text:
# Comparing initial-data sets for binary black holes

## I Introduction

Numerical evolutions of black holes have been improved slowly but steadily over the last few years and now first attempts are being made to extract physical information from these evolutions. Most notably one wants to predict the gravitational radiation emitted during black hole coalescence Baker et al. (2001); Alcubierre et al. (2001); Baker et al. (2002).

The quality of the initial data will be crucial to the success of the predictions of the gravitational wave forms. Unphysical gravitational radiation present in the initial data will contribute to the gravitational waves computed in an evolution and might overwhelm the true gravitational wave signature of the physical process under consideration. Therefore an important question is how to control the gravitational wave content of initial-data sets, and how to specify astrophysically relevant initial data with the appropriate gravitational wave content, for e.g. two black holes orbiting each other. Unfortunately, assessing and controlling the gravitational wave content of initial-data sets is not well understood at all.

The mere construction of an initial-data set alone is fairly involved, since every initial-data set must satisfy a rather complicated set of four partial differential equations, the so-called constraint equations of general relativity. The question of how to solve these equations, and how to specify initial data representing binary black holes in particular, has received considerable attention.

We consider in this paper three different approaches that transform the constraint equations into elliptic equations: The conformal transverse-traceless (TT) decompositionYork, Jr. (1979), the physical TT decomposition Murchadha and York, Jr. (1974a, b); Murchadha and York, Jr. (1976) and the conformal thin sandwich decompositionYork, Jr. (1999). These decompositions split the variables on the initial-data surface into various pieces in such a way that the constraint equations determine some of the pieces while not restricting the others. After these freely specifiable pieces are chosen, the constraint equations are solved and the results are combined with the freely specifiable pieces to yield a valid initial-data set.

Any reasonable choice for the freely specifiable pieces will lead to a valid initial-data set. Furthermore, any one of these decompositions can generate any desired initial-data set, given the correct choices of the freely specifiable pieces. However, it is not clear what choices of freely specifiable pieces lead to initial-data sets with the desired properties.

The decompositions we consider here lead to four coupled nonlinear elliptic partial differential equations. Since such equations are difficult to solve, the early approach to constructing initial data was pragmatic: One used the conformal TT decomposition with additional restrictions on the freely specifiable pieces, most notably conformal flatness and maximal slicing. These assumptions decouple the constraints and allow for analytical solutions to the momentum constraints, the so-called Bowen-York extrinsic curvatureBowen (1979); Bowen and York, Jr. (1980); Kulkarni et al. (1983). All that remains is to solve a single elliptic equation, the Hamiltonian constraint. This approach has been used in several variationsThornburg (1987); Cook et al. (1993); Brandt and Brügmann (1997).

However, these numerical simplifications come at a cost. The freely specifiable pieces have been restricted to a small subset of all possible choices. One therefore can generate only a subset of all possible initial-data sets, one that might not contain the desired astrophysically relevant initial-data sets.

Over the last few years there have been additional developments: Post-Newtonian results have indicated that binary black hole metrics are not conformally flatRieth (1997); Damour et al. (2000). With certain restrictions on the slicing, it has also been shown that a single stationary spinning black hole cannot be represented with a conformally flat spatial metric Monroe (1976); Garat and Price (2000). In Pfeiffer et al. (2000), it was shown that conformally flat initial data sets for spinning binary black holes contain an unphysical contamination. Moreover, computations in spherical symmetryLousto and Price (1998) indicated that initial-data sets depend strongly on the choice of the extrinsic curvature and that the use of the Bowen-York extrinsic curvature might be problematic.

Therefore it is necessary to move beyond conformally flat initial data and to explore different choices for the extrinsic curvature. Matzner et alMatzner et al. (1999) proposed a non-flat conformal metric based on the superposition of two Kerr-Schild metrics; a solution based on this proposal was obtained in Marronetti and Matzner (2000). This work demonstrated the existence of solutions to the 3D set of equations, but did not examine the data sets in any detail. Refs. Gourgoulhon et al. (2002); Grandclément et al. (2002) obtained solutions to a similar set of equations during the computation of quasi-circular orbits of binary black holes. However, these works assumed conformal flatness.

In this paper we present a code capable of solving the three above-mentioned decompositions of the constraint equations for arbitrary choices of the freely specifiable pieces. This code is based on spectral methods which have been used successfully for several astrophysical problems (see e.g. Bonazzola et al. (1993, 1998); Kidder et al. (2000); Gourgoulhon et al. (2001); Kidder et al. (2001); Ansorg et al. (2002); Grandclément et al. (2002)). Our code is described in detail in a separate paperPfeiffer et al. (2002).

We compute solutions of the different decompositions for the non-flat conformal metric proposed in Ref. Matzner et al. (1999). Each decomposition has certain choices for the freely specifiable pieces and boundary conditions that seem “natural” and which we use in our solutions. We compare the computed initial-data sets with each other and with the “standard” conformally-flat solution using the Bowen-York extrinsic curvature. Our major results confirm that

1. the different decompositions generate different physical initial-data sets for seemingly similar choices for the freely specifiable pieces.
2. the choice of extrinsic curvature is critical.

The first result is certainly not unexpected, but each of these factors can cause relative differences of several per cent in gauge-invariant quantities like the ADM-energy.

We also find that the conformal TT/physical TT decompositions generate initial-data sets with ADM-energies $`23`$% higher than data sets of the conformal thin sandwich decomposition. We demonstrate that this higher ADM-energy is related to the choice of the freely specifiable part of the extrinsic curvature. In addition, we find that the solutions depend significantly on the boundary conditions used.

The paper is organized as follows. In the next section we describe the three decompositions. Section III explains how we choose the freely specifiable data within each decomposition. In section IV we describe and test our elliptic solver. Section V presents our results, which we discuss in section VI.

## II Decompositions of Einstein’s equations and the constraint equations

### II.1 3+1 Decomposition

In this paper we use the standard 3+1 decomposition of Einstein’s equations. We foliate the spacetime with $`t=\text{const}`$ hypersurfaces and write the four-dimensional metric as

$${}_{}{}^{(4)}ds^2=N^2dt^2+\gamma _{ij}(dx^i+N^idt)(dx^j+N^jdt),$$
(1)
where $`\gamma _{ij}`$ represents the induced 3-metric on the hypersurfaces, and $`N`$ and $`N^i`$ represent the lapse function and the shift vector, respectively. We define the extrinsic curvature $`K_{ij}`$ on the slice by

$$𝐊=\frac{1}{2}_n{}_{}{}^{(4)}𝐠$$
(2)
where $`{}_{}{}^{(4)}𝐠`$ is the space-time metric, $`n`$ the unit normal to the hypersurface, and $``$ denotes the projection operator into the $`t=\text{const}`$ slice. Einstein’s equations divide into constraint equations, which constrain the data $`(\gamma _{ij},K^{ij})`$ on each hypersurface, and into evolution equations, which determine how the data $`(\gamma _{ij},K^{ij})`$ evolve from one hypersurface to the next. The constraint equations are

$`R+K^2K_{ij}K^{ij}`$ $`=16\pi G\rho `$ (3)
$`_j\left(K^{ij}\gamma ^{ij}K\right)`$ $`=8\pi Gj^i.`$ (4)
Eq. (3) is called the Hamiltonian constraint, and Eq. (4) is referred to as the momentum constraint. $`K=\gamma _{ij}K^{ij}`$ is the trace of the extrinsic curvature, $``$ and $`R`$ denote the three dimensional covariant derivative operator and the Ricci scalar compatible with $`\gamma _{ij}`$. $`\rho `$ and $`j^i`$ are the energy and momentum density, respectively. Both vanish for the vacuum spacetimes considered here.

The evolution equation for $`\gamma _{ij}`$ is

$$_t\gamma _{ij}=2NK_{ij}+_iN_j+_jN_i,$$
(5)
which follows from Eq. (2). There is a similar albeit longer equation for $`_tK_{ij}`$ which we will not need in this paper. The choices of $`N`$ and $`N^i`$ are arbitrary. One can in principle use any lapse and shift in the evolution off the initial-data surface, although some choices of lapse and shift are better suited to numerical implementation than others.

Later in this paper we will often refer to the trace-free piece of Eq. (5). Denote the tracefree piece of a tensor by $`\text{TF}(.)`$, and define $`\gamma det\gamma _{ij}`$. From Eq. (5) and the fact that $`\delta \mathrm{ln}\gamma =\gamma ^{kl}\delta \gamma _{kl}`$, it follows that

$$\text{TF}(_t\gamma _{ij})=\gamma ^{1/3}_t\left(\gamma ^{1/3}\gamma _{ij}\right)=2NA_{ij}+(𝕃N)_{ij}.$$
(6)
Here $`A_{ij}=K_{ij}\frac{1}{3}\gamma _{ij}K`$ denotes the trace-free extrinsic curvature, and

$$(𝕃N)^{ij}^iN^j+^jN^i\frac{2}{3}\gamma ^{ij}_kN^k.$$
(7)
$`𝕃`$ always acts on a vector, so the ’N’ in $`(𝕃N)^{ij}`$ denotes the shift vector $`N^i`$ and not the lapse $`N`$.

### II.2 Decomposition of the constraint equations

Equations (3) and (4) constrain four degrees of freedom of the 12 quantities $`(\gamma _{ij},K^{ij})`$. However, it is not immediately clear which pieces of $`\gamma _{ij}`$ and $`K^{ij}`$ are constrained and which pieces can be chosen at will. Several decompositions have been developed to divide the 12 degrees of freedom into freely specifiable and constrained pieces. We will now review some properties of the three decompositions we consider in this paper.

All three decompositions follow the York-Lichnerowicz approach and use a conformal transformation on the physical 3-metric $`\gamma _{ij}`$,

$$\gamma _{ij}=\psi ^4\stackrel{~}{\gamma }_{ij}.$$
(8)
$`\psi `$ is called the conformal factor, $`\stackrel{~}{\gamma }_{ij}`$ the background metric or conformal metric. We will denote all conformal quantities with a tilde. In particular, $`\stackrel{~}{}`$ is the covariant derivative operator associated with $`\stackrel{~}{\gamma }_{ij}`$, and $`\stackrel{~}{R}_{ij}`$ and $`\stackrel{~}{R}`$ are the Ricci tensor and Ricci scalar of $`\stackrel{~}{\gamma }_{ij}`$.

The extrinsic curvature is split into its trace and trace-free part,

$$K^{ij}=A^{ij}+\frac{1}{3}\gamma ^{ij}K.$$
(9)
The three decompositions of the constraint equations we discuss in this paper differ in how $`A^{ij}`$ is decomposed. For each decomposition, we discuss next the relevant equations, and describe how we choose the quantities one has to specify before solving the equations. We use the conventions of Cook (2000).

#### II.2.1 Conformal TT Decomposition

In this decomposition one first conformally transforms the traceless extrinsic curvature,

$$A^{ij}=\psi ^{10}\stackrel{~}{A}^{ij},$$
(10)
and then applies a TT decomposition with respect to the background metric $`\stackrel{~}{\gamma }_{ij}`$:

$$\stackrel{~}{A}^{ij}=\stackrel{~}{A}_{TT}^{ij}+(\stackrel{~}{𝕃}X)^{ij}.$$
(11)
The operator $`\stackrel{~}{𝕃}`$ is defined by Eq. (7) but using the conformal metric $`\stackrel{~}{\gamma }_{ij}`$ and derivatives associated with $`\stackrel{~}{\gamma }_{ij}`$. $`\stackrel{~}{A}_{TT}^{ij}`$ is transverse with respect to the conformal metric, $`\stackrel{~}{}_j\stackrel{~}{A}_{TT}^{ij}=0`$, and is traceless.

Substituting Eqs. (10) and (11) into the momentum constraint (4), one finds that it reduces to an elliptic equation for $`X^i`$, whereas $`\stackrel{~}{A}_{TT}^{ij}`$ is unconstrained.

In order to specify the transverse-traceless tensor $`\stackrel{~}{A}_{TT}^{ij}`$ one usually has to construct it from a general symmetric trace-free tensor $`\stackrel{~}{M}^{ij}`$ by subtracting the longitudinal piece. As described in Cook (2000) one can incorporate the construction of $`\stackrel{~}{A}_{TT}^{ij}`$ from $`\stackrel{~}{M}^{ij}`$ into the momentum constraint, arriving at the following equations:

$$\stackrel{~}{}^2\psi \frac{1}{8}\psi \stackrel{~}{R}\frac{1}{12}\psi ^5K^2+\frac{1}{8}\psi ^7\stackrel{~}{A}_{ij}\stackrel{~}{A}^{ij}=2\pi G\psi ^5\rho ,$$
(12)
$$\stackrel{~}{\mathrm{\Delta }}_𝕃V^i\frac{2}{3}\psi ^6\stackrel{~}{}^iK+\stackrel{~}{}_j\stackrel{~}{M}^{ij}=8\pi G\psi ^{10}j^i,$$
(13)
where $`\stackrel{~}{A}^{ij}`$ and the operator $`\stackrel{~}{\mathrm{\Delta }}_𝕃`$ are defined by

$$\stackrel{~}{A}^{ij}=(\stackrel{~}{𝕃}V)^{ij}+\stackrel{~}{M}^{ij}$$
(14)
and

$$\stackrel{~}{\mathrm{\Delta }}_𝕃V^i\stackrel{~}{}_j(\stackrel{~}{𝕃}V)^{ij}.$$
(15)
After solving these equations for $`\psi `$ and $`V^i`$, one obtains the physical metric $`\gamma _{ij}`$ from (8) and the extrinsic curvature from

$$K^{ij}=\psi ^{10}\stackrel{~}{A}^{ij}+\frac{1}{3}\psi ^4\stackrel{~}{\gamma }^{ij}K.$$
(16)
We will refer to Eqs. (12) and (13) together with (14), (16) and (8) as the conformal TT equations. In these equations we are free to specify the background metric $`\stackrel{~}{\gamma }_{ij}`$, the trace of the extrinsic curvature $`K`$, and a symmetric traceless tensor $`\stackrel{~}{M}^{ij}`$. The solution $`V^i`$ will contain a contribution that removes the longitudinal piece from $`\stackrel{~}{M}^{ij}`$ and the piece that solves the momentum constraint if $`\stackrel{~}{M}^{ij}`$ were transverse-traceless.

This decomposition has been the most important in the past, since if one chooses a constant $`K`$ and if one considers vacuum spacetimes then the momentum constraint (13) decouples from the Hamiltonian constraint (12). Moreover, if one assumes conformal flatness and $`\stackrel{~}{M}^{ij}=0`$, it is possible to write down analytic solutions to Eq. (13), the so-called Bowen-York extrinsic curvature. In that case one has to deal with only one elliptic equation for $`\psi `$. The Bowen-York extrinsic curvature can represent multiple black holes with arbitrary momenta and spins. One can fix boundary conditions for $`\psi `$ by requiring that the initial-data slice be inversion symmetric at both throatsMisner (1963); Cook (1991). In that case one has to modify the extrinsic curvature using a method of images. We will include initial-data sets obtained with this approach below, where we will refer to them as inversion symmetric initial data.

Reasonable choices for the freely specifiable pieces $`\stackrel{~}{\gamma }_{ij}`$, $`K`$, $`\stackrel{~}{M}^{ij}`$ will lead to an initial-data set $`(\gamma _{ij},K^{ij})`$ that satisfies the constraint equations. How should we choose all these functions in order to obtain a desired physical configuration, say a binary black hole with given linear momenta and spins for the individual holes? We can gain insight into this question by considering how the conformal TT decompositions can recover a known solution.

Suppose we have a known solution $`(\gamma _{0ij},K_0^{ij})`$ of the constraint equations. Denote the trace and trace-free parts of this extrinsic curvature by $`K_0`$ and $`A_0^{ij}`$, respectively. If we set

$$\stackrel{~}{\gamma }_{ij}=\gamma _{0ij},K=K_0,\stackrel{~}{M}^{ij}=A_0^{ij}$$
(17)
then

$$\psi =1,V^i=0$$
(18)
trivially solve Eqs. (12-13). Note that we have to set $`\stackrel{~}{M}^{ij}`$ equal to the trace-free part of the extrinsic curvature.

Now suppose we have a guess for a metric and an extrinsic curvature, which —most likely— will not satisfy the constraint equations (3) and (4). Set $`\stackrel{~}{\gamma }_{ij}`$ to the guess for the metric, and set $`K`$ and $`\stackrel{~}{M}^{ij}`$ to the trace and trace-free piece of the guess of the extrinsic curvature. By solving the conformal TT equations we can compute $`(\gamma _{ij},K^{ij})`$ that satisfy the constraint equations. If our initial guess is “close” to a true solution, we will have $`\psi 1`$ and $`V^i0`$, so that $`\gamma _{ij}`$ and $`K^{ij}`$ will be close to our initial guess.

Thus one can guess a metric and extrinsic curvature as well as possible and then solve the conformal TT equations to obtain corrected quantities that satisfy the constraint equations.

An artifact of the conformal TT decomposition is that one has no direct handle on the transverse traceless piece with respect to the physical metric. For any vector $`X^i`$,

$$(𝕃X)^{ij}=\psi ^4(\stackrel{~}{𝕃}X)^{ij}.$$
(19)
Thus, Eqs. (10) and (11) imply

$$A^{ij}=\psi ^{10}\stackrel{~}{A}_{TT}^{ij}+\psi ^6(𝕃X)^{ij}.$$
(20)
For any symmetric traceless tensor $`S^{ij}`$

$$_jS^{ij}=\psi ^{10}\stackrel{~}{}_j\left(\psi ^{10}S^{ij}\right).$$
(21)
Therefore the first term on the right hand side of Eq. (20) is transverse-traceless with respect to the physical metric,

$$_j\left(\psi ^{10}\stackrel{~}{A}_{TT}^{ij}\right)=0.$$
(22)
However, the second term on the right hand side of Eq. (20) is conformally weighted. Therefore, Eq. (20) does not represent the usual TT decomposition.

#### II.2.2 Physical TT Decomposition

In this case one decomposes the physical traceless extrinsic curvature directly:

$$A^{ij}=A_{TT}^{ij}+(𝕃X)^{ij}.$$
(23)
As above in the conformal TT decomposition, the momentum constraint becomes an elliptic equation for $`X^i`$. We can again incorporate the construction of the symmetric transverse traceless tensor $`A_{TT}^{ij}`$ from a general symmetric tensor $`\stackrel{~}{M}^{ij}`$ into the momentum constraint. Then one obtains the physical TT equations:

$$\stackrel{~}{}^2\psi \frac{1}{8}\psi \stackrel{~}{R}\frac{1}{12}\psi ^5K^2+\frac{1}{8}\psi ^5\stackrel{~}{A}_{ij}\stackrel{~}{A}^{ij}=2\pi G\psi ^5\rho ,$$
(24)
$$\stackrel{~}{\mathrm{\Delta }}_𝕃V^i+6(\stackrel{~}{𝕃}V)^{ij}\stackrel{~}{}_j\mathrm{ln}\psi \frac{2}{3}\stackrel{~}{}^iK+\psi ^6\stackrel{~}{}_j\stackrel{~}{M}^{ij}=8\pi G\psi ^4j^i,$$
(25)
where $`\stackrel{~}{A}^{ij}`$ is defined by

$$\stackrel{~}{A}^{ij}=(\stackrel{~}{𝕃}V)^{ij}+\psi ^6\stackrel{~}{M}^{ij}.$$
(26)
When we have solved (24) and (25) for $`\psi `$ and $`V^i`$, the physical metric is given by (8), and the extrinsic curvature is

$$K^{ij}=\psi ^4\left(\stackrel{~}{A}^{ij}+\frac{1}{3}\stackrel{~}{\gamma }^{ij}K\right).$$
(27)
We are free to specify the background metric $`\stackrel{~}{\gamma }_{ij}`$, the trace of the extrinsic curvature $`K`$, and a symmetric traceless tensor $`\stackrel{~}{M}^{ij}`$. As with the conformal TT equations, the solution $`V^i`$ will contain a contribution that removes the longitudinal piece from $`\stackrel{~}{M}^{ij}`$ and a piece that solves the momentum constraint if $`\stackrel{~}{M}^{ij}`$ were transverse-traceless.

These equations can be used in the same way as the conformal TT equations. Guess a metric and extrinsic curvature, set $`\stackrel{~}{\gamma }_{ij}`$ to the guess for the metric, and $`K`$ and $`\stackrel{~}{M}^{ij}`$ to the trace and trace-free pieces of the guess for the extrinsic curvature. Then solve the physical TT equations to obtain a corrected metric $`\gamma _{ij}`$ and a corrected extrinsic curvature $`K^{ij}`$ that satisfy the constraint equations.

The transverse traceless piece of $`K^{ij}`$ (with respect to $`\gamma _{ij}`$) will be the transverse traceless piece of $`\psi ^{10}\stackrel{~}{M}^{ij}`$. One can also easily rewrite the physical TT equations such that $`\psi ^{10}\stackrel{~}{M}^{ij}`$ can be freely chosen instead of $`\stackrel{~}{M}^{ij}`$. So, in this decomposition we can directly control the TT piece of the physical extrinsic curvature. We have chosen to follow Cook (2000) since it seems somewhat more natural to specify two conformal quantities, $`\stackrel{~}{\gamma }_{ij}`$ and $`\stackrel{~}{M}^{ij}`$ than to specify one conformal and one physical quantity.

#### II.2.3 Conformal thin sandwich decomposition

The conformal and physical TT decompositions rely on a tensor splitting to decompose the trace-free part of the extrinsic curvature. In contrast, the conformal thin sandwich decomposition simply defines $`A^{ij}`$ by Eq. (10) and the decomposition

$$\stackrel{~}{A}^{ij}\frac{1}{2\stackrel{~}{\alpha }}\left((\stackrel{~}{𝕃}\beta )^{ij}\stackrel{~}{u}^{ij}\right),$$
(28)
where $`\stackrel{~}{u}^{ij}`$ is symmetric and tracefree. Eq. (28) is motivated by Eq. (6): If one evolves an initial-data set with $`A^{ij}`$ of the form (28) using as lapse and shift

$`N`$ $`=\psi ^6\stackrel{~}{\alpha },`$ (29)
$`N^i`$ $`=\beta ^i,`$
then

$$\text{TF}(_t\gamma _{ij})=\psi ^4\stackrel{~}{u}_{ij}.$$
(30)
Therefore, the decomposition (28) is closely related to the kinematical quantities in an evolution. Although $`\stackrel{~}{\alpha }`$ and $`\beta ^i`$ are introduced in the context of initial data, one usually refers to them as the “conformal lapse” and “shift”. While the form of Eq. (28) is similar in form to the conformal and physical TT decompositions, there are differences. In particular, $`\stackrel{~}{u}^{ij}`$ is not divergenceless.

Within the conformal thin sandwich decomposition, the constraint equations take the form:

$$\stackrel{~}{}^2\psi \frac{1}{8}\psi \stackrel{~}{R}\frac{1}{12}\psi ^5K^2+\frac{1}{8}\psi ^7\stackrel{~}{A}_{ij}\stackrel{~}{A}^{ij}=2\pi G\psi ^5\rho $$
(31)
$$\stackrel{~}{\mathrm{\Delta }}_𝕃\beta ^i(\stackrel{~}{𝕃}\beta )^{ij}\stackrel{~}{}_j\mathrm{ln}\stackrel{~}{\alpha }\frac{4}{3}\stackrel{~}{\alpha }\psi ^6\stackrel{~}{}^iK$$

$$\stackrel{~}{\alpha }\stackrel{~}{}_j\left(\frac{1}{\stackrel{~}{\alpha }}\stackrel{~}{u}^{ij}\right)=16\pi G\stackrel{~}{\alpha }\psi ^{10}j^i$$
(32)
Having solved Eqs. (31) and (32) for $`\psi `$ and the vector $`\beta ^i`$, one obtains the physical metric from (8) and the extrinsic curvature from

$$K^{ij}=\psi ^{10}\stackrel{~}{A}^{ij}+\frac{1}{3}\psi ^4\stackrel{~}{\gamma }^{ij}K.$$
(33)
In this decomposition we are free to specify a conformal metric $`\stackrel{~}{\gamma }_{ij}`$, the trace of the extrinsic curvature $`K`$, a symmetric trace-free tensor $`\stackrel{~}{u}^{ij}`$ and a function $`\stackrel{~}{\alpha }`$.

It seems that the conformal thin sandwich decomposition contains additional degrees of freedom in the form of the function $`\stackrel{~}{\alpha }`$ and three additional unconstrained components of $`\stackrel{~}{u}^{ij}`$. This is not the case. The longitudinal piece of $`\stackrel{~}{u}^{ij}`$ corresponds to the gauge choice of the actual shift vector used in an evolution. Thus $`\stackrel{~}{u}^{ij}`$ really only contributes two degrees of freedom, just like $`\stackrel{~}{M}^{ij}`$ in the conformal and physical TT decompositions. Furthermore, we can reach any reasonable physical solution $`(\gamma _{ij},K^{ij})`$ with any reasonable choice of $`\stackrel{~}{\alpha }`$; each choice of $`\stackrel{~}{\alpha }`$ simply defines a new decomposition. A forthcoming article by YorkYork, Jr. (2002) will elaborate on these ideas. Note that for $`\stackrel{~}{\alpha }=1/2`$ we recover the conformal TT decomposition.

Let us now turn to the question of how one should pick the freely specifiable data in the conformal thin sandwich approach. We motivate our prescription again by considering how to recover a known spacetime: Assume we are given a full four-dimensional spacetime with 3+1 quantities $`\gamma _{0ij}`$, $`K_0^{ij}`$, $`N_0^i`$ and $`N_0`$. Further assume the spacetime is stationary and the slicing is such that $`_t\gamma _{ij}=_tK_{ij}=0`$. An example for such a situation is a Kerr black hole in Kerr-Schild or Boyer-Lindquist coordinates.

Using $`_t\gamma _{0ij}=0`$ in Eq. (6) yields a relation for the trace-free extrinsic curvature

$$A_0^{ij}=\frac{1}{2N_0}(𝕃N_0)^{ij}.$$
(34)
This is a decomposition of the form (28) with $`\stackrel{~}{u}^{ij}=0`$. Therefore, if we choose the freely specifiable data for the conformal thin sandwich equations as

$`\stackrel{~}{\gamma }_{ij}`$ $`=\gamma _{0ij},`$ $`\stackrel{~}{\alpha }`$ $`=N_0,`$ (35)
$`K`$ $`=K_0,`$ $`\stackrel{~}{u}^{ij}`$ $`=0,`$
and if we use appropriate boundary conditions, then the solution of the conformal thin sandwich equations will be $`\psi =1`$ and $`\beta ^i=N_0^i`$. As part of the solution, we obtain the shift vector needed for an evolution to produce $`\text{TF}(_t\gamma _{ij})=0`$. Not needing a guess for the trace-free extrinsic curvature, and having the solution $`\beta ^i`$ automatically provide an initial shift for evolution, make the conformal thin sandwich equations very attractive.

In order to generate initial-data slices that permit an evolution with zero time derivative of the conformal metric — a highly desirable feature for quasi-equilibrium data, or for a situation with holes momentarily at rest — one can proceed as follows: Set $`\stackrel{~}{\gamma }_{ij}`$ and $`K`$ to the guess for the metric and trace of extrinsic curvature, respectively. Set $`\stackrel{~}{\alpha }`$ to the lapse function that one is going to use in the evolution, and set $`\stackrel{~}{u}^{ij}=0`$. If these guesses are good, the conformal factor $`\psi `$ will be close to 1, and $`N=\psi ^6\stackrel{~}{\alpha }`$ as well as $`N^i=\beta ^i`$ give us the actual lapse function and shift vector to use in the evolution.

## III Choices for the freely specifiable data

### III.1 Kerr-Schild coordinates

We base our choice for the freely specifiable data on a superposition of two Kerr black holes in Kerr-Schild coordinates. In this section we describe the Kerr-Schild solution and collect necessary equations. We also describe how we compute the 3-metric, lapse, shift and extrinsic curvature for a boosted black hole with arbitrary spin.

A Kerr-Schild metric is given by

$$g_{\mu \nu }=\eta _{\mu \nu }+2Hl_\mu l_\nu ,$$
(36)
where $`\eta _{\mu \nu }`$ is the Minkowski metric, and $`l_\mu `$ is a null-vector with respect to both the full metric and the Minkowski metric: $`g^{\mu \nu }l_\mu l_\nu =\eta ^{\mu \nu }l_\mu l_\nu =0`$. The 3-metric, lapse and shift are

$`\gamma _{ij}`$ $`=\delta _{ij}+2Hl_il_j,`$ (37)
$`N`$ $`=(1+2Hl^tl^t)^{1/2},`$ (38)
$`N^i`$ $`={\displaystyle \frac{2Hl^tl^i}{1+2Hl^tl^t}}.`$ (39)
For a black hole at rest at the origin with mass $`M`$ and angular momentum $`M\stackrel{}{a}`$, one has

$`H`$ $`={\displaystyle \frac{Mr^3}{r^4+(\stackrel{}{a}\stackrel{}{x})^2}},`$ (40)
$`l_\mu ^{\text{rest}}`$ $`=(1,\stackrel{}{l}_{\text{rest}}),`$ (41)
$`\stackrel{}{l}_{\text{rest}}`$ $`={\displaystyle \frac{r\stackrel{}{x}\stackrel{}{a}\times \stackrel{}{x}+(\stackrel{}{a}\stackrel{}{x})\stackrel{}{a}/r}{r^2+a^2}},`$ (42)
with

$$r^2=\frac{\stackrel{}{x}^2\stackrel{}{a}^2}{2}+\left(\frac{(\stackrel{}{x}^2\stackrel{}{a}^2)^2}{4}+(\stackrel{}{a}\stackrel{}{x})^2\right)^{1/2}.$$
(43)
For a nonrotating black hole with $`\stackrel{}{a}=0`$, $`H`$ has a pole at the origin, whereas for rotating black holes, $`r`$ has a ring singularity. We will therefore have to excise from the computational domain a region close to the center of the Kerr-Schild black hole.

Under a boost, a Kerr-Schild coordinate system transforms into a Kerr-Schild coordinate system. Applying a Lorentz transformation with boost velocity $`v^i`$ to $`l_\mu ^{\text{rest}}`$, we obtain the null-vector $`l_\mu `$ of the boosted Kerr-Schild coordinate system. Eqs. (37-39) give then the boosted 3-metric, lapse, and shift. Since all time-dependence is in the uniform motion, evolution with lapse $`N`$ and shift $`N^i`$ yields $`_t\gamma _{ij}=v^k_k\gamma _{ij}`$, and from Eq. (5) one can compute the extrinsic curvature

$$K_{ij}=\frac{1}{2N}\left(v^k_k\gamma _{ij}+_iN_j+_jN_i\right).$$
(44)
If this initial-data set is evolved with the shift $`N^i`$, the black hole will move through the coordinate space with velocity $`v^i`$. However, if the evolution uses the shift vector $`N^i+v^i`$, the coordinates will move with the black hole, and the hole will be at rest in coordinate space. The spacetime is nonetheless different from a Kerr black hole at rest. The ADM-momentum will be $`P_{ADM}^i=\gamma Mv^i`$, where $`M`$ is the rest-mass of the hole and $`\gamma =(1\stackrel{}{v}^2)^{1/2}`$.

### III.2 Freely specifiable pieces

We want to generate initial data for a spacetime containing two black holes with masses $`M_{A,B}`$, velocities $`\stackrel{}{v}_{A,B}`$ and spins $`M_A\stackrel{}{a}_A`$ and $`M_B\stackrel{}{a}_B`$.

We follow the proposal of Matzner et al Matzner et al. (1999); Marronetti and Matzner (2000) and base our choices for the freely specifiable choices on two Kerr-Schild coordinate systems describing two individual black holes. The first black hole with label A has an associated Kerr-Schild coordinate system with metric

$$\gamma _{Aij}=\delta _{ij}+2H_Al_{Ai}l_{Aj},$$
(45)
and with an extrinsic curvature $`K_{Aij}`$, a lapse $`N_A`$ and a shift $`N_A^i`$. The trace of the extrinsic curvature is $`K_A`$. All these quantities can be computed as described in the previous section, III.1. The second black hole has a similar set of associated quantities which are labeled with the letter B.

For all three decompositions, we need to choose a conformal metric and the trace of the extrinsic curvature. We choose

$$\stackrel{~}{\gamma }_{ij}=\delta _{ij}+2H_Al_{Ai}l_{Aj}+2H_Bl_{Bi}l_{Bj}$$
(46)
$$K=K_A+K_B$$
(47)
The metric is singular at the center of each hole. Therefore we have to excise spheres around the center of each hole from the computational domain. We now specify for each decomposition the remaining freely specifiable pieces and boundary conditions.

#### III.2.1 Conformal TT and physical TT decompositions

For the conformal TT and physical TT decompositions we will be solving for a correction to our guesses. As guess for the trace-free extrinsic curvature, we use a superposition

$$\stackrel{~}{M}^{ij}=\left(K_{Ak}^{(i}+K_{Bk}^{(i}\frac{1}{3}\delta _k^{(i}(K_A+K_B)\right)\stackrel{~}{\gamma }^{j)k}.$$
(48)
$`\stackrel{~}{M}^{ij}`$ is symmetric and trace-free with respect to the conformal metric, $`\stackrel{~}{\gamma }_{ij}\stackrel{~}{M}^{ij}=0`$. Solving for a correction only, we expect that $`\psi 1`$ and $`V^i0`$, hence we use Dirichlet boundary conditions

$$\psi =1,V^i=0.$$
(49)
#### III.2.2 Conformal thin sandwich

For conformal thin sandwich, we restrict the discussion to either two black holes at rest, or in a quasi-circular orbit in corotating coordinates. In these cases, one expects small or even vanishing time-derivatives, $`_t\gamma _{ij}0`$, and so Eq. (30) yields the simple choice

$$\stackrel{~}{u}^{ij}=0.$$
(50)
The conformal 3-metric and the trace of the extrinsic curvature are still given by Eqs. (46) and (47). Orbiting black holes in a corotating frame will not move in coordinate space, therefore we do not boost the individual Kerr-Schild metrics in this decomposition: $`v_{A/B}^i=0`$. The lapse functions $`N_{A/B}`$ and the shifts $`N_{A/B}^i`$ are also for unboosted Kerr-Schild black holes.

We use Dirichlet boundary conditions:

$`\psi `$ $`=1`$ all boundaries (51a)
$`\beta ^i`$ $`=N_A^i`$ sphere inside hole A (51b)
$`\beta ^i`$ $`=N_B^i`$ sphere inside hole B (51c)
$`\beta ^i`$ $`=\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{r}`$ outer boundary (51d)
Eq. (51d) ensures that we are in a corotating reference frame; the cross-product is performed in flat space, and $`\stackrel{}{\mathrm{\Omega }}=0`$ corresponds to two black holes at rest. Close to the holes we force the shift to be the shift of a single black hole in the hope that this choice will produce a hole that is at rest in coordinate space.

For the conformal lapse we use

$$\stackrel{~}{\alpha }=N_A+N_B1$$
(52)
or

$$\stackrel{~}{\alpha }=N_AN_B.$$
(53)
The first choice of $`\stackrel{~}{\alpha }`$ follows the philosophy of adding quantities of each individual hole. However, $`\stackrel{~}{\alpha }`$ of Eq. (52) becomes negative sufficiently close to the center of each hole and is therefore a bad choice if the excised spheres are small. The choice (53) does not change sign and has at large distances the same behavior (same $`1/r`$ term) as (52).

## IV Numerical Implementation

We implemented an elliptic solver that can solve all three decompositions we described above in complete generality. The solver uses domain decomposition and can handle nontrivial topologies. It is based on pseudospectral collocation, that is, it expresses the solution in each subdomain as an expansion in basis functions. This elliptic solver is described in detail in a separate paper Pfeiffer et al. (2002).

From the computational domain we excise two spheres containing the singularities of the Kerr-Schild metric close to the center of each hole. Around each of the excised spheres we place a spherical shell. These shells are patched together with $`5\times 3\times 3=45`$ rectangular blocks, with the two blocks at the location of the spheres removed. Around these 43 blocks, another spherical shell is placed that extends far out, typically to an outer radius of $`10^7M`$. In the rectangular blocks, we expand in Chebyshev polynomials, while in the spheres we use Chebyshev polynomials radially and spherical harmonics for the angular variables. This setup is depicted in Fig. 1.

The domain decomposition in Fig. 1 is fairly complicated. Even if the shells were made as large as possible, they do not cover the full computational domain when the excised spheres are close together. Thus additional subdomains are needed in any case. Choosing the 43 cubes as depicted allows for relatively small inner shells and for a relatively large inner radius of the outer shell. Thus each shell covers a region of the computational domain in which the angular variations of the solution are fairly low, allowing for comparatively few angular basis-functions.

The code can handle a general conformal metric. In principle, the user needs to specify only $`\stackrel{~}{\gamma }_{ij}`$. Then the code computes $`\stackrel{~}{\gamma }^{ij}`$, and —using numerical derivatives— the Christoffel symbols, Ricci tensor and Riemann scalar. For the special case of the Kerr-Schild metric of a single black hole and the superposition of two Kerr-Schild metrics, Eq. (46), we compute first derivatives analytically and use numerical derivatives only to compute the Riemann tensor.

The solver implements Eqs. (12) and (13) for the conformal TT decomposition, Eqs. (24) and (25) for the physical TT decomposition, and Eqs. (31) and (32) for the conformal thin sandwich decomposition.

After solving for $`(\psi ,V^i)`$ \[conformal TT and physical TT\], or $`(\psi ,\beta ^i)`$ \[thin sandwich\] we compute the physical metric $`\gamma _{ij}`$ and the physical extrinsic curvature $`K^{ij}`$ of the solution. Utilizing these physical quantities $`(\gamma _{ij},K^{ij})`$, we implement several analysis tools. We evaluate the constraints in the form of Eq. (3) and (4), we compute ADM-quantities and we search for apparent horizons. Note that these analysis tools are completely independent of the particular decomposition; they rely only on $`\gamma _{ij}`$ and $`K^{ij}`$.

Next we present tests ensuring that the various systems of equations are solved correctly. We also include tests of the analysis tools showing that we can indeed compute constraints, ADM-quantities and apparent horizons with good accuracy.

### IV.1 Testing the conformal TT and physical TT decompositions

We can test the solver by conformally distorting a known solution of the constraint equations. Given a solution to the constraint equations $`(\gamma _{0ij},K_0^{ij})`$ pick functions

$$\mathrm{\Psi }>0,W^i$$
(54)
and set

$`\stackrel{~}{\gamma }_{ij}`$ $`=\mathrm{\Psi }^4\gamma _{0ij},`$ (55)
$`K`$ $`=K_0,`$ (56)
and

$`\stackrel{~}{M}^{ij}`$ $`=\mathrm{\Psi }^{10}\left(K_0^{ij}{\displaystyle \frac{1}{3}}\gamma _0^{ij}K_0\right)\mathrm{\Psi }^4(𝕃_0W)^{ij}`$ (57)
for conformal TT or
$`\stackrel{~}{M}^{ij}`$ $`=\mathrm{\Psi }^{10}\left(K_0^{ij}{\displaystyle \frac{1}{3}}\gamma _0^{ij}K_0(𝕃_0W)^{ij}\right)`$ (58)
for physical TT. With these freely specifiable pieces and appropriate boundary conditions, a solution of the conformal TT equations (12), (13) or the physical TT equations (24), (25) will be

$`\psi `$ $`=\mathrm{\Psi }`$ (59)
$`V^i`$ $`=W^i.`$ (60)
From Eq. (8) we recover our initial metric $`\gamma _{0ij}`$, and from Eq. (16) \[conformal TT\] or Eq. (27) \[physical TT\] we recover the extrinsic curvature $`K_0^{ij}`$.

In our tests we used the particular choices

$$\mathrm{\Psi }=1+\frac{8(r2)}{36+x^2+0.9y^2+1.3(z1)^2}$$
(61)
$$W^i=\frac{50(r2)}{(6^4+r^4)}(y,x,1).$$
(62)
These functions are plotted in Fig. 2. $`\mathrm{\Psi }`$ varies between 0.8 and 1.5, $`W^i`$ varies between $`\pm 0.5`$, and both take their maximum values around distance $`7`$ from the center of the hole. We used for $`(\stackrel{~}{\gamma }_{0ij},K_0^{ij})`$ a single, boosted, spinning black hole in Kerr-Schild coordinates.

Figure 3 shows results of testing the conformal TT decomposition on a single spherical shell. The numerical solution $`(\psi ,V^i)`$ converges to the analytic solutions $`(\mathrm{\Psi },W^i)`$ exponentially with the number of basis functions as expected for a properly constructed spectral method. Moreover, the reconstructed metric and extrinsic curvature satisfy the constraints.

Now we test the solver for the physical TT decomposition, and demonstrate that we can correctly deal with multiple domains. In this example the computational domain is covered by an inner spherical shell extending for $`1.5Mr10M`$. This shell is surrounded by 26 rectangular blocks that overlap with the shell and extend out to $`x,y,z=\pm 25M`$. Finally another spherical shell covers the region $`20Mr10^6M`$. As can be seen in Fig. 4, the solution converges again exponentially.

For realistic cases we do not know the analytic solution and therefore need a measure of the error. Our major tool will be the change in results between different resolution. In particular we consider the $`L_2`$ norm of the point-wise differences of the solution at some resolution and at the next lower resolution. This diagnostic is labeled by circles in Fig. 4. Since the solution converges exponentially, these circles essentially give the error of the lower of the two resolutions.

In addition to testing the equations, this example tests domain decomposition and the integration routines for the ADM quantities. The ADM quantities are computed by the standard integrals at infinity in Cartesian coordinates,

$`E_{ADM}`$ $`={\displaystyle \frac{1}{16\pi }}{\displaystyle _{\mathrm{}}}\left(\gamma _{ij,j}\gamma _{jj,i}\right)d^2S_i,`$ (63)
$`J_{(\xi )}`$ $`={\displaystyle \frac{1}{8\pi }}{\displaystyle _{\mathrm{}}}\left(K^{ij}\gamma ^{ij}K\right)\xi _jd^2S_i.`$ (64)
For the $`x`$-component of the linear ADM-momentum, $`\xi =\widehat{e}_x`$ in Eq. (64). The choice $`\xi =x\widehat{e}_yy\widehat{e}_x`$ yields the $`z`$-component of the ADM-like angular momentum as defined by York York, Jr. (1979). Since the space is asymptotically flat there is no distinction between upper and lower indices in Eqs. (63) and (64). Note that Eq. (63) reduces to the familiar monopole term

$$\frac{1}{2\pi }_{\mathrm{}}_r\psi dA$$
(65)
only for quasi-isotropic coordinates. Our outer domain is large, but since it does not extend to infinity, we extrapolate $`r\mathrm{}`$.

For a Kerr black hole with mass $`M`$ and spin $`\stackrel{}{a}`$, that is boosted to velocity $`\stackrel{}{v}`$, the ADM-quantities will be

$`E_{ADM}`$ $`=\gamma M,`$ (66)
$`\stackrel{}{P}_{ADM}`$ $`=\gamma M\stackrel{}{v},`$ (67)
$`\stackrel{}{J}_{ADM}`$ $`=\left[\gamma \stackrel{}{a}(\gamma 1){\displaystyle \frac{(\stackrel{}{a}\stackrel{}{v})\stackrel{}{v}}{\stackrel{}{v}^2}}\right]M,`$ (68)
where $`\gamma =(1\stackrel{}{v}^2)^{1/2}`$ denotes the Lorentz factor. Eq. (68) reflects the fact that under a boost, the component of the angular momentum perpendicular to the boost-direction is multiplied by $`\gamma `$.

The example in Fig. 4 uses $`\stackrel{}{v}=(0.2,0.3,0.4)`$, and $`\stackrel{}{a}=(1/4,1/4,1/6)M`$. The evaluation of the angular momentum $`J_z`$ seems to be less accurate since our current procedure to extrapolate to infinity magnifies roundoff. We plan to improve this in a future version of the code. Until then we seem to be limited to an accuracy of $`10^6`$.

### IV.2 Testing conformal thin sandwich equations

The previous decompositions could be tested with a conformally distorted known solution. In order to test the conformal thin sandwich decomposition we need to find an analytic decomposition of the form (28). To do this, we start with a stationary solution of Einstein’s equations and boost it with uniform velocity $`v^i`$. Denote the metric, extrinsic curvature, lapse and shift of this boosted spacetime by $`\stackrel{~}{\gamma }_{0ij}`$, $`K_0^{ij}=A_0^{ij}+\frac{1}{3}\gamma _0^{ij}K_0`$, $`N_0`$ and $`N_0^i`$, respectively. Since we boosted the static solution, we will not find $`_t\gamma _{ij}=0`$ if we evolve it with the shift $`N_0^i`$. However, all time-dependence of this spacetime is due to the uniform motion, so in the comoving reference frame specified by the shift $`N_0^i+v^i`$, we will find $`_t\gamma _{ij}=0`$. In this case, Eq. (6) yields

$$A_0^{ij}=\frac{1}{2N_0}(𝕃(N_0+v))^{ij}.$$
(69)
If we choose $`\stackrel{~}{\alpha }=N_0`$ and $`\stackrel{~}{u}^{ij}=0`$, the thin sandwich equations (31) and (32) will thus be solved by $`\psi =1`$ and $`\beta ^i=N_0^i+v^i`$. Similar to the conformal TT and physical TT decomposition above, we can also conformally distort the metric $`\gamma _{0ij}`$. Furthermore, we can consider nonvanishing $`\stackrel{~}{u}^{ij}`$. We arrive at the following method to test the solver for the conformal thin sandwich equations:

Given a boosted version of a stationary solution with shift $`N_0^i`$, lapse $`N_0`$, 3-metric $`\gamma _{0ij}`$, trace of extrinsic curvature $`K_0`$, and boost-velocity $`v^i`$. Pick any functions

$$\mathrm{\Psi }>0$$
(70)
$$W^i$$
(71)
and set

$`\stackrel{~}{\gamma }_{ij}`$ $`=\mathrm{\Psi }^4\gamma _{0ij}`$ (72)
$`K`$ $`=K_0`$ (73)
$`\stackrel{~}{\alpha }`$ $`=\mathrm{\Psi }^6N_0`$ (74)
$`\stackrel{~}{u}^{ij}`$ $`=\mathrm{\Psi }^4(𝕃_0W)^{ij}`$ (75)
Then a solution to the thin sandwich equations (31-32) will be

$`\psi `$ $`=\mathrm{\Psi }`$ (76)
$`\beta ^i`$ $`=N_0^i+v^i+W^i`$ (77)
assuming boundary conditions respecting this solution.

Figure 5 shows results of this test for a single spherical shell and a Kerr black hole with $`\stackrel{}{v}=(0.2,0.3,0.1)`$, $`\stackrel{}{a}=(0.4,0.3,0.1)M`$. The solution converges to the expected analytical result exponentially. In addition, apparent horizon searches were performed. For the numerically found apparent horizons, the apparent horizon area $`A_{AH}`$ as well as the apparent horizon mass

$$M_{AH}=\sqrt{\frac{A_{AH}}{16\pi }}$$
(78)
were computed. The figure compares $`M_{AH}`$ to the expected value

$$M\left(\frac{1}{2}+\frac{1}{2}\sqrt{1\frac{\stackrel{}{a}^2}{M^2}}\right)^{1/2}.$$
(79)
As described in Baumgarte et al. (1996); Pfeiffer et al. (2000), the apparent horizon finder expands the apparent horizon surface in spherical harmonics up to a fixed order $`L`$. For fixed $`L`$, the error in the apparent horizon mass is dominated by discretization error of the elliptic solver at low resolution $`N`$. As $`N`$ is increased, the discretization error of the elliptic solver falls below the error due to finite $`L`$. Then the error in $`M_{AH}`$ becomes independent of $`N`$. Since the expansion in spherical harmonics is spectral, the achievable resolution increases exponentially with $`L`$. Note that for exponential convergence it is necessary to position the rays in the apparent horizon finder at the abscissas of Gauss-Legendre integration.

### IV.3 Convergence of binary black hole solutions

Figure 6 present the convergence of the solver in the binary black hole case. In this particular example, the conformal TT equations were solved for two black holes at rest with coordinate separation of $`10M`$. The computational domain is structured as in Fig. 1. The excised spheres have radius $`r_{exc}=2M`$, the inner spherical shells extend to radius $`4M`$. The rectangular blocks cover space up to $`x,y,z=\pm 25M`$, and the outer spherical shell extending from inner radius $`20M`$ to an outer radius of $`R=10^7M`$.

We do not use fall-off boundary conditions at the outer boundary; we simply set $`\psi =1`$ there. This limits the computations presented in this paper to an accuracy of order $`1/R10^7`$. Figure 6 shows that even for the next to highest resolution ($`N80`$) the solution will be limited by the outer boundary condition. All results presented in the following section are obtained at resolutions around $`N80`$. If the need arises to obtain solutions with higher accuracy, one can easily change to a fall-off boundary condition, or just move the outer boundary further out.

## V Results

The purpose of this paper is to compare the initial-data sets generated by different decompositions using simple choices for the freely specifiable pieces in each decomposition. We solve

$``$ ConfTT: Conformal TT equations (12) and (13) with freely specifiable pieces and boundary conditions given by Eqs. (46), (47), (48) and (49).

$``$ PhysTT: Physical TT equations (24) and (25) with same freely specifiable pieces and boundary conditions as ConfTT.

$``$ CTS: Conformal thin sandwich equations (31) and (32) with freely specifiable pieces and boundary conditions given by Eqs. (46), (47), (50) and (51). The lapse $`\stackrel{~}{\alpha }`$ is given by either Eq. (52), or by Eq. (53).

We will apply the terms “ConfTT”, “PhysTT” and “CTS” only to these particular choices of decomposition, freely specifiable pieces and boundary conditions. When referring to different freely specifiable pieces, or a decomposition in general, we will not use these shortcuts. If we need to distinguish between the two choices of $`\stackrel{~}{\alpha }`$ in CTS, we will use “CTS-add” for the additive lapse Eq. (52) and “CTS-mult”for the multiplicative lapse Eq. (53). Below in section V.3 we will also introduce as a forth term “mConfTT”.

### V.1 Binary black hole at rest

We first examine the simplest possible configuration: Two black holes at rest with equal mass, zero spin, and with some fixed proper separation between the apparent horizons of the holes. We solve

* ConfTT
* PhysTT
* CTS (with both choices of $`\stackrel{~}{\alpha }`$).

In the comparisons, we also include inversion symmetric conformally flat initial data obtained with the conformal-imaging formalism.

We excised spheres with radius $`r_{exc}=2M`$, which is close to the event horizon for an individual Eddington-Finkelstein black hole. This results in the boundary conditions being imposed close to, but within the apparent horizons of the black holes. The centers of the excised spheres have a coordinate separation of $`d=10M`$.

We now discuss the solutions. The conformal factor $`\psi `$ is very close to $`1`$ for each of the three decompositions. It deviates from $`1`$ by less than 0.02, indicating that a conformal metric based of a superposition of Kerr-Schild metrics does not deviate far from the constraint surface.

Figure 7 presents a plot of the conformal factor along the axis through the centers of the holes. One sees that $`\psi `$ is close to $`1`$; however, between the holes ConfTT and CTS force $`\psi `$ in opposite directions. For CTS, $`\psi >1`$ between the holes, for ConfTT, $`\psi <1`$! The contour plots in Fig. 8 also show this striking difference between the decompositions.

The result of PhysTT was found to be almost identical with ConfTT. This is reasonable, since these two decompositions differ only in that in one case the TT decomposition is with respect to the conformal metric, and in the other case the TT decomposition it is with respect to the physical metric. Since $`\psi 1`$, the conformal metric is almost identical to the physical metric, and only minor differences arise. In the following we will often use ConfTT/PhysTT when referring to both data sets.

We performed apparent horizon searches for these cases. For all data sets, the apparent horizon is outside the sphere with radius $`2M`$, that is outside the coordinate location for the apparent horizon in a single hole spacetime. For ConfTT/PhysTT, the radius of the apparent horizon surface is $`2.05M`$, for CTS it is $`2.15M`$. We computed the apparent horizon area $`A_{AH}`$, the apparent horizon mass

$$M_{AH}=\sqrt{\frac{A_{AH}}{16\pi }}$$
(80)
of either hole, and the combined mass of both holes,

$$m=2M_{AH}.$$
(81)
There is no rigorous definition of the mass of an individual black hole in a binary black hole spacetime, and Eq. (80) represents the true mass on an individual black hole only in the limit of wide separation of the black holes. A hard upper limit on the possible gravitational radiation emitted to infinity during the coalescence process of a binary will be

$$E_{MPRC}=E_{ADM}\sqrt{\frac{2A_{AH}}{16\pi }},$$
(82)
where $`2A_{AH}`$ is the combined apparent horizon area of both holes. Thus, $`E_{MPRC}`$ represents the maximum possible radiation content (MPRC) of the initial data. This, of course, makes the unlikely assumption that the binary radiates away all of its angular momentum.

We also compute the proper separation $`\mathrm{}`$ between the apparent horizon surfaces along the straight line connecting the centers of the excised spheres. In order to compare different data sets we consider the dimensionless quantities $`\mathrm{}/E_{ADM}`$, $`E_{ADM}/m`$ and $`E_{MPRC}/E_{ADM}`$. We will also use $`E_b/\mu `$ which will be defined shortly.

Table 1 summarizes these quantities for all three decompositions. It also includes results for inversion symmetric initial data, which for black holes at rest reduces to the Misner dataMisner (1963)<sup>1</sup><sup>1</sup>1Although the Misner solution can be obtained analytically, we found it more convenient to solve the Hamiltonian constraint numerically. The configuration in Table 1 corresponds to a separation $`\beta =12`$ in terms of Cook (1991).. The results in Table 1 are intended to represent nearly the same physical configuration.

From Table 1, one finds that the black holes have roughly the same dimensionless proper separation. However, the scaled ADM-energy $`E_{ADM}/m`$ differs by as much as 4.7% between the different data sets. $`E_{MPRC}/E_{ADM}`$, which does not depend on any notion of individual black hole masses at all, differs by 16% between the different data sets.

The inversion symmetric data has lowest $`E_{ADM}/m`$ and $`E_{MPRC}/E_{ADM}`$, CTS has somewhat larger values, and ConfTT/PhysTT lead to the biggest values for $`E_{ADM}/m`$ and $`E_{MPRC}/E_{ADM}`$. This indicates that, relative to the sizes of the black holes, ConfTT/PhysTT and CTS probably contain some excess energy.

A slightly different argument uses the binding energy which is defined as

$$\frac{E_b}{\mu }\frac{E_{ADM}2M_{AH}}{\mu },$$
(83)
where $`\mu =M_{AH}/2`$ is the reduced mass. Two Newtonian point masses at rest satisfy

$$\frac{E_b}{\mu }=\frac{1}{\mathrm{}/m}.$$
(84)
From Table 1 we see that for ConfTT/PhysTT, $`|E_b/\mu |>(l/m)^1`$, and for CTS, $`|E_b/\mu |(l/m)^1`$. Since gravity in general relativity is typically stronger than Newtonian gravity, we find again that CTS and ConfTT/PhysTT contain too much energy relative to the black hole masses, ConfTT/PhysTT having even more than CTS.

The proper separation between the apparent horizons $`\mathrm{}/m`$ is about 4% smaller for CTS than for ConfTT/PhysTT. By Eq. (84) this should lead to a relative difference in binding energy of the same order of magnitude. Since $`E_b/\mu `$ differs by almost a factor of two between the different decompositions, the differences in $`\mathrm{}/m`$ play only a minor role.

### V.2 Configurations with angular momentum

Now we consider configurations which are approximating two black holes in orbit around each other. The conformal metric is still a superposition of two Kerr Schild metrics. The black holes are located along the $`x`$-axis with a coordinate separation $`d`$. For ConfTT/PhysTT, we boost the individual holes to some velocity $`\pm v\widehat{e}_y`$ along the $`y`$-axis. For CTS we go to a co-rotating frame with an angular frequency $`\stackrel{}{\mathrm{\Omega }}=\mathrm{\Omega }\widehat{e}_z`$. Thus, for each decomposition we have a two parameter family of solutions, the parameters being $`(d,v)`$ for ConfTT and PhysTT, and $`(d,\mathrm{\Omega })`$ for CTS.

By symmetry, these configuration will have an ADM angular momentum parallel to the $`z`$-axis which we denote by $`J`$. In order to compare solutions among each other, and against the conformally flat inversion symmetric data sets, we adjust the parameters $`(d,v)`$ and $`(d,\mathrm{\Omega })`$, such that each initial-data set has angular momentum $`J/\mu m=2.976`$ and a proper separation between the apparent horizons of $`l/m=4.880`$. In Ref. Cook (1994), these values were determined to be the angular momentum and proper separation of a binary black hole at the innermost stable circular orbit.

Table 2 lists the parameters corresponding to this situation as well as results for each initial-data set<sup>3</sup><sup>3</sup>3 Because of the Lorentz contraction, the apparent horizons for ConfTT/PhysTT intersect the sphere with radius $`2`$. In order to have the full apparent horizon inside the computational domain, the radius of the excised spheres was reduced to 1.9 for these data sets.. As with the configuration with black holes at rest, we find again that ConfTT/PhysTT and CTS lead to different ADM-energies. Now, $`E_{ADM}/m`$ and $`E_{MPRC}/E_{ADM}`$ differ by $`0.02`$ and $`0.013`$, respectively, between CTS and ConfTT/PhysTT. However, in contrast to the cases where the black holes at rest, now CTS and the inversion symmetric data set have very similar values for $`E_{ADM}/m`$ and $`E_{MPRC}/E_{ADM}`$.

### V.3 Reconciling conformal TT and thin sandwich

We now investigate further the difference between ConfTT/PhysTT and CTS. Since the resulting initial-data sets for PhysTT and ConfTT are very similar, we restrict our discussion to ConfTT.

#### V.3.1 Motivation

The construction of binary black hole data for the ConfTT/PhysTT cases produces an extrinsic curvature that almost certainly contains a significant TT component. It would be interesting to know how significant this component is to the value of the various physical parameters we are comparing. Ideally, we would like to completely eliminate the TT component and see what effect this has on the resulting data sets. Unfortunately, this is a difficult, if not impossible, task.

The TT component of a symmetric tensor $`M^{ij}`$ is defined as

$$M_{TT}^{ij}M^{ij}(𝕃Y)^{ij},$$
(85)
where the vector $`Y^i`$ is obtained by solving an elliptic equation of the form

$$\mathrm{\Delta }_𝕃Y^i=_jM^{ij}.$$
(86)
The problem resides in the fact that the meaning of the TT component depends of the boundary conditions used in solving (86).

For the ConfTT/PhysTT cases we are actually solving for a vector $`V^i`$ that is a linear combination of two components, one that solves an equation of the form of (86) to obtain the TT component of $`\stackrel{~}{M}^{ij}`$ and one that solves the momentum constraint. But by imposing inner-boundary conditions of $`V^i=0`$, we don’t specify the boundary conditions on either part independently. Nor is it clear what these boundary conditions should be. Since we cannot explicitly construct the TT component of the extrinsic curvature, we cannot eliminate it. Although it is not ideal, there is an alternative we can consider that does provide some insight into the importance of the initial choice of $`\stackrel{~}{M}^{ij}`$.

#### V.3.2 Black holes at rest

Consider the following numerical experiment for two black holes at rest: Given $`\stackrel{~}{M}^{ij}`$ from Eq. (48), make a transverse traceless decomposition by setting

$$2N\stackrel{~}{M}^{ij}=\stackrel{~}{M}_{TT}^{ij}+(\stackrel{~}{𝕃}Y)^{ij}$$
(87)
where $`\stackrel{~}{}_j\stackrel{~}{M}_{TT}^{ij}=0`$ and $`N=N_A+N_B1`$. Notice that we are decomposing $`2N\stackrel{~}{M}^{ij}`$, not $`\stackrel{~}{M}^{ij}`$. Taking the divergence of Eq. (87) one finds

$$\stackrel{~}{\mathrm{\Delta }}_𝕃Y^i=\stackrel{~}{}_j\left(2N\stackrel{~}{M}^{ij}\right).$$
(88)
The decomposition chosen in Eq. (87) is motivated by the conformal thin sandwich decomposition. With this decomposition we can, in fact, use the shift vector $`N^i`$ to fix boundary conditions on $`Y^i`$, just as it was used to fix the boundary conditions in Eqs. (51b51d). For the black holes at rest in this case, we have $`\mathrm{\Omega }=0`$. After solving Eq. (88) for $`Y^i`$, we can construct a new conformal extrinsic curvature by

$$\stackrel{~}{M}{}_{}{}^{}{}_{}{}^{ij}=\frac{1}{2N}(\stackrel{~}{𝕃}Y)^{ij}$$
(89)
which is similar to what would result if we could eliminate $`\stackrel{~}{M}_{TT}^{ij}`$ from $`\stackrel{~}{M}^{ij}`$. Using $`\stackrel{~}{M}^{ij}`$ in place of $`\stackrel{~}{M}^{ij}`$, we can again solve the conformal TT equations. The result of this modified conformal TT decomposition “mConfTT” is striking: Figure 7 shows that mConfTT generates a conformal factor $`\psi `$ that is very similar to $`\psi `$ of CTS. mConfTT is also included in Table 1 where it can be seen that the quantities $`E_{ADM}/m`$ and $`E_{MPRC}/E_{ADM}`$ differ only slightly between mConfTT and CTS.

The fact that modification of the extrinsic curvature changes the ADM-energy by such a large amount underlines the importance of a careful choice for the extrinsic curvature $`\stackrel{~}{M}^{ij}`$ in ConfTT/PhysTT. The extremely good agreement between CTS and mConfTT is probably caused by our procedure to determine $`\stackrel{~}{M}_{}^{ij}{}_{}{}^{}`$. We force the extrinsic curvature of mConfTT into the form Eq. (89). This is precisely the form of the extrinsic curvature in CTS, Eq. (28), even using the same function $`N`$ and the same boundary conditions on the vectors $`Y^i`$ and $`\beta ^i`$.

#### V.3.3 Black holes with angular momentum

We now apply the modified conformal TT decomposition to the orbiting configurations of section V.2. In the corotating frame, the black holes are at rest, and we start with $`\stackrel{~}{\gamma }_{ij}`$ and $`\stackrel{~}{M}^{ij}`$ of two black holes at rest with coordinate separation $`d=11.860`$. We now solve Eq. (88) with

$$N=N_A+N_B1$$
(90)
and corotating boundary conditions on $`Y^i`$ \[cf. Eqs. (51b)–(51d)\]:

$`Y^i`$ $`=N_A^i`$ sphere inside hole A, (91a)
$`Y^i`$ $`=N_B^i`$ sphere inside hole B, (91b)
$`Y^i`$ $`=\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{r}`$ outer boundary. (91c)
$`N_{A/B}`$ and $`N_{A/B}^i`$ are lapse and shift of individual Kerr-Schild black holes at rest. $`\stackrel{~}{M}_{}^{ij}{}_{}{}^{}`$ is again constructed by Eq. (89) and used in solving the conformal TT equations.

Results from this procedure are included in Table 2. Again, mConfTT generates results very close to CTS. $`E_{ADM}/m`$ changes by 1.8% of the total mass between ConfTT and mConfTT, again highlighting the importance of the extrinsic curvature.

### V.4 Dependence on the size of the excised spheres

The framework presented in this paper requires the excision of the singularities at the centers of each hole<sup>4</sup><sup>4</sup>4Marronetti and MatznerMarronetti and Matzner (2000) effectively excised the centers, too, by using “blending functions”.. So far we have used $`r_{exc}=2M`$ or $`r_{exc}=1.9M`$ in order to impose boundary conditions close to the apparent horizons, but different choices can be made. Indeed, one might expect that the boundary conditions (49) and (51) become “better” farther inside the apparent horizon, where the metric and extrinsic curvature of that black hole dominate the superposed metric $`\stackrel{~}{\gamma }_{ij}`$ and superposed extrinsic curvature $`\stackrel{~}{M}^{ij}`$.

In order to test this assumption, we solve the constraint equations for two black holes at rest for different radii $`r_{exc}`$. We find that for all three decompositions, the data sets depend strongly on the radius of the excised spheres.

Figure 9 presents plots of the conformal factor $`\psi `$ and $`V^x`$ for ConfTT with different $`r_{exc}`$. There is no clear sign of convergence of $`\psi `$ as $`r_{exc}0`$. For $`r_{exc}=0.2M`$, the conformal factor $`\psi `$ even oscillates close to the excised sphere. Table 3 displays various quantities for the ConfTT decomposition for different $`r_{exc}`$. As $`r_{exc}`$ varies between $`2.0M`$ and $`0.1M`$, the ADM-energy varies between $`2.065`$ and $`2.106`$, whereas the apparent horizon area changes by nearly 4%. The apparent horizons move around somewhat as $`r_{exc}`$ changes. Figure 10 shows the location of the apparent horizons for different $`r_{exc}`$.

For CTS-add (with $`\stackrel{~}{\alpha }=N_A+N_B1`$), the initial-data sets seem to diverge as $`r_{exc}`$ is decreased. This has to be expected, since this choice for $`\stackrel{~}{\alpha }`$ changes sign if the excised spheres become sufficiently small. Changing to $`\stackrel{~}{\alpha }=N_AN_B`$ so that the lapse does not change sign reduces this divergent behavior. Von Neumann boundary conditions on $`\psi `$ at the excised spheres,

$$\frac{\psi }{r}=0,$$
(92)
lead to an increase in $`A_{AH}`$ especially for large excised spheres. This combination of lapse $`\stackrel{~}{\alpha }`$ and boundary conditions exhibits the smallest variations in $`E_{ADM}/m`$; cuts through $`\psi `$, $`\beta ^x`$ and through the apparent horizons are shown in Figs. 11 and 12. From the three examined combinations of lapse and boundary conditions, the one shown behaves best, but there is still no convincing sign of convergence.

Table 4 presents ADM-energies and apparent horizon areas and masses for CTS with different $`r_{exc}`$ and different choices of lapse and boundary condition. From the unscaled ADM-energy $`E_{ADM}`$ it is apparent that $`\stackrel{~}{\alpha }=N_A+N_B1`$ diverges most strongly. Note that between all choices of lapse, boundary conditions and $`r_{exc}`$, the unscaled quantities $`E_{ADM}`$, $`M_{AH}`$, and $`\mathrm{}`$ exhibit a much broader variation than the scaled quantities $`E_{ADM}/m`$ and $`\mathrm{}/E_{ADM}`$.

## VI Discussion

Our results clearly show that different decompositions lead to different initial-data sets, even when seemingly similar choices for the freely specifiable pieces are used. From Tables 1 and 2, one sees that $`E_{ADM}/m`$ changes by as much as 0.029 between ConfTT/PhysTT and CTS. The difference between ConfTT/PhysTT and the inversion symmetric data is even larger, 0.047. These numbers seem to be small; however, current evolutions of binary data usually find the total energy emitted in gravitational radiation $`E_{GW}/m`$ to be between 0.01 and 0.03 Baker et al. (2001); Alcubierre et al. (2001); Brandt et al. (2000), which is the same order of magnitude as the changes in $`E_{ADM}/m`$ we find. This means that, in principle, most of the energy radiated in these simulations could originate from “spurious” energy in the system and not from the dynamics of the binary system we are interested in.

These findings highlight the fact that current binary black hole initial data sets are inadequate for the task of accurately describing realistic binary systems. We see that the choices of the conformal 3-geometry $`\stackrel{~}{\gamma }_{ij}`$ and the freely specifiable portions of the extrinsic curvature, embedded in $`\stackrel{~}{M}^{ij}`$, influence the content of the initial data at a significant level. Furthermore, the results suggest that small changes in the free data associated with the extrinsic curvature are more significant than small changes in the choice of $`\stackrel{~}{\gamma }_{ij}`$.<sup>5</sup><sup>5</sup>5Following submission of this paper, a preprint by Damour et al.Damour et al. (2002) has appeared that lends support to our idea that the extrinsic curvature plays a key role in constructing quasi-equilibrium binary black hole initial data. This assertion is supported by the fact that $`E_{ADM}/m`$ is consistently larger for the ConfTT solutions than for the CTS solutions but the two approaches can be made to produce quite consistent results by using the modified extrinsic curvature of the mConfTT method. All of these decompositions use the same non-flat conformal metric, but differ in the extrinsic curvature. On the other hand, results for the conformally-flat inversion-symmetric data agree rather well with the results from the CTS method when we consider orbiting black holes. For black holes at rest, CTS differs from the inversion symmetric data, which seems to contradict our conclusion. However, this difference is likely due to the time-symmetric nature of the inversion symmetric data, which is especially adapted to the time-symmetry of the particular configuration of “two black holes at rest”.

Improved binary black hole initial data will require choices for the freely specifiable data that are physically motivated, rather than chosen for computational convenience. The same is true for the boundary conditions used in solving the constraints. The boundary conditions used in this paper carry the implicit assumptions that the approximate metric and extrinsic curvature are correct at the excision boundaries and that the value of the single-hole Kerr-Schild shift at the excision boundary is correct in a multi-hole situation. This is clearly not true, but we might hope that the impact of the error in this choice would diminish as we decrease the radius of the excision boundary. However, our results presented in Tables 3 and 4 do not support this conjecture. Examining the change in $`E_{ADM}/m`$ as we vary $`r_{exc}`$ shows only a small change, but more importantly, it shows no sign of converging as we decrease $`r_{exc}`$. The effects of changing $`r_{exc}`$ are much more significant for $`\mathrm{}/E_{ADM}`$, changing its value by as much as 10% in the case of CTS-mult for the range of values considered. Furthermore, as with the energy, we see no sign of convergence in $`\mathrm{}/E_{ADM}`$ as $`r_{exc}`$ decreases. Interestingly, although the solutions show no sign of convergence as we shrink the excision radius, we do find that the dimensionless quantities $`E_{ADM}/m`$ and $`\mathrm{}/E_{ADM}`$ do become independent of the choice of the inner-boundary condition as $`r_{exc}`$ decreases. This can be seen in comparing the result in Table 4 for the cases using $`\psi =1`$ and $`\psi /r=0`$ as inner-boundary conditions. Additional tests, not reported in this paper, further support this assertion.

## VII Conclusion

Using a new elliptic solver capable of solving the initial-value problem of general relativity for any of three different decompositions and any choice for the freely specifiable data, we have examined data sets representing binary black hole spacetimes. We find that the choices for the freely specifiable data currently in use are inadequate for the task of simulating the gravitational radiation produced in astrophysically realistic situations. In particular, we studied the results of using a superposition of two Kerr-Schild black holes to fix the freely specifiable data and compared them to the results obtained from conformally flat initial data.

Although the new Kerr-Schild based data provide a valuable point of comparison, it is not clear that the data produced are significantly superior to previous conformally-flat data. What is clear is that the choice of the freely specifiable data will be very important in constructing astrophysically realistic binary black hole initial data. Progress will require that these data, and the boundary conditions needed to solve the constraints, must be chosen based on physical grounds rather than computational convenience.

How can better initial data be achieved and how can the quality of initial data be measured? We believe that the conformal thin sandwich decomposition will be especially useful. Genuine radiative degrees of freedom cannot in principle be recognized on a single time slice. The conformal thin sandwich method uses in effect two nearby surfaces, giving it a potential advantage over other methods. Also, it avoids much of the uncertainty related to specifying a conformal extrinsic curvature. Moreover, the conformal thin sandwich approach is especially well suited for the most interesting configurations, a black-hole binary in a quasi-equilibrium orbit. In this case time derivatives of all quantities are small and the choice $`\stackrel{~}{u}^{ij}=0`$ is physically motivated. One should exploit the condition of quasi-equilibrium as fully as possible, i.e. one should use the conformal thin sandwich approach together with the constant $`K`$ equation, $`_tK=0`$. The latter yields another elliptic equation for the lapse which removes the arbitrariness inherent in choosing a conformal lapse $`\stackrel{~}{\alpha }`$. One will also need more physical boundary conditions. Work in this direction was begun inGourgoulhon et al. (2002); Grandclément et al. (2002) and refined in Cook (2002).

Ultimately, the gravitational wave content of an initial-data set can be determined only by long term evolutions. One must compute an initial-data set representing a binary black hole in quasi-circular orbit and evolve it. Then one must repeat this process with an initial-data set representing the same binary black hole, say, one orbital period earlier, and evolve that data set, too. If both evolutions lead to the same gravitational waves (modulo time offset) then one can be confident that the gravitational radiation is indeed astrophysically realistic. This approach has recently been used for the first time in conjunction with conformally flat puncture data, where it proved remarkably successful Baker et al. (2002).

###### Acknowledgements.

We thank Lawrence Kidder, Mark Scheel, and James York for helpful discussions. This work was supported in part by NSF grants PHY-9800737 and PHY-9900672 to Cornell University, and by NSF grant PHY-9988581 to Wake Forest University. Computations were performed on the IBM SP2 of the Department of Physics, Wake Forest University, with support from an IBM SUR grant.