### 5.1 Integrating null geodesics forwards in time

The first generation of event-horizon finders were based directly on Hawking’s original definition of an event horizon: an event is within the black hole region of spacetime if and only if there is no future-pointing “escape route” null geodesic from to ; the event horizon is the boundary of the black hole region.

That is, as described by Hughes et al. [88], we numerically integrate the null geodesic equation

(where is an affine parameter) forwards in time from a set of starting events and check which events have “escaping” geodesics. For analytical or semi-analytical studies like that of Bishop [31], this is an excellent algorithm.

For numerical work it is straightforward to rewrite the null geodesic equation (8) as a coupled system of two first-order equations, giving the time evolution of photon positions and 3-momenta in terms of the 3 + 1 geometry variables , , , and their spatial derivatives. These can then be time-integrated by standard numerical algorithms. However, in practice several factors complicate this algorithm.

We typically only know the 3 + 1 geometry variables on a discrete lattice of spacetime grid points, and we only know the 3 + 1 geometry variables themselves, not their spatial derivatives. Therefore we must numerically differentiate the field variables, then interpolate the field variables and their spacetime derivatives to each integration point along each null geodesic. This is straightforward to implement, but the numerical differentiation tends to amplify any numerical noise that may be present in the field variables.

Another complicating factor is that the numerical computations generally only span a finite region of spacetime, so it is not entirely obvious whether or not a given geodesic will eventually reach . However, if the final numerically-generated slice contains an apparent horizon, we can use this as an approximation: Any geodesic which is inside this apparent horizon will definitely not reach , while any other geodesic may be assumed to eventually reach if its momentum is directed away from the apparent horizon. If the final slice (or at least its strong-field region) is approximately stationary, the error from this approximation should be small. I discuss this stationarity assumption further in Section 5.3.1.

#### 5.1.1 Spherically-symmetric spacetimes

In spherical symmetry this algorithm works well and has been used by a number of researchers. For example, Shapiro and Teukolsky [141142143144] used it to study event horizons in a variety of dynamical evolutions of spherically symmetric collapse systems. Figure 2 shows an example of the event and apparent horizons in one of these simulations.

#### 5.1.2 Non-spherically-symmetric spacetimes

In a non-spherically-symmetric spacetime, several factors make this algorithm very inefficient:

• Many trial events must be tried to accurately resolve the event horizon’s shape. (Hughes et al. [88] describe a 2-stage adaptive numerical algorithm for choosing the trial events so as to accurately locate the event horizon as efficiently as possible.)
• At each trial event we must try many different trial-geodesic starting directions to see if any of the geodesics escape to (or our numerical approximation to it). Hughes et al. [88] report needing only 48 geodesics per trial event in several nonrotating axisymmetric spacetimes, but about 750 geodesics per trial event in rotating axisymmetric spacetimes, with up to 3000 geodesics per trial event in some regions of the spacetimes.
• Finally, each individual geodesic integration requires many (short) time steps for an accurate integration, particularly in the strong-field region near the event horizon.

Because of these limitations, for non-spherically-symmetric spacetimes the “integrate null geodesics forwards” algorithm has generally been supplanted by the more efficient algorithms I describe in the following.