Table Of ContentHigh-accuracy waveforms for binary black hole inspiral, merger, and ringdown
Mark A. Scheel,1 Michael Boyle,1 Tony Chu,1 Lawrence E. Kidder,2 Keith D. Matthews,1 and Harald P. Pfeiffer1
1Theoretical Astrophysics 130-33, California Institute of Technology, Pasadena, California 91125
2Center for Radiophysics and Space Research, Cornell University, Ithaca, New York 14853
(Dated: January 20, 2009)
The first spectral numerical simulations of 16 orbits, merger, and ringdown of an equal-mass
nonspinningbinaryblackholesystemarepresented. Gravitationalwaveformsfromthesesimulations
haveaccumulatednumericalphaseerrorsthroughringdownof<∼0.1radianwhenmeasuredfromthe
beginningofthesimulation,and<∼0.02radianwhenwaveformsaretimeandphaseshiftedtoagree
atthepeakamplitude. Thewaveformseenbyanobserveratinfinityisdeterminedfromwaveforms
computed at finite radii by an extrapolation process accurate to <∼ 0.01 radian in phase. The
phasedifference between this waveform at infinity and the waveform measured at a finite radius of
r=100M isabouthalfaradian. TheratiooffinalmasstoinitialmassisM /M =0.95162±0.00002,
9 f
and thefinal black hole spin is S /M2 =0.68646±0.00004.
0 f f
0
2 PACSnumbers: 04.25.D-,04.25.dg,04.30.-w,04.30.Db,02.70.Hm
n
a
I. INTRODUCTION detectiontemplates,itismostlikelyinadequateforLIGO
J
parameterestimationandisfarfromwhatisrequiredfor
0
LISA data analysis [42].
2 Beginning with the groundbreaking binary black hole
evolutions of Pretorius [1] and the development of the One approach to increasing the accuracy and effi-
] ciency of simulations is to adopt more efficient numer-
c moving puncture method [2, 3], it has recently become
ical methods. In particular, a class of numerical tech-
q possible to solveEinstein’s equationsnumericallyfor the
niques known as spectral methods holds much promise.
- inspiral, merger, and ringdown of two black holes in a
r For smooth solutions, the errors produced by spectral
g binary orbit. Already these simulations have provided
methods decrease exponentially as computational re-
[ tests of post-Newtonian approximations [4, 5, 6, 7, 8, 9,
sources are increased, whereas the errors of finite dif-
10, 11, 12, 13, 14], have allowed initial exploration of
2 ference methods, the methods used by the majority of
v the orbital dynamics of spinning binaries [15, 16, 17, 18,
binaryblackholesimulations,decreasepolynomially. In-
7 19, 20], have determined the recoil velocity of the final
deed, spectral methods have been used to produce very
6 black hole when the masses are unequal [21, 22, 23, 24],
accurate initial data for binary black holes and neutron
7 andhaveledto the discoveryofdramaticallylargerecoil
1 stars[43,44,45,46,47,48,49,50,51,52,53,54,55,56],
velocity from certain spin configurations [18, 25, 26, 27,
. and they have been used to produce the longest and
0 28, 29, 30, 31, 32, 33, 34, 35, 36, 37].
most accurate binary black hole inspiral simulation to
1
Waveforms from these numerical simulations are im- date [9, 57].
8
0 portant for gravitational-wave detectors such as LIGO However, a key difficulty with time-dependent spec-
: and LISA. This is not only because detected waveforms tral binary black hole simulations has been handling the
v
can be compared with numerical models to measure as- merger of the two holes. For example, the spectral sim-
i
X trophysicalpropertiesofthesourcesofgravitationalradi- ulations described in [9, 12, 57] are very accurate and
r ation,butalsobecausethedetectionprobabilityitselfcan efficient, but they follow only the inspiral of the two
a beincreasedviathetechniqueofmatchedfiltering[38],in black holes, and fail just before the holes merge. This
whichnoisydataareconvolvedwithnumericaltemplates is sufficient for some applications, such as comparing
to enhance the signal. post-Newtonian formulae with numerical results during
However, binary black hole simulations are time con- the inspiralandfinding accurateanalytic templates that
suming: a single simulation following approximately 10 match the numerical inspiral waveforms [9, 12], but for
orbits, merger, and ringdown typically requires a few most purposes the merger is the most crucial part of the
weeks of runtime on approximately 50 or 100 processors process: for instance the gravitational-wave emission is
ofaparallelsupercomputer,andtypicallysuchasimula- the strongest during merger, and details of the merger
tion produces waveforms of only modest accuracy. This determine the recoil velocity of the final black hole.
largecomputationalexpenseprecludes,forexample,pro- In this paper we present a spectral binary black hole
ducing a full template bank ofnumericalwaveformscov- simulation that follows sixteen orbits of the binary plus
eringtheentireparameterspaceofblackholemassesand mergerandringdownofthemergedblackhole. InSec.II
spins. Hencetherehasbeenmuchinterestinconstruction we describe the equations,gauge conditions,and numer-
of phenomenological analytical waveforms [7, 39, 40, 41] icalmethodsweusetosolveEinstein’sequations;inpar-
that can be computed quickly and are calibrated by a ticular,Secs.IIC andIIDdescribe changesto ourgauge
small number of numerical simulations. While the accu- conditions that allow simulation of the merger, and our
racyoftypicalsimulationsissufficientforcreatingLIGO methodforextendingtheevolutionthroughringdown. In
2
Sec. III we discuss extraction of the gravitational wave- only internally in the code as a means of treating mov-
formfromthesimulation,includingtheprocessofextrap- ing holes with a fixed domain. Therefore all coordinate
olatingthe waveformto infinity. Sec.III alsoincludes an quantities (e.g. black hole trajectories, wave-extraction
estimate of the uncertainty in the waveform from sev- radii) mentioned in this paper are inertial-frame values
eral sources. Finally, in Sec. IV we discuss outstanding unless explicitly stated otherwise.
difficulties and future improvements. As described in [9], the mapping between inertial and
comovingcoordinatesfor the inspiral, expressedin polar
coordinates relative to the center of mass of the system,
II. SOLUTION OF EINSTEIN’S EQUATIONS is
r′2
A. Initial data r = a(t)+(1 a(t)) r′, (1)
(cid:20) − R′2(cid:21)
0
′
The initial data describe two nonspinning black holes, θ = θ , (2)
each with Christodoulou mass M/2, in quasicircular or- φ = φ′+b(t), (3)
bit with low eccentricity. The initial data are exactly as
describedin Ref.[9]. Briefly,initial data areconstructed wherea(t)andb(t)arefunctionsoftime,andR′ isacon-
0
within the conformal thin sandwich formalism [58, 59] stantusuallychosentoberoughlytheradiusoftheouter
using a pseudospectral elliptic solver [49]. We employ boundary in comoving coordinates. Here primes denote
quasiequilibrium boundary conditions [50, 60] on spher- the comoving coordinates. For the choice R′ = , the
0 ∞
ical excision boundaries, choose conformal flatness and mapping is simply a rotationby b(t) plus anoverallcon-
maximal slicing, and use Eq. (33a) of Ref. [53] as the traction given by a(t). The functions a(t) and b(t) are
lapse boundary condition. The spins of the black holes determined by a dynamical control system as described
aremadeverysmall( 10−7)viaanappropriatechoiceof in Ref. [57]. This control system dynamically adjusts
∼
the tangentialshift at the excision surfaces,as described a(t)andb(t)sothatthe centersofthe apparenthorizons
in[53]. Finally,the initialorbitaleccentricityis tunedto remain stationary in the comoving frame. Note that the
averysmallvalue( 5 10−5)usingtheiterativeproce- outerboundaryofthecomputationaldomainisatafixed
∼ × ′
dure described in Ref. [9], which is an improved version comoving radius R , so the inertial-coordinate radius
max
of the procedure of Ref. [61]. of the outer boundary R (t) is a function of time.
max
Thegaugefreedominthegeneralizedharmonicsystem
is fixed via a freely specifiable gauge source function H
a
B. Evolution of the inspiral phase that satisfies the constraint
0= Γ b+H , (4)
Theevolutionofthefirst 15binaryorbitsisidentical a ab a
C ≡
∼
to the simulation presented in Ref. [9]. We describe it
where Γa are the spacetime Christoffel symbols. To
here briefly in order to facilitate the presentation of our bc
choose this gauge source function, we first define a new
methodforcontinuingthe evolutionthroughmergerand
quantityH˜ thathasthefollowingtwoproperties: (1)H˜
ringdown, which is described in Secs. IIC and IID. a a
transforms like a tensor, and (2) in inertial coordinates
The Einstein evolution equations are solved with the
H˜ = H . We choose H so that the constraint equa-
pseudospectral evolution code described in Ref. [57]. a a a
This code evolves a first-order representation [62] of the tion (4) is satisfied initially, and we demand that H˜a′ is
generalized harmonic system [63, 64, 65]. We handle constant in the moving frame, i.e., that ∂t′H˜a′ =0.
the singularitiesby excising the blackhole interiorsfrom
the computational domain. Our outer boundary condi-
tions[62,66,67]aredesignedtopreventtheinfluxofun- C. Extending inspiral runs through merger
physical constraint violations [68, 69, 70, 71, 72, 73, 74]
and undesired incoming gravitational radiation [75, 76], Iftheinspiralrunsdescribedaboveareallowedtocon-
while allowing the outgoing gravitational radiation to tinue without any modification of the algorithm, then
pass freely through the boundary. as the binary approaches merger, the horizons of the
We employ the dual-frame method described in black holes become extremely distorted and the dynam-
Ref. [57]: we solve the equations in an “inertial frame” ical fields begin to develop sharp (but numerically con-
that is asymptotically Minkowski, but our domain de- vergent) features near each hole. These features grow
composition is fixed in a “comoving frame” that rotates rapidly in time, eventually halting the simulation before
with respect to the inertial frame and also shrinks with merger. This is due to a gauge effect: The gauge condi-
respect to the inertial frame as the holes approach each tion used during the inspiral, namely fixing H in time
a
other. The positions of the holes are fixed in the comov- in the comoving frame, was chosen based on the idea
ing frame; we account for the motion of the holes by dy- that each black hole is in quasiequilibriumin this frame.
namicallyadjustingthecoordinatemappingbetweenthe Oncetheblackholesbegintointeractstrongly,thisgauge
two frames. Note that the comoving frame is referenced condition no longer allows the coordinates to sufficiently
3
reacttothe changinggeometry,andcoordinatesingular- gauge does not change too rapidly. We choose
ities develop.
f(x,t)=g(x,t) = (2 e−(t−tg)/σ1)
Therefore we must modify our gauge conditions in or- −
der to handle merger. Because the inspiral gauge works (1 e−(t−tg)2/σ22)e−r′2/σ32, (8)
× −
sowellbefore merger,wechooseto remainin thatgauge
where r′ is the coordinate radius in comoving coordi-
until some time t = t , and then we change (smoothly)
g nates, and the constants are σ 17.5M, σ 15M,
to a new gauge. 1 ∼ 2 ∼
and σ 40M. Here M is the sum of the initial
We have experimented with several gauge condi- 3 ∼
Christodoulou masses of the two holes.
tions[77],butsofarthesimplestgaugechoicethatworks,
Equation (5) is a second-order hyperbolic equation,
and the one used in the simulations presented here, is
whichweevolveinfirst-orderformbydefiningnewfields
based on the gauge treatment of Pretorius [1, 65, 78]:
ΠH and ΦH, representing (up to the addition of con-
We promote the gauge source function H to an inde- a ia
a straints) the appropriate time and space derivatives of
pendent dynamical field that satisfies
H , respectively:
a
c H =Q (x,t,ψ )+ξ tb∂ H , (5) ΠH = tb∂ H , (9)
∇ ∇c a a ab 2 b a a − b a
ΦH = ∂ H . (10)
ia i a
where c isthecurvedspacescalar waveoperator(i.e.
c
each co∇m∇ponent of H is evolved as a scalar), ψ is the Therepresentationofwaveequationsofthistypeinfirst-
a ab
spacetime metric, and ta is the timelike unit normal to orderformiswellunderstood,seee.g.,Refs.[62, 80]; the
the hypersurface. The driving function Q is result for Eq. (5) is
a
∂ H = NΠH +NkΦH, (11)
1 N t a − a ka
Qt = f(x,t)ξ1 N−η , (6) ∂tΠHa = Nk∂kΠHa −Ngki∂kΦHia−γ2HNk∂kHa
N + γHNkΦH +N(Γkj gkj∂ N)ΦH
Q = g(x,t)ξ i. (7) 2 ka j − j ka
i 3N2 + NKΠH +Q , (12)
a a
∂ ΦH = Nk∂ ΦH N∂ ΠH +γHN∂ H
HereN andNiarethelapsefunctionandtheshiftvector, t ia k ia− i a 2 i a
η, ξ1, ξ2, and ξ3 are constants, and f(x,t) and g(x,t) − ΠHa ∂iN +ΦHka∂iNk−γ2HNΦHia, (13)
areprescribedfunctionsofthespacetimecoordinates(we
where g is the spatial metric and K is the trace of the
describe our choices for these objects below). ij
extrinsic curvature. We choose the constraint-damping
Equation (5) is a damped, driven wave equation with
parameter γH to be γH =4/M.
damping parameter ξ and driving function Q . The 2 2
2 a Theseequationsaresymmetrichyperbolic,andrequire
driving term Q in Eq. (6) was introduced by Preto-
t boundary conditions on all incoming characteristicfields
rius [1, 65] to drive the lapse function towards unity so
atallboundaries. The characteristicfieldsforEqs.(11)–
as to prevent it from becoming small. The driving term
(13) in the direction of a unit spacelike covector n are
i
Q is new; it drives the shift vector towards zero near
i
the horizons. This causes the horizons to expand in co- UH± = ΠH niΦH γHH , (14)
a a ± ia− 2 a
ordinate space, and has the effect of smoothing out the ZH1 = H , (15)
dynamical fields near the horizon and preventing gauge a a
ZH2 = (δk n nk)ΦH. (16)
singularities from developing. A different gauge choice ia i − i ka
that causes similar coordinate expansion of the horizons The (coordinate) characteristic speeds for UH±, ZH1,
a a
was introduced in Ref. [79]. Care must be taken so that and ZH2 are N n Ni, 0, and n Ni, respectively.
thehorizonsdonotexpandtooquicklyrelativetotheex- At tihae exci±sion−bouindaries all c−hariacteristic fields are
cision boundaries; otherwise the characteristic fields will outgoing (i.e. into the holes) or nonpropagating, so no
fail to be purely outgoing (into the holes) at the exci- boundaryconditionsarenecessaryandnoneareimposed.
sionboundaries,andexcisionwillfail. We findthatwith At the outer boundary, we must impose boundary con-
appropriate choices of ξ1, ξ3, f(x,t), and g(x,t) as de- ditions on UH− and ZH2. Define
a ia
scribed below, the horizons expand gradually and not
too rapidly. Dt(UaH±) ≡ ∂tΠHa ±ni∂tΦHia−γ2H∂tHa, (17)
ξ2F=or1t0h,earnudnξs3p=res0e.n4t.edThheerfeuwncetciohnososfe(xη,=t)4a,nξd1g=(x0,.1t), DDt((ZZaHH12)) ≡ ∂(δtkHa,n nk)∂ ΦH, ((1189))
in Eqs. (6) and (7) are chosen based on two criteria: the t ia ≡ i − i t ka
first is that the driving terms Q are nonzero only near where the time derivatives on the right-hand side are
a
the blackholes where they are needed; if these terms are evaluated using Eqs. (11)–(13). Then we impose the fol-
nonzerointhe wave-extractionzonethey leadto compli- lowing boundary conditions:
cated gauge dynamics in this region, making waveform
∂ UH− = γHD (ZH1), (20)
extractiondifficult. Thesecondcriterionisthatthedriv- t a − 2 t a
ing terms are turned on in a gradualmanner so that the ∂ ZH2 = D (ZH2)+2n Nknj∂ ΦH . (21)
t ia t ia k [i j]a
4
Equation (20) is the outgoing-wave boundary condition in detail.
described in detail in Ref. [80]. Equation (21) ensures Our new computational domain contains only a single
that violations of the artificial constraint C ΦH excised region, and is much simpler than the one used
ia ≡ ia −
∂ H = 0 do not enter the domain through the bound- until merger. It consists only of nested spherical-shell
i a
ary;it is the directanalogueofthe constraint-preserving subdomains that extend from a new excision boundary
boundaryconditionweapplytotheanalogousvariablein R′′ , chosen to be slightly inside the common apparent
min
the generalizedharmonicformulationofEinstein’s equa- horizon, to an outer boundary R′′ that coincides with
max
tions, Eq. (65) of Ref. [62]. the outer boundary of the old domain.
Note that Eqs. (11)–(13) involve only first derivatives Tounderstandhowwechooseournewcomovingframe,
of the spacetime metric, and similarly, the generalized first recall that in the dual-frame technique [57], the co-
harmonicEinsteinequationsinvolveonlyfirstderivatives moving frame is the one in which the computational do-
of H . Therefore, adding Eqs. (11)–(13) to the system main is fixed, the inertial frame is the one in which the
a
does not change the hyperbolicity or characteristicfields coordinates are Minkowski-like at infinity, and the two
ofthegeneralizedharmonicEinsteinequations,sowecan frames are related by a mapping that is chosen so that
impose the sameboundary conditionsonthe generalized thecomputationaldomaintracksthemotionoftheblack
harmonic variables as we do during the inspiral, as de- holes. Let xa represent the inertial coordinates (which
scribed in Refs. [57, 66]. are the same before and after merger), let x′a represent
Equations(11)–(13)requireasinitialdatathevaluesof the old comoving coordinates, and let x′′a represent the
H andΠH att=t . Thesequantities canbe computed new comoving coordinates. The mapping between x′a
a a g
fromthe gaugechoiceusedduringthe inspiralfort t , and xa is given by Eqs. (1)–(3). The mapping between
g
so we choose them to be continuous at t=t . ≤ x′′a and xa is chosen to be
g
Note that Eqs. (11)–(13) and the boundary condi-
tions (20) and (21) are written in the inertial coordinate r = r˜ 1+sin2(πr˜/2R′′ )
(cid:20) max
system. The equations are actually solvedin the comov-
ing coordinate system using the dual-frame method de- R′ R′3
A(t) max +(1 A(t)) max 1 (,22)
scribed in Ref. [57]. ×(cid:18) R′′ − R′′ R′2 − (cid:19)(cid:21)
max max 0
With the modifications to the gauge conditions de-
ℓmax ℓ
scribed here, the evolution of the binary can be tracked r˜ = r′′ q(r′′) λ (t)Y (θ′′,φ′′), (23)
ℓm ℓm
up until (and shortly after) the formation of a common −
Xℓ=0mX=−ℓ
horizon that encompasses both black holes. Because of ′′
θ = θ , (24)
the more rapiddynamics and the distortions ofthe hori-
′′
zonsduringthe merger,wetypicallyincreasethe numer- φ = φ +B(t), (25)
ical resolution slightly when we make these changes to
whereR′ istheouterboundaryofthepremergercom-
the gauge conditions (this is the difference between the max
putational domain in the old comoving coordinates, and
first and second entry in the N column in Table I).
pts q(r′′), A(t), B(t), and λ (t) are functions we will now
Afterthecommonhorizonforms,theproblemreducesto ℓm
discuss.
evolving a single highly distorted dynamical black hole,
First we describe the angular map: The function B(t)
rather than two separate black holes. We change the al-
ischosensothatthenewcomovingframeinitiallyrotates
gorithm to take advantage of this, as described in the
withrespecttotheinertialframe,butthisrotationslows
next section.
to a halt after a short time. In particular,
B(t)=B +(B +B (t t ))e−(t−tm)/τB, (26)
0 1 2 m
D. Evolution from merger through ringdown −
where the constants B , B , and B are chosen so that
0 1 2
We make three main changes to our evolution algo- B(t) matches smoothly onto b(t) from Eq. (3): B(tm)=
rithmoncewedetectacommonapparenthorizon. First, b(t ), B˙(t ) = b˙(t ), and B¨(t ) = ¨b(t ). Here t
m m m m m m
because there is now only one black hole and not two, is the time at which we transition to the new domain
we interpolate all variables onto a new computational decomposition. The constant τ is chosen to be on the
B
domain that contains only a single excised region. Sec- order of 20M.
ond, we choose a new comoving coordinate system (and The radial map is a composition of two individual
a correspondingmapping to inertialcoordinates)sothat maps: Eqs. (22) and (23). The purpose of Eq. (22)
the new excision boundary tracks the shape of the (dis- is to match the outer boundary of the new domain
torted, rotating, pulsating) apparent horizon in the in- smoothly onto that of the old domain, while far from
ertial frame, and so that the outer boundary behaves the outer boundary Eq. (22) approaches the identity.
smoothly in time. Third, we modify the gauge condi- We have found that without the use of Eq. (22), the
tions so that the shift vector is no longer driven towards (inertial-coordinate) location of the boundary changes
zero, allowing the solution to eventually relax to a time- nonsmoothlyatt=t ,therebygeneratingaspuriousin-
m
independent state. We now discuss these three changes going gauge pulse that spoils waveform extraction. The
5
function A(t) is Run Rm′ ax Rm′′ax R0′ Npts CPU-h CPU-h/T
A(t)=A +(A +A (t t ))e−(t−tm)/τA, (27) 30c1/N4 462 462 698 (573,593,573) 8,800 2.0
0 1 2 m
− 30c1/N5 462 462 698 (623,663,633) 15,000 3.4
where the constants A0, A1, and A2 are chosen so that 30c1/N6 462 462 698 (673,733,703) 23,000 5.3
A(t) matches smoothly onto a(t) from Eq. (1): A(t )=
m 30c2/N6 722 96 ∞ (713,763,633) 25,000 5.7
a(t ),A˙(t )=a˙(t ),andA¨(t )=a¨(t ). Theconstant
m m m m m
τA is chosen to be on the order of 5M. TABLE I: Outer boundary parameters, collocation points,
The other piece of the radial map, Eq. (23), is chosen and CPU usage for several zero-spin binary black hole evo-
so that the apparent horizon is nearly spherical in the lutions. The first column identifies the inspiral run in the
new comoving coordinates x′′a. The function q(r′′) is nomenclature of Ref. [9]. Npts is the approximate number
of collocation points used to cover the entire computational
q(r′′)=e−(r′′−R′A′H)3/σq3, (28) domain. The three values for Npts are those for the inspiral,
where R′′ is the radius of the apparent horizon in co- merger, and ringdown portions of the simulation, which are
AH described in Sections IIB, IIC, and IID, respectively. The
moving coordinates, and σq is a constant of order 20M. outer boundary parameters R′ , R′′ and R′, as well as
This function q(r′′) ensures that the piece of the radial max max 0
runtimesT,areinunitsoftheinitialChristodouloumassM
map represented by Eq. (23) acts only in the vicinity of ofthesystem,whichprovidesanaturaltimeandlengthscale.
the merged hole and not in the exterior wave-extraction
region.
Wenowdiscussthechoiceofthefunctionsλ (t)that
ℓm
appear in Eq. (23). Given the known location of the where σ4 = 7M. The modification of g(x,t) turns off
apparent horizon in inertial coordinates, the λ (t) de- the term in the gauge evolution equations that drives
ℓm
termine the shape of the apparent horizon in comoving the shift to zero near the holes. Before merger, it is ad-
coordinates. At t = t , we choose these quantities so vantageous to have the shift driven to zero so that the
m
that the apparent horizon is spherical (up to spherical horizonsexpandincoordinatespaceandsothatgrowing
harmoniccomponentℓ=ℓ )incomovingcoordinates: gauge modes remain inside the common horizon. After
max
thatis,ifthecomoving-coordinateradiusoftheapparent merger,however,it is no longer desirable for the horizon
horizon as a function of angles is written as to expand, since this would prevent the solution from
eventually settling down to a time-independent state in
ℓmax ℓ
′′ ′′ ′′ ′′ ′′ which the horizonis stationarywith respect to the coor-
r (θ ,φ ) Q (t)Y (θ ,φ ), (29)
AH ≡ ℓm ℓm dinates.
Xℓ=0 mX=−ℓ
Tosummarize,thestepsinvolvedinthetransitionfrom
then for 1 ℓ ℓmax we choose λℓm(tm) so that evolving a binary black hole spacetime to evolving a
≤ ≤
Qℓm(tm) = 0. In addition, we choose λ00(tm) = 0; this merged single black hole spacetime are as follows: (1)
determines RA′′H. For t > tm, λℓm(t) are determined Find the common apparent horizon in the inertial frame
by a dynamical feedback control system identical to the at several times near t = t . (2) Solve for the λ (t )
m ℓm m
one described in Ref. [57], which adjusts these functions that make the horizon spherical in the comoving frame,
so that the apparenthorizonis drivento a sphere (up to and simultaneously solve for R′′ . (3) Choose the in-
AH
sphericalharmoniccomponentℓ=ℓmax)incomovingco- ner boundary of the new computational domain R′′ to
min
ordinates. This dynamical feedback control allows us to be slightly less than R′′ , and choose the outer bound-
AH
freelychoosethefirstandsecondtime derivativesofλℓm ary Rm′′ax [for sufficiently small a(tm) it is necessary to
at t = tm. Simply choosing these to be zero causes the choose R′′ <R′ so that the mapping (22) is invert-
max max
control system to oscillate wildly before settling down, ible]. At this point the computational domain and the
and unless the time step is very small, these oscillations mapping (22)– (25) have been determined. (4) Interpo-
are large enough that the excision boundary crosses the late all dynamical variables from the old computational
horizon and our excision algorithm fails. So instead, we domain onto the new one. This interpolation is done via
obtainthetimederivativesofλℓmbyfindingtheapparent thespectralexpansionintheolddomain,soitintroduces
horizon at several times surrounding t = tm, comput- no additional error. (5) Modify the gauge source evolu-
ing λℓm at these times, and finite-differencing in time. tion equations so that the shift is no longer driven to
For the equal-mass zero-spin merger presented here, in zero. (6) Continue the evolution on the new computa-
Eq.(23)itsuffices to sumonlyoverevenℓandm andto tional domain. All of these steps can be automated.
choose ℓ =6.
max
The last change we make before continuing the simu-
lation past merger is to modify the functions f(x,t) and
E. Properties of the numerical solution
g(x,t), which before merger were given by Eq. (8), to
f(x,t) = (2 e−(t−tg)/σ1) In Table I we list outer boundary parameters, resolu-
−
(1 e−(t−tg)2/σ22)e−r′′2/σ32, (30) tions, and run times of several runs we have done using
× − the algorithm described above. Three of these runs are
g(x,t) = f(x,t)e−(t−tm)2/σ42, (31) identical except for numerical resolution, and the fourth
6
-2
3960 10 -2
4000 30c1 10
|| C ||
t=t
m -4
t/M -4 10
10
30c2
3000 t=t 3900 3920 3940 3960
30c1 3920 g
-6
10
t/M 240 260
r/M
2000 || C ||
-4
10 N4
3960 (unnormalized) N5
30c2
N6
-6
t=t 10
1000 m
t/M
-8
10
t=t
3920 g -10
10
0 0 1000 2000 3000 4000
0 250 500 750 100 150 200
r/M r/M t/M
FIG. 1: Spacetime diagram showing the spacetime volume FIG. 2: Constraint violations of run 30c1. The top panel
simulatedbythenumericalevolutionslistedinTableI. Each shows the L2 norm of all constraints, normalized by the L2
curve represents the worldline of the outer boundary for a norm of the spatial gradients of all dynamical fields. The
particularsimulation. Themagnifiedviewsontherightshow bottom panel shows the same data, but without the normal-
that the outer boundary moves smoothly near merger. The izationfactor. TheL2normsaretakenovertheportionofthe
transition times tg =3917M and tm =3940M are indicated computational volume that lies outside apparent horizons.
on theright panels.
At t=3917M (t=3927M for resolutionN4), the gauge
conditions are changed (cf. Sec. IIC) and the resolution
is performedona differentdomainwith a differentouter
is increased slightly (compare the first and second entry
boundarylocation. Asdiscussedabove,theouterbound-
in the N column in Table I). Because of the change
ary of our simulation varies in time because of the dual- pts
ofresolution,the constraints droprapidly by almosttwo
frameapproachweusetofollowtheblackholes. Figure1
orders of magnitude, but then they begin to grow again.
is a spacetime diagram illustrating the region of space-
The transition to a single-hole evolution (cf. Sec. IID)
time being evolved in our simulation.
occurs at t = 3940M (t = 3948M for resolution N4).
We do not explicitly enforce either the Einstein con-
At this time the constraint norm drops by about an or-
straintsorthesecondaryconstraintsthatarisefromwrit-
derofmagnitude because the regioninwhichthe largest
ing the system in first-order form. Therefore, examin-
constraint violations occur—the interior of the common
ing how well these constraints are satisfied provides a
horizon—is newly excised.
useful consistency check. Figure 2 shows the constraint
After the binary proceeds through inspiral, merger,
violations for run 30c1. The top panel shows the L2
and ringdown, it settles down to a final stationary black
norm of all the constraint fields of our first-order gener-
hole. In our simulation this final state is not expressed
alized harmonic system, normalized by the L2 norm of
inanystandardcoordinatesystemusedto describe Kerr
thespatialgradientsofthedynamicalfields(seeEq.(71)
spacetime, but nevertheless the final mass and spin of
of Ref. [62]). The bottom panel shows the same quan-
the hole can be determined. The area A of the apparent
tity, but without the normalization factor [i.e., just the
horizon provides the irreducible mass of the final black
numerator of Eq. (71) of Ref. [62]]. The L2 norms are
hole,
takenoverthe portionofthe computationalvolumethat
liesoutsideapparenthorizons. Atearlytimes,t<500M,
M = A/16π, (32)
irr
theconstraintsconvergeratherslowlywithresolutionbe-
p
cause the junk radiationcontains high frequencies. Con- which we find to be M /M =0.88433 0.00001,where
irr
vergenceismorerapidduringthe smoothinspiralphase, M isthesumoftheinitialirreduciblem±assesoftheblack
after the junk radiation has exited through the outer holes. TheuncertaintyinM /M isdeterminedfromthe
irr
boundary. differencebetweenruns30c1/N6,30c1/N5,and30c2/N6,
The constraints increase as the holes approach each so it includes only uncertainties due to numerical reso-
other and the solution becomes increasingly distorted. lution and outer boundary location. We have verified
7
Initial orbital eccentricity: e ∼ 5×10−5 Penrose scalar Ψ4, following the same procedure as in
Initial spin of each hole: Si/M2 <∼ 10−7 Refs. [61, 86]. To summarize, we compute
Time of evolution: T/M = 4330
Final Christodoulou mass: M /M = 0.95162±0.00002 Ψ4 = Cαµβνℓµℓνm¯αm¯β, (34)
f −
Final spin: S /M2 = 0.68646±0.00004
f f
where
TABLE II: Physical parameters describing the equal-mass
1
nonspinning binary black hole evolutions presented here. ℓµ = (tµ rµ), (35a)
The dimensionful quantity M is the initial sum of the √2 −
Christodoulou masses of the black holes. Uncertainty esti- µ
1 ∂ 1 ∂
matesincludenumericaluncertaintiesandtheeffectsofvary- mµ = +i . (35b)
ing theouter boundary location. √2r (cid:18)∂θ sinθ∂φ(cid:19)
Here (r,θ,φ) denote the standard spherical coordinates
that the uncertainty due to the finite resolution of our intheinertialframe,tµ isthetimelikeunitnormaltothe
apparent horizon finder is negligible. spatialhypersurface,andrµ istheoutward-pointingunit
The final spin S of the black hole can be computed normal to the extraction sphere. We then expand Ψ in
f 4
by integrating a quasilocal angular momentum density termsofspin-weightedsphericalharmonicsofweight 2:
−
over the final apparent horizon [81, 82]. Our imple-
mentation of this method is described in detail in Ap- Ψ4(t,r,θ,φ)= Ψl4m(t,r)−2Ylm(θ,φ), (36)
pendix A of [55]. Furthermore, an alternative method of
Xlm
computing the final spin, which is based on evaluating
the extremal values of the 2-dimensional scalar curva- where the Ψlm are expansion coefficients defined by this
tureontheapparenthorizonandcomparingthesevalues 4
equation.
to those obtained analytically for a Kerr black hole, is
Note that our choice of mµ is not exactly null nor
also described in [55]. Using these measures, we deter-
exactly of unit magnitude at finite r, as is required by
mine the dimensionless spin of the final black hole to
the standard definition. The resulting Ψlm computed at
be S /M2 = 0.68646 0.00004, where the uncertainty 4
f f ± finite r will therefore disagree with the waveforms ob-
is dominated by the difference between runs 30c1/N6
served at infinity. Our definition does, however, agree
and 30c1/N5 rather than by the differences between dif- with the standard definition of Ψlm as r . Because
ferent methods of measuring the spin. Here M is the 4 → ∞
f we extrapolate the extracted waves to find the asymp-
Christodoulou mass of the final black hole,
totic radiation field (see Sec. IIIC), these tetrad effects
S2 shouldnotplayarole: RelativeerrorsinΨl4m introduced
M2 =M2 + f . (33) by using the simple coordinate tetrad fall off like pow-
f irr 4M2
irr ers of M/r, and thus should vanish after extrapolating
We find that the ratio of the final to initial black hole to obtain the asymptotic behavior. More careful treat-
mass is M /M =0.95162 0.00002. The mass and spin ment of the extraction method—such as those discussed
f
ofthe finalhole areconsist±entwith thosefound by other in Refs. [87, 88, 89]—may improve the quality of extrap-
groups [2, 4, 83, 84, 85]. Physical parameters describing olationandwouldbe interestingtoexploreinthe future.
the evolutions are summarized in Table II. Inthis paper,we focusonthe dominant(l,m)=(2,2)
mode. Following common practice (see e.g. [84, 85]), we
split the extracted waveform into real phase φ and real
III. COMPUTATION OF THE WAVEFORM amplitude A, defined by
The numerical solution of Einstein’s equations ob- Ψ22(r,t)=A(r,t)e−iφ(r,t). (37)
4
tained using the methods described above yields the
spacetime metric and its first derivatives at all points The gravitational-wavefrequency is given by
inthecomputationaldomain. Inthissectionwedescribe
how this solution is used to compute the key quantity dφ
ω = (38)
relevant for gravitational-wave observations: the gravi- dt
tational waveform as seen by an observer infinitely far
from the source. The minus sign in Eq. (37) is chosen so that the phase
increases in time and ω is positive.
The (l,m) = (2,2) waveform, extracted at a single
A. Waveform extraction radius for run 30c1/N6, is shown in Fig. 3. The short
pulse at t 200M is caused by imperfect initial data
∼
Gravitationalwavesare extractedfrom the simulation thatarenotpreciselyinequilibrium;thispulseisusually
on a sphere of coordinate radius r using the Newman- referred to as “junk radiation”.
8
0.001
0.05 1
22ΨM)4 0 0 ans) N4-N6
e(r di
R a10-2 N5-N6
-0.05 r
-0.001 (
0 1000 2000 3000 4000 4100 4200 φ
t/M t/M ∆ 30c2 - 30c1
-4
FIG.3: Gravitationalwaveformextractedatfiniteradiusr= 10
225M, for thecase30c1/N6 in TableI. Theleft panelzooms 1
in on theinspiral waveform, and theright panel zooms in on
themerger and ringdown.
A -2 N4-N6
10
/
A
B. Convergence of extracted waveforms
∆ N5-N6
-4
Inthissectionweexaminetheconvergenceofthegrav- 10 30c2 - 30c1
itational waveforms extracted at fixed radius, without
extrapolation to infinity. This allows us to study the
0 1000 2000 3000 4000
behavior of our code without the complications of ex-
t/M
trapolation. The extrapolationprocessandthe resulting
extrapolated waveforms are discussed in Sec. IIIC.
FIG. 4: Convergence of waveforms with numerical resolution
Figure 4 shows the convergence of the gravitational-
andouterboundarylocation. Shownarephaseandamplitude
wave phase φ and amplitude A with numerical resolu- differencesbetweennumericalwaveformsΨ22computedusing
4
tion. Forthisplot,thewaveformwasextractedatafixed different numerical resolutions. Shown also is the difference
inertial-coordinate radius of r = 60M. This fairly small betweenourhighest-resolution waveformsusingtwodifferent
extraction radius was chosen to allow a comparison of outer boundary locations. All waveforms are extracted at
the simulations30c1and30c2. Eachsolidline inthe top r = 60M, and no time shifting or phase shifting is done to
panel showsthe absolute difference between φ computed align waveforms.
at some particular resolution and φ computed from our
highest-resolution run, labeled 30c1/N6 in Table I. The
solid curves in the bottom panel similarly show the rela-
tive amplitude differences. When subtracting results at
different resolutions, no time or phase adjustment has sible only for extractionradiir<75M. Comparingruns
beenperformed. Thenoiseatearlytimesisdueto“junk 30c1 and 30c2 provides an esti∼mate of the uncertainty
radiation”generatedneart=0. Whilemostofthisradi-
in the waveform due to outer boundary effects such as
ation leaves through the outer boundary after one cross-
imperfect boundary conditions that might reflect outgo-
ing time, some remains visible for a few crossing times1.
ing waves. From Fig. 4 we estimate this uncertainty to
The plots show that the phase difference accumulated
be 0.03 radiansin phase and half a percent in amplitude
over 16 orbits plus merger and ringdown is less than 0.1
(when no time shift is applied).
radians for our medium resolution, and the relative am-
plitudedifferencesarelessthan0.015;thesenumberscan
Figure 5 is the same as Fig. 4 except eachwaveformis
betakenasanestimateofthenumericaltruncationerror
time shiftedandphaseshiftedsothatthe maximumam-
of our medium resolution run.
plitude of the wave occurs at the same time and phase.
Also shown as a dotted curve in each panel of Fig. 4
This type of comparison is relevant for analysis of data
is the difference between our highest-resolution run,
from gravitational-wave detectors: when comparing ex-
30c1/N6, and a similar run but with a different outer
perimental data with numericaldetection templates, the
boundary location, 30c2/N6. The 30c2 run initially has
template will be shifted in both time and phase to best
amoredistantouterboundarythan30c1,butduringthe
matchthedata. Forthistypeofcomparison,Fig.5shows
inspiraltheouterboundarymovesrapidlyinward,asseen
thatthenumericaltruncationerrorofourmediumresolu-
in Fig. 1, so that extraction of the full waveform is pos-
tionrunislessthan0.01radiansinphaseand0.1percent
in amplitude for t>1000M. At earlier times, the errors
aresomewhatlargerandaredominated by residualjunk
radiation. Ouruncertaintyduetoouterboundaryeffects
1 The junk radiation at early times is discussed in more detail
is similar to that in Fig. 4: about 0.02 radians in phase
in Ref. [9] (specifically, just before Eq. (9) and in the third
and half a percent in amplitude. Boundary effects are
paragraphofSec. IIE),whichpresentstheexactsamewaveform
asshownherebutwithoutmergerandringdown. most prominent during the ringdown.
9
-1 and the tortoise-coordinate radius [90] is
10 N4-N6
)
ns r∗ =r +2M ln rareal 1 . (40)
a areal ADM (cid:18)2M − (cid:19)
i ADM
ad10-2 N5-N6
r Here MADM is the ADM mass of the initial data,
(
φ 30c2 - 30c1 and rareal = A/4π, where A is the measured (time-
∆ dependent) arpea of the extraction sphere. If we were
10-3 to choose ts to be simply the coordinate time t, then
the retarded time coordinate u would fail to be null,
largely because the lapse function in our simulation is
-2
10 time-dependentanddiffersfromtheSchwarzschildvalue.
N4-N6 We attempt to account for this by assuming that our
A -3 backgroundspacetimecoordinatesareSchwarzschild,but
10
/ with g replaced by N2 , where N is the (time-
A tt − avg avg
∆ -4 N5-N6 dependent) averagevalue ofthe lapsefunction measured
10 on the extraction sphere. Under these assumptions, it
30c2 - 30c1
can be shown that the one-form
-5
10
N
avg ∗
0 1000 2000 3000 4000 dt dr (41)
1 2M /r −
ADM areal
t/M −
p
isnull,soweequatethisone-formwithduandthusdefine
FIG. 5: Convergence of waveforms with numerical resolution
and outer boundary location. Same as Fig. 4 except wave- t Navg
t = dt. (42)
formsaretime-shiftedandphase-shiftedsothatthemaximum s Z 1 2M /r
0 ADM areal
amplitude occurs at thesame time and phase. −
p
We show below (cf. Fig. 9) that choosing Eq. (42) in-
steadof t =t significantly increases the accuracyof our
s
C. Extrapolation of waveforms to infinity extrapolation procedure during merger and ringdown.
Having computed the retarded time at each extrac-
tion radius, we now consider the extracted waveformsas
Ournumericalsimulationscoveronlyafinitespacetime
functions of retardedtime u andextractionradius r ,
volume, as shown in Fig. 1, so it is necessary to extract areal
i.e. Ψ22(u,r ). At each value of u, we have the phase
our numerical waveforms at a finite distance from the 4 areal
and amplitude of Ψ22 at several extraction radii r .
source. However, gravitational-wave detectors measure 4 areal
Therefore at eachvalue of u, we fit phase and amplitude
waveforms as seen by an observer infinitely far from the
separately to a polynomial in 1/r :
source. Accordingly,afterextractingwaveformsatmulti- areal
plefiniteradii,weextrapolatethesewaveformstoinfinite n φ (u)
radius using a procedure similar to that described in [9]. φ(u,r )=φ (u)+ (k) , (43)
areal (0) rk
This extrapolation procedure is intended to remove not kX=1 areal
onlynear-fieldeffectsthatareabsentatinfinity, butalso n A (u)
(k)
gaugeeffects thatcanbe causedby the time-dependence r A(u,r)=A (u)+ . (44)
areal (0) rk
of the lapse function or the nonoptimal choice of tetrad Xk=1 areal
for computing Ψ .
4
Thephaseandamplitudeofthedesiredasymptoticwave-
TheextractionproceduredescribedinSec.IIIAyields
form are thus given by the leading-order term of the ap-
a set of waveforms Ψ22(t,r), with each waveform ex-
4 propriate polynomial, as a function of retarded time:
tracted at a different radius. To extrapolate to infinite
radiuswemustcomparewaveformsatdifferentradii,but
φ(u)=φ (u), (45)
(0)
thesewaveformsmustbeoffsetintimebythelight-travel
r A(u)=A (u). (46)
time between adjacent radii. To account for this time areal (0)
shift, for each extraction radius we compute Ψ22(u,r),
4 Figure 6 shows phase and amplitude differences be-
where u is the retarded time at that radius. Assuming
tween extrapolated waveforms that are computed us-
for simplicity that the background spacetime is nearly
ing different values of polynomial order n in Eqs. (43)
Schwarzschild, we compute the retarded time u using
and (44). For the extrapolation we use waveforms ex-
tractedatradii75M,85M,100M,110M,130M,140M,
∗
u ts r , (39) 150M, 160M, 170M, 180M, 190M, 200M, 210M, and
≡ −
225M. FromFig. 6 it is clearthat increasingn increases
where t is some approximation of Schwarzschild time, the accuracyof the extrapolationin smooth regions,but
s
10
n=2
n=4
0.02 0.01
) )
s s
n n=2 n=4 n
a a
i i
d 0 d 0
a a
r r
( (
n=1
φ φ
∆-0.02 n=3 ∆-0.01 n=1
n=3
3950 4000 4050
(t -r*)/M
0.05
s
A n=2 FIG. 7: Late-time phase convergence of extrapolation to in-
A/ 0 finity. Same as the top panel of Fig. 6, except zoomed to
∆ n=3 n=4 late times. The peak amplitude of the waveform occurs at
ts−r∗ =3954M.
-0.05
n=1
N5
0 1000 2000 3000 4000
)0.01 N6 n=2
(t -r*)/M s
n
s a
i
d 0
a
FIG. 6: Convergence of extrapolation to infinity for extrap- r
(
olation of order n. For each n, plotted is the extrapolated
φ
waveform from run 30c1/N6 using order n+1 minusthe ex- ∆-0.01
n=1
trapolated waveform using order n. The top panel shows
phase differences, the bottom panel shows amplitude differ-
ences. No shifting in time or phase has been done for this 3950 4000 4050
comparison. Increasing n increases accuracy in smooth re- (t -r*)/M
gions butalso amplifies noise. s
FIG.8: Effectofnumericalresolution onextrapolation toin-
finity. Thesolidcurvesareidenticaltothe“n=1”and“n=2”
also amplifies any noise present in the waveform. Our curvesfromFig.7. Thedottedcurvesarethesamequantities
preferredchoice,n=3,givesaphase errorof0.005radi- computed using thelower resolution run 30c1/N5.
ans and a relative amplitude error of 0.003 during most
of the inspiral, and a phase error of 0.01 radians and a
relative amplitude error of 0.01 in the ringdown. The effects. Such gauge effects might be reduced by improv-
junk radiation epoch ts r∗ < 1000M has moderately ing the gauge conditions in the numerical simulation or
larger errors than the ri−ngdow∼n. If we were to choose by adopting more sophisticated wave extraction and ex-
instead n = 4, we would gain higher accuracy in the trapolation algorithms that better compensate for dy-
smooth regions at the expense of increased noise in the namically varying gauge fields.
junkradiationepochandslightlylargererrorsduringthe Indeed, we have already made a first attempt at cor-
merger and ringdown. recting for a time-dependent lapse function by using t
s
Figure 7 is the same as the top panel of Fig. 6, ex- from Eq. (42) to compute the retarded time. Figure 9
cept zoomed to late times. Note that during merger illustrates the importance of this correction. Figures 7
andringdown,theextrapolationproceduredoesnotcon- and 9 differ only in the choice of ts used to compute the
verge with increasing extrapolation order n: the phase retarded time: In Fig. 7, ts is obtained from Eq. (42),
differences are slightly larger for larger n. This lack of and in Fig. 9, ts is simply the coordinate time t. Us-
convergencesuggeststhatthenonextrapolatednumerical ing the naive choice ts =t clearly results in much larger
waveform contains some small contamination that does phasedifferencesthatdivergewithincreasingnandgrow
not obey the fitting formulae, Eqs. (43) and (44). Fig- in time.
ure 8 shows the n=1 and n=2 convergence curves from In Fig. 10 we examine the difference between extrapo-
Fig. 7, but computed for two different numerical resolu- latedwaveformsandwaveformsthathavebeenextracted
tions, 30c1/N5 and 30c1/N6. The N5 and N6 lines are at a finite radius. We compare our preferred wave-
very close to each other in this figure, indicating that form, 30c1/N6 extrapolated to infinity using n=3, ver-
the lackofconvergencewith extrapolationordern isnot sus nonextrapolated waveforms and versus extrapolated
dominated by insufficient numerical resolution. We sus- waveformswithdifferentvaluesofn. Becausetheextrap-
pect that the main contribution is instead due to gauge olated and nonextrapolated waveforms differ by overall
