Vacuum and black hole evolutions

6.2 Vacuum and black hole evolutions

If one wants to simulate the most interesting astrophysical sources of gravitational radiation, one must have a code able to follow, in a stable manner, gravitational waves themselves on a background spacetime. It has been observed by all numerical relativity groups that the stability of a numerical code, which solves Einstein’s field equations, not only depends on the numerical algorithm, but also on the particular formulation of the equations. Successes in the simulations of binary systems of compact objects in general relativity (see Section 6.3) are also due to numerous studies and advances in the formulations of Einstein’s equations. The methods known at present that work for the numerical evolution of binaries are the generalized harmonic coordinate method [82 , 86 , 177

] and the BSSN method (for Baumgarte–Shapiro–Shibata–Nakamura [25

, 195

]). In particular, these two formulations of the field equations have the important property that constraint violating modes (see discussion following) stay at a reasonable level during the evolution. Use of the generalized harmonic gauge requires constraint damping terms; and a particular method suited for harmonic evolution, which was proposed by Gundlach et al. [111 ], enabled Pretorius to evolve black hole spacetimes [177

, 176

It is, therefore, a crucial step to devise such a stable formulation, and more particularly with spectral methods, because they add very little numerical dissipation and thus, even the smallest instability is not dissipated away and can grow to an unacceptable level. The situation becomes even more complicated with the setup of an artificial numerical boundary at a finite distance from the source, needing appropriate boundary conditions to control the physical wave content, and possibly to limit the growth of unstable modes. All these points have been extensively studied since 2000 by the Caltech/Cornell groups and their pseudospectral collocation code SpEC [125 , 127 , 188 , 187 , 138 , 120 , 126 , 137 , 49 ]; they were followed in 2004 by the Meudon group [37 ] and in 2006 by Tichy [219 ].

Next, it is necessary to be able to evolve black holes. Successful simulations of black hole binaries have been performed using the black-hole puncture technique [55 , 18 ]. Unfortunately, the dynamic part of Einstein fields are not regular at the puncture points and it seems difficult to regularize them so as to avoid any Gibbs-like phenomenon using spectral methods. Therefore, punctures are not generally used for spectral implementations; instead the excision technique is employed, removing part of the coordinate space inside the apparent horizon. There is no need for boundary conditions on this new artificial boundary, provided that one uses a free-evolution scheme [188 ], solving only hyperbolic equations. In the considered scheme, and for hydrodynamic equations as well, one does not need to impose any boundary condition, nor do any special treatment on the excision boundary, contrary to finite difference techniques, where one must construct special one-sided differencing stencils. On the other hand, with a constrained scheme, elliptic-type equations are to be solved [37 ] and, as for initial data (see Sections 5.3 and 5.6), boundary conditions must be provided, e.g., on the apparent horizon, from the dynamic horizon formalism [105 ].

Finally, good outer boundary conditions, which are at the same time mathematically well posed, consistent with the constraints and prevent as much as possible reflections of outgoing waves, must be devised. In that respect, quite complete boundary conditions have been obtained by Buchman and Sarbach [53 ].

6.2.1 Formulation and boundary conditions

Several formulations have been proposed in the literature for the numerical solution of Einstein’s equations using spectral methods. The standard one is the 3+1 (a.k.a. Arnowitt–Deser–Misner – ADM) formalism of general relativity [14 , 229 ] (for a comprehensive introduction, see the lecture notes by Gourgoulhon [97 ]), which has been reformulated into the BSSN [25 , 195 ] for better stability. But first, let us mention an alternative characteristic approach based on expanding null hypersurfaces foliated by metric two-spheres developed by Bartnik [20 ]. This formalism allows for a simple analysis of the characteristic structure of the equations and uses the standard “edth” ( $∂$ ) operator on $2 S$ to express angular derivatives. Therefore, Bartnik and Norton [21 ] use spin-weighted spherical harmonics (see Section 3.2.2) to numerically describe metric fields.

Coming back to the 3+1 formalism, Einstein’s equations split into two subsets of equations. First, the dynamic equations specifying the way the gravitational field evolves from one timeslice to the next; then, the constraint equations, which must be fulfilled on each timeslice. Still, it is well known that for the Einstein system, as well as for Maxwell’s equations of electromagnetism, if the constraints are verified on the initial timeslice, then the dynamic equations guarantee that they shall be verified in the future of that timeslice. Unfortunately, when numerically doing such free evolution, i.e. solving only for the dynamic equations, small violations of the constraints due to round-off errors appear to grow exponentially (for an illustration with spectral methods, see, e.g., [188 , 219 ]). The opposite strategy is to discard some of the evolution equations, keeping the equations for the two physical dynamic degrees of freedom of the gravitational field, and to solve for the four constraint equations: this is a constrained evolution [37 ].

The advantages of the free evolution schemes are that they usually allow one to write Einstein’s equations in the form of a strongly- or symmetric-hyperbolic system, for which there are many mathematical theorems of existence and well-posedness. In addition, it is possible to analyze such systems in terms of characteristics, which can give very clear and easy-to-implement boundary conditions [126 ]. Using finite-difference numerical methods, it is also very CPU-time consuming to solve for constraint equations, which are elliptic in type, but this is not the case with spectral methods. On the other hand, constrained evolution schemes have, by definition, the advantage of not being subject to constraint-violation modes. Besides, the equations describing stationary spacetimes are usually elliptic and are naturally recovered when taking the steady-state limit of such schemes. Finally, elliptic PDEs usually do not exhibit instabilities and are known to be well posed. To be more precise, constrained evolution using spectral methods has been implemented by the Meudon group [37 ], within the framework of the BSSN formulation. Free-evolution schemes have been used by Tichy [219 ] (with the BSSN formulation) and by the Caltech/Cornell group, which has developed their Kidder–Scheel–Teukolsky (KST) scheme [127 ] and have later used the Generalized-Harmonic (GH) scheme [137 ]. The KST scheme is, in fact, a 12-parameter family of hyperbolic formulations of Einstein’s equations, which can be fine tuned in order to stabilize the evolution of, e.g., black hole spacetimes [188 ].

Even when doing so, constraint-violating modes grow exponentially and three ways of controlling their growth have been studied by the Caltech/Cornell group. First, the addition of multiples of the constraints to the evolution system in order to minimize this growth. The parameters linked with these additional terms are then adjusted to control the evolution of the constraint norm. This generalized version of the dynamic constraint control method used by Tiglio et al. [221 ] has been presented by Lindblom et al. [138 ] and tested on a particular representation of the Maxwell equations. Second, Lindblom et al. [138 ] devised constraint preserving boundary conditions from those of Calabrese et al. [54 ], where the idea was to get maximally dissipative boundary conditions on the constraint evolution equations [138 , 126 ]. This second option appeared to be more efficient, but still did not completely eliminate the instabilities. Finally, bulk constraint violations cannot be controlled by constraint-preserving boundary conditions alone, so Holst et al. [120 ] derived techniques to project at each timestep the solution of the dynamic equations onto the constraint submanifold of solutions. This method necessitates the solution of a covariant inhomogeneous Helmholtz equation to determine the optimal projection. Nevertheless, the most efficient technique seems to be the use of the GH formulation, which also incorporates multiples of the constraints, thus exponentially suppressing bulk constraint violation, together with constraint-preserving boundary conditions [137 ].

Boundary conditions are not only important for the control of the constraint-violation modes in free evolutions. Because they cannot be imposed at spatial infinity (see Section 3.1.2), they must be completely transparent to gravitational waves and prevent any physical wave from entering the computational domain. A first study of interest for numerical relativity was done by Novak and Bonazzola [158 ], in which gravitational waves are considered in the wave zone, as perturbations of flat spacetime. The specificity of gravitational waves is that they start at the quadrupole level ( $ℓ = 2$ ) in terms of spherical harmonics expansion. Standard radiative boundary conditions (known as Sommerfeld boundary conditions [202 ]) being accurate only for the $ℓ = 0$ component, a generalization of these boundary conditions has been done to include quadrupolar terms [158 ]. They strongly rely on the spectral decomposition of the radiative field in terms of spherical harmonics and on spherical coordinates. More specific boundary conditions for the Einstein system, in order to control the influx of the radiative part of the Weyl tensor, have been devised by Kidder et al. [126 ] for the KST formulation, generalizing earlier work by Stewart [208 ] and Calabrese et al. [54 ]. They were then adapted to the GH formulation by Lindblom et al. [137 ]. Rinne [180 ] has studied the well-posedness of the initial-boundary–value problem of the GH formulation of Einstein’s equations. He has considered first-order boundary conditions, which essentially control the incoming gravitational radiation through the incoming fields of the Weyl tensor. He has also tested the stability of the whole system with a collocation pseudospectral code simulating a Minkowski or Schwarzschild perturbed spacetime. Generalizing previous works, a hierarchy of absorbing boundary conditions has been introduced by Buchman and Sarbach [53 ], which have then been implemented in the Caltech/Cornell SpEC code by Ruiz et al. [182 ], together with new sets of absorbing and constraint-preserving conditions in the generalized harmonic gauge. Ruiz et al. have shown that their second-order conditions can control the incoming gravitational radiation, up to a certain point. In addition, they have analyzed the well-posedness of the initial-boundary–value problems arising from these boundary conditions. Rinne et al. [181 ] have compared various methods for treating outer boundary conditions. They have used the SpEC code to estimate the reflections caused by the artificial boundary, as well as the constraint violation it can generate.

6.2.2 Gauges and wave evolution

The final ingredient before performing a numerical simulation of the dynamic Einstein system is the gauge choice. For example, the analytical study of the linearized gravitational wave in vacuum has been done with the harmonic gauge, for which the coordinates $μ {x }$ verify the scalar covariant wave equation

This is the definition of the form $H μ$ , where $gμν$ is the metric and $∇ σ$ the associated covariant derivative. Recent work by the Caltech/Cornell group uses the GH formulation in which the gauge choice is achieved through the specification of $H μ$ as an arbitrary function of ${x μ}$ and $g μν$ , which can be set, for instance, to its initial value [189

]. Still, it is with the KST formulation, and with the lapse and shift set from the analytic values, that Boyle et al. [49

] have submitted the Caltech/Cornell SpEC code to the “Mexico City tests” [1 ]. These are a series of basic numerical relativity code tests to verify their accuracy and stability, including small amplitude linear plane wave, gauge wave and Gowdy spacetime evolutions. These tests have been passed by the Caltech-Cornell code using a Fourier basis for all three Cartesian coordinates and a fourth-order Runge-Kutta timestepping scheme. In the particular case of the linear plane wave, Boyle et al. [49 ] exhibit the proper error behavior, which increases as the square of the wave amplitude, because all nonlinear terms are neglected in this test. The authors have also shown that the use of filtering of the spherical harmonics coefficients was very effective in reducing nonlinear aliasing instabilities. Gauge drivers for the GH formulation of Einstein’s equations have been devised by Lindblom et al. [136 ]. They provide an elegant way of imposing gauge conditions that preserve hyperbolicity for many standard gauge conditions. These drivers have been tested with the SpEC code.

Within the constrained formulation of Einstein’s equations, the Meudon group has introduced a generalization of the Dirac gauge to any type of spatial coordinates [37 ]. Considering the conformal 3+1 decomposition of Einstein’s equations, the Dirac gauge requires that the conformal three-metric $&tidle;γij$ (such that $ij det &tidle;γ = 1$ ) be divergence-free with respect to the flat three-metric (defined as the asymptotic structure of the three-metric and with the associated covariant derivative $¯D$ )

The time coordinate is set by the standard maximal slicing condition. These conditions turn out to be dynamic gauge conditions: the lapse and the shift are determined through the solution of elliptic PDEs at each timestep. With this setting, Bonazzola et al. have studied the propagation of a three-dimensional gravitational wave, i.e. the solution of the fully nonlinear Einstein equations in vacuum. Their multidomain spectral code based on the Lorene library [99

] was able to follow the wave using spherical coordinates, including the (coordinate) singular origin, and to let it out of the numerical grid with transparent boundary conditions [158 ]. Evolution was performed with a second-order semi-implicit Crank–Nicolson time scheme, and the elliptic system of constraint equations was solved iteratively. Since only two evolution equations were solved (out of six), the others were used as error indicators and proved the awaited second-order time convergence. A preliminary analysis of the mathematical structure of the evolution part of this formalism done by Cordero et al. [64 ] has shown that the choice of Dirac’s gauge for the spatial coordinates guarantees the strongly hyperbolic character of that system as a system of conservation laws.

6.2.3 Black hole spacetimes

As stated at the beginning of Section 6.2, the detailed strategy to perform numerical simulations of black hole spacetimes depends on the chosen formulation. With the characteristic approach, Bartnik and Norton [21 ] modeled gravitational waves propagating on a black hole spacetime in spherical coordinates, but with a null coordinate $z = t − r$ . Interestingly, they combined a spectral decomposition on spin-weighted spherical harmonics for the angular coordinates and an eighth-order scheme using spline convolution to calculate derivatives in the $r$ or $z$ direction. Integration in these directions was done with a fourth or eighth-order Runge–Kutta method. For the spectral part, they had to use Orszag’s 2/3 rule [56 ] for antialiasing. This code achieved a global accuracy of 10^–5 and was able to evolve the black hole spacetime up to $z = 55M$ . More recently, Tichy has evolved a Schwarzschild black hole in Kerr–Schild coordinates in the BSSN formulation, up to $t ≃ 100M$ [219 ]. He used spherical coordinates in a shell-like domain, excising the interior of the black hole. The expansion functions are Chebyshev polynomials for the radial direction, and Fourier series for the angular ones.

Most successful simulations in this domain have been performed by the Caltech/Cornell group, who seem to be able to stably evolve forever not only a Schwarzschild, but also a Kerr black hole perturbed by a gravitational wave pulse [137 ], using their GH formulation with constraint-damping and constraint-preserving boundary conditions. However, several attempts had been reported by this group before, starting with the spherically-symmetric evolution of a Schwarzschild black hole by Kidder et al. [128 ]. Problems had arisen when trying three-dimensional simulations of such physical systems with the new parameterized KST formalism [127 ]. Using spherical coordinates in a shell-like domain, the authors decomposed the fields (or Cartesian components for tensor fields) on a Chebyshev radial base and scalar spherical harmonics. The integration in time was done using a fourth-order Runge–Kutta scheme and the gauge variables were assumed to keep their analytical initial values. The evolution was limited by the growth of constraint-violating modes at $t ∼ 1000M$ . With a fine-tuning of the parameters of the KST scheme, Scheel et al. [188 ] have been able to extend the lifetime of the numerical simulations to about $8000M$ . On the other hand, when studying the behavior of a dynamic scalar field on a fixed Kerr background, Scheel et al. [187 ] managed to get nice results on the late time decay of this scalar field. They had to eliminate high-frequency numerical instabilities, with a filter on the spherical harmonics basis, following again Orszag’s 2/3 rule [56 ] and truncating the last third of the coefficients. It is interesting to note that no filtering was necessary on the radial (Chebyshev) basis functions. A more complicated filtering rule has been applied by Kidder et al. [126 ] when trying to limit the growth of constraint-violation in three-dimensional numerical evolutions of black hole spacetimes with appropriate boundary conditions. They have set to zero the spherical harmonics terms with $ℓ ≥ ℓmax − 3$ in the tensor spherical harmonics expansion of the dynamic fields. The stable evolutions reported by Lindblom et al. [137 ], thus, might be due to the following three ingredients:

GH formalism, exponentially suppressing all small short-wavelength constraint violations,
constraint-preserving boundary conditions,
the filtering of spherical harmonics spectral coefficients.