Black hole puncture initial data with realistic gravitational wave content

TL;DR

This paper introduces improved initial data for non-spinning black-hole binaries that incorporate realistic gravitational wave content, enhancing the accuracy of numerical simulations and bridging the gap between inspiral and merger waveforms.

Contribution

The authors explicitly calculate integral terms to include realistic gravitational waves in initial data, improving accuracy in the far zone for black-hole binary simulations.

Findings

Improved initial data includes realistic gravitational waves.

Agreement of metric with quadrupole-approximation waveforms.

Enables direct connection between inspiral and merger waveforms.

Abstract

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

Figures

Equations

$$g_{ij}$$

$${\pi }_i^i$$

$${\varphi }^{(2)}=4G\sum _{A=1}^N\frac{m_A}{r_A},$$

$$h_{ij}^{\mathrm{TT}}=-{\delta }_{ij}^{\mathrm{TT}kl}{\square }_{ret}^{-1}{s}_{kl},$$

$$h_{ij}^{\mathrm{TT}}$$

$$s_{kl}=16\pi G\sum _A\frac{p_{Ak}{p}_{Al}}{m_A}\delta (x-x_A)+\frac{1}{4}{\varphi }_{,k}^{(2)}{\varphi }_{,l}^{(2)}$$

$${\bar{s}}_{kl}$$

$$h_{ij,(div)}^{\mathrm{TT}}=-{\delta }_{ij}^{\mathrm{TT}kl}({\square }_{ret}^{-1}-{\Delta }^{-1})(s_{kl}-{\bar{s}}_{kl})\sim O(v/c{)}^5$$

$$h_{ij}^{\mathrm{TT}}=h_{ij}^{\mathrm{TT}(\mathrm{NZ})}-{\delta }_{ij}^{\mathrm{TT}kl}({\square }_{ret}^{-1}-{\Delta }^{-1}){\bar{s}}_{kl}+h_{ij,(div)}^{\mathrm{TT}},$$

$$h_{ij}^{\mathrm{TT}}$$

$$h_{ij}^{\mathrm{TT}(\mathrm{NZ})}=h_{ij}^{\mathrm{TT}(4)}+h_{ij}^{\mathrm{TT}(5)}+O(v/c{)}^6.$$

$$h^{\mathrm{TT}(4)ij}$$

$$h_{+}^{(\mathrm{NZ})}$$

$$h_{\times }^{(\mathrm{NZ})}$$

$$h_{+}$$

$$h_{\times }$$

$$\Phi (\tau )=\Phi (t_c)-\frac{1}{\nu }{\Theta }^{5/8},$$

$$\Omega (\tau )=\frac{1}{8GM}{\Theta }^{-3/8},$$

$$r_{12}$$

$$p$$

$$h_{ij}^{\mathrm{TT}}$$

$$H_{ij}^{\mathrm{TT}A}[\vec{u}]$$

$$[u_i{u}_j{]}^{\mathrm{TT}}$$

$$t-t_A^r-r_A(t_A^r)=0,$$

$$H_{\mathrm{TT}A}^{ij}[\vec{u}]$$

$$H_{\mathrm{TT}A}^{ij}[\vec{u};t]$$

$$H_{\mathrm{TT}A}^{ij}[\vec{u};t_A^r]$$

$$H_{\mathrm{TT}A}^{ij}[\vec{u};t_A^r\to t]$$

$$h_{ij}^{\mathrm{TT}(pot,now)}$$

$$h_{ij}^{\mathrm{TT}(pot,4)}+h_{ij}^{\mathrm{TT}(pot,now)}$$

$$\frac{r_{12}(\Omega )}{GM}$$

$$\frac{p(\Omega )}{M\nu }$$

$$K^{ij}$$

$$\mathcal{S}\approx 1-3\frac{{\psi }_0{\psi }_4}{{\psi }_2^2}+O({ϵ}^3),$$

Full text

I Introduction

Post-Newtonian (PN) methods have played a fundamental role in our understanding of the astrophysical implications of Einstein’s theory of general relativity. Most importantly, they have been used to confirm that the radiation of gravitational waves accounts for energy loss in known binary pulsar configurations. They have also been used to create templates for the gravitational waves emitted from compact binaries which might be detected by ground-based gravitational wave observatories, such as LIGO Vogt (1992); Abbott et al. (2004), and the NASA/ESA planned space-based mission, LISA Bender et al. (1998); Danzmann and Rudiger (2003). However, PN methods have not been extensively used to provide initial data for binary evolution in numerical relativity, nor, until recently (see Buonanno et al. (2007); Berti et al. (2007)), have they been extensively studied so that their limitations could be well identified and the results of numerical relativity independently confirmed.

Until the end of 2004, the field of numerical relativity had been struggling to compute even a single orbit for a black-hole binary (BHB). Although debate occurred on the advantages of one type of initial data over another, the primary focus within the numerical relativity community was on code refinement which would lead to more stable evolution. Astrophysical realism was very much a secondary issue. However, this situation has radically changed in the last few years with the introduction of two essentially independent, but equally successful techniques: the generalized harmonic gauge (GHG) method developed by Pretorius  Pretorius (2005) and the “moving puncture” approach, independently developed by the UTB and NASA Goddard groups Campanelli et al. (2006a); Baker et al. (2006a). Originally introduced by Brandt & Brügmann Brandt and Brügmann (1997) in the context of initial data, the puncture method explicitly factored out the singular part of the metric. When used in numerical evolution in which the punctures remained fixed on the numerical grid, it resulted in distortions of the coordinate system and instabilities in the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) Shibata and Nakamura (1995); Baumgarte and Shapiro (1999) evolution scheme. The revolutionary idea behind the moving puncture approach was precisely, not to factor out the singular part of the metric, but rather evolve it together with the regular part, allowing the punctures to move freely across the grid with a suitable choice of the gauge.

A golden age for numerical relativity is now emerging, in which multiple groups are using different computer codes to evolve BHBs for several orbits before plunge and merger Brügmann et al. (2004); Campanelli et al. (2006b); Pretorius (2006); Baker et al. (2006b); Brügmann et al. (2006); Scheel et al. (2006); Marronetti et al. (2007); Tichy (2006); Pfeiffer et al. (2007). Comparison of the numerical results obtained from these various codes has taken place  Baker et al. (2007a); Thornburg et al. (2007); NRw , and comparison with PN inspiral waveforms has also been carried out with encouraging success Buonanno et al. (2007); Berti et al. (2007); Baker et al. (2007b, 2006c). The application of successful numerical relativity tools to study some important astrophysical properties (e.g. precession, recoil, spin-orbit coupling, elliptical orbits, etc) of spinning and/or unequal mass-black hole systems is currently producing extremely interesting new results  Campanelli (2005); Herrmann et al. (2006); Baker et al. (2006d); Campanelli et al. (2006c, d, 2007a); Gonzalez et al. (2007a); Herrmann et al. (2007); Campanelli et al. (2007b); Koppitz et al. (2007); Gonzalez et al. (2007b); Choi et al. (2007); Baker et al. (2007c); Pretorius and Khurana (2007); Campanelli et al. (2007c); Tichy and Marronetti (2007). It now seems that the primary obstacle to further progress is simply one of computing power. In this new situation, it is perhaps time to return to the question of what initial data will best describe an astrophysical BHB.

To date, the best-motivated description of pre-merger BHBs has been supplied by PN methods. We might expect, then, that a PN-based approach would give us the most astrophysically correct initial data from which to run full numerical simulations. In practice, PN results are frequently obtained in a form ill-adapted to numerical evolution. PN analysis often deals with the full four-metric, in harmonic coordinates; numerical evolutions frequently use ADM-type coordinates, with a canonical decomposition of the four-metric into a spatial metric and extrinsic curvature.

Fortunately, many PN results have been translated into the language of ADM by Ohta, Damour, Schäfer and collaborators. Explicit results for 2.5PN BHB data in the near zone were given by Schäfer Schäfer (1985) and Jaranowski & Schäfer (JS) Jaranowski and Schäfer (1998), and these were implemented numerically by Tichy *et al. *​Tichy et al. (2003a). Their insight was that the ADM-transverse-traceless (TT) gauge used by Schäfer was well-adapted to a puncture approach. To facilitate comparison with this earlier work Tichy et al. (2003a), we continue to use the results of Schäfer and co-workers, anticipating that higher-order PN results should eventually become available in a useful form.

The initial data provided previously by Tichy *et al. *​already include PN information. They are accurate up to order $(v/c)^{5}$ in the near zone ($r\ll\lambda$), but the accuracy drops to order $(v/c)^{3}$ in the far zone ($r\gg\lambda$) [here $\lambda\sim\pi\sqrt{r_{12}^{3}/G(m_{1}+m_{2})}$ is the gravitational wavelength]. These data were incomplete in the sense that they did not include the correct TT radiative piece in the metric, and thus did not contain realistic gravitational waves.

In this paper, we revisit the PN data problem in ADM-TT coordinates, with the aim of supplying Numerical Relativity with initial BHB data that extend as far as necessary, and contain realistic gravitational waves. To do this, we have evaluated the “missing pieces” of Schäfer’s TT metric for the case of two non-spinning particles. We have analyzed the near- and far-zone behavior of these data, and incorporated them numerically in the Cactus Cac framework. In principle, the most accurate PN metric available could be used at this step, but it is not currently available in ADM-TT form.

The remainder of this paper is laid out as follows. In Section II, we summarize the results of Schäfer (1985) Schäfer (1985), and Jaranowski & Schäfer (1997) Jaranowski and Schäfer (1998) and their application by Tichy *et al. *​(2003) Tichy et al. (2003a), to the production of puncture data for numerical evolution. In Section III, we describe briefly the additional terms necessary to complete $h^{{\rm TT}}$ to order $(v/c)^{4}$, deferring details to the Appendix. In Section IV, we study the full data both analytically and numerically. Section V summarizes our results, and lays the groundwork for numerical evolution of these data, to be presented in a subsequent article.

II ADM-TT Gauge in Post-Newtonian Data

The “ADM-TT” gauge Ohta et al. (1974); Schäfer (1985) is a 3+1 split of data where the three-metric differs from conformal flatness precisely by a TT radiative part:

[TABLE]

The fields $\phi$, $\pi^{ij}$ and $h^{{\rm TT}}_{ij}$ can all be expanded in a post-Newtonian series. Solving the constraint equations of 3+1 general relativity in this gauge, Schäfer (1985); Jaranowski and Schäfer (1998) obtained explicit expressions valid up to $O(v/c)^{5}$ in the near zone, incorporating an arbitrary number of spinless point particles, with arbitrary masses $m_{A}$. For $N$ particles, the lowest-order contribution to the conformal factor is111We explicitly include the gravitational constant $G$ in all expressions here, as the standard convention $G=1$ used in Numerical Relativity differs from the convention $16\pi G=1$ employed by Schäfer (1985); Jaranowski and Schäfer (1998).:

[TABLE]

where $r_{A}=\sqrt{\vec{x}-\vec{x}_{A}}$ is the distance from the field point to the location of particle $A$.

In principle $h^{{\rm TT}}_{ij}$ is computed from

[TABLE]

where $\Box_{ret}^{-1}$ is the (flat space) inverse d’Alembertian (with a “no-incoming-radiation” condition Fock (1964)), $s_{kl}$ is a non-local source term and $\delta^{{\rm TT}\,kl}_{ij}$ is the TT-projection operator. In order to compute $h^{{\rm TT}}_{ij}$ we first rewrite Eq. (4) as

[TABLE]

Note that the near-zone approximation $h^{{\rm TT}\,({\rm NZ})}_{ij}$ of $h^{{\rm TT}}_{ij}$ has already been computed in Schäfer (1985) up to order $O(v/c)^{4}$ (see also Eq. 12 below). The last term in Eq. (5) is difficult to compute because

[TABLE]

is a non-local source. However, we can approximate $s_{kl}$ by

[TABLE]

and show that

[TABLE]

in the near zone. Furthermore, outside the near zone $h^{{\rm TT}}_{ij,(div)}\sim 1/r^{2}$, so that $h^{{\rm TT}}_{ij,(div)}$ falls off much faster than rest of $h^{{\rm TT}}_{ij}$, which falls off like $1/r$. Hence

[TABLE]

where $h^{{\rm TT}}_{ij,(div)}$ can be neglected if we only keep terms up to $O(v/c)^{4}$ generally, and $O(1/r)$ at infinity.

The full expression for $h^{{\rm TT}}_{ij}$ for $N$ interacting point particles from Eq. (4.3) of Schäfer (1985) is:

[TABLE]

The first term in (10), $h^{{\rm TT}\,({\rm NZ})}_{ij}$ can be expanded in $v/c$ as

[TABLE]

The leading order term at $O(v/c)^{4}$, is given explicitly by Eq. (A20) of Jaranowski and Schäfer (1998):

[TABLE]

where $s_{AB}\equiv r_{A}+r_{B}+r_{AB}$. The other two terms in Eq. (10) can be shown to be small in the near zone ($r\ll\lambda$, where the characteristic wavelength $\lambda\sim 100M$ for $r_{AB}\sim 10M$). However, $h^{{\rm TT}\,({\rm NZ})}_{ij}$ is only a valid approximation to $h^{{\rm TT}}_{ij}$ in the near zone, and becomes highly inaccurate when used further afield.

Setting aside these far-field issues, Tichy *et al. *​Tichy et al. (2003a) applied Schäfer’s formulation, in the context of a black-hole binary system, to construct initial data that are accurate up to $O(v/c)^{5}$ in the near zone. They noted that the ADM-TT decomposition was well-adapted to the use of a puncture approach to handle black-hole singularities. This approach is essentially an extension of the method introduced in Brandt and Brügmann (1997). It allows a simple numerical treatment of the black holes without the need for excision.

The PN-based puncture data of Tichy *et al. *​ have not been used for numerical evolutions. This is in part because these data, just like standard puncture data Brandt and Brügmann (1997); Tichy et al. (2003b); Tichy and Brügmann (2004); Ansorg et al. (2004), do not contain realistic gravitational waves in the far zone: $h^{{\rm TT}\,({\rm NZ})}_{ij}$ does not even vaguely agree with the 2PN approximation to the waveform amplitude nor with the quadrupole approximation to the waveform phase for realistic inspiral.

To illustrate this, we restrict to the case of two point sources, and compute the “plus” and “cross” polarizations of the near-zone approximation for $h^{{\rm TT}}_{ij}$:

[TABLE]

For comparison, the corresponding polarizations of the quadrupole approximation for the gravitational-wave strain are given by (paraphrasing Eq. (3.4) of Finn and Chernoff (1993)):

[TABLE]

where $\mathcal{M}\equiv\nu^{3/5}\,M$ is the “chirp mass” of the binary, given in terms of the total PN mass of the system $M=m_{1}+m_{2}$, and the symmetric mass ratio $\nu=m_{1}m_{2}/M^{2}$. The angle $\theta$ is the “inclination angle of orbital angular momentum to the line of sight toward the detector”; that is, just the polar angle to the field point, when the binary moves in the $x$-$y$ plane. $\Phi_{{\rm GW}}$ and $f_{{\rm GW}}$ are the phase and frequency of the radiation at time $t$, exactly twice the orbital phase $\Phi(t-r)$ and orbital frequency $\Omega(t-r)/2\pi$.

The lowest-order PN prediction for radiation-reaction effects yields a simple inspiral of the binary over time, with orbital phasing given by

[TABLE]

where $\Theta\equiv\nu\,(t_{c}-\tau)/5\,GM$, $M$ and $\nu$ are given below (16), and $t_{c}$ is a nominal “coalescence time”. To evaluate (13-14), we need the transverse momentum $p$ corresponding to the desired separation $r_{12}$. The simplest expression for this is the classical Keplerian relation, which we give parameterized by $\Omega(\tau)$:

[TABLE]

In Fig. 1 we compare the plus polarization of the two waveforms (13) and (15) at a field point $r=100M$, $\theta=\pi/4$, $\phi=0$, for a binary in the $x$-$y$ plane, with initial separation $r_{12}=10M$. The orbital frequency of the binary is related to the separation $r_{12}$ and momenta $p$ entering (13) by (19-20). To this level of approximation, the binary has a nominal PN coalescence time $t_{c}\approx 780M$. As might have been anticipated, both phase and amplitude of $h^{{\rm TT}\,(4)}_{ij}$ are wrong outside the near zone. This means that the data constructed from $h^{{\rm TT}\,(4)}_{ij}$ have the wrong wave content, but nevertheless these data are still accurate up to order $(v/c)^{3}$ in the far zone.

It is evident from the present-time dependence of (12) that it cannot actually contain any of the past history of an inspiralling binary. We would expect that a correct “wave-like” contribution should depend rather on the retarded time of each contributing point source. It seems evident that the correct behavior must, in fact, be contained in the as-yet unevaluated parts of (10). The requisite evaluation is what we undertake in the next section.

III Completing the Evaluation of $h^{{\rm TT}}_{ij}$

To move forward, we will simplify (10) and (12) to the case of only two particles. Then (10) reduces to:

[TABLE]

where

[TABLE]

Here, the “TT projection” is effected using the operator $P_{i}^{j}\mathrel{\mathop{:}}=\delta_{i}^{j}-k_{i}\,k^{j}/k^{2}$. For an arbitrary spatial vector $\vec{u}$,

[TABLE]

Details on the evaluation of these terms are presented in Appendix A. After calculation, we write the result as a sum of terms evaluated at the present field-point time $t$, the retarded time $t_{A}^{r}$ defined by

[TABLE]

and integrals between $t_{A}^{r}$ and $t$,

[TABLE]

where the three parts are given by:

[TABLE]

In Fig. 2, we show the retarded times calculated for each particle, as measured at points along the $x$ axis, for the same orbit as in Fig. 1. We also show the corresponding retarded times for a binary in an exactly circular orbit. Since the small-scale oscillatory effect of the finite orbital radius would be lost by the overall linear trend, we have multiplied by the orbital radius.

III.1 Reconciling with Jaranowski & Schäfer’s $h^{{\rm TT}\,(4)}_{ij}$

From the derivation above it is clear that $h^{{\rm TT}}_{ij}$ includes retardation effects, so it will not depend solely on the present time. We might even expect that all “present-time” contributions should vanish individually, or should cancel out. It can be seen easily from (27) that the “$t$” part of the second and third terms of Eq. (22) exactly cancel out the “kinetic” part (first line) of Eq. (12). Thus, we can simply remove that line in Eq. (12), and use the “$t^{r}$” part instead. One may similarly inquire whether the “$t$” parts of the fourth and fifth terms of Eq. (22) above,

[TABLE]

also cancel the remaining, “potential” parts of Eq. (12). The answer is “not completely”; expanding in powers of $1/r$, we find:

[TABLE]

where $W\equiv\sin\theta\,\cos(\phi-\Phi(t))$, and $\Phi(t)$ is the orbital phase of particle 1 at the present time $t$. That is, the “new” contribution cancels the $1/r$ and $1/r^{2}$ pieces of $h^{{\rm TT}\,(4)}_{ij}$ entirely. In the far zone the result is thus smaller than the $h^{{\rm TT}}_{ij,(div)}$ term which we are ignoring everywhere, since it is small both in the near and the far zone Schäfer (1985).

We note here two general properties of the contributions to the full $h^{{\rm TT}}_{ij}$.

In the near zone $h^{{\rm TT}\,(4)}_{ij}$ is the dominant term since all other terms arise from $(\Box_{ret}^{-1}-\Delta^{-1})s_{kl}$. Thus all other terms must cancel within the accuracy of the near-zone approximation. 2. 2.

$h^{{\rm TT}\,(4)}_{ij}$ is wrong far from the sources; thus, the new corrections should “cancel” $h^{{\rm TT}\,(4)}_{ij}$ entirely, far from sources. Note, however, that while $h_{ij}=-\Box_{ret}^{-1}s_{kl}$ depends only on retarded time, its TT-projection $h^{{\rm TT}}_{ij}=\delta^{{\rm TT}\,kl}_{ij}h_{kl}$ has a more complicated causal structure; E.g. the finite time integral comes from applying the TT-projection. [Proof: Even if we had a source given exactly by $\bar{s}_{kl}$, $h^{{\rm TT}\,(4)}_{ij}$ would depend only the present time, $h_{ij}$ would depend only on retarded time, and $h^{{\rm TT}}_{ij}$ would (as we have computed) contain a finite time integral term.]

Additionally, the full $h^{{\rm TT}}_{ij}$ agrees well with quadrupole predictions, which we demonstrate in Section IV.

IV Numerical Results and Invariants

IV.1 Phasing and Post-Keplerian Relations

It has been known for some time (see for example Cutler et al. (1993)) that gravitational wave phase plays an even more important part in source identification than does wave amplitude. In PN work, phase and amplitude are estimated somewhat separately; the amplitude requires knowledge of the time-dependent multipoles, used in developing the the full metric, while the phase can be relatively simply approximated from the orbital equations of motion, taking into account the gravitational wave flux at infinity to evolve the orbital parameters Tichy et al. (2000).

The quadrupole waveform introduced for the comparison in Fig. 1 had an amplitude accurate to $O(v/c)^{4}$ and the simplest available time evolution for the phase. Waveform phase is a direct consequence of orbital phase. To lowest order, we could have assumed a binary moving in a circular orbit (of zero eccentricity) since, up to 2PN order, we can have circular orbits, where the linear momentum, $p$, of each particle is related to the separation $r_{12}$ by, say, Eq. (24) of Tichy et al. (2003a). Nevertheless, circular orbits are physically unrealistic – since radiation reaction will lead to inspiral and merger of the particles – and Eqs. (17-18) already include leading-order radiation-reaction effects. Moreover, the phase errors that would accrue from using purely circular orbits would be larger, the further from the sources we tried to compute them.

The calculations of section III lead to waveform amplitudes that are accurate at $O(v/c)^{4}$ everywhere. However, we desire that our initial-data wave content already encode the phase as accurately as possible. Highly accurate phase for our initial data (via $h^{{\rm TT}}$), and hence in the leading edge of the waveforms we would extract from numerical evolution, is critical for parameter estimation following a detection.

For demonstrative purposes, in this section, we will restrict ourselves to the simplest phasing relations consistent with radiation-reaction inspiral as given by Eqs. (17-18), while using higher-order PN expressions than Eqs. (19 -20) for relating the orbit to the phase. For example, from Schäfer and Wex (1993), we have found to second PN (beyond leading) order:

[TABLE]

and we note that higher-order equivalents of these can be computed from Memmesheimer et al. (2004).

In the numerical construction of initial data, the primary input is the coordinate separation of the holes. In placing the punctures on the numerical grid, the separation must be maintained exactly. To ensure this, we invert Eq. (32) to obtain the exact $\Omega_{r}$ corresponding to our desired $r_{12}$. Then we use Eq. (18) with $t=0$ to find the coalescence time $t_{c}$ that yields this $\Omega_{r}$. Once we have obtained $t_{c}$, we then find the orbital phase $\Phi$ and frequency $\Omega$ at any source time $\tau$ directly from Eqs. (17-18), and the corresponding separation $r_{12}$ and momentum $p$ from Eqs. (32-33), or their higher-order equivalents.

In Fig. 3, we show a representative component of the retarded-time part of $h^{{\rm TT}}_{ij}$ for both circular and leading-order inspiral orbits. For both orbits, we use the extended Keplerian relations (32) and (33); otherwise the orbital configuration is that of Fig. 1. The coalescence time is now $t_{c}\sim 1100M$. We can see that the cumulative wavelength error of the circular-orbit assumption becomes very large at large distances from the sources.

This demonstrates that using inspiral orbits instead of circular orbits will significantly enhance the phase accuracy of the initial data, even though circular orbits are in principle sufficient when we include terms only up to $O(v/c)^{4}$ as done in this work. From now on we use only inspiral orbits.

Next, we compare our full waveform $h^{{\rm TT}}_{ij}$ (expressed as the combinations $h_{+}$ and $h_{\times}$) at an intermediate-field position ($r=100M$, $\theta=\pi/4$, $\phi=0$) to the lowest-order quadrupole result. In Fig. 4, the orbital configuration is the same as for Fig. 1.

As one can see, both the $+$ and $\times$ polarizations of our $h^{{\rm TT}}_{ij}$ agree very well with quadrupole results, as they should. We demonstrate the near- and intermediate-zone behavior of the new data on the initial spatial slice in Fig. 5. The quadrupole and full solutions agree very well outside $\sim 100M$. However, the full solution’s phase and amplitude approach the NZ solution closer to the sources.

IV.2 Numerical Implementation

After having confirmed that we have a PN three-metric $g_{ij}$ that is accurate up to errors of order $O(v/c)^{5}$, and that correctly approaches the quadrupole limit outside the near zone, we are now ready to construct initial data for numerical evolutions. In order to do so, we need the intrinsic curvature $K_{ij}$, which can be computed as in Tichy *et al. *​Tichy et al. (2003a) from the conjugate momentum. The difference is that here we use the full $\dot{h}^{{\rm TT}}_{ij}$ instead of the near-zone approximation $\dot{h}^{{\rm TT}\,(4)}_{ij}$ to obtain the conjugate momentum Schäfer (1985). The result is

[TABLE]

where the error term comes from neglecting terms like $h^{{\rm TT}}_{ij,(div)}$ at $O(v/c)^{5}$ in ${h}^{{\rm TT}}_{ij}$, and where $\psi_{PN}$, $\tilde{\pi}^{ij}_{(3)}$ and $\phi_{(2)}$ can be found in Tichy *et al. *​Tichy et al. (2003a). An additional difference is that the time derivative of $h^{{\rm TT}}_{ij}$ is evaluated numerically in this work. Note that the results for $g_{ij}$ are accurate up to $O(v/c)^{4}$, while the results for $K_{ij}$ are accurate up $O(v/c)^{5}$, because $K_{ij}$ contains an additional time derivative Tichy et al. (2003a); Yunes et al. (2006); Yunes and Tichy (2006).

Next we show the violations of the Hamiltonian and momentum constraints computed from $g_{ij}$ and $K_{ij}$, as functions of the binary separation $r_{12}$. As we can see in both panels of Fig. 6, the constraints become smaller for larger separations, because the post-Newtonian approximation gets better. Note that, as in Tichy et al. (2003a), the constraint violation remains finite everywhere, and is largest near each black hole.

IV.3 Curvature Invariants and Asymptotic Flatness

In analysis of both initial and evolved data, it is often instructive to investigate the behavior of scalar curvature invariants, as these give some idea of the far-field properties of our solution. We expect, for an asymptotically flat space-time, that in the far field, the speciality index $\mathcal{S}\equiv 27\mathcal{J}^{2}/\mathcal{I}^{3}$ will be close to unity. This can be seen from the following arguments. Let us choose a tetrad such that the Weyl tensor components $\psi_{1}$ and $\psi_{3}$ are both zero. Further, we assume that in the far field $\psi_{0}$ and $\psi_{4}$ are both perturbations of order $\epsilon$ off a Kerr background. Then

[TABLE]

which is indeed close to one. Note however, that this argument only works if the components of the Weyl tensor obey the peeling theorem, such that $\psi_{2}\sim O(r^{-3})$, $\psi_{0}\sim O(r^{-5})$ and $\psi_{4}\sim O(r^{-1})$. In particular, if $\psi_{0}$ falls off more slowly than $O(r^{-5})$, $S$ will grow for large $r$. Now observe that $\psi_{0}\sim O(r^{-5})\sim M^{3}/r^{5}$ is formally of $O(v/c)^{6}$. Thus, in order to see the expected behavior of $\mathcal{S}\approx 1$ in the far-field we need to go to $O(v/c)^{6}$. If we only go to $O(v/c)^{4}$ (as done in this work) $\psi_{0}$ consists of uncontrolled remainders only, which should in principle be dropped. When we numerically compute $\mathcal{S}$ we find that for our data, $\mathcal{S}$ deviates further and further from unity for large distances from the binary. This reflects the fact that the so-called “incoming” Weyl scalar $\psi_{0}$ only falls off as $1/r^{3}$, due to uncontrolled remainders at $O(v/c)^{6}$, which arise from a mixing of the background with the TT waveform.

V Discussion and Future Work

Exploring and validating PN inspiral waveforms is crucially important for gravitational-wave detection and for our theoretical understanding of black-hole binaries. Our goal has been to provide a step forward in this understanding by building a direct interface between the PN approach and numerical evolution, along the lines initially outlined in Ref. Tichy et al. (2003a). In this paper we have essentially completed the calculation of the transverse-traceless part of the ADM-TT metric to $O(v/c)^{4}$ provided in Tichy et al. (2003a), yielding data that, on the initial Cauchy slice, will describe the space-time into the far-field. We have incorporated this formulation into a numerical initial-data routine adapted to the “puncture” topology that has been so successful recently, and have explored these data’s numerical properties on the initial slice.

Our next step is to evolve these data with moving punctures, and investigate how the explicit incorporation of post-Newtonian waveforms in the initial data affects both the ensuing slow binary inspiral of the sources and the release of radiation from the system. We note especially that our data are non-conformally flat beyond $O(v/c)^{3}$. We expect our data to incorporate smaller unphysical initial distortions in the black holes than is possible with conformal flatness, and hence less spurious gravitational radiation during the numerical evolution. We see this as a very positive step toward providing further validation of numerical relativity results for multiple orbit simulations, since it permits comparison with PN results where they are expected to be reliable. Our initial data will also allow us to fully evaluate the validity of PN results for merging binaries by enabling comparison with the most accurate numerical relativity results.

We expect that further development of these data will certainly involve the use of more accurate orbital phasing information than the leading order given by Eqs. (17-18). This information is available in radiative coordinates (see, e.g. Eq. (6.29) of Blanchet (1996)) appropriate for far-field evaluation of the gravitational radiative modes; it may be possible to produce them in ADM-TT coordinates via a contact transformation, or by direct calculation (see, e.g. Damour et al. (2005)). For initial separations similar to the fiducial test case of this paper, $r_{12}\!=\!10M$, the order necessary for clean matching of the initial wave content with the new radiation generated in evolution should not be particularly high Baker et al. (2006c). As noted, the Keplerian relations Eqs. (32-33) can easily be extended to higher PN order.

The data presented already allow for arbitrary initial mass ratios $\nu$; this introduces the possibility of significant gravitational radiation in odd-$l$ multipoles, together with associated phenomena, such as in-plane recoil “kicks”. An interesting future development of these data will be the inclusion of spin angular momenta on the pre-merger holes. This will open our initial-data prescription to describing an even richer spectrum of binary radiation.

Acknowledgements.

We would like to thank L. Blanchet and G. Schäfer for generous assistance and helpful discussion. M.C., B.K. and B.W. gratefully acknowledge the support of the NASA Center for Gravitational Wave Astronomy (NAG5-13396). M.C. and B.K. also acknowledge the NSF for financial support under grants PHY-0354867 and PHY-0722315. B.K. also acknowledges support from the NASA Postdoctoral Program at the Oak Ridge Associated Universities. The work of W.T. was supported by NSF grant PHY-0555644. W.T. also acknowledges partial support from the NCSA under Grant PHY-060040T. The work of B.W. was also supported by NSF grants PHY-0245024 and PHY-0555484.

Appendix A Details of Integral Calculation

Here we present some more details of the calculations that lead to the three contributions to Eq. (23): Eqs. (27-29). Inserting Eq. (24) in the general integral (23), we can write $H^{ij}_{{\rm TT}\,A}[\vec{u}]$ as a combination of scalar and tensor terms:

[TABLE]

where the “$I$” integrals are defined as:

[TABLE]

Here $T\equiv t-\tau$, and $\vec{r}_{A}\equiv\vec{x}-\vec{x}_{A}$. We have also taken our integration coordinates such that $\vec{r}_{A}$ lies in the $z$ direction, so that the dummy momentum vector $\vec{k}$ satisfies

[TABLE]

Define the unit orthogonal vectors $\hat{n}_{A}\equiv(0,0,1)\;\;,\;\;\hat{\ell}\equiv(\cos\phi,\sin\phi,0)$. Then we can write

[TABLE]

We can also define a projector tensor onto $\hat{\ell}$:

[TABLE]

A.1 Angular integration

We will neglect the $A$ subscript for now, until it becomes relevant again. To calculate the integrals (37-39), we begin with the $\phi$ integration. The only $\phi$ dependence comes from the $\vec{\ell}$ parts of the $\vec{k}$ terms. It can be seen from elementary trigonometric integrals that:

[TABLE]

We use these to calculate the $\phi$ integrals for $I^{a\,b}_{A}$ and $I^{a\,b\,c\,d}_{A}$. Define $w\equiv\cos\theta$. Then

[TABLE]

So the next integrals will differ in their $\theta$ dependence, contained in the powers of $w$ above. The $\theta$ integrals will contain the following basic types:

[TABLE]

Now $I^{a\,b}$ and $I^{a\,b\,c\,d}$ can be written as the linear combinations:

[TABLE]

$I$ here can be expressed in terms of $g_{0}(a)$ above:

[TABLE]

The $1/2$ factor is because we moved to integrating $k$ over the whole real line instead of the positive half-line (this is permissible as $g_{n}(a)$ is an even function of $a$). $K$ and $L$ are defined analogously to $I$, but with extra even powers of $\cos\theta=w$:

[TABLE]

A.2 Momentum integration

Now we address the $k$ integrals, defined as:

[TABLE]

where we collect the positive exponents in the $g_{n}$ in the integrand of $f^{+}_{n}(k)$, and the negative exponents in $f^{-}_{n}(k)$:

[TABLE]

We calculate this as the sum of contour integrals of the “plus” and “minus” integrands (necessary, as the opposite signs require different contours). Each of these has poles at $k=0$, $k=k_{+}\equiv\omega+i\,\epsilon$, and $k=k_{-}\equiv-\omega-i\,\epsilon$ (the first of these is from the $g_{n}$). We integrate the “plus” integrands anticlockwise around the contour $C_{1}$, and the “minus” integrands anticlockwise around the contour $C_{2}$ (see Fig. 7); taking the limit $|k|\rightarrow\infty$, the contribution from the curved segments vanishes, and the residue theorem gives us:

[TABLE]

Calculating the residues, we find the values of each of the $J_{n}$:

[TABLE]

A.3 Frequency integration

Now we perform the $\omega$ integration. Inserting the results (51-53) into (47-49) respectively, we see that each of $I$, $K$ and $L$ contains a delta function, which we can extract:

[TABLE]

where the new terms on the right-hand side come from the $J_{n}$ above, grouped by exponential, as that is what determines the contours chosen during integration (see Fig. 7):

[TABLE]

Now the residues are as follows (taking the $\epsilon\rightarrow 0$ limit):

[TABLE]

The only pole is at $\omega=-i\,\epsilon$, so if we can close the contour in the upper half-plane, we will get zero.

• •

For $T<0$, both the “a” and “b” integrals can be closed in $C_{1}$. Result: zero contribution.

• •

For $0<T<r$, the “a” integrals can be closed in $C_{1}$, but the “b” integrals must be closed in $C_{2}$. Result: “b” contribution.

• •

For $T>r$, both the “a” and “b” integrals must be closed in $C_{2}$. But then the “a” and “b” residues cancel out. Result: zero contribution.

Thus the only interesting contribution happens in the interval $0<T<r\Leftrightarrow t-r(\tau)<\tau<t$. In this case, the final integrals yield

[TABLE]

leading to the final result for $K$ and $L$:

[TABLE]

We use these to calculate the $I^{i\,j}$ and $I^{i\,j\,k\,l}$:

[TABLE]

A.4 Time integration

The final integrations will be over the source time $\tau$. The “crossing times” for the two $\Theta$ functions are $\tau=t$ and $\tau=t^{r}$, where $t$ is the present field time, and $t^{r}$ the corresponding retarded time defined by (25). Now taking a general function $y(\tau)$, we find that

[TABLE]

These can now be substituted into the general integral (36). We write the result as a sum of terms at the present field-point time $t$, the retarded time $t_{A}^{r}$, and interval terms between them,

[TABLE]

[TABLE]

[TABLE]

[TABLE]