Table of contents

Traditionally, frontiers represent a treacherous terrain to venture into, where hidden obstacles are present and uncharted territories lie ahead. At the same time, frontiers are also a place where new perspectives can be appreciated and have often been the cradle of new and thriving developments. With this in mind and inspired by this spirit, the Numerical Relativity Group at the Albert Einstein Institute (AEI) organized a `New Frontiers in Numerical Relativity' meeting on 17–21 July 2006 at the AEI campus in Potsdam, Germany.

It is an interesting historical remark that the suggestion of the meeting was first made in the late summer of 2005 and thus at a time that for many reasons has been a turning point in the recent history of numerical relativity. A few months earlier (April 2005) in fact, F Pretorius had announced the first multi-orbit simulations of binary black holes and computed the waveforms from the inspiral, merger and ring-down (`Numerical Relativity', Banff International Research Station, Banff, Canada, 16–21 April 2005). At that time, the work of Pretorius served as an important boost to the research in this field and although no other group has yet adopted the techniques he employed, his results provided the numerical relativity community with clear evidence that the binary black hole problem could be solved. A few months later (November 2005), equally striking results were presented by the NASA Goddard and Texas/Brownsville groups, who also reported, independently, multi-orbit evolutions of binary black holes using numerical techniques and formulations of the Einstein equations which were markedly distinct from those suggested by Pretorius (`Numerical Relativity 2005', Goddard Space Flight Centre, Greenbelt, MD, USA, 2–4 November 2005).

A few months later other groups were able to repeat the same simulations and obtain equivalent results, testifying that the community as a whole had reached comparable levels of maturity in both the numerical techniques and the mathematical methods needed for successful solution of the Einstein equations for binary black holes. Clearly, an important frontier, and actually a long-awaited one, was finally open and the `gold rush' was just about to begin by the time the `New Frontiers in Numerical Relativity' meeting started its sessions in July 2006.

And so, almost 20 years since the almost homonymous meeting held at Urbana–Champaign (`Frontiers in Numerical Relativity', University of Illinois, IL, USA, 1988), the `New Frontiers in Numerical Relativity' meeting at the AEI saw the enthusiastic participation of a great part of the community, with 127 participants present (in 1988 they were 55) and with a large majority being represented by students and postdocs, a reassuring sign of good health for the community. Faithful to the title of the conference, the programme was dedicated to the many and diversified `frontiers' in numerical relativity and organized so as to have few talks with ample time dedicated to discussions.

Overall, the talks presented at the meeting covered all of the most salient aspects of numerical relativity: from the formulation of the Einstein equations, over to the initial-value problem in general relativity, from the evolution of vacuum and non-vacuum spacetimes, to multiblock adaptive mesh-refinement techniques, from boundary conditions and perturbative methods, to relativistic fluids and plasmas. The contributions in this special issue represent a selection of that research, but also include invited papers from authors who were not present at the meeting but were pursuing research at the forefronts of numerical relativity.

In addition to the more traditional sessions, the `New Frontiers in Numerical Relativity' meeting also hosted a less traditional session, dedicated to an `unconstrained' discussion which covered some of the most controversial issues that emerged during the conference. During this session, chaired by E Seidel, a lively discussion took place in the non-trivial attempt of marking the new frontiers on the map of numerical relativity. The transcript of this discussion is an integral part of this issue and it is available, along with the audio recording, in the online version only. We believe they embody an important part of the development of this field and, like a good bottle of wine, it will be interesting to read them again once sufficiently aged.

As a concluding remark we note that it is almost one year since the `New Frontiers in Numerical Relativity' meeting and dozens of excellent papers have been published or posted on the preprint archive. Some of the scientific results obtained over these months, especially those revolving around binary black holes, were simply unimaginable a few years ago and represent an indisputable evidence that the research in numerical relativity has never been as exciting as it is now.

These results have already had an impact in astrophysics and the community interested in the analysis of gravitational-wave data, thus opening new and different frontiers in numerical relativity. Interestingly, all of this is happening while ground-based gravitational wave detectors in the US and Europe are operating at a sensitivity such that gravitational radiation may soon be directly detected.

While much still needs to be understood and improved, the gold rush towards the new frontiers of numerical relativity does not yet show any sign of being close to a rapid end.

We present a multi-domain spectral method to compute initial data of binary systems in general relativity. By utilizing adapted conformal coordinates, the vacuum region exterior to the gravitational sources is divided into two subdomains within which the spectral expansion of the field quantities is carried out. If a component of the binary is a neutron star, a further subdomain covering the star's interior is added. As such, the method can be used to construct arbitrary initial data corresponding to binary black holes, binary neutron stars or mixed binaries. In particular, it is possible to describe a black hole component by the puncture ansatz as well as through an excision technique. First examples are given for binary black hole excision data that fulfil the requirements of the quasi-stationary framework, which combines the conformal thin sandwich formulation of the constraint equations with the isolated horizon conditions for black holes in quasi-equilibrium. These numerical solutions were obtained to extremely high accuracy with moderate computational effort. Moreover, the method proves to be applicable even when tending towards limiting cases such as large mass ratios of the binary components.

It is well known that Bowen–York initial data contain spurious radiation. Although this 'junk' radiation has been seen to be small for non-spinning black-hole binaries in circular orbit, its magnitude increases when the black holes are given spin. It is possible to reduce the spurious radiation by applying the puncture approach to multiple Kerr black holes, as we demonstrate for examples of head-on collisions of equal-mass black-hole binaries.

This is a particularly exciting time for gravitational wave physics. Ground-based gravitational wave detectors are now operating at a sensitivity such that gravitational radiation may soon be directly detected, and recently several groups have independently made significant breakthroughs that have finally enabled numerical relativists to solve the Einstein field equations for coalescing black-hole binaries, a key source of gravitational radiation. The numerical relativity community is now in the position to begin providing simulated merger waveforms for use by the data analysis community, and it is therefore very important that we provide ways to validate the results produced by various numerical approaches. Here, we present a simple comparison of the waveforms produced by two very different, but equally successful approaches—the generalized harmonic gauge and the moving puncture methods. We compare waveforms of equal-mass black hole mergers with minimal or vanishing spins. The results show exceptional agreement for the final burst of radiation, with some differences attributable to small spins on the black holes in one case.

We present results from fully nonlinear simulations of unequal mass binary black holes plunging from close separations well inside the innermost stable circular orbit with mass ratios q ≡ M₁/M₂ = {1, 0.85, 0.78, 0.55, 0.32}, or equivalently, with reduced mass parameters η ≡ M₁M₂/(M₁ + M₂)² = {0.25, 0.248, 0.246, 0.229, 0.183}. For each case, the initial binary orbital parameters are chosen from the Cook–Baumgarte equal-mass ISCO configuration. We show waveforms of the dominant ℓ = 2, 3 modes and compute estimates of energy and angular momentum radiated. For the plunges from the close separations considered, we measure kick velocities from gravitational radiation recoil in the range 25–82 km s⁻¹. Due to the initial close separations our kick velocity estimates should be understood as a lower bound. The close configurations considered are also likely to contain significant eccentricities influencing the recoil velocity.

Binary black hole simulations have traditionally been computationally very expensive: current simulations are performed in supercomputers involving dozens if not hundreds of processors, thus systematic studies of the parameter space of binary black hole encounters still seem prohibitive with current technology. Here we show how the multi-layered refinement level code BAM can be used on dual processor workstations to simulate certain binary black hole systems. BAM, based on the moving punctures method, provides grid structures composed of boxes of increasing resolution near the centre of the grid. In the case of binaries, the highest resolution boxes are placed around each black hole and they track them in their orbits until the final merger when a single set of levels surrounds the black hole remnant. This is particularly useful when simulating spinning black holes since the gravitational fields gradients are larger. We present simulations of binaries with equal mass black holes with spins parallel to the binary axis and intrinsic magnitude of S/m² = 0.75. Our results compare favourably to those of previous simulations of this particular system. We show that the moving punctures method produces stable simulations at maximum spatial resolutions up to M/160 and for durations of up to the equivalent of 20 orbital periods.

Binary black hole simulations starting from quasi-circular (i.e., zero radial velocity) initial data have orbits with small but nonzero orbital eccentricities. In this paper, the quasi-equilibrium initial-data method is extended to allow nonzero 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 (∼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 Phys. Rev.136 1224–32). The gravitational waveforms, which contain ∼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 describe recent numerical simulations of the merger of a class of equal mass, non-spinning, eccentric binary black hole systems in general relativity. We show that with appropriate fine tuning of the initial conditions one can reach a region of parameter space we denote the threshold of immediate merger. Here, the binary enters a phase of close interaction in a near-circular orbit, stays there for an amount of time proportional to the logarithmic distance from the threshold in parameter space, then either separates or merges to form a single Kerr black hole. To gain a better understanding of this phenomenon, we study an analogous problem in the evolution of equatorial geodesics about a central Kerr black hole. A similar threshold of capture exists for appropriate classes of initial conditions, and tuning to threshold the geodesics approach one of the unstable circular geodesics of the Kerr spacetime. Remarkably, with a natural mapping of the parameters of the geodesic to that of the equal mass system, the scaling exponents describing the whirl phase of each system turn out to be quite similar. Armed with this lone piece of evidence that an approximate correspondence might exist between near-threshold evolution of geodesics and generic binary mergers, we illustrate how this information can be used to estimate the cross section and energy emitted in the ultra-relativistic black hole scattering problem. This could eventually be of use in providing estimates for the related problem of parton collisions at the large hadron collider in extra dimension scenarios where black holes are produced.

We discuss the transition from quasi-circular inspiral to plunge of a system of two nonrotating black holes of masses m₁ and m₂ in the extreme-mass-ratio limit m₁m₂ ≪ (m₁ + m₂)². In the spirit of the effective one body (EOB) approach to the general relativistic dynamics of binary systems, the dynamics of the two black hole system is represented in terms of an effective particle of mass μ ≡ m₁m₂/(m₁ + m₂) moving in a (quasi-)Schwarzschild background of mass M ≡ m₁ + m₂ and submitted to an radiation reaction force defined by Padé resumming high-order post-Newtonian results. We then complete this approach by numerically computing, in the style of Regge–Wheeler–Zerilli, the gravitational radiation emitted by such a particle. Several tests of the numerical procedure are presented. We focus on gravitational waveforms and the related energy and angular momentum losses. We view this work as a contribution to the matching between analytical and numerical methods within an EOB-type framework.

We present our latest results of simulation for merger of black hole (BH)–neutron star (NS) binaries in full general relativity which is performed preparing a quasicircular state as the initial condition. The BH is modelled by a moving puncture with no spin and the NS by the Γ-law equation of state with Γ = 2 and a corotating velocity field as a first step. The mass of the BH is chosen to be ≈3.2M_⊙ or 4.0M_⊙, and the rest mass of the NS ≈1.4M_⊙ with a relatively large radius of the NS ≈13–14 km. The NS is tidally disrupted near the innermost stable orbit, but ∼80–90% of the material is swallowed into the BH and the resulting disc mass is not very large as ∼0.3M_⊙ even for a small BH mass ∼3.2M_⊙. The result indicates that the system composed of BH and a massive disc of ∼M_⊙ is not formed from nonspinning BH–NS binaries irrespective of the BH mass, although a disc of mass ∼0.1M_⊙ is a possible outcome for this relatively small BH mass range as ∼3–4M_⊙. Our results indicate that the merger of low-mass BH and NS may form a central engine of short-gamma-ray bursts.

We present results from the first 2 + 1 and 3 + 1 simulations of the collapse of rotating stellar iron cores in general relativity employing a finite-temperature equation of state and an approximate treatment of deleptonization during collapse. We compare full 3 + 1 and conformally-flat spacetime evolution methods and find that the conformally-flat treatment is sufficiently accurate for the core-collapse supernova problem. We focus on the gravitational wave (GW) emission from rotating collapse, core bounce and early postbounce phases. Our results indicate that the GW signature of these phases is much more generic than previously estimated. In addition, we track the growth of a nonaxisymmetric instability of dominant m = 1 character in two of our models that leads to prolonged narrow-band GW emission at ∼920–930 Hz over several tens of milliseconds.

We present a new general relativistic hydrodynamics code specifically designed to study magneto-rotational, relativistic, stellar core collapse. The code is an extension of an existing (and thoroughly tested) hydrodynamics code, which has been applied in the recent past to study relativistic rotational core collapse. It is based on the conformally-flat approximation of Einstein's field equations and conservative formulations for the magneto-hydrodynamics equations. As a first step towards magneto-rotational core collapse simulations the code assumes a passive (test) magnetic field. The paper is focused on the description of the technical details of the numerical implementation, with emphasis on the magnetic field module. A number of code tests are presented and discussed, along with a representative core collapse simulation.

We present new results on dynamical instabilities in rapidly rotating relativistic stars. In particular, using numerical simulations in full general relativity, we analyse the effects that the stellar compactness has on the threshold for the onset of the dynamical bar-mode instability, as well as on the appearance of other dynamical instabilities. By using an extrapolation technique developed and tested in our previous study (Baiotti L et al 2007 Phys. Rev. D 75 044023), we explicitly determine the threshold for a wide range of compactnesses using four sequences of models of constant baryonic mass comprising a total of 59 stellar models. Our calculation of the threshold is in good agreement with the Newtonian prediction and improves the previous post-Newtonian estimates. In addition, we find that for stars with sufficiently large mass and compactness, the m = 3 deformation is the fastest growing one. For all of the models considered, the non-axisymmetric instability is suppressed on a dynamical timescale with an m = 1 deformation dominating the final stages of the instability. These results, together with those presented in Baiotti L et al (2007 Phys. Rev. D 75 044023), suggest that an m = 1 deformation represents a general and late-time feature of non-axisymmetric dynamical instabilities both in full general relativity and in Newtonian gravity.

We provide details and present additional results on the numerical study of the gravitational-wave emission from the collapse of neutron stars to rotating black holes in three dimensions. More specifically, we concentrate on the advantages and disadvantages of the use of the excision technique and on how alternative approaches to that of excision can be successfully employed. Furthermore, as a first step towards source characterization, we present a systematic discussion of the influence that rotation and different perturbations have on the waveforms and hence on the energy emitted in gravitational waves.

The capacity to model magnetohydrodynamical (MHD) flows in dynamical, strongly curved spacetimes significantly extends the reach of numerical relativity in addressing many problems at the forefront of theoretical astrophysics. We have developed and tested an evolution code for the coupled Einstein–Maxwell-MHD equations which combines a BSSN solver with a high resolution shock capturing scheme. As one application, we evolve magnetized, differentially rotating neutron stars under the influence of a small seed magnetic field. Of particular significance is the behaviour found for hypermassive neutron stars (HMNSs), which have rest masses greater than the mass limit allowed by uniform rotation for a given equation of state. The remnant of a binary neutron star merger is likely to be a HMNS. We find that magnetic braking and the magnetorotational instability lead to the collapse of HMNSs and the formation of rotating black holes surrounded by massive, hot accretion tori and collimated magnetic field lines. Such tori radiate strongly in neutrinos, and the resulting neutrino–antineutrino annihilation (possibly in concert with energy extraction by MHD effects) could provide enough energy to power short-hard gamma-ray bursts. To explore the range of outcomes, we also evolve differentially rotating neutron stars with lower masses and angular momenta than the HMNS models. Instead of collapsing, the non-hypermassive models form nearly uniformly rotating central objects which, in cases with significant angular momentum, are surrounded by massive tori.

A flexible spectral code for the study of general relativistic magnetohydrodynamics is presented. Aiming at investigating the physics of slowly rotating magnetized compact stars, this new code makes use of various physically motivated approximations. Among them, the relativistic anelastic approximation is a key ingredient of the current version of the code. In this paper, we mainly outline the method, putting emphasis on algorithmic techniques that enable to benefit as much as possible from the non-dissipative character of spectral methods, showing also a potential astrophysical application and providing a few illustrative tests.

The accurate modelling of astrophysical scenarios involving compact objects and magnetic fields, such as the collapse of rotating magnetized stars to black holes or the phenomenology of γ-ray bursts, requires the solution of the Einstein equations together with those of general-relativistic magnetohydrodynamics. We present a new numerical code developed to solve the full set of general-relativistic magnetohydrodynamics equations in a dynamical and arbitrary spacetime with high-resolution shock-capturing techniques on domains with adaptive mesh refinements. After a discussion of the equations solved and of the techniques employed, we present a series of testbeds carried out to validate the code and assess its accuracy. Such tests range from the solution of relativistic Riemann problems in flat spacetime, over to the stationary accretion onto a Schwarzschild black hole and up to the evolution of oscillating magnetized stars in equilibrium and constructed as consistent solutions of the coupled Einstein–Maxwell equations.

The radio source Sagittarius A* (Sgr A*) is believed to be a hot, inhomogeneous, magnetized plasma flowing near the event horizon of the 3.6 × 10⁶ M_⊙ black hole at the galactic centre. At a distance of 8 kpc (≃ 2.5 × 10²² cm) the black hole would be among the largest black holes as judged by angular size. Recent observations are consistent with the idea that the millimetre and sub-millimetre photons are dominated by optically thin, thermal synchrotron emission. Anticipating future Very Long Baseline Interferometry (VLBI) observations of Sgr A* at these wavelengths, we present here the first dynamically self-consistent models of millimetre and sub-millimetre emission from Sgr A* based on general relativistic numerical simulations of the accretion flow. Angle-dependent spectra are calculated assuming a thermal distribution of electrons at the baryonic temperature dictated by the simulation and the accretion rate, which acts as a free parameter in our model. The effects of varying model parameters (black hole spin and inclination of the spin to the line of sight) and source variability on the spectrum are shown. We find that the accretion rate value needed to match our calculated millimetre flux to the observed flux is consistent with constraints on the accretion rate inferred from detections of the rotation measure. We also describe the relativistic jet that is launched along the black hole spin axis by the accretion disc and evolves to scales of ∼10³GMc⁻², where M is the mass of the black hole.

We describe an explicit in time, finite-difference code designed to simulate black holes by using the excision method. The code is based upon the harmonic formulation of the Einstein equations and incorporates several features regarding the well-posedness and numerical stability of the initial-boundary problem for the quasilinear wave equation. After a discussion of the equations solved and of the techniques employed, we present a series of testbeds carried out to validate the code. Such tests range from the evolution of isolated black holes to the head-on collision of two black holes and then to a binary black hole inspiral and merger. Besides assessing the accuracy of the code, the inspiral and merger test has revealed that the marginally trapped surfaces contained within the common apparent horizon of the merged black hole can touch and even intersect. This novel feature in the dynamics of the marginally trapped surfaces is unexpected but consistent with theorems on the properties of these surfaces.

This paper covers some of the current techniques and issues involved in performing numerical simulations of the formation of singularities.

In a recent article, we constructed a hierarchy of outer boundary conditions for Einstein's field equations with the property that, for a spherical outer boundary, it is perfectly absorbing for linearized gravitational radiation up to a given angular momentum number L. In this paper, we generalize so that it can be applied to fairly general foliations of spacetime by space-like hypersurfaces and general outer boundary shapes and further, we improve in two steps: (i) we give a local boundary condition which is perfectly absorbing including first-order contributions in 2M/R of curvature corrections for quadrupolar waves (where M is the mass of the spacetime and R is a typical radius of the outer boundary) and which significantly reduces spurious reflections due to backscatter, and (ii) we give a non-local boundary condition which is exact when first-order corrections in 2M/R for both curvature and backscatter are considered, for quadrupolar radiation.

The characteristic approach to numerical relativity is a useful tool in evolving gravitational systems. In the past this has been implemented using two patches of stereographic angular coordinates. In other applications, a six-patch angular coordinate system has proved effective. Here we investigate the use of a six-patch system in characteristic numerical relativity, by comparing an existing two-patch implementation (using second-order finite differencing throughout) with a new six-patch implementation (using either second- or fourth-order finite differencing for the angular derivatives). We compare these different codes by monitoring the Einstein constraint equations, numerically evaluated independently from the evolution. We find that, compared to the (second-order) two-patch code at equivalent resolutions, the errors of the second-order six-patch code are smaller by a factor of about 2, and the errors of the fourth-order six-patch code are smaller by a factor of nearly 50.

We present a method for extracting gravitational waves from numerical spacetimes which generalizes and refines one of the standard methods based on the Regge–Wheeler–Zerilli perturbation formalism. At the analytical level, this generalization allows a much more general class of slicing conditions for the background geometry, and is thus not restricted to Schwarzschild-like coordinates. At the numerical level, our approach uses high-order multi-block methods, which improve both the accuracy of our simulations and of our extraction procedure. In particular, the latter is simplified since there is no need for interpolation, and we can afford to extract accurate waves at large radii with only little additional computational effort. We then present fully nonlinear three-dimensional numerical evolutions of a distorted Schwarzschild black hole in Kerr–Schild coordinates with an odd parity perturbation and analyse the improvement that we gain from our generalized wave extraction, comparing our new method to the standard one. In particular, we analyse in detail the quasinormal frequencies of the extracted waves, using both methods. We do so by comparing the extracted waves with one-dimensional high resolution solutions of the corresponding generalized Regge–Wheeler equation. We explicitly see that the errors in the waveforms extracted with the standard method at fixed, finite extraction radii do not converge to zero with increasing resolution. We find that even with observers as far out as R = 80M—which is larger than what is commonly used in state-of-the-art simulations—the assumption in the standard method that the background is close to having Schwarzschild-like coordinates increases the error in the extracted waves considerably. Furthermore, those errors are dominated by the extraction method itself and not by the accuracy of our simulations. For extraction radii between 20M and 80M and for the resolutions that we use in this paper, our new method decreases the errors in the extracted waves, compared to the standard method, by between one and three orders of magnitude. In a general scenario, for example a collision of compact objects, there is no precise definition of gravitational radiation at a finite distance, and gravitational wave extraction methods at such distances are thus inherently approximate. The results of this paper bring up the possibility that different choices in the wave extraction procedure at a fixed and finite distance may result in relative differences in the waveforms which are actually larger than the numerical errors in the solution.

Time evolutions of certain spherically symmetric systems are investigated where simple explicit second order finite difference methods are not applicable. Due to a compactified space coordinate, efficiency and long-term numerical stability require at least fourth order accuracy for both the massive Klein–Gordon field and the SU(2) Yang–Mills–Higgs system. Moreover, adaptive mesh refinement (AMR) has a crucial role in dealing with high frequency oscillations that appear as an initial disturbance is radiated away. The incompatibility of AMR with fully fourth order accuracy is discussed and a solution is presented. Finally, compactification is compared to standard spherical coordinates and truncated grids in terms of efficiency.