{'fraction_non_alphanumeric': 0.027556850005600986, 'fraction_numerical': 0.007766700272581308, 'mean_word_length': 4.855268911237429, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
3+1 General Relativity in Hyperbolic Form

This paper focuses on the imposition of boundary conditions for numerical relativity simulations of black holes. This issue is used to motivate the discussion of a new hyperbolic formulation of 3+1 general relativity. The paper will appear in the Proceedings of the Les Houches School on Astrophysical Sources of Gravitational Radiation, 1995, edited by J.-A. Marck and J.-P. Lasota to be published by Springer-Verlag.


INTRODUCTION
Numerical relativity represents the only currently viable method for obtaining solutions to Einstein's equations for highly dynamical and strong field sources of gravitational radiation. The most astrophysically interesting example is probably the final stages of binary inspiral and coalescence. Partly motivated by the prospect of observations with the next generation of gravitational wave detectors, a multi-institutional "Grand Challenge" effort [1] is underway in the US aimed at solving the full Einstein equations numerically for coalescing black hole binaries. In addition to the tremendous computational resource demands of this problem, it has been realized for some time that, unfortunately, the standard 3+1 formulation of Einstein's theory as a Cauchy problem (cf. [2]) is somewhat deficient because of the difficulty of imposing boundary conditions which maintain numerical accuracy (and in some cases the physical correctness) of the solution.
The problem of boundary conditions is most dramatically evident in the study of black hole spacetimes. Inside black holes spacelike slices either a) run into singularities causing termination of simulations, b) freeze their evolution necessitating the commitment of more and more computational resources to the astrophysically irrelevant black hole throat as the simulation progresses. An obvious solution to this problem is to excise the interior of the black hole from the computational domain. Since it is impossible to identify the event-horizon dynamically during the course of a simulation, a possible alternative is to use the apparent horizon (which can be located on a single timeslice) and always lies inside the true horizon (assuming cosmic censorship holds).
Numerical techniques based on the notion of causal differencing (cf. Seidel/Suen in this volume and [3]) have been proposed for dealing with apparent horizon boundary conditions. In practice, it seems clear that the success of these techniques is crucially dependent on the form of the Einstein equations used [4,5]. For spacetimes including gravitational radiation a purely hyperbolic evolution system is imperative because boundary conditions for the full set of constraint equations are not available on the apparent horizon. Furthermore, a purely hyperbolic evolution scheme with "simple" characteristics, where the only nonzero propagation speed is the speed of light, enables one to ignore entirely the spacetime inside the apparent horizon without concern that any unphysical (gauge) field quantities should be escaping the horizon.
Outer boundary conditions for simulations on spacelike slices of asymptotically flat spacetimes are another important issue for the computation of gravitational waveforms. Since it is not feasible to simulate out to spatial infinity where there is no radiation, it is important to have boundary conditions that allow radiation to pass cleanly off the mesh. If an outgoing boundary condition is applied to the wrong variable, spurious radiation is produced which can contaminate the computed gravitational waveform. Additionally, for some problems it is necessary to put the outer boundary at such a small radius from the isolated source that backscatter of radiation off curvature is significant. This source of incoming radiation then needs to be built into the outer boundary conditions. The usual approach is to match the interior numerical solution onto perturbation theory for the exterior region [6,7,8,9,10]. Work is also underway to allow the connection of interior Cauchy solutions to exterior numerical solutions on characteristic hypersurfaces [11]. Both approaches benefit greatly from the use of a hyperbolic evolution scheme with simple characteristics for the interior solution. Outer boundary values can be assigned without the necessity of forming complicated admixtures of gauge and physical data.
In this paper we motivate and discuss a remarkable new hyperbolic formulation of general relativity [12,13] which may be thought of as a natural extension of the usual 3+1 split of spacetime. This formulation preserves complete spatial covariance by means of an arbitrary shift vector. The standard 3+1 treatment [14,2], is gauge covariant in this sense but not hyperbolic. Naturally, our formulation does require a condition on the time slicing to deal with the timereparametrization invariance of the theory. This is physically intuitive; for example, we believe that a complete understanding of the generality of slicing conditions is a necessary first step towards addressing the problem of time in quantum gravity. We expect many other applications of this formulation. Early indications are that it will provide a powerful new approach to perturbation theory and approximation schemes for general relativity.
The plan of this paper is as follows. First we will motivate the derivation of wave equations for general relativity by considering the vastly simpler case of a scalar wave and show how causal boundary conditions can be implemented.
We then turn our attention to electromagnetism to demonstrate the general procedure for gauge theories for producing wave equations. We then take general relativity in 3+1 form and apply the same method to obtain an explicitly hyperbolic formulation. This formulation is then written in first-order symmetric hyperbolic form ideal for numerical implementation. Finally we discuss perturbative reductions of these equations and their use in outer boundary conditions and radiation extraction.

BOUNDARY CONDITIONS FOR THE SCALAR WAVE EQUA-TION
Consider the simple scalar wave equation in flat space: Boundary conditions for this equation are straightforward to state and implement because the equation is manifestly hyperbolic. Since the equation is linear, it is clear how to impose outgoing wave boundary conditions at the edge of a numerical domain. It is also possible to employ inner "no boundary" or causal boundary conditions at the edge of an expanding null-surface (analogous to a black hole). Not surprisingly, since the equation is purely hyperbolic, stable and accurate solutions can be obtained by merely ignoring the causally disconnected region inside the null boundary.
We have performed tests of this notion using the simple scalar wave equation (and straightforward nonlinear extensions) in 2+1 dimensions and Cartesian coordinates: An arbitrary point in the grid is used as the origin of an expanding spherical wavefront that itself is used as the inner boundary. Outgoing wave boundary conditions are imposed at the rectangular outer boundary. Equation (2) is integrated using a predictor-corrector scheme. The second spatial derivatives are computed using centered differences except at the boundaries. At every timestep, the expanding inner boundary is located, boundary points identified and special finite difference operators constructed as shown in Figure 1. Enough points are used so that the 2nd-order error term can be set equal to the 2ndorder error in the centered difference scheme used in the mesh interior. We also simulate the effect of a shift-vector or grid velocity by allowing the coordinate system to be time dependent. This is accomplished by a highorder interpolation step following each time-step update. Since this shift-vector can be larger than the speed of propagation, it is possible for grid points to "emerge" from inside the horizon. These points are filled with data from outside using high-order extrapolation. This code has been extensively tested with grid velocities up to 5 times the propagation speed. It is stable, secondorder convergent, and equal in accuracy to a comparison code using standard boundary points boundary conditions. Nonlinear source terms seem to present no problems for this scheme. A similar, equally successful, algorithm has been developed in the context of a flux-conservative first-order version of (2) and a Lax-Wendroff evolution scheme [15].

DEVELOPMENT OF HYPERBOLIC SYSTEM
In the previous section we have demonstrated that hyperbolic wave equations are very amenable to imposition of causal boundary conditions. Here we discuss the construction of analogous wave equations for considerably more complicated (and gauge dependent) theories.

Electromagnetism
As a first example of a gauge theory, consider Maxwell's equations in flat space, written here in 3+1 form [16]. The dynamical equations are where A i is the vector potential, E i is the electric field, J i is the current source, ∇ is a 3D flat space covariant derivative, and φ is the gauge variable (scalar potential). These equations are supplemented by the initial value constraint with ρ the charge density, and a gauge condition on A i which requires the computation of φ. To produce a wave equation, one approach is to take a time derivative of the A i evolution equation and substitute the E i evolution equation. To produce a D'Alembertian operator, it is necessary to apply, for example, a transversality condition on A i which in turn imposes a radiation gauge condition on φ: . We have obtained a wave equation for A i in the "Coulomb gauge." Alternatively, we could employ the Lorentz gauge: An alternative, gauge-covariant, approach is to take a time-derivative of the evolution equation for the electric field: and use the constraint (5) to eliminate the first term yielding the wave equation Interestingly, A i doesn't appear in (7); the dynamics of electromagnetism have been cleanly separated from the gauge-dependent evolution of the vector and scalar potentials.

General Relativity
Consider a globally hyperbolic manifold of topology Σ × R with metric g µν . A foliation of this spacetime is defined by a closed 1-form ω = ∇ α t where t is this coordinate time function and ω has normalization ||ω|| = −N 2 with N the lapse function. The four-dimensional line-element associated with g µν may be decomposed in the general ADM [14] form as where β i is the shift vector. It is convenient to introduce the non-coordinate co-frame, with corresponding dual (convective) derivatives In this non-coordinate basis the ADM metric takes the simple form: Note that where the C's are the structure functions of the co-frame, dθ α = − 1 2 C βγ α θ β ∧ θ γ . Instead of using the usual time-congruence ∂/∂t = ∂ 0 + β k ∂ k which follows the spatial coordinates to define our evolution direction, we define a more natural time derivative for evolution [2] where L β is the Lie derivative in a time slice Σ along the shift vector. In combination with the lapse as N −1∂ 0 , this is the derivative with respect to proper time along the normal to Σ, and it always lies inside the light cone, in contrast to ∂/∂t. In addition, it has the useful property that it commutes with the spatial coordinate derivatives, [∂ 0 , ∂ i ] = 0. This time-derivative is particularly appropriate to our form of the equations as we have subtracted out the momentum constraints which, in the Hamiltonian formulation, turn out to be generated by the shift-vector. The dynamical variables in the standard 3+1 decomposition are the 3-metric g ij and the extrinsic curvature of the slice Σ as defined by the relation∂ In four dimensions, we can write Einstein's equation as Using the moving basis defined above, (14) can be split up into constraints and evolution equations. The time-time part of the Einstein equation leads to the Hamiltonian constraint where H ≡ K i i andR is the trace of the 3-dimensional Ricci tensorR ij (barred quantities are always spatial in our notation with indices running from 1-3). The time-space parts of Einstein's equation yield the momentum constraints: The purely spatial parts of Einstein's equation give us the evolution of the extrinsic curvature: The standard 3+1 formulation of general relativity consists of the evolution equations (13) and (17) with initial data (g ij , K ij ) satisfying the constraints (15) and (16). These equations are supplemented by equations for the sources (if any) and for the kinematical variables β i and N . The "slicing" equation for N often determined by a condition on H using the trace of (17).
To derive a wave equation for general relativity, one could, for example, follow the classic method of eliminating from R µν the unwanted second derivatives of g µν by using the full spacetime harmonic condition Γ µ = 0 [17,18,19]. One would then obtain a non-geometric D'Alembertian g αβ ∂ α ∂ β g µν . This procedure is analogous to using the Lorentz gauge in Maxwell's theory.
Instead, we shall follow the spatially gauge-covariant analog of the procedure that led to (6) and (7) for Maxwell's equations. The spatial metric g ij is analogous to A i , the shift β k to φ, and the extrinsic curvature K ij of Σ to E i . The lapse, on the other hand, is a quantity found only in time-reparametrization invariant theories. We take a time derivative of the equations of motion and subtract spatial gradients of the momentum constraints, thus obtaining a new quantity Ω ij∂ In terms of the dynamical gravitational variables, Ω ij may be expressed as where ✷ = −N −1∂ 0 N −1∂ 0 +∇ k∇ k is the physical wave operator for arbitrary β k . It is constructed from second proper time-derivatives and second covariant spatial-derivatives. The source term is given by and contains no second-derivatives of the extrinsic curvature. The slicing term S ij is given by For Ω ij to produce a wave equation, S ij must be equal to a functional involving fewer than second derivatives of K ij . The most obvious way to guarantee this is to use the harmonic condition (cf. [20] when β k = 0) (This can easily be generalized by adding an ordinary well behaved function f (t, x) to the right hand side. The slicing generality compatible with hyperbolic evolution schemes is the subject of Ref. [21].) Imposing (22) for all time amounts to imposing an equation of motion for N . The shift is completely arbitrary: it can be given as a function of space and time or solved for on each timestep based on some condition on the evolved variables. This "System I" is purely hyperbolic. Initial data consists of 3-metric and extrinsic curvature satisfying the constraints and proper time-derivative of the extrinsic curvature satisfying (17). The initial value for the lapse must also be specified. An alternative to the harmonic condition is to specify the trace of the extrinsic curvature as a known function for all time H = h(x, t). This also eliminates the second derivatives of unknown functions in S ij and provides a time-dependent elliptic equation, to solve for N on each timestep. The shift vector is still freely specifiable. This "System II" is mixed hyperbolic/elliptic. It is possible to prove [12] using the doubly contracted Bianchi identity, that Systems I and II are completely equivalent to Einstein's theory. These systems are summarized in Table I. It is easily seen that for first-order perturbations of static backgrounds, the evolution equation for the 3-metric (13) and the wave equation for the extrinsic curvature (19) become completely decoupled. We will explore this idea further in the section on perturbative reduction. This situation is analogous to the separation of the equation for A i from the wave equation for E i we saw in linear electromagnetism. It is also consistent with physical intuition about the separation of transverse wave motion from longitudinal fields in general relativity. Locally, the 3-metric provides a background on top of which the extrinsic curvature propagates.

First-order hyperbolic form for System I
In order to reduce System I to hyperbolic form it is necessary to define some new variables. (Here we restrict ourselves to the vacuum case and the simple harmonic slicing condition f (t, x) = 0.) We introduce a i = N −1∇ i N -the acceleration of the local Eulerian observers (those at rest in the time slices)its derivatives a 0i = N −1∂ 0 a i and a ji =∇ j a i = a ij , as well as time and space derivatives of the extrinsic curvaturê and M kij =∇ k K ij . System I can now be cast in flux-conservative first-order symmetric hyperbolic form [12,13]. The 49 unknowns of the first-order system are g ij , N , K ij L ij , M kij , a i , a ji and a 0i , and the equations are (13), (22), (17), (24, and∂ where J ij is computed using (20). To reduce this system to completely firstorder form, the 3-dimensional Riemann curvature appearing in J ij is expressed in terms of the 3-dimensional Ricci curvature using the identitȳ The 3D Ricci tensor is then eliminated using (17) rewritten in terms of firstorder variables:R The 4D Ricci tensor is computed from sources using (14). To demonstrate strictly first-order form it is also necessary to make Christoffel symbols part of the system by introducing a background metric as in Ref. [12]. We stress again that the shift vector is not one of the unknown fields; the form of the equations is completely independent of β k . With this form of the equations it is possible to read off the characteristic speeds of the different fields and verify one's physical expectations about the propagating degrees of freedom. Since there are no spatial derivatives on the left-hand-side of their evolution equations, we see that g ij , K ij , N , and a i all propagate with zero speed with respect to the Eulerian observers: they glide up the the normal to the foliation driven by the dynamical sources. Only the derivatives of the extrinsic curvature (L ij , M ijk ) and the derivatives of the acceleration (a oi , a ij ) propagate with the speed of light. These quantities all appear in the components of the spacetime Riemann tensor and thus represent tidal fields.
Along with A. Anderson (UNC) we have performed numerical tests of this form of the equations on the simple dynamical problem of even and odd-parity cylindrical waves. A Lax-Wendroff scheme is easily coded for the 27 equations necessary to describe this system with complete spatial gauge freedom. We find that a stable and accurate evolution can be computed which is comparable with that obtained by solving the usual 3+1 equations in fully harmonic coordinates (see [22] for this version of the equations).

PERTURBATIVE REDUCTION AND OUTER BOUNDARY CONDITIONS
Here we sketch the reduction of System I for perturbations of the Schwarzschild metric. Full details are given in Ref. [10]. This reduction has helped us elucidate many aspects of the full theory and provides new insight into the nature of gauge-invariant perturbation theory. It also allows us to define a framework for both radiation extraction and outer boundary conditions based on Schwarzschild perturbation theory. Such a framework can be used in conjunction with numerical simulations.

First-order perturbation theory of Schwarzschild
For first order perturbations of static Schwarzschild, we make the following decomposition: Tildes denote background values. The background metric and lapse take their standard static Schwarzschild values and the background extrinsic curvature and shift are zero. Unless otherwise noted, covariant derivatives are with respect to the background metric. The evolution of the 3-metric (13) reduces to Notice that this equation is entirely gauge dependent: the arbitrary choice of shift v i translates into arbitrary distortion of metric perturbations. The harmonic condition for the lapse splits into: The wave equation for the extrinsic curvature (19) reduces to Here κ = K i i andR ij andR ijkl are background, spatial Ricci and Riemann tensors respectively. This equation is entirely decoupled from the evolution of the 3-metric perturbation (37). It could be directly evolved in an exterior region as a perturbative version of the full evolution equations. (Note that for flat space,Ñ = 1, (40) reduces to ✷κ ij = 0 which has radiative solutions corresponding to first time-derivatives of the usual transverse-traceless metric perturbations for gravitational waves.) Alternatively, one can perform a decomposition of (40) in terms of tensor spherical harmonics and produce scalar wave equations for the different ℓ, m mode combinations. Here, for simplicity, we restrict attention to the slicing independent odd-parity perturbations. For odd-parity perturbations, κ ij is decomposed with two tensor spherical harmonics and two amplitude functions. The component κ rφ is expressed in terms of the amplitude function a × and angular functions as κ rφ = a × (t, r) sin θ∂ θ Y ℓm .
Taking the r-φ component of (40) and utilizing the φ component of the momentum constraint, a 1-dimensional scalar wave equation purely in terms of the amplitude function a × is formed: To form the standard Regge-Wheeler equation for odd-parity perturbations of Schwarzschild (cf. [23]) one takes a time-derivative of this equation using which is the perturbative reduction of (17) for odd-parity perturbations. Here the covariant derivatives are with respect to the perturbed background and the 3D Ricci tensor is computed from the perturbed metric. The variable ∂ t a × satisfies the usual Regge-Wheeker equation. We note that no work has been required to construct spatial gauge invariants. These come "for free" in our spatially covariant wave-equation.
The situation for even parity perturbations is somewhat more involved because the lapse perturbations (39) couple into our wave equation via the trace of κ ij and the harmonic slicing condition at the extrinsic curvature level. The same basic procedure holds, however. Using a tensor spherical-harmonic decomposition and the radial component of the momentum constraint, coupled 1D wave equations are formed for projections of κ rr and of κ. Connection to the Zerilli equation can be made by taking a time-derivative of these equations. The usual gauge invariant perturbation equations for Schwarzschild spacetime [23] are seen, not surprisingly, to represent curvature evolution.

Radiation extraction and outer boundary conditions
Perturbation theory has proven to be a powerful tool for extracting physical information from numerically generated spacetimes (cf. [6,7,8,9]). The basic scheme is to match the full nonlinear interior solution to perturbation theory along the timelike cylinder representing the boundary of the computational domain. This idea is illustrated in Figure 2.
The main steps in this scheme are to 1) construct perturbatively gaugeinvariant quantities from evolved code variables, 2) propagate these gaugeinvariants to large radius to remove near-zone effects, 3) use information from step 2 to construct code variables at the edge of the mesh, thus providing outer boundary conditions. The construction of gauge-invariants makes it possible to use the same extraction procedure in conjunction with numerical simulations using different choices of spatial gauges. This follows closely the conceptual picture for the calculation of gravitational radiation from isolated sources laid out by Thorne [24]. The calculation of the strong field and dynamical source is performed by the numerical simulation. The overlap/matching region is in the nondynamical near zone region (within a typical wavelength of the source). The goal is to compute waveforms in the wave zone beyond which the geometric optics approximation can be used to propagate the waves. In the procedure shown in Fig. 2, the waveform is read off the perturbative variables at the outer boundary of the exterior evolution. Effects of backscatter off background curvature between the outer boundary of the interior nonlinear solution and the outer boundary of the exterior perturbative solution have been taken into account in both the waveform and the boundary conditions imposed on the interior solution.
As should be clear from the discussion in the previous section, the new hyperbolic formulation elucidates the process of attaching the standard 3+1 variables g ij , K ij etc. onto perturbation theory. This subject is explored fully both for weak-field and Schwarzschild perturbation theory in Ref. [10]. In the weak field case, the exterior perturbative evolution can be done analytically. For Schwarzschild perturbations this requires a straightforward numerical integration using the same coordinate time steps as the interior evolution. Boundary data for the exterior equations is computed via multipolar projections of the components of K ij and L ij . The Schwarzschild mass is found from the ADM . During the course of the evolution, boundary data for the perturbative evolution is read off from the nonlinear solution at the inner boundary. Data from the perturbative solution is used in turn to construct outer boundary data for the nonlinear simulation. Approximate asymptotic waveforms are read off the perturbative solution at large radius.
surface integral performed near the edge of the interior mesh.
To produce boundary data for the interior simulation, the components of K ij and L ij are reformed from the perturbative variables using the momentum constraint equations. For System I, the lapse is determined by the harmonic slicing condition which ties it to the trace of the extrinsic curvature. So the lapse at the outer boundary is set directly from the exterior evolution. (If an elliptic slicing condition is used, then the lapse at the outer boundary will be known given boundary values for the extrinsic curvature and an imposed condition on (23).) The boundary condition on the 3-metric is determined using (13), the known boundary values for the extrinsic curvature and the lapse, and the chosen boundary condition on the shift vector components.
As mentioned in the previous section, an alternative procedure for the exterior evolution is to simply integrate Eq. (40) in the exterior on a 3D finite difference grid. This can be accomplished either using a Cauchy or characteristic formulation of the equation and a spherical polar topology numerical mesh for computational efficiency. Since the coordinate singularity at r = 0 will not be part of the evolution domain, the usual difficulties with numerical instabilities will be avoided. In addition, it will be sufficient to perform adaptive mesh refinement, if desired, in only the radial direction. Imposition of boundary values is trivial for the extrinsic curvature and proceeds exactly as above for the other variables.