INTRODUCTION
============

The physics of compact objects is entering a particularly exciting phase, as new instruments can now yield unprecedented observations. For example, there is evidence that the Rossi X-ray Timing Explorer has identified the innermost stable circular orbit around an accreting neutron star [@zsss98]. Also, the new generation of gravitational wave detectors under construction, including LIGO, VIRGO, GEO and TAMA, promise to detect, for the first time, gravitational radiation directly (see, e.g., [@t95]).

In order to learn from these observations (and, in the case of the gravitational wave detectors, to dramatically increase the likelihood of detection), one has to predict the observed signal from theoretical modeling. The most promising candidates for detection by the gravitational wave laser interferometers are the coalescences of black hole and neutron star binaries. Simulating such mergers requires self-consistent, numerical solutions to Einstein’s field equations in 3 spatial dimensions, which is extremely challenging. While several groups, including two “Grand Challenge Alliances” [@gc], have launched efforts to simulate the coalescence of compact objects (see also [@on97; @wmm96]), the problem is far from being solved.

Before Einstein’s field equations can be solved numerically, they have to be cast into a suitable initial value form. Most commonly, this is done via the standard 3+1 decomposition of Arnowitt, Deser and Misner (ADM, [@adm62]). In this formulation, the gravitational fields are described in terms of spatial quantities (the spatial metric and the extrinsic curvature), which satisfy some initial constraints and can then be integrated forward in time. The resulting “$\dot g - \dot
K$” equations are straightforward, but do not satisfy any known hyperbolicity condition, which, as it has been argued, may cause stability problems in numerical implementations. Therefore, several alternative, hyperbolic formulations of Einstein’s equations have been proposed [@fr94; @bmss95; @aacy96; @pe96; @f96; @acy98]. Most of these formulations, however, also have disadvantages. Several of them introduce a large number of new, first order variables, which take up large amounts of memory in numerical applications and require many additional equations. Some of these formulations require taking derivatives of the original equations, which may introduce further inaccuracies, in particular if matter sources are present. It has been widely debated if such hyperbolic formulations have computational advantages [@texas95]; their performance has yet to be compared directly with that of the original ADM equations. Accordingly, it is not yet clear if or how much the numerical behavior of the ADM equations suffers from their non-hyperbolicity.

In this paper, we demonstrate by means of a numerical experiment and a direct comparison that the standard implementation of the ADM system of equations, consisting of evolution equations for the bare metric and extrinsic curvature variables, is more susceptible to numerical instabilities than a modified form of the equations based on a conformal decomposition as suggested by Shibata and Nakamura [@sn95]. We will refer to the standard, “$\dot g - \dot
K$” form of the equations as “System I” (see Section \[sys1\] below). We follow Shibata and Nakamura and modify these original ADM equations by factoring out a conformal factor and introducing a spatial field of connection functions (“System II”, see Section \[sys2\] below). The conformal decomposition separates “radiative” variables from “nonradiative” ones in the spirit of the “York-Lichnerowicz” split [@l44; @y71]. With the help of the connection functions, the Ricci tensor becomes an elliptic operator acting on the components of the conformal metric. The evolution equations can therefore be reduced to a set of wave equations for the conformal metric components, which are coupled to the evolution equations for the connection functions. These wave equations reflect the hyperbolic nature of general relativity, and can also be implemented numerically in a straight-forward and stable manner.

We evolve low amplitude gravitational waves in pure vacuum spacetimes, and directly compare Systems I and II for both geodesic slicing and harmonic slicing. We find that System II is not only more appealing mathematically, but performs far better numerically than System I. In particular, we can evolve low amplitude waves in a stable fashion for hundreds of light travel timescales with System II, while the evolution crashes at an early time in System I, independent of gauge choice. We present these results in part to alert developers of 3+1 general relativity codes, many of whom currently employ System I, that a better set of equations may exist for numerical implementation.

The paper is organized as follows. In Section \[sec2\], we present the basic equations of both Systems I and II. We briefly discuss our numerical implementation in Section \[sec3\], and present numerical results in Section \[sec4\]. In Section \[sec5\], we summarize and discuss some of the implications of our findings.

BASIC EQUATIONS {#sec2}
===============

System I {#sys1}
--------

We write the metric in the form $$ds^2 = - \alpha^2 dt^2 + \gamma_{ij} (dx^i + \beta^i dt)(dx^j + \beta^j dt),$$ where $\alpha$ is the lapse function, $\beta^i$ is the shift vector, and $\gamma_{ij}$ is the spatial metric. Throughout this paper, Latin indices are spatial indices and run from 1 to 3, whereas Greek indices are spacetime indices and run from 0 to 3. The extrinsic curvature $K_{ij}$ can be defined by the equation $$\label{gdot1}
\frac{d}{dt} \gamma_{ij} = - 2 \alpha K_{ij},$$ where $$\frac{d}{dt} = \frac{\partial}{\partial t} - {\cal L}_{\beta}$$ and where ${\cal L}_{\beta}$ denotes the Lie derivative with respect to $\beta^i$.

The Einstein equations can then be split into the Hamiltonian constraint $$\label{ham1}
R - K_{ij}K^{ij} + K^2 = 2 \rho,$$ the momentum constraint $$\label{mom1}
D_j K^{j}_{~i} - D_i K = S_i,$$ and the evolution equation for the extrinsic curvature $$\label{Kdot1}
\frac{d}{dt} K_{ij} = - D_i D_j \alpha + \alpha ( R_{ij} 
        - 2 K_{il} K^l_{~j} + K K_{ij} - M_{ij} )$$ Here $D_i$ is the covariant derivative associated with $\gamma_{ij}$, $R_{ij}$ is the three-dimensional Ricci tensor $$\begin{aligned}
 \label{ricci}
R_{ij} & = & \frac{1}{2} \gamma^{kl} 
    \Big( \gamma_{kj,il} + \gamma_{il,kj} 
        - \gamma_{kl,ij} - \gamma_{ij,kl} \Big) \\[1mm]
    & & + \gamma^{kl} \Big( \Gamma^m_{il} \Gamma_{mkj}
    - \Gamma^m_{ij} \Gamma_{mkl} \Big), \nonumber\end{aligned}$$ and $R$ is its trace $R = \gamma^{ij} R_{ij}$. We have also introduced the matter sources $\rho$, $S_i$ and $S_{ij}$, which are projections of the stress-energy tensor with respect to the unit normal vector $n_{\alpha}$, $$\begin{aligned}
\rho & = & n_{\alpha} n_{\beta} T^{\alpha \beta}, \nonumber \\[1mm]
S_i  & = & - \gamma_{i\alpha} n_{\beta} T^{\alpha \beta}, \\[1mm]
S_{ij} & = & \gamma_{i \alpha} \gamma_{j \beta} T^{\alpha \beta}, \nonumber\end{aligned}$$ and have abbreviated $$M_{ij} \equiv S_{ij} + \frac{1}{2} \gamma_{ij}(\rho - S),$$ where $S$ is the trace of $S_{ij}$, $S = \gamma^{ij} S_{ij}$.

The evolution equations (\[gdot1\]) and (\[Kdot1\]) together with the constraint equations (\[ham1\]) and (\[mom1\]) are equivalent to the Einstein equations, and are commonly referred to as the ADM form of the gravitational field equations [@adm62; @footnote1]. We will call these equations System I. This system is widely used in numerical relativity calculations (e.g. [@aetal98; @cetal98]), even though its mathematical structure is not simple to characterize and may not be ideal for computation. In particular, the Ricci tensor (\[ricci\]) is not an elliptic operator: while the last one of the four terms involving second derivatives, $\gamma^{kl}\gamma_{ij,kl}$, is an elliptic operator acting on the components of the metric, the elliptic nature of the whole operator is spoiled by the other three terms involving second derivatives. Accordingly, the system as a whole does not satify any known hyperbolicity condition (see also the discussion in [@f96]). Therefore, to establish existence and uniqueness of solutions to Einstein’s equations, most mathematical analyses rely either on particular coordinate choices or on different formulations.

System II {#sys2}
---------

Instead of evolving the metric $\gamma_{ij}$ and the extrinsic curvature $K_{ij}$, we can evolve a conformal factor and the trace of the extrinsic curvature separately (“York-Lichnerowicz split” [@l44; @y71]). Such a split is very appealing from both a theoretical and computational point of view, and has been widely applied in numerical axisymmetric (2+1) calculations (see, e.g., [@e84]). More recently, Shibata and Nakamura [@sn95] applied a similar technique in a three-dimensional (3+1) calculation. Adopting their notation, we write the conformal metric as $$\tilde \gamma_{ij} = e^{- 4 \phi} \gamma_{ij}$$ and choose $$e^{4 \phi} = \gamma^{1/3} \equiv \det(\gamma_{ij})^{1/3},$$ so that the determinant of $\tilde \gamma_{ij}$ is unity. We also write the trace-free part of the extrinsic curvature $K_{ij}$ as $$A_{ij} = K_{ij} - \frac{1}{3} \gamma_{ij} K,$$ where $K = \gamma^{ij} K_{ij}$. It turns out to be convenient to introduce $$\tilde A_{ij} = e^{- 4 \phi} A_{ij}.$$ We will raise and lower indices of $\tilde A_{ij}$ with the conformal metric $\tilde \gamma_{ij}$, so that $\tilde A^{ij} = e^{4 \phi} A^{ij}$ (see [@sn95]).

Taking the trace of the evolution equations (\[gdot1\]) and (\[Kdot1\]) with respect to the physical metric $\gamma_{ij}$, we find [@footnote2] $$\label{phidot2}
\frac{d}{dt} \phi = - \frac{1}{6} \alpha K$$ and $$\label{Kdot2}
\frac{d}{dt} K = - \gamma^{ij} D_j D_i \alpha + 
    \alpha(\tilde A_{ij} \tilde A^{ij}
    + \frac{1}{3} K^2) + \frac{1}{2} \alpha (\rho + S),$$ where we have used the Hamiltonian constraint (\[ham1\]) to eliminate the Ricci scalar from the last equation. The tracefree parts of the two evolution equations yield $$\label{gdot2}
\frac{d}{dt} \tilde \gamma_{ij} =  
    - 2 \alpha \tilde A_{ij}.$$ and $$\begin{aligned}
 \label{Adot2}
\frac{d}{dt} \tilde A_{ij} & = & e^{- 4 \phi} \left( 
    - ( D_i D_j \alpha )^{TF}  +
    \alpha ( R_{ij}^{TF} - S_{ij}^{TF} ) \right)
    \nonumber \\[1mm]
    & & + \alpha (K \tilde A_{ij} - 2 \tilde A_{il} \tilde A^l_{~j}).\end{aligned}$$ In the last equation, the superscript $TF$ denotes the trace-free part of a tensor, e.g. $R_{ij}^{TF} = R_{ij} - \gamma_{ij} R/3$. Note that the trace $R$ could again be eliminated with the Hamiltonian constraint (\[ham1\]). Note also that $\tilde \gamma_{ij}$ and $\tilde A_{ij}$ are tensor densities of weight $-2/3$, so that their Lie derivative is, for example, $${\cal L}_{\beta} \tilde A_{ij} = \beta^k \partial_k \tilde A_{ij}
    + \tilde A_{ik} \partial_j \beta^k 
    + \tilde A_{kj} \partial_i \beta^k
    - \frac{2}{3} \tilde A_{ij} \partial_k \beta^k.$$

The Ricci tensor $R_{ij}$ in (\[Adot2\]) can be written as the sum $$R_{ij} = \tilde R_{ij} + R_{ij}^{\phi}.$$ Here $R_{ij}^{\phi}$ is $$\begin{aligned}
R^{\phi}_{ij} & = & - 2 \tilde D_i \tilde D_j \phi - 
    2 \tilde \gamma_{ij} \tilde D^l \tilde D_l \phi \nonumber \\[1mm]
    & & + 4 (\tilde D_i \phi)(\tilde D_j \phi)
    - 4 \tilde \gamma_{ij} (\tilde D^l \phi) (\tilde D_l \phi),\end{aligned}$$ where $\tilde D_i$ is the derivative operator associated with $\tilde
\gamma_{ij}$, and $\tilde D^i = \tilde \gamma^{ij} \tilde D_j$.

The “tilde” Ricci tensor $\tilde R_{ij}$ is the Ricci tensor associated with $\tilde \gamma_{ij}$, and could be computed by inserting $\tilde \gamma_{ij}$ into equation (\[ricci\]). However, we can bring the Ricci tensor into a manifestly elliptic form by introducing the “conformal connection functions” $$\label{cgsf}
\tilde \Gamma^i \equiv \tilde \gamma^{jk} \tilde \Gamma^{i}_{jk}
    = - \tilde \gamma^{ij}_{~~,j},$$ where the $\tilde \Gamma^{i}_{jk}$ are the connection coefficients associated with $\tilde \gamma_{ij}$, and where the last equality holds because $\tilde \gamma = 1$. In terms of these, the Ricci tensor can be written [@footnote3] $$\begin{aligned}
\tilde R_{ij} & = & - \frac{1}{2} \tilde \gamma^{lm}
    \tilde \gamma_{ij,lm} 
    + \tilde \gamma_{k(i} \partial_{j)} \tilde \Gamma^k
    + \tilde \Gamma^k \tilde \Gamma_{(ij)k}  + \nonumber \\[1mm]
    & & \tilde \gamma^{lm} \left( 2 \tilde \Gamma^k_{l(i} 
    \tilde \Gamma_{j)km} + \tilde \Gamma^k_{im} \tilde \Gamma_{klj} 
    \right).\end{aligned}$$ The principal part of this operator, $\tilde \gamma^{lm} \tilde
\gamma_{ij,lm}$, is that of a Laplace operator acting on the components of the metric $\tilde \gamma_{ij}$. It is obviously elliptic and diagonally dominant (as long as the metric is diagonally dominant). All the other second derivatives of the metric appearing in (\[ricci\]) have been absorbed in the derivatives of the connection functions. At least in appropriately chosen coordinate systems (for example $\beta^i = 0$), equations (\[gdot2\]) and (\[Adot2\]) therefore reduce to a coupled set of nonlinear, inhomogeneous wave equations for the conformal metric $\tilde
\gamma_{ij}$, in which the gauge terms $K$ and $\tilde \Gamma^i$, the conformal factor $\exp(\phi)$, and the matter terms $M_{ij}$ appear as sources. Wave equations not only reflect the hyperbolic nature of general relativity, but can also be implemented numerically in a straight-forward and stable manner. The same method has often been used to reduce the four-dimensional Ricci tensor $R_{\alpha\beta}$ [@c62] and to bring Einstein’s equations into a symmetric hyperbolic form [@fm72].

Note that the connection functions $\tilde \Gamma^i$ are pure gauge quantities in the sense that they could be chosen, for example, to vanish by a suitable choice of spatial coordinates (“conformal three-harmonic coordinates”, compare [@sy78]). The $\tilde
\Gamma^i$ would then play the role of “conformal gauge source functions” (compare [@c62; @fm72]). Here, however, we impose the gauge by choosing the shift $\beta^i$, and evolve the $\tilde
\Gamma^i$ with equation (\[Gammadot2\]) below. Similarly, $K$ is a pure gauge variable (and could be chosen to vanish by imposing maximal time slicing).

An evolution equation for the $\tilde \Gamma^i$ can be derived by permuting a time derivative with the space derivative in (\[cgsf\]) $$\begin{aligned}
\frac{\partial}{\partial t} \tilde \Gamma^i
& = & - \frac{\partial}{\partial x^j} \Big( 2 \alpha \tilde A^{ij} 
    \nonumber \\[1mm]
& & - 2 \tilde \gamma^{m(j} \beta^{i)}_{~,m}
    + \frac{2}{3} \tilde \gamma^{ij} \beta^l_{~,l} 
    + \beta^l \tilde \gamma^{ij}_{~~,l} \Big).\end{aligned}$$ It turns out to be essential for the numerical stability of the system to eliminate the divergence of $\tilde A^{ij}$ with the help of the momentum constraint (\[mom1\]), which yields $$\begin{aligned}
 \label{Gammadot2}
\frac{\partial}{\partial t} \tilde \Gamma^i
& = & - 2 \tilde A^{ij} \alpha_{,j} + 2 \alpha \Big(
    \tilde \Gamma^i_{jk} \tilde A^{kj} -  \nonumber \\[1mm]
& & \frac{2}{3} \tilde \gamma^{ij} K_{,j}
    - \tilde \gamma^{ij} S_j + 6 \tilde A^{ij} \phi_{,j} \Big) +
    \nonumber \\[1mm]
& & \frac{\partial}{\partial x^j} \Big(
    \beta^l \tilde \gamma^{ij}_{~~,l} 
    - 2 \tilde \gamma^{m(j} \beta^{i)}_{~,m}
    + \frac{2}{3} \tilde \gamma^{ij} \beta^l_{~,l} \Big).\end{aligned}$$

We now consider $\phi$, $K$, $\tilde \gamma_{ij}$, $\tilde A_{ij}$ and $\tilde \Gamma^i$ as fundamental variables. These can be evolved with the evolution equations (\[phidot2\]), (\[Kdot2\]), (\[gdot2\]), (\[Adot2\]), and (\[Gammadot2\]), which we call System II. Note that obviouly not all these variables are independent. In particular, the determinant of $\tilde \gamma_{ij}$ has to be unity, and the trace of $\tilde A_{ij}$ has to vanish. These conditions can either be used to reduce the number of evolved quantities, or, alternatively, all quantities can be evolved and the conditions can be used as a numerical check (which is what we do in our implementation).

NUMERICAL IMPLEMENTATION {#sec3}
========================

In order to compare the properties of Systems I and II, we implemented them numerically in an identical environment. We integrate the evolution equations with a two-level, iterative Crank-Nicholson method. The iteration is truncated after a certain accuracy has been achieved. However, we iterate at least twice, so that the scheme is second order accurate.

The gridpoints on the outer boundaries are updated with a Sommerfeld condition. We assume that, on the outer boundaries, the fundamental variables behave like outgoing, radial waves $$Q(t,r) = \frac{G(\alpha t - e^{2 \phi} r)}{r}.$$ Here $Q$ is any of the fundamental variables (except for the diagonal components of $\tilde \gamma_{ij}$, for which the radiative part is $Q
= \tilde \gamma_{ii} - 1$), and $G$ can be found by following the characteristic back to the previous timestep and interpolating the corresponding variable to that point (see also [@sn95]). We found that a linear interpolation is adequate for our purposes.

We impose octant symmetry in order to minimize the number of gridpoints, and impose corresponding symmetry boundary conditions on the symmetry plains. Unless noted otherwise, the calculations presented in this paper were performed on grids of $(32)^3$ gridpoints, and used a Courant factor of $1/4$. The code has been implemented in a parallel environment on SGI Power ChallengeArray and SGI CRAY Origin2000 computer systems at NCSA using DAGH [@dagh] software for parallel processing.

RESULTS {#sec4}
=======

Initial Data
------------

For initial data, we choose a linearized wave solution (which is then evolved with the full nonlinear systems I and II). Following Teukolsky [@t82], we construct a time-symmetric, even-parity $L=2$, $M = 0$ solution. The coefficients $A, B$ and $C$ (see equation (6) in [@t82]) are derived from a function $$F(t,r) = {\cal A}\, (t \pm r) \, \exp(- \lambda (t \pm r)^2).$$ Unless noted otherwise, we present results for an amplitude ${\cal A} = 10^{-3}$ and a wavelength $\lambda = 1$. The outer boundary conditions are imposed at $x, y, z = 4$.

We evolve these initial data for zero shift $$\beta^i = 0,$$ and compare the performance of Systems I and II for both geodesic and harmonic slicing.

Geodesic Slicing
----------------

In geodesic slicing, the lapse is unity $$\alpha = 1.$$ Since the acceleration of normal observers satisfies $a_a = D_a \ln \alpha = 0$, these observers follow geodesics. The energy content of even a small, linear wave packet will therefore focus these observers, and even after the wave has dispersed, the observers will continue to coast towards each other. Since $\beta^i = 0$, normal observers are identical to coordinate observers, hence geodesic slicing will ultimately lead to the formation of a coordinate singularity even for arbitrarily small waves.

The timescale for the formation of this singularity can be estimated from equation (\[Kdot2\]) with $\alpha = 1$ and $\beta^i = 0$. The $\tilde A_{ij}$, which can be associated with the gravitational waves, will cause $K$ to increase to some finite value, say $K_0$ at time $t_0$, even if $K$ was zero initially. After roughly a light crossing time, the waves will have dispersed, and the further evolution of $K$ is described by $\partial_t K \sim K^2/3$, or $$\label{K_approx}
K \sim \frac{3 K_0}{3 - K_0(t - t_0)}$$ (see [@sn95]). Obviously, the coordinate singularity forms at $t \sim 3/K_0 + t_0$ as a result of the nonlinear evolution.

We can now evolve the wave initial data with Systems I and II and compare how well they reproduce the formation of the coordinate singularity.

In Figure 1, we show $K$ at the origin ($x = y = z = 0$) as a function of time both for System I (dashed line) and System II (solid line). We also plot the approximate analytic solution (\[K\_approx\]) as a dotted line, which we have matched to the System I solution with values $K_0 = 0.00518$ and $t_0 = 10$. For these values, equation (\[K\_approx\]) predicts that the coordinate singularity appears at $t \sim 590$. In the insert, we show a blow-up of System II for early times. It can be seen very clearly how the initial wave content lets $K$ grow from zero to the “seed” value $K_0$. Once the waves have dispersed, System II approximately follows the solution (\[K\_approx\]) up to fairly late times. System I, on the other hand, crashes long before the coordinate singularity appears.

In Figure 2, we compare the extrinsic curvature component $K_{zz}$ evaluated at the origin. The noise around $t \sim 8$, which is present in the evolutions of both systems, is caused by reflections of the initial wave off the outer boundaries. It is obvious from these plots that System II evolves the equations stably to a fairly late time, at which the integration eventually becomes inaccurate as the coordinate singularity approaches. We stopped this calculation when the iterative Crank-Nicholson scheme no longer converged after a certain maximum number of iterations. It is also obvious that System I performs extremely poorly, and crashes at a very early time, well before the coordinate singularity.

It is important to realize that the poor performance of System I is [*not*]{} an artifact of our numerical implementation. For example, the ADM code currently being used by the Black Hole Grand Challenge Alliance, is based on the equations of System I, and also crashes after a very similar time [@r98] (see also [@aetal98], where a run with a much smaller initial amplitude nevertheless crashes earlier than our System II). This shows that the code’s crashing is intrinsic to the equations and slicing, and not to our numerical implementation.

Harmonic Slicing
----------------

Since geodesic slicing is known to develop coordinate singularities for generic, nontrivial initial data, it is obviously not a very good slicing condition. We therefore also compare the two Systems using harmonic slicing. In harmonic slicing, the coordinate time $t$ is a harmonic function of the coordinates $\nabla^{\alpha} \nabla_{\alpha} t = 0$, which is equivalent to the condition $$\Gamma^0 \equiv g^{\alpha\beta} \Gamma^0_{\alpha\beta} = 0$$ (where the $\Gamma^{\alpha}_{\beta\gamma}$ are the connection coefficients associated with the four-dimensional metric $g_{\alpha\beta}$). For $\beta^i = 0$, the above condition reduces to $$\partial_t \alpha = - \alpha^2 K.$$ Inserting (\[phidot2\]), this can be written as $$\partial_t (\alpha e^{-6 \phi}) = 0
\mbox{~~~or~~~}
\alpha = C(x^i) e^{6 \phi},$$ where $C(x^i)$ is a constant of integration, which depends on the spatial coordinates only. In practice, we choose $C(x^i) = 1$.

In Figure 3, we show results for the same initial data as in the last section. Obviously, both Systems do much better for this slicing condition. System I crashes much later than in geodesic slicing (after about 40 light crossing times, as opposed to about 10 for geodesic slicing), but it still crashes. System II, on the other hand, did not crash after even over 100 light crossing times. We never encountered a growing instability that caused the code to crash.

SUMMARY AND CONCLUSION {#sec5}
======================

We numerically implement two different formulations of Einstein’s field equations and compare their performance for the evolution of linear wave initial data. System I is the standard set of ADM equations for the evolution of $\gamma_{ij}$ and $K_{ij}$. In System II, we conformally decompose the equations and introduce connection functions. The conformal decomposition naturally splits “radiative” variables from “nonradiative” ones, and the connection functions are used to bring the Ricci tensor into an elliptic form. These changes are appealing mathematically, but also have a striking numerical consequence: System II performs far better than System I.

It is interesting to note that most earlier axisymmetric codes (e.g. [@e84]) also relied on a decomposition similar to that of System II. Much care was taken to identify radiative variables and to integrate those variables as opposed to the raw metric components. It is surprising that this experience was abandoned in the development of most 3+1 codes, which integrate equations equivalent to System I. These codes have been partly successful [@cetal98], but obvious problems remain, as for example the inability to integrate low amplitude waves for arbitrarily long times. While efforts have been undertaken to stabilize such codes with the help of appropriate outer boundary conditions [@aetal98; @rar98], our findings point to the equations themselves as the fundamental cause of the problem, and not to the outer boundaries. Obviously, boundary conditions as employed in the perturbative approach in [@aetal98; @rar98] or in the characteristic approach in [@bglswi96] are still needed for accuracy – but our results clearly suggest that they are not needed for stability [@footnote4].

Some of the recently proposed hyperbolic systems are very appealing in that they bring the equations in a first order, symmetric hyperbolic form, and that all characteristics are physical (i.e., are either at rest with respect to normal observers or travel with the speed of light) [@aacy96; @acy98]. These properties may be very advantageous for numerical implementations, in particular at the boundaries (both outer boundaries and, in the case of black hole evolutions, inner “apparent horizon” boundaries). Some of these systems have also been implemented numerically, and show stability properties very similar to our System II [@cs98]. Our System II, on the other hand, uses fewer variables than most of the hyperbolic formulations, and does not take derivatives of the equations, which may be advantageous especially when matter sources are present. This suggests that a system similar to System II may be a good choice for evolving interior solutions and matter sources, while one may want to match to one of the hyperbolic formulations for a better treatment of the boundaries.

The mathematical structure of System II is more appealing than that of System I, and these improvements are reflected in the numerical behavior. We therefore conclude that the mathematical structure has a very deep impact on the numerical behavior, and that the ability to finite difference the standard “$\dot g - \dot K$” ADM equations may not be sufficient to warrant a stable evolution.

It is a pleasure to thank A. M. Abrahams, L. Rezzolla, J. W. York and M. Shibata for many very useful conversations. We would also like to thank H. Friedrich for very valuable comments, and S. A. Hughes for a careful checking of our code. Calculations were performed on SGI CRAY Origin2000 computer systems at the National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign. This work was supported by NSF Grant AST 96-18524 and NASA Grant NAG 5-3420 at Illinois, and the NSF Binary Black Hole Grand Challenge Grant Nos. NSF PHYS 93-18152, NSF PHY 93-10083 and ASC 93-18152 (ARPA supplemented).

W. Zhang, Z. P. Smale, T. E. Stohmayer and J. H. Swank, Astrophys. J. [**500**]{}, L171 (1998).

K. Thorne, in [*Proceedings of the Seventeenth Texas Symposium on Relativistic Astrophysics and Cosmology*]{}, edited by H. Böhringer, G. E. Morfill and J. E. Trümper (Annals of the New York Academy of Sciences, Vol. 759, New York, 1995).

Information about the Binary Black Hole Grand Challenge can be found at [www.npac.syr.edu/projects/bh/]{}, and about the Binary Neutron Star Grand Challenge at [wugrav.wustl.edu/Relativ/nsgc.html]{}

K. Oohara and T. Nakamura, in [*Relativistic Gravitation and Gravitational Radiation*]{}, edited by J.-A. Marck and J.-P. Lasota (Cambridge University Press, Cambridge, 1997).

J. R. Wilson and G. J. Mathews, Phys. Rev. Lett. [**75**]{}, 4161 (1995); J. R. Wilson, G. J. Mathews and P. Marronetti, Phys. Rev. D [**54**]{}, 1317 (1996).

R. Arnowitt, S. Deser and C. W. Misner, in [*Gravitation: An Introduction to Current Research*]{}, edited by L. Witten (Wiley, New York, 1962).

S. Frittelli and O. Reula, Commun. Math. Phys. [**166**]{}, 221 (1994).

C. Bona, J. Massó, E. Seidel and J. Stela, Phys. Rev. Lett. [**75**]{}, 600 (1995).

A. Abrahams, A. Anderson, Y. Choquet-Bruhat and J. W. York, Jr., Phys. Rev. Lett. [**75**]{}, 3377 (1996).

M. H. P. M. van Putten and D. M. Eardley, Phys. Rev. D [**53**]{}, 3056 (1996).

H. Friedrich, Class. Quantum Gravit. [**13**]{}, 1451 (1996).

A. Anderson, Y. Choquet-Bruhat and J. W. York, Jr., to appear in Topol. Methods in Nonlinear Analysis (also gr-qc/9710041).

For example at the Third Texas Workshop on 3-dimensional Numerical Relativity of the Binary Black Hole Grand Challenge, held at Austin, Texas, 1995 (Proceedings can be obtained from Richard Matzner, Center for Relativity, University of Texas at Austin, Texas).

M. Shibata and T. Nakamura, Phys. Rev. D [**52**]{}, 5428 (1995).

A. Lichnerowicz, J. Math. Pure Appl. [**23**]{}, 37 (1944).

J. W. York, Jr., Phys. Rev. Lett. [**26**]{}, 1656 (1971).

Note, however, that Arnowitt, Deser and Misner [@adm62] wrote these equations in terms of the conjugate momenta $\pi_{ij}$ instead of the extrinsic curvature $K_{ij}$.

A. M. Abrahams [*et. al.*]{} (The Binary Black Hole Grand Challenge Alliance) Phys. Rev. Lett. [**80**]{}, 1812 (1998);

G. B. Cook [*et. al.*]{} (The Binary Black Hole Grand Challenge Alliance) Phys. Rev. Lett. [**80**]{}, 2512 (1998).

For example: J. M. Bardeen and T. Piran, Phys. Rep. [**96**]{}, 205 (1983); C. R. Evans, PhD thesis, University of Texas at Austin (1984); A. M. Abrahams and C. R. Evans, Phys. Rev. D [**37**]{}, 318 (1988); A. M. Abrahams, G. B. Cook, S. L. Shapiro and S. A. Teukolsky, Phys. Rev. D [**49**]{}, 5153 (1994).

Note also that ${\cal L}_{\beta} \phi = \beta^i \partial_i \phi +
    \partial_i \beta^i / 6$.

Shibata and Nakamura [@sn95] use a similar auxiliary variable $F_i = \tilde \gamma_{ij,j}$ to eliminate some second derivatives from the Ricci tensor.

T. De Donder, [*La gravifique einsteinienne*]{} (Gauthier-Villars, Paris, 1921); C. Lanczos, Phys. Z. [**23**]{}, 537 (1922); Y. Choquet-Bruhat, in [*Gravitation: An Introduction to Current Research*]{}, edited by L. Witten (Wiley, New York, 1962).

A. Fischer and J. Marsden, Comm. Math. Phys., [**28**]{}, 1 (1972).

L. Smarr and J. W. York, Jr., Phys. Rev. D [**17**]{}, 1945 (1978).

M. Parashar and J. C. Brown, in [*Proceedings of the International Conference for High Performance Computing*]{}, eds. S. Sahni, V. K. Prasanna and V. P. Bhatkar (Tata McGraw-Hill, New York, 1995), also [www.caip.rutgers.edu/$\sim$parashar/DAGH/]{}

S. A. Teukolsky, Phys. Rev. D [**26**]{}, 745 (1982).

L. Rezzolla, talk presented at the Binary Black Hole Grand Challenge Alliance’s meeting at the University of Pittsburgh, April 1998.

M. E. Rupright, A. M. Abrahams and L. Rezzolla, Phys. Rev. D [**58**]{}, 044005 (1998); L. Rezzolla, A. M. Abrahams, R. A. Matzner, M. Rupright and S. L. Shapiro, 1998, submitted.

N. Bishop, R. Gomez, L. Lehner, B. Szilagyi, J. Winicour and R. Isaacson, Phys. Rev. Lett. [**76**]{}, 4303 (1996).

Alternatively, the outer boundary conditions can be completely removed by a conformal rescaling; see, for example, P. Hübner, Phys. Rev. D [**53**]{}, 701 (1996) and P. Hübner, 1998, submitted (also gr-qc/9804065).

G. B. Cook, M. S. Scheel, private communication.
