Reducing phase error in long numerical binary black hole evolutions with sixth order finite differencing
Sascha Husa, Jose A. Gonzalez, Mark Hannam, Bernd Bruegmann, Ulrich Sperhake
aa r X i v : . [ g r- q c ] J un Reducing phase error in long numerical binary blackhole evolutions with sixth order finite differencing
Sascha Husa, Jos´e A. Gonz´alez, Mark Hannam, BerndBr¨ugmann, Ulrich Sperhake
Friedrich Schiller University Jena, Max-Wien-Platz 1, 07743 Jena, GermanyE-mail: [email protected]
Abstract.
We describe a modification of a fourth-order accurate “movingpuncture” evolution code, where by replacing spatial fourth-order accuratedifferencing operators in the bulk of the grid by a specific choice of sixth-orderaccurate stencils we gain significant improvements in accuracy. We illustrate theperformance of the modified algorithm with an equal-mass simulation coveringnine orbits.PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf
1. Introduction
In the last two years the numerical solution of the general relativistic two-body problemhas made a giant leap forward with a series of breakthroughs in 2005 [1, 2, 3]. Morethan forty years after Hahn and Lindquist started the numerical investigation ofcolliding black holes [4], the field has now passed several crucial milestones towardsimulating general inspiral situations, such as simulations of unequal-mass binaries andcalculations of the gravitational recoil effect and the evolution of black-hole binarieswith spin [5, 6, 7, 8, 9, 10]. The latter have recently lead to the spectacular findingthat extremely large recoils are possible for spinning black holes [11, 12, 13].In order to fulfill numerical relativity’s promise of providing useful informationto the gravitational wave data analysis community, it is desirable to perform longnumerical inspiral evolutions that allow us to cleanly match fully general relativisticand post-Newtonian waveforms with error bars, and thus produce “complete”waveforms, which contain large numbers of gravitational-wave cycles from the inspiralphase, as well as simulating the merger and ringdown phases. Such simulations will benecessary for a sufficiently dense sample of the black-hole binary configuration space.Comparisons with post-Newtonian results have already started, and severalgroups have published results showing good agreement of various aspects of non-spinning simulations with post-Newtonian predictions (see e.g. [14, 15, 16, 17, 18, 19,20]). Precise error estimates and detailed coverage of the unequal-mass and spinningcases are however missing. One serious technical problem is that performing such longevolutions with good accuracy in the phase is still computationally expensive — atleast for standard “moving puncture” finite-difference codes [21, 22, 23, 16, 13, 10].In order to overcome phase inaccuracies in long evolutions, an alternative route is inary black hole evolutions with sixth order finite differencing ” grid configuration) and 27 (“76 ” grid) CPU hoursper M [26] (we will quote time and length units of the total black-hole mass M , see[16]). The highest resolution run presented here performs at 5.5 CPU hours per M on an Intel Woodcrest 2.66 GHz dual core processor. Both codes show satisfactoryaccuracy for long evolutions, and we certainly expect both codes to undergo furtheroptimizations. Long evolutions that show fourth order convergence for most of thesimulation have also been presented by the Goddard group [27], using impressivelylow spatial resolution, but no details are given on code timing.We have previously reported accurate evolutions for approximately two orbits(initial separation of D = 6 . M ) in [16], with an error in the merger time of 0.2 % ata computational cost of 505 CPU hours (1.44 CPU hours/M), and we have reported amerger time error of 0.5% for D = 8 M simulations [28]. At larger initial separationsthe number of orbits is a steep function of the separation, and phase accuracy rapidlydecreases. Our fourth order code would thus require resolutions which we find hardto tolerate for performing large parameter studies in the style of [6]. In the contextof finite differencing it is natural to consider higher order methods. For example, italready turned out to be important to move from second-order finite differencing tofourth-order finite differencing as the feasible evolution time for puncture evolutionsincreased.However, using higher order finite differencing for the types of codes used tosimulate black-hole binary inspiral is not entirely trivial: for the moving-puncturemethod, which is currently employed in the majority of the current codes, thecontinuum equations become singular at the location of the “puncture”, and onemight worry about the robustness of finite-difference schemes. Furthermore, currentmesh refinement algorithms in the field are based on the use of buffer zones, whosenumber depends on the stencil width of the finite differencing scheme. In three spatialdimensions high-order finite-difference schemes with wide stencils easily lead to adrastic decrease in performance caused by the additional computational load due tothe extra buffer points. This gain in accuracy pertains in particular to the phase ofthe evolution.In the present paper we report on a first step to significantly improve the accuracyof current finite-difference codes to evolve black-hole binaries by using sixth orderaccurate finite differencing operators in the bulk of the grid. We combine the sixthorder accurate derivative operators with fourth-order accurate dissipation operators[setting r = 3 in Eq. (3)] and Runge-Kutta time integrators, and aggressively reducethe number of AMR buffer zones compared to the number we would theoreticallyrequire for sixth order convergence. The penalty in computational cost is a rathermoderate 30% compared with our fourth-order code. (We have compared the averagespeed over 100 M of evolution time for our largest grid configuration, which produces inary black hole evolutions with sixth order finite differencing
3a 29% increase in computational cost; in our experience this is typical for our currentcode).
2. Summary of the “moving-puncture method” as implemented in theBAM code
There is a large freedom in writing the Einstein equations as a system of partialdifferential equations, and much research has gone into finding optimal choices. In thiswork we employ the currently most popular choice, the BSSN system [29, 30, 31, 32],which so far is the only system for which results have been reported of long-timeblack-hole binary simulations that do not rely on black hole excision as has been usedfor example in [1, 24, 25].We use the BSSN system together with the 1+log and gamma freezing coordinategauges [33, 34, 31] as described in [16] (choosing in particular the parameter η inthe gamma freezing shift condition as η = 2 as we have done previously). Thesegauge conditions allow the “punctures” to move across the grid (“moving puncture”approach [2, 3]) and allow an effective softening of the singularity in the metricassociated with an internal asymptotic region [35, 36, 37], which had been prohibitedby the traditional “fixed punctures” approach. The BSSN system is based on aconformal decomposition of the spatial geometry, writing the physical spatial metric as g ij = χ − ˜ g ij (following [2]). The blowup of the metric at the “punctures” is absorbedinto the conformal factor χ , which vanishes at the “puncture”.For our numerical evolutions we use the BAM [38, 39] code, which is designed tosolve partial differential equations on structured meshes, in particular a coupled systemof (typically hyperbolic) evolution equations and elliptic equations. The complexityof the equations is addressed by using a Mathematica package integrated into thecode, which produces C-code from Mathematica expressions in tensor notation.Using such a system as we do in BAM, or as has been discussed in detail for theCactus environment in [40] drastically simplifies the modification of complex codesfor black-hole binary simulations, as was required to adapt codes from the “fixedpuncture” to the “moving puncture” paradigm, or in the present case to implementthe improved numerical algorithms discussed here. The structure of the BAM codehas also made it straightforward to implement higher order finite differencing methodsThe computational domain is decomposed into rectangular boxes, following standarddomain-decomposition algorithms, and is parallelized with MPI [41]. Our meshrefinement algorithm is based on the standard Berger-Oliger algorithm, but withadditional buffer zones, along the lines of [42, 43] as described in [16] and summarizedin the next section. We essentially use a fixed-mesh-refinement strategy, with innerlevel refinement boxes following the motion of the black holes. Typically we use about10 refinement levels (refining the grid spacing by factors of 2), roughly half of whichfollow the movement of the black holes.In order to represent black holes in the initial data, we use the so-called “puncturemethod” [44]. For these data it is well understood how to write the constraintequations in a form suitable for numerical solution [45, 46]. Following the approachof [44] our initial data sets are chosen to be conformally flat with Bowen-Yorkextrinsic curvature. The momentum parameter in the Bowen-York extrinsic curvatureis determined from a quasi-equilibrium condition at third post-Newtonian order asdescribed in [16]. The elliptic constraint equations are solved in BAM with the pseudo-spectral collocation code described in [47]. AMR data are then obtained by barycentric inary black hole evolutions with sixth order finite differencing
3. Sixth order finite differencing
Our numerical evolution algorithm is based on a method-of-lines approach using finitedifferencing in space and explicit fourth-order accurate Runge-Kutta time stepping(with a fixed time step). We apply sixth-order accurate polynomial interpolation inspace between different refinement levels so that all spatial operations of the AMRmethod (i.e. restriction and prolongation) are sixth-order accurate, such that thesecond derivatives of interpolated values are at least fourth-order accurate. Althoughthe time stepping used for evolution is also fourth-order accurate through the Runge-Kutta integrator, there arises the additional issue of how to provide boundary valuesfor the intermediate time-levels of the Berger-Oliger algorithm that are not alignedin time with a coarser level (otherwise spatial interpolation can be used). Usinghigher than third-order interpolation has lead to spurious noise at mesh refinementboundaries as described in [16]. We therefore use third-order interpolation intime, which introduces a second-order error within the Berger-Oliger time-steppingscheme [16], which however is not noticeable in typical runs as we have checked byrunning with uniform (as opposed to Berger-Oliger non-uniform) time-stepping. Insummary, if the outer boundary is placed sufficiently far away and if time-interpolationerrors at refinement boundaries are small, then fourth-order convergence can beobserved.A relatively straightforward modification of the standard Berger-Oliger scheme isto replace the single-point refinement boundary by a buffer zone consisting of severalpoints, e.g., [42, 43, 48]. For a sufficiently large number of buffer zones (the product ofnumber of points in the stencil toward the mesh refinement boundary and the numberof source evaluations during a full time step of the coarse grid), no time interpolation isrequired and excellent results have been reported for this scheme [43] (note that specialmethods like [49] seem to achieve similar performance). For example, our fourth-orderRunge Kutta scheme requires four source evaluations, and if the lop-sided stencil withthree points in one direction, f ′ ( x ) = − f − − f + 18 f − f + f h − f (5) ( x ) h + O ( h ) (1)is used, then the numerical domain of dependence for a given point has a radius of 12points. Here and in the following we use the notation f j = f ( x + j ∆ x ). Therefore, itis possible to provide 12 buffer points at the refinement boundary and to perform oneRK4 time step with size three stencils that does not require any boundary updates.Only after the time step is completed, the buffer zones have to be repopulated. Inthe context of Berger-Oliger AMR, the buffer update is based on interpolation fromthe coarser levels. Since every second time step at level l coincides in time with level l −
1, one can provide 24 buffer points, perform two time steps, and then update thebuffer by interpolation in space. With 12 buffer points, one can interpolate in time toobtain data for the buffer points at intermediate time levels.To use fewer than 12 buffer points, we can interpolate into all buffer pointsbefore starting an RK4 update as described, and then evolve all points except the inary black hole evolutions with sixth order finite differencing / ≈ .
44 and 128 / ≈ .
69. For clarity, we alwaysquote grid sizes without buffer points, because this is the number of points owned by aparticular grid. Experimentally we have found that for the fourth-order case using justsix buffer points leads to very small differences compared to 12 buffer points, but evensmaller buffer zones lead to noticeable differences. For simulations with fourth-orderaccurate derivative operators, we have therefore chosen a standard setup of RK4 withdissipation and lop-sided advection stencils, 6 buffer points, quadratic interpolation intime, and Berger-Oliger time-stepping on all but the outermost grids following [16].As is common in numerical relativity, we use symmetric finite-difference stencilsfor all spatial derivatives but the advection terms associated with the shift vector,where we use lop-sided upwind stencils, see e.g. [50] for the fourth-order accuratecase.
In finite-difference codes targeted at smooth solutions of nonlinear hyperbolicequations, it is common practice to add artificial dissipation terms to all right-hand-sides of the time evolution equations, schematically written as ∂ t u → ∂ t u + Q u . (2)Such dissipation terms are very efficient at suppressing very high-frequency waves,which are not part of the physical solution. This may be necessary for numericalstability [51], but also to reduce numerical noise generated at mesh-refinementboundaries. As has become rather common in numerical relativity, we follow [51]and choose an operator ( Q ) of order 2 r as Q = σ ( − h ) r − ( D + ) r ( D − ) r / r , (3)for consistency with a 2 r − σ a parameter regulating thestrength of the dissipation. As we have done in the past with our fourth-order codewe choose the factor σ as σ = 0 . σ = 0 . inary black hole evolutions with sixth order finite differencing For the sixth-order case we find that several choices for the advection term stencilsyield stable evolutions, but the lop-sided upwind stencil which is closest to thesymmetric case yields (probably not surprisingly) by far the best accuracy, i.e. weuse f ′ ( x ) = 2 f − − f − − f + 80 f − f + 8 f − f h + 1105 d f ( x ) dx h + O ( h ) . (4)Alternative asymmetric choices would be f ′ ( x ) = − f − − f + 150 f − f + 50 f − f + 2 f h − d f ( x ) dx h + O ( h ) ,f ′ ( x ) = − f + 360 f − f + 400 f − f + 72 f − f )60 h + 17 d f ( x ) dx h + O ( h ) , for the stencils that deviate more from the symmetric choice. We can see that thefirst choice has the smallest leading error term. The symmetric stencil has an evensmaller error term, f ′ ( x ) = − f + f − − f − f − ) + 45( f − f − )60 h − d f ( x ) dx h + O ( h ) , but does not show equally robust results, as is common for solving advection equations.For non-advection derivative terms we again use the standard symmetric stencil,similarly for second derivatives in one direction we use the symmetric stencil f ′′ ( x ) = − f + 270( f + f − ) − f + f − ) + 2( f + f − ))180 h − f (8) ( x ) h + O ( h ) . For mixed derivatives, we use the stencils which result from a product of the symmetricsixth order accurate first derivative operators.
4. Results for long equal-mass evolutions
All runs are carried out with the symmetry ( x, y, z ) → ( − x, − y, z ) and ( x, y, z ) → ( x, y, − z ), reducing the computational cost by a factor of four. The Courant factor C = ∆ t/h l is kept constant, and is set to C = 1 / D = 12 M , the horizon mass for each individual hole is chosen inary black hole evolutions with sixth order finite differencing Table 1.
Grid setups used for convergence test simulations. The notation in the“Run” column is the same as we have used in [16]. The quantities h min and h max (rounded to three digits) denote the finest and coarsest grid spacings, and r max isthe location of the outer boundary (rounded to 4 digits), and all are in units of M .Also specified are the numbers of processors used, maximal memory requirementin GByte (to be precise, we quote the resident size of the program, i.e., the physicalmemory a task has used), and average speed in M /hour for the Kepler cluster atthe University of Jena (using Intel dual Woodcrest CPUs running at 2.66 GHz).The number in bold are used to indicate individual simulations throughout thispaper. Run h min h max r max procs. mem. (GByte) M/hour χ η =2 [5 × : 5 ×
96 : 6] 1 / . χ η =2 [5 × : 5 ×
112 : 6] 1 / . χ η =2 [5 × : 5 ×
128 : 6] 1 / . χ η =2 [5 × : 5 ×
144 : 6] 1 / . χ η =2 [5 × : 5 ×
160 : 6] 1 / . (cid:144) M - - - r h (cid:144) M Figure 1.
Left panel: Coordinate tracks of the puncture location of one black holefor simulations { , , } . Only in the last few orbits are differences betweenthe three runs discernable. Right panel: the waveform plotted as the real part of rh , as defined in [20]. as m i = 0 . M , which corresponds to puncture mass parameters of m p = 0 . M . Theinitial momenta are obtained as p = ± . M from a 3PN-accurate quasicircularity-condition as in [16].Our algorithm for gravitational wave extraction in terms of the Newman-Penrosescalar Ψ has been described in [16]. It is useful to write the signal in terms of atime-dependent amplitude and phase asΨ = A ( t ) e iϕ ( t ) , and define the gravitational wave frequency as ω = ˙ ϕ .The coordinate tracks of the puncture locations are shown in figure 1 forthe simulations for which we have obtained sixth-order convergence in the phase, inary black hole evolutions with sixth order finite differencing (cid:144) M0123456 D (cid:144) M
600 800 1000 1200 1400 1600 1800t (cid:144)
M255075100125150175200 Φ @ r a d i a n D Figure 2.
Left panel: coordinate distance of the black holes for the fourth orderversion of the 48-configuration and sixth-order simulations { , , } in the orderof increasing merger time. Right panel: the gravitational wave phase for the sameruns. The 72 simulation would not be distinguishable from the 80 simulation onthe scale shown here. Figure 3.
Convergence test for the gravitational wave phase. Plotted are thedifference between the 72 and 80 runs, and the difference between the 64 and 72runs rescaled for sixth-order convergence. Also shown is the convergence factordivided by 6, which shows a “glitch” around the time that the phase increasesvery sharply, and the error estimate after performing Richardson extrapolation.The left panel shows a linear scaling, the right panel shows the same plotwith a logarithmic scaling to emphasize the slow but clean exponential growth δϕ = 0 . . t/M of the phase error at intermediate times. { , , } . Only during the last two orbits are differences between the three runsdistinguishable. The orbital tracks show roughly 9 orbits before merger, and for thegravitational wave we obtain roughly 26 cycles before the ringdown signal becomes toonoisy at t ≈ M . The right panel of figure 1 shows the real part of the l = 2 , m = 2mode of the rescaled strain r h ( Ψ = ¨ h , compare [20]).Figure 2 shows the coordinate distance and gravitational wave phase for the blackholes for the { , } simulations (the 72 simulation would not be distinguishable fromthe 80 run on the scale shown here). We find that lower resolutions merge earlier, andthis is systematic for all the runs we have performed. For the { , , , } runs weobtain “merger times” of t = (1746 . , . , . , . M , and the maximum ofthe radiation power is reached at times t = (1818 . , . , . , . M . “Mergertime” here is understood only as a rough indicator of the time of merger, which isused for convergence tests, and is defined as the time when the coordinate distanceof the punctures drops below 0 . M . For the { , , } runs the merger times showconvergence at order 5 .
55, the radiation peak times show convergence at order 5 . inary black hole evolutions with sixth order finite differencing - - -
500 0t (cid:144)
M0.000010.00010.0010.01 È r Y È @ M - D - - - - - - -
250 0t (cid:144)
M0.10.20.30.40.50.6 Ω M Figure 4.
The l = 2 , m = 2 mode of the wave signal is split into the absolute valueof Ψ , (left panel) and the wave frequency ω (right panel). Both panels showthe simulations { , , } , aligned in time to coincide at the peak of | Ψ | . Thecurves are clipped at early times, where they are very noisy due to the smallnessof the signal and finite differencing error in computing the wave frequency fromthe phase. (cid:144) M - - D H È r Y È L @ M - D H L - H LHH L - H LL(cid:144) Φ GW @ radian D - - D È r Y È @ M - D H L - H LHH L - H LL(cid:144)
Figure 5.
Convergence plot for the wave amplitude | ψ | in the l = 2 , m = 2mode. Both panels show the difference between the 72 and 80 runs and differencesbetween the 64 and 72 runs rescaled for sixth-order convergence. In the left paneldata at different resolutions are compared at the same coordinate time, whichleads to a seeming loss of convergence near the radiation peak, which is due tothe relatively large phase error. In the right panel the data are compared atthe same value of the gravitational wave phase, which restores clean sixth orderconvergence. Richardson extrapolation with convergence order 5 . ≈ M for both times. Note that oscillations, which are probably mainly due to eccentricity,can clearly be seen in the black hole distance, and also in the wave amplitude shownin figure 4. A method to reduce eccentricity will be discussed in [52].We have obtained roughly sixth order convergence for the gravitational wavephase between about t = 1000 M and t = 1800 M , as shown in figure 3. At earliertimes the convergence factor becomes very noisy due to the smallness of the signal.Shortly before the merger the convergence factor “glitches” to a value of roughly 7.This problem can also clearly be seen in the convergence of the radiation frequencyand amplitude as shown in figures 6 and 5. This “glitch” appears when the frequencyand phase increase very sharply, and small phase errors have a large effect. Thelogarithmic scaling version of the convergence plot in figure 3 shows a slow and ratherclean exponential growth for the phase error at intermediate times, a nonlinear fitfor 300 M ≤ t ≤ M yields δϕ = 0 . . t/M ) for the phase error. This inary black hole evolutions with sixth order finite differencing
80 85 90 95 100 105 110 Φ GW @ radian D - - - - - D s e p a r a ti on @ M D H L - H LHH L - H LL(cid:144) (cid:144)
M0.00010.0010.010.11 ∆ Φ @ w a v ec y c l e s D Figure 6.
Left panel: Reparametrisation of the black hole distance bygravitational wave phase yields clean sixth order convergence until the time ofmerger (which occurs roughly at the peak of the error) – the difference betweenthe 72 and 80 runs and differences between the 64 and 72 runs rescaled for sixth-order convergence lie essentially on top of each other. Right panel: The differencein wave phase between the 4th and 6th order versions of the lowest resolution 48configuration shows the fourth-order algorithm “falling behind”. observation provides one way to optimize numerical methods in the inspiral phase,without evolving all the way to the merger. In the last stage of the inspiral thephase error grows very rapidly. We have noted this previously in a different contextin [16], where we have compared different methods to provide quasicircular inspiraldata. Small changes on the order of 1% of the initial momenta have lead to drasticchanges of ≈ M in merger time. In order to clarify the convergence behaviourof our code, we have applied a new technique to check for convergence in situationswhere the numerical error is dominated by phase shift: We first perform a convergenceanalysis for the dependence of the gravitational (or orbital) phase on code time, andthen perform a standard convergence test on a quantity like the puncture separation(Fig. (6)) or wave amplitude (Fig. (5)) regarding the functional dependence of thisquantity on the phase.An important question aside from convergence is how the sixth order and fourthorder algorithms compare in absolute numbers. For this purpose we have re-run the 48configuration with our standard fourth order algorithm. We find that already at thislow resolution the sixth order algorithm is superior, as shown in figures 2 and 6, whileat higher resolutions the larger convergence factor increases the gain in accuracy.
5. Conclusions
We have described a minimal extension of the fourth-order accurate evolutionalgorithm described in [16], where by replacing spatial fourth-order accuratedifferencing operators in the bulk of the grid by sixth-order accurate stencils, wegain drastic improvements in accuracy for the phase in long simulations of equal-mass inspiral. The crucial technical point regarding the choice of sixth-order accuratefinite difference operators has been a specific choice for the advection stencil, which isused to discretize Lie derivative terms with respect to the shift vector in the Einsteinequations. Using this method we have demonstrated evolutions of about nine orbits or1800 M with a phase error of approximately 4 M in the merger time, requiring ≈ inary black hole evolutions with sixth order finite differencing Acknowledgments
References [1] Frans Pretorius. Evolution of binary black hole spacetimes.
Phys. Rev. Lett. , 95:121101, 2005.[2] Manuela Campanelli, Carlos O. Lousto, Pedro Marronetti, and Yosef Zlochower. Accurateevolutions of orbiting black-hole binaries without excision.
Phys. Rev. Lett. , 96:111101, 2006.[3] John G. Baker, Joan Centrella, Dae-Il Choi, Michael Koppitz, and James van Meter.Gravitational wave extraction from an inspiraling configuration of merging black holes.
Phys.Rev. Lett. , 96:111102, 2006.[4] Susan G. Hahn and Richard W. Lindquist. The two body problem in geometrodynamics.
Ann.Phys. , 29:304–331, 1964.[5] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, J. van Meter, and M. C. Miller. Getting akick out of numerical relativity.
Astrophys. J , 2007. astro-ph/0603204.[6] J. A. Gonz´alez, U. Sperhake, B. Br¨ugmann, M. Hannam, and S. Husa. The maximum kick fromnonspinning black-hole binary inspiral.
Phys. Rev. Lett. , 2006. gr-qc/0610154.[7] M. Campanelli, C. O. Lousto, and Y. Zlochower. Gravitational radiation from spinning-black-hole binaries: The orbital hang up.
Phys. Rev. D , 74:041501, 2006. gr-qc/0604012.[8] M. Campanelli, C. O. Lousto, and Y. Zlochower. Spin-orbit interactions in black-hole binaries.
Phys. Rev. D , 74:084023, 2006. gr-qc/0608275.[9] F. Herrmann, I. Hinder, D. Shoemaker, P. Laguna, and R. Matzner. Gravitational recoil fromspinning binary black hole mergers. 2007. gr-qc/0701143.[10] Michael Koppitz, Denis Pollney, Christian Reisswig, Luciano Rezzolla, Jonathan Thornburg,Peter Diener, and Erik Schnetter. Getting a kick from equal-mass binary black hole mergers.2007. gr-qc/0701163.[11] J. A. Gonzalez, M. D. Hannam, U. Sperhake, B. Br¨ugmann, and S. Husa. Supermassive kicksfor spinning black holes. 2007. gr-qc/0702052.[12] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Merritt. Large Merger Re-coils and Spin Flips From Generic Black-Hole Binaries. 2007. Final version, .[13] F. Herrmann, I. Hinder, D. Shoemaker, and P. Laguna. Presented at “NRm3PN Meeting”, St.Louis, February 8-11, 2007.[14] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter. Binary black hole mergerdynamics and waveforms.
Phys. Rev. D , 73:104002, 2006. gr-qc/0602026.[15] J. G. Baker, J. R. van Meter, S. T. McWilliams, J. Centrella, and B. J. Kelly. Consistency ofpost-newtonian waveforms with numerical relativity. 2006. gr-qc/0612024.[16] Bernd Br¨ugmann, Jos´e A. Gonz´alez, Mark Hannam, Sascha Husa, Ulrich Sperhake, andWolfgang Tichy. Calibration of moving puncture simulations. 2006. gr-qc/0610128. inary black hole evolutions with sixth order finite differencing [17] Alessandra Buonanno, Gregory B. Cook, and Frans Pretorius. Inspiral, merger and ring-downof equal mass black-holes binaries. 2006. gr-qc/0610122.[18] Emanuele Berti et al. Inspiral, merger and ringdown of unequal mass black hole binaries: Amultipolar analysis. 2007. gr-qc/0703053.[19] Yi Pan et al. A data-analysis driven comparison of analytic and numerical coalescing binarywaveforms: Nonspinning case. 2007. arXiv:0704.1964 [gr-qc].[20] P. Ajith et al. Phenomenological template family for black-hole coalescence waveforms. 2007.[21] Manuela Campanelli, C. O. Lousto, and Y. Zlochower. The last orbit of binary black holes. Phys. Rev. D , 73:061501(R), 2006.[22] John G. Baker, Joan Centrella, Dae-Il Choi, Michael Koppitz, and James van Meter. Binaryblack hole merger dynamics and waveforms.
Phys. Rev. D , 73:104002, 2006.[23] U. Sperhake. Binary black-hole evolutions of excision and puncture data. 2006. gr-qc/0606079.[24] Mark A. Scheel, Harald P. Pfeiffer, Lee Lindblom, Lawrence E. Kidder, Oliver Rinne, andSaul A. Teukolsky. Solving Einstein’s equations with dual coordinate frames.
Phys. Rev. D ,74:104006, 2006.[25] Harald P. Pfeiffer, Duncan Brown, Lawrence E. Kidder, Lee Lindblom, Geoffrey Lovelance, andMark A. Scheel. Reducing orbital eccentricity in binary black hole simulations. 2007.[26] Mark Scheel. Binary black hole evolutions using spectral methods. Talk at the NFNRConference, Potsdam, Germany, July 17-21 2006.[27] John G. Baker et al. Binary black hole late inspiral: Simulations for gravitational waveobservations. 2006. gr-qc/0612117.[28] Bernd Br¨ugmann, Jos´e A. Gonz´alez, Mark Hannam, Sascha Husa, and Ulrich Sperhake.Presented at LSC Meeting, Boston, November 4-6, 2006.[29] Thomas W. Baumgarte and Stuart L. Shapiro. On the numerical integration of Einstein’s fieldequations.
Phys. Rev. D , 59:024007, 1999.[30] Masaru Shibata and Takashi Nakamura. Evolution of three-dimensional gravitational waves:Harmonic slicing case.
Phys. Rev. D , 52:5428, 1995.[31] Miguel Alcubierre, Bernd Br¨ugmann, Peter Diener, Michael Koppitz, Denis Pollney, EdwardSeidel, and Ryoji Takahashi. Gauge conditions for long-term numerical black hole evolutionswithout excision.
Phys. Rev. D , 67:084023, 2003.[32] Carsten Gundlach and Jose M. Martin-Garcia. Hyperbolicity of second-order in space systemsof evolution equations.
Class. Quantum Grav. , 23:S387–S404, 2006.[33] Miguel Alcubierre and Bernd Br¨ugmann. Simple excision of a black hole in 3+1 numericalrelativity.
Phys. Rev. D , 63:104006, 2001.[34] John Baker, Bernd Br¨ugmann, Manuela Campanelli, Carlos O. Lousto, and Ryoji Takahashi.Plunge waveforms from inspiralling binary black holes.
Phys. Rev. Lett. , 87:121103, 2001.[35] Mark Hannam, Sascha Husa, Denis Pollney, Bernd Br¨ugmann, and Niall ´O Murchadha.Geometry and regularity of moving punctures. 2006. gr-qc/0606099.[36] Mark Hannam, Sascha Husa, Niall ´O Murchadha, Bernd Br¨ugmann, Jos´e A. Gonz´alez, andUlrich Sperhake. Where do moving punctures go?
Journal of Physics: Conference series ,2007. in press.[37] J. David Brown. Puncture evolution of schwarzschild black holes. 2007. arXiv:0705.1359 [gr-qc].[38] B. Br¨ugmann. Binary black hole mergers in 3D numerical relativity.
Int. J. Mod. Phys. D ,8:85, 1999.[39] Bernd Br¨ugmann, Wolfgang Tichy, and Nina Jansen. Numerical simulation of orbiting blackholes.
Phys. Rev. Lett. , 92:211101, 2004.[40] Sascha Husa, Ian Hinder, and Christiane Lechner. Kranc: a Mathematica application to generatenumerical codes for tensorial evolution equations.
Comput. Phys. Comm. , 174:983–1004,2006.[41] W. Gropp, E. Lusk, N. Doss, and A. Skjellum. A high-performance, portable implementation ofthe MPI message passing interface standard.
Parallel Computing , 22(6):789–828, September1996.[42] Erik Schnetter, Scott H. Hawley, and Ian Hawke. Evolutions in 3D numerical relativity usingfixed mesh refinement.
Class. Quantum Grav. , 21(6):1465–1488, 21 March 2004.[43] Luis Lehner, Steven L. Liebling, and Oscar Reula. AMR, stability and higher accuracy. 2005.[44] S. Brandt and B. Br¨ugmann. A simple construction of initial data for multiple black holes.
Phys. Rev. Lett. , 78(19):3606–3609, 1997.[45] R. Beig and N. O’Murchadha.
Class. Quantum Grav. , 11:419, 1994.[46] S. Dain and H. Friedrich. Asymptotically flat initial data with prescribed regularity at infinity.
Comm. Math. Phys. , 222:569, 2001. gr-qc/0102047.[47] Marcus Ansorg, Bernd Br¨ugmann, and Wolfgang Tichy. A single-domain spectral method for inary black hole evolutions with sixth order finite differencing black hole puncture data. Phys. Rev. D , 70:064011, 2004.[48] P. Csizmadia. Testing a new mesh refinement code in the evolution of a spherically symmetricKlein-Gordon field.
Int. J. Mod. Phys. D , 15:107–119, 2006.[49] John G. Baker and James R. Meter. Reducing reflections from mesh refinement interfaces innumerical relativity.
Phys. Rev. D , 72:104010, 2005.[50] Y. Zlochower, J. G. Baker, M. Campanelli, and C. O. Lousto. Accurate black hole evolutionsby fourth-order numerical relativity.
Phys. Rev. D , 72:024021, 2005.[51] Heinz Otto Kreiss and Joseph Oliger.
Methods for the approximate solution of time dependentproblems . GARP publication series No. 10, Geneva, 1973.[52] Sascha Husa, Mark Hannam, Jos´e A. Gonz´alez, Ulrich Sperhake, and Bernd Br¨ugmann.Reducing eccentricity in black-hole binary evolutions with initial parameters from post-newtonian inspiral. In preparation.[53] Mark Hannam, Sascha Husa, Jos´e A. Gonz´alez, Ulrich Sperhake, and Bernd Br¨ugmann. 2007.In preparation.[54] Gioel Calabrese, Ian Hinder, and Sascha Husa. Numerical stability for finite differenceapproximations of Einstein’s equations.