Table Of ContentA Simple method to set up low eccentricity initial data for moving puncture
simulations
Wolfgang Tichy and Pedro Marronetti
Department of Physics, Florida Atlantic University, Boca Raton, FL 33431, USA
We introduce two new eccentricity measures to analyze numerical simulations. Unlike earlier
definitionstheseeccentricitymeasuresdonotinvolveanyfreeparameterswhichmakesthemeasyto
use. Weshowhowrelativelyinexpensivegridsetupscanbeusedtoestimatetheeccentricityduring
the early inspiral phase. Furthermore, we compare standard puncture data and post-Newtonian
data in ADMTT gauge. We find that both use different coordinates. Thus low eccentricity initial
momentum parameters for a certain separation measured in ADMTT coordinates are hard to use
in puncturedata, because it is not known how theseparation in puncturecoordinates is related to
the separation in ADMTT coordinates. As a remedy we provide a simple approach which allows
us to iterate the momentum parameters until our numerical simulations result in acceptably low
1
eccentricities.
1
0
PACSnumbers: 04.25.dg,04.25.Nx,04.20.Ex,04.30.Db,
2
n
a I. INTRODUCTION whichconsidereccentricitiesoforderafewtimes10−3ac-
J ceptably small. There have been several previous works
4 that have considered eccentricities for puncture initial
Currently several gravitationalwave detectors such as data [20–22]. However, the most successful approach in
]
c LIGO [1, 2], Virgo [3, 4] or GEO [5] are already operat- termsofachievingloweccentricitieswasimplementedfor
q ing, while severalothers are in the planning or construc- excision type initial data [23–25]. The method discussed
- tion phase [6]. One of the most promising sources for inthisworkaimsatloweringtheeccentricityforthekind
r
g these detectors are the inspirals and mergers of binary of puncture initial data that is routinely used with the
[ black holes. In order to make predictions about the fi- moving puncture approach.
nal phase of such inspirals and mergers, fully non-linear Throughout we will use units where G = c = 1. The
2
v numerical simulations of the Einstein Equations are re- black hole masses are denoted by m1 and m2. We also
6 quired. To numerically evolve the Einstein equations, at introducethetotalmassM =m1+m2,thereducedmass
3 least two ingredients are necessary. First we need a spe- µ=m1m2/M and ν =µ/M.
9 cificformulationoftheevolutionequations. Andsecond, The paper is organized as follows. Sec. II introduces
2 to start such simulations initial data are needed. As the and compares several eccentricity measures. In Sec. III
.
0 first ingredient most groups nowadays use the BSSNOK we describe grid setups that can be used in numerical
1 formulation[7–9]oftheevolutionequations. Thisformu- simulations aimed at measuring the eccentricity. Sec. IV
0 lationisusuallyevolvedusingfinitedifferencingmethods, discusses a simple method to pick initial momentum pa-
1
but forsingleblackholes there havebeen someattempts rameters. ThisisfollowedbySec.Vwhichdescribeshow
:
v to use spectral methods [10–12]. For binary black holes to iterate these parameters to arrive at a reduced eccen-
i the BSSNOK system is usually used together with the tricity. We conclude with a discussion of our results in
X
moving puncture approach [13, 14]. This approach so Sec. VI.
r
a farworksonlywithfinite differencing methods sincecer-
tain evolved variables are not smooth inside the black
holes (at the punctures). Almost all simulations using II. DEFINING ECCENTRICITY FOR INSPIRAL
the BSSNOK formulationto date use standardpuncture ORBITS
data [15, 16] as initial data. These initial data are very
flexibleinthatthey containfreeparametersforthe posi- Real binary black hole orbits can never be circular.
tion,momentumandspinofeachblackholeandthusone Theyalwaysfollowspirals. Sowhenweareaimingforlow
can setup practically any kind of orbit. Note, however, eccentricityinitialdata,wereallywantdatathatresultin
that the emission of gravitational waves tends to circu- trajectoriesthatspiralinsmoothlywithoutoscillationsin
larize the orbits [17, 18]. Thus for realistic binary black the black hole separation. Of course this issue is further
holesystemsthathavebeeninspiralingalreadyforalong complicated by the fact that trajectories are coordinate
time,weexpectthetwoblackholestobeinquasi-circular dependent.
orbits around each other with a radius which shrinks on There are severalearliereccentricity definitions for in-
a timescale much larger than the orbital timescale. This spiral orbits in the literature [20–25]. All of them define
means that the initial data should be such that the or- eccentricity as a deviation from an underlying smooth,
bit has no or at least very small eccentricity. For our secular trendin some specific quantity that is associated
purposes here we follow the NRAR (Numerical Relativ- with the orbits. In [20] the frequency of the dominant
ity - Analytical Relativity) collaboration [19] guidelines l=m=2modeofthegravitationalwavesemittedisfit-
2
ted to a fourth order monotonic polynomial, and the de- binary. The idea is to determine the separation in a
viation of the frequency from this fit is used to compute more gauge invariant way from the amplitude of Ψ in-
4
the eccentricity. This approach works for non-spinning stead of using the gauge dependent coordinate separa-
binaries. Essentiallythesamemethodisalsousedin[21], tion. In [26] it is shown that for a non-precessing binary
but instead of the gravitationalwave frequency [21] uses in the quadrupole approximation the amplitude of the
the orbital frequency and also the coordinate separation l = m = 2 spin-weighted Spherical Harmonic mode is
to obtain two eccentricity measures. These same mea- given by
sures are also used in [22]. The approaches in [23, 24] fit
a linear function plus a sine function to the coordinate |C22|=32 π/5ν(Mω)8/3, (7)
separation and also the proper separation. The eccen-
where ω is the orbital anpgular velocity. Using Kepler’s
tricity can then be obtained from the amplitude of the
law we can define a separation
fitted sine function. In [25] the same fitting approach as
in [21] is used, but the fitted quantities are coordinate r =M1/3ω−2/3 =M[|C |/(32 π/5ν)]−1/4 (8)
separation,properseparationandalso orbitalfrequency. 22 22
Allthe approachesbasedinorbitalparametersshouldin which is directly related to the amplpitude of |C | of the
22
principle alsoworkforsystemswithspin. Also notethat
l = m = 2 mode of Ψ . Replacing the coordinate sepa-
4
alltheseeccentricitydefinitionsarechosensuchthatthey
ration in Eq. (1) by r we define
22
result in the correct value for Newtonian orbits.
Below we will introduce two new eccentricity defini- ∆r (t)−∆r (t)
22,max 22,min
e (t)= , (9)
tions and compare them to the earlier definition based 22 2r
22,av
on fitting the orbital frequency [21, 25].
The firsteccentricitydefinition is basedonthe coordi- which is an eccentricity definition that can be computed
nate separation of the two black holes. It is given by from Ψ4 alone. Note that this definition needs to be
extended for the case of precessing orbits, since in that
e (t)= ∆rmax(t)−∆rmin(t) (1) case |C22| will oscillate even for spherical orbits (i.e. or-
r 2r bits with r = const). The extension could be achieved
av
byinsteadusinga|C |thatiscomputedinacoordinate
22
where the average separation, and the maximum and
system where the z-axis points along the instantaneous
minimum deviation from a smoothed value rs are given orbital angular momentum.
by
Wehavealsotestedaneccentricitydefinitionbasedon
the coordinate angular velocity ω. Here the eccentricity
t+T/2
′ ′ is defined by [21, 25]
r = r(t)dt/T (2)
av
Zt−T/2
ω(t)−ω (t)
∆r (t) = max [r(t′)−r (t′,t)] (3) e (t)= fit , (10)
max s ω
t′∈[t−T/2,t+T/2] 2ωfit(t)
′ ′
∆r (t) = min [r(t)−r (t,t)]. (4)
min t′∈[t−T/2,t+T/2] s where ω(t) is simply the coordinate angularvelocityand
ω (t) is a polynomial fit of order 5 to ω(t) over a time
fit
Here the period T is defined using Kepler’s law interval corresponding to several complete orbits. The
hope is that the fit will smooth out oscillations so that
T =2π(r3/M)1/2. (5) e (t)∝ω(t)−ω (t)becomes ameasureofhowmuchω
ω fit
oscillates. The actual eccentricity is the maximum mag-
Notice thatthe actualorbitalperiodmaybe slightlydif-
nitude of e (t).
ω
ferent, but this estimate suffices to get an approximate
Inordertocomparethesethreeeccentricitydefinitions
eccentricity measure. The smoothed value r (t′,t) is ob-
s we now present the eccentricities from an actual numer-
tained from
ical simulation. As we can see in Fig. 1 all three eccen-
r(t+T/2)−r(t−T/2) tricity definitions agree well for 300M ≤ t ≤ 3000M.
′ ′
rs(t,t)=r(t)+ T (t −t), (6) They yield a value of about 0.02. Notice that er and e22
are direct eccentricity measures, while in the case of e
ω
butdifferentsmoothingsarepossible(e.g.,byperforming theeccentricitycorrespondstothemaximumvalueofthe
a least-squares fit of e.g. a linear or quadratic function magnitude of e . The eccentricity definition e which is
ω r
to r(t) in the interval [t+T/2,t−T/2]). Essentially the calculated from the separation has certain problems for
definition in Eq. (1) measures how much the coordinate t< T/2∼135M. These occur because we have no data
separation oscillates over the time T. For Newtonian for t < 0, which is needed in the average over one com-
orbits it coincides with the usual eccentricity definition plete period (centered around t) and also for the fitting
forellipticorbits. Fororbitswhoseradiusshrinkslinearly in Eq. (6). These same problems also affect e . In e ,
22 22
in time (without any oscillations) e (t) is zero. however, they are exacerbated by the initial junk radi-
r
Another similar eccentricity measure can be obtained ation that dominates |C | until about 150M (see mid-
22
using the gravitational wave signal of the inspiraling dle panel). We see that e does not completely settle
22
3
12 NOKformalism[7–9]inthevariationknownasthe“mov-
10
M 8 ing punctures” method [13, 14]. The particulars of our
r/ 46 BSSNOK implementation can be found in [28, 29]. For
2
completeness we note that lapse and shift evolve accord-
0
0.0010 ing to
0.0008
|C|2200..00000046 (∂t−βi∂i)α = −2αK,
0.0002 3
0.0000 (∂ −βk∂ )βi = Bi,
t k
4
entricity 000...000024 eeer22 (∂t−βk∂k)Bi = (∂t−βk∂k)Γ˜i−ηBi. (11)
c ω
ec-0.02 Theshiftdriverparameterissetη =2/M inallourruns.
0 500 1000 1500 2000 2500 3000 TheBAMcodeisbasedonamethodoflinesapproach
t/M
using sixth order finite differencing in space and explicit
fourth order Runge-Kutta time stepping. The time step
FIG. 1: This plot shows results from a numerical run (with
size is chosen such that the Courant factor is either 0.25
parameters from row 6 of table I. The upper panel shows
the coordinate separation r. The middle panel depicts the or0.5. Forefficiency,Berger-Oligertypemeshrefinement
magnitude of the l=m=2 mode |C22| of Ψ4 extracted at a is used [30]. The numerical domain is represented by a
separationof70M fromthecenterofmass. Inthelowerpanel hierarchy of nested Cartesian boxes. The hierarchy con-
we plot the three different eccentricities obtained with the sistsofL+1levelsofrefinement,indexedbyl =0,...,L.
definitions given in Eqs. (1), (9) and (10). Note that er and ArefinementlevelconsistsofoneortwoCartesianboxes
er directly measure eccentricity, while for eω the eccentricity with a constant grid-spacing h = h /2l on level l. We
l 0
corresponds to the maximum values of |eω(t)|. haveusedhereL=10to11forthenumberofrefinement
levels,withthelevels0through5eachconsistingofasin-
gle fixed box centeredon the origin(the center of mass).
down until about 700M. Notice also that during stan-
On each of the finer levels 6 through L, we initially use
dard moving puncture evolutions the coordinates adjust
two sets of moving boxes centered on each black hole.
quite rapidly initially. This means that any eccentricity
When the black holes get close enough that two of these
definition that is based on the coordinate separation or
boxes start touching, they are replaced by a single box.
the coordinate angular velocity will not be completely
The position of each hole is tracked by integrating the
reliableduringthe first100M orso. Thusitcomesasno
shift vector. We use this same set up but with different
surprisethattheeccentricitydefinitione hasadifferent
ω resolutionsdependingonthe purposeofeachsimulation.
frequency and amplitude in the beginning.
For an accurate simulation of the inspiral and merger
The curves for e and e in Fig. 1 oscillate even at
r 22 oftwonon-spinningequalmassblackholeswe mightuse
later times. This oscillation is due to the fact that the
L=10 with a resolution h =M/96 on the finest level
10
period T which we get from Kepler’s law is not exactly
using 144 points on the fixed levels and 72 points on the
equaltotheactualorbitalperiod. Themagnitudeofthis
moving levels. The notation we use to describe this grid
oscillationcanbe usedas anerrorestimate ofour eccen-
setup for this simulation is:
tricitymeasures. Foraconservativeeccentricityestimate
we can use the maximum values of e and e .
r 22 [5×72,6×144][M/h =96,OB =768M][C =0.25]
10
Let us point out that e and e are easier to com-
r 22 (12)
pute than e . The latter depends on a polynomial fit
ω which indicates that we have 5 moving levels with 72
to the measured orbital angular velocity. The problem
points in each box and 6 fixed levels with 144 points
is that this fit has to be done over a certain time inter-
each. TheresolutionisgivenbyM/h =96onthefinest
10
val, that must terminate well before the merger. Thus
level, which results in an outer boundary at 768M. The
it requires a certain amount of fine tuning and human
Courantfactorhere is chosento be 0.25,whichimplies a
intervention. On the other hand e and e can be com-
r 22 time step of dt =0.25h =M/384 on the finest level.
10 10
puted at any time in a very simple way. Recall however,
Ifthe blackholeshavespinsand/orunequalmasseseven
that e is more sensitive to the initial junk radiation.
22 more resolution is needed. For example for a mass ratio
Thus for short runs used to probe the eccentricity of a
of 3 and a dimensionless spin magnitude of 0.6 on the
configuration we usually just use e .
r larger hole we might use a setup described by:
[6×72,6×144][M/h =192,OB =768M][C =0.25]
11
III. NUMERICAL EVOLUTIONS (13)
this setup has twice the resolution on the finest level, so
Thenumericalresultsdiscussedinthispaperhavebeen that the resolution in terms of the individual masses is
obtained with the BAM code [27–29]. As already men- nowgivenbym /h =48forthesmallerandm /h =
1 11 2 11
tioned,thegravitationalfieldsareevolvedusingtheBSS- 144 for the larger hole. Such runs are quite expensive.
4
0.0025 2 3 δ pknl
[6x72, 6x144][M/h =192, OB=768M][C=0.25] Kij = ψ−10 p(inj)− kl A A(δij −ni nj )
11 r2 A A 2 A A
0.0020 [5x48, 6x96][M/h10=85.3, OB=576M][C=0.5] AX=1n A (cid:20) (cid:21)
6
+ ǫ Sl nmδk(inj) . (16)
r2 klm A A A
A
0.0015 h io
er Herepi andSi arethemomentumandspinparameters
A A
0.0010 of black hole A, while rA and niA denote the distance
andnormalvectormeasuredfromholeA. Theconformal
factor is
0.0005
m m
ψ =1+ b1 + b2 +u, (17)
2r 2r
0.0000 1 2
0 200 400 600 800 1000 1200 1400
t/M
where m and m denote the black hole bare mass pa-
b1 b2
rameters. The scalaru is computed by numerically solv-
FIG. 2: This plot shows the eccentricity er for two different ing the Hamiltonian constraint. These initial data are
grid setups. Both start with the parameters from row 7 of
very flexible since the free parameters for the position,
tableI. Thesolidlineshowserforahighresolutiongrid,while
thebrokenlineshowsthemeasureder ifweusealessaccurate momentum and spin of each black hole can be chosen
andcomputerintensivegridsetup. Bothyieldsimilar orbital freely. Thus one can setup practically any kind of orbit.
eccentricities. Note, however, that our goal is to set up data for black
holes that are in quasi-circular orbit. This means that
we need to choose our momentum parameters such that
Until mergerthey take about one month ona Cray XT5
the eccentricity is as small as possible. Since we start
supercomputer like NICS Kraken if we use 96 cores.
our evolutions in a frame where the center of mass is at
However, if our objective is to simply measure the or-
rest, both black holes have momenta that are equal in
bital eccentricity of a binary (characterized by certain
magnitude but opposite in direction. Thus we have to
initial parameters), it is sufficient to evolve for only a
choose only two parameters: the tangential and radial
few orbits. We have found that such an evolution does
component of the momentum of one of the black holes.
not need great accuracy. So if the goal is to simply de-
To complete the definition of the initial data, we also
terminetheinitialeccentricityweusethefollowingsetup
need to specify initial values for the lapse α and shift
[5×48,6×96][M/h =85.3,OB =576M][C =0.5] vector βi. At time t=0 we use
10
(14)
With this grid setup it takes only two days to evolve up α = ψ−2,
t=800M if we use 48coresonNICS CrayXT5 Kraken. βi = 0. (18)
InFig. 2weseethattheeccentricitye fromthischeaper
r
run agrees quite well with the e from a more expensive
r
run. As explained above no eccentricity estimate is ac-
A. Finding the momentum parameters
curateatthe startofarun. We needto waituntilabout
500M for e to settle down to a regular oscillation. As
r
mentionedbeforeweusethemaximaofe asaconserva- There have been various attempts to guess appropri-
r
tive eccentricity estimate. Thus the eccentricities of the ate momentum parameters. Some have used the quasi-
cheap and expensive runs are 0.0013 and 0.0010. Hence equilibrium approach [31–33], but most are based on
weneedtoevolveuntilatleastabout800M togetareli- post-Newtonian (PN) approximations (see e.g. [29, 34,
able estimate. From the oscillations in the curves we see 35]). The latter have been taken to their extreme in [21]
that the errors in both eccentricity estimates are about and [22] who integrate PN equations in the ADMTT
0.0005. gauge [36] in the hope that the thus obtained momen-
tum parameters will lead to orbits with less eccentricity.
However,asmentionedin[22],itis notcertainthatnon-
IV. CHOOSING PUNCTURE PARAMETERS eccentric PN parameters in ADMTT gauge should pro-
duce non-eccentric orbits in full General Relativity. We
want to point out here, that standard puncture data are
In order to start our simulations we need initial data
inconsistentwithPNtheorybeyond(v/c)3 [46] [34,37–
forbinaryblackholeswitharbitraryspinsandmassesat
39], even in ADMTT gauge! The reasons for this incon-
somegiveninitialseparationr. Since wewill employthe
sistency are as follows. First, the 3-metric of puncture
movingpuncturesapproachinourevolutions,wewilluse
data is always conformally flat (see Eq.(15)), while the
standard puncture initial data [15]. Thus the 3-metric
PN 3-metric contains deviations from conformal flatness
and extrinsic curvature are given by
at order (v/c)4 [37, 40]. This is true for both harmonic
g = ψ4δ (15) and ADMTT gauge. Second, the conformal factor in
ij ij
5
ADMTT gauge is given by [37] p =µ|r˙|. (21)
r
E1 E2 In all cases discussed in this paper we put the two punc-
ψ =1+ + , (19)
PN
2r 2r tures on the y-axis at y = m r/M and y = −m r/M.
1 2 1 2 2 1
The two initial black hole momenta are then p =
where E is the energy of each particle used to model 1
A
(−p ,−p ,0) and p =(p ,p ,0).
the black holes. If we compare Eqs. (19) and (17) we t r 2 t r
see that the conformal factor in ADMTT gauge is not
identical to the conformal factor used in puncture initial
B. Determining the bare mass parameters
data. One difference is that the ADMTT ψ contains
PN
the particle energies E while the puncture ψ contains
A
the bare masses mbA. These two only agree for infinite The momenta p1, p2 and the spins S1, S2 for a co-
separation. Furthermore, the puncture ψ contains the ordinate separation r directly enter the Bowen-York ex-
additional piece u that is obtained by numerically solv- trinsic curvature of standard puncture data. Note, how-
ing the Hamiltonian constraint, while PN data violates ever, that the bare mass parameters mb1 and mb2 which
theHamiltonianconstraint. Onemightliketoarguethat appear in the construction of standard puncture data
the PN data should agree with puncture data up to the are not equal to the individual black hole masses. As
higherordertermsthatareneglectedinthe PNapproxi- in [29, 34, 42, 43] we obtain the bare masses from the
mation. Thisishowevernottrueclosetotheblackholes, condition that the ADM masses
because PN theory is not valid in regions of very strong mADM =m (1+u )+ mb1mb2, (22)
gravity. Thus near the black holes the ADMTT metric A bA A 2r
differsfromthepuncturemetricbytermsthatareoflow
measured at the punctures [44] should be equal to m
PN order. 1
and m , where u is the value of u at puncture A. As
From the above explanation it is clear that the coor- 2 A
in[16,33,44]weassumethattheADMmassesmeasured
dinate separation in ADMTT gauge does not have the
ateachpunctureareagoodapproximationfortheinitial
same physical meaning as the coordinate separation in
individual black hole masses. Numerically this condition
puncture initial data. Thus if by some method we arrive
is implemented as a root finder in the initial data solver
atthemomentumparametersneededataparticularsep-
that picks m and m such that the ADM masses at
arationr(inADMTTcoordinates),weshouldnotsimply b1 b2
the punctures are equal to m and m .
usethesesamemomentumparametersatthesamepunc- 1 2
ture separation.
So if (for some configurations) the more complicated
V. REDUCING THE ECCENTRICITY
approachesin[21]or[22]leadtolesseccentricorbitsthan
e.g. the simpler approaches in [29, 34, 35] this should
Now that we know how to generate initial data on
be considered a coincidence, since it is not due to the
approximately circular orbits for arbitrary spins, masses
inclusion of higher order terms or due to the integration
and separations,it is time to present some numerical re-
of post-Newtonian equations of motion. In fact, we find
sults to demonstrate what eccentricities we get and how
that (as anticipated by [22]) for many configurationsthe
we can reduce them. For our purposes here, we consider
approaches in [21] or [22] do not lead to less eccentricity
the eccentricity small enough if it is 0.003 or less, which
than the ones in [29, 34, 35].
agrees with the NRAR target of an eccentricity of a few
For these reasons we take a different approach here.
times 10−3. The first example we consider is an equal
We will take a simple PN formula to obtain a reason-
mass binary without spin. The momentum parameters
able guess for both the tangentialand radialmomentum
are picked according to Eqs. (20) and (21). The values
components p and p . Then we numerically evolve the
t r
of all the initial parameters are given in the first line of
resulting initial data for a short time (using the efficient
table I. As we can see, the eccentricity is about 0.005 in
grid setup in Eq. (14)) to see how eccentric they really
this case. In order to test how it varies with p we have
are. Afterwardweadjustp toreducetheeccentricity. In t
t
also performed a run where we have increased p by
thiswaywecanobtainloweccentricityorbitsforanycon- t
figuration. In order to come up with a reasonable guess M M 3 M 4
for pt and pr we use 2PN accurate expressions of Kid- ∆pt =µ 11.29 −92.37 . (23)
r r r
der[41]inharmonicgauge. Specifically,wefreely choose r " (cid:18) (cid:19) (cid:18) (cid:19) #
the two masses m , m , the six spin components of S
and S and a sepa1ratio2n r. Next, choosing Lˆ in the z1- The results of this increase are shown in the second line
2 N of table I and yield an eccentricity that is reduced by
direction, we use Eqs. (2.8) and (4.7) of [41] to compute
more than a factor of 3. The expression for Eq. (23)
the total orbital angular momentum L, and Eq. (4.12)
comesfromthefittingoftwoofourolderequalmassruns
of [41] to compute r˙. We rotate L, S and S so that L
1 2 that resulted in low eccentricity. Thus non-zero spins or
pointsinthez-direction. Thenweobtainthemomentum
unequal masses are not taken into account by this fit.
in the xy-plane as
Therefore Eq. (23) certainly does not present the opti-
p =|L|/r (20) mum momentum correction. In fact we have found that
t
6
row Mr mMb1 mMb2 mM1 mM2 |mS121|,θ1,φ1 |mS222|,θ2,φ2 10M2pt 10M4pr 103er
1 11.9718 0.488255 0.488255 1/2 1/2 0 0 8.5018 3.88 5.0
2 11.9718 0.488249 0.488249 1/2 1/2 0 0 8.5168 3.88 1.5
3 11.9694 0.404597 0.404687 1/2 1/2 0.6,60,0 0.6,120,90 8.5013 3.87 4.0
4 11.9694 0.404592 0.404681 1/2 1/2 0.6,60,0 0.6,120,90 8.5163 3.87 2.0
5 11.9694 0.404588 0.404677 1/2 1/2 0.6,60,0 0.6,120,90 8.5280 3.87 1.5
6 11.4546 0.2239 0.6145 1/4 3/4 0.4,0,0 0.6,0,0 6.4416 2.21 20
7 10.7271 0.223341 0.614005 1/4 3/4 0.4,0,0 0.6,0,0 6.6396 2.50 1.0
8 12.0815 0.224174 0.614794 1/4 3/4 0.4,180,0 0.6,180,0 6.5900 2.50 1.0
9 11.7173 0.223891 0.614572 1/4 3/4 0.4,180,0 0.6,180,0 6.7497 4.08 5.0
10 10.9876 0.223317 0.614121 1/4 3/4 0.4,0,0 0.6,180,0 6.9756 3.19 3.0
11 11.5010 0.300017 0.543886 1/3 2/3 0.4,60,0 0.6,60,0 7.6504 3.18 2.0
12 11.5033 0.300020 0.543926 1/3 2/3 0.4,60,0 0.6,60,90 7.6509 3.19 1.8
13 11.6182 0.300023 0.543980 1/3 2/3 0.4,60,0 0.6,120,90 7.7896 3.51 2.1
14 11.5502 0.310095 0.656229 1/3 2/3 0.3,0,0 0 7.7196 3.32 1.3
TABLE I: Initial data parameters. The black holes have coordinate separation r. We give both bare masses mb1, mb2 as well
as physical masses m1, m2. The puncturesare located on they-axis at y1 =m2r/M and y2 =−m1r/M. The spins are given
in termsoftheirmagnitudesandtheusualpolar angles ofspherical coordinates measured indegrees. Thelinear momentaare
(∓pt,∓pr,0). The last column shows theresulting eccentricity.
adding ∆pt to pt does notalwaysreduce the eccentricity 12.0
sincethe PNestimateofEq.(20)forp is sometimestoo
t
large and sometimes too small for generic orbits. The
expression in Eq. (23) simply gives us a rough estimate 11.5
by how much we might have to raise or lowerp in order
t
to reduce the eccentricity.
M11.0
Soin orderto reallyreduce the eccentricityweusually r/
start one run on the coarse grid described by Eq. (14).
This run normally uses p given by Eq. (20). We then 10.5
t
look at the coordinate separation r(t) for this run. Usu-
r(0)=12.08M, p=.0659M, p=.00025M
allyonecantellwhethertheinitialtangentialmomentum t r
p was too large or too small. We then start a new run 10.0 r(0)=11.72M, pt=.0675M, pr=.00041M
t
whereweeither increaseordecreasep by∆p . This run 0 200 400 600 800 1000 1200
t t
t/M
then gives a new r(t) and e . From the two resulting
r
simulationsone canthen extrapolateto zeroeccentricity
FIG.3: Thisplot shows r(t)for thetwonumericalrunswith
to obtain a more refined tangential momentum parame-
parametersfrom rows 8(solid line)and9(brokenline) ofta-
ter. This procedure is illustrated in rows 3,4, and 5 of
bleI. Theparametersforthesolidlinewerepickedaccording
table I. In row 3 we have used the p =8.5013×10−2M
t to the method explained in this paper while the broken line
from Eq. (20) and obtained an eccentricity of 0.004. In is theresult of setting theparameters according to themuch
row 4 we have increased this pt by ∆pt which leads to a more complicated method in [22]. As we can see the eccen-
decrease of the eccentricity to 0.002. A further increase tricity difference of about a factor of 5 is quite noticeable in
to p = 8.5280× 10−2M (as in row 5) then yields an ther(t) curves.
t
eccentricity of 0.0015.
Row6oftableIshowstheinitialdatathatwereusedto
produceFig.1. Inordertoproduceaneccentricitywhich
from the method introduced in [22] for the same mass
is clearly visible in the r(t) curve we have used a p that
t ratio and spin configuration. We can see that it is about
isdeliberatelychosenmuchlargerthanwhatispredicted
5 times larger, which supports our argument that inte-
by Eq. (20). In row 7 we show the same configuration
gratingPNequationsofmotiondoesnotnecessarilylead
(atsomewhatcloserseparation),butthistimewechoose
to better results. The results from these two runs are
p according to Eq. (20), we see that in this case the
t also shown in Fig. 3. The solid line corresponds to row
eccentricity is already small enough to satisfy e.g. the
8 (i.e. our approach) and the broken line is generated
NRAR guidelines.
using the approach in [22]. As we can see our method
Inrow8wehaveusedourmethodtoproducemomen- does give a noticeably less eccentric r(t) curve. Notice
tumparametersforasimilarconfiguration,butthistime that the initial dip in the coordinate separation r(t), is
both spins point in the negative z-direction. Again we due tothe aforementionedinitialcoordinateadjustment.
can reach the NRAR target of an eccentricity of order Itdoesnotmeanthetwoholesreallyplungetowardeach
0.001. In row 9 we also show the eccentricity resulting other initially.
7
ThelastfiverowsoftableIgiveafewmoreexamplesof gular velocity to a polynomial of some low degree. This
parameters and resulting eccentricities. They show that degreeisessentiallyanotherfreeparameterandisusually
we can obtain low eccentricities for generic mass ratios chosen to be 4 or 5. Note, however, that all eccentricity
and spin orientations. definitions are ambiguous to a certain extend since the
Notice that our eccentricity reduction method is simi- entire concept of eccentricity is only rigorously defined
lar in spirit to the method described in [23, 24], in that for periodic orbits. All eccentricity definitions for inspi-
we also use short numerical runs to adjust certain pa- ralorbitsdependonhowwesplitafunctionoftime(like
rameters. There are, however, important differences. the separation) into a smooth and an oscillating piece.
In [23, 24] excision type initial data, not punctures are Thus we do not claim that our definitions e and e
r 22
used. These data are constructed using an extension are less ambiguous than earlier ones. For example our
to the conformal thin sandwich formalism [45]. Empir- definitions depend on the period T which we compute
ically it turns out that the tangential momentum that from Kepler’s law, but other choices are possible (e.g.
is achieved with the conformal thin sandwich method is an estimate of the actual orbital period of the numerical
quite close to what is needed for low eccentricity data. simulation at each time). Since our definitions e and
r
The main reason for eccentricity is the absence of a ra- e do not require us to choose extra parameterssuch as
22
dial momentum in the original conformal thin sandwich fitting intervals they are easier to use.
method. In[23]amethodisdevelopedthataddsanarbi- We also show that certain low resolution grid setups
trary radial velocity parameter to the initial data. This (whicharerelativelyinexpensive)canbeusedtoestimate
radial velocity parameter together with the tangential the initial eccentricity using our measure e . This gives
r
momentum are then adjusted to reduce the eccentricity. us a relatively efficient way to measure the eccentricity
Themethodintroducedin[23,24]iscapableofproducing of any initial data.
eccentricities of order of only 10−5. For standard punc- Furthermore, we have explained why the coordinates
ture data the method in[23, 24] cannotbe used directly. of standard puncture data are not the same as in PN
One problem is that standard moving puncture simula- calculations in ADMTT or harmonic gauge. Thus even
tions start with the lapse and shift given in Eqs. (18). if one arrives at a highly accurate estimate of the mo-
Thusthecoordinatesusedarenotwelladaptedtoquasi- mentum parametersneeded for low eccentricity orbits in
equilibrium. Hence oscillations in the hole separation r ADMTT gauge, there is no easy way to incorporate this
are not due to realeccentricity alone but also due to the knowledgeintosimulationsstartingfromstandardpunc-
fact that the coordinates are still evolving as well. This turedata. Indeedwefindthatusingaccurateparameters
problem is quite visible in Fig. 3. The broken line shows in ADMTT gauge [21, 22] in standard puncture initial
oscillationsthat have more thanone frequency. Thus we data can lead to relatively large eccentricities.
cannot fit the curve very well with a straight line plus
We provide a simpler approach which starts from mo-
a single sine function as in [23, 24]. To summarize, re-
mentumparametersusingtherelativelyshortexpressions
ducing the eccentricity for puncture data is harder than
given by Kidder [41]. After measuring the resulting ec-
for the extended conformal thin sandwich data in [23].
centricitye inashortinexpensivenumericalsimulation,
r
Therefore, our results in both final eccentricities and re-
these parameters can then be refined by changing p by
t
duction of eccentricity per iteration are not as good as
an amount of order ∆p (see Eq. (23)). In this way we
t
in [23, 24].
canalwaysarriveateccentricitiesthatarelowenoughfor
In the cases we have studied so far, we have found
the purposes of the NRAR collaboration [19]. In order
it unnecessary to adjust p (away from the PN value in
r to achieve this we usually do not need more than three
Eq.(21))inorderto reachaneccentricityoforder0.001.
numerical runs.
It is much more important to choose an appropriate p .
t
However, we expect that adjusting p will be necessary
r
to reach even lower eccentricities.
Acknowledgments
VI. DISCUSSION It is a pleasure to thank Doreen Mu¨ller for use-
ful discussions about the approach in [22]. This work
Wehaveintroducedthetwoneweccentricitymeasures was supported by NSF grant PHY-0855315. Compu-
e and e . Both are easy to compute since their calcu- tational resources were provided by the Ranger cluster
r 22
lationdoes notinvolveanyfreeparametersunlike earlier at the Texas Advanced Computing Center (allocation
definitions akin to e (t). To compute eccentricity mea- TG-PHY090095)andtheKrakencluster(allocationTG-
ω
sures such as e (t) one needs to specify a time interval PHY100051)attheNationalInstituteforComputational
ω
duringtheinspiralphaseoverwhichwefittheorbitalan- Sciences.
[1] B. Abbott et al. (The LIGO Scientific Collaboration) [2] LIGO, http://www.ligo.caltech.edu.
(2007), arXiv:0711.3041 [gr-qc], 0711.3041.
8
[3] F.Acerneseetal.,J.Opt.A:PureAppl.Opt.10,064009 D75, 124018 (2007), gr-qc/0610122.
(2008). [27] B.Bru¨gmann,W.Tichy,andN.Jansen,Phys.Rev.Lett.
[4] VIRGO, http://www.virgo.infn.it. 92, 211101 (2004), gr-qc/0312112.
[5] GEO600, http://www.geo600.org. [28] B. Bru¨gmann, J. Gonz´alez, M. Hannam, S. Husa,
[6] B. Schutz,Class. Quantum Grav. 16, A131 (1999). U. Sperhake, and W. Tichy, Phys. Rev. D77, 024027
[7] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. (2008), gr-qc/0610128.
Phys.Suppl.90, 1 (1987). [29] P.Marronetti,W.Tichy,B.Bru¨gmann,J.Gonz´alez,and
[8] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 U.Sperhake,Phys.Rev.D77,064010(2008),0709.2160.
(1995). [30] M. J. Berger and J. Oliger, Journal of Computational
[9] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D59, Physics 53, 484 (1984).
024007 (1998), gr-qc/9810065. [31] T. W. Baumgarte, Phys. Rev. D 62, 024018 (2000), gr-
[10] W. Tichy, Phys. Rev. D74, 084005 (2006), gr- qc/0004050.
qc/0609087. [32] W. Tichy, B. Bru¨gmann, and P. Laguna, Phys. Rev.
[11] W. Tichy, Class. Quant. Grav. 26, 175018 (2009), D68, 064008 (2003), gr-qc/0306020.
0908.0620. [33] W. Tichy and B. Bru¨gmann, Phys. Rev. D 69, 024006
[12] W. Tichy,Phys.Rev. D80, 104034 (2009), 0911.0973. (2004), gr-qc/0307027.
[13] M. Campanelli, C. O. Lousto, P. Marronetti, and [34] P. Marronetti, W. Tichy, B. Bru¨gmann, J. Gonz´alez,
Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), gr- M. Hannam, S. Husa, and U. Sperhake, Class. Quant.
qc/0511048. Grav. 24, S43 (2007), gr-qc/0701123.
[14] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and [35] B. Bru¨gmann, J. A. Gonz´alez, M. Hannam, S. Husa,
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), gr- U. Sperhake, and W. Tichy, Phys. Rev. D77, 024027
qc/0511103. (2008), gr-qc/0610128.
[15] S. Brandt and B. Bru¨gmann, Phys. Rev. Lett. 78, 3606 [36] G. Sch¨afer, Annalsof Physics 161, 81 (1985).
(1997), gr-qc/9703066. [37] W. Tichy, B. Bru¨gmann, M. Campanelli, and P. Diener,
[16] M. Ansorg, B. Bru¨gmann, and W. Tichy, Phys. Rev. Phys. Rev.D67, 064008 (2003), gr-qc/0207011.
D70, 064011 (2004), gr-qc/0404056. [38] N.YunesandW.Tichy,Phys.Rev.D74,064013(2006),
[17] P.C.PetersandJ.Mathews,Phys.Rev.131,435(1963). gr-qc/0601046, gr-qc/0601046.
[18] P.C. Peters, Phys. Rev.136, B1224 (1964). [39] N. Yunes, W. Tichy, B. J. Owen, and B. Bru¨gmann,
[19] NRAR, https://www.ninja- Phys. Rev.D74, 104011 (2006), gr-qc/0503011.
project.org/doku.php?id=nrar:home. [40] N.K.Johnson-McDaniel, N.Yunes,W.Tichy,andB.J.
[20] J.G.Baker, J. R.vanMeter, S.T.McWilliams, J.Cen- Owen, Phys. Rev.D80, 124039 (2009), 0907.0891.
trella, and B. J. Kelly, Phys. Rev. Lett. 99, 181101 [41] L. E. Kidder, Phys. Rev. D52, 821 (1995), gr-
(2007), gr-qc/0612024. qc/9506022.
[21] S. Husa, M. Hannam, J. A. Gonzalez, U. Sperhake, [42] W. Tichy and P. Marronetti, Phys. Rev. D76, 061502
and B. Bruegmann, Phys. Rev. D77, 044037 (2008), (2007), gr-qc/0703075.
0706.0904. [43] W. Tichy and P. Marronetti, Phys. Rev. D78, 081501
[22] B. Walther, B. Bruegmann, and D. Mueller, Phys. Rev. (2008), 0807.2985.
D79, 124040 (2009), 0901.0993. [44] W. Tichy, B. Bru¨gmann, and P. Laguna, Phys. Rev.
[23] H.P. Pfeiffer et al., Class. Quant.Grav.24, S59 (2007), D68, 064008 (2003), gr-qc/0306020.
gr-qc/0702106. [45] J. W. York,Phys.Rev.Lett. 82, 1350 (1999).
[24] M.Boyle, D.A.Brown,L.E.Kidder,A.H.Mrou,H.P. [46] ThestandardPNjargonforatermoforder(v/c)ninthe
Pfeiffer, M. A.Scheel, G. B. Cook, and S. A. Teukolsky, equations ofmotion is nPN term.Noticethat according
2
Phys.Rev.D76, 124038 (2007), 0710.0158. to the virial theorem v2 ∼M/r, so that nPN terms are
2
[25] A. H. Mroue, H. P. Pfeiffer, L. E. Kidder, and S. A. n
also O(M/r)2.
Teukolsky (2010), 1004.4697.
[26] A.Buonanno, G. B. Cook, and F. Pretorius, Phys. Rev.
