Table of contents

Two meeetings gave rise to this special issue on numerical relativity: the workshop 'Numerical relativity' at the Banff International Research Station on 16–21 April 2005 and the conference 'New directions in numerical relativity' which was held at Southampton University on the 18 and 19 August 2005 as a satellite meeting of the Newton Institute Programme 'Global problems in mathematical relativity'. This edition contains contributions drawn from these two meetings.

Looking back, 2005 will be remembered as the year in which key advances were made on a number of fronts which allowed significant progress in the binary black hole merger problem: at the Banff meeting, Frans Pretorius announced the first multi-orbit simulations, using a generalization of harmonic coordinates in which Friedrich's gauge source functions have been promoted to dynamical variables.

Then, at the 'Numerical Relativity 2005' meeting held on 2–4 November 2005 at NASA–Goddard, the NASA–Goddard and Texas/Brownsville groups independently (in back-to-back talks!) announced multi-orbit simulations with waveforms using the Baumgarte–Shapiro–Shibata–Nakamura 3+1 formulation with improved hyperbolic lapse and shift drivers, and representing the black holes as wormholes ('punctures') moving through the grid.

These highlights were made possible by previous progress. Particularly important is the implementation of adaptive mesh refinement in general relativity in two and three dimensions, which not only allows for improved accuracy, but reduces the amount of time taken by 3D simulations, thus allowing systematic testing and improvement of 3D codes.

In addition, the community is now much more aware of the importance of well-posedness of the continuum problem and the stability of the numerical methods, and some formal investigations of these matters have caused practical improvements. The same applies for the role of gauge choices and boundary conditions.

Beyond the binary black hole problem, more incremental but steady progress is being made in neutron star merger and collapse simulations. The general relativity side of the simulations is just coming under control, and groups are now implementing more realistic matter models and, in particular, magneto-hydrodynamics. Vacuum and matter simulations will remain closely related. The investigation of critical collapse continues in 1D and 2D and remains a challenge for 3D codes. Additionally, work continues towards incorporating more refined numerical techniques for improved accuracy, stability and robustness.

Some, but not all, of these developments are reflected in the papers collected in this issue. They were invited mainly from participants of the two meetings above and make up this special issue which is the first of several that Classical and Quantum Gravity will publish. It is clear that many significant results are to be expected in the near future and will be reflected in this series.

Building on previous work on the critical behaviour in gravitational collapse of the self-gravitating SU(2) σ-field and using high precision numerical methods we uncover a fine structure hidden in a narrow window of parameter space. We argue that this numerical finding has a natural explanation within a dynamical system framework of critical collapse.

The large scale binary black hole effort in numerical relativity has led to increasing distinction between numerical and mathematical relativity. This paper discusses this situation and gives some examples of successful interactions between numerical and mathematical models in general relativity.

Computational methods are essential to provide waveforms from coalescing black holes, which are expected to produce strong signals for the gravitational wave observatories being developed. Although partial simulations of the coalescence have been reported, scientifically useful waveforms have so far not been delivered. The goal of the AppleswithApples (AwA) Alliance is to design, coordinate and document standardized code tests for comparing numerical relativity codes. The first round of AwA tests has now been completed and the results are being analysed. These initial tests are based upon the periodic boundary conditions designed to isolate performance of the main evolution code. Here we describe and carry out an additional test with periodic boundary conditions which deals with an essential feature of the black hole excision problem, namely a non-vanishing shift. The test is a shifted version of the existing AwA gauge wave test. We show how a shift introduces an exponentially growing instability which violates the constraints of a standard harmonic formulation of Einstein's equations. We analyse the Cauchy problem in a harmonic gauge and discuss particular options for suppressing instabilities in the gauge wave tests. We implement these techniques in a finite difference evolution algorithm and present test results. Although our application here is limited to a model problem, the techniques should benefit the simulation of black holes using harmonic evolution codes.

We present strongly stable semi-discrete finite difference approximations to the quarter space problem (x > 0, t > 0) for the first order in time, second order in space wave equation with a shift term. We consider space-like (pure outflow) and time-like boundaries, with either second- or fourth-order accuracy. These discrete boundary conditions suggest a general prescription for boundary conditions in finite difference codes approximating first order in time, second order in space hyperbolic problems, such as those that appear in numerical relativity. As an example we construct boundary conditions for the Nagy–Ortiz–Reula formulation of the Einstein equations coupled to a scalar field in spherical symmetry.

A major obstacle in the numerical simulation of general relativistic spacetimes is the fact that coordinates have to be specified in order to obtain a well-defined numerical evolution. While the choice of coordinates has no impact on the geometry, it does influence the evolution in ways which are rather difficult to predict. For this reason it might be useful to devise numerical procedures which are manifestly coordinate invariant. In this paper we present an alternative way to obtain discrete versions of the Einstein equations. We formulate the Einstein equations as an exterior system in terms of differential forms. Such a formulation can be used to interpret the variables and equations in a discrete setting as discrete differential forms. We describe the basic equations and discuss the accuracy of the discrete equations.

A possible definition of strong/symmetric hyperbolicity for a second-order system of evolution equations is that it admits a reduction to first order which is strongly/symmetric hyperbolic. We investigate the general system that admits a reduction to first order and give necessary and sufficient criteria for strong/symmetric hyperbolicity of the reduction in terms of the principal part of the original second-order system. An alternative definition of strong hyperbolicity is based on the existence of a complete set of characteristic variables, and an alternative definition of symmetric hyperbolicity is based on the existence of a conserved (up to lower-order terms) energy. Both these definitions are made without any explicit reduction. Finally, strong hyperbolicity can be defined through a pseudo-differential reduction to first order. We prove that both definitions of symmetric hyperbolicity are equivalent and that all three definitions of strong hyperbolicity are equivalent (in three space dimensions). We show how to impose maximally dissipative boundary conditions on any symmetric hyperbolic second-order system. We prove that if the second-order system is strongly hyperbolic, any closed constraint evolution system associated with it is also strongly hyperbolic, and that the characteristic variables of the constraint system are derivatives of a subset of the characteristic variables of the main system, with the same speeds.

In the harmonic description of general relativity, the principal part of the Einstein equations reduces to a constrained system of ten curved space wave equations for the components of the spacetime metric. We use the pseudo- differential theory of systems which are strongly well posed in the generalized sense to establish the well posedness of constraint-preserving boundary conditions for this system when treated in a second-order differential form. The boundary conditions are of a generalized Sommerfeld type that is benevolent for numerical calculation.

Efforts to achieve better accuracy in numerical relativity have so far focused either on implementing second-order accurate adaptive mesh refinement or on defining higher order accurate differences and update schemes. Here, we argue for the combination, that is a higher order accurate adaptive scheme. This combines the power that adaptive gridding techniques provide to resolve fine scales (in addition to a more efficient use of resources) together with the higher accuracy furnished by higher order schemes when the solution is adequately resolved. To define a convenient higher order adaptive mesh refinement scheme, we discuss a few different modifications of the standard, second-order accurate approach of Berger and Oliger. Applying each of these methods to a simple model problem, we find these options have unstable modes. However, a novel approach to dealing with the grid boundaries introduced by the adaptivity appears stable and quite promising for the use of high order operators within an adaptive framework.

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.

For the case of an electric or scalar charge travelling along a timelike geodesic in a curved, background spacetime, we apply a Green's theorem argument to transform the tail contribution to the particle's self-force at a point p along its trajectory to a form which involves only an integral over the past light cone from p. One potential advantage of this reformulation is that the tail contribution to the fundamental solution for a (self-adjoint) tensor wave equation in an arbitrary spacetime is (at least within so-called causal domains) explicitly computable on the past light cone of an arbitrary point through the integration of certain linear transport equations defined along the null generators of the cone. The conventional approach to computing the self-force requires this tail field throughout the interior of the past light cone and extracts it from the particle's (numerically generated) total field through an intricate mode-by-mode numerical decomposition. By contrast, our approach requires that the particle's total field be paired with the (explicitly computable) tail field on the light cone itself and then integrated over this cone. We thus avoid the need for the aforementioned mode-by-mode decomposition of the particle's total field into direct ⊕ tail contributions. We speculate that in some circumstances it may be sufficient to truncate the particle's total field, as needed in our calculation, to its (also explicitly computable) direct Liénard Wiechert approximation.

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, we show 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 Ψ₀ yields a priori estimates on the metric, connection and curvature fields. These estimates are expected to be useful for obtaining 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, an a priori specified shift vector 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.

A new algorithm for solving the general relativistic MHD equations is described in this paper. We design our scheme to incorporate black hole excision with smooth boundaries, and to simplify solving the combined Einstein and MHD equations with AMR. The fluid equations are solved using a finite difference convex ENO method. Excision is implemented using overlapping grids. Elliptic and hyperbolic divergence cleaning techniques allow for maximum flexibility in choosing coordinate systems, and we compare both methods for a standard problem. Numerical results of standard test problems are presented in two-dimensional flat space using excision, overlapping grids and elliptic and hyperbolic divergence cleaning.

A numerical solution scheme for the Einstein field equations based on generalized harmonic coordinates is described, focusing on details which are not provided before in the literature and which are of particular relevance to the binary black hole problem. This includes demonstrations of the effectiveness of constraint damping, and how the time slicing can be controlled through the use of a source function evolution equation. In addition, some results from an ongoing study of binary black hole coalescence, where the black holes are formed via scalar field collapse, are shown. Scalar fields offer a convenient route to exploring certain aspects of black hole interactions, and one interesting though tentative suggestion from this early study is that behaviour reminiscent of 'zoom-whirl' orbits in particle trajectories is also present in the merger of equal mass, non-spinning binaries, with appropriately fine-tuned initial conditions.

We describe a generic infrastructure for time evolution simulations in numerical relativity using multiple grid patches. After a motivation of this approach, we discuss the relative advantages of global and patch-local tensor bases. We describe both our multi-patch infrastructure and our time evolution scheme, and comment on adaptive time integrators and parallelization. We also describe various patch system topologies that provide spherical outer and/or multiple inner boundaries. We employ penalty inter-patch boundary conditions, and we demonstrate the stability and accuracy of our three-dimensional implementation. We solve both a scalar wave equation on a stationary rotating black hole background and the full Einstein equations. For the scalar wave equation, we compare the effects of global and patch-local tensor bases, different finite differencing operators and the effect of artificial dissipation onto stability and accuracy. We show that multi-patch systems can directly compete with the so-called fixed mesh refinement approach; however, one can also combine both. For the Einstein equations, we show that using multiple grid patches with penalty boundary conditions leads to a robustly stable system. We also show long-term stable and accurate evolutions of a one-dimensional nonlinear gauge wave. Finally, we evolve weak gravitational waves in three dimensions and extract accurate waveforms, taking advantage of the spherical shape of our grid lines.

We introduce a computational framework which avoids solving explicitly hydrodynamic equations and is suitable for studying the pre-merger evolution of black hole–neutron star binary systems. The essence of the method consists of constructing a neutron star model with a black hole companion and freezing the internal degrees of freedom of the neutron star during the course of the evolution of the spacetime geometry. We present the main ingredients of the framework, from the formulation of the problem to the appropriate computational techniques to study these binary systems. In addition, we present numerical results of the construction of initial data sets and evolutions that demonstrate the feasibility of this approach.

Results from helically symmetric scalar-field models and first results from a convergent helically symmetric binary neutron-star code are reported here; these are models stationary in the rotating frame of a source with constant angular velocity Ω. In the scalar-field models and the neutron-star code, helical symmetry leads to a system of mixed elliptic–hyperbolic character. The scalar-field models involve nonlinear terms of the form ψ³, (∇ψ)² and ψ□ψ that mimic nonlinear terms of the Einstein equation. Convergence is strikingly different for different signs of each nonlinear term; it is typically insensitive to the iterative method used, and it improves with an outer boundary in the near zone. In the neutron-star code, one has no control on the sign of the source, and convergence has been achieved only for an outer boundary less than ∼1 wavelength from the source or for a code that imposes helical symmetry only inside a near zone of that size. The inaccuracy of helically symmetric solutions with appropriate boundary conditions should be comparable to the inaccuracy of a waveless formalism that neglects gravitational waves, and the (near zone) solutions we obtain for waveless and helically symmetric BNS codes with the same boundary conditions nearly coincide.

In computational relativity, critical behaviour near the black hole threshold has been studied numerically for several models in the last decade. In this paper we present a spatial Galerkin method suitable for finding numerical solutions of the Einstein–Dirac equations in spherically symmetric spacetime (in polar/areal coordinates). The method features exact conservation of the total electric charge and allows for a spatial mesh adaption based on physical arclength. Numerical experiments confirm excellent robustness and convergence properties of our approach. Hence, this new algorithm is well suited for studying critical behaviour.