# Accurate black hole evolutions by fourth-order numerical relativity

###### Abstract

We present techniques for successfully performing numerical relativity simulations of binary black holes with fourth-order accuracy. Our simulations are based on a new coding framework which currently supports higher order finite differencing for the BSSN formulation of Einstein’s equations, but which is designed to be readily applicable to a broad class of formulations. We apply our techniques to a standard set of numerical relativity test problems, demonstrating the fourth-order accuracy of the solutions. Finally we apply our approach to binary black hole head-on collisions, calculating the waveforms of gravitational radiation generated and demonstrating significant improvements in waveform accuracy over second-order methods with typically achievable numerical resolution.

###### pacs:

04.25.Dm, 04.25.Nx, 04.30.Db, 04.70.Bw## I Introduction

Gravitational wave astronomy will soon provide astrophysicists with a new tool to observe and analyze some of the most energetic phenomena in our universe. Collisions of compact objects, such as neutron stars, stellar-mass black holes, and super-massive black holes, should produce characteristic gravitational wave signals that are observable to cosmological distances. In particular binary black hole systems are among the most promising sources of gravitational waves for both the current generation of ground-based detectors, such as LIGO Vogt (1992), and for the next generation of space-based detectors, such as LISA LIS ; Danzmann and Rudiger (2003). While early ground-based detectors will need to use templates from theoretically derived waveforms in order to extract the weak signal from the much larger noise, space-based detectors will require accurate models of the gravitational wave sources to extract the characteristic information (e.g., mass, spin) from the detected gravitational waves. Thus there is a critical need for accurate gravitational waveforms from realistic simulations of merging black holes.

Numerical studies of Einstein’s equations play an essential role in studying highly dynamic systems, such as binary black holes, that have little or no symmetry. These are not only computationally very demanding, requiring high-performance massive parallel computers, but also mathematically and numerically very challenging. Despite these difficulties a great deal of progress has been made in the last few years Alcubierre (2004), and it is now possible to follow binary black hole evolutions during the last moments of the merger phase Brandt et al. (2000); Alcubierre et al. (2001a, 2004a) and perhaps even through a complete orbit Brügmann et al. (2004). The calculation of the gravitational radiation emission from merging binary black holes has also been made possible through the use of the Lazarus approach, which bridges numerical relativity and perturbative techniques to extract approximate gravitational waveforms Baker et al. (2001, 2002a, 2002b, 2004); Campanelli (2004). These calculations are in good agreement with recent numerical relativity evolutions Alcubierre et al. (2004a).

At present, one of the most serious limiting factor in numerical relativity has been the inability to extract accurate and reliable results from stable three-dimensional numerical simulations of coalescing binary black hole spacetimes. Considerable resources have been dedicated to finding numerically stable formulations of the Einstein evolution equations Nakamura et al. (1987); Shibata and Nakamura (1995); Friedrich (1996); Frittelli and Reula (1996); Baumgarte and Shapiro (1999); Anderson and York (1999); Alcubierre et al. (2001a); Kidder et al. (2001); Sarbach and Tiglio (2002); Shinkai and Yoneda (2003); Lindblom and Scheel (2003); Bona and Palenzuela (2004). This effort has led to stable evolutions of single black hole spacetimes Alcubierre et al. (2001b); Yo et al. (2002); Alcubierre et al. (2003a); Sperhake et al. (2004). However, there is a growing realization that solving a well-posed formulation of the Einstein equations with second-order accurate finite differencing is not sufficient to produce accurate simulations of binary black hole spacetimes. Unfortunately, black hole simulations require large domains with high resolution near the horizons, and the unigrid second-order accurate codes developed in the past are not sufficient, given the foreseeable constraints on memory and CPU speed, to produce accurate waveforms from merging black holes Miller (2005). The Numerical Relativity community is currently pursuing several different approaches to produce accurate evolutions. Among them are the use of finite elements, spectral methods Pfeiffer et al. (2003); Kidder et al. (2000), adaptive mesh finite-difference codes Imbiriba et al. (2004); Fiske et al. (2005); Sperhake et al. (2005); Schnetter et al. (2004); Pretorius and Lehner (2004), and higher order finite-difference codes Calabrese et al. (2004). In addition, fourth-order accurate algorithms have recently been designed to evolve perturbations of nonrotating Lousto (2005) and rotating black holes Pazos-Avalos and Lousto (2004). Ideally, the next generation of finite-difference codes should combine higher order finite-differencing in combination with mesh refinement techniques. However, considerable progress toward highly accurate evolutions can be achieved using a unigrid higher-order finite-difference code along with a coordinate system, such as ‘Fisheye’ Baker et al. (2002a), that concentrates the gridpoints in the central region containing the black holes. In this paper we present a successful approach to performing fourth-order accurate numerical relativity simulations.

Our simulations are based on a new numerical relativity evolution code framework, LazEv, to solve the full non-linear Einstein’s equations in 3D. The LazEv framework is designed to be modular. The code consists of a ‘Method of Lines’ time integrator along with several Mathematica scripts Brügmann et al. that convert tensorial equations into finite-difference algorithms (of arbitrary order) in C. This modularity allows us to quickly implement new formulations and gauge choices as needed. We have currently implemented the ADM formulation as well as several ‘flavors’ of the BSSN formulation Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1999). In Sec. II we describe the LazEv framework in detail. In Sec. III we enumerate many of the basic techniques utilized in numerical relativity simulations which we have implemented within LazEv, providing the foundation for detailed testing and application in black hole simulations.

In Sec. IV we test our fourth-order techniques on standard numerical relativity testbed problems Alcubierre et al. (2004b) demonstrating fourth-order accuracy. Finally, in Sec V we apply these methods to study head-on binary black hole collisions. Although head-on collisions have limited astrophysical relevance, there has been recent interest in using these spacetimes to develop and test numerical techniques for evolving more generic binary black hole spacetimes Alcubierre et al. (2004c); Fiske et al. (2005); Sperhake et al. (2005). In this paper we demonstrate that our fourth-order techniques can realize significant improvements in gravitational radiation waveform accuracies at typically achievable numerical resolutions.

## Ii The LazEv Evolution Framework

The LazEv evolution framework is a foundation for rapid implementation of numerical relativity evolution system routines, supporting higher order finite differencing techniques with ‘method of lines’ (MoL) integration Gustafsson et al. (1995). Our numerical code uses the Cactus Computational Toolkit cactus_web for parallelization and IO. We have constructed the code so that it is compatible with the initial data and analysis thorns provided with Cactus. Specific evolution routines are produced using Mathematica scripts Brügmann et al. that convert tensorial equations into finite difference algorithms in C. These scripts provide several types of consistency checks that prevent many types of potential coding errors in setting up the tensorial equations. The Mathematica scripts support generic centered spatial finite differencing, with optional upwinded differencing (and other off-center differencing) provided by external macros. Currently we use second and fourth-order spatial finite differencing.

The LazEv framework applies the MoL technique for solving first-order-in-time, hyperbolic PDE’s, in which the PDE’s are converted into coupled ODE’s for all variables at every gridpoint. This is achieved by choosing a discrete numerical grid and finite difference stencils for the spatial derivatives (the finite different stencils couple the fields at neighboring gridpoints), and then integrating the time derivatives of the fields at all gridpoints. This time integration can be carried out by standard ODE integrators.

The LazEv MoL integrator provides a generic framework for integrating hyperbolic PDE’s using Runge-Kutta or ICN style time integrations. The integrator itself has no knowledge of the system that is being evolved. The MoL integrator provides internal timebins pre-ministep, ministep, dissipation, and post-ministep, which are called during each step of the Runge-Kutta or ICN integration (there are several ministeps in each full timestep). The evolution system is chosen by registering routines with MoL to be evaluated during each of these bins. For example, in the BSSN system (see Sec. III.1) we register routines to calculate the time derivative of all the BSSN variables in the ministep timebin, a routine to rescale the BSSN conformal metric to unit determinant in post-ministep, and a routine to subtract off the numerical trace of the BSSN trace-free, conformal extrinsic curvature in post-ministep. In the ADM system, on the other hand, we register routines only in the ministep bin. Additional evolution systems may also register routines during these steps. For example, the gauge evolution systems, which are independent of the main evolution systems, also register routines in ministep and post-ministep.

Currently we have implemented ICN (second-order) and Runge-Kutta (second, third, and fourth-order). In some cases we use Kreiss-Oliger Kreiss and Oliger (1973) dissipation of the form

(1) |

where is one of the evolution equations, is the gridspacing in the th direction, and are the forward and backwards differencing operators (in the th direction), and is the order of the finite differencing used to evaluate .

We use several different styles of finite difference stencils, depending on the order of accuracy desired and the type of derivative considered. Our standard choices for fourth-order accurate evolutions are:

(2) | |||||

(3) | |||||

(4) | |||||

Note that the mixed derivative is obtained by applying the x-derivative and y-derivatives sequentially (the order is irrelevant).

We do not use standard centered differencing for the advection terms (i.e. terms of form ). For these terms we use the following upwinded stencils:

(5) | |||||

(6) | |||||

In addition we use lower-order centered differencing on the planes adjacent to the boundaries. We have two choices for how the lower-order stencils are constructed. We can choose to use second-order centered differencing at all points on these planes and for all directions, or we can use second-order centered differencing only in the direction perpendicular to the boundary. When we use this latter choice we construct mixed second derivatives by applying the fourth-order accurate centered first derivative operator (in a direction tangent to the boundary) to the second-order centered derivative operator (perpendicular to the boundary).

Although we plan to implement excision boundary conditions soon, we currently use the puncture approach (see Ref. Alcubierre et al. (2004c) for a comparison of methods), along with singularity avoiding slicings, to evolve black hole space-times. When using punctures, we found that we needed to either use second-order stencils (but still fourth-order Runge-Kutta time integration) in regions inside the apparent horizons, or use a second-order upwinded stencil for the advection terms over the entire computational domain. The former method produces significantly better waveforms. See Sec. V.3 for further details.

## Iii Formulation

Within our LazEv framework we have implemented support for several standard options in formulating numerical relativity modeling problems. Such a formulation includes selection of an Einstein evolution system and a choice of gauge for the evolving spacetime. In this section we discuss the formulation options presently available within our LazEv framework. Of these we have specific realizations which are appropriate for the fourth-order applications discussed in the subsequent sections.

### iii.1 Evolution Systems

Many alternative formulations of Einstein’s evolution equations have been considered so far (see Ref. Alcubierre et al. (2001c) and references therein). Within the most popular Cauchy or approach, one can currently choose among the first order symmetric hyperbolic formulations of the evolution equations which are mathematically very attractive Kidder et al. (2001); Sarbach and Tiglio (2002); Lindblom and Scheel (2003); Bona and Palenzuela (2004), and various flavors of spatially second-order hyperbolic formulations which are numerically more tractable Nagy et al. (2004); Gundlach and Martin-Garcia (2004). While the former have the advantage of possessing a well-posed continuum limit in the presence of maximally dissipative boundary conditions, it is not yet clear how one can handle dozens of free parameters Kidder et al. (2001); Sarbach and Tiglio (2002) and dynamical gauge choices for the lapse and shift. On the other hand, some second-order hyperbolic formulations have proven to be empirically much more robust than others Shibata and Nakamura (1995); Baumgarte and Shapiro (1999); Pretorius (2005); Kreiss and Ortiz ((2002). In this paper we will focus on the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1999) system, a strongly hyperbolic system Sarbach et al. (2002); Beyer and Sarbach (2004) that has been shown to have some attractive stability properties Alcubierre et al. (2000); Sarbach et al. (2002); Beyer and Sarbach (2004).

The BSSN system is an extension of the standard ADM Arnowitt et al. (1962) system with better numerical properties than the original system. In the ADM system the Einstein Equations are split into six evolution equations for the metric (along with six auxiliary evolution equations for the extrinsic curvature)

(7) | |||||

(8) | |||||

and four constraints equations

(9) | |||||

(10) |

Here is the operator , is the Lie derivative with respect to the shift vector , is the covariant derivative associated with the 3-metric , is the three-dimensional Ricci tensor, the Ricci scalar, and is the trace of .

The BSSN system of equations is obtained from the standard ADM equations by the following substitutions,

(11) | |||||

(12) |

where and . Three additional variables

(13) |

are also introduced. The BSSN variables (, , , , and ) obey the following evolution equations Alcubierre et al. (2003b)

(14) | |||||

(15) | |||||

(16) | |||||

(17) | |||||

(18) | |||||

where indicates that only the trace-free part of the tensor is used and is given by

(19) | |||||

(20) | |||||

and is the covariant derivative with respect to . is replaced by in Eq’s (14) - (20) wherever it is not differentiated. Note that Eq. (18) gives the derivative of rather than the derivative. The Lie derivatives of the non-tensorial quantities (, , and ) are given by

(21) | |||||

(22) | |||||

(23) | |||||

In addition to the evolution equations, the BSSN variables must also obey the following seven differential constraint equations

(24) | |||||

(25) | |||||

(26) |

and the following two algebraic constraint equations

(27) | |||||

(28) |

We monitor the differential constraints but do not enforce them. However, we enforce the algebraic constraints by rescaling the evolved and subtracting off the trace of the evolved at every timestep.

We can evaluate in two ways. We could either calculate and find its trace numerically and subtract off the numerical trace, or we can use the Hamiltonian constraint Eq. (24) to calculate and subtract that from . Our tests showed very little difference between the two approaches even when the Hamiltonian constraint violations were relatively large. Unless otherwise specified we use the latter method.

### iii.2 Gauge Choices

We have implemented three basic types of lapse conditions: (i) maximal slicing, (ii) “Bona-Massó” type lapses Bona et al. (1995); Alcubierre et al. (2001b); Balakrishna et al. (1996); Alcubierre et al. (2003b), (ii) and densitized lapses. The maximal slicing condition has the form

(29) |

and is implemented using a modified version of the Cactus ‘Maximal’ thorn cactus_web in combination with a fourth-order accurate version of Bernd Brügmann’s ‘BAM_Elliptic’ Brandt and Brügmann (1997); cactus_web thorn. This fourth-order accurate version of ‘BAM_Elliptic’ was developed at UTB by Mark Hannam. The “K-driver” lapses have the form

(30) |

where

(31) |

and is some specified function (usually zero). We have also implemented modifications to this general form by replacing with , adding to the right-hand-side of Eq. (30), and multiplying by factors of in Eq. (30), where is the conformal factor of the puncture data. Equation (30) includes harmonic slicing (), 1+log slicing (), and ‘shock-avoiding’ slicing () Alcubierre (2003).

We have implemented ‘Gamma-driver’ shift conditions Alcubierre et al. (2003b) as well as static corotating shifts. The two versions of the ‘Gamma-driver’ shifts that we implemented are

(32) |

and

(33) |

where and are some prescribed functions over space, and is a function of and . Unless otherwise specified, we take .

### iii.3 Boundaries

Thus far we have implemented fairly simple boundary conditions. For fourth-order upwinded stencils we alternatively use centered finite differencing or second-order upwinding at the second point from the boundary, and second-order centered differencing at the first point from the boundary. The boundary points are filled using radiative boundary conditions Alcubierre et al. (2003b). However, when the shift is zero, or sufficiently small, we evolve on the boundary using the zero-shift form of Eq. (14). We have observed that the boundary algorithm has very little effect on the convergence of the component of the waveforms in Sec. V.

## Iv Apples With Apples Tests

We apply the LazEv framework to a standard set of numerical relativity tests known as ‘Apples with Apples’ tests Alcubierre et al. (2004b). These tests consists of evolving spacetimes with topology with the advantage that there are no boundaries. The results from these testbeds indicate that LazEv is stable and fourth-order convergent with our choices of fourth-order stencils.

### iv.1 Gauge wave test

For this test we evolve the metric

(34) |

where

and is the wavelength. The ADM variables for this metric have the form

(35) | |||||

(36) | |||||

(37) | |||||

(38) |

where all remaining ADM variables (including ) are zero. Note that obeys the harmonic slicing condition

(39) |

In our simulations is ‘live’ and we evolve it using the harmonic slicing condition. A two-dimensional test is obtained by the coordinate transformation

producing non-trivial , , , and .

In order to obtain stable runs we needed to include dissipation of the form given in Eq. (1), with dissipation coefficient .

These one and two-dimensional problems were evolved using a three-dimensional grid with periodic boundary conditions in all directions. We chose a wavelength () of 1 and constructed the grid so that it contained points (6 points for ghost-zones) in the non-trivial directions and 9 points in the trivial directions (the stencil requires 3 ghost-zones), where is the gridsize and is an integer.

Our first test is a weak gauge wave with amplitude and resolutions , where , , . Figure 1 shows the norm of the error in versus time (rescaled by a factor of ). Note that fourth-order convergence (as evident by the overlap of the rescaled curves) implies that the error for the case is 256 times smaller than the case. This relationship breaks down near crossing times, though the higher resolution runs remain in the convergence regime longer. Using the criterion that the numerical solution is sufficiently accurate if the norm of the error is smaller than of the amplitude , we find that the fourth-order code produces accurate results to , , and for the , , and runs respectively.

Figure 2 shows the norm of (rescaled by a factor of ). Again, the higher resolution runs remain in the convergence regime longer. Using the criterion that the numerical solution is sufficiently accurate if the norm of the Hamiltonian constraint is smaller than of the amplitude , we find that the fourth-order code produces accurate results to , , and for the , , and runs respectively. The Hamiltonian constrain violation is a stricter measure of the quality of the results than the error in because the Hamiltonian constraint involves second derivatives of the metric (which are harder to model).

Our next test is a stronger gauge wave. Figure 3 shows the norm of versus time (rescaled by a factor of ) for the gauge-wave test with amplitude and resolutions , where , , . Note that in this higher amplitude case the runs remain convergent for 80 crossing times (approximately one-tenth as long as the runs). Using the criterion that the numerical solution is sufficiently accurate if the norm of the Hamiltonian constraint is smaller than of the amplitude , we find that the fourth-order code produces accurate results to , , and for the , , and runs respectively. Hence these runs are accurate for roughly one-tenth the time that the are accurate.

Though we have applied our fully three-dimensional code to the last two test problems, the dynamics in these cases are non-trivial only in one (numerical grid) direction. Our last gauge wave test is non-trivial in two grid directions. Figure 4 shows the norm of versus time (rescaled by a factor of ) for the two-dimensional gauge-wave test with amplitude and resolutions , where , , . These two-dimensional runs are less stable than the corresponding one-dimensional run, crashing before 100 crossing times. Convergence is, nonetheless, maintained for about 55 crossing times, longer for the higher resolution runs. Using the criterion that the numerical solution is sufficiently accurate if the norm of the Hamiltonian constraint is smaller than of the amplitude , we find that the fourth-order code produces accurate results to , , and for the , , and runs respectively. Note that these two-dimensional runs are accurate roughly two-thirds as long as the corresponding one-dimensional runs.

### iv.2 Gowdy wave test

In this section we present results for the ‘Polarized Gowdy Wave’ test. The ‘Polarized Gowdy Wave’ metric is given by

(40) |

where

(41) | |||||

(42) | |||||

We use this metric to obtain initial data and evolve backwards in time using the harmonic slicing condition. The above metric only varies as a function of and but we can obtain a two-dimensional problem by introducing new coordinates , , .

Here again we use a full three-dimensional grid to solve one and two-dimensional problems. We use periodic boundary conditions in all directions, and our grid consists of points (plus 6 for ghost-zones) in the non-trivial directions, and 9 points in the trivial directions. Additionally, we needed to add dissipation of the form given in Eq. (1), with dissipation coefficient , to stabilize the runs.

Figures 5 and 6 show the fourth-order convergence of the Hamiltonian constraint and the component of the metric of the one-dimensional Gowdy Wave test for 1000 crossing times. Note that these convergence plots imply that the errors in and in the run are 256 times smaller than these errors in the run. Unlike the previous gauge-wave test, here the amplitude of the metric functions are damped. So our criterion for an acceptable value of the Hamiltonian constrain violation is time dependent. In this case we will use the criterion that the error norm, divided by the norm of , must be smaller than . At the end of the run the of (the function, not the error) is of order 1. Both the Hamiltonian constraint and the error in are smaller than this quality criterion for all resolutions.

Figure 7 shows the fourth-order convergence of the Hamiltonian constraint for the two-dimensional Gowdy Wave test. The early time lack of convergence is due to roundoff effects, and, although the curve does not lie on the and curves, the overlap of the and curves confirm that the code is fourth-order convergent with sufficient resolution.

Figure 8 shows the fourth-order convergence of for the two-dimensional Gowdy Wave test. Unlike , clearly demonstrates fourth-order convergence at the low () resolution.

We conclude from the gauge wave and Gowdy wave results that our BSSN code is stable, accurate, and convergent.

## V Head-on Binary Black-Hole Collisions

In this section we present results for head-on collisions of two equal-mass Misner-Wheeler-Brill-Lindquist (MWBL) black holes Brill and Lindquist (1963); Misner and Wheeler (1957). These results extend our previous tests of LazEv to more interesting nonlinear spacetimes containing binary black holes.

The MWBL data represent conformally flat slices of multiple black hole space-times with punctures. The data are parametrized by puncture masses and puncture positions (in the conformal space). The ADM mass is given by the sum of , and the ADM linear momentum and angular momentum are zero. The data have the form

(43) | |||||

(44) |

where

(45) |

and .

### v.1 Setup

We used MWBL octant-symmetric data consisting of two equal mass black holes aligned along the z-axis, allowing us to evolve the data using the octant. In addition, we use a ‘Transition Fisheye’ transformation Baker et al. (2002a) to enlarge the physical domain without sacrificing resolution near the punctures. The ‘Transition Fisheye’ transformation is a smooth radial transformation from an inner resolution fixed by the Cactus gridspacing () to an outer resolution that is generally lower than the inner resolution. This transformation is parametrized by an outer ‘de-resolution parameter’ (the effective grid-spacing at the edge of the grid is ), a transition width parameter ‘’, and the center of the transition region ‘’. The transformation has the form

(46) |

where is the physical radius corresponding to the coordinate radius and is given by

(47) |

For these runs we set the mass parameter of the two holes to , and the puncture positions to and . For most of the runs, the computational grid extended to , which corresponds to a physical size of . We also performed some short runs where the computational grid extended to , which corresponds to a physical size of . These latter runs were used to determine the effect of the boundaries on the waveforms, constraint violations, and horizon mass. The Fisheye parameters used in these runs were , , . We performed runs at resolutions of with gridsizes of gridpoints along one octant, where . The gridsizes were chosen so that the punctures lie halfway between gridpoints. For the runs with the boundary at we used the same resolutions but chose gridsizes of gridpoint along one octant, where .

We evolved the data with the standard 1+log lapse, where the initial value of (see Eq. 31) is zero, and a modified 1+log lapse, where the initial value of is

(48) |

We choose 1 for the initial value of the lapse. When using the modified lapse we set , and . This modified lapse, unlike the original 1+log lapse, collapses at the puncture in the continuum limit. Analytically the standard 1+log lapse should retain its initial value at the puncture throughout the entire evolution. Since we start with a lapse of one, the lapse should not collapse at the puncture. However, the lapse will collapse near the puncture, and a lack of resolution drives the lapse to zero at the puncture as well. However, as the grid is refined this artificial collapse at the puncture is delayed, leading to short wavelength features and blow-ups near the puncture. The initial value of in Eq. (48) forces the lapse to collapse near the puncture in the continuum limit. Unfortunately this modification introduces its own problems as described below (see Fig. 21).

We used the second form of the ‘Gamma-Driver’ shift Eq. (III.2) to evolve the shift. The initial value of the shift and its time derivative were zero. We used both a constant () and a spatially varying equal to at infinity and at the punctures, as suggested by Diener Diener . The spatially varying was of the form Alcubierre et al. (2004c)

(49) |

where and are parameters specifying at the puncture and infinity respectively, and is the puncture data conformal factor. The function in the ‘Gamma-Driver’ shift was set to .

Pure fourth-order runs proved to be very unstable near the punctures both with and without upwinded stencils, and Kreiss-Oliger dissipation further destabilized these runs (near the punctures). We had success in stabilizing these runs with two techniques: (i) using second-order upwinding of advection terms over the entire grid, and (ii) a localized (within the apparent horizon) order reduction (LOR) to all second-order accurate spatial differentiation (including second-order upwinding of advection terms). This latter method provides significantly more accurate waveforms. When using the LOR method we found that lower-order accuracy is needed both at the punctures and at the origin. However, we only reduce the order of accuracy near the origin after a common apparent horizon forms (typically after the common apparent horizon forms). We found that Kreiss-Oliger dissipation can be used to remove late time instabilities if the dissipation coefficient is set to zero in the LOR region. However, we did not use dissipation in the runs presented below.

We implement LOR in the following way. First we calculate the time derivatives of all variables using the standard fourth-order stencils in Eq’s (2) - (6). Then we overwrite all points within some given coordinate-ellipsoid with the time derivatives calculated using the standard second-order stencils (with second-order upwinded advection terms). The dissipation operator inside the LOR region can be either the same one chosen for the rest of the grid or a lower-order operator (in either case the dissipation coefficient inside the LOR region can be different from the dissipation coefficient used for the rest of the grid). Only the spatial finite differencing is changed; the time integrator is the same for all gridpoints. The size of the ellipsoid is chosen at runtime and up to eleven ellipsoids may be used. The user is free to choose when each ellipsoid is activated. For example, in the head-on binary black hole collision case, we specify a small ellipsoid for each of the individual apparent horizons and a larger ellipsoid for the common apparent horizon, where the larger ellipsoid is activated only after the common apparent horizons forms. The sizes of these ellipsoids can be determined by evolving with second-order accuracy over the entire grid and finding the sizes of the apparent horizons versus time.

For the head-on-collision runs the LOR regions consisted of spheres of radius centered on the punctures activated at , as well as an ellipsoid with semi-axes activated at (the punctures were located on the z-axis). For comparison, the common apparent horizon has a minimum radius of and a maximum radius of at . The individual apparent horizons are approximately spherical with radii at .

We used radiative boundary conditions for all variables except (which were evolved via on the boundary). Table 1 gives the radiative boundary condition parameters for all variables. In Table 1 ‘’ is the Fisheye de-resolution parameter. We enforce the algebraic constraints Eq’s (27),(28) prior to applying the boundary conditions. Hence, the boundary points do not satisfy these constraints exactly.

Variables | Asymptote | speed |
---|---|---|

0 | ||

0 | ||

0 | ||

0 |

We updated the Zorro thorn of the Lazarus Toolkit Baker et al. (2002a) to compute the Weyl scalars to fourth-order accuracy and to make it compatible with octant, bitant, and quadrant symmetries in any direction as well as rotation symmetry about the z-axis. In addition we added a spherical harmonic decomposition routine to generate the and modes shown below.

### v.2 Second-Order Accurate Results

We first confirmed that our code can reproduce the second-order accurate waveforms published in Alcubierre et al. (2003b). For these runs we used standard 1+log lapse (i.e. ) and Gamma-Driver shift with . We ran with gridsizes of gridpoints and grid resolutions of for . We calculated as well as its and components using Zorro Baker et al. (2002a). Note that the and modes of are purely real for this test.

Figure 9 shows the mode of at for these three resolutions along with the Richardson extrapolation of these data. In addition Fig. 9 shows the differences and . These two difference curves overlap reasonably well indicating that the waveforms are second-order accurate. However, the phase drift between the two curves makes an evaluation of the exact order of convergence difficult. Figure 10 shows the convergence rate for this mode given by

The measured convergence rate oscillates wildly but averages to about 2.

Strictly speaking is not in the radiation zone so given above does not represent the asymptotic waveform. However, the code’s performance in calculating the mode of at is indicative of its performance in calculating this mode at large . Additionally, boundary errors contaminate the modes at large (see Sec. V.3). Thus by extracting at we can maximize the amount of non-trivial, physically correct quasinormal oscillations in the waveforms.

Figure 11 shows the norm of the Hamiltonian constraint for these three runs. Note that the constraints have not been rescaled. From the figure we can see that the constraints tend to get bigger as the resolution is increased. The source of these constraint violations is located between the puncture and origin. High frequency features develop there, leading to strong constraint violations. However, as seen in Fig. 12 these large constraint violations do not leak out of the apparent horizon. Thus, although there are significant problems inside the horizon, the region outside the horizon remains uncontaminated.

Our calculations of the gravitational waveforms from the Newman-Penrose scalar , which contains higher-order spatial derivatives of the metric, do not reveal any additional degree of noise and distortion with respect the waveforms obtained from the Zerilli-Moncrief formalism Alcubierre et al. (2004c); Sperhake et al. (2005).

### v.3 Fourth-Order Accurate Results

Purely fourth-order runs of puncture data proved to be unstable with more than one black hole. We found that we could stabilize the evolution by reducing the spatial discretization to second-order inside the apparent horizons, and we successfully ran the MWBL head-on collision data with fourth-order discretization, fourth-order Runge Kutta time integration, and second-order LOR regions inside the apparent horizons. For these runs we used the spatially varying form of the Gamma-driver shift given by Eq. (III.2) and Eq. (49). Figure 13 demonstrates the convergence of the component of at . Note that the system is not strictly fourth-order accurate. There appears to be an additional phase discrepancy between the two differences. However, the amplitude of the differences appears to be falling at a rate consistent with better than fourth-order accuracy. The run crashed at due to an instability near the origin (see Fig. 19). This unstable mode was triggered by advection terms (the collapse of the lapse near the origin ensures that these are the only terms which can lead to a blowup). The fields which blew up most strongly were the . This blow-up does not appear to be related to the quadratic blow-up in at the puncture discussed later. That term led to a blow up proportional to (in ) near the puncture, which masked the unstable behavior near the origin. Both blow-ups can be modified by various choices in the gauge conditions and it is likely that a better choice of gauge will lead to longer fourth-order evolutions.

A fit of the mode of to the quasinormal form gave an exponential damping factor of and frequency of . The exponential damping factor agrees to within 6% of the expected value of , and the frequency agrees to within of the expected value of Nollert and Schmidt (1992).

In Fig. 14 we show the Richardson extrapolated value of the mode of , calculated using the second-order accurate results, along with the values obtained from second and fourth-order evolutions with LOR. Note that the medium resolution fourth-order run outperforms the high resolution second-order run. Figure 15 shows a magnified view of Fig. 14 (note that the legends in the two figures are different).

Figure 16 shows the convergence of the mode of at calculated using fourth-order accuracy with LOR. The plot shows the differences and , as well the mode itself for . Note that the latter rescaled difference is smaller than the former. This indicates that the component of the waveform is converging faster than fourth-order. In addition, note that the error in the waveform is comparable to the amplitude of the waveform.

Figure 17 shows a magnified view of the mode of at for the , , and fourth-order runs as well as the second-order run. Note that the second-order results are again inferior to the fourth-order (with LOR) results.

The modes, unlike the modes, contain significant contamination from the boundary at late times. For example, when extracting at we find that the amplitude of the mode of (which is non-zero due to boundary errors) is approximately that of the mode at . This contamination is stronger when the extraction sphere is further out. However, the waveform extracted at shows significantly less contamination. The reason for this is that the incoming error is delayed by about while the outgoing mode is advanced by (for a net gain of in reliable data). In addition the outgoing mode, which has a falloff, is correspondingly larger at this smaller radius. This boundary contamination problem can be mitigated by pushing the outer boundary further out. This can be achieved within the existing LazEv framework using a stronger fisheye transformation or adding more gridpoints.

Although the waveforms from the fourth-order runs converge as expected the Hamiltonian constraint does not. High frequency features near the puncture and origin lead to large Hamiltonian constraint violation. At the points surrounding the puncture these high frequency features induce a quadratic blow-up in for the run. This quadratic behavior, which is confined to the nearest neighbor points to the puncture, leads to a blow-up of in the Hamiltonian constraint due to terms proportional to . This blow-up in is entirely localized to the puncture and does not affect points outside. In the following figures we demonstrate that the large constraint violations inside the apparent horizon do not leak out and that the constraint violations outside the apparent horizon converge to zero when boundary effects are removed.

Figure 18 shows the Hamiltonian constraint along the z-axis (the points surrounding the puncture have been removed) at for the second and fourth-order runs for . The Hamiltonian constraint inside the horizon is as much as times larger for the fourth-order run. However, outside the horizon the constraint violations are very similar. Note that we do not expect the constraints to converge to zero in this case because our boundary conditions are not constraint preserving and boundary constraint violations have contaminated the solution at this time (the boundaries were at in physical coordinates). The fourth-order run with LOR crashed due to an instability near the origin. Figure 19 shows the unstable mode in .

In Fig. 20 we show the norms of the Hamiltonian constraint , momentum constraint , and the BSSN constraints (note that the and components of the both momentum and BSSN constraints are equal). We restricted these norms to the region outside (the horizon is at ) and inside , where is the physical radius. The outer boundaries were located at in physical coordinates. These runs required double the number of gridpoints (along each direction) as the standard runs for the same resolutions. The reasonable overlap of the and rescaled curves for of evolution indicates that the norms are approximately fourth-order convergent (the convergence is lost near due to constraint violating modes propagating in from the boundary). In addition, the amplitude of the constraint violations is , which is nine order of magnitude smaller than the maximum Hamiltonian violation inside the horizon. We conclude from this figure that the second-order errors introduced by the LOR differencing, as well as the extremely large Hamiltonian constraint violation inside the apparent horizon, do not propagate outside of the apparent horizon.

The extreme Hamiltonian violation near the puncture can be removed using the modified 1+log lapse. However, this modification tends to destabilize the runs at late times. Figure 21 shows the norm of the Hamiltonian constraint violation for the fourth-order run with LOR. The modified lapse produces a constraint violation smaller than 1 for but also causes the run to crash at . The standard 1+log lapse produces a “stable” evolution to (see above) but has a constraint violation proportional to at the puncture. If short term evolutions, like the ones required for the Lazarus techniques, are required, then the modified 1+log is preferred.

Figure 22 shows the horizon mass versus time, as calculated by the AHFinderDirect thorn Thornburg (2004), for the and runs with LOR. The plot also shows the horizon mass of a run with double the standard number of gridpoints in each direction and outer boundary at in physical coordinates. The relatively sharp increase in the horizon mass at is due to contamination from the boundary. When the boundary is moved out to (from ) this increase is delayed by roughly . Also note that the slope of the increase is reduced when the boundary is further out. The measured convergence rate of the horizon mass (based on runs with the boundary at ) exceeded 3.8 for of evolution.

The horizon-mass versus time plot shows that the late time boundary contamination contains significant erroneous gravitational radiation. To determine when the calculated gravitational waveforms are good approximations to the true waveforms, we need to confirm that (i) the waveforms converge and do not change significantly when the resolution is further increased, (ii) the waveforms do not change significantly when the boundary effects are removed (which we achieved by causally disconnecting the boundaries from the observer), and (iii) the Hamiltonian constraint violation (inside the extraction region) is much smaller than the waveform amplitude. We have shown that the waveforms converge, and that the Hamiltonian constraint violation converges in the extraction region (see Fig. 20). Here we show the effects of the boundary on the Hamiltonian constraint and the waveforms. Figure 23 shows Hamiltonian constraint violation along the axis at time for the standard run and a run with twice the standard number of gridpoints. The boundary in this latter case was at (note that the -axis plots the fisheye coordinate and that corresponds to the physical coordinate and corresponds to ). The figure shows the Hamiltonian constraint outside the apparent horizon. Note that at this time the boundary errors for the larger run have just reached . The Hamiltonian constraint, in the range to , is 100 times smaller for the run with the larger boundary. Even at , when the boundary errors have contaminated the entire grid, the Hamiltonian constraint is 20 times smaller for the larger run in this range.

As seen in Fig. 24, the effect of reducing the boundary noise is not readily apparent in the mode of . However, as seen in Fig. 25, the effect is readily observable in the mode. We conclude, based on these two figures, that the mode is represented accurately to when the boundaries are located at , while the mode is represented accurately only to for the same boundary location. The mode is more accurate because the boundary contamination contains relatively high frequencies which are filtered more effectively by the angular integration. Note that we placed the observer at and that the waveforms become inaccurate sooner when the extraction sphere is placed at larger radii.