Initial Data for Numerical Relativity

Abstract

Initial data are the starting point for any numerical simulation. In the case of numerical relativity, Einstein’s equations constrain our choices of these initial data. We will examine several of the formalisms used for specifying Cauchy initial data in the 3 + 1 decomposition of Einstein’s equations. We will then explore how these formalisms have been used in constructing initial data for spacetimes containing black holes and neutron stars. In the topics discussed, emphasis is placed on those issues that are important for obtaining astrophysically realistic initial data for compact binary coalescence.

1 Introduction

The goal of numerical relativity is to study spacetimes that cannot be studied by analytic means. The focus is therefore primarily on dynamical systems. Numerical relativity has been applied in many areas: cosmological models, critical phenomena, perturbed black holes and neutron stars, and the coalescence of black holes and neutron stars, for example. In any of these cases, Einstein’s equations can be formulated in several ways that allow us to evolve the dynamics. While Cauchy methods have received a majority of the attention, characteristic and Reggi calculus based methods have also been used. All of these methods begin with a snapshot of the gravitational fields on some hypersurface, the initial data, and evolve these data to neighboring hypersurfaces.

The focus of this review is on the initial data needed for Cauchy evolutions of Einstein’s equations. These initial data cannot be freely specified in their entirety. Rather they are subject to certain constraints which must be satisfied. Because of the nonlinearity of Einstein’s equations, there is no unique way of choosing which pieces of the initial data can be freely specified and which are constrained. In Section 2, I will look at the various formalisms that exist for expressing the initial-value equations. I will use the 3 + 1 (or ADM) decomposition of Einstein’s equations throughout and I begin the review with a brief introduction of this in Section 2.1. In Section 2.2, Section 2.3, and Section 2.4, I will explore the most important and widely used decompositions that have been developed.

In the remainder of the review, I will focus on initial data for black holes and neutron stars. Section 3 deals with black-hole initial data, which has received considerable attention over the years. The evolution of black-hole spacetimes is of particular importance because it allows for the study of pure geometro-dynamics. The spacetime is either pure vacuum or any matter can be hidden behind an event horizon. Thus, black-hole spacetimes allow us to study the dynamics of gravity alone. In Section 3.1, I begin with a review of some of the classic analytic solutions of the initial-data equations. Section 3.2 covers current methods for constructing general multi-black-hole initial data sets and also explores some of the limitations of these data sets. Section 3.3 deals with a class of black-hole initial data that has received little attention until recently and which I feel may play an important role in future work on constructing black-hole initial data. Finally, Section 3.4 examines the issue of constructing black-hole initial data for quasiequilibrium binaries.

This last topic is of extreme importance. Several laser interferometer gravitational wave detectors will become operational in the near future (cf. Ref. [89] for a Living Review on this topic) and the coalescence of black-hole binaries is considered to be one of the strongest candidates for detection by the earliest generation of detectors. The chances of detecting these events and then unraveling the information contained in the gravitational wave signals will be greatly increased if we have accurate numerical simulations of black-hole binary coalescence. While post-Newtonian techniques can be used to simulate the inspiral of a compact binary system, the final plunge and coalescence must be simulated numerically. Astrophysically realistic initial data will be needed before these simulations can provide reliable results.

The final plunge and coalescence of neutron-star binaries must also be simulated numerically. Section 4 examines the construction of neutron-star initial data. The major difference between black-hole and neutron-star initial data is the need to deal with the neutron star’s matter. In this section, I will deal with the neutron star matter in general, without considering any particular equation of state. In Section 4.1, I look at the issue of hydrostatic equilibrium of matter. Section 4.2 takes a brief look at constructing equilibrium initial data for isolated neutron stars. Finally, Section 4.3 examines the issues involved in constructing quasiequilibrium neutron-star binary initial data. In particular, we look at recent advances in constructing irrotational fluid models.

This article is far from a complete review of all the work that has been done on initial data for numerical relativity. I offer my apologies to all whose work I have not included. In particular, I have not covered any of the issues associated with existence and uniqueness of solutions to the constraint equations, nor dealt extensively with the issue of asymptotic falloff rates. I welcome correspondence on topics that you feel should be covered in this review, references that you think should be included, and the inevitable typographical errors.

1.1 Conventions

I will use a 4-metric g_μν with signature (-, +, +, +). Following the MTW [79] conventions, I define the Riemann tensor as

$${R^\mu }_{\nu \alpha \beta } \equiv {\partial _\alpha }{\Gamma ^\mu }_{\nu \beta } - {\partial _\beta }{\Gamma ^\mu }_{\nu \alpha } + {\Gamma ^\mu }_{\sigma \alpha }{\Gamma ^\sigma }_{\nu \beta } - {\Gamma ^\mu }_{\sigma \beta }{\Gamma ^\sigma }_{\nu \alpha }.$$

((1))

The Ricci tensor is defined as

$${R_{\mu \nu }} \equiv {R^\sigma }_{\mu \sigma \nu },$$

((2))

and the Einstein tensor as

$${G_{\mu \nu }} \equiv {R_{\mu \nu }} - {\textstyle{1 \over 2}}{g_{\mu \nu }}R.$$

((3))

A spatial 3-metric will be written as γ_ij and the Riemann and Ricci tensors associated with it will be defined by (1) and (2). Four-dimensional indices will be denoted by Greek letters, while 3-dimensional indices will be denoted by Latin letters.

To avoid confusion, the covariant derivative and Ricci tensor associated with γ_ij will be written with over-bars - ̄∇_j and ̄R_ij. I will also frequently deal with an auxiliary 3-dimensional space with a metric that is conformally related to the metric γ_ij of the physical space. The metric for this space will be denoted ̃γ_ij. The covariant derivative and Ricci tensor associated with this metric will be written with tildes - ̄∇_j and ̄R_ij. Other quantities that have a conformal relationship to quantities in the physical space will also be written with a tilde over them.

2 The Initial-Value Equations

Einstein’s equations, G_μν = 8πGTμν, represent ten independent equations. Since there are ten equations and ten independent components of the 4-metric g_μν, it seems that we have the same number of equations as unknowns. From the definition of the Einstein tensor (3), we see that these ten equations are linear in the second derivatives and quadratic in the first derivatives of the metric. We might expect that these ten second-order equations represent evolution equations for the ten components of the metric. However, a close inspection of the equations reveals that only six of the ten involve second time-derivatives of the metric. The remaining four equations are not evolution equations. Instead, they are constraint equations. The full system of equations is still well posed, however, because of the Bianchi identities

$${\nabla _\nu }{G^{\mu \nu }} \equiv 0.$$

((4))

The four constraint equations appear as a result of the general covariance of Einstein’s theory, which gives us the freedom to apply general coordinate transformations to each of the four coordinates and leave the interval

$${\rm{d}}{s^2} = {g_{\mu \nu }}{\rm{d}}{x^\mu }{\rm{d}}{x^\nu }$$

((5))

unchanged.

If we consider Einstein’s equations as a Cauchy problem¹, we find that the ten equations separate into a set of four constraint or initial-value equations, and six evolution or dynamical equations. If the four initial-value equations are satisfied on some spacelike hypersurface, which we can label with t = 0, then the Bianchi identities (4) guarantee that the evolution equations preserve the constraints on neighboring spacelike hypersurfaces.

2.1 Initial Data

In the Cauchy formulation of Einstein’s equations, we begin by foliating the 4-dimensional manifold as a set of spacelike, 3-dimensional hypersurfaces (or slices) {Σ}. These slices are labeled by a parameter t or, more simply, each slice of the 4-dimensional manifold is a t = constant hypersurface. Following the standard 3 + 1 decomposition [6, 111, 37], we let n^μ be the future-pointing timelike unit normal to the slice, with

$${n^\mu } \equiv - \alpha {\nabla ^\mu }t.$$

((6))

Here, α is called the lapse function (frequently denoted N in the literature). The scalar lapse function sets the proper interval measured by observers as they move between slices on a path that is normal to the hypersurface (so-called normal observers):

$${\text{d}}s{|_{along\,\,{n^\mu }}} = \alpha dt.$$

((7))

Of course, there is no reason that observers must move along a path normal to the hypersurface. In general, we can define the time vector as

$${t^\mu } \equiv \alpha {n^\mu } + {\beta ^\mu },$$

((8))

where

$${\beta ^\mu }{n_\mu } \equiv 0.$$

((9))

Here, α^μ is called the shift vector (frequently denoted N^μ in the literature).

Because of (9), α^μ has only three independent components and is a spatial vector, tangent to the hypersurface on which it resides. At this point, it is convenient to introduce a coordinate system adapted to the foliation {Σ}. Let xⁱ be the spatial coordinates in the slice. The fourth coordinate, t, is the parameter labeling each slice. With this adapted coordinate system, we find that 3-dimensional coordinate values remain constant as we move between slices along the t^μ direction (8). The four parameters, α and βⁱ, are a manifestation of the 4-dimensional coordinate invariance, or gauge freedom, in Einstein’s theory. If we let γ_ij represent the metric of the spacelike hypersurfaces, then we can rewrite the interval (5) as

$${\rm{d}}{s^2} = - {\alpha ^2}{\rm{d}}{t^2} + {\gamma _{ij}}({\rm{d}}{x^i} + {\beta ^i}{\rm{d}}t)({\rm{d}}{x^j} + {\beta ^j}{\rm{d}}t).$$

((10))

In the Cauchy formulation of Einstein’s equations, γ_ij is regarded as the fundamental variable and values for its components must be given as part of a well-posed initial-value problem. Since Einstein’s equations are second order, we must also specify something like a time derivative of the metric. For this, we use the second fundamental form, or extrinsic curvature, of the slice, K_ij, defined² by

$${K_{ij}} \equiv - \frac{1}{2}{{\mathcal L}_n}{\gamma _{ij}},$$

((11))

where ${{\mathcal L}_n}$ denotes the Lie derivative along the n^μ direction.

Together, γ_ij and K_ij are the minimal set of initial data that must be specified for a Cauchy evolution of Einstein’s equations. The metric γ_ij on a hypersurface is induced on that surface by the 4-metric g_μν. This means that the values γ_ij receives depend on how Σ is embedded in the full spacetime. In order for the foliation of slices {Σ} to fit into the higher-dimensional space, they must satisfy the Gauss-Codazzi-Ricci conditions. Combining these conditions with Einstein’s equations, and using (10), the six evolution equations become

$$\begin{array}{l} {\partial _t}{K_{ij}} = \alpha [{{\bar R}_{ij}} - 2{K_{i\ell }}K_j^\ell + K{K_{ij}} - 8\pi G{S_{ij}} + 4\pi G{\gamma _{ij}}(S - \rho )]\\ \;\;\;\;\;\;\;\;\;\;\;\; - {{\bar \nabla }_i}{{\bar \nabla }_j}\alpha + {\beta ^\ell }{{\bar \nabla }_\ell }{K_{ij}} + {K_{i\ell }}{{\bar \nabla }_j}{\beta ^\ell } + {K_{j\ell }}{{\bar \nabla }_i}{\beta ^\ell }. \end{array}$$

((12))

Here, ̄∇_i is the spatial covariant derivative compatible with γ_ij, ̄R_ij is the Ricci tensor associated with γ_ij, K = K ⁱ _i , ρ is the matter energy density, S_ij is the matter stress tensor, and S = S ⁱ _i .³ We have also used the fact that in our adapted coordinate system, ${{\mathcal L}_t} \equiv {\partial _t}$. The set of second-order evolution equations is completed by rewriting the definition of the extrinsic curvature (11) as

$${\partial _t}{\gamma _{ij}} = - 2\alpha {K_{ij}} + {\bar \nabla _i}{\beta _j} + {\bar \nabla _j}{\beta _i}.$$

((13))

Equations (12) and (13) are a first-order representation of a complete set of evolution equations for given initial data γ_ij and K_ij. However, the data cannot be freely specified in their entirety. The four constraint equations, following the same procedure outlined above for the evolution equations, become

$$\bar R + {K^2} - {K_{ij}}{K^{ij}} = 16\pi G\rho $$

((14))

and

$${{\bar \nabla }_j}\;({K^{ij}} - {\gamma ^{ij}}K) = 8\pi G{j^i}.$$

((15))

Here, ̄R = ̄R ⁱ _i and jⁱ is the matter momentum density. Equation (14) is referred to as the Hamiltonian or scalar constraint, while (15) are referred to as the momentum or vector constraints. Valid initial data for the evolution equations (12) and (13) must satisfy this set of constraints. And, as mentioned earlier, the Bianchi identities (4) guarantee that the evolution equations will preserve the constraints on future slices of the evolution.

As we will see below, the Hamiltonian constraint (14) most naturally constrains the 3-metric γ_ij, while the momentum constraints (15) naturally constrain the extrinsic curvature K_ij. Taking the constraints into consideration, it seems that the 3-metric has five degrees of freedom remaining, while the extrinsic curvature has three. But we know that the gravitational field in Einstein’s theory has two dynamical degrees of freedom, so we expect that both γ_ij and K_ij should each have only two free components. The answer to this problem is, once again, the coordinate invariance of Einstein’s theory. This seems strange at first, because we have already used the coordinate invariance of the theory to narrow our scope from the ten components of the 4-metric to the six components of the 3-metric. However, the lapse and shift do not completely specify the coordinate gauge. Rather, they specify how an initial choice of gauge will evolve with the foliation. The metric on a given hypersurface retains full 3-dimensional coordinate invariance, reducing the number of freely specifiable components to two [107, 108]. There also remains one degree of gauge freedom associated with the time coordinate which must be fixed. Each hypersurface represents a t = const. slice of the spacetime, so how the initial hypersurface is embedded in the full 4-dimensional manifold represents our temporal gauge choice. There is no unique way to specify this choice, but it is often convenient to let the trace of the extrinsic curvature K represent this temporal gauge choice [108]. Thus, we find that we are108 allowed to choose freely five components of the 3-metric and three components of the extrinsic curvature. However, only two of the components for each field represent dynamical degrees of freedom, the remainder are gauge degrees of freedom.

The four constraint equations, (14) and (15), represent conditions which the 3-metric and extrinsic curvature must satisfy. But, they do not specify which components (or combination of components) are constrained and which are freely specifiable. In the weak field limit where Einstein’s equations can be linearized, there are clear ways to determine which components are dynamic, which are constrained, and which are gauge. However, in the full nonlinear theory, there is no unique decomposition. In this case, one must choose a method for decomposing the constraint equations. The goal is to transform the equations into standard elliptic forms which can be solved given appropriate boundary conditions [80, 82, 111]. Each different decomposition yields a unique set of elliptic equations to be solved and a unique set of freely specifiable parameters which must be fixed somehow. Seemingly similar sets of assumptions applied to different decompositions can lead to physically different initial conditions.

2.2 York-Lichnerowicz Conformal Decompositions

For general initial-data configurations, the most widely used class of constraint decompositions are the York-Lichnerowicz conformal decompositions. At their heart are a conformal decomposition of the metric and certain components of the extrinsic curvature, together with a transverse-traceless decomposition of the extrinsic curvature.

First, the metric is decomposed into a conformal factor ψ multiplying an auxiliary 3-metric [68, 107, 108]:

$${\gamma _{ij}} \equiv {\psi ^4}{\tilde \gamma _{ij}}.$$

((16))

The auxiliary 3-metric ̃γ_ij is often called the conformal or background 3-metric, and it carries five degrees of freedom. Its natural definition is given by

$${\tilde \gamma _{ij}} = {\gamma ^{ - 1/3}}{\gamma _{ij}},$$

((17))

leaving ̃γ = 1, but we are free to choose any normalization for ̃γ. Using (16), we can rewrite the Hamiltonian constraint (14) as

$${\tilde \nabla ^2}\psi - \frac{1}{8}\psi \tilde R - \frac{1}{8}{\psi ^5}{K^2} + \frac{1}{8}{\psi ^5}{K_{ij}}{K^{ij}} = - 2\pi G{\psi ^5}\rho ,$$

((18))

where ̃∇² = ̃∇ⁱ̃∇_i is the scalar Laplace operator, and ̃∇_i and ̃R are the covariant derivative and Ricci scalar associated with ̃γ_ij. Equation (18) is a quasilinear elliptic equation for the conformal factor ψ, and we see that the Hamiltonian constraint naturally constrains the 3-metric.

The conformal decomposition of the Hamiltonian constraint was proposed by Lichnerowicz. But, the key to the full decomposition is the treatment of the extrinsic curvature introduced by York [109, 110]. This begins by splitting the extrinsic curvature into its trace and tracefree parts,

$${K_{ij}} \equiv {A_{ij}} + \frac{1}{3}{\gamma _{ij}}K.$$

((19))

The decomposition proceeds by using the fact that we can covariantly split any symmetric tracefree tensor as follows:

$${S^{ij}} \equiv {({\rm{\mathbb{L}}}X)^{ij}} + {T^{ij}}.$$

((20))

Here, T^ij is a symmetric, transverse-traceless tensor (i. e., ∇_jT^ij = 0 and T ⁱ _i = 0) and

$${({\rm{\mathbb{L}}}X)^{ij}} \equiv {\nabla ^i}{X^j} + {\nabla ^j}{X^i} - \frac{2}{3}{\gamma ^{ij}}{\nabla _\ell }{X^\ell }.$$

((21))

After separating out the transverse-traceless portion of ${{\mathcal S}^{ij}}$, what remains, $({\Bbb LX)^{ij}}$, is referred to as its “longitudinal” part. We now want to apply this transverse-traceless decomposition to the tracefree part of the extrinsic curvature A_ij. However, the conformal decomposition of the metric leaves us with at least two ways to proceed.

The goal of the decomposition is to produce a coupled set of elliptic equations to be solved with some prescribed boundary conditions. We have already reduced the Hamiltonian constraint to an elliptic equation being solved on a background space in terms of differential operators that are compatible with the conformal 3-metric. In the end, we want to reduce the momentum constraints to a set of elliptic equations based on differential operators that are compatible with the same conformal 3-metric. But, the longitudinal operator (21) can be defined with respect to any metric space. In particular, it is natural to consider decompositions with respect to both the physical and conformal 3-metrics.

2.2.1 Conformal Transverse-Traceless Decomposition

Let us first consider decomposing A^ij with respect to the conformal 3-metric [111, 116]. As we will see, when certain assumptions are made, this decomposition has the advantage of producing a simpler set of elliptic equations that must be solved. The first step is to define the conformal tracefree extrinsic curvature ̃A^ij by

$${A^{ij}} \equiv {\psi ^{ - 10}}{\tilde A^{ij}}\;\;\;\;\;{\rm{or}}\;\;\;\;\;\;{A_{ij}} \equiv {\psi ^{ - 2}}{\tilde A_{ij}}.$$

((22))

Next, the transverse-traceless decomposition is applied to the conformal extrinsic curvature,

$${\tilde A^{ij}} \equiv {(\tilde {\rm{\mathbb{L}}}X)^{ij}} + {\tilde Q^{ij}}.$$

((23))

Note that the longitudinal operator ${\tilde {\mathbb{L}}}$ and the symmetric, transverse-tracefree tensor ̃Q^ij are both defined with respect to covariant derivatives compatible with ̃γ_ij.

Applying equations (16), (19), (21), (22), and (23) to the momentum constraints (15), we find that they simplify to

$${\tilde \Delta _{\rm{\mathbb{L}}}}{X^i} = \frac{2}{3}{\psi ^6}{\tilde \nabla ^i}K + 8\pi G{\psi ^{10}}{j^i},$$

((24))

where

$$\tilde \Delta _\mathbb{L} X^i \equiv \tilde \nabla _j (\tilde {\mathbb{L}}X)^{ij} = \tilde \nabla ^2 X^i + \tfrac{1} {3}\tilde \nabla ^i (\tilde \nabla _j X^j ) + \tilde R_j^i X^j ,$$

((25))

and we have used the fact that

$${\bar \nabla _j}{{\mathcal S}^{ij}} = {\psi ^{ - 10}}{\tilde \nabla _j}({\psi ^{10}}{{\mathcal S}^{ij}})$$

((26))

for any symmetric tracefree tensor ${{\mathcal S}^{ij}},$.

In deriving equation (24), we have also used the fact that ̃Q^ij is transverse (i. e. ̃v∇_j̃Q^ij = 0). However, in general, we will not know if a given symmetric tracefree tensor, say ̃M^ij, is transverse. By using (20) we can obtain its transverse-traceless part ̃Q^ij via

$$\tilde Q^{ij} \equiv \tilde M^{ij} - (\tilde {\mathbb{L}}Y)^{ij} ,$$

((27))

and using the fact that if ̃Q^ij is transverse, we find

$${\tilde \nabla _j}{\tilde Q^{ij}} \equiv 0 = {\tilde \nabla _j}{\tilde M^{ij}} - {\tilde \Delta _{\rm{\mathbb{L}}}}{Y^i}.$$

((28))

Thus, Eqs. (27) and (28) give us a general way of constructing the required symmetric transverse-traceless tensor from a general symmetric traceless tensor.

Using the linearity of ${\tilde {\mathbb{L}}}$, we can rewrite (23) as

$${\tilde A^{ij}} = {(\tilde {\rm{\mathbb{L}}}V)^{ij}} + {\tilde M^{ij}},$$

((29))

where

$${V^i} \equiv {X^i} - {Y^i}.$$

((30))

Similarly, using the linearity of ${\tilde \Delta _\Bbb L}$, we can rewrite (24) as

$${\tilde \Delta _{\rm{\mathbb{L}}}}{V^i} = \frac{2}{3}{\psi ^6}{\tilde \nabla ^i}K - {\tilde \nabla _j}{\tilde M^{ij}} + 8\pi G{\psi ^{10}}{j^i}.$$

((31))

By solving directly for Vⁱ, we can combine the steps of decomposing ̃M^ij with that of solving the momentum constraints.

After applying (19) and (22) to the Hamiltonian constraint (18), we obtain the following full decomposition, which I will list together here for convenience:

$$\begin{array}{*{20}c} {\quad \quad \quad \quad \quad \quad \quad \gamma _{ij} = \psi ^4 \tilde \gamma _{ij} ,} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {\text{ }}K^{ij} = \psi ^{ - 10} \tilde A^{ij} + \frac{1} {3}\psi ^{ - 4} \tilde \gamma ^{ij} K,} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad {\text{ }}\tilde A^{ij} = ({\tilde{\mathbb{L}}}V)^{ij} + \tilde M^{ij} ,} \\ {\quad \quad \quad \quad {\text{ }}\;\;\;\;\;\;\;\;\;\tilde \Delta _\mathbb{L} V^i - \frac{2} {3}\psi ^6 \tilde \nabla ^i K = - \tilde \nabla _j \tilde M^{ij} + 8\pi G\psi ^{10} j^i ,} \\ {\tilde \nabla ^2 \psi - \frac{1} {8}\psi \tilde R - \frac{1} {{12}}\psi ^5 K^2 + \frac{1} {8}\psi ^{ - 7} \tilde A_{ij} \tilde A^{ij} = - 2\pi G\psi ^5 \rho .} \\ \end{array} $$

((32))

In the decomposition given by (32), we are free to specify a symmetric tensor ̃γ_ij as the conformal 3-metric, a symmetric tracefree tensor ̃M^ij, and a scalar function K. Then, with given matter energy and momentum densities, ρ and jⁱ, and appropriate boundary conditions, the coupled set of constraint equations for ψ and Vⁱ are solved. Finally, given the solutions, we can construct the physical initial data, γ_ij and K^ij.

The decomposition outlined above has the interesting property that if we choose K to be constant and if the momentum density vanishes⁴, then the momentum constraint equations fully decouple from the Hamiltonian constraint. As we will see later, this simplification has proven to be useful.

2.2.2 Physical Transverse-Traceless Decomposition

Alternatively, we can decompose A^ij with respect to the physical 3-metric [81, 82, 83]. We decompose the extrinsic curvature as

$$A^{ij} \equiv ({\bar {\mathbb{L}}W})^{ij} + Q^{ij} .$$

((33))

In this case, the longitudinal operator ${\bar {\mathbb{L}}}$ and the symmetric transverse-tracefree tensor Q^ij are both defined with respect to covariant derivatives compatible with γ_ij.

Applying equations (16), (19), (21), (33), and (26) to the momentum constraint (15), we find that it simplifies to

$${\tilde \Delta _{\rm{\mathbb{L}}}}{W^i} + 6{(\tilde {\rm{\mathbb{L}}}W)^{ij}}{\tilde \nabla _j}\ln \psi = \frac{2}{3}{\tilde \nabla ^i}K + 8\pi G{\psi ^4}{j^i},$$

((34))

where we have used the fact that

$$({\bar {\mathbb{L}}W})^{ij} = \psi ^{ - 4} ({\tilde {\mathbb{L}}W})^{ij} .$$

((35))

As in the previous section, we will obtain the symmetric transverse-traceless tensor Q^ij from a general symmetric tracefree tensor ̃M^ij by using (20). In this case, we take

$$Q^{ij} \equiv \psi ^{ - 10} \tilde M^{ij} - ({\bar {\mathbb{L}}Z})^{ij} ,$$

((36))

and use the fact that Q^ij is transverse, to obtain

$${\tilde \Delta _{\mathbb{L}}} Z^i + 6({\tilde{ \mathbb{L}}}Z)^{ij} {\tilde \nabla} _j \ln \psi = \psi ^{ - 6} {\tilde \nabla} _j \tilde M^{ij} .$$

((37))

Again, we can define

$${V^i} \equiv {W^i} - {Z^i},$$

((38))

and use the linearity of ${\tilde {\mathbb{L}}}$ and ${{\tilde \Delta }_\Bbb L}$ to combine the process of obtaining the transverse-traceless part of M^ij and solving the momentum constraints. We obtain the following full decomposition, which I will list together here for convenience⁵:

$$\begin{array}{l} \;{\gamma _{ij}} = {\psi ^4}{{\tilde \gamma }_{ij}},\\ {K^{ij}} = {\psi ^{ - 4}}({{\tilde A}^{ij}} + \frac{1}{3}{{\tilde \gamma }^{ij}}K), \end{array}$$

In the decomposition given by (39), we are again free to specify a symmetric tensor ̃γ_ij as the conformal 3-metric, a symmetric tracefree tensor ̃M^ij, and a scalar function K. Then, with given matter energy and momentum densities, ρ and jⁱ, and appropriate boundary conditions, the coupled set of constraint equations for ψ and Vⁱ are solved. Finally, given the solutions, we can construct the physical initial data, γ_ij and K^ij.

Notice that, while very similar to the decomposition from Section 2.2.1, the sets of equations are distinctly different. In general, if we make the same choices for the freely specifiable data in both decompositions (i. e., we choose ̃γ_ij, ̃M^ij, and K the same), we will produce two different sets of initial data. Both will be equally valid solutions of the constraint equations, but they will have distinct physical properties.

There is at least one exception to this. Assume we have a valid set of initial data γ_ij and K_ij, which satisfies the constraint equations (14) and (15). For any everywhere-positive function Ψ we define our freely specifiable data as follows:

$$\begin{array}{*{20}c} {\quad \quad {\text{ }}\tilde A^{ij} = ({\tilde {\mathbb{L}}}V)^{ij} + \psi ^{ - 6} {\tilde M^{ij}} ,} \\ {{\tilde \Delta} _{\mathbb{L}} V^i + 6({\tilde {\mathbb{L}}}V)^{ij} \tilde \nabla _j \ln \psi = \frac{2} {3}\tilde \nabla ^i K - \psi ^{ - 6} {\tilde \nabla} _j \tilde M^{ij} + 8\pi G\psi ^4 j^i ,} \\ {\tilde \nabla ^2 \psi - \frac{1} {8}\psi \tilde R - \frac{1} {{12}}\psi ^5 K^2 + \frac{1} {8}\psi ^5 \tilde A_{ij} \tilde A^{ij} = - 2\pi G\psi ^5 \rho .\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ \end{array} $$

((39))

Then the solution to both sets of equations, assuming we use correct boundary conditions, will be ψ = Ψ and Vⁱ = 0, which yields the original data as the solution for each decomposition.

2.3 Conformal Thin-Sandwich Decomposition

The two general initial-value decompositions outlined in Section 2.2.1 and Section 2.2.2 require identical freely specified data (̃γ_ij, ̃M^ij, and K), yet they usually produce different physical initial data. One shortcoming of these approaches is that they provide no direct insight into how to choose the freely specifiable data. All of the data are determined by schemes that involve only a single spacelike hypersurface. The resulting constraint equations are independent of the kinematical variables α and βⁱ that govern how the coordinates move through spacetime, and thus there is no connection to dynamics. York’s conformal thin-sandwich decomposition [115] takes a different approach by considering the evolution of the metric between two neighboring hypersurfaces (the thin sandwich). This decomposition is very similar to an approach originated by Wilson [104, 46], but is somewhat more general. Perhaps the most attractive feature of this decomposition is the insight it yields into the choice of the freely specifiable data.

The decomposition begins with the standard conformal decomposition of the 3-metric (16). However, we next make use of the evolution equation for the metric (13) in order to connect the 3-metrics on the two neighboring hypersurfaces. Label the two slices by t and t′, with t′ = t + δt, then = ′γ_ij = γ_ij + (∂_tγ_ij) δt. We would like to specify how the 3-metric evolves, but we do not have full freedom to do this. We know we can freely specify only the conformal 3-metric, and similarly, we are free to specify only the evolution of the conformal 3-metric. We make the following definitions:

$$\begin{array}{*{20}{l}} {\;\,{{\tilde \gamma }_{ij}} \equiv {\Psi ^{ - 4}}{\gamma _{ij}},}\\ {{{\tilde M}^{ij}} \equiv {\Psi ^{10}}({K^{ij}} - \frac{1}{3}{\gamma ^{ij}}K),}\\ {\;\;K \equiv K_i^i.} \end{array}$$

((40))

$${u_{ij}} \equiv {\gamma ^{1/3}}{\partial _t}({\gamma ^{ - 1/3}}{\gamma _{ij}}),$$

((41))

$${\tilde u_{ij}} \equiv {\partial _t}{\tilde \gamma _{ij}},$$

((42))

and

$${\tilde \gamma ^{ij}}{\tilde u_{ij}} \equiv 0.$$

((43))

The latter definition is made for convenience, so that we can treat ψ, ̃γ_ij, and ũ_ij as regular scalars and tensors instead of as scalar- and tensor-densities within this thin-sandwich formalism.

The conformal scaling of u_ij follows directly from (16), (41), (42), (43), and the identity that, for any small perturbation, δ ln γ = γ^ijδγ_ij. The result is

$${u_{ij}} = {\psi ^4}{\tilde u_{ij}},$$

((44))

which relies on the useful intermediate result that⁶

$${\partial _t}\tilde \gamma = 0.$$

((45))

Equation (41) represents the tracefree part of the evolution of the 3-metric, so (13) becomes

$$u^{ij} = - 2\alpha A^{ij} + ({\bar {\mathbb{L}}}\beta )^{ij} .$$

((46))

Using the conformal scalings (22), (35), and (44), we obtain

$$\tilde A^{ij} = \frac{{\psi ^6 }} {{2\alpha }}\left( {({\tilde {\mathbb{L}}}\beta )^{ij} - \tilde u^{ij} } \right).$$

((47))

York has pointed out that it is natural to use the following conformal rescal-ing of the lapse:

$$\alpha = {\psi ^6}\tilde \alpha .$$

((48))

This rescaling follows naturally from the “slicing function” that replaces the usual lapse $(\alpha = \sqrt \gamma \tilde \alpha )$ which has been critical in solving several problems [4]. It also results in the natural conformal scaling (22) postulated for the tracefree part of the extrinsic curvature. Substituting (48) into (47) yields what is taken as the definition of the tracefree part of the conformal extrinsic curvature,

$$\tilde A^{ij} \equiv \frac{1} {{2\tilde \alpha }}\left( {({\tilde {\mathbb{L}}}\beta )^{ij} - \tilde u^{ij} } \right).$$

((49))

Because the tracefree extrinsic curvature satisfies the normal conformal scaling, the Hamiltonian constraint will take on the same form as in (32). However, the momentum constraint will have a very different form. Combining equations (16), (19), (21), (22), and (49) with the momentum constraint (15), we find that it simplifies to

$${\tilde \Delta _{\rm{\mathbb{L}}}}{\beta ^i} - {(\tilde {\rm{\mathbb{L}}}\beta )^{ij}}{\tilde \nabla _j}\ln \tilde \alpha = \frac{4}{3}\tilde \alpha {\psi ^6}{\tilde \nabla ^i}K + \tilde \alpha {\tilde \nabla _j}(\frac{1}{{\tilde \alpha }}{\tilde u^{ij}}) + 16\pi G\tilde \alpha {\psi ^{10}}{j^i}.$$

((50))

Let us, for convenience, group together all the equations that constitute the conformal thin-sandwich decomposition:

$$\begin{array}{*{20}c} {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {\text{ }}\gamma _{ij} = \psi ^4 \tilde \gamma _{ij} ,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ {K^{ij} = \psi ^{ - 10} \tilde A^{ij} + \frac{1} {3}\psi ^{ - 4} \tilde \gamma ^{ij} K,} \\ {\tilde A^{ij} = \frac{1} {{2\tilde \alpha }}\left( {({\tilde{\mathbb{L}}}\beta )^{ij} - \tilde u^{ij} )} \right),} \\ {{\tilde \Delta _{\mathbb{L}}} \beta ^i - ({\tilde{\mathbb{L}}}\beta )^{ij} \tilde \nabla _j \ln \tilde \alpha - \frac{4} {3}\tilde \alpha \psi ^6 \tilde \nabla ^i K = \tilde \alpha \tilde \nabla _j \left( {\frac{1} {{\tilde \alpha }}\tilde u^{ij} } \right) + 16\pi G\tilde \alpha \psi ^{10} j^i ,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ {\tilde \nabla ^2 \psi - \frac{1} {8}\psi \tilde R - \frac{1} {{12}}\psi ^5 K^2 + \frac{1} {8}\psi ^{ - 7} \tilde A_{ij} \tilde A^{ij} = - 2\pi G\psi ^5 \rho .\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ \end{array} $$

((51))

In this decomposition (51), we are free to specify a symmetric tensor ̃γ_ij as the conformal 3-metric, a symmetric tracefree tensor ũ^ij, a scalar function K, and the scalar function ̃α. Solving this set of equations with appropriate boundary conditions yields initial data γ_ij and K_ij on a single hypersurface. However, we also know the following: If we chose to use the shift vector obtained from solving (50) and the lapse from (48) via our choice of ̃α and our solution to the Hamiltonian constraint, then the rate of change of the physical 3-metric is given by

$$\begin{array}{l} {\partial _t}{\gamma _{ij}} = {u_{ij}} + \frac{2}{3}{\gamma _{ij}}\left( {{{\bar \nabla }_k}{\beta ^k} - \alpha K} \right)\\ \;\quad \quad = {\psi ^4}\left[ {{{\tilde u}_{ij}} + \frac{2}{3}{{\tilde \gamma }_{ij}}\left( {{{\tilde \nabla }_k}{\beta ^k} + 6{\beta ^k}{{\tilde \nabla }_k}\ln \psi - {\psi ^6}\tilde \alpha K} \right)} \right]. \end{array}$$

((52))

This direct information about the consequences of our choices for the freely specifiable data is something not present in the previous decompositions. As we will see later, this framework has been used to construct initial data that are in quasiequilibrium.

2.4 Stationary Solutions

When there is sufficient symmetry present, it is possible to construct initial data that are in true equilibrium. These solutions possess at least two Killing vectors, one that is timelike at large distances and one that is spatial, representing an azimuthal symmetry. When these symmetries are present, solving for the initial data produces a global solution of Einstein’s equations and the solution is said to be stationary. The familiar Kerr-Neumann solution for rotating black holes is an example of a stationary solution in vacuum. Stationary configurations supported by matter are also possible, but the matter sources must also satisfy the Killing symmetries, in which case the matter is said to be in hydrostatic equilibrium [8].

The basic approach for finding stationary solutions begins by simplifying the metric to take into account the symmetries. Many different forms have been used for the metric (cf. Refs. [8, 9, 22, 34, 62, 21]). I will use a decomposition that makes comparison with the previous decompositions straightforward. First, define the interval as

$${\rm{d}}{s^2} = - {\psi ^{ - 4}}{\rm{d}}{t^2} + {\psi ^4}[{A^2}({\rm{d}}{r^2} + {r^2}{\rm{d}}{\theta ^2}) + {B^2}{r^2}{\sin ^2}\theta {({\rm{d}}\phi + {\beta ^\phi }{\rm{d}}t)^2}].$$

((53))

This form of the metric can describe any stationary spacetime. Notice that the lapse is related to the conformal factor by

$$\alpha = {\psi ^{ - 2}},$$

((54))

and that the shift vector has only one component

$${\beta ^i} = (0,\;0,\:{\beta ^\phi }).$$

((55))

I have used the usual conformal decomposition of the 3-metric (16) and have written the conformal 3-metric with two parameters as

$${\tilde \gamma _{ij}} = \left( {\begin{array}{*{20}{c}} {{A^2}}&0&0\\ 0&{{A^2}{r^2}}&0\\ 0&0&{{B^2}{r^2}{{\sin }^2}\theta } \end{array}} \right).$$

((56))

The four functions ψ, β^ϕ, A, and B are functions of r and θ only.

The equations necessary to solve for these four functions are derived from the constraint equations (14) and (15), and the evolution equations (12) and (13). For the evolution equations, we use the fact that ∂_tγ_ij = 0 and ∂_tK_ij = 0. The metric evolution equation (13) defines the extrinsic curvature in terms of derivatives of the shift

$${K_{ij}} \equiv \frac{1}{{2\alpha }}({\bar \nabla _i}{\beta _j} + {\bar \nabla _j}{\beta _i}).$$

((57))

With the given metric and shift, we find that K = 0 and the divergence of the shift also vanishes. This means we can write the tracefree part of the extrinsic curvature as

$$A^{ij} = \psi ^{ - 10} \tilde A^{ij} = \frac{1} {2}\psi ^{ - 2} ({\tilde{\mathbb{L}}}\beta )^{ij} .$$

((58))

We find that the Hamiltonian and momentum constraints take on the forms given by the conformal thin-sandwich decomposition (51) with ũ^ij = K ≡ 0 and ̃α = ψ^-8. Only one of the momentum constraint equations is non-trivial, and we find that the constraints yield elliptic equations for ψ and β^ϕ. What remains unspecified as yet are A and B (i. e., the conformal 3-metric).

The conformal 3-metric is determined by the evolution equations for the traceless part of the extrinsic curvature. Of these five equations, one can be written as an elliptic equation for B, and two yield complementary equations that can each be solved by quadrature for A. The remaining equations are redundant as a result of the Bianchi identities.

Of course, the clean separation of the equations I have suggested above is an illusion. All four equations must be solved simultaneously, and clever combinations of the four metric quantities can greatly simplify the task of solving the system of equations. This accounts for the numerous different systems used for solving for stationary solutions.

3 Black Hole Initial Data

In this section, we will look at Cauchy initial data that represent one or more black holes in an asymptotically flat spacetime. The majority of these will be either vacuum solutions or solutions of the Einstein-Maxwell equations. With no matter to support the gravitational field, we find that we must usually use a spacetime with a non-trivial topology, although a black hole can be supported by a compact gravitational wave [2]. It is certainly possible to construct black-hole solutions supported by matter [90, 91, 5], but it is often desirable to avoid the complications of matter sources.

This raises a point about solutions of Einstein’s equations which we have not yet mentioned. When constructing solutions of Einstein’s initial-value equations, we are free to specify the topology of the initial-data hypersurface. Einstein’s equations of general relativity place no constraints on the topology of the spacetime they describe or of spacelike hypersurfaces that foliate it. For astrophysical black holes (i. e., black holes in an asymptotically flat spacetime), the freedom in the choice of the topology has relatively minor consequences. The primary effects of different topology choices are hidden within the black hole’s event horizon.

In the sections below, we will explore many of the existing black-hole solutions and the schemes for generating them.

3.1 Classic Solutions

3.1.1 Schwarzschild

Of course, the simplest black-hole solution is the Schwarzschild solution. It represents a static spacetime containing a single black hole that connects two causally disconnected, asymptotically flat universes. There are actually many different coordinate representations of the Schwarzschild solutions. The simplest representations are time-symmetric (K_ij= 0), and so exist on a “maximally embedded” spacelike hypersurface (K = 0). These choices fix the foliation {Σ}. Spherical symmetry fixes two of the three spatial gauge choices. If we choose an “areal-radial coordinate”⁷, then the interval is written as

$${\rm{d}}{s^2} = - \left( {1 - \frac{{2M}}{r}} \right){\rm{d}}{t^2} + {\left( {1 - \frac{{2M}}{r}} \right)^{ - 1}}{\rm{d}}{r^2} + {r^2}({\rm{d}}{\theta ^{\rm{2}}} + {\sin ^2}\theta {\rm{d}}{\phi ^{\rm{2}}}).$$

((59))

If we choose an isotropic radial coordinate, then the interval is written as

$${\rm{d}}{s^2} = - {\left( {\frac{{1 - {\textstyle{M \over {2\tilde r}}}}}{{1 + {\textstyle{M \over {2\tilde r}}}}}} \right)^2}{\rm{d}}{t^2} + {\left( {1 - \frac{M}{{2\tilde r}}} \right)^4}({\rm{d}}{\tilde r^2} + {\tilde r^2}{\rm{d}}{\theta ^2} + {\tilde r^2}{\sin ^2}\theta {\rm{d}}{\phi ^2}).$$

((60))

In both (59) and (60), M represents the mass of the black hole as measured at spacelike infinity. Both of these solutions exist on the same foliation of t = const. slices. But, notice that the 3-geometry of the slice associated with (60) is conformally flat, while the 3-geometry associated with (59) is not.

The solution given in (60) is easily generated by any of the methods in Section 2.2 or Section 2.3. By choosing a time-symmetric initial-data hypersurface, we immediately get K_ij = 0, which eliminates the need to solve the momentum constraints. If we choose the conformal 3-geometry to be given by a flat metric (in spherical coordinates in this case), then the vacuum Hamiltonian constraint (18) becomes

$${\tilde \nabla ^2}\psi = 0,$$

((61))

where ̃∇² is the flat-space Laplace operator. For the solution ψ to yield an asymptotically flat physical 3-metric, we have the boundary condition that ψ(̃r → ∞) = 1. The simplest solution of this equation is

$$\psi = 1 + \frac{M}{{2\tilde r}},$$

((62))

where we have chosen the remaining integration constant to give a mass at infinity of M.

We now have full Cauchy initial data representing a single black hole. If we want to generate a full solution of Einstein’s equations, we must choose a lapse and a shift vector and integrate the evolution equations (12) and (13). In this case, a reasonable approach for specifying the lapse is to demand that the time derivative of K vanish. For the case of K = 0, this yields the so-called maximal slicing equation which, for the current situation, takes the form

$${\tilde \nabla ^2}(\alpha \psi ) = 0.$$

((63))

If we choose boundary conditions so that the lapse is frozen on the event horizon (α(̃r = M/2) = 0) and goes to one at infinity, we find that the solution is

$$\alpha = \frac{{1 - \frac{M}{{2\tilde r}}}}{{1 + \frac{M}{{2\tilde r}}}}.$$

((64))

If we now choose βⁱ = 0, we find that the left-hand sides of the evolution equations (12) and (13) vanish identically, and we have found the static solution of Einstein’s equations given in (60). We can, of course, recover the usual Schwarzschild coordinate solution (59) by using the purely spatial coordinate transformation $r = \tilde r{\left( {1 + \frac{M}{{2\tilde r}}} \right)^2}$.

It is interesting to examine the differences in these two representations of the Schwarzschild solution. The isotropic radial coordinate representation is well behaved everywhere except, it seems, at ̃r = 0. However, even here, the solution is well behaved. The 3-geometry is invariant under the coordinate transformation

$$\tilde r \to {\left( {\frac{M}{2}} \right)^2}\frac{1}{{r'}}.$$

((65))

The event horizon at $\tilde <Emphasis Type="Italic">r</Emphasis> \to {\left( {\frac{M}{2}} \right)^2}\frac{1}{{r'}}.$ is a fixed-point set of the isometry condition (65) which identifies points in two causally disconnected, asymptotically flat universes. We see that ̃r = 0 is simply an image of infinity in the other universe [31].

Given our choice for the lapse (64), which is frozen on the event horizon, we find that the solution can cover only the exterior of the black hole. To cover any of the interior with the lapse pinned to zero at the horizon would require we use a slice that is not spacelike everywhere. This is exactly what happens when the usual Schwarzschild areal-radial coordinate is used. At the event horizon, r = 2M, there is a coordinate singularity, and inside this radius the t = const. hypersurface is no longer spacelike. It is impossible to perform a Cauchy evolution interior to the event horizon using the areal-radial coordinate and the given time slicing.

We find that a Cauchy evolution, using the usual Schwarzschild time slicing that is frozen at the horizon, is capable of evolving only the region exterior to the black hole’s event horizon. Portions of the interior of the black hole can be covered by an evolution that begins with data on a standard Schwarzschild time slice, but the result is not a time-independent solution. As we will see later, there are other slicings of the Schwarzschild spacetime that cover the interior of the black hole and yield time-independent solutions.

3.1.2 Time-Symmetric Multi-Hole Solutions

As we saw in Section 3.1.1, the simplest approach for generating initial data is to assume time symmetry and let the conformal 3-geometry be flat. We could have generated the exterior solution for Schwarzschild with the areal-radial coordinate by a clever choice for the conformal 3-geometry. The Hamiltonian constraint is still linear, even if ̃γ_ij is not flat. But there is no obvious motivation for the correct ̃γ_ij that would yield the desired solution.

One approach for generating a time-symmetric multi-hole solution is straightforward. Brill-Lindquist initial data [31, 69] again assume a flat conformal 3-geometry, and the only non-trivial constraint equation is the Hamiltonian constraint, which again takes the form given in (61). But this time, we use the linearity of the Hamiltonian constraint and choose the solution to be a superposition of solutions with the form of (62). More precisely, we choose the solution

$$\psi = 1 + \sum\limits_{\sigma = 1}^N {\frac{{{\mu _\sigma }}}{{2|{\rm{x}} - {{\rm{C}}_\sigma }|}}} .$$

((66))

Here, |x - C_σ| is a coordinate distance from the point C_σ in the Euclidean conformal space, and the μ_σare constants related to the masses of the black holes. Assuming the points C_σ are sufficiently far apart, this solution of the initial-value equations represents N black holes momentarily at rest in “our” asymptotically flat universe. As was the case for a single black hole, each singular point in the solution, x = C_σ, represents infinity in a different, causally disconnected universe. In fact, each black hole connects “our” universe to a different universe, so that there are N +1 asymptotically flat hypersurfaces connected together at the throats of N black holes. While we started with a 3-dimensional Euclidean manifold, the requirement that we delete the singular points, C_σ, results in a manifold that is not simply connected. This solution is often referred to as having a topology with N +1 “sheets”.

Brill-Lindquist initial data are very similar to the Schwarzschild initial data in isotropic coordinates except for one major difference: The solution does not represent two identical universes that have been joined together. The coordinate transformation (65) can still be used to show that each pole in the solution corresponds to infinity on an asymptotically flat hypersurface. But, the solution has N +1 different asymptotically flat universes connected together, not two, and “our” universe containing N black holes cannot be isometric to any of the other universes which contain just one. Interestingly, Misner [78] found that it is possible to construct a solution of the vacuum, time-symmetric Hamiltonian constraint (61) that has two isometric asymptotically flat hypersurfaces connected by N black holes. The case of two black holes that satisfies this isometry condition (often called inversion symmetry) is usually referred to as “Misner initial data”. It has an analytic representation in terms of an infinite series expansion. The construction is tedious, and there are several representations of the solution [93, 38]⁸. Like the Brill-Lindquist data, the non-simply connected topology of the full manifold is represented on a Euclidean 3-manifold by the presence of singular points that must be removed. Brill-Lindquist data have N singular points, each representing an image of infinity as seen through the throat of the black hole connecting our universe to that hole’s other universe. Misner’s data contain an infinite number of singular points for each black hole, each representing an image of one of two asymptotic infinities. The result is seen to represent two identical asymptotically flat universes joined by N black holes. The two universes join together at the throats of the N black holes, with each throat being a coordinate 2-sphere in the conformal space. Data at any point in one universe are related to data at the corresponding point in the alternate universe via the same isometry condition (65) found in the Schwarzschild case. And, as in the Schwarzschild case, the 2-sphere throat of each black hole forms a fixed-point set of the isometry condition.

Both the Brill-Lindquist data and Misner’s data represent N black holes at a moment of time-symmetry (i. e., all the black holes are momentarily at rest). Both are conformally flat and the difference in the topology of the two solutions is hidden from an observer outside the black holes. Yet solutions where the holes are chosen to have the same size and separation yield similar, but physically distinct, solutions [36, 3].

3.2 General Multi-Hole Solutions

Time-symmetric black-hole solutions of the constraint equations such as those described in Section 3.1 are useful as test cases because they have analytic representations. However, they have very little physical relevance. General time-asymmetric solutions are needed to represent black holes that are moving and spinning⁹. A few approaches for generating general multi-hole solutions have been explored and, below, we will look at the two approaches which are direct generalizations of the Misner and Brill-Lindquist data. Generalizing these two approaches was the most natural first step toward constructing general multi-hole solutions.

The first approach to be developed generalized the Misner approach [78]. It was attractive because an isometry condition relating the two asymptotically flat universes provides two useful things. First, because the two universes are identical, finding a solution in one universe means that you have the full solution. Second, because the throats are fixed-point sets of the isometry, we can construct boundary conditions on any quantity there. This allows us to excise the region interior to the spherical throats from one of the Euclidean background spaces and solve for the initial data in the remaining volume.

The generalization of the Misner approach seemed preferable to trying to solve the constraints on N +1 Euclidean manifolds stitched together smoothly at the throats of N black holes. However, Brandt and Brügmann [27] realized it was possible to factor out analytically the behavior of the singular points in the Euclidean manifold of the N + 1 sheeted approach. Referred to as the “puncture” method, this approach allows us to rewrite the constraint equations for functions on an N + 1 sheeted manifold as constraint equations for new functions on a simple Euclidean manifold.

Another approach tried, which we will not discuss in detail, avoided the issue of the topology of the initial-data slice entirely. Developed by Thornburg [100], this approach was based on the idea that only the domain exterior to the apparent horizon of a black hole is relevant. The equation describing the location of an apparent horizon can be rewritten in a form that can be used as a boundary condition for the conformal factor in the Hamiltonian constraint equation. Thus, given a compatible solution to the momentum constraints, this boundary condition can be used to construct a solution of the Hamiltonian constraint in the domain exterior to the apparent horizons of any black holes, with no reference at all to the topology of the full manifold.

3.2.1 Bowen-York Data

The generalization of the Misner approach was developed by York and his collaborators [25, 23, 24, 66, 113, 65, 26, 38]. The approach is often called the “conformal-imaging” method, and the data are usually referred to as “Bowen-York” data. This approach begins with a set of simplifying assumptions that is common to all three of the approaches described above. These assumptions are that

$$\begin{array}{*{20}{l}} {\,\:K = 0,\,\,\,\,\,\,\quad maximal\,slicing,}\\ {\:{{\tilde \gamma }_{ij}} = {f_{ij}},\,\,\,\,\,\;\;conformal\,flatness,}\\ {\psi {|_\infty } = 1,\,\,\,\,\,\,\quad asymptotic\,flatness.} \end{array}$$

((67))

Here, f_ij represents a flat metric in any suitable coordinate system. The assumption of conformal flatness means that the differential operators in the constraints are the familiar flat-space operators. More importantly, if we use the conformal transverse-traceless decomposition (32), we find that in vacuum the momentum constraints completely decouple from the Hamiltonian constraint.

The importance of this last property stems from the fact that York and Bowen were able to find analytic solutions of this version of the momentum constraints, solutions that represent a black hole with both linear momentum and spin [25, 23, 24]. If we choose ̃M^ij = 0, then the momentum constraints (31) become

$${\tilde \nabla ^2}{V^i} + \frac{1}{3}{\tilde \nabla ^i}({\tilde \nabla _j}{V^j}) = 0.$$

((68))

A solution of this equation is

$$V^{i}=-\displaystyle \frac{1}{4r}[7P^{i}+n^{i}n_{j}P^{j}]+\frac{1}{r^{2}}\epsilon^{ijk}n_{j}S_{k}.$$

((69))

Here, Pⁱ and Sⁱ are vector parameters, r is a coordinate radius, and nⁱ is the outward-pointing unit normal of a sphere in the flat conformal space (nⁱ = xⁱ/r). ∊^ijk is the 3-dimensional Levi-Civita tensor.

This solution of the momentum constraints yields the tracefree part of the extrinsic curvature,

$$\begin{array}{l} {{\tilde A}_{ij}} = \frac{3}{{2{r^2}}}[{P_i}{n_j} + {P_j}{n_i} - ({f_i}_j - {n_i}{n_j}){P^k}{n_k}]\\ \;\;\;\;\;\;\;\; + \frac{3}{{{r^3}}}[{ \epsilon _{ki\ell }}{S^\ell }{n^k}{n_j} + { \epsilon _{kj\ell }}{S^\ell }{n^k}{n_i}]. \end{array}$$

((70))

Remarkably, using this solution (70) and the assumptions in (67), we can determine the physical values for the linear and angular momentum of any initial data we can construct. The momentum contained in an asymptotically flat initial-data hypersurface can be calculated from the integral [112]

$${\Pi ^i}\xi _{(k)}^i = \frac{1}{{8\pi }}\oint_\infty {\left( {K_i^j - \delta _i^jK} \right)} \;\xi _{(k)}^i{{\rm{d}}^2}{S_j},$$

((71))

where ξ ⁱ_{( k)} is a Killing vector of the 3-metric γ_ij¹⁰. Since we are not likely to have true Killing vectors, we make use of the asymptotic translational and rotational Killing vectors of the flat conformal space. We find from (71), (70), and (67) that the physical linear momentum of the initial-data hypersurface is Pⁱ and the physical angular momentum of the slice is Sⁱ. Furthermore, because the momentum constraints are linear, we can add any number of solutions of the form of Eq. (70) to represent a collection of linear and angular momentum sources. The total physical linear momentum of the initial-data slice will simply be the vector sum of the individual linear momenta. The total physical angular momentum cannot be obtained by simply summing the individual spins because this neglects the orbital angular momentum of the various sources. However, the total angular momentum can still be computed without having to solve the Hamiltonian constraint [114].

The Bowen-York solution for the extrinsic curvature is the starting point for all the general multi-hole initial-data sets we have discussed in Section 3.2. However, this solution is not inversion symmetric. That is, it does not satisfy the isometry condition that any field must satisfy to exist on a two-sheeted manifold like that of Misner’s solution. Fortunately, there is a method of images, similar to that used in electrostatics but applicable to tensors, that can be used to make any tensor inversion symmetric [78, 25, 66, 65, 113].

For the conformal extrinsic curvature of a single black hole, there are two inversion-symmetric solutions [25],

$$\begin{array}{l} \tilde A_{ij}^ \pm = \frac{3}{{2{r^2}}}\left[ {{P_i}{n_j} + {P_j}{n_i} - ({f_i}_j - {n_i}{n_j}){P^k}{n_k}} \right]\\ \;\;\;\;\;\;\;\;\; \mp \frac{{3{a^2}}}{{2{r^4}}}\left[ {{P_i}{n_j} + {P_j}{n_i} + ({f_i}_j - 5{n_i}{n_j}){P^k}{n_k}} \right]\\ \;\;\;\;\;\;\;\;\; + \frac{3}{{{r^3}}}\left[ {{ \epsilon _{ki\ell }}{S^\ell }{n^k}{n_j} + { \epsilon _{kj\ell }}{S^\ell }{n^k}{n_i}} \right]. \end{array}$$

((72))

Here, a is the radius of the coordinate 2-sphere that is the throat of the black hole. Of course, this coordinate 2-sphere is the fixed-point set of the isometry and is the surface on which we can impose boundary conditions. Notice that this radius enters the solutions only when we make it inversion symmetric.

When the extrinsic curvature represents more than one black hole, the process for making the solution inversion symmetric is rather complex and results in an infinite-series solution. However, in most cases of interest, the solution converges rapidly and it is straightforward to evaluate the solution numerically [38].

Given an inversion-symmetric conformal extrinsic curvature, it is possible to find an inversion-symmetric solution of the Hamiltonian constraint [25]. Given our assumptions (67), the Hamiltonian constraint becomes

$${\tilde \nabla ^2}\psi + \frac{1}{8}{\psi ^{ - 7}}{\tilde A_{ij}}{\tilde A^{ij}} = 0.$$

((73))

The isometry condition imposes a condition on the conformal factor at the throat of each hole. This condition takes the form [25]

$$n_\sigma ^i{\tilde \nabla _i}{\left. \psi \right|_{{a_\sigma }}} = - {\left. {\frac{\psi }{{2{r_\sigma }}}} \right|_{{a_\sigma }}},$$

((74))

where n ⁱ _σ is the outward-pointing unit-normal vector to the σ^th throat and a_σ is the coordinate radius of that throat. This condition can be used as a boundary condition when solving (73) in the region exterior to the throats.

In addition to boundary conditions on the throats, a boundary condition on the outer boundary of the domain is needed before the quasilinear elliptic equation in (73) can be solved as a well-posed boundary-value problem. This final boundary condition comes from the fact that we want an asymptotically flat solution. This implies that the solution behaves as

$$\psi = 1 + \frac{E}{{2r}} + {\mathcal O}({r^{ - 2}}),$$

((75))

where E is the total ADM energy content of the initial-data hypersurface. Equation (75) can be used to construct appropriate boundary conditions either at infinity or at a large, but finite, distance from the black holes [116].

3.2.2 Puncture Data

The generalization of the Brill-Lindquist data developed by Brandt and Brügmann [27] begins with the same set of assumptions (67) as the conformal-imaging approach outlined in Section 3.2.1. We immediately have Eq. (70) from the solution of the momentum constraints, and we must solve the Hamiltonian constraint, which again takes the form of Eq. (73). At this point, however, the method of solution differs from the conformal-imaging approach.

Based on the time-symmetric solution, it is reasonable to assume that the conformal factor will take the form

$$\psi = \frac{1}{\chi } + u,\:\;\;\;\;\;\;\;\frac{1}{\chi } \equiv \sum\limits_{\sigma = 1}^N {\frac{{{\mu _\sigma }}}{{2|{\rm{x}} - {{\rm{C}}_\sigma }|}}} .$$

((76))

If u is sufficiently smooth, (76) implies that the manifold will have the topology of N +1 asymptotically flat regions just as in the Brill-Lindquist solution. In this case, asymptotic flatness requires that $u = 1 + {\mathcal O}({r^{ - 1}})$.

Substituting (76) into the Hamiltonian constraint (73) yields

$${\tilde \nabla ^2}u + \eta {(1 + \chi u)^{ - 7}} = 0,$$

((77))

where

$$\eta = \frac{1}{8}{\chi ^7}{\tilde A_{ij}}{\tilde A^{ij}}.$$

((78))

Near each singular point, or “puncture”, we find that χ ≈ |x - C_σ|. From (70), we see that AjAjj behaves no worse than Ã_ijÃ^ij behaves no worse than |x - C_σ|^-6, so η vanishes at the punctures at least as fast as |x - C_σ|.

With this behavior, Brandt and Brügmann [27] have shown the existence and uniqueness of C² solutions of the modified Hamiltonian constraint (77). The resulting scheme for constructing multiple black hole initial data is very simple. The mass and position of each black hole are parameterized by μ_σ and C_σ, respectively. Their linear momenta and spin are parameterized by P_σ and S_σ in the conformal extrinsic curvature (70) used for each hole. Finally, the solution for u is found on a simple Euclidean manifold, with no need for any inner boundaries to avoid singularities. This is a great simplification over the conformal-imaging approach, where proper handling of the inner boundary is the most difficult aspect of solving the Hamiltonian constraint numerically [41].

3.2.3 Problems with These Data

Both the conformal-imaging and puncture methods for generating multiple black hole initial data allow for completely general configurations of the relative sizes of the black holes, as well as their linear and angular momenta. This does not mean that these schemes allow for the generation of all desired black-hole initial data. The two schemes rely on specific assumptions about the freely specifiable gravitational data. In particular, they assume K = 0, ̃M^ij = 0, and, most importantly, that the 3-geometry is conformally flat.

These choices for the freely specifiable data are not always commensurate with the desired physical solution. For example, if we choose to use either method to construct a single spinning black hole, we will not obtain the Kerr solution. The Kerr-Newman solution can be written in terms of a quasi-isotropic radial coordinate on a K = 0 time slice [28]. Let r denote the usual Boyer-Lindquist radial coordinate and make the standard definitions

$${\rho ^2} \equiv {r^2} + {a^2}{\cos ^2}\theta \;\;\;\;{\rm{and}}\;\;\;\;\Delta \equiv {r^2} - 2Mr + {a^2} + {Q^2}.$$

((79))

A quasi-isotropic radial coordinate ̃r can be defined via

$$r = \tilde r\left( {1 + \frac{{M + \sqrt {{a^2} + {Q^2}} }}{{2\tilde r}}} \right)\left( {1 + \frac{{M - \sqrt {{a^2} + {Q^2}} }}{{2\tilde r}}} \right).$$

((80))

The interval then becomes

$${\rm{d}}{s^2} = - {\alpha ^2}{\rm{d}}{t^2} + {\psi ^4}[{e^{2\mu /3}}({\rm{d}}{\tilde r^2} + {\tilde r^2}{\rm{d}}{\theta ^2}) + {\tilde r^2}{\sin ^2}\theta {e^{ - 4\mu /3}}{({\rm{d}}\phi + {\beta ^\phi }{\rm{d}}t)^2}],$$

((81))

with

$$\begin{array}{l} {\alpha ^2} = \frac{{{\rho ^2}\Delta }}{{{{({r^2} + {a^2})}^2} - \Delta {a^2}{{\sin }^2}\theta }},\\ {\beta ^\phi } = - a\frac{{{r^2} + {a^2} - \Delta }}{{{{({r^2} + {a^2})}^2} - \Delta {a^2}{{\sin }^2}\theta }},\\ {\psi ^4} = \frac{{{\rho ^{2/3}}{{({{({r^2} + {a^2})}^2} - \Delta {a^2}{{\sin }^2}\theta )}^{1/3}}}}{{{{\tilde r}^2}}},\\ {e^{2\mu }} = \frac{{{\rho ^4}}}{{{{({r^2} + {a^2})}^2} - \Delta {a^2}{{\sin }^2}\theta }}. \end{array}$$

((82))

We see immediately that the 3-geometry associated with a t = const. hypersurface of (81) is not conformally flat. In fact, Garat and Price [54] have shown that in general there is no spatial slicing of the Kerr spacetime that is axisymmetric, conformally flat, and smoothly goes to the Schwarzschild solution as the spin parameter a → 0.

Since the Kerr solution is stationary, the inescapable conclusion is that conformally flat initial data for a single rotating black hole must also contain some nonvanishing dynamical component. When we evolve such data, the system will emit gravitational radiation and eventually settle down to the Kerr geometry [25, 29]. But, it cannot be the Kerr geometry initially, and it is unlikely that the spurious gravitational radiation content of the initial data has any desirable physical properties. Conformally flat initial data for spinning holes contain some amount of unphysical “junk” radiation. A similar conclusion is reached for conformally flat data for a single black hole with linear momentum [112].

The choice of a conformally flat 3-geometry was originally made for convenience. Combined with the choice of maximal slicing, these simplifying assumptions allowed for an analytic solution of the momentum constraints which vastly simplified the process of constructing black-hole initial data. Yet there has been much concern about the possible adverse physical effects that these choices (especially the choice of conformal flatness) will have in trying to study black-hole spacetimes [40, 70, 87, 60, 86]. While these conformally flat data sets may still be useful for tests of black-hole evolution codes, it is becoming widely accepted that the unphysical initial radiation will significantly contaminate any gravitational waveforms extracted from evolutions of these data. In short, these data are not astrophysically realistic.

The various initial-data decompositions outlined in Section 2.2 and Section 2.3 are all capable of producing completely general black-hole initial data sets. The only limitation of these schemes is our understanding of what choices to make for the freely specifiable data and the boundary conditions to apply when solving the sets of elliptic equations. All these choices will have a critical impact on the astrophysical significance of the data produced. It is also important to remember that similar choices for the freely specifiable data will result in physically different solutions when applied to the different schemes.

The first studies of black-hole initial data that are not conformally flat were carried out by Abrahams et al. [1]. They looked at the superposition of a gravitational wave and a black hole. Using a form of the conformal metric that allows for so-called Brill waves, they constructed time-symmetric initial data that were not conformally flat and yet satisfied the isometry condition (65) used in the conformal-imaging method of Section 3.2.1. These data were further generalized by Brandt and Seidel [30] to include rotating black holes with a superimposed gravitational wave. In this case, the data are no longer time-symmetric yet they satisfy a generalized form of the isometry condition so that the solution is still represented on two isometric, asymptotically flat hypersurfaces.

Matzner et al. [77] have begun to move beyond conformally flat initial data for binary black holes Their proposal is to use boosted versions of the Kerr metric written in the Kerr-Schild form to represent each black hole. Thus, an isolated black hole will have no spurious radiation content in the initial data. To construct solutions with multiple black holes, they propose, essentially, to use a linear combination of the single-hole solutions. The resulting metric can be used as the conformal 3-metric, the trace of the resulting extrinsic curvature can be used for K, and the tracefree part of the resulting extrinsic curvature can be used for ̃M^ij. Their scheme uses York’s conformal transverse-traceless decomposition outlined in Section 2.2.1, with the boundary conditions of ψ = 1 and Vⁱ = 0 on the horizons of the black holes and conditions appropriate for asymptotic flatness at large distances from the holes. A related method has been proposed by Bishop et al. [15], but their approach is much different and outside the current scope of this review.

The approach outlined by Matzner should certainly yield “cleaner” data than the conformally flat data currently available. For the task of specifying data for astrophysical black-hole binaries in nearly circular orbits, it is still true that these new data will not contain the correct initial gravitational wave content. Because the black holes are in orbit, they must be producing a continuous wave-train of gravitational radiation. This radiation will not be included in the method proposed by Matzner et al. Also, it is clear that the boundary conditions being used do not correctly account for the tidal distortion of each black hole by its companion. When the black holes are sufficiently far apart, the radiation from the orbital motion can be computed using post-Newtonian techniques. One possibility for producing astrophysically realistic, binary black-hole initial data is to use information from these post-Newtonian calculations to obtain better guesses for ̃γ_ij, ̃M^ij, and K.

3.3 Horizon-Penetrating Solutions

We noted in Section 3.1.1 that the time-independent maximal slicing of Schwarzschild with isotropic coordinates covers only the exterior of the black hole. This is because the time independence of this gauge requires that the lapse vanish on the horizon. It is possible to evolve into the black hole’s interior when starting from initial data constructed in this gauge, but it requires a choice for the lapse that yields a time-dependent solution [94]. The time dependence of such a solution is purely gauge, of course, since the spacetime is static.

It is possible to cover all, or part, of the interior of a single black hole with a time-independent slicing. However, doing so seems to require that we give up the maximal-slicing condition. To cover the interior of the black hole, we need a slicing that passes smoothly through the event horizon. A convenient way to generate such solutions is to begin with the metric in standard ingoing-null coordinates. If we want to consider a rotating and charged black hole, then we use the Kerr-Newman geometry in Kerr coordinates:

$$\begin{array}{l} d{s^2} = - \left( {1 - \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}} \right){\rm{d}}{{\tilde V}^2} + 2{\rm{d}}\tilde V{\rm{d}}r - 2\frac{{2Mr - {Q^2}}}{{{\rho ^2}}}a{\sin ^2}\theta \;{\rm{d}}\tilde V{\rm{d}}\tilde \phi \\ \;\;\;\;\;\;\;\; + {\rho ^2}{\rm{d}}{\theta ^2} - 2a{\sin ^2}\theta \;{\rm{d}}r{\rm{d}}\tilde \phi + \frac{1}{{{\rho ^2}}}\left[ {{{({r^2} + {a^2})}^2} - \Delta {a^2}{{\sin }^2}\theta } \right]{\sin ^2}\theta \;{\rm{d}}{{\tilde \phi }^2}, \end{array}$$

((83))

where ρ and Δ are defined by Eq. (79), and ̃V is the ingoing-null coordinate. This metric is regular at $r = {r_ \pm } \equiv M \pm \sqrt {{M^2} - {a^2} - {Q^2}}$, where r₊ and r_- are the locations of the event horizon and the Cauchy horizon, respectively.

This metric can be put into a form suitable for producing time-independent Cauchy initial data by making coordinate transformations of the general form

$${\rm{d}}t = {\rm{d}}\tilde V + f(r){\rm{d}}r,\:\;\;\;\;\;\;{\rm{d}}\phi = {\rm{d}}\tilde \phi + g(r){\rm{d}}r,$$

((84))

where f and g are suitably chosen functions of the radial coordinate r. There are a few particularly significant solutions for the general Kerr-Newman geometry, and I will outline these below, listing the nonzero components of the metric in the standard 3+1 format.

3.3.1 Kerr-Schild Coordinates

A spherical coordinate version of the standard Kerr-Schild coordinate system is obtained from (83) by using the coordinate choice (cf. Refs. [73, 42])

$$t = \tilde V - r\;\;\;\;{\rm{and}}\;\;\;\;\;\phi = \tilde \phi .$$

((85))

The nonzero components of the lapse, shift, and 3-metric are then given by:

$${\alpha ^{ - 2}} = 1 + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}},$$

((86))

$${\beta ^r} = {\alpha ^2}\frac{{2Mr - {Q^2}}}{{{\rho ^2}}},$$

((87))

$${\gamma _{rr}} = 1 + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}},$$

((88))

$${\gamma _{r\phi }} = - \left[ {1 + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}} \right]\;a{\sin ^2}\theta ,$$

((89))

$${\gamma _{\theta \theta }} = {\rho ^2},$$

((90))

$${\gamma _{\phi \phi }} = \left[ {{r^2} + {a^2} + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}{a^2}{{\sin }^2}\theta } \right]{\sin ^2}\theta .$$

((91))

Cartesian coordinate components can be obtained from these via the standard Kerr-Schild coordinate transformations [79]

$$x + iy = (r + ia){e^{i\phi }}\sin \theta \;\;\;\;{\rm{and}}\;\;\;\;z = r\cos \theta .$$

((92))

This yields the implicit definition of r from

$${r^4} - {r^2}({x^2} + {y^2} + {z^2} - {a^2}) - {a^2}{z^2} = 0,$$

((93))

with r > 0 and r = 0 on the disk described by z = 0 and x² + y² ≤ a².

3.3.2 Harmonic Coordinates

Harmonic time slicing is integral to some hyperbolic formulations of general relativity, and a time-independent harmonic slicing of the Kerr-Newman geometry does exist [17, 42]. The harmonic time slicing condition is □t = 0, which can be written

$$\frac{1}{{\sqrt { - g} }}{\partial _\mu }\left( {\sqrt { - g} {g^{0\mu }}} \right) = 0.$$

((94))

This equation is satisfied by using the coordinate choice

$$t = \tilde V - r + 2M\ln \left| {\frac{{2M}}{{r - {r_ - }}}} \right|,\:\;\;\;\;\;\;\phi = \tilde \phi .$$

((95))

The nonzero components of the lapse, shift, and 3-metric are then given by:

$${\alpha ^{ - 2}} = 1 + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}\left( {\frac{{r + {r_ + }}}{{r - {r_ - }}}} \right) + \frac{{r_ + ^2 + {a^2}}}{{{\rho ^2}}}\left( {\frac{{2M}}{{r - {r_ - }}}} \right),$$

((96))

$${\beta ^r} = {\alpha ^2}\frac{{r_ + ^2 + {a^2}}}{{{\rho ^2}}},$$

((97))

$${\beta ^\phi } = - {\alpha ^2}\frac{a}{{{\rho ^2}}}\left( {\frac{{2M}}{{r - {r_ - }}}} \right),$$

((98))

$${\gamma _{rr}} = \left[ {2 - \left( {1 - \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}} \right)\frac{{r + {r_ + }}}{{r - {r_ - }}}} \right]\frac{{r + {r_ + }}}{{r - {r_ - }}},$$

((99))

$${\gamma _{r\phi }} = - \left[ {1 + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}\left( {\frac{{r + {r_ + }}}{{r - {r_ - }}}} \right)} \right]\;a{\sin ^2}\theta ,$$

((100))

$${\gamma _{\theta \theta }} = {\rho ^2},$$

((101))

$${\gamma _{\phi \phi }} = \left[ {{r^2} + {a^2} + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}{a^2}{{\sin }^2}\theta } \right]\;{\sin ^2}\theta .$$

((102))

Cartesian coordinate components can be obtained from these via the standard Kerr-Schild coordinate transformations (92) and (93). However, for the harmonic slicing, the t = const. hypersurface is spacelike only outside the Cauchy horizon at r > r_-.

Fully harmonic coordinates (□x^μ = 0) can be defined when Cartesian spatial coordinates are used by employing a variation of the standard Kerr-Schild coordinate transformations [42]

$$x + iy = (r - m + ia){e^{i\phi }}\sin \theta \;\;\;\;{\rm{and}}\;\;\;\;z = (r - m)\cos \theta .$$

((103))

This yields the implicit definition of r from

$${(r - m)^4} - {(r - m)^2}({x^2} + {y^2} + {z^2} - {a^2}) - {a^2}{z^2} = 0.$$

((104))

Fully harmonic coordinates are useful because applying a boost to a harmonically sliced black hole yields a solution that satisfies (94) only if the black hole is written in fully harmonic coordinates. In this case, the boosted solution also satisfies the fully harmonic coordinate conditions.

3.3.3 Generalized Painlevé-Gullstrand Coordinates

The Painlevé-Gullstrand gauge choice for the Schwarzschild geometry has been rediscovered many times because of its simple form (cf. Refs. [85, 59, 88, 64, 67]). It is another time-independent solution, but the 3-geometry is completely flat (not simply conformally flat). The lapse is one in this gauge, and all of the information regarding the curvature of spacetime is contained in the shift. The Painlevé-Gullstrand gauge also has an intuitive physical interpretation [74]. An observer starting at rest at infinity and freely falling will trace out a world line that is everywhere orthogonal to the t = const. hypersurfaces in Painlevé-Gullstrand coordinates.

A generalization of the Painlevé-Gullstrand gauge derived by Doran [48] includes the Kerr spacetime, and the extension of this solution to the full Kerr-Newman spacetime is trivial. In the limit that a and Q vanish, this solution reduces to the Painlevé-Gullstrand gauge. The coordinate transformation is written most easily as

$${\rm{d}}t = {\rm{d}}\tilde V - \frac{{{\rm{d}}r}}{{1 + \sqrt {\frac{{2Mr - {Q^2}}}{{{r^2} + {a^2}}}} }},\;\;\;\;\;\;{\rm{d}}\phi = {\rm{d}}\tilde \phi - \frac{a}{{{r^2} + {a^2}}}\frac{{{\rm{d}}r}}{{1 + \sqrt {\frac{{2Mr - {Q^2}}}{{{r^2} + {a^2}}}} }}.$$

((105))

The nonzero components of the lapse, shift, and 3-metric are then given by:

$${\alpha ^{ - 2}} = 1,$$

((106))

$${\beta ^r} = {\alpha ^2}\sqrt {\frac{{{r^2} + {a^2}}}{{{\rho ^2}}}} \sqrt {\frac{{2Mr - {Q^2}}}{{{\rho ^2}}}} ,$$

((107))

$${\gamma _{rr}} = \frac{{{\rho ^2}}}{{{r^2} + {a^2}}},$$

((108))

$${\gamma _{r\phi }} = - \left[ {\sqrt {\frac{{{\rho ^2}}}{{{r^2} + {a^2}}}} \sqrt {\frac{{2Mr - {Q^2}}}{{{\rho ^2}}}} } \right]\;a{\sin ^2}\theta ,$$

((109))

)

$${\gamma _{\theta \theta }} = {\rho ^2},$$

((110))

$${\gamma _{\phi \phi }} = \left[ {{r^2} + {a^2} + \frac{{2Mr - {Q^2}}}{{{\rho ^2}}}{a^2}{{\sin }^2}\theta } \right]\;{\sin ^2}\theta .$$

((111))

Notice that the lapse remains one, but the 3-geometry is no longer flat when the black hole is spinning. Cartesian coordinate components can be obtained from these via the standard Kerr-Schild coordinate transformations (92) and (93). Like the Kerr-Schild time slicing, a t = const. slice of the generalized Painlevé-Gullstrand gauge remains spacelike for all r ≥ 0.

3.4 Quasicircular Binary Data

One of the primary driving forces behind the development of black-hole initial data has been the two-body problem of general relativity: the inspiral and coalescence of a pair of black holes. This problem is of fundamental importance. Not only is the relativistic two-body problem the most fundamental dynamical problem of general relativity, it is also considered one of the most likely candidates for observation with the upcoming generation of gravitational wave laser interferometers. Because of the circularizing effects of gravitational radiation damping, we expect the orbits of most tight binary systems to have small eccentricities. It is therefore desirable to have a method that can discern which data sets, within the very large parameter space of binary black-hole initial-data sets, correspond to black-hole binaries in a nearly circular (quasicircular) orbit.

Currently, only one approach has been developed for locating quasicircular orbits in a parameter space of binary black-hole initial data [39]. It is based on the fact that minimizing the energy of a binary system while keeping the orbital angular momentum fixed will yield a circular orbit in Newtonian gravity. The idea does not hold strictly for general relativistic binaries since they emit gravitational radiation and cannot be in equilibrium. However, for orbits outside the innermost stable circular orbit, the gravitational radiation reaction time scale is much longer than the orbital period. Thus it is a good approximation to treat such binaries as an equilibrium system. Called an “effective potential method”, this approach was used originally to find the quasicircular orbits and innermost stable circular orbit (ISCO) for equal-sized nonspinning black holes [39]. In this work, the initial data for binary black holes were computed using the conformal-imaging approach outlined in Section 3.2.1. The approach was also applied to binary black-hole data computed using the puncture method [10], where similar results were found. Configurations containing a pair of equal-sized black holes with spin also have been examined [86]. In this case, the spins of the black holes are equal in magnitude, but are aligned either parallel to, or anti-parallel to, the direction of the orbital angular momentum.

The approach defines an “effective potential” based on the binding energy of the binary. The binding energy is defined as

$${E_{\rm{b}}} \equiv {E_{{\rm{ADM}}}} - {M_1} - {M_2},$$

((112))

where E_ADM is the total ADM energy of the system measured at infinity, and M₁ and M₂ are the masses of the individual black holes. Quasicircular orbital configurations are obtained by minimizing the effective potential (defined as the nondimensional binding energy E_b/μ (where μ = M₁M₂/(M₁ + M₂)) as a function of separation, while keeping the ratio of the masses of the black holes M₁/M₂, the spins of the black holes S₁/M ² ₁ and S₂/M ² ₂ , and the total angular momentum J/(M₁M₂) constant.

This approach is limited primarily by the ambiguity in defining the individual masses of the black holes, M₁ and M₂. There is no rigorous definition for the mass of an individual hole in a binary configuration and some approximation must be made here. There is also no rigorous definition for the individual spins of holes in a binary configuration. The problem of defining the individual masses becomes particularly pronounced when the holes are very close together (see Ref. [86]). The limiting choice of conformal flatness for the 3-geometry also has proven to be problematic. The effects of this choice have been clearly seen in the case of quasicircular orbits of spinning black holes [86], but it is also believed to be a serious problem for any binary configuration because binary configurations are not conformally flat at the second post-Newtonian order [87].

To date, the results of the effective-potential method have not matched well to the best result from post-Newtonian approximations [47]¹¹. It will be interesting to see if the results from post-Newtonian approximations and numerical initial-data sets converge, especially when the approximation of conformal flatness is eliminated.

4 Neutron-Star Initial Data

The construction of initial data for neutron stars requires that the state of the neutron-star matter be specified before the gravitational data can be determined. Of course, a solution of the gravitational constraint equations can be found in principal for any given energy density ρ and momentum density jⁱ. But, with neutron-star solutions, we are usually interested in situations where the matter is in (or nearly in) hydrostatic equilibrium and the gravitational fields are also in (or nearly in) equilibrium.

4.1 Hydrostatic Equilibrium

For a neutron star to be in true equilibrium, the spacetime must be stationary as discussed in Section 2.4. This means that the spacetime possesses both “temporal” and “angular” Killing vectors (cf. Ref. [35]). If the matter is also to be in equilibrium, then the 4-velocity of the matter u^μ must be a linear combination of these two Killing vectors. If we use coordinates as defined in (53) with the angular Killing vector in the ϕ direction, then

$${u^\mu } = {u^t}[1,\;0,\;0,\:\Omega ].$$

((113))

Here, u^t and Ω are functions of r and θ only. Ω is the angular velocity of the matter as measured at infinity.

It is common to define v as the relative velocity between the matter and a normal observer (often called a zero angular momentum observer) so that

$$\frac{1}{{\sqrt {1 - {v^2}} }} = - {n_\mu }{u^\mu } = \alpha {u^t}.$$

((114))

The velocity v is then fixed by the normalization condition u^μu_μ = -1.

If we assume that the matter source is a perfect fluid, then the stress-energy tensor is given by

$${T^{\mu \nu }} = (\varepsilon + P){u^\mu }{u^\nu } + P{g^{\mu \nu }},$$

((115))

where ε and P are the total energy density and pressure, respectively, as measured in the rest frame of the fluid. The vanishing of the divergence of the stress-energy tensor yields the equation of hydrostatic equilibrium (often referred to as the relativistic Bernoulli equation). In differential form, this is

$${\rm{d}}P - (\varepsilon + P)({\rm{d}}\ln {u^t} - {u^t}{u_\phi }{\rm{d}}\Omega ) = 0.$$

((116))

If the fluid is barytropic¹², then we can define the relativistic enthalpy as

$$h(P) \equiv \exp \left[ {\int_0^P {\frac{{dP}}{{\varepsilon + P}}} } \right],$$

((117))

and rewrite the relativistic Bernoulli equation as

$$\ln h(P) = \ln {h_0} + \ln {u^t} - \ln u_0^t - \int_{{\Omega _0}}^\Omega {{u^t}} {u_\phi }d\Omega .$$

((118))

The constants h₀, u ^t₀ , and Ω₀ are the values their respective quantities have at some reference point, often taken to be the surface of the neutron star at the axis of rotation. When uniform rotation is assumed (dΩ = 0), Eq. (118) is rather easy to solve. The case of differential rotation is somewhat more complicated. An integrability condition of (116) requires that u^tu_ϕ be expressible as a function of Ω, so

$${u^t}{u_\phi } \equiv F(\Omega ).$$

((119))

F(Ω) is a specifiable function of Ω which determines the rotation law that the neutron star must obey [35].

4.2 Isolated Neutron Stars

The simplest models of isolated neutron stars are static (i. e., nonrotating), spherically symmetric models that can be constructed, given a suitable equation of state, by solving the Oppenheimer-Volkoff (OV) equations [84]:

$$\begin{array}{l} \frac{{dP}}{{dr}} = - \frac{{(\varepsilon + P)(m + 4\pi {r^3}P)}}{{r(r - 2m)}},\\ \frac{{dm}}{{dr}} = 4\pi {r^2}\varepsilon ,\\ \frac{{d\alpha }}{{dr}} = \frac{{\alpha (m + 4\pi {r^3}\varepsilon )}}{{r(r - 2m)}}, \end{array}$$

((120))

for 0 ≤ r ≤ R, where R is the radius of the surface of the star. Here, r is an areal radius and m(r) is the mass inside radius r. Exterior to the surface of the star, the metric is the standard Schwarzschild metric as in Eq. (59) with M = m(R). Interior to the surface of the star, the metric is

$${\rm{d}}{s^2} = - {\alpha ^2}{\rm{d}}{t^2} + {\left( {1 - \frac{{2m(r)}}{r}} \right)^{ - 1}}{\rm{d}}{r^2} + {r^2}({\rm{d}}{\theta ^2} + {\sin ^2}\theta {\rm{d}}{\phi ^2}).$$

((121))

The boundary conditions are that m(0) = 0, ε(0) is some chosen constant ε₀, and $\alpha (R) = \sqrt {1 - 2M/R}$. The solutions of this equation form a one-parameter family, parameterized by ε₀ which determines how relativistic the system is. A method for solving these equations in both the areal coordinate r and an isotropic radial coordinate ̃r can be found in Ref. [32].

More generally, isolated neutron stars will be rotating. If the neutron stars are uniformly rotating, then, for any given equation of state, the solutions form a two-parameter family. These models can be parameterized by their central density, which determines how relativistic they are, and by the amount of rotation. If the models are allowed to have differential rotation, then some rotation law must be chosen.

To construct a neutron-star model, the equations for a stationary solution of Einstein’s equations outlined in Section 2.4 must be solved self-consistently with the equations for hydrostatic equilibrium of the matter outlined above in Section 4.1. The equations that must be solved depend on the form of the metric chosen, and numerous formalisms and numerical schemes have been used. An incomplete list of references to work on constructing neutron-star models include [103, 22, 34, 32, 35, 33, 52, 62, 63, 43, 45, 44, 98, 21, 56, 16, 19]. Further review information on neutron-star models can be found in Refs. [97, 51, 50].

4.3 Neutron-Star Binaries

Neutron-star binaries, like any relativistic binary system, cannot exist in a true equilibrium configuration since they must emit gravitational radiation. But, as is true for black-hole binary data, for orbits outside the innermost stable circular orbit, the gravitational radiation reaction time scale is much longer than the orbital period and it is a reasonable approximation to consider the stars to be in a quasiequilibrium state.

A binary configuration obviously lacks the azimuthal symmetry that was assumed in the discussions of stationary solutions of Einstein’s equations in Section 2.4 and hydrostatic equilibrium in Section 4.1. Fortunately, the condition of hydrostatic equilibrium requires only the presence of a single, timelike Killing vector. With the assumption that gravitational radiation is negligible, we can assume that the matter is in some equilibrium state as viewed from the reference frame that is rotating along with the binary. That is, if the binary has a constant orbital angular velocity of Ω, then the time vector in the rotating frame t′^μ is a Killing vector and it is related to the time vector in the rest frame of the binary t^μ by

$${t'^\mu } = {t^\mu } + \Omega {\xi ^\mu },$$

((122))

where ξ^μ is a generator of rotations about the rotation axis¹³. Bonazzola et al. [18] refer to t*^μ as a helicoidal Killing vector.

Two equilibrium states for the matter have been explored in the literature. The simplest case is that of co-rotation, where the 4-velocity of the matter is proportional to t′^μ. In this case, the matter is at rest in the frame of reference rotating with the binary system, the corotating reference frame. The second equilibrium state is that of counter-rotation, where there is no rotation in the rest frame of the binary. We will explore these two cases further below.

Stationarity of the gravitational field, unlike hydrostatic equilibrium, requires the presence of separate timelike and azimuthal Killing vectors. For the case of a binary system, there is no unique definition of quasiequilibrium. The earliest results on constructing quasiequilibrium solutions of Einstein’s equations stem from work by Wilson and Mathews [104, 105, 106], and others have explored similar schemes [11, 18, 12]. Although written in slightly different forms, the system of equations for the gravitational fields in all of these schemes are fundamentally identical. While they were developed before the conformal thin-sandwich decomposition of Section 2.3, the conformal thin-sandwich decomposition (see Eq. (51)) offers the easiest way to interpret this approach. We consider ourselves to be in the corotating reference frame so that our time vector is t′^μ. To make the transition back to the rest frame of the binary as easy as possible, we write the shift vector of our 3+1 decomposition as

$${B^i} = {\beta ^i} + \Omega {\xi ^i},$$

((123))

so that βⁱ remains as the shift vector of the 3+1 decomposition made with respect to the rest frame of the binary system.

The primary assumptions are that the conformal 3-metric γ_ij is flat, the initial-data slice is maximal so that K = 0, and ũ^ij = 0. We see from (42) that the last assumption implies that the conformal 3-geometry is instantaneously stationary as seen in the corotating reference frame. The final choice that must be made is for the conformally rescaled lapse ̃μ. An elliptic equation for the lapse can be obtained by demanding that the trace of the extrinsic curvature K also be instantaneously stationary in the corotating reference frame. This is the so-called maximal slicing condition on the lapse. For the particular assumptions we have made here, this equation can be written

$${\tilde \nabla ^2}(\tilde \alpha {\psi ^7}) = (\tilde \alpha {\psi ^7})\left[ {\frac{7}{8}\psi {}^{ - 8}{{\tilde A}_{ij}}{{\tilde A}^{ij}} + 2\pi G{\psi ^4}(\rho + 2S)} \right].$$

((124))

It is interesting to note that, for a conformally flat 3-geometry,

$$({\tilde{\mathbb{L}}B)^{ij}} = ({\tilde{\mathbb{L}}\beta })^{ij} ,$$

((125))

so Ω does not appear in the equations for the gravitational fields except in the matter terms and possibly in boundary conditions. Equations (51) and (124) can be solved for the gravitational fields on an initial-data hypersurface, given values for the matter terms and appropriate boundary conditions.

If we choose the matter so that it is in hydrostatic equilibrium with respect to the pseudo-Killing vector t′^μ, then these equations for the gravitational fields will yield data that are in quasiequilibrium in the sense that ̃γ_ij and K are both instantaneously stationary.

For corotating binaries, the matter is at rest in the corotating reference frame of the binary. It is rigidly rotating and hydrostatic equilibrium is specified by solving the relativistic Bernoulli equation (118), with dΩ = 0, self-consistently with the equations for the gravitational fields.

For counterrotating binaries, the matter is not rotating in the rest frame of the binary. Counterrotating equilibrium configurations can be obtained by assuming the matter to have irrotational flow [99, 92]. As long as the flow is isentropic, we can express the enthalpy (117) as¹⁴

$$h = \frac{{\varepsilon + P}}{{{\rho _0}}},$$

((126))

where ρ₀ is the rest-mass density. For irrotational flow, the vorticity of the fluid (cf. Ref. [99]) is zero. Combining this with the Euler equation, we find that the 4-velocity of the fluid can be expressed as

$$h{u_\mu } = {\nabla _\mu }\varphi ,$$

((127))

where φ is the velocity potential, or flow field. Equation (127), together with the normalization condition u^μu_μ= -1, automatically satisfies the Euler equation and we are left with the continuity equation which must be satisfied,

$${\nabla _\mu }({\rho _0}{u^\mu }) = 0.$$

((128))

The continuity equation (128), the normalization condition, and Eq. (127) yield

$${\nabla ^2}\varphi = - ({\nabla ^\mu }\varphi ){\nabla _\mu }\ln ({\rho _0}/h)\;\;\;\;{\rm{with}}\;\;\;\;h = \sqrt { - ({\nabla ^\nu }\varphi )({\nabla _\nu }\varphi )} .$$

((129))

Stationarity (or quasistationarity) in the corotating reference frame requires that

$$h{u_\mu }{t^{\prime \mu }} = - C,$$

((130))

where C is a positive constant. Now, in terms of the 3+1 decomposition with Bⁱ as the shift vector (see Eq. (123)), we find that the Bernoulli equation is written as

$${h^2} = - ({\bar \nabla ^i}\varphi ){\bar \nabla _i}\varphi + \frac{1}{{{\alpha ^2}}}{(C + {B^j}{\bar \nabla _j}\varphi )^2}.$$

((131))

The flow field φ must satisfy

$$\begin{array}{*{20}{l}} {{{\bar \nabla }^2}\varphi - {B^i}{{\bar \nabla }_i}\left( {\frac{\lambda }{{{\alpha ^2}}}} \right) - \frac{\lambda }{\alpha }K = - \left( {{{\bar \nabla }^i}\varphi - \frac{\lambda }{{{\alpha ^2}}}{B^i}} \right){{\bar \nabla }_i}\ln \left( {\frac{{\alpha {\rho _0}}}{h}} \right)\,,}\\ {\;\:\:\:\:\:\:\:\:\:\:\:\:\:\;\;\;\;\;\;\;\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\;\lambda \equiv C + {B^i}{{\bar \nabla }_i}\varphi ,} \end{array}$$

((132))

subject to the boundary condition¹⁵ at the surface of the flow,

$$\left( {{{\bar \nabla }^i}\varphi - \frac{\lambda }{{{\alpha ^2}}}{B^i}} \right){\bar \nabla _i}\mathop {\left. {{\rho _0}} \right|}\nolimits_{{\rm{surf}}} = 0.$$

((133))

Solving Eqs. (131) and (132) self-consistently with the equations for the gravitational fields yields a counterrotating binary in hydrostatic equilibrium.

As mentioned above, the earliest work on neutron-star binaries was carried out by Wilson and Mathews [104, 105]. Wilson et al. [106] describe their approach for generating initial data for equilibrium neutron-star binaries. In these early works, the equation of hydrostatic equilibrium was not used. Rather, an initial guess for the density profile was chosen and the full hydrodynamic system was evolved with viscous damping until equilibrium was reached. During each step of the hydrodynamic evolution, the equations for the gravitational fields were resolved. The resulting data represented neither strictly co- nor counter-rotating binary neutron stars. This work led to the controversial result [75] that each neutron star in the binary may become radially unstable and collapse prior to the merger of the pair of stars. While an error was found in this work [49, 76] with the result that the signature of collapse is significantly weaker, the controversy has not yet been completely resolved.

The first use of corotating hydrostatic equilibrium with the Wilson-Mathews approach for specifying the gravitational fields was by Cook et al. [46] for the test case of an isolated neutron star. This approach was then used to study corotating neutron-star binaries by Baumgarte et al. [11, 12] and by Marronetti et al. [71]. Interestingly, turning-point methods for detecting secular instabilities [95, 96, 53] can be applied to the case of corotating binaries [13].

Corotating binary configurations are relatively easy to construct. However, it is believed that the viscosity of neutron-star matter is not large enough to allow for synchronization of the spin with the orbit [61, 14]. But if the initial spins of the neutron stars are not too large, close binaries should be well approximated by irrotational models. Bonazzola et al. [18] (as corrected by Asada [7]) developed the first approach for constructing counterrotating binary configurations. However, simpler formulations of irrotational flow were developed independently by Teukolsky [99] and Shibata [92], and Gourgoulhon [55] showed that all three approaches were equivalent. Numerical solutions of the equations for irrotational flow coupled to the equations for the gravitational fields are more difficult to construct than those for corotation because of the boundary condition (133) on the flow field that must be applied on the surface of each neutron star. This boundary condition is particularly difficult to implement because the location of the surface of the star is not known a priori, and will move as the equations are being solved. The first models of irrotational binary neutron stars were constructed by Bonazzola et al. [20], Marronetti et al. [72], and Uryū and collaborators [101, 102]. A description of the numerical methods used by Bonazzola et al. can be found in Refs. [57, 58].