Numerical Evolution of Black Holes with a Hyperbolic Formulation of General Relativity

We describe a numerical code that solves Einstein’s equations for a Schwarzschild black hole in spherical symmetry, using a hyperbolic formulation introduced by Choquet-Bruhat and York. This is the ﬁrst time this formulation has been used to evolve a numerical spacetime containing a black hole. We excise the hole from the computational grid in order to avoid the central singularity. We describe in detail a causal diﬀerencing method that should allow one to stably evolve a hyperbolic system of equations in three spatial dimensions with an arbitrary shift vector, to second-order accuracy in both space and time. We demonstrate the success of this method in the spherically symmetric case.

We describe a numerical code that solves Einstein's equations for a Schwarzschild black hole in spherical symmetry, using a hyperbolic formulation introduced by Choquet-Bruhat and York. This is the first time this formulation has been used to evolve a numerical spacetime containing a black hole. We excise the hole from the computational grid in order to avoid the central singularity. We describe in detail a causal differencing method that should allow one to stably evolve a hyperbolic system of equations in three spatial dimensions with an arbitrary shift vector, to second-order accuracy in both space and time. We demonstrate the success of this method in the spherically symmetric case.

I. INTRODUCTION
A key goal of numerical relativity is to determine the gravitational radiation produced by the inspiral and coalescence of two black holes in a decaying binary orbit. The importance of solving this problem is heightened by the possibility that gravitational wave detectors such as LIGO, VIRGO, and GEO may observe gravitational waveforms from binary black hole coalescence within the next decade. Not only could a comparison of measured waveforms with numerical simulations provide a crucial strong-field test of general relativity, but in addition, accurate templates produced by these simulations could significantly increase the sensitivity of waveform measurements [1] and reduce the uncertainties in astrophysical parameters derived from gravitational wave data.
Although constructing a numerical simulation of binary black hole coalescence is a difficult and unsolved problem, several recent advances have brought us closer to a solution. One key advance is the development of so-called apparent horizon boundary condition (AHBC) methods [2][3][4][5][6][7], which treat black holes by evolving only the regions of spacetime that lie outside apparent horizons. These methods take advantage of the fact that information cannot emerge from within the apparent horizon of a black hole (which, assuming cosmic censorship, is contained within the event horizon). Without AHBC schemes, the spacetime singularity that inevitably forms inside a black hole eventually causes numerical simulations to terminate, typically on a time scale of order 10 to 100M , where M is the mass of the system. For this reason, AHBC methods may be crucial for solving the binary black hole problem, where one hopes to evolve two holes long enough to see them orbit, coalesce, and eventually settle into a final state containing a single (Kerr) black hole.
Another promising development is the construction of manifestly hyperbolic formulations of Einstein's equations [8][9][10][11][12][13]. Not only do these formulations offer insights into the mathematical structure of general relativity, but they may also be better-suited for numerical solution than the usual ADM [14] equations, which are not manifestly hyperbolic. The reason for this is twofold: First, there exists an extensive literature concerning stable and efficient numerical methods for solving hyperbolic systems of equations [15]. Some of these methods have been successfully applied to hyperbolic formulations of general relativity by Bona and Masso [6]. The second reason is that a non-hyperbolic set of equations can present a fundamental difficulty for black hole simulations that employ an AHBC approach.
To understand this difficulty, consider a non-hyperbolic set of equations, or even a hyperbolic set that has characteristics lying outside the local light cone. Suppose such a system is to be solved on a domain that includes a black hole. Although the physics guarantees that nothing can emerge from the hole, the equations do not know this, and nonphysical (gauge) modes can propagate outward through the apparent horizon. Solving such a system of equations on a restricted domain that excludes the interior of the hole is mathematically well-posed only if appropriate boundary conditions are imposed on the horizon. While it should be possible to impose explicit horizon boundary conditions to fix the coordinate system [16], it is unclear which boundary conditions are appropriate for dynamical variables, particularly in the general three-dimensional case in which one may have a nonspherical horizon and a significant amount of gravitational radiation. Now consider a hyperbolic set of equations with characteristics that never lie outside the local light cone. In this case, future-pointing characteristics inside the apparent horizon (which must be an outgoing non-timelike surface) can never intersect the horizon itself, so that quantities at or outside the horizon cannot depend on the interior region. Consequently, one can solve this set of equations on a restricted domain that excludes the interior of the hole, and one can do so without imposing boundary conditions on the horizon.
In this paper, we concentrate on the hyperbolic formulation of Einstein's equations originally proposed by Choquet-Bruhat and York [9,10], hereafter referred to as the CBY system. This formulation has several advantageous features. First, the characteristics of this set of equations are extremely simple: they lie either along the light cone or along the normal to the current time slice. This guarantees that no information, not even gauge information, can propagate acausally. Second, the equations admit an arbitrary shift vector, and the characteristics are independent of the choice of shift. Finally, the fundamental variables in this formulation are spatially covariant three-dimensional tensors. These tensors directly measure spacetime curvature, and from them one can form all components of the spacetime Weyl tensor.
While we believe that the CBY formalism holds considerable promise for numerical simulations, particularly when combined with AHBC methods, there is no experience with solving this particular set of equations on a computer. Therefore, before expending the significant effort required to implement the CBY equations in a full three-dimensional code, it is important to demonstrate that such an approach is feasible. Accordingly, we have developed a numerical code that solves the spherically symmetric Einstein equations using the CBY formalism. We evolve a Schwarzschild spacetime using a causal differencing scheme, and we use an AHBC method to avoid the central singularity. This is the first time the CBY equations have been used to evolve a numerical spacetime containing a black hole.
Our code runs in parallel, and is rigorously second order convergent. As described in detail elsewhere [17], it can run for times in excess of 1000M provided certain constraints are regularly enforced. Our code is based on the DAGH [18] software package originally developed for the Binary Black Hole Grand Challenge Alliance. It is written in C++, and uses Fortran-90 numerical kernels. The DAGH system contains support for adaptive mesh refinement, but we have not yet taken advantage of this feature.
Our code provides an important demonstration that the CBY formulation works well in numerical simulations. It allows us to study the details of implementing the CBY equations in a simple setting, and it provides a testing ground for apparent horizon excising schemes and causal differencing algorithms. It also serves as an important check on a code that solves the CBY equations in three spatial dimensions, a code that is currently under development and will be described elsewhere.
We employ results and methods specific to spherical symmetry as little as possible, so that our techniques are readily generalizable to the three-dimensional case with Cartesian coordinates. For example, we do not use logarithmic radial coordinates [2,4] despite the great advantages they provide for spherically symmetric codes. Likewise, although there are many shift conditions that work well with AHBC methods in spherical symmetry [3], we consider only those that can be applied in the general three-dimensional case.
Recent progress has also been made in spherical symmetry by choosing a global coordinate system that has desirable properties near the horizon [7]. The spatial gauge used by [7] depends on the concept of an areal radius, and is thus applicable only in spherical symmetry. We forego such a coordinate choice in favor of a more general approach.
In section II, we summarize the CBY formulation of Einstein's equations, and we specialize this formulation to the case of spherical symmetry. In section III, we present a key ingredient of our code: a causal differencing scheme that is second-order convergent and has a stability criterion that is independent of the shift vector. We describe this scheme in detail for the general case of three spatial dimensions, and then apply it to the spherically symmetric case. In section IV we describe the AHBC method employed at the inner boundary of our grid, and the shift conditions that we use in order to implement this method. These shift conditions not only allow us to control the motion of the apparent horizon through the grid, but they also prevent coordinate singularities that may result from differential stretching or compression of the grid in the remainder of the spacetime. In section V we discuss the boundary condition that we impose at the outer boundary of our grid. This boundary is ideally at spatial infinity, but in practice it is placed at a large but finite radius. In section VI we present the results of rigorous convergence tests using our code. In section VII, we close with a short discussion of our results.

A. The CBY Formalism
Here we summarize the fundamental equations and variables used in the CBY representation of general relativity. For details of the CBY formulation and a derivation of the equations, see [9].
We write the metric in the usual 3+1 form where N is the lapse function, β i is the shift vector, and g ij is the three-metric on a spatial hypersurface of constant t.
Define the variables Here D is the three-dimensional covariant derivative compatible with the three-metric g ij , the time derivative operator is∂ and £ denotes a Lie derivative. The quantity K ij is the usual extrinsic curvature. If we assume that the time coordinate satisfies the harmonic slicing condition then the vacuum evolution equations take the form where H is the trace of K ij , X (ij) ≡ (X ij + X ji )/2, the quantities Γ k ij are the connection coefficients associated with the spatial covariant derivative D k , and Eq. (5d) is the harmonic slicing condition (4). There are considerably more variables and equations in the CBY formalism than in the usual ADM formalism. However, the form of the equations is much simpler in the CBY case. While the right-hand sides of Eqs. (5) contain many terms, these terms consist solely of algebraic combinations of the dynamical variables and involve no derivatives. Eqs. (5f-5i) are tensor wave equations whose characteristics are along the light cones. Eqs. (5a-5e) are even simplerthey drag information normal to the surfaces of constant t, that is, along zero-velocity (with respect to the normal) characteristics. There are no other characteristics in the system.
Although we have eliminated some gauge freedom in the choice of lapse function by imposing the harmonic time slicing condition (4), the shift β i is unspecified and completely arbitrary. The shift is not a dynamical variable in this formalism, in the sense that it obeys no evolution equation, and that it appears in Eqs. (5) only through the time derivative operator∂ 0 . Instead, the shift is an auxiliary gauge variable that may be freely chosen on each time slice, and may even change discontinuously from one slice to the next.
In vacuum, the constraint equations include whereR ij is the three-dimensional Ricci tensor. Eqs. (6a) and (6b) are the familiar Hamiltonian and momentum constraints rewritten in terms of the CBY variables, and Eq. (6c) is a result of harmonic time slicing. Eq. (6d) is the familiar ADM evolution equation for K ij , which in the CBY picture becomes a constraint on L ij . Eqs. (2c), (2d), (2f), and the usual relation between Γ k ij and derivatives of g ij are also constraints that must be satisfied at all times. All constraints are preserved by the evolution equations.

B. Spherical Symmetry
We write the spherically symmetric three-metric in the general form where (r, θ, φ) are the usual spherical coordinates. Spherical symmetry reduces the number of dynamical variables (including the connection coefficients) from 67 to 16. Define where the subscript T denotes "transverse". Then the vacuum evolution equations (5) take the form where The constraints (6) become The additional constraints (2c), (2d), (2f), and the usual relation between Γ k ij and derivatives of g ij take the form

III. CAUSAL EVOLUTION METHOD
Here we present the causal differencing method we use to evolve Eqs. (5) from one spacelike hypersurface to the next. Straightforward differencing schemes typically become unstable for large shifts, which are needed for the implementation of AHBC methods. Our method is second-order accurate and has a stability criterion that is independent of the shift vector. We emphasize that our method is not specific to Eqs. (5), but can be used to handle advective terms in any system of first-order evolution equations. Similar causal differencing methods have been used by [2,3] and [19] in order to treat large shifts in the standard ADM formulation of Einstein's equations.
Our causal differencing method is independent of the actual form of the shift vector. We only place two restrictions on the shift: First, it must be a smooth functional of the dynamical variables and of the space and time coordinates. Second, it can be computed to second-order accuracy, given second-order values for these variables and coordinates. The particular prescription that we use for computing the shift will be discussed in section IV. In this section we assume only that such a prescription exists.
Although the code described in this paper assumes spherical symmetry, our causal evolution method is general. In this section we first describe our method for the general case of three spatial dimensions plus time, and then we specialize to spherical symmetry.
A. Overview of the Method Figure 1 shows an initial spatial hypersurface labeled by t = t 0 , and a subsequent spatial hypersurface labeled by t = t 0 + ∆t. We wish to evolve quantities defined at the point (t 0 , x i ) to the point (t 0 + ∆t, x i ), that is, along the vector t shown in Figure 1. Here t and x i are the coordinates defined in Eq. (1), and the vector t is given by where n is the unit normal to the hypersurface t = t 0 . A large shift vector β tends to cause stability problems in most numerical schemes. Some schemes, including many implicit ones, are unconditionally unstable whenever t is non-timelike, as is the case in Figure 1. Other schemes can be made stable for an arbitrary shift, but only at the expense of a very small time step ∆t.
Spacetime diagram illustrating the relation t = t + β and showing both the (t, x i ) and the (t,x i ) coordinate systems. The vector t must always lie within the light cone. This is not true for the vector t for a sufficiently large shift β.
In order to construct a differencing scheme that works for an arbitrary shift, we introduce an auxiliary coordinate system. First define a new timelike vector Then define new coordinates (t, and such that the spatial coordinatesx i coincide with x i at t = t 0 . The new coordinates and their relationship to the vectors t and t are shown in Figure 1. Partial derivatives with respect to the new coordinates (t,x i ) are given by Our method works by breaking up each time step into two substeps, as illustrated in Figure 2: First, we evolve quantities along the vector t , that is, we evolve using the (t,x i ) coordinate system, from the points on the slice t = t 0 in Figure 2 to the points on the slice t = t 0 + ∆t that are labeled by circles. We then complete the timestep by interpolating from the points labeled by circles to those labeled by dots. These two substeps will be considered separately in sections III C and III D below.
B. Transforming into the (t,x i ) system Each evolution equation in the system (5) can be written in the form Two-dimensional spacetime diagrams illustrating the two coordinate systems used in our causal differencing scheme. Solid and dotted arrows denote the normal vector t and the time vector t at each grid point on the initial time slice. Solid dots represent grid points with particular values of x i , and circles represent grid points with particular values ofx i . The dots and circles coincide at t = t0, but not at t = t0 + ∆t. This is because x i is constant along t whilex i is constant along t . Diagrams A and B represent the same spacetime. The only difference is that A is drawn in the computational coordinate system (t, x i ), where the time axis lies along t, and B is drawn in the normal coordinate system (t,x i ), where the time axis lies along t and is normal to the spatial slices.∂ where T , the quantity being evolved, is not necessarily a scalar, but may be a coordinate-dependent object. We wish to rewrite Eq. (19), which is defined in the computational coordinate system (t, x i ), in terms of the normal coordinate system (t,x i ).
If we consider all quantities to be defined in the (t, x i ) basis, we can rewrite the∂ 0 operator in the (t,x i ) coordinate system as follows:∂ so that Eq. (19) becomes whereL Here we have used Eqs. (3), (13), and (17). The first two lines of Eq. (20) are coordinate independent, but in the third line we have assumed the (t, x i ) basis in order to write £ t = ∂/∂t . In the fourth line we have separated the Lie derivative along β into two pieces: The advective piece, β i ∂/∂x i , is the one responsible for the instability that often arises when one tries to evolve along t with a large shift vector. This piece is eventually absorbed into the time derivative ∂/∂t. The remaining piece,L β , when operating on some quantity T , describes the change in T induced by the change in basis vectors along β. The operatorL β vanishes when operating on a scalar. Furthermore,L β T does not actually contain any derivatives of T , but only contains derivatives of β. Therefore, no spatial derivatives of T appear in Eq. (21).
Note that Eq. (21) is not the same as Eq. (19) transformed into the (t,x i ) basis. By splitting the Lie derivative along β into two pieces, we have derived an equation that describes the evolution of T defined with respect to the (t, x i ) basis along a path of constantx i . The coordinate system used to define tensor components, (t, x i ), is different from the coordinate system used to label spacetime points during the evolution, (t,x i ).
Note also that we have introduced an additional auxiliary variable: the Jacobian ∂x j /∂x i that appears in Equation (21). From Eqs. (17) and (18), we can find the rate of change of the Jacobian along the vector t : To compute the Jacobian, we evolve this equation along with Eq. (21). Because we setx i = x i at the beginning of each time step t = t 0 , the initial value of ∂x j /∂x i on each time step is the Kroeneker delta δ j i . Finally, we require derivatives of β, which appear in the operatorL β and also on the right hand side of Eq. (24). Assuming that we have some prescription for choosing the shift given the values of x i , t, and the dynamical variables at each grid point, we simply use this prescription to compute β, and then we obtain its derivatives either analytically (if the shift is analytic) or by finite differencing. C.
Step 1: Evolve along t The first step of our causal differencing method is to evolve Eq. (21) and the auxiliary equation (24) along the vector t , fromt = t 0 tot = t 0 + ∆t. Although this can be done using any standard differencing scheme, the algorithm described in this section assumes a scheme with two time levels; a three-level scheme (such as leapfrog) requires a slight modification of the algorithm. We use a Macormack predictor-corrector method, which we illustrate with a simple wave equation in spherical symmetry, written in first-order form in the (t,r) coordinate system: Here R and S are arbitrary functions of P , Q,r, andt. Given a discrete set of uniformly spaced grid pointsr i , we denote P and Q at grid point i and timet =t n by P n i and Q n i respectively. To compute P and Q at timẽ t =t n+1 =t n + ∆t, we first compute initial guessesP andQ using the "predictor" step where ∆r =r i+1 −r i . The quantitiesP n+1 i andQ n+1 i are first order accurate in both space and time. Notice that the finite difference approximation to the spatial derivative is one-sided.
Once we have the predicted quantities, we then compute P and Q at timet =t n+1 =t n + ∆t to second order in both space and time using the "corrector" step The spatial derivatives are taken in the opposite direction to the predictor step. This ensures that the first order error terms introduced by the one-sided derivatives in the corrector step cancel those produced by the one-sided derivatives in the predictor step. This cancellation would be spoiled by substituting Q n+1 forQ n+1 in the spatial derivative term of Eq. (29) or by substituting P n+1 forP n+1 in the spatial derivative term of Eq. (30). However, it is irrelevant for accuracy or stability whether the right hand sides R and S in the corrector step are computed using predicted values of P and Q as in Eq. (29) or using corrected values of P and predicted values of Q as in Eq. (30). We use corrected values in the right-hand sides whenever possible because it minimizes the memory required for storing temporary variables on the computer. The above Macormack scheme is stable (disregarding the boundaries, which will be discussed in section V) whenever ∆t < ∆r. For a wave equation with characteristic speed v, the stability condition for the above scheme is the familiar Courant condition v∆t < ∆r.
To obtain a corrected value for a particular variable, we require predicted values for all quantities appearing in the equation for that variable. For our prototypical equation (21), in order to compute a corrected value for T , we must have predicted values not only for T , Q, S i , and R, but also for the Jacobian ∂x j /∂x i and for β and its derivatives (which appear in the operatorL β ). The predicted value of the Jacobian is obtained by evolving Eq. (24) to first order using the Macormack predictor step. The predicted value of the shift is obtained by using our shift prescription to compute β from the predicted values of the dynamical variables.

D. Step 2: Interpolation
One entire numerical time step should take variables defined on discrete points with particular values of x i , and compute quantities at the same values of x i but at a later time. In Figure 2, this corresponds to taking values defined at the discrete points on the slice t = t n and computing quantities at the solid dots on the slice t = t n + ∆t. However, when we evolve along t as described in section III C, we compute quantities at the points that correspond to the circles on the slice t = t n + ∆t in Figure 2.
Our causal differencing method therefore requires a second step, namely, interpolation of our dynamical variables from the (t,x i ) coordinate system back into the (t, x i ) coordinate system, that is, from the circled points to the dotted points on the slice t = t n + ∆t in Figure 2. The values of x i at the dotted points and the values ofx i at the circled points are known; both are the same as the values of x i at the appropriate grid points on the initial slice. To perform the interpolation, we must also know either the values of x i at the circled points or the values ofx i at the dotted ones.
From Eq. (17), the change in the coordinate x i along the vector t is given by Therefore, if we evolve Eq. (31) along with Eq. (21) using the Macormack scheme, we obtain the value of x i at each of the circled grid points at t = t n + ∆t in Figure 2. This allows us to interpolate quantities from the circles to the dots, working in the x i coordinate system. Note, however, that the circled grid points are, in general, not uniformly spaced in x i , as shown in Figure 2A. Instead, the dotted grid points are uniformly spaced in x i . While this poses no difficulty for the spherically symmetric case discussed in this paper, in the general three-dimensional case interpolating from an arbitrary set of points onto a uniform grid is a nontrivial numerical problem that cannot be treated very efficiently. Much easier and much less costly in terms of computer time is to interpolate from a uniform grid to an arbitrary set of points.
To handle this difficulty, notice that if we evolve along the vector t instead of along the vector t , the coordinates x i remain constant and the coordinatesx i vary. The change inx i along the vector t is given by where we have used Eq. (17). Therefore, if we evolve Eq. (32) along the vector t using the Macormack scheme, we obtain the value ofx i at each of the dotted grid points at t = t n + ∆t in Figure 2. This allows us to interpolate quantities from the circled points to the dotted ones, working in thex i coordinate system. The circled points are uniformly spaced inx i , as shown in Figure 2B. We thus interpolate from a uniform grid (inx i ) to an arbitrary set of points. This can be done relatively easily in three spatial dimensions. There is a subtlety in evolving Eq. (32) along the vector t using the Macormack scheme: In order to implement the corrector step of Eq. (32), we require predicted values of the quantities β j ∂x i /∂x j at the dotted points, but these quantities are only known at the circled points. However, from the predictor step of Eq. (32), we already have predicted values ofx i at the dotted points. Therefore, we use these values to interpolate the predicted values of β j ∂x i /∂x j from the circled points to the dotted ones. In this case, we are again interpolating from a uniform grid to an arbitrary set of points.

E. Implementation in Spherical Symmetry
In spherical symmetry, we solve equations (9) on a numerical grid that extends from some value r = r min just inside the apparent horizon to a large radius r = r max . We denote our numerical grid points by r i , where i runs from i min to i max , corresponding to the innermost and outermost grid points. We use a zone-centered grid, so that the innermost and outermost grid points do not correspond to r min and r max . Instead, we set where ∆r = r max − r min i max − i min + 1 .
When evolving along t , the spherically symmetric vacuum evolution equations (9) take the form ∂ ∂t a r = N a 0r + a r ∂r ∂r where the right hand sides (9q-9s) are unchanged, and we have included theL β terms explicitly. Note the second derivative of β r in the equation for Γ r rr . This term results from applying £ β to a non-tensorial quantity. Because β r is not an unknown in the system of equations (35), but is instead an auxiliary variable that may be, for example, given analytically as a function of the coordinates, this second derivative should not spoil the hyperbolicity of the system. In spherical symmetry, Equations (24) Eqs. (36) provide values for ∂r/∂r to be used in the corrector step, as well as coordinate information for the interpolation step. We use cubic interpolation for causal differencing. This is accurate to fourth order in ∆r. Quadratic interpolation would also suffice, but linear interpolation would not yield a scheme that was second order convergent in time. The reason is that linear interpolation makes an O(∆r) 2 error on each time step, so that after N time steps the error is of order ∆r. This is because for a fixed total time t total that we wish to evolve, the Courant condition requires N ∼ t total (∆r) −1 .

IV. SHIFT VECTOR AND INNER BOUNDARY CONDITION
In this section we describe how we choose a shift vector β, and how this choice affects how we handle the inner boundary of our computational domain.
Because both the CBY formalism and the causal differencing scheme discussed in section III place no restrictions on β, we are free to choose any shift we wish. Although setting β = 0 is the simplest choice, it is often useful to employ a nonzero shift vector. One technique is to use the shift to simplify the form of the Einstein equations, for example, to eliminate particular components of the spatial metric [20]. The disadvantage of this approach is that it involves an actual change in the equations being solved. Instead of Eqs. (5), one would be solving a different (but physically equivalent) system of equations that would include the shift as a dynamical variable, and would no longer be hyperbolic.
We instead use the shift for a different purpose: to allow us to truncate our computational domain just inside the apparent horizon, so that we evolve only the exterior region. We thus avoid the spacetime singularity inside the black hole. If we were to attempt such a truncation with β = 0, numerical grid points originally located just outside the apparent horizon would soon fall in, and any grid points located inside the apparent horizon would eventually encounter the singularity. This is because for zero shift, numerical grid points follow the world lines of normal observers, and these world lines are necessarily timelike.

A. Shift at the Inner Boundary
We force the inner boundary of our grid, r = r min , to hover within a grid spacing of the apparent horizon. To accomplish this for a static black hole, we choose β near the horizon to point along the outward normal to the horizon, and we choose its magnitude so that the local coordinate speed of light in that direction is zero. The horizon is then approximately stationary with respect to the spatial coordinates. For the spherically symmetric case, the local coordinate speed of light in the outgoing radial direction is given by so we set β r at the horizon equal to N/A. We track the apparent horizon at each time step, and retain only the grid points that lie on the outside. This is done via a masking algorithm that labels grid points outside the apparent horizon as valid, and those inside as invalid. Invalid points are never used in the computation, just as if those points did not exist. Because we use a cell-centered grid, the inner boundary r = r min is located half a grid spacing inside the innermost valid grid point.
Numerical errors typically cause the horizon to drift through the grid, even if one tries to lock it in place by forcing β r = N/A. While this drift seems to cause minimal difficulty with either the stability or accuracy of the code, for coarse-resolution runs it produces a small but distracting gauge pulse each time the horizon crosses a grid zone. Horizon drift can be eliminated by introducing a feedback mechanism [21] that adjusts the magnitude of the shift to compensate. We do this as follows: If we wish to force the horizon to remain at some radius r 0 , but we find the horizon is actually at radius r AH , then we set where β r , N , and A refer to values at the horizon location r AH . This feedback mechanism is not necessary for sufficiently fine grid spacing. While locating an apparent horizon on a numerically generated time slice is a difficult problem in multidimensions [22][23][24][25], it is trivial in spherical symmetry. The marginal outer trapped surface equation where s i is the spatial unit normal to the surface, reduces to for a spherically symmetric system. Here we have used our variables defined in Eqs. (7) and (8). We find the apparent horizon by first evaluating ϑ(r) at each grid point, and then by using 3-point interpolation to locate the outermost root that satisfies ϑ ′ (r) > 0. When locating the apparent horizon after the Macormack predictor step, we must be sure to use the value of r and notr in evaluating the function ϑ(r) at each grid point. In principle, one can also use a shift condition at the apparent horizon to move a black hole through the numerical grid. This could be accomplished by adjusting the shift so that grid points on one side of the hole fall into the horizon, and grid points on the other side emerge from it. A similar idea could in principle be applied to rotating black holes or systems with more than one black hole.
B. Inner Boundary Condition on Eqs. (9) Treating the inner boundary correctly is a primary motivation for using a hyperbolic formulation of the Einstein equations. The spherically symmetric CBY equations have only three characteristics at each spacetime point. These lie along the ingoing and outgoing null rays, and along the normal vector t . A boundary condition is required only at a point that cannot obtain information from one or more of the characteristics passing through it. For example, the outer boundary cannot obtain information from the ingoing null characteristics that it intersects, because these characteristics originate from outside the computational domain, where we have no data. To update quantities at the outer boundary, this information must be provided by a boundary condition. Similarly, the inner boundary ordinarily requires a boundary condition because it cannot obtain information from the outgoing null characteristics that it intersects. However, if the inner boundary follows an outgoing spacelike or null trajectory, then each point on the inner boundary can obtain information from all three characteristics passing through it, so a boundary condition is not required. This is the essence of an AHBC scheme for treating the inner boundary: by forcing the inner boundary to move along with the apparent horizon, we force it along an outgoing non-timelike path. Therefore, there is no need to impose an explicit inner boundary condition. Regardless of the mathematical formulation of Einstein's equations being used, general relativity tells us that when the inner boundary follows an outgoing spacelike or null path, information with physical content cannot penetrate this boundary from the inside, since this information cannot propagate outside the light cone. A key advantage of a hyperbolic formulation with only simple (non-spacelike) characteristics is that in such a formulation, this statement applies to gauge information as well.
The way we solve Eqs. (9) without imposing an explicit condition at the inner boundary is by simply ignoring the innermost grid point during the interpolation step of our causal differencing scheme. In the case where the inner boundary r = r min moves with respect to the (t,r) coordinate a distance less than ∆r during each time step, the interpolation becomes an extrapolation at the innermost point. This is always the case in our simulation because of the Courant limit: Because the inner boundary, which moves along with the apparent horizon, has a velocity with respect to the (t,r) coordinate system of approximately N/A, the inner boundary can never move farther than ∆r during a time step without violating the Courant condition ∆t < (A/N )∆r. One could avoid this extrapolation by using an implicit differencing scheme [19,4] to get around the Courant condition, but such a scheme requires much more computer time than explicit schemes, especially in the multidimensional case.

C. Shift in the Remainder of the Spacetime
Once a shift criterion at the apparent horizon has been chosen, one must then determine the shift in the remainder of the spacetime. At spatial infinity, one presumably would like the shift to approach zero, so that the spacetime metric components approach Minkowski values. Or perhaps, in the case of spacetimes with nonzero angular momentum, one would like the asymptotic shift to describe a co-rotating frame. However, given the shift at infinity and at the apparent horizon, it is not clear how to choose the shift elsewhere.
One possibility is to use a parameterized analytic function whose parameters are set so that the shift behaves appropriately near the apparent horizon and far from the black hole. For example, when evolving a single black hole in spherical symmetry we have tried the Gaussian form where C, r c , and w are chosen so that β r is equal to N/A at the apparent horizon, ∂β r /∂r is equal to ∂(N/A)/∂r at the apparent horizon, and β r is smaller than some threshold at the outer boundary of the grid. Although this choice results in a second order convergent evolution, we find that the grid points tend to compress or stretch where the shift gradients are large, and eventually coordinate singularities develop that cause the simulation to terminate. Similarly, one can choose so that the shift is equal to N/A far inside some arbitrary radius r = r c , and zero far outside r = r c . In this case, the grid points become compressed near r = r c as one evolves in time, and again the simulation terminates. It therefore appears that in addition to a prescription for specifying the shift at the apparent horizon and at infinity, one must impose some additional restriction on the shift that ensures that it will not induce any coordinate pathologies elsewhere. Such a restriction can be provided by two different elliptic shift conditions that were introduced for the very purpose of minimizing coordinate strain caused by a shift vector. The first is the minimal distortion condition [26], which can be written in the form The minimal distortion shift minimizes the average change of shape of a spatial volume element as it is dragged from one time slice to the next. A related choice is the minimal strain condition [26], which can be written in the form The minimal strain shift minimizes the average change in the three metric g ij as one evolves from one time slice to the next. It differs from the minimal distortion shift in that it takes into account the change in size of spatial volume elements as well as their change in shape. The downside of these shift conditions is that they require one to solve elliptic equations. This can be costly in terms of computer time, especially in three dimensions. It may be possible to use a parameterized analytic function to mimic one of these conditions, or it may suffice to use an approximate solution. However, since it appears that these conditions give us a useful shift vector, we will adopt them in the spherically symmetric case, where the computational burden is not so severe.
Both the minimal distortion and minimal strain conditions work well in spherical symmetry, as shown in section VI. They prevent grid points from becoming locally compressed or stretched to the point where coordinate singularities form. Using the variables defined in (2) and (8), we find that the minimal distortion equation (43) We have written both equations in ther coordinate system and included the factors ∂r/∂r and ∂ 2r /∂r 2 because the shift must be computed not only after the corrector step of the Macormack scheme whenr = r, but also after the predictor step, whenr = r. Eqs. (45) and (46) require boundary conditions at both ends of the numerical grid. We impose the condition ∂r ∂r at the outer boundary r = r max so that the spacetime is asymptotically Minkowski, and we impose at the apparent horizon, so that the apparent horizon is stationary with respect to the coordinates. Both the minimal distortion and minimal strain equations can be written in the general form We solve this equation using the usual three-point finite difference approximation In order to retain second-order accuracy, we must be careful always to impose the boundary condition (47) at the point r = r max , which may be different from the outer boundary of the grid,r = r max . We write Eq. (47) in the second-order accurate finite-difference form where λ ≡r max −r i ∆r .
Herer max is the value ofr corresponding to the point r = r max , and Eqs. (51) and (52) are evaluated at i = i max . Combining Eqs. (50) and (51) to eliminate the fictitious grid point i max + 1 yields a boundary condition for β r imax in terms of β r imax−1 . To impose the boundary condition (48) at the horizon, we use the finite-difference expression where λ ≡r AH −r i ∆r .
Herer AH is ther-coordinate location of the apparent horizon, and Eqs. (53) and (54) are evaluated at the value of i such that i − 1 is the innermost grid point that lies outside the horizon. Combining Eqs. (50) and (53) to eliminate the grid point at i + 1 yields a boundary condition for β r i−1 in terms of β r i . We currently use two methods to solve the linear system of equations that results from Eqs. (50), (51), and (53). We use a standard tridiagonal algorithm [27] when running in serial, or when running in parallel using a small number of processors. Although very efficient, this is a serial algorithm, so it begins to dominate the computation time as the number of processors increases. For a large number of processors, we instead use a multigrid technique [27,28] which is parallelized with the help of the DAGH system.

V. OUTER BOUNDARY
The only variables that require a boundary condition at the outer boundary r = r max are the quantities L rr , L T , M rrr , M rT , a 0r , and a rr , which propagate along the ingoing and outgoing null characteristics. All other variables evolve along the characteristics parallel to the vector t , and are never differentiated with respect tor in Eqs. (35).
To apply a boundary condition on these variables, we assume that for large r, the function f for which we need a boundary condition behaves like for some integer n and constants C j . We then impose r ∂ 2 ∂r 2 (r n f ) + 2 ∂ ∂r (r n f ) = 0 (56) at the outer boundary. This boundary condition truncates the series (55) after the second term, that is, it uses the correct values of C 1 and C 2 , but sets C j to zero for j > 2. We set the value of n for each variable according to its analytic falloff rate for the static solution.
We call this boundary condition an "extended Robin" method because of its similarity to the familiar Robin boundary condition which one often imposes on a function that behaves like Because Eq. (56) involves approximating the large-r behavior of the function f , it introduces errors that are no longer second-order in the grid discretization, but fall off rapidly as one increases the radius of the outer boundary. These errors result solely from the fact that Eq. (56) is not exactly satisfied by the initial data. To obtain a secondorder convergent scheme, we effectively eliminate these approximation errors by increasing r max until the second-order error terms dominate.
Because we use a cell-centered grid, the outer boundary is located at half a grid spacing beyond the outermost grid point: r max = r imax +1/2 . Therefore, to ensure second-order convergence in ∆r, we must impose Eq. (56) not at the outermost grid point r imax , but at the outer boundary r max . This is necessary because r imax depends on the grid spacing ∆r.
We impose Eq. (56) by using the following finite-difference approximations for the first and second derivatives: Substituting these expressions in Eq. (56) with i = i max yields an equation that can be solved for the quantity r n f at the fictitious grid point i = i max + 1. Using this quantity, we can now evaluate the usual one-sided first derivative operator in Eqs. (27)(28) at the outermost grid point i = i max . This boundary condition needs to be applied only during the predictor step of the Macormack scheme discussed in section III C. This is because the corrector step uses a backwards one-sided first derivative operator that is welldefined for the outermost grid point. Therefore, we need not worry about whether we use r orr in Eqs. (56) and (59), because the normal and computational coordinate systems coincide at the beginning of the predictor step.
The above boundary condition works well in the case of a spherically symmetric, static solution, but it does not properly handle wave propagation. This includes not only physical gravitational waves, which are absent in spherical symmetry, but also wavelike gauge modes, which are present in our numerical solution. While there exist many methods for imposing wavelike outer boundary conditions, we find for this one-dimensional problem that it suffices to simply move the outer boundary far from the origin.

VI. DIAGNOSTIC TESTS
In this section we present results of rigorous convergence tests performed with our code. In all cases below, we evolve a Schwarzschild black hole in spherical symmetry. We begin each simulation with time-symmetric initial data corresponding to the v = 0 plane of the Kruskal diagram. The initial three-metric is written in isotropic coordinates: where M is the mass of the black hole, so that the initial values of A, B, Γ r rr , Γ rT , K rr , K T , M rrr , and M rT are given by  but not explicitly given in Eqs. (73). These error terms vanish faster than (∆r) 2 as ∆r → 0, as indicated by the convergence of the plots toward each other as the grid resolution is increased.  Figure 5 shows the second-order error in the quantity r 2 β r , for the same case shown in Figures 3 and 4. This demonstrates that our solution of the minimal distortion equation is second-order convergent, and that it satisfies the outer boundary condition to O(∆r) 2 . Higher-order error terms are also present in this figure. They converge to zero as the grid resolution is increased.

B. Convergence of Constraint Equations
Even if a numerical code is second-order convergent in the self-consistent manner described above, it still may not converge to the analytic solution. To test whether this is the case, we evaluate the left-hand sides of the constraint equations (10-11), which are not explicitly solved in the code except on the initial time slice. If we are indeed solving Einstein's equations, these quantities should vanish like (∆r) 2 . Figure 6 shows left-hand side of Eq. (10b) for the same case shown in Figures 3-5, and for five different grid resolutions. For each finer grid discretization, this quantity is multiplied by a factor of 4, so that for strict secondorder convergence, the five plots should coincide. Because they do coincide except for small r where higher-order error terms are significant, we see that the momentum constraint is satisfied to O(∆r) 2 . The same is true for the other constraints.
There is an additional test we can use to determine whether we are solving the correct equations: For a static black hole, the total mass inside each radius r, which is a well-defined locally measurable quantity in spherical symmetry, should be conserved. The degree to which mass conservation is violated provides a good indicator for the overall accuracy of the code. An invariant expression for the mass inside a spherical surface of symmetry is [29] M where A is the area of the surface. In terms of our variables, this expression reduces to Figure 7 shows the quantity at time t = 5.625M as calculated from Eq. (75), for five grid discretizations. This quantity is multiplied by factors of 4 for each finer grid discretization. The mass of the system is conserved to O(∆r) 2 .

VII. DISCUSSION
While the CBY formulation of Einstein's equations at first glance looks more complicated than the usual ADM formulation, in most respects it is actually much simpler, particularly from a numerical point of view. Each equation in the CBY system is either a wave equation with an advective term, or a simple advective equation. In thex i coordinate system that is used for causal differencing, the equations are even simpler: They are either wave equations or firstorder ordinary differential equations in time. Admittedly, the right-hand sides of these equations are complicated and nonlinear, but because these right-hand sides contain no derivatives, they do not make numerical solution of the equations any more difficult. The large number of evolution equations in the CBY formalism is also not a serious drawback because the equations are all of the same form, so they can all be solved by the same method.
Unlike a previous work [4] in which we forced the horizon to sit at a particular grid point, here we allow the horizon to lie at some arbitrary location on the grid. This is a closer approximation to what we expect in the threedimensional case, where one would most likely use Cartesian coordinates to describe a spherical hole. Likewise, in the same work we imposed explicit boundary conditions on the apparent horizon in order to solve elliptic constraint equations outside the hole. Here we do not impose explicit boundary conditions (except for the shift via Eq. (48)) because in the three-dimensional case it is not only difficult to impose such conditions on a nonspherical boundary, but also the number of boundary conditions needed for solving all of the constraints in three dimensions exceeds the number of boundary conditions available at the horizon.
A key milestone in three-dimensional black hole simulations is the ability to stably move a hole through the numerical grid. This is arguably a necessary precursor to simulating binary orbits. Techniques for moving holes through the grid could be tested in spherical symmetry using our code, and we plan such tests in the future. To preserve spherical symmetry, the motion in this case would be the expansion or contraction (in coordinate space) of the hole rather than the translation of the hole's center, but since the location of the center is irrelevant for AHBC methods, many of the same methods should apply to both cases.
The code as described here tends to lose accuracy at times greater than 10 or 20 M , where M is the mass of the black hole. Because this behavior is relatively independent of numerical details such as grid resolution, we believe that it may be due to either gauge modes or unphysical rapidly-growing solutions of the evolution equations that do not satisfy the constraints. One way of solving this problem is to enforce several of the constraint equations after each time step. With this modification, we can integrate past t = 1000M . However, the details of constraint-violating solutions, gauge modes, and methods for dealing with them are beyond the scope of this paper and will be dealt with elsewhere [17].