### 5.3 Integrating null surfaces backwards in time

Anninos et al. [7], Libson et al. [103], and Walker [162] introduced the important concept of explicitly (numerically) finding the event horizon as a null surface in spacetime. They observed that if we parameterize the event horizon with any level-set function satisfying the basic level-set definition (3), then the condition for the surface to be null is just
Applying a 3 + 1 decomposition to this then gives a quadratic equation which can be solved to find the time evolution of the level-set function,

Alternatively, assuming the event horizon in each slice to be a Strahlkörper in the manner of Section 2.2, we can define a suitable level-set function by Equation (7). Substituting this definition into Equation (10) then gives an explicit evolution equation for the horizon shape function,

Surfaces near the event horizon share the same “attraction” property discussed in Section 5.2 for geodesics near the event horizon. Thus by integrating either surface representation (10) or (11) backwards in time, we can refine an initial guess into a very accurate approximation to the event horizon.

In contrast to the null geodesic equation (8), neither Equation (10) nor Equation (11) contain any derivatives of the 4-metric (or equivalently the 3 + 1 geometry variables). This makes it much easier to integrate these latter equations accurately. This formulation of the event-horizon finding problem also completely eliminates the tangential-drifting problem discussed in Section 5.2, since the level-set function only parameterizes motion normal to the surface.

#### 5.3.1 Error bounds: Integrating a pair of surfaces

For a practical algorithm, it is useful to integrate a pair of trial null surfaces backwards: an “inner-bound” one which starts (and thus always remains) inside the event horizon and an “outer-bound” one which starts (and thus always remains) outside the event horizon. If the final slice contains an apparent horizon then any 2-surface inside this can serve as our inner-bound surface. However, choosing an outer-bound surface is more difficult.

It is this desire for a reliable outer bound on the event horizon position that motivates our requirement (Section 4) for the final slice (or at least its strong-field region) to be approximately stationary: In the absence of time-dependent equations of state or external perturbations entering the system, this requirement ensures that, for example, any surface substantially outside the apparent horizon can serve as an outer-bound surface.

Assuming we have an inner- and an outer-bound surface on the final slice, the spacing between these two surfaces after some period of backwards integration then gives an error bound for the computed event horizon position. Equivalently, a necessary (and, if there are no other numerical problems, sufficient) condition for the event-horizon finding algorithm to be accurate is that the backwards integration must have proceeded far enough for the spacing between the two trial surfaces to be “small”. For a reasonable definition of “small”, this typically takes at least of backwards integration, with or more providing much higher accuracy.

In some cases it is difficult to obtain a long enough span of numerical data for this backwards integration. For example, in some simulations of binary black hole collisions, the evolution becomes unstable and crashes soon after a common apparent horizon forms. This means that we cannot compute an accurate event horizon for the most interesting region of the spacetime, that which is close to the black-hole merger. There is no good solution to this problem except for the obvious one of developing a stable (or less-unstable) simulation that can be continued for a longer time.

#### 5.3.2 Explicit Strahlkörper surface representation

The initial implementations of the “integrate null surface backwards” algorithm by Anninos et al. [7], Libson et al. [103], and Walker [162] were based on the explicit Strahlkörper surface integration formula (11), further restricted to axisymmetry.

For a single black hole the coordinate choice is straightforward. For the two-black-hole case, the authors used topologically cylindrical coordinates , where the two black holes collide along the axisymmetry () axis. Based on the symmetry of the problem, they then assumed that the event horizon shape could be written in the form

in each slice.

This spacetime’s event horizon has the now-classic “pair of pants” shape, with a non-differentiable cusp along the “inseam” (the axis ) where new generators join the surface. The authors tried two ways of treating this cusp numerically:

• Since the cusp’s location is known a priori, it can be treated as a special case in the angular finite differencing, using one-sided numerical derivatives as necessary.
• Alternatively, in 1994 Thorne suggested calculating the union of the event horizon and all its null generators (including those which have not yet joined the surface). This “surface” has a complicated topology (it self-intersects along the cusp), but it is smooth everywhere. This is illustrated by Figure 4, which shows a cross-section of this surface in a single slice, for a head-on binary black hole collision. For comparison, Figure 5 shows a perspective view of part of the event horizon and some of its generators, for a similar head-on binary black hole collision.

Caveny et al. [4446] implemented the “integrate null surfaces backwards” algorithm for fully generic numerically-computed spacetimes using the explicit Strahlkörper surface integration formula (11). To handle moving black holes, they recentered each black hole’s Strahlkörper parameterization (4) on the black hole’s coordinate centroid at each time step.

For single-black-hole test cases (Kerr spacetime in various coordinates), they report typical accuracies of a few percent in determining the event horizon position and area. For binary-black-hole test cases (Kastor–Traschen extremal-charge black hole coalescence with a cosmological constant), they detect black hole coalescence (which appears as a bifurcation in the backwards time integration) by the “necking off” of the surface. Figure 6 shows an example of their results.

#### 5.3.3 Level-set parameterization

Caveny et al. [4445] and Diener [60] (independently) implemented the “integrate null surfaces backwards” algorithm for fully generic numerically-computed spacetimes, using the level-set function integration formula (10). Here the level-set function is initialized on the final slice of the evolution and evolved backwards in time using Equation (10) on (conceptually) the entire numerical grid. (In practice, only a smaller box containing the event horizon need be evolved.)

This surface parameterization has the advantage that the event-horizon topology and (non-)smoothness are completely unconstrained, allowing the numerical study of configurations such as toroidal event horizons (discussed in Section 4). It is also convenient that the level-set function is defined on the same numerical grid as the spacetime geometry, so that no interpolation is needed for the evolution.

The major problem with this algorithm is that during the backwards evolution, spatial gradients in tend to steepen into a jump discontinuity at the event horizon, eventually causing numerical difficulty.

Caveny et al. [4445] deal with this problem by adding an artificial viscosity (i.e. diffusion) term to the level-set function evolution equation, smoothing out the jump discontinuity in . That is, instead of Equation (10), they actually evolve via

where is the right hand side of Equation (10) and is a generic 2nd order linear (elliptic) spatial differential operator, and is a (small) dissipation constant. This scheme works, but the numerical viscosity does seem to lead to significant errors (several percent) in their computed event-horizon positions and areas, and even failure to converge to the correct solution for some test cases (e.g. rapidly-spinning Kerr black holes).

Alternatively, Diener [60] developed a technique of periodically reinitializing the level-set function to approximately the signed distance from the event horizon. To do this, he periodically evolves

in an unphysical “pseudo-time” until an approximate steady state has been achieved. He reports that this works well in most circumstances but can significantly distort the computed event horizon if the isosurface (the current approximation to the event horizon) is only a few grid points thick in any direction, as typically occurs just around the time of a topology change in the isosurface. He avoids this problem by estimating the minimum thickness of this isosurface and, if it is below a threshold, deferring the reinitialization.

In various tests on analytical data, Diener [60] found this event-horizon finder, EHFinder, to be robust and highly accurate, typically locating the event horizon to much less than 1% of the 3-dimensional grid spacing. As an example of results obtained with EHFinder, Figure 7 shows two views of the numerically-computed event horizon for a spiraling binary black hole collision. As another example, Figure 8 shows the numerically-computed event and apparent horizons in the collapse of a rapidly rotating neutron star to a Kerr black hole. (The apparent horizons were computed using the AHFinderDirect code described in Section 8.5.7.)

EHFinder is implemented as a freely available module (“thorn”) in the Cactus computational toolkit (see Table 2). It originally worked only with the PUGH unigrid driver, but work is ongoing [61] to enhance it to work with the Carpet mesh-refinement driver [134131].