Anninos et al.  and Libson et al.  observed that while this instability is a problem for the “integrate null geodesics forwards in time” algorithm (it forces that algorithm to take quite short time steps when integrating the geodesics), we can turn it to our advantage by integrating the geodesics backwards in time: The geodesics will now converge on to the horizon10.
This event-horizon finding algorithm thus integrates a large number of such (future-pointing outgoing) null geodesics backwards in time, starting on the final numerically-generated slice. As the backwards integration proceeds, even geodesics which started far from the event horizon will quickly converge to it. This can be seen, for example, in Figures 2 and 3.
Unfortunately, this convergence property holds only for outgoing geodesics. In spherical symmetry the distinction between outgoing and ingoing geodesics is trivial but, as described by Libson et al. ,
[…] for the general 3D case, when the two tangential directions of the EH are also considered, the situation becomes more complicated. Here normal and tangential are meant in the 3D spatial, not spacetime, sense. Whether or not a trajectory can eventually be “attracted” to the EH, and how long it takes for it to become “attracted,” depends on the photon’s starting direction of motion. We note that even for a photon which is already exactly on the EH at a certain instant, if its velocity at that point has some component tangential to the EH surface as generated by, say, numerical inaccuracy in integration, the photon will move outside of the EH when traced backward in time. For a small tangential velocity, the photon will eventually return to the EH [… but] the position to which it returns will not be the original position.
This kind of tangential drifting is undesirable not just because it introduces inaccuracy in the location of the EH, but more importantly, because it can lead to spurious dynamics of the “EH” thus found. Neighboring generators may cross, leading to numerically artificial caustic points […].
Libson et al.  also observed:
Another consequence of the second order nature of the geodesic equation is that not just the positions but also the directions must be specified in starting the backward integration. Neighboring photons must have their starting direction well correlated in order to avoid tangential drifting across one another.
Libson et al.  give examples of the numerical difficulties that can result from these difficulties and conclude that this event-horizon finding algorithm
[…] is still quite demanding in finding an accurate history of the EH, although the difficulties are much milder than those arising from the instability of integrating forward in time.
Because of these difficulties, this algorithm has generally been supplanted by the “backwards surface” algorithm I describe next.
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License.