The non-linear perturbation of a black hole by gravitational waves. I. The Bondi-Sachs mass loss

With the work by Bondi [6], Sachs [31], Newman and Penrose [23, 24, 25] in the 1960’s it has become clear that Einstein’s theory admits gravitational waves. It has also become apparent that, just as in Maxwell’s theory, the notion of radiation is well defined only at infinitely large distances from a source. Penrose’s proposal for the geometric treatment of asymptotically flat space-times [26] demonstrated that the notion of infinity can be precisely captured by the idea of conformal compactification which stipulates that asymptotically flat space-times are those for which one can conformally rescale the metric in such a way that one can attach a conformal boundary, commonly called null infinity and denoted by $\mathscr{I}$ (for more details see [9, 18, 35]). The mass loss formula proved by Bondi et.al. [6, 31], Newman and Penrose [23, 25] is a balance law between asymptotic quantities, i.e., it holds at null infinity $\mathscr{I}$, relating the decrease of energy-momentum contained inside the space-time at a certain instant over time to the gravitational flux radiated away from the source.

In the new era of gravitational physics which started with the detection of gravitational waves [1] the fact that wave forms can be computed accurately only at infinity poses a problem, at least in principle. Current commonly employed frameworks for numerically solving the Einstein equations operate on a finite part of the space-time and thus approximation and limiting procedures must be used (see [5] for an overview.) Nevertheless, these methods have been very successful in practice, in particular for calculating outgoing gravitational waveforms from compact binaries for comparison with observational data from LIGO.

There do exist methods by which the asymptotic quantities can be computed directly without the need for approximation or limiting procedures. They all make use, to varying degrees, of compactification. One of the earliest methods [19, 20, 5] is based on characteristic hypersurfaces which are compactified in order to bring the points at infinity to a finite location. Another method, recently increasingly popular, is based on the hyperboloidal initial value problem in which the space-time is foliated by space-like hypersurfaces which extend to null infinity. These hyperboloidal hypersurfaces are spatially compactified [30]. The method of Cauchy-Characteristic matching [4] tries to combine the characteristic method with a Cauchy evolution of the Einstein equations with a matching procedure across a time-like interface. The first two methods have been implemented successfully but have certain drawbacks. The characteristic approach suffers from the fact that the null hypersurfaces may develop caustics in strong fields while the equations used in the hyperboloidal method are singular on $\mathscr{I}$.

Friedrich’s conformal field equations [13, 12, 14] are a system of equations which generalise Einstein’s equations by implementing Penrose’s proposal of conformal compactification [26]. In this way, the conformal field equations address an asymptotically simple space-time including its conformal boundary (and beyond). Recently, an Initial Boundary Value Problem (IBVP) framework for the Generalized Conformal Field Equations (GCFE) [14] was presented and verified to be well-posed (at least numerically) for a variety of different initial and boundary conditions [3]. In particular, $\mathscr{I}^{+}$ was successfully incorporated in the computational domain for the numerical evolution of a Schwarzschild space-time perturbed by gravitational waves. To do so, the equations were written completely in the space-spinor formalism [33] and Newman and Penrose’s $\eth$-calculus [27, 28] for spin-weighted spherical harmonics on the sphere was employed to allow for fast and accurate pseudo-spectral methods to be utilized [21]. A method for prescribing boundary conditions was put forward which guarantees constraint propagation and was numerically verified. It was shown that in axisymmetry and using a simple ${}_{2}Y_{20}$ mode for the ingoing gravitational radiation, the numerical evolution was stable up to and beyond $\mathscr{I}^{+}$.

Since the computational domain contains at least a portion of the conformal boundary, one can compute the global energy-momentum quantities there which is what we will focus on in this work. Eventually, each time-slice of this numerical evolution will intersect null infinity in a cut. On each cut, one can calculate the asymptotic quantities such as the components of the Bondi-Sachs momentum 4-covector. However there is a major snag: In the literature (e.g. [28]) judicial choices of the coordinates, frame and conformal factor are made to simplify the mathematical analysis as much as possible. These choices, which are generally made based on geometric considerations, are in general inconsistent with what is required to implement a stable numerical scheme. So then, how does one who has run a numerical evolution that includes $\mathscr{I}$ in the computational domain actually compute the correct quantities? In a recent paper [10] we discuss this question and derive expressions for the Bondi-Sachs energy-momentum and the gravitational flux which are valid on arbitrary cuts and for any choice of coordinates, frame and conformal factor.

The aim of this paper is to apply these formulae in a concrete setting and to numerically test them for non-linear gravitational perturbations of the Schwarzschild space-time using the IBVP framework for the GCFE as the numerical evolution scheme. Since this scenario was described in detail in [3] we refer the reader there for thorough checks of correctness, such as constraint convergence tests etc.

The layout of the paper is as follows: Section 2 provides a summary the IBVP framework for the GCFE whilst in Section 4 we describe the specifics of the system to be evolved. Section 3 describes how we make use of the formula for the Bondi components given in [10] and Section 7 presents the numerical results. The paper is concluded with a brief discussion and a summary in Section 8. We use the conventions of [3] throughout.

Here we briefly summarize the most relevant aspects of the GCFE and the IBVP framework for them and point the reader to [15] and [3] respectively for a more comprehensive review.

Given a space-time $(\widetilde{M},\widetilde{g})$, thought of as the “physical” space-time, Einstein’s vacuum field equation with vanishing cosmological constant is given by $\widetilde{R}_{ab}=0$. We are interested in the properties of solutions of this equation “at infinity”. In order to access the points at infinity one introduces an appropriate conformal factor $\Theta$ and defines the “conformal” metric as $g=\Theta^{2}\widetilde{g}$ and the conformal boundary is then defined as the set of points for which $\Theta=0$ and $\mathrm{d}\Theta\neq 0$. This condition defines a regular hypersurface $\mathscr{I}$ in the conformal space-time $(M,g_{ab})$ which can be shown to be null in the case of vanishing cosmological constant.

The metric conformal field equations are a regular extension of the vacuum equations incorporating the conformal rescaling of the physical metric. They express the physical content of the Einstein equations entirely in terms of the geometry on the conformal manifold. The relevant quantities are $\{g_{ab},\Theta,s,P_{ab},K_{abc}{}^{d}\}$ where $s:=\frac{1}{4}\nabla^{c}\nabla_{c}\Theta-\frac{1}{24}R\Theta$, $P_{ab}$ is the Schouten tensor and $K_{abc}{}^{d}:=\Theta^{-1}C_{abc}{}^{d}$ is the gravitational field tensor. In this form of the conformal field equations the conformal factor $\Theta$ is taken to be an unknown of the system, and hence the location of the conformal boundary is not known before a numerical evolution is performed, it must be determined during the evolution.

The GCFE arise by utilizing the full conformal freedom available. This is accomplished by freeing the connection from being compatible with a metric and instead demanding only that it be compatible with the conformal class of the physical metric i.e., by using a Weyl connection $\widehat{\nabla}$. This allows for the use of conformal geodesics, curves which are defined with respect to a Weyl connection in a similar way to how geodesics are defined in terms of a metric connection. The conformal Gauß gauge can be introduced in analogy to the standard Gauß gauge but being adapted to time-like conformal geodesics instead of time-like geodesics. This gives rise to an additional term in the geodesic deviation equation, which helps avoiding the development of caustics, and one can in fact obtain a semi-global covering of the Schwarzschild space-time in this gauge [16].

A startling consequence of exploiting the full conformal freedom through introducing a Weyl connection and the conformal Gauß gauge is that the conformal factor can now be explicitly determined from initial data alone. This means that the location of $\mathscr{I}$ is known before the evolution proceeds. In addition, the evolution equations simplify significantly which has important consequences for the design of boundary conditions. We will now briefly describe the system of equations and the variables we use for the problem at hand.

In complete analogy to other approaches in numerical relativity we perform a $3+1$ decomposition of the space-time and the variables. The difference in our approach compared to others is due to the use of spinors to decompose the variables into irreducible parts. This is motivated mainly by the fact that the Bianchi equation for the rescaled Weyl spinor takes a very efficient form when written in the spinor formalism compared to a tensorial formulation. We use the space-spinor formalism which implements the $3+1$ decomposition in a straightforward way. Most tensorial quantities are written in a spinorial representation with respect to a spin-frame $(o^{A},\iota^{A})$ and its complex conjugate. The resulting $SL(2,\mathbb{C})$ spinors are then transformed to $SU(2)$ spinors i.e., primed indices are transformed to unprimed ones using the Hermitian form $t_{AA^{\prime}}$ (normalised by $t_{AA^{\prime}}t^{AA^{\prime}}=2$) on spin space defined by the time-like tangent vector to the congruence of time-like conformal geodesics used for the Gauß gauge. Finally, the resulting unprimed spinors are then decomposed into irreducible parts.

In a similar way, we split the covariant derivative $\nabla_{a}$ at every point of the space-time manifold into a part which is parallel to $t^{a}=t^{AA^{\prime}}$, and a part orthogonal to it. We define $\partial$ as the covariant derivative operator in the tangential direction which annihilates $t^{a}$ and $\partial_{AB}=\partial_{(AB)}$ is a covariant derivative in the orthogonal direction, also annihilating $t^{a}$, so that

We define frame components $c^{\mu}_{AB}$ as the action of $\partial_{AB}$ on the coordinates

and the $\gamma_{ABC}$ as the action of $\partial_{AB}$ on the spin-frame

where the analogous quantities defined with the $\partial$ are fixed as $c^{\mu}:=\partial x^{\mu}=\sqrt{2}\delta^{\mu}{}_{0}$ and $\gamma_{A}:=\partial o_{A}=0$.

Written in the form of [3] the GCFE are given as a system of differential equations for the unknowns

where $c^{\mu}_{AB}$ are frame components, $\gamma_{ABC}$ are connection coefficients, $K_{ABCD}$ can be thought of as the extrinsic curvature of hypersurfaces of constant time, $f_{AB}$ is defined by the action of the Weyl connection $\widehat{\nabla}$ on the conformal metric $g$ as $\widehat{\nabla}_{c}g_{ab}=-2f_{c}g_{ab}$, $P_{ABCD}$ is the Schouten spinor with respect to the Weyl connection and $\psi_{ABCD}:=\Theta^{-1}\Psi_{ABCD}$ is the rescaled Weyl spinor, which describes the gravitational field.

These are evolved with the symmetric hyperbolic evolution system

Here, $\hat{\psi}_{ABCD}$ denotes the Hermitian conjugate of $\psi_{ABCD}$ defined in the space-spinor formalism by complex conjugation followed by transvection with appropriately many $t^{A^{\prime}}{}_{A}$. It can be seen from this system that there are no equations for the conformal factor $\Theta$ and the 1-form $h_{a}$. In fact, the conformal factor and the components of $h_{a}$ can be obtained explicitly from the Gauß gauge

where the underlined quantities are computed solely with initial data at $t=0$. The corresponding constraint equations are given in the appendix.

As indicated above this system of evolution equations shows a remarkable structure: all equations except for (5g) are simply transport equations along the time-like conformal geodesics. Only the equation for the gravitational field features a derivative in spatial directions. In fact, it is easy to show that this subsystem for the rescaled Weyl spinor is symmetric hyperbolic. This fact is mirrored in the constraint propagation system (see [3] for these equations) in which the constraints for $\psi_{ABCD}$ are propagated with a symmetric hyperbolic system of equations as well, while all the other constraint quantities are simply advected along the time-like conformal geodesics.

Having defined the system, we need to take care of imposing boundary conditions for the components of $\psi_{ABCD}$. This has to be dealt with very carefully, and we follow the maximally dissipative approach as in [17], which was successfully numerically implemented also in [11]. This approach looks at how an “energy” changes over time, and the boundary conditions are chosen so that no energy propagates in from the boundaries. In short, using that the evolution system for $\psi_{ABCD}$ is symmetric hyperbolic we can diagonalize it at every boundary point and find the five distinct characteristic variables $\{\widetilde{\psi}_{0},\widetilde{\psi}_{1},\widetilde{\psi}_{2},\widetilde{\psi}_{3},\widetilde{\psi}_{4}\}$ with speeds $\{-\alpha,-\beta,0,\beta,\alpha\}$, $\alpha>\beta>0$ perpendicular to the boundary. These speeds correspond to $\widetilde{\psi}_{0}$ and $\widetilde{\psi}_{4}$ travelling along light-cones and $\widetilde{\psi}_{1}$ and $\widetilde{\psi}_{3}$ travelling in a time-like fashion. The diagonalization procedure yields a transformation between the $\psi_{i}$ and the diagonalized $\widetilde{\psi}_{i}$. One can prescribe free data for those components of $\widetilde{\psi}_{i}$ that travel inwards. These resulting fields are then translated back into boundary conditions for the components of $\psi_{i}$.

The structure of the Weyl system (5g) shows that at a given boundary there will be two ingoing modes, two outgoing modes and one mode travelling tangential to the time-like boundary. The boundary conditions must be such that they do not inject any constraint violating modes into the domain because these would cause the violation of the constraints across the entire domain and hence render the computation useless.

Thus we must also perform a characteristic analysis of the constraint propagation system. The only constraint equation of interest is

whose propagation equation has the principal part

The characteristics of the components of this constraint propagation equation for the components $\{G_{0},G_{1},G_{2}\}$ are $\{-\beta,0,\beta\}$. It is of no surprise that there is a correspondence between the characteristics of the evolution and subsidiary systems, see for example [2]. These equations have one ingoing mode which is the one that needs to be suppressed at the boundary. The vanishing of this mode leads to a condition on the ingoing modes for the $\psi_{i}$ which eliminates one degree of freedom in the choice of boundary data. Thus, at every boundary there is exactly one degree of freedom, i.e., one complex valued function, that can be specified freely on the boundary. It corresponds exactly to the ingoing spin-2 component of the Weyl spinor with respect to the boundary. We do not elucidate this procedure further here and again refer to [3] for a full account.

In this section we summarize the results obtained in [10] and put forward the analytical expressions required to compute the Bondi-Sachs energy on arbitrary cuts of $\mathscr{I}^{+}$.

First, we need to transform to a Bondi frame, say $\{L^{a},N^{a},M^{a},\overline{M}^{a}\}$, which is a null tetrad such that $N^{a}$ points along the null generators¹¹1Note, that here we divert in notation with [28] where $N^{a}=-\nabla^{a}\Theta$., and that $M^{a}$ is tangent to a given cut $C$. These conditions are expressed as equations on $\mathscr{I}^{+}$

where $A$ is a positive non-zero scalar on $\mathscr{I}^{+}$. The inclusion of this factor is necessary for maintaining the GHP formalism. However, even if one would fix the choice of a frame in some way then $A$ would show up as a normalisation factor. This frame can be thought of as a “physical detector” on $\mathscr{I}^{+}$, which can “detect” the radiation coming from the physical space-time.

Next we consider the Bondi-Sachs momentum 4-covector components, henceforth called Bondi components, where each component is defined as an integral over the cut $C$. The work in [10] results in a generalised expression for these components and represents them in a conformally invariant GHP (cGHP) formalism. As there are many quantities in these equations that need an explanation, we first present them all and then unpack their meanings. The Bondi components are obtained from integration over $C$ through

where

is closely related to Bondi’s news function and $U$ is an asymptotic translation obtained from the equations

To unpack all this, we note the presence of the spin-coefficients $\rho^{\prime},\sigma,\tau$ from the Newman-Penrose formalism, the $\eth$ and $\eth^{\prime}$ operators from the associated GHP formalism and their conformally invariant counterparts $\eth_{c},\overline{\eth}_{c}$ from the cGHP formalism. The quantity $\Phi_{20}$ is a component of the spinor $\Phi_{ABA^{\prime}B^{\prime}}$ which essentially represent the trace-free Ricci tensor and $\psi_{2}$ is the spin-zero component of the gravitational spinor. Finally, $K$ is related to the Gauß curvature $k$ of the cut which is given on $C$ by $k:=2K=2(\Phi_{11}+\Lambda-\rho\rho^{\prime})$, which is a real quantity, where $\Lambda$ is the scalar curvature and $\Phi_{11}$ is another component of $\Phi_{ABA^{\prime}B^{\prime}}$. We note that the three terms in the parentheses in Eq. (11) are together referred to as the mass aspect.

Eq. (11) is manifestly conformally invariant, assumes no conditions on the conformal factor $\Theta$, which is used to compactify the space-time, or spin-coefficients and reduces to the original definition when the specialization in [28] is taken. $U$ represents an infinitesimal translation. The translations form a subgroup of the super-translation group, which itself is a subgroup of the BMS group. The quantity $\mathcal{Q}$ arises naturally from action of the commutator $[\eth_{c},\overline{\eth}_{c}]$ on properly weighted quantities and the quantity $\mathcal{R}$ is uniquely determined by $\mathcal{Q}$ and is referred to as the co-curvature.

If the cut was a unit sphere and if we were in standard polar coordinates, then the equation for $U$ would reduce to $\eth^{2}U=0$ and it would have four independent solutions, namely the first four spherical harmonics $Y_{lm}$ ($l=0,\;m=0$ and $l=1,\;m=-1,0,1$.) These four choices for $U$ correspond to asymptotic time and space translations, and when normalized with a certain Lorentzian metric lead to the energy and momenta, respectively, when integrated against the mass aspect as given by Eq. (11) in a Bondi frame.

It is useful to present the Bondi-Sachs mass-loss

where $\mathcal{N}:=N+\bar{\mathcal{R}}$ and $\mathscr{I}_{1}^{2}$ is the part of null infinity between two cuts $C_{1}$ and $C_{2}$ on which the corresponding component is evaluated. When $U$ represents an appropriately normalized asymptotic time translation, Eq. (14) relates the difference in Bondi energy at two different cuts of $\mathscr{I}^{+}$ to the “energy flux” between them due to the outgoing gravitational radiation. This will be numerically checked in Sec. 7.3.

In this paper we consider only the Bondi energy, where the associated $U$ is equal to the conformal factor $\Omega$ that scales the metric of a given cut $C$ of $\mathscr{I}^{+}$ to the unit 2-sphere metric [10].

In the previous section we outlined how to solve an IBVP for the GCFE, taking proper care of the boundary conditions. We now discuss our choice of initial data. We choose to evolve Schwarzschild initial data from a space-like initial hypersurface as laid out in [16] and following the previously implemented numerical evolution in [3]. This choice of initial data, evolved with the GCFE, was analytically shown [16] to remain smooth and free of degeneracies up to and including null infinity. Even with gravitational waves impinging on this space-time region from the time-like outer boundary, we have demonstrated in [3] that the evolution of these initial data remains valid up to $i^{+}$, where the code breaks down due to the physical singularity there.

Even though the code is written to handle a general 3D problem for efficiency we now assume that the space-time is axisymmetric to reduce the problem to $2+1$ dimensions. This has the consequence that the spectral coefficients of a function in the spin-weighted spherical harmonic basis ${}_{s}Y_{lm}$ vanish when $m=0$.

We assume that the initial hypersurface can be foliated into concentric spheres and we adapt the coordinate system to this choice by labeling each sphere with a radial coordinate $r$ and introducing standard polar coordinates $(\theta,\phi)$ on each sphere. We define complex derivative operators on each sphere by

By deeming these to be complex null vectors we have implicitly introduced a fiducial metric on each sphere which makes it into a unit-sphere. This has nothing to do with the actual metric on each sphere induced from the space-time metric and exists for the sole purpose of having available a framework for using a numerical $\eth$-formalism based on the unit-sphere.

This structure on the initial hypersurface is carried along by the evolution along the conformal geodesics which are parameterised by the coordinate $t$. Thus, we use as a basis in each tangent space the vectors $\{\partial_{t},\partial_{r},\widetilde{\delta},\widetilde{\delta}^{\prime}\}$ with dual basis $\{\text{d}t,\text{d}r,\widetilde{\mathbf{m}},\overline{\widetilde{\mathbf{m}}}\}$, where $\{\widetilde{\delta},\,\widetilde{\delta}^{\prime}\}$ are tangent to the spheres of constant $(t,r)$. We assume the range of radial values lie in an interval so that the left and right time-like boundaries at a fixed $t$ will be spheres of constant $r$, denoted by $r_{l}$ and $r_{r}$ respectively with $r_{l}<r_{r}$.

The purpose of this section is to numerically investigate the response of the Schwarzschild black hole to the gravitational wave and how the Bondi energy along $\mathscr{I}^{+}$ is affected. We evolve as close as possible to $i^{+}$, at which point we expect the evolution to break down since, in the exact Schwarzschild space-time, it is a singular point. In fact, along the time-like conformal geodesic that goes through $i^{+}$, $\psi_{2}$ takes the form [16]

which diverges exactly at time-like infinity located at $t=2/(\sqrt{5}m)$.

Unless otherwise stated, we choose the initial mass of the black hole $m=0.5$ and solve the equations using a $r$-resolution of $401$ points, a $\theta$-resolution of $33$ points with an associated $l_{max}$ of $22$, a time-step of $\Delta t=c\,\Delta r$ with CFL number $c=0.5$. The location of $i^{+}$ for the exact Schwarzschild space-time can be calculated a priori using results of Friedrich [16]. We have $t_{i^{+}}=2/(\sqrt{5}m)\approx 1.789$ and $r_{i^{+}}=(3+\sqrt{5})m/4\approx 0.655$. These should be reasonable indications for the location of time-like infinity also in the perturbed cases.

In this paper we presented a framework for studying the interaction between gravitational waves hitting a black hole. This is achieved by setting up and solving an initial boundary value problem for the general conformal field equations. We described the effect of a gravitational wave impinging on an initially spherically symmetric (Schwarzschild) black hole. The advantage of using the GCFE together with the conformal Gauß gauge is that the time-like boundary eventually reaches and intersects null infinity so that the full information of the outgoing radiation becomes available. This allows us to evaluate the expressions for the Bondi-Sachs energy-momentum and numerically verify the validity of the mass loss formula directly on $\mathscr{I}$.

After the gravitational wave is sent into the space-time we observe a “periodic” behaviour in the components of the Weyl spinor (and other field quantities as well). It is tempting to relate these to the quasinormal modes of a linearly perturbed Schwarzschild black hole which have been extensively studied previously. However, this is not straightforward since those satisfy different boundary conditions from what we find in our simulation. While the linear modes are forced to vanish on the horizon, our excitation definitely persist on and even inside the horizon. This raises the question as to whether there is a difference in the characteristic frequencies between the linear and our non-linear modes and how would one find out?

The quasinormal modes are commonly computed by solving the radial Regge-Wheeler-Zerilli equation [29, 36] for a master function representing a linear perturbation of the Schwarzschild metric. This is then presented in Schwarzschild time along curves of constant Schwarzschild radius. However, in our situation we have neither, as our non-linear gravitational perturbations couple to the geometry, causing a physical deviation from the Schwarzschild metric. This makes it difficult to make physically meaningful comparisons.

Due to the use of the conformal Gauß gauge which has the property that the relationship between the physical (proper or retarded) times and the conformal time parameter is ultimately exponential we see the “ringing” currently only for three periods. This is definitely not enough to make any long-term comparisons with the linear regime. However, we pointed out that rescaling the conformal parameter could be used to extend the simulation, in principle indefinitely. Another possibility would be to abolish the conformal Gauß gauge altogether and use a different choice. However, this would entail non-trivial changes to the system of equations. In particular, we would lose the simplicity of the GCFE with the clear splitting into wave-type equations for the Weyl curvature and advection equations for the remaining variables. This would have the consequence that the IBVP would become much more complicated and one would have to discuss much more complicated boundary conditions.

The framework as presented here seems to be ideal to study the interaction of black holes and gravitational waves in a very clean setting. The possibility of setting up initially unperturbed situations and injecting the perturbation from the boundary opens the way to study many more questions of principle. For example, our next task will be to study whether we can “kick” a black hole by injecting gravitational waves with a linear momentum into the space-time. What would the response be? It should be straightforward also to replace the initial Schwarzschild black hole with a Kerr or Reissner-Nordström black hole. This would allow us to study phenomena such as super-radiance, the effect of gravitational perturbations on the inner Cauchy horizons and the stability of the solutions in general.