We show that by adding suitable lower-order terms to the Z4 formulation of the Einstein equations, all constraint violations except constant modes are damped. This makes the Z4 formulation a particularly simple example of a λ-system as suggested by Brodbeck et al. We also show that the Einstein equations in harmonic coordinates can be obtained from the Z4 formulation by a change of variables that leaves the implied constraint evolution system unchanged. Therefore the same method can be used to damp all constraints in the Einstein equations in harmonic gauge.

In the Cauchy problem of general relativity one considers initial data that satisfies certain constraints. The evolution equations guarantee that the evolved variables will satisfy the constraints at later instants of time. This is only true within the domain of dependence of the initial data. If one wishes to consider situations where the evolutions are studied for longer intervals than the size of the domain of dependence, as is usually the case in three dimensional numerical relativity, one needs to give boundary data. The boundary data should be specified in such a way that the constraints are satisfied everywhere, at all times. In this paper we address this problem for the case of general relativity linearized around Minkowski space using the generalized Einstein-Christoffel symmetric hyperbolic system of evolution equations. We study the evolution equations for the constraints, specify boundary conditions for them that make them well posed and further choose these boundary conditions in such a way that the evolution equations for the metric variables are also well posed. We also consider the case of a manifold with a non-smooth boundary, as is the usual case of the cubic boxes commonly used in numerical relativity. The techniques discussed should be applicable to more general cases, as linearizations around more complicated backgrounds, and may be used to establish well posedness in the full non-linear case.

This is the first paper in a series aimed to implement boundary conditions consistent with the constraints' propagation in 3D unconstrained numerical relativity. Here we consider spherically symmetric black hole spacetimes in vacuum or with a minimally coupled scalar field, within the Einstein-Christoffel (EC) symmetric hyperbolic formulation of Einstein's equations. By exploiting the characteristic propagation of the main variables and constraints, we are able to single out the only free modes at the outer boundary for these problems. In the vacuum case a single free mode exists which corresponds to a gauge freedom, while in the matter case an extra mode exists which is associated with the scalar field. We make use of the fact that the EC formulation has no superluminal characteristic speeds to excise the singularity.We present a second-order, finite difference discretization to treat these scenarios, where we implement these constraint-preserving boundary conditions, and are able to evolve the system for essentially unlimited times. As a test of the robustness of our approach, we allow large pulses of gauge and scalar field to enter the domain through the outer boundary. We reproduce expected results, such as trivial (in the physical sense) evolution in the vacuum case (even in gauge-dynamical simulations), and the tail decay for the scalar field.

We extend the notion of numerical stability of finite difference approximations to include hyperbolic systems that are first order in time and second order in space, such as those that appear in numerical relativity and, more generally, in Hamiltonian formulations of field theories. By analyzing the symbol of the second order system, we obtain necessary and sufficient conditions for stability in a discrete norm containing one-sided difference operators. We prove stability for certain toy models and the linearized Nagy-Ortiz-Reula formulation of Einstein's equations.We also find that, unlike in the fully first order case, standard discretizations of some well-posed problems lead to unstable schemes and that the Courant limits are not always simply related to the characteristic speeds of the continuum problem. Finally, we propose methods for testing stability for second order in space hyperbolic systems.

We discuss finite difference techniques for hyperbolic equations in non-trivial domains, as those that arise when simulating black hole spacetimes. In particular, we construct dissipative and difference operators that satisfy the summation by parts property in domains with excised multiple cubic regions. This property can be used to derive semi-discrete energy estimates for the associated initial-boundary value problem which in turn can be used to prove numerical stability. I. INTRODUCTIONAn important problem in astrophysics is to model in a detailed way the collision of two black holes [1,2]. This requires numerically integrating Einstein's field equations and extracting from the simulations the relevant physical information. Unfortunately, it is difficult to obtain numerical solutions of these equations in generic three-dimensional settings, especially for long term simulations. Obstacles to this goal are encountered in both the analytical and numerical arenas. In the analytical one, the formulation of a well posed initial-boundary value problem is not completely understood. This includes the definition of proper initial and boundary conditions and the equations determining the future evolution of the fields. In the numerical arena one seeks to define numerical techniques that allow for long term accurate evolutions. This requires the construction of appropriate discrete operators to implement the initial-boundary value problem. To date, despite considerable advances in both fronts [23], the challenge of simulating generic three-dimensional black hole systems remains unattained.This article intends to provide some initial steps for setting up numerical techniques suitable to address the numerical stability of equations like the ones in question, by extending and devising finite difference techniques to tackle first order symmetric hyperbolic problems in non-trivial domains, with numerical stability being guaranteed in the linear case. Furthermore, via a local argument, one can assert that these methods should be useful in evolving smooth solutions of quasi-linear symmetric hyperbolic equations as well, as is the case of the full, non-linear Einstein vacuum equations when appropriately written [5]. Although the main motivation of this work is to present techniques for the simulation of Einstein's equations on domains with excised regions, the techniques here presented are readily applicable to any symmetric hyperbolic problem in such domains. Applications of these techniques in a variety of scenarios will be presented elsewhere [6,7,8,9,10].This work is organized in the following way. In sections II and III we review some of the issues involved in obtaining stable numerical schemes through the energy method, to set the stage for the specialized discussions that follow. In section IV, the main new results of the article are presented: we derive three dimensional difference operators satisfying summation by parts for non-trivial domains which enable one to obtain energy estimates. We further introduce dissipative ...

We use rigorous techniques from numerical analysis of hyperbolic equations in bounded domains to construct stable finite-difference schemes for Numerical Relativity, in particular for their use in black hole excision. As an application, we present 3D simulations of a scalar field propagating in a Schwarzschild black hole background.The numerical implementation of Einstein's equations represents a daunting task. The involved nature of the equations themselves, and a number of technical difficulties (related to the necessarily finite computational domain, the limitations in relative computational power and the need to deal with singularities, constraints and gauge freedom) imply a significant challenge.Numerical solutions of Einstein's equations involve solving a nonlinear set of partial differential equations on a bounded domain, and formally constitute an initialboundary-value problem (IBVP). Constructing stable and long term well behaved numerical approximations for such systems with boundaries is highly non-trivial. Here the term numerical stability is used in the sense that is equivalent, through Lax's theorem, to convergence (that is, the property that the numerical solution will approach the continuum one when resolution is increased). Delicate complications arise due to corner and edges at outer (and possible inner) boundaries, all of which introduce subtleties for a stable implementation.Recently, however, several sophisticated tools of rigorous numerical analysis have been developed for systematically constructing stable numerical schemes for IBVP's. They employ a discrete form of well-posedness and thus are stable by construction, at least for linear systems. At this time, these tools have essentially not been used by the numerical relativity community. The purpose of this paper is to outline their use, in particular for black hole excision.An IBVP consists of three ingredients: a partial differential equation, initial and boundary data. It is wellposed if a solution exists, is unique, and depends continuously on the initial and boundary data. It is well known that stable finite difference schemes approximating an IBVP can only be constructed for well-posed systems. While problems in general relativity are typically nonlinear we consider here the linear IBVP, as stability in the linearized case is a necessary condition for stability in the full nonlinear system and these methods may also be applied to nonlinear equations. Consider the linear IBVP on a domain [0, ∞) × Ωwhere u is a vector-valued function, w + and w − are incoming and outgoing modes, and S is sufficiently small (maximally dissipative boundary conditions [3]). The system is assumed to be symmetric hyperbolic, i.e., there exists a symmetric positive definite matrix H(t, x), the symmetrizer, satisfying HA i = (HA i ) T . The usual proof of well-posedness proceeds by defining an "energy", E = Ω u T Hu d 3 x , and by showing that E can be bounded as a function of the initial and boundary data. Analogously, a way to construct stable numerical s...

It is often the case in numerical relativity that schemes that are known to be convergent for well posed systems are used in evolutions of weakly hyperbolic (WH) formulations of Einstein's equations. Here we explicitly show that with several of the discretizations that have been used through out the years, this procedure leads to non-convergent schemes. That is, arbitrarily small initial errors are amplified without bound when resolution is increased, independently of the amount of numerical dissipation introduced. The lack of convergence introduced by this instability can be particularly subtle, in the sense that it can be missed by several convergence tests, especially in 3+1 dimensional codes. We propose tests and methods to analyze convergence that may help detect these situations.Convergence is a central element of any numerical simulation. It refers to the property that if one refines the simulation by adding more points to the grid, numerical errors should diminish. In the limit of zero spacing, they should go to zero and one should get the exact solution. Having a convergent code is a key element for numerical simulations to have predictive power: although one will always be limited in practice to a finite number of points in the grid, one can extrapolate the results for more and more refined simulations and have a very good approximation to the true (continuum) results. Codes that do not converge produce answers that, even if they remain finite (at least for a while), have little predictive power: there generically is no way to know if the solutions found approximate the desired continuum solution.In this paper we want to emphasize that discretization schemes that yield convergent code for strongly hyperbolic (SH) systems of equations do not necessarily do so for weakly hyperbolic (WH) or completely ill posed (CIP) systems. The relevance of this observation is that several formulations of the Einstein evolution equations commonly used in numerical relativity, including, for example, the Arnowitt-Deser-Misner (ADM) [1] and Baumgarte-Shapiro-Shibata-Nakamura (BSSN) [2] formulations with fixed lapse and shift, are not SH.We concretely prove that a simple system of linear WH equations with constant coefficients, when discretized with either the iterated Crank-Nicholson (ICN) with fixed number of iterations or second order Runge-Kutta (2RK) methods, leads to unconditionally unstable codes, even if numerical dissipation is explicitly added. This simple system of equations is related to the equations one encounters in the WH formulations of general relativity. It therefore strongly suggests that these formulations should produce code that is not convergent. The lack of convergence is of a particularly pernicious nature, since in the WH case it might grow slower than the "usual" von-Neumann numerical instability. It is such that if one tests the code with stationary solutions (for instance simulating a single black hole) and performs convergence tests, one could easily be confused into believing one has conve...

It is expected that the realization of a convergent and long-term stable numerical code for the simulation of a black hole inspiral collision will depend greatly upon the construction of stable algorithms capable of handling smooth and, most likely, time dependent boundaries. After deriving single grid, energy conserving discretizations for axisymmetric systems containing the axis of symmetry, we present a new excision method for moving black holes using multiple overlapping coordinate patches, such that each boundary is fixed with respect to at least one coordinate system. This multiple coordinate structure eliminates all need for extrapolation, a commonly used procedure for moving boundaries in numerical relativity.We demonstrate this excision method by evolving a massless Klein-Gordon scalar field around a boosted Schwarzschild black hole in axisymmetry. The excision boundary is defined by a spherical coordinate system co-moving with the black hole. Our numerical experiments indicate that arbitrarily high boost velocities can be used without observing any sign of instability.