Numerical Relativity Papers Published by Members of the SXS Collaboration

We present a code for solving the coupled Einstein-hydrodynamics equations to evolve relativistic, self-gravitating fluids. The Einstein field equations are solved in generalized harmonic coordinates on one grid using pseudospectral methods, while the fluids are evolved on another grid using shock-capturing finite difference or finite volume techniques. We show that the code accurately evolves equilibrium stars and accretion flows. Then we simulate an equal-mass nonspinning black hole-neutron star binary, evolving through the final four orbits of inspiral, through the merger, to the final stationary black hole. The gravitational waveform can be reliably extracted from the simulation.

We test the resummation techniques used in developing Pad\'e and Effective One Body (EOB) waveforms for gravitational wave detection. Convergence tests show that Pad\'e approximants of the gravitational wave energy flux do not accelerate the convergence of the standard Taylor approximants even in the test mass limit, and there is no reason why Pad\'e transformations should help in estimating parameters better in data analysis. Moreover, adding a pole to the flux seems unnecessary in the construction of these Pad\'e-approximated flux formulas. Pad\'e approximants may be useful in suggesting the form of fitting formulas. We compare a 15-orbit numerical waveform of the Caltech-Cornell group to the suggested Pad\'e waveforms of Damour et al. in the equal mass, nonspinning quasi-circular case. The comparison suggests that the Pad\'e waveforms do not agree better with the numerical waveform than the standard Taylor based waveforms. Based on this result, we design a simple EOB model by modifiying the ET EOB model of Buonanno et al., using the Taylor series of the flux with an unknown parameter at the fourth post-Newtonian order that we fit for. This simple EOB model generates a waveform having a phase difference of only 0.002 radians with the numerical waveform, much smaller than 0.04 radians the phase uncertainty in the numerical data itself. An EOB Hamiltonian can make use of a Pad\'e transformation in its construction, but this is the only place Pad\'e transformations seem useful.

Expressions for the gravitational wave (GW) energy flux and center-of-mass energy of a compact binary are integral building blocks of post-Newtonian (PN) waveforms. In this paper, we compute the GW energy flux and GW frequency derivative from a highly accurate numerical simulation of an equal-mass, non-spinning black hole binary. We also estimate the (derivative of the) center-of-mass energy from the simulation by assuming energy balance. We compare these quantities with the predictions of various PN approximants (adiabatic Taylor and Pade models; non-adiabatic effective-one-body (EOB) models). We find that Pade summation of the energy flux does not accelerate the convergence of the flux series; nevertheless, the Pade flux is markedly closer to the numerical result for the whole range of the simulation (about 30 GW cycles). Taylor and Pade models overestimate the increase in flux and frequency derivative close to merger, whereas EOB models reproduce more faithfully the shape of and are closer to the numerical flux, frequency derivative and derivative of energy. We also compare the GW phase of the numerical simulation with Pade and EOB models. Matching numerical and untuned 3.5 PN order waveforms, we find that the phase difference accumulated until $M \omega = 0.1$ is -0.12 radians for Pade approximants, and 0.50 (-0.28) radians for an EOB approximant with Keplerian (non-Keplerian) flux. We fit free parameters within the EOB models to minimize the phase difference, and discover degeneracies among these parameters. By tuning pseudo 4PN order coefficients in the radial potential or in the flux, or, if present, the location of the pole in the flux, we find that the accumulated phase difference can be reduced - if desired - to much less than the estimated numerical phase error (0.04 radians).

We present a new numerical scheme to solve the initial value problem for black hole-neutron star binaries. This method takes advantage of the flexibility and fast convergence of a multidomain spectral representation of the initial data to construct high-accuracy solutions at a relatively low computational cost. We provide convergence tests of the method for both isolated neutron stars and irrotational binaries. In the second case, we show that we can resolve the small inconsistencies that are part of the quasi-equilibrium formulation, and that these inconsistencies are significantly smaller than observed in previous works. The possibility of generating a wide variety of initial data is also demonstrated through two new configurations inspired by results from binary black holes. First, we show that choosing a modified Kerr-Schild conformal metric instead of a flat conformal metric allows for the construction of quasi-equilibrium binaries with a spinning black hole. Second, we construct binaries in low-eccentricity orbits, which are a better approximation to astrophysical binaries than quasi-equilibrium systems.

Numerical simulations of 15 orbits of an equal-mass binary black-hole system are presented. Gravitational waveforms from these simulations, covering more than 30 cycles and ending about 1.5 cycles before merger, are compared with those from quasicircular zero-spin post-Newtonian (PN) formulae. The cumulative phase uncertainty of these comparisons is about 0.05 radians, dominated by effects arising from the small residual spins of the black holes and the small residual orbital eccentricity in the simulations. Matching numerical results to PN waveforms early in the run yields excellent agreement (within 0.05 radians) over the first ~15 cycles, thus validating the numerical simulation and establishing a regime where PN theory is accurate. In the last 15 cycles to merger, however, generic time-domain Taylor approximants build up phase differences of several radians. But, apparently by coincidence, one specific post-Newtonian approximant, TaylorT4 at 3.5PN order, agrees much better with the numerical simulations, with accumulated phase differences of less than 0.05 radians over the 30-cycle waveform. Gravitational-wave amplitude comparisons are also done between numerical simulations and post-Newtonian, and the agreement depends on the post-Newtonian order of the amplitude expansion: the amplitude difference is about 6%–7% for zeroth order and becomes smaller for increasing order. A newly derived 3.0PN amplitude correction improves agreement significantly (<1% amplitude difference throughout most of the run, increasing to 4% near merger) over the previously known 2.5PN amplitude terms.

We investigate the late time behavior of a scalar field on a fixed Kerr background using a 2+1 dimensional pseudospectral evolution code. We compare evolutions of pure axisymmetric multipoles in both Kerr-Schild and Boyer-Lindquist coordinates. We find that the late-time power-law decay rate depends upon the slicing of the background, confirming previous theoretical predictions for those decay rates. The accuracy of the numerical evolutions is sufficient to decide unambiguously between competing claims in the literature.

The generalized harmonic representation of Einstein's equation is manifestly hyperbolic for a large class of gauge conditions. Unfortunately most of the useful gauges developed over the past several decades by the numerical relativity community are incompatible with the hyperbolicity of the equations in this form. This paper presents a new method of imposing gauge conditions that preserves hyperbolicity for a much wider class of conditions, including as special cases many of the standard ones used in numerical relativity: e.g., K-freezing, Gamma-freezing, Bona-Masso slicing, conformal Gamma-drivers, etc. Analytical and numerical results are presented which test the stability and the effectiveness of this new gauge driver evolution system.

The increasing sophistication and accuracy of numerical simulations of compact binaries (especially binary black holes) presents the opportunity to test the regime in which post-Newtonian (PN) predictions for the emitted gravitational waves are accurate. In order to confront numerical results with those of post-Newtonian theory, it is convenient to compare multipolar decompositions of the two waveforms. It is pointed out here that the individual modes can be computed to higher post-Newtonian order by examining the radiative multipole moments of the system, rather than by decomposing the 2.5PN polarization waveforms. In particular, the dominant (l = 2, m = 2) mode can be computed to 3PN order. Individual modes are computed to as high a post-Newtonian order as possible given previous post-Newtonian results.

We present a straightforward approach for estimating the final black hole spin of a binary black hole coalescence with arbitrary initial masses and spins. Making some simple assumptions, we estimate the final angular momentum to be the sum of the individual spins plus the orbital angular momentum of a test particle orbiting at the last stable orbit around a Kerr black hole with a spin parameter of the final black hole. The formula we obtain is able to reproduce with reasonable accuracy the results from available numerical simulations, but, more importantly, it can be used to investigate what configurations might give rise to interesting dynamics. In particular, we discuss scenarios which might give rise to a ``flip'' in the direction of the total angular momentum of the system. By studying the dependence of the final spin upon the mass ratio and initial spins we find that our simple approach suggests that it is not possible to spin-up a black hole to extremal values through merger scenarios irrespective of the mass ratio of the objects involved.

We analyze Einstein's vacuum field equations in generalized harmonic coordinates on a compact spatial domain with boundaries. We specify a class of boundary conditions which is constraint-preserving and sufficiently general to include recent proposals for reducing the amount of spurious reflections of gravitational radiation. In particular, our class comprises the boundary conditions recently proposed by Kreiss and Winicour, a geometric modification thereof, the freezing-Psi0 boundary condition and the hierarchy of absorbing boundary conditions introduced by Buchman and Sarbach. Using the recent technique developed by Kreiss and Winicour based on an appropriate reduction to a pseudo-differential first order system, we prove well posedness of the resulting initial-boundary value problem in the frozen coefficient approximation. In view of the theory of pseudo-differential operators it is expected that the full nonlinear problem is also well posed. Furthermore, we implement some of our boundary conditions numerically and study their effectiveness in a test problem consisting of a perturbed Schwarzschild black hole.

In this Comment we compute the contributions of the radiation reaction force in the 2.5 post-Newtonian (PN) gravitational wave polarizations for compact binaries in circular orbits. (i) We point out and correct an inconsistency in the derivation of Arun, Blanchet, Iyer, and Qusailah. (ii) We prove that all contributions from radiation reaction in the 2.5PN waveform are actually negligible since they can be absorbed into a modification of the orbital phase at the 5PN order.

Various methods of treating outer boundaries in numerical relativity are compared using a simple test problem: a Schwarzschild black hole with an outgoing gravitational wave perturbation. Numerical solutions computed using different boundary treatments are compared to a `reference' numerical solution obtained by placing the outer boundary at a very large radius. For each boundary treatment, the full solutions including constraint violations and extracted gravitational waves are compared to those of the reference solution, thereby assessing the reflections caused by the artificial boundary. These tests use a first-order generalized harmonic formulation of the Einstein equations. Constraint-preserving boundary conditions for this system are reviewed, and an improved boundary condition on the gauge degrees of freedom is presented. Alternate boundary conditions evaluated here include freezing the incoming characteristic fields, Sommerfeld boundary conditions, and the constraint-preserving boundary conditions of Kreiss and Winicour. Rather different approaches to boundary treatments, such as sponge layers and spatial compactification, are also tested. Overall the best treatment found here combines boundary conditions that preserve the constraints, freeze the Newman-Penrose scalar Psi_0, and control gauge reflections.

A new constraint suppressing formulation of the Einstein evolution equations is presented, generalizing the five-parameter first-order system due to Kidder, Scheel and Teukolsky (KST). The auxiliary fields, introduced to make the KST system first-order, are given modified evolution equations designed to drive constraint violations toward zero. The algebraic structure of the new system is investigated, showing that the modifications preserve the hyperbolicity of the fundamental and constraint evolution equations. The evolution of the constraints for pertubations of flat spacetime is completely analyzed, and all finite-wavelength constraint modes are shown to decay exponentially when certain adjustable parameters satisfy appropriate inequalities. Numerical simulations of a single Schwarzschild black hole are presented, demonstrating the effectiveness of the new constraint-damping modifications.

We present a new code for solving the coupled Einstein-hydrodynamics equations to evolve relativistic, self-gravitating fluids. The Einstein field equations are solved on one grid using pseudospectral methods, while the fluids are evolved on another grid by finite differencing. We discuss implementation details, such as the communication between the grids and the treatment of stellar surfaces, and present code tests.

Binary black hole simulations starting from quasi-circular (i.e., zero radial velocity) initial data have orbits with small but non-zero orbital eccentricities. In this paper the quasi-equilibrium initial-data method is extended to allow non-zero radial velocities to be specified in binary black hole initial data. New low-eccentricity initial data are obtained by adjusting the orbital frequency and radial velocities to minimize the orbital eccentricity, and the resulting ($\sim 5$ orbit) evolutions are compared with those of quasi-circular initial data. Evolutions of the quasi-circular data clearly show eccentric orbits, with eccentricity that decays over time. The precise decay rate depends on the definition of eccentricity; if defined in terms of variations in the orbital frequency, the decay rate agrees well with the prediction of Peters (1964). The gravitational waveforms, which contain $\sim 8$ cycles in the dominant l=m=2 mode, are largely unaffected by the eccentricity of the quasi-circular initial data. The overlap between the dominant mode in the quasi-circular evolution and the same mode in the low-eccentricity evolution is about 0.99.

We study the appearance of multiple solutions to certain decompositions of Einstein's constraint equations. Pfeiffer and York recently reported the existence of two branches of solutions for identical background data in the extended conformal thin-sandwich decomposition. We show that the Hamiltonian constraint alone, when expressed in a certain way, admits two branches of solutions with properties very similar to those found by Pfeiffer and York. We construct these two branches analytically for a constant-density star in spherical symmetry, but argue that this behavior is more general. In the case of the Hamiltonian constraint this non-uniqueness is well known to be related to the sign of one particular term, and we argue that the extended conformal thin-sandwich equations contain a similar term that causes the breakdown of uniqueness.

We study the numerical implementation of a set of boundary conditions derived from the isolated horizon formalism, and which characterize a black hole whose horizon is in quasi-equilibrium. More precisely, we enforce these geometrical prescriptions as inner boundary conditions on an excised sphere, in the numerical resolution of the Conformal Thin Sandwich equations. As main results, we firstly establish the consistency of including in the set of boundary conditions a "constant surface gravity" prescription, interpretable as a lapse boundary condition, and secondly we assess how the prescriptions presented recently by Dain et al. for guaranteeing the well-posedness of the Conformal Transverse Traceless equations with quasi-equilibrium horizon conditions extend to the Conformal Thin Sandwich elliptic system. As a consequence of the latter analysis, we discuss the freedom of prescribing the expansion associated with the ingoing null normal at the horizon.

The accuracy and stability of the Caltech-Cornell pseudospectral code is evaluated using the KST representation of the Einstein evolution equations. The basic ``Mexico City Tests'' widely adopted by the numerical relativity community are adapted here for codes based on spectral methods. Exponential convergence of the spectral code is established, apparently limited only by numerical roundoff error. A general expression for the growth of errors due to finite machine precision is derived, and it is shown that this limit is achieved here for the linear plane-wave test. All of these tests are found to be stable, except for simulations of high amplitude gauge waves with nontrivial shift.

We construct exact solutions to the Bianchi equations on a flat spacetime background. When the constraints are satisfied, these solutions represent in- and outgoing linearized gravitational radiation. We then consider the Bianchi equations on a subset of flat spacetime of the form [0,T] x B_R, where B_R is a ball of radius R, and analyze different kinds of boundary conditions on \partial B_R. Our main results are: i) We give an explicit analytic example showing that boundary conditions obtained from freezing the incoming characteristic fields to their initial values are not compatible with the constraints. ii) With the help of the exact solutions constructed, we determine the amount of artificial reflection of gravitational radiation from constraint-preserving boundary conditions which freeze the Weyl scalar Psi_0 to its initial value. For monochromatic radiation with wave number k and arbitrary angular momentum number l >= 2, the amount of reflection decays as 1/(kR)^4 for large kR. iii) For each L >= 2, we construct new local constraint-preserving boundary conditions which perfectly absorb linearized radiation with l <= L. (iv) We generalize our analysis to a weakly curved background of mass M, and compute first order corrections in M/R to the reflection coefficients for quadrupolar odd-parity radiation. For our new boundary condition with L=2, the reflection coefficient is smaller than the one for the freezing Psi_0 boundary condition by a factor of M/R for kR > 1.04. Implications of these results for numerical simulations of binary black holes on finite domains are discussed.

A method is introduced for solving Einstein's equations using two distinct coordinate systems. The coordinate basis vectors associated with one system are used to project out components of the metric and other fields, in analogy with the way fields are projected onto an orthonormal tetrad basis. These field components are then determined as functions of a second independent coordinate system. The transformation to the second coordinate system can be thought of as a mapping from the original ``inertial'' coordinate system to the computational domain. This dual-coordinate method is used to perform stable numerical evolutions of a black-hole spacetime using the generalized harmonic form of Einstein's equations in coordinates that rotate with respect to the inertial frame at infinity; such evolutions are found to be generically unstable using a single rotating coordinate frame. The dual-coordinate method is also used here to evolve binary black-hole spacetimes for several orbits. The great flexibility of this method allows comoving coordinates to be adjusted with a feedback control system that keeps the excision boundaries of the holes within their respective apparent horizons.

Pointlike objects cause many of the divergences that afflict physical theories. For instance, the gravitational binding energy of a point particle in Newtonian mechanics is infinite. In general relativity, the analog of a point particle is a black hole and the notion of binding energy must be replaced by quasilocal energy. The quasilocal energy (QLE) derived by York, and elaborated by Brown and York, is finite outside the horizon but it was not considered how to evaluate it inside the horizon. We present a prescription for finding the QLE inside a horizon, and show that it is finite at the singularity for a variety of types of black hole. The energy is typically concentrated just inside the horizon, not at the central singularity.

We present a modification to the Berger and Oliger adaptive mesh refinement algorithm designed to solve systems of coupled, non-linear, hyperbolic and elliptic partial differential equations. Such systems typically arise during constrained evolution of the field equations of general relativity. The novel aspect of this algorithm is a technique of "extrapolation and delayed solution" used to deal with the non-local nature of the solution of the elliptic equations, driven by dynamical sources, within the usual Berger and Oliger time-stepping framework. We show empirical results demonstrating the effectiveness of this technique in axisymmetric gravitational collapse simulations. We also describe several other details of the code, including truncation error estimation using a self-shadow hierarchy, and the refinement-boundary interpolation operators that are used to help suppress spurious high-frequency solution components ("noise").

We construct approximate analytical solutions to the constraint equations of general relativity for binary black holes of arbitrary mass ratio in quasicircular orbit. We adopt the puncture method of Brandt and Bruegmann to solve the constraint equations in the transverse-traceless decomposition and consider perturbations of Schwarzschild black holes caused by boosts and the presence of a binary companion. A superposition of these two perturbations then yields approximate, but fully analytic binary black hole initial data that are accurate to first order in the inverse of the binary separation and the square of the black holes' momenta. Even close to the innermost stable circular orbit, the perturbative treatment introduces errors that are remarkably small and only somewhat larger than the errors caused by the underlying assumptions of puncture data.

The construction of initial data for black-hole binaries usually involves the choice of free parameters that define the spins of the black holes and essentially the eccentricity of the orbit. Such parameters must be chosen carefully to yield initial data with the desired physical properties. In this paper, we examine these choices in detail for the quasiequilibrium method coupled to apparent-horizon/quasiequilibrium boundary conditions. First, we compare two independent criteria for choosing the orbital frequency, the "Komar-mass condition" and the "effective-potential method," and find excellent agreement. Second, we implement quasi-local measures of the spin of the individual holes, calibrate these with corotating binaries, and revisit the construction of non-spinning black hole binaries. Higher-order effects, beyond those considered in earlier work, turn out to be important. Without those, supposedly non-spinning black holes have appreciable quasi-local spin; furthermore, the Komar-mass condition and effective potential method agree only when these higher-order effects are taken into account. We compute a new sequence of quasi-circular orbits for non-spinning black-hole binaries, and determine the innermost stable circular orbit of this sequence.

A new representation of the Einstein evolution equations is presented that is first order, linearly degenerate, and symmetric hyperbolic. This new system uses the generalized harmonic method to specify the coordinates, and exponentially suppresses all small short-wavelength constraint violations. Physical and constraint-preserving boundary conditions are derived for this system, and numerical tests that demonstrate the effectiveness of the constraint suppression properties and the constraint-preserving boundary conditions are presented.

We discuss the initial-boundary value problem of General Relativity. Previous considerations for a toy model problem in electrodynamics motivate the introduction of a variational principle for the lapse with several attractive properties. In particular, it is argued that the resulting elliptic gauge condition for the lapse together with a suitable condition for the shift and constraint-preserving boundary conditions controlling the Weyl scalar Psi_0 are expected to yield a well posed initial-boundary value problem for metric formulations of Einstein's field equations which are commonly used in numerical relativity. To present a simple and explicit example we consider the 3+1 decomposition introduced by York of the field equations on a cubic domain with two periodic directions and prove in the weak field limit that our gauge condition for the lapse and our boundary conditions lead to a well posed problem. The method discussed here is quite general and should also yield well posed problems for different ways of writing the evolution equations, including first order symmetric hyperbolic or mixed first-order second-order formulations. Well posed initial-boundary value formulations for the linearization about arbitrary stationary configurations will be presented elsewhere.

This paper is concerned with the initial-boundary value problem for the Einstein equations in a first-order generalized harmonic formulation. We impose boundary conditions that preserve the constraints and control the incoming gravitational radiation by prescribing data for the incoming fields of the Weyl tensor. High-frequency perturbations about any given spacetime (including a shift vector with subluminal normal component) are analyzed using the Fourier-Laplace technique. We show that the system is boundary-stable. In addition, we develop a criterion that can be used to detect weak instabilities with polynomial time dependence, and we show that our system does not suffer from such instabilities. A numerical robust stability test supports our claim that the initial-boundary value problem is most likely to be well-posed even if nonzero initial and source data are included.

Based on the recent understanding of the role of the densitized lapse function in Einstein's equations and of the proper way to pose the thin sandwich problem, a slight readjustment of the minimal distortion shift gauge in the 3+1 approach to the dynamics of general relativity allows this shift vector to serve as the vector potential for the longitudinal part of the extrinsic curvature tensor in the new approach to the initial value problem, thus extending the initial value decomposition of gravitational variables to play a role in the evolution as well. The new shift vector globally minimizes the changes in the conformal 3-metric with respect to the spacetime measure rather than the spatial measure on the time coordinate hypersurfaces, as the old shift vector did.

Combining deeper insight of Einstein's equations with sophisticated numerical techniques promises the ability to construct accurate numerical implementations of these equations. We illustrate this in two examples, the numerical evolution of ``bubble'' and single black hole spacetimes. The former is chosen to demonstrate how accurate numerical solutions can answer open questions and even reveal unexpected phenomena. The latter illustrates some of the difficulties encountered in three-dimensional black hole simulations, and presents some possible remedies.

The need to smoothly cover a computational domain of interest generically requires the adoption of several grids. To solve the problem of interest under this grid-structure one must ensure the suitable transfer of information among the different grids involved. In this work we discuss a technique that allows one to construct finite difference schemes of arbitrary high order which are guaranteed to satisfy linear numerical and strict stability. The technique relies on the use of difference operators satisfying summation by parts and {\it penalty techniques} to transfer information between the grids. This allows the derivation of semidiscrete energy estimates for problems admitting such estimates at the continuum. We analyze several aspects of this technique when used in conjuction with high order schemes and illustrate its use in one, two and three dimensional numerical relativity model problems with non-trivial topologies, including truly spherical black hole excision.

A first order symmetric hyperbolic tetrad formulation of the Einstein equations developed by Estabrook and Wahlquist and put into a form suitable for numerical relativity by Buchman and Bardeen (the WEBB formulation) is adapted to explicit spherical symmetry and tested for accuracy and stability in the evolution of spherically symmetric black holes (the Schwarzschild geometry). The lapse and shift which specify the evolution of the coordinates relative to the tetrad congruence are reset at frequent time intervals to keep the constant-time hypersurfaces nearly orthogonal to the tetrad congruence and the spatial coordinate satisfying a kind of minimal rate of strain condition. By arranging through initial conditions that the constant-time hypersurfaces are asymptotically hyperbolic, we simplify the boundary value problem and improve stability of the evolution. Results are obtained for both tetrad gauges (``Nester'' and ``Lorentz'') of the WEBB formalism using finite difference numerical methods. We are able to obtain stable unconstrained evolution with the Nester gauge for certain initial conditions, but not with the Lorentz gauge.

Outer boundary conditions for strongly and symmetric hyperbolic formulations of 3D Einstein's field equations with a live gauge condition are discussed. The boundary conditions have the property that they ensure constraint propagation and control in a sense made precise in this article the physical degrees of freedom at the boundary. We use Fourier-Laplace transformation techniques to find necessary conditions for the well posedness of the resulting initial-boundary value problem and integrate the resulting three-dimensional nonlinear equations using a finite-differencing code. We obtain a set of constraint-preserving boundary conditions which pass a robust numerical stability test. We explicitly compare these new boundary conditions to standard, maximally dissipative ones through Brill wave evolutions. Our numerical results explicitly show that in the latter case the constraint variables, describing the violation of the constraints, do not converge to zero when resolution is increased while for the new boundary conditions, the constraint variables do decrease as resolution is increased. As an application, we inject pulses of ``gravitational radiation'' through the boundaries of an initially flat spacetime domain, with enough amplitude to generate strong fields and induce large curvature scalars, showing that our boundary conditions are robust enough to handle nonlinear dynamics. We expect our boundary conditions to be useful for improving the accuracy and stability of current binary black hole and binary neutron star simulations, for a successful implementation of characteristic or perturbative matching techniques, and other applications. We also discuss limitations of our approach and possible future directions.

In this article, dedicated to one of the best specialist of the FOSH systems, we couple the Bianchi equations with the equations satisfied by the dynamical acceleration of a charged fluid and the derivatives of the associated Maxwell field.

The most promising way to compute the gravitational waves emitted by binary black holes (BBHs) in their last dozen orbits, where post-Newtonian techniques fail, is a quasistationary approximation introduced by Detweiler and being pursued by Price and others. In this approximation the outgoing gravitational waves at infinity and downgoing gravitational waves at the holes' horizons are replaced by standing waves so as to guarantee that the spacetime has a helical Killing vector field. Because the horizon generators will not, in general, be tidally locked to the holes' orbital motion, the standing waves will destroy the horizons, converting the black holes into naked singularities that resemble black holes down to near the horizon radius. This paper uses a spherically symmetric, scalar-field model problem to explore in detail the following BBH issues: (i) The destruction of a horizon by the standing waves. (ii) The accuracy with which the resulting naked singularity resembles a black hole. (iii) The conversion of the standing-wave spacetime (with a destroyed horizon) into a spacetime with downgoing waves by the addition of a ``radiation-reaction field''. (iv) The accuracy with which the resulting downgoing waves agree with the downgoing waves of a true black-hole spacetime (with horizon). The model problem used to study these issues consists of a Schwarzschild black hole endowed with spherical standing waves of a scalar field. It is found that the spacetime metric of the singular, standing-wave spacetime, and its radiation-reaction-field-constructed downgoing waves are quite close to those for a Schwarzschild black hole with downgoing waves -- sufficiently close to make the BBH quasistationary approximation look promising for non-tidally-locked black holes.

We describe early success in the evolution of binary black hole spacetimes with a numerical code based on a generalization of harmonic coordinates. Indications are that with sufficient resolution this scheme is capable of evolving binary systems for enough time to extract information about the orbit, merger and gravitational waves emitted during the event. As an example we show results from the evolution of a binary composed of two equal mass, non-spinning black holes, through a single plunge-orbit, merger and ring down. The resultant black hole is estimated to be a Kerr black hole with angular momentum parameter a~0.70. At present, lack of resolution far from the binary prevents an accurate estimate of the energy emitted, though a rough calculation suggests on the order of 5% of the initial rest mass of the system is radiated as gravitational waves during the final orbit and ringdown.

This thesis is concerned with formulations of the Einstein equations in axisymmetric spacetimes which are suitable for numerical evolutions. We develop two evolution systems based on the (2+1)+1 formalism. The first is a (partially) constrained scheme with elliptic gauge conditions arising from maximal slicing and conformal flatness. The second is a strongly hyperbolic first-order formulation obtained by combining the (2+1)+1 formalism with the Z4 formalism. A careful study of the behaviour of regular axisymmetric tensor fields enables us to cast the equations in a form that is well-behaved on the axis. Further topics include (non)uniqueness of solutions to the elliptic equations arising in constrained schemes, and comparisons between various boundary conditions used in numerical relativity. The numerical implementation is applied to adaptive evolutions of nonlinear Brill waves, including twist.

The conformal thin sandwich (CTS) equations are a set of four of the Einstein equations, which generalize the Laplace-Poisson equation of Newton's theory. We examine numerically solutions of the CTS equations describing perturbed Minkowski space, and find only one solution. However, we find {\em two} distinct solutions, one even containing a black hole, when the lapse is determined by a fifth elliptic equation through specification of the mean curvature. While the relationship of the two systems and their solutions is a fundamental property of general relativity, this fairly simple example of an elliptic system with non-unique solutions is also of broader interest.

The conformal method for constructing initial data for Einstein's equations is presented in both the Hamiltonian and Lagrangian picture (extrinsic curvature decomposition and conformal thin sandwich formalism, respectively), and advantages due to the recent introduction of a weight-function in the extrinsic curvature decomposition are discussed. I then describe recent progress in numerical techniques to solve the resulting elliptic equations, and explore innovative approaches toward the construction of astrophysically realistic initial data for binary black hole simulations.

In many numerical implementations of the Cauchy formulation of Einstein's field equations one encounters artificial boundaries which raises the issue of specifying boundary conditions. Such conditions have to be chosen carefully. In particular, they should be compatible with the constraints, yield a well posed initial-boundary value formulation and incorporate some physically desirable properties like, for instance, minimizing reflections of gravitational radiation. Motivated by the problem in General Relativity, we analyze a model problem, consisting of a formulation of Maxwell's equations on a spatially compact region of spacetime with timelike boundaries. The form in which the equations are written is such that their structure is very similar to the Einstein-Christoffel symmetric hyperbolic formulations of Einstein's field equations. For this model problem, we specify a family of Sommerfeld-type constraint-preserving boundary conditions and show that the resulting initial-boundary value formulations are well posed. We expect that these results can be generalized to the Einstein-Christoffel formulations of General Relativity, at least in the case of linearizations about a stationary background.

The periodic standing wave (PSW) method for the binary inspiral of black holes and neutron stars computes exact numerical solutions for periodic standing wave spacetimes and then extracts approximate solutions of the physical problem, with outgoing waves. The method requires solution of a boundary value problem with a mixed (hyperbolic and elliptic) character. We present here a new numerical method for such problems, based on three innovations: (i) a coordinate system adapted to the geometry of the problem, (ii) an expansion in multipole moments of these coordinates and a filtering out of higher moments, and (iii) the replacement of the continuum multipole moments with their analogs for a discrete grid. We illustrate the efficiency and accuracy of this method with nonlinear scalar model problems. Finally, we take advantage of the ability of this method to handle highly nonlinear models to demonstrate that the outgoing approximations extracted from the standing wave solutions are highly accurate even in the presence of strong nonlinearities.

This paper is concerned exclusively with axisymmetric spacetimes. We want to develop reductions of Einstein's equations which are suitable for numerical evolutions. We first make a Kaluza-Klein type dimensional reduction followed by an ADM reduction on the Lorentzian 3-space, the (2+1)+1 formalism. We include also the Z4 extension of Einstein's equations adapted to this formalism. Our gauge choice is based on a generalized harmonic gauge condition. We consider vacuum and perfect fluid sources. We use these ingredients to construct a strongly hyperbolic first-order evolution system and exhibit its characteristic structure. This enables us to construct constraint-preserving stable outer boundary conditions. We use cylindrical polar coordinates and so we provide a careful discussion of the coordinate singularity on axis. By choosing our dependent variables appropriately we are able to produce an evolution system in which each and every term is manifestly regular on axis.

New boundary conditions are constructed and tested numerically for a general first-order form of the Einstein evolution system. These conditions prevent constraint violations from entering the computational domain through timelike boundaries, allow the simulation of isolated systems by preventing physical gravitational waves from entering the computational domain, and are designed to be compatible with the fixed-gauge evolutions used here. These new boundary conditions are shown to be effective in limiting the growth of constraints in 3D nonlinear numerical evolutions of dynamical black-hole spacetimes.

We use the numerical solution describing the evolution of a perturbed black string presented in Choptuik et al. (2003) to elucidate the intrinsic behavior of the horizon. It is found that by the end of the simulation, the affine parameter on the horizon has become very large and the expansion and shear of the horizon in turn very small. This suggests the possibility that the horizon might pinch off in infinite affine parameter.

In calculations of the inspiral of binary black holes an intermediate approximation is needed that can bridge the post-Newtonian methods of the early inspiral and the numerical relativity computations of the final plunge. We describe here the periodic standing wave approximation: A numerical solution is found to the problem of a periodic rotating binary with helically symmetric standing wave fields, and from this solution an approximation is extracted for the physically relevant problem of inspiral with outgoing waves. The approximation underlying this approach has been recently confirmed with innovative numerical methods applied to nonlinear model problems.

A new numerical scheme to solve the Einstein field equations based upon the generalized harmonic decomposition of the Ricci tensor is introduced. The source functions driving the wave equations that define generalized harmonic coordinates are treated as independent functions, and encode the coordinate freedom of solutions. Techniques are discussed to impose particular gauge conditions through a specification of the source functions. A 3D, free evolution, finite difference code implementing this system of equations with a scalar field matter source is described. The second-order-in-space-and-time partial differential equations are discretized directly without the use first order auxiliary terms, limiting the number of independent functions to fifteen--ten metric quantities, four source functions and the scalar field. This also limits the number of constraint equations, which can only be enforced to within truncation error in a numerical free evolution, to four. The coordinate system is compactified to spatial infinity in order to impose physically motivated, constraint-preserving outer boundary conditions. A variant of the Cartoon method for efficiently simulating axisymmetric spacetimes with a Cartesian code is described that does not use interpolation, and is easier to incorporate into existing adaptive mesh refinement packages. Preliminary test simulations of vacuum black hole evolution and black hole formation via scalar field collapse are described, suggesting that this method may be useful for studying many spacetimes of interest.

A method is presented to construct initial data for Einstein's equations as a superposition of a gravitational wave perturbation on an arbitrary stationary background spacetime. The method combines the conformal thin sandwich formalism with linear gravitational waves, and allows detailed control over characteristics of the superposed gravitational wave like shape, location and propagation direction. It is furthermore fully covariant with respect to spatial coordinate changes and allows for very large amplitude of the gravitational wave.

We have modelled the nonlinear development of the secular bar-mode instability that is driven by gravitational radiation-reaction (GRR) forces in rotating neutron stars. In the absence of any competing viscous effects, an initially uniformly rotating, axisymmetric n=1/2 polytropic star with a ratio of rotational to gravitational potential energy T/|W| = 0.181 is driven by GRR forces to a bar-like structure, as predicted by linear theory. The pattern frequency of the bar slows to nearly zero, that is, the bar becomes almost stationary as viewed from an inertial frame of reference as GRR removes energy and angular momentum from the star. In this ``Dedekind-like'' state, rotational energy is stored as motion of the fluid in highly noncircular orbits inside the bar. However, in less than 10 dynamical times after its formation, the bar loses its initially coherent structure as the ordered flow inside the bar is disrupted by what appears to be a purely hydrodynamical, short-wavelength, ``shearing'' type instability. The gravitational waveforms generated by such an event are determined, and an estimate of the detectability of these waves is presented.

The complete form of the constraints following from their conformal structure is extended so as to include constant mean curvature and other mean curvature foliations. This step is demonstrated using the momentum phase space approach. This approach yields equations of exactly the same form as the extended conformal thin sandwich approach. In solving the equations, it is never necessary actually to perform a tensor decomposition.

We define and extensively test a set of boundary conditions that can be applied at black hole excision surfaces when the Hamiltonian and momentum constraints of general relativity are solved within the conformal thin-sandwich formalism. These boundary conditions have been designed to result in black holes that are in quasiequilibrium and are completely general in the sense that they can be applied with any conformal three-geometry and slicing condition. Furthermore, we show that they retain precisely the freedom to specify an arbitrary spin on each black hole. Interestingly, we have been unable to find a boundary condition on the lapse that can be derived from a quasiequilibrium condition. Rather, we find evidence that the lapse boundary condition is part of the initial temporal gauge choice. To test these boundary conditions, we have extensively explored the case of a single black hole and the case of a binary system of equal-mass black holes, including the computation of quasi-circular orbits and the determination of the inner-most stable circular orbit. Our tests show that the boundary conditions work well.

Techniques are developed for projecting the solutions of symmetric hyperbolic evolution systems onto the constraint submanifold (the constraint-satisfying subset of the dynamical field space). These optimal projections map a field configuration to the ``nearest'' configuration in the constraint submanifold, where distances between configurations are measured with the natural metric on the space of dynamical fields. The construction and use of these projections is illustrated for a new representation of the scalar field equation that exhibits both bulk and boundary generated constraint violations. Numerical simulations on a black-hole background show that bulk constraint violations cannot be controlled by constraint-preserving boundary conditions alone, but are effectively controlled by constraint projection. Simulations also show that constraint violations entering through boundaries cannot be controlled by constraint projection alone, but are controlled by constraint-preserving boundary conditions. Numerical solutions to the pathological scalar field system are shown to converge to solutions of a standard representation of the scalar field equation when constraint projection and constraint-preserving boundary conditions are used together.

We report a new critical solution found at the threshold of axisymmetric gravitational collapse of a complex scalar field with angular momentum. To carry angular momentum the scalar field cannot be axisymmetric; however, its azimuthal dependence is defined so that the resulting stress energy tensor and spacetime metric are axisymmetric. The critical solution found is non-spherical, discretely self-similar with an echoing exponent of 0.42 (+- 4%), and exhibits a scaling exponent of 0.11 (+- 10%) in near critical collapse. Our simulations suggest that the solution is universal (within the imposed symmetry class), modulo a family-dependent constant phase in the complex plane.

The use of adaptive mesh refinement (AMR) techniques is crucial for accurate and efficient simulation of higher dimensional spacetimes. In this work we develop an adaptive algorithm tailored to the integration of finite difference discretizations of wave-like equations using characteristic coordinates. We demonstrate the algorithm by constructing a code implementing the Einstein-Klein-Gordon system of equations in spherical symmetry. We discuss how the algorithm can trivially be generalized to higher dimensional systems, and suggest a method that can be used to parallelize a characteristic code.

Motivated by the need to control the exponential growth of constraint violations in numerical solutions of the Einstein evolution equations, two methods are studied here for controlling this growth in general hyperbolic evolution systems. The first method adjusts the evolution equations dynamically, by adding multiples of the constraints, in a way designed to minimize this growth. The second method imposes special constraint preserving boundary conditions on the incoming components of the dynamical fields. The efficacy of these methods is tested by using them to control the growth of constraints in fully dynamical 3D numerical solutions of a particular representation of the Maxwell equations that is subject to constraint violations. The constraint preserving boundary conditions are found to be much more effective than active constraint control in the case of this Maxwell system.

The initial value problem is introduced after a thorough review of the essential geometry. The initial value equations are put into elliptic form using both conformal transformations and a treatment of the extrinsic curvature introduced recently. This use of the metric and the extrinsic curvature is manifestly equivalent to the author's conformal thin sandwich formulation. Therefore, the reformulation of the constraints as an elliptic system by use of conformal techniques is complete.

We investigate the behavior of a dynamical scalar field on a fixed Kerr background in Kerr-Schild coordinates using a 3+1 dimensional spectral evolution code, and we measure the power-law tail decay that occurs at late times. We compare evolutions of initial data proportional to f(r) Y_lm(theta,phi) where Y_lm is a spherical harmonic and (r,theta,phi) are Kerr-Schild coordinates, to that of initial data proportional to f(r_BL) Y_lm(theta_BL,phi), where (r_BL,theta_BL) are Boyer-Lindquist coordinates. We find that although these two cases are initially almost identical, the evolution can be quite different at intermediate times; however, at late times the power-law decay rates are equal.

We present results from a numerical study of critical gravitational collapse of axisymmetric distributions of massless scalar field energy. We find threshold behavior that can be described by the spherically symmetric critical solution with axisymmetric perturbations. However, we see indications of a growing, non-spherical mode about the spherically symmetric critical solution. The effect of this instability is that the small asymmetry present in what would otherwise be a spherically symmetric self-similar solution grows. This growth continues until a bifurcation occurs and two distinct regions form on the axis, each resembling the spherically symmetric self-similar solution. The existence of a non-spherical unstable mode is in conflict with previous perturbative results, and we therefore discuss whether such a mode exists in the continuum limit, or whether we are instead seeing a marginally stable mode that is rendered unstable by numerical approximation.

Black strings, one class of higher dimensional analogues of black holes, were shown to be unstable to long wavelength perturbations by Gregory and Laflamme in 1992, via a linear analysis. We revisit the problem through numerical solution of the full equations of motion, and focus on trying to determine the end-state of a perturbed, unstable black string. Our preliminary results show that such a spacetime tends towards a solution resembling a sequence of spherical black holes connected by thin black strings, at least at intermediate times. However, our code fails then, primarily due to large gradients that develop in metric functions, as the coordinate system we use is not well adapted to the nature of the unfolding solution. We are thus unable to determine how close the solution we see is to the final end-state, though we do observe rich dynamical behavior of the system in the intermediate stages.

We discuss the initial value problem of general relativity in its recently unified Lagrangian and Hamiltonian pictures and present a multi-domain pseudo-spectral collocation method to solve the resulting coupled nonlinear partial differential equations. Using this code, we explore several approaches to construct initial data sets containing one or two black holes: We compute quasi-circular orbits for spinning equal mass black holes and unequal mass (nonspinning) black holes using the effective potential method with Bowen-York extrinsic curvature. We compare initial data sets resulting from different decompositions, and from different choices of the conformal metric with each other. Furthermore, we use the quasi-equilibrium method to construct initial data for single black holes and for binary black holes in quasi-circular orbits. We investigate these binary black hole data sets and examine the limits of large mass-ratio and wide separation. Finally, we propose a new method for constructing spacetimes with superposed gravitational waves of possibly very large amplitude.

The Einstein evolution equations have been written in a number of symmetric hyperbolic forms when the gauge fields--the densitized lapse and the shift--are taken to be fixed functions of the coordinates. Extended systems of evolution equations are constructed here by adding the gauge degrees of freedom to the set of dynamical fields, thus forming symmetric hyperbolic systems for the combined evolution of the gravitational and the gauge fields. The associated characteristic speeds can be made causal (i.e. less than or equal to the speed of light) by adjusting 14 free parameters in these new systems. And 21 additional free parameters are available, for example to optimize the stability of numerical evolutions. The gauge evolution equations in these systems are generalizations of the ``K-driver'' and ``Gamma-driver'' conditions that have been used with some success in numerical black hole evolutions.

We present a new solver for coupled nonlinear elliptic partial differential equations (PDEs). The solver is based on pseudo-spectral collocation with domain decomposition and can handle one- to three-dimensional problems. It has three distinct features. First, the combined problem of solving the PDE, satisfying the boundary conditions, and matching between different subdomains is cast into one set of equations readily accessible to standard linear and nonlinear solvers. Second, touching as well as overlapping subdomains are supported; both rectangular blocks with Chebyshev basis functions as well as spherical shells with an expansion in spherical harmonics are implemented. Third, the code is very flexible: The domain decomposition as well as the distribution of collocation points in each domain can be chosen at run time, and the solver is easily adaptable to new PDEs. The code has been used to solve the equations of the initial value problem of general relativity and should be useful in many other problems. We compare the new method to finite difference codes and find it superior in both runtime and accuracy, at least for the smooth problems considered here.

We present a new numerical code designed to solve the Einstein field equations for axisymmetric spacetimes. The long term goal of this project is to construct a code that will be capable of studying many problems of interest in axisymmetry, including gravitational collapse, critical phenomena, investigations of cosmic censorship, and head-on black hole collisions. Our objective here is to detail the (2+1)+1 formalism we use to arrive at the corresponding system of equations and the numerical methods we use to solve them. We are able to obtain stable evolution, despite the singular nature of the coordinate system on the axis, by enforcing appropriate regularity conditions on all variables and by adding numerical dissipation to hyperbolic equations.

The Einstein initial-value equations in the extrinsic curvature (Hamiltonian) representation and conformal thin sandwich (Lagrangian) representation are brought into complete conformity by the use of a decomposition of symmetric tensors which involves a weight function. In stationary spacetimes, there is a natural choice of the weight function such that the transverse traceless part of the extrinsic curvature (or canonical momentum) vanishes.

Three dimensional (3D) numerical evolutions of static black holes with excision are presented. These evolutions extend to about 8000M, where M is the mass of the black hole. This degree of stability is achieved by using growth-rate estimates to guide the fine tuning of the parameters in a multiparameter family of symmetric hyperbolic representations of the Einstein evolution equations. These evolutions were performed using a fixed gauge in order to separate the intrinsic stability of the evolution equations from the effects of stability-enhancing gauge choices.

The Einstein evolution equations may be written in a variety of equivalent analytical forms, but numerical solutions of these different formulations display a wide range of growth rates for constraint violations. For symmetric hyperbolic formulations of the equations, an exact expression for the growth rate is derived using an energy norm. This expression agrees with the growth rate determined by numerical solution of the equations. An approximate method for estimating the growth rate is also derived. This estimate can be evaluated algebraically from the initial data, and is shown to exhibit qualitatively the same dependence as the numerically-determined rate on the parameters that specify the formulation of the equations. This simple rate estimate therefore provides a useful tool for finding the most well-behaved forms of the evolution equations.

First-order hyperbolic systems are promising as a basis for numerical integration of Einstein's equations. In previous work, the lapse and shift have typically not been considered part of the hyperbolic system and have been prescribed independently. This can be expensive computationally, especially if the prescription involves solving elliptic equations. Therefore, including the lapse and shift in the hyperbolic system could be advantageous for numerical work. In this paper, two first-order symmetrizable hyperbolic systems are presented that include the lapse and shift as dynamical fields and have only physical characteristic speeds.

We compare the results of constructing binary black hole initial data with three different decompositions of the constraint equations of general relativity. For each decomposition we compute the initial data using a superposition of two Kerr-Schild black holes to fix the freely specifiable data. We find that these initial-data sets differ significantly, with the ADM energy varying by as much as 5% of the total mass. We find that all initial-data sets currently used for evolutions might contain unphysical gravitational radiation of the order of several percent of the total mass. This is comparable to the amount of gravitational-wave energy observed during the evolved collision. More astrophysically realistic initial data will require more careful choices of the freely specifiable data and boundary conditions for both the metric and extrinsic curvature. However, we find that the choice of extrinsic curvature affects the resulting data sets more strongly than the choice of conformal metric.

We write a first order symmetric hyperbolic system coupling the Riemann with the dynamical a