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 levelset 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 4metric (or equivalently the 3 + 1 geometry variables). This makes it much easier to integrate these latter equations accurately^{11}. This formulation of the eventhorizon finding problem also completely eliminates the tangentialdrifting problem discussed in Section 5.2, since the levelset function only parameterizes motion normal to the surface.
For a practical algorithm, it is useful to integrate a pair of trial null surfaces backwards: an “innerbound” one which starts (and thus always remains) inside the event horizon and an “outerbound” one which starts (and thus always remains) outside the event horizon. If the final slice contains an apparent horizon then any 2surface inside this can serve as our innerbound surface. However, choosing an outerbound 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 strongfield region) to be approximately stationary: In the absence of timedependent 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 outerbound surface.
Assuming we have an inner and an outerbound 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 eventhorizon 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 blackhole merger. There is no good solution to this problem except for the obvious one of developing a stable (or lessunstable) simulation that can be continued for a longer time.
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^{12}.
For a single black hole the coordinate choice is straightforward. For the twoblackhole 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 nowclassic “pair of pants” shape, with a nondifferentiable cusp along the “inseam” (the axis ) where new generators join the surface. The authors tried two ways of treating this cusp numerically:


Caveny et al. [44, 46] implemented the “integrate null surfaces backwards” algorithm for fully generic numericallycomputed 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 singleblackhole test cases (Kerr spacetime in various coordinates), they report typical accuracies of a few percent in determining the event horizon position and area. For binaryblackhole test cases (Kastor–Traschen extremalcharge 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.

Caveny et al. [44, 45] and Diener [60] (independently) implemented the “integrate null surfaces backwards” algorithm for fully generic numericallycomputed spacetimes, using the levelset function integration formula (10). Here the levelset 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 eventhorizon 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 levelset 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^{14}, eventually causing numerical difficulty.
Caveny et al. [44, 45] deal with this problem by adding an artificial viscosity (i.e. diffusion) term to the levelset 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 eventhorizon positions and areas^{15}, and even failure to converge to the correct solution for some test cases (e.g. rapidlyspinning Kerr black holes).Alternatively, Diener [60] developed a technique of periodically reinitializing the levelset function to approximately the signed distance from the event horizon. To do this, he periodically evolves
in an unphysical “pseudotime” 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 eventhorizon finder, EHFinder, to be robust and highly accurate, typically locating the event horizon to much less than 1% of the 3dimensional grid spacing. As an example of results obtained with EHFinder, Figure 7 shows two views of the numericallycomputed event horizon for a spiraling binary black hole collision. As another example, Figure 8 shows the numericallycomputed 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 meshrefinement driver [134, 131].


http://www.livingreviews.org/lrr20073 
This work is licensed under a Creative Commons AttributionNoncommercialNo Derivative Works 2.0 Germany License. Problems/comments to 