1 Part I Introduction

Systems with strong gravitational fields, particularly systems which may contain event horizons and/or apparent horizons, are a major focus of numerical relativity. The usual output of a numerical relativity simulation is some (approximate, discrete) representation of the spacetime geometry (the 4-metric and possibly its derivatives) and any matter fields, but not any explicit information about the existence, precise location, or other properties of any event/apparent horizons. To gain this information, we must explicitly find the horizons from the numerically-computed spacetime geometry. The subject of this review is numerical algorithms and codes for doing this, focusing on calculations done using the 3 + 1 ADM formalism [14, 163]. Baumgarte and Shapiro [27, Section 6] have also recently reviewed event and apparent horizon finding algorithms. The scope of this review is limited to the finding of event/apparent horizons and omits any but the briefest mention of the many uses of this information in gaining physical understanding of numerically-computed spacetimes.

In this review I distinguish between a numerical algorithm (an abstract description of a mathematical computation; also often known as a “method” or “scheme”), and a computer code (a “horizon finder”, a specific piece of computer software which implements a horizon finding algorithm or algorithms). My main focus is on the algorithms, but I also mention specific codes where they are freely available to other researchers.

In this review I have tried to cover all the major horizon finding algorithms and codes, and to accurately credit the earliest publication of important ideas. However, in a field as large and active as numerical relativity, it is not unlikely that I have overlooked and/or misdescribed some important research. I apologise to anyone whose work I’ve slighted, and I ask readers to help make this a truly “living” review by sending me corrections, updates, and/or pointers to additional work (either their own or others) that I should discuss in future revisions of this review.

The general outline of this review is as follows: In the remainder of Part I, I define notation and terminology (Section 1), discuss how 2-surfaces should be parameterized (Section 2), and outline some of the software-engineering issues that arise in modern numerical relativity codes (Section 3). I then discuss numerical algorithms and codes for finding event horizons (Part II) and apparent horizons (Part III). Finally, in the appendices I briefly outline some of the excellent numerical algorithms/codes available for two standard problems in numerical analysis, the solution of a single nonlinear algebraic equation (Appendix A) and the time integration of a system of ordinary differential equations (Appendix B).

2 Notation and Terminology

Except as noted below, I generally follow the sign and notation conventions of Wald [160]. I assume that all spacetimes are globally hyperbolic, and for event-horizon finding I further assume asymptotic flatness; in this latter context \({{\mathcal J}^ +}\) is future null infinity. I use the Penrose abstract-index notation, with summation over all repeated indices. 4-indices abc range over all spacetime coordinates {xa}, and 3-indices ijk range over the spatial coordinates {xi} in a spacelike slice t = constant. The spacetime coordinates are thus xa = (t, xi).

Indices uvw range over generic angular coordinates (θ, ϕ) on S2 or on a horizon surface. Note that these coordinates are conceptually distinct from the 3-dimensional spatial coordinates xi. Depending on the context, (θ, ϕ) may or may not have the usual polar-spherical topology. Indices ijk label angular grid points on S2 or on a horizon surface. These are 2-dimensional indices: a single such index uniquely specifies an angular grid point. δIJ is the Kronecker delta on the space of these indices or, equivalently, on surface grid points.

For any indices p and q, p and pq are the coordinate partial derivatives ∂/∂xp and 2/∂xp∂xq respectively; for any coordinates µ and ν, u and µν are the coordinate partial derivatives ∂/∂/µ and 2/∂µ∂ν respectively. Δ is the flat-space angular Laplacian operator on S2, while Δx refers to a finite-difference grid spacing in some variable x.

gab is the spacetime 4-metric, and gab the inverse spacetime 4-metric; these are used to raise and lower 4-indices. \(\Gamma _{ab}^c\) are the 4-Christoffel symbols. \({{\mathcal L}_\upsilon}\) is the Lie derivative along the 4-vector field va.

I use the 3 + 1 “ADM” formalism first introduced by Arnowitt, Deser, and Misner [14]; York [163] gives a general overview of this formalism as it is used in numerical relativity. gij is the 3-metric defined in a slice, and gij is the inverse 3-metric; these are used to raise and lower 3-indices. ∇i is the associated 3-covariant derivative operator, and \(\Gamma _{ij}^k\) are the 3-Christoffel symbols. α and βi are the 3 + 1 lapse function and shift vector respectively, so the spacetime line element is

$$d{s^2} = {g_{ab}}d{x^a}d{x^b}$$
(1)
$$= - ({\alpha ^2} - {\beta _i}{\beta ^i})d{t^2} + 2{\beta _i}d{x^i}dt + {g_{ij}}d{x^i}d{x^j}.$$
(2)

As is common in 3 + 1 numerical relativity, I follow the sign convention of Misner, Thorne, and Wheeler [112] and York [163] in defining the extrinsic curvature of the slice as \({K_{ij}} = - {1 \over 2}{{\mathcal L}_n}{g_{ij}} = - {\nabla _i}{n_j}\), where na is the future-pointing unit normal to the slice. (In contrast, Wald [160] omits the minus signs from this definition.) \(K\equiv {K_i^i}\) is the trace of the extrinsic curvature Kij. mADM is the ADM mass of an asymptotically flat slice.

I often write a differential operator as F = F (y, uy, uvy; gij, kgij, Kij), where the “;” notation means that F is a (generally nonlinear) algebraic function of the variable y and its 1st and 2nd angular derivatives, and that F also depends on the coefficients gij, kgij, and Kij at the apparent horizon position.

There are three common types of spacetimes/slices where numerical event or apparent horizon finding is of interest: spherically-symmetric spacetimes/slices, axisymmetric spacetimes/slices, and spacetimes/slices with no continuous spatial symmetries (no spacelike Killing vectors). I refer to the latter as “fully generic” spacetimes/slices.

In this review I use the abbreviations “ODE” for ordinary differential equation, “PDE” for partial differential equation, “CE surface” for constant-expansion surface, and “MOTS” for marginally outer trapped surface. Names in Small Capitals refer to horizon finders and other computer software.

When discussing iterative numerical algorithms, it is often convenient to use the concept of an algorithm’s “radius of convergence”. Suppose the solution space within which the algorithm is iterating is S. Then given some norm ‖ · ‖ on S, the algorithm’s radius of convergence about a solution sS is defined as the smallest r > 0 such that the algorithm will converge to the correct solution s for any initial guess g with ‖gs‖ ≤ r. We only rarely know the exact radius of convergence of an algorithm, but practical experience often provides a rough estimateFootnote 1.

3 2-Surface Parameterizations

3.1 Level-set-function parameterizations

The most general way to parameterize a 2-surface in a slice is to define a scalar “level-set function” F on some neighborhood of the surface, with the surface itself then being defined as the level set

$$F = 0\quad {\rm{on}}\;{\rm{the}}\;{\rm{surface}}.$$
(3)

Assuming the surface to be orientable, it is conventional to choose F so that F > 0 (F < 0) outside (inside) the surface. The choice of level-set function for a given surface is non-unique, but in general this is not a problem.

This parameterization is valid for any surface topology including time-dependent topologies. The 2-surface itself can then be found by a standard isosurface-finding algorithm such as the marching-cubes algorithm [105]. (This algorithm is widely used in computer graphics and is implemented in a number of widely-available software libraries.)

3.2 Strahlkörper parameterizations

Most apparent horizon finders, and some event-horizon finders, assume that each connected component of the apparent (event) horizon has S2 topology. With the exception of toroidal event horizons (discussed in Section 4), this is generally a reasonable assumption.

To parameterize an S2 surface’s shape, it is common to further assume that we are given (or can compute) some “local coordinate origin” point inside the surface such that the surface’s 3-coordinate shape relative to that point is a “Strahlkörper” (literally “ray body”, or more commonly “star-shaped region”), defined by Minkowski [138, Page 108] as

a region in n-D Euclidean space containing the origin and whose surface, as seen from the origin, exhibits only one point in any direction.

The Strahlkörper assumption is a significant restriction on the horizon’s coordinate shape (and the choice of the local coordinate origin). For example, it rules out the coordinate shape and local coordinate origin illustrated in Figure 1: a horizon with such a coordinate shape about the local coordinate origin could not be found by any horizon finder which assumes a Strahlkörper surface parameterization.

Figure 1
figure 1

This figure shows a cross-section of a coordinate shape (the thick curve) which is not a Strahlkörper about the local coordinate origin shown (the large dot). The dashed line shows a ray from the local coordinate origin, which intersects the surface in more than one point.

For event-horizon finding, algorithms and codes are now available which allow an arbitrary horizon topology with no Strahlkörper assumption (see the discussion in Section 5.3.3 for details). For apparent horizon finding, the flow algorithms discussed in Section 8.7 theoretically allow any surface shape, although many implementations still make the Strahlkörper assumption. Removing this assumption for other apparent horizon finding algorithms might be a fruitful area for further research.

Given the Strahlkörper assumption, the surface can be explicitly parameterized as

$$r = h(\theta, \phi),$$
(4)

where r is the Euclidean distance from the local coordinate origin to a surface point, (θ, ϕ) are generic angular coordinates on the horizon surface (or equivalently on S2), and the “horizon shape function” h : S2 → ℜ+ is a positive real-valued function on the domain of angular coordinates defining the surface shape. Given the choice of local coordinate origin, there is clearly a one-to-one mapping between Strahlkörper 2-surfaces and horizon shape functions.

There are two common ways to discretize a horizon shape function:

  • Spectral representation

    Here we expand the horizon shape function h in an infinite series in some (typically orthonormal) set of basis functions such as spherical harmonics Yℓm or symmetric trace-free tensorsFootnote 2,

    $$h(\theta, \phi) = \sum\limits_{\ell, m} {{a_{\ell m}}{Y_{\ell m}}(\theta, \phi).}$$
    (5)

    This series can then be truncated at some finite order max, and the Ncoeff = max(max+1) coefficients {aℓm} used to represent (discretely approximate) the horizon shape. For reasonable accuracy, max is typically on the order of 8 to 12.

  • Finite difference representation

    Here we choose some finite grid of angular coordinates {(θK, ϕK)}, K = 1, 2, 3, …, Nang on S2 (or equivalently on the surface)Footnote 3, and represent (discretely approximate) the surface shape by the Nang values

    $$\{h({\theta _{\rm{K}}},{\phi _{\rm{K}}})\} \quad {\rm{K = 1,2,3,}} \ldots, {N_{{\rm{ang}}}}.$$
    (6)

    For reasonable accuracy, Nang is typically on the order of a few thousand.

It is sometimes useful to explicitly construct a level-set function describing a given Strahlkörper. A common choice here is

$$F \equiv r - h(\theta, \phi).$$
(7)

3.3 Finite-element parameterizations

Another way to parameterize a 2-surface is via finite elements where the surface is modelled as a triangulated mesh, i.e. as a set of interlinked “vertices” (points in the slice, represented by their spatial coordinates {xi}), “edges” (represented by ordered pairs of vertices), and faces. Typically only triangular faces are used (represented as oriented triples of vertices).

A key benefit of this representation is that it allows an arbitrary topology for the surface. However, determining the actual surface topology (e.g. testing for whether or not the surface self-intersects) is somewhat complicated.

This representation is similar to that of Regge calculus [128, 72]Footnote 4, and can similarly be expected to show 2nd order convergence with the surface resolution.

4 Software-Engineering Issues

Historically, numerical relativists wrote their own codes from scratch. As these became more complex, many researchers changed to working on “group codes” with multiple contributors.

4.1 Software libraries and toolkits

More recently, particularly in work on fully generic spacetimes, where all three spatial dimensions must be treated numerically, there has been a strong trend towards the use of higher-level software libraries and modular “computational toolkits” such as Cactus [74] (http://www.cactuscode.org). These have a substantial learning overhead, but can allow researchers to work much more productively by focusing more on numerical relativity instead of computer-science and softwareengineering issues such as parameter-file parsing, parallelization, I/O, etc.

A particularly important area for such software infrastructure is mesh refinementFootnote 5. This is essential to much current numerical-relativity research but is moderately difficult to implement even in only one spatial dimension, and much harder in multiple spatial dimensions. There are now a number of software libraries providing multi-dimensional mesh-refinement infrastructure (sometimes combined with parallelization), such as those listed in Table 1. The Cactus toolkit can be used in either unigrid or mesh-refinement modes, the latter using a “mesh-refinement driver” such as PAGH or Carpet [134, 131] (http://www.carpetcode.org).

Table 1 This table lists some software toolkits for multi-dimensional mesh refinement. All these toolkits also provide parallelization.

In this review I point out event and apparent horizon finders which have been written in particular frameworks and comment on whether they work with mesh refinement.

4.2 Code reuse and sharing

Another important issue is that of code reuse and sharing. It is common for codes to be shared within a research group but relatively uncommon for them to be shared between different (competing) research groups. Even apart from concerns about competitive advantage, without a modular structure and clear documentation it is difficult to reuse another group’s code. The use of a common computational toolkit can greatly simplify such reuse.

If such reuse can be accomplished, it becomes much easier for other researchers to build on existing work rather than having to “reinvent the wheel”. As well as the obvious ease of reusing existing code that (hopefully!) already works and has been thoroughly debugged and tested, there is another — less obvious — benefit of code sharing: It greatly eases the replication of past work, which is essential as a foundation for new development. That is, without access to another researcher’s code, it can be surprisingly difficult to replicate her results because the success or failure of a numerical algorithm frequently depends on subtle implementation details not described in even the most complete of published papers.

Event and apparent horizon finders are excellent candidates for software reuse: Many numerical-relativity researchers can benefit from using them, and they have a relatively simple interface to an underlying numerical-relativity simulation. Even if a standard computational toolkit is not used, this relatively simple interface makes it fairly easy to port an event or apparent horizon finder to a different code.

Table 2 lists event and apparent horizon finders which are freely available to any researcher.

Table 2 This table lists event and apparent horizon finders which are freely available to any researcher, along with the cvs repositories or web pages from which they may be obtained.

4.3 Using multiple event/apparent horizon finders

It is useful to have multiple event or apparent horizon finders available: Their strengths and weaknesses may complement each other, and the extent of agreement or disagreement between their results can help to estimate the numerical accuracy. For example, Figure 11 shows a comparison between the irreducible masses of apparent horizons in a binary black hole coalescence simulation (Alcubierre et al. [5], [Figure 4b]), as computed by two different apparent horizon finders in the Cactus toolkit, AHFinder and AHFinderDireot. In this case the two agree to within about 2% for the individual horizons and 0.5% for the common horizon.

5 Part II Finding Event Horizons

6 Introduction

The black hole region of an asymptotically-flat spacetime is defined [81, 82] as the set of events from which no future-pointing null geodesic can reach future null infinity (\({{\mathcal J}^ +}\)). The event horizon is defined as the boundary of the black hole region. The event horizon is a null surface in spacetime with (in the words of Hawking and Ellis [82, Page 319]) “a number of nice properties” for studying the causal stucture of spacetime.

The event horizon is a global property of an entire spacetime and is defined nonlocally in time: The event horizon in a slice is defined in terms of (and cannot be computed without knowing) the full future development of that slice.

In practice, to find an event horizon in a numerically-computed spacetime, we typically instrument a numerical evolution code to write out data files of the 4-metric. After the evolution (or at least the strong-field region) has reached an approximately-stationary final state, we then compute a numerical approximation to the event horizon in a separate post-processing pass, using the 4-metric data files as inputs.

As a null surface, the event horizon is necessarily continuous. In theory it need not be anywhere differentiableFootnote 6, but in practice this behavior rarely occursFootnote 7: The event horizon is generally smooth except for possibly a finite set of “cusps” where new generators join the surface; the surface normal has a jump discontinuity across each cusp. (The classic example of such a cusp is the “inseam” of the “pair of pants” event horizon illustrated in Figures 4 and 5.)

A black hole is defined as a connected component of the black hole region in a 3 + 1 slice. The boundary of a black hole (the event horizon) in a slice is a 2-dimensional set of events. Usually this has 2-sphere (S2) topology. However, numerically simulating rotating dust collapse, Abrahams et al. [1] found that in some cases the event horizon in a slice may be toroidal in topology. Lehner et al. [99], and Husa and Winicour [91] have used null (characteristic) algorithms to give a general analysis of the event horizon’s topology in black hole collisions; they find that there is generically a (possibly brief) toroidal phase before the final 2-spherical state is reached. Lehner et al. [100] later calculated movies showing this behavior for several asymmetric black hole collisions.

7 Algorithms and Codes for Finding Event Horizons

There are three basic event-horizon finding algorithms:

  • Integrate null geodesics forwards in time (Section 5.1).

  • Integrate null geodesics backwards in time (Section 5.2).

  • Integrate null surfaces backwards in time (Section 5.3).

I describe these in detail in the following.

7.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 \({\mathcal P}\) is within the black hole region of spacetime if and only if there is no future-pointing “escape route” null geodesic from \({\mathcal P}\) to \({{\mathcal J}^ +}\); 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

$${{{d^2}{x^a}} \over {d{\lambda ^2}}} + \Gamma _{bc}^a{{d{x^b}} \over {d\lambda}}{{d{x^c}} \over {d\lambda}} = 0$$
(8)

(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 α, βi, gij, and their spatial derivatives. These can then be time-integrated by standard numerical algorithmsFootnote 8. 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 implementFootnote 9, 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 \({{\mathcal J}^ +}\). 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 \({{\mathcal J}^ +}\), while any other geodesic may be assumed to eventually reach \({{\mathcal J}^ +}\) 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.

7.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 [141, 142, 143, 144] 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.

Figure 2
figure 2

This figure shows part of a simulation of the spherically symmetric collapse of a model stellar core (\(a\,\Gamma = {5 \over 3}\) polytrope) to a black hole. The event horizon (shown by the dashed line) was computed using the “integrate null geodesics forwards” algorithm described in Section 5.1 ; solid lines show outgoing null geodesics. The apparent horizon (the boundary of the trapped region, shown shaded) was computed using the zero-finding algorithm discussed in Section 8.1 . The dotted lines show the world lines of Lagrangian matter tracers and are labeled by the fraction of baryons interior to them. Figure reprinted with permission from [142]. © 1980 by the American Astronomical Society.

7.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 \({{\mathcal J}^ +}\) (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.

7.2 Integrating null geodesics backwards in time

It is well-known that future-pointing outgoing null geodesics near the event horizon tend to diverge exponentially in time away from the event horizon. Figure 3 illustrates this behavior for Schwarzschild spacetime, but the behavior is actually quite generic.

Figure 3
figure 3

This figure shows a number of light cones and future-pointing outgoing null geodesics in a neighborhood of the event horizon in Schwarzschild spacetime, plotted in ingoing Eddington-Finkelstein coordinates (t, r). (These coordinates are defined by the conditions that t + r is an ingoing null coordinate, while r is an areal radial coordinate.) Note that for clarity the horizontal scale is expanded relative to the vertical scale, so the light cones open by more than ±45°. All the geodesics start out close together near the event horizon; they diverge away from each other exponentially in time (here with an e-folding time of 4m near the horizon). Equivalently, they converge towards each other if integrated backwards in time (downwards on the page).

Anninos et al. [7] and Libson et al. [103] 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 horizonFootnote 10.

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. [103],

[…] 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. [103] 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. [103] 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.

7.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 F satisfying the basic level-set definition (3), then the condition for the surface F = 0 to be null is just

$${g^{ab}}{\partial _a}F{\partial _b}F = 0.$$
(9)

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,

$${\partial _t}F = {{- {g^{ti}}{\partial _i}F + \sqrt {{{({g^{ti}}{\partial _i}F)}^2} - {g^{tt}}{g^{ij}}{\partial _i}f{\partial _j}f}} \over {{g^{tt}}}}.$$
(10)

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 F by Equation (7). Substituting this definition into Equation (10) then gives an explicit evolution equation for the horizon shape function,

$${\partial _t}h = {{- {g^{tr}} + {g^{ru}}{\partial _u}h + \sqrt {{{({g^{tr}} - {g^{tu}}{\partial _u}h)}^2} - {g^{tt}}({g^{rr}} - 2{g^{ru}}{\partial _u}h + {g^{uv}}{\partial _u}h{\partial _v}h)}} \over {{g^{tt}}}}.$$
(11)

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 accuratelyFootnote 11. 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.

7.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 15mADM of backwards integration, with 20mADM 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.

7.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 axisymmetryFootnote 12.

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

$$\rho = h(z)$$
(12)

in each t = constant slice.

This spacetime’s event horizon has the now-classic “pair of pants” shape, with a non-differentiable cusp along the “inseam” (the z axis ρ = 0) 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)Footnote 13. 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.

Figure 4
figure 4

This figure shows a view of the numerically-computed event horizon in a single slice, together with the locus of the event horizon’s generators that have not yet joined the event horizon in this slice, for a head-on binary black hole collision. Notice how the event horizon is non-differentiable at the cusp where the new generators join it. Figure reprinted with permission from [103]. © 1996 by the American Physical Society.

Figure 5
figure 5

This figure shows a perspective view of the numerically-computed event horizon, together with some of its generators, for the head-on binary black hole collision discussed in detail by Matzner et al. [108]. Figure courtesy of Edward Seidel.

Caveny et al. [44, 46] 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.

Figure 6
figure 6

This figure shows the cross-section of the numerically-computed event horizon in each of five different slices, for the head-on collision of two extremal Kastor-Traschen black holes. Figure reprinted with permission from [46]. © 2003 by the American Physical Society.

7.3.3 Level-set parameterization

Caveny et al. [44, 45] 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 F 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 F 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 F tend to steepen into a jump discontinuity at the event horizonFootnote 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 level-set function evolution equation, smoothing out the jump discontinuity in F. That is, instead of Equation (10), they actually evolve F via

$${\partial _t}F = {\rm{rh}}{{\rm{s}}_{{\rm{Eq}}{\rm{.(10)}}}} + {\varepsilon ^2}{\nabla ^2}F,$$
(13)

where rhs Eq. (10) is the right hand side of Equation (10) and ∇2 is a generic 2nd order linear (elliptic) spatial differential operator, and ε > 0 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 areasFootnote 15, 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

$${\partial _\lambda}F = - {F \over {\sqrt {{F^2} + 1}}}(\vert \nabla F\vert - 1)$$
(14)

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 F = 0 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 AHFinderDireot code described in Section 8.5.7.)

Figure 7
figure 7

This figure shows two views of the numerically-computed event horizon’s cross-section in the orbital plane for a spiraling binary black hole collision. The two orbital-plane dimensions are shown horizontally; time runs upwards. The initial data was constructed to have an approximate helical Killing vector, corresponding to black holes in approximately circular orbits (the D =18 case of Grandclément et al. [78]), with a proper separation of the apparent horizons of 6.9 m. (The growth of the individual event horizons by roughly a factor of 3 in the first part of the evolution is an artifact of the coordinate choice — the black holes are actually in a quasi-equilibrium state.) Figure courtesy of Peter Diener.

Figure 8
figure 8

This figure shows the polar and equatorial radii of the event horizon (solid lines) and apparent horizon (dashed lines) in a numerical simulation of the collapse of a rapidly rotating neutron star to form a Kerr black hole. The dotted line shows the equatorial radius of the stellar surface. These results are from the D4 simulation of Baiotti et al. [21]. Notice how the event horizon grows from zero size while the apparent horizon first appears at a finite size and grows in a spacelike manner. Notice also that both surfaces are flattened due to the rotation. Figure reprinted with permission from [21]. © 2005 by the American Physical Society.

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 [134, 131].

8 Summary of Algorithms/Codes for Finding Event Horizons

In spherical symmetry the “integrate null geodesics forwards” algorithm (Section 5.1) can be used, although the “integrate null geodesics backwards” and “integrate null surfaces backwards” algorithms (Sections 5.2 and 5.3 respectively) are more efficient.

In non-spherically-symmetric spacetimes the “integrate null surfaces backwards” algorithm (Section 5.3) is clearly the best algorithm known: It is efficient, accurate, and fairly easy to implement. For generic spacetimes, Diener’s event-horizon finder EHFinder [60] is particularly notable as a freely available implementation of this algorithm as a module (“thorn”) in the widely-used Cactus computational toolkit (see Table 2).

9 Part III Finding Apparent Horizons

10 Introduction

10.1 Definition

Given a (spacelike) 3 + 1 slice, a “trapped surface” is defined as a smooth closed 2-surface in the slice whose future-pointing outgoing null geodesics have negative expansion Θ. The “trapped region” in the slice is then defined as the union of all trapped surfaces, and the “apparent horizon” is defined as the outer boundary of the trapped region.

While mathematically elegant, this definition is not convenient for numerically finding apparent horizons. Instead, an alternate definition can be used: A MOTS is defined as a smooth (differentiable) closed orientable 2-surface in the slice whose future-pointing outgoing null geodesics have zero expansion Θ.Footnote 16 There may be multiple MOTSs in a slice, either nested within each other or intersectingFootnote 17. An apparent horizon is then defined as an outermost MOTS in a slice, i.e. a MOTS not contained in any other MOTS. Kriele and Hayward [98] have shown that subject to certain technical conditions, this definition is equivalent to the “outer boundary of the trapped region” one.

Notice that the apparent horizon is defined locally in time (it can be computed using only Cauchy data on a spacelike slice), but (because of the requirement that it be closed) non-locally in spaceFootnote 18. Hawking and Ellis [82] discuss the general properties of MOTSs and apparent horizons in more detail.

Except for flow algorithms (Section 8.7), all numerical “apparent horizon” finding algorithms and codes actually find MOTSs, and hereinafter I generally follow the common (albeit sloppy) practice in numerical relativity of blurring the distinction between an MOTS and an apparent horizon.

10.2 General properties

Given certain technical assumptions (including energy conditions), the existence of any trapped surface (and hence any apparent horizon) implies that the slice contains a black holeFootnote 19. (The converse of this statement is not true: An arbitrary (spacelike) slice through a black hole need not contain any apparent horizonFootnote 20.) However, if an apparent horizon does exist, it necessarily coincides with, or is contained in, an event horizon. In a stationary spacetime the event and apparent horizons coincide.

It is this relation to the event horizon which makes apparent horizons valuable for numerical computation: An apparent horizon provides a useful approximation to the event horizon in a slice, but unlike the event horizon, an apparent horizon is defined locally in time and so can be computed “on the fly” during a numerical evolution.

Given a family of spacelike 3 + 1 slices which foliate part of spacetime, the union of the slices’ apparent horizons (assuming they exist) forms a world-tubeFootnote 21. This world-tube is necessarily either null or spacelike. If it is null, this world-tube is slicing-independent (choosing a different family of slices gives the same world-tube, at least so long as each slice still intersects the world-tube in a surface with 2-sphere topology). However, if the world-tube is spacelike, it is slicing-dependent: Choosing a different family of slices will in general give a different world-tubeFootnote 22.

10.3 Trapping, isolated, and dynamical horizons

Hayward [83] introduced the important concept of a “trapping horizon” (roughly speaking an apparent horizon world-tube where the expansion becomes negative if the surface is deformed in the inward null direction) along with several useful variants. Ashtekar, Beetle, and Fairhurst [16], and Ashtekar and Krishnan [18] later defined the related concepts of an “isolated horizon”, essentially an apparent horizon world-tube which is null, and a “dynamical horizon”, essentially an apparent horizon world-tube which is spacelike.

These world-tubes obey a variety of local and global conservation laws, and have many applications in analyzing numerically-computed spacetimes. See the references cited above and also Dreyer et al. [63], Ashtekar and Krishnan [19, 20], Gourgoulhon and Jaramillo [76], Booth [36], and Schnetter, Krishnan, and Beyer [137] for further discussions, including applications to numerical relativity.

10.4 Description in terms of the 3 + 1 variables

In terms of the 3 + 1 variables, a MOTS (and thus an apparent horizon) satisfies the conditionFootnote 23

$$\Theta \equiv {\nabla _i}{s^i} + {K_{ij}}{s^i}{s^j} - K = 0,$$
(15)

where si is the outward-pointing unit 3-vector normal to the surfaceFootnote 24. Assuming the Strahlkörper surface parameterization (4), Equation (15) can be rewritten in terms of angular 1st and 2nd derivatives of the horizon shape function h,

$$\Theta \equiv \Theta (h,{\partial _u}h,{\partial _{uv}}h;{g_{ij}},{\partial _k}{g_{ij}},{K_{ij}}) = 0,$$
(16)

where Θ is a complicated nonlinear algebraic function of the arguments shown. (Shibata [146] and Thornburg [153, 156] give the Θ(h, uh, uvh) function explicitly.)

10.5 Geometry interpolation

Θ depends on the slice geometry variables gij, kgij, and Kij at the horizon positionFootnote 25. In practice these variables are usually only known on the (3-dimensional) numerical grid of the underlying numerical-relativity simulationFootnote 26, so they must be interpolated to the horizon position and, more generally, to the position of each intermediate-iterate trial shape the apparent horizon finding algorithm tries in the process of (hopefully) converging to the horizon position.

Moreover, usually the underlying simulation gives only gij and Kij, so gij must be numerically differentiated to obtain kgij. As discussed by Thornburg [156, Section 6.1], it is somewhat more efficient to combine the numerical differentiation and interpolation operations, essentially doing the differentiation inside the interpolatorFootnote 27.

Thornburg [156, Section 6.1] argues that for an elliptic-PDE algorithm (Section 8.5), for best convergence of the nonlinear elliptic solver the interpolated geometry variables should be smooth (differentiable) functions of the trial horizon surface position. He argues that that the usual Lagrange polynomial interpolation does not suffice here (in some cases his Newton’s-method iteration failed to converge) because this interpolation gives results which are only piecewise differentiableFootnote 28. He uses Hermite polynomial interpolation to avoid this problem. Cook and Abrahams [51], and Pfeiffer et al. [124] use bicubic spline interpolation; most other researchers either do not describe their interpolation scheme or use Lagrange polynomial interpolation.

10.6 Criteria for assessing algorithms

Ideally, an apparent horizon finder should have several attributes:

  • Robust: The algorithm/code should find an (the) apparent horizon in a wide range of numerically-computed slices, without requiring extensive tuning of initial guesses, iteration parameters, etc. This is often relatively easy to achieve for “tracking” the time evolution of an existing apparent horizon (where the most recent previously-found apparent horizon provides an excellent initial guess for the new apparent horizon position) but may be difficult for detecting the appearance of a new (outermost) apparent horizon in an evolution, or for initial-data or other studies where there is no “previous time step”.

  • Accurate: The algorithm/code should find an (the) apparent horizon to high accuracy and should not report spurious “solutions” (“solutions” which are not actually good approximations to apparent horizons or, at least, to MOTSs).

  • Efficient: The algorithm/code should be efficient in terms of its memory use and CPU time; in practice CPU time is generally the major constraint. It is often desirable to find apparent horizons at each time step (or, at least, at frequent intervals) during a numerical evolution. For this to be practical the apparent horizon finder must be very fast.

In practice, no apparent horizon finder is perfect in all these dimensions, so trade-offs are inevitable, particularly when ease of programming is considered.

10.7 Local versus global algorithms

Apparent horizon finding algorithms can usefully be divided into two broad classes:

  • Local algorithms are those whose convergence is only guaranteed in some (functional) neighborhood of a solution. These algorithms require a “good” initial guess in order to find the apparent horizon. Most apparent horizon finding algorithms are local.

  • Global algorithms are those which can (in theory, ignoring finite-step-size and other numerical effects) converge to the apparent horizon independent of any initial guess. Flow algorithms (Section 8.7) are the only truely global algorithms. Zero-finding in spherical symmetry (Section 8.1) and shooting in axisymmetry (Section 8.2) are “almost global” algorithms: They require only 1-dimensional searches, which (as discussed in Appendix A) can be programmed to be very robust and efficient. In many cases horizon pretracking (Section 8.6) can semi-automatically find an initial guess for a local algorithm, essentially making the local algorithm behave like an “almost global” one.

One might wonder why local algorithms are ever used, given the apparently superior robustness (guaranteed convergence independent of any initial guess) of global algorithms. There are two basic reasons:

  • In practice, local algorithms are much faster than global ones, particularly when “tracking” the time evolution of an existing apparent horizon.

  • Due to finite-step-size and other numerical effects, in practice even “global” algorithms may fail to converge to an apparent horizon. (That is, the algorithms may sometimes fail to find an apparent horizon even when one exists in the slice.)

11 Algorithms and Codes for Finding Apparent Horizons

Many researchers have studied the apparent horizon finding problem, and there are a large number of different apparent horizon finding algorithms and codes. Almost all of these require (assume) that any apparent horizon to be found is a Strahlkörper (Section 2) about some local coordinate origin; both finite-difference and spectral parameterizations of the Strahlkörper are common.

For slices with continuous symmetries, special algorithms are sometimes used:

  • Zero-Finding in Spherical Symmetry (Section 8.1)

    In spherical symmetry the apparent horizon equation (16) becomes a 1-dimensional nonlinear algebraic equation, which can be solved by zero-finding.

  • The Shooting Algorithm in Axisymmetry (Section 8.2)

    In axisymmetry the apparent horizon equation (16) becomes a nonlinear 2-point boundary value ODE, which can be solved by a shooting algorithm.

Alternatively, all the algorithms described below for generic slices are also applicable to axisymmetric slices and can take advantage of the axisymmetry to simplify the implementation and boost performance.

For fully generic slices, there are several broad categories of apparent horizon finding algorithms and codes:

  • Minimization Algorithms (Section 8.3)

    These algorithms define a scalar norm on Θ over the space of possible trial surfaces. A generalpurpose scalar-function-minimization routine is then used to search trial-surface-shape space for a minimum of this norm (which should give Θ = 0).

  • Spectral Integral-Iteration Algorithms (Section 8.4)

    These algorithms expand the (Strahlkörper) apparent horizon shape function in a spherical-harmonic basis, use the orthogonality of spherical harmonics to write the apparent horizon equation as a set of integral equations for the spectral coefficients, and solve these equations using a functional-iteration algorithm.

  • Elliptic-PDE Algorithms (Section 8.5)

    These algorithms write the apparent horizon equation (16) as a nonlinear elliptic (boundary-value) PDE for the horizon shape and solve this PDE using (typically) standard elliptic-PDE numerical algorithms.

  • Horizon Pretracking (Section 8.6)

    Horizon pretracking solves a slightly more general problem than apparent horizon finding: Roughly speaking, the determination of the smallest E ≥ 0 such that the equation Θ = E has a solution, and the determination of that solution. By monitoring the time evolution of E and of the surfaces satisfying this condition, we can determine — before it appears — approximately where (in space) and when (in time) a new MOTS will appear in a dynamic numerically-evolving spacetime. Horizon pretracking is implemented as a 1-dimensional (binary) search using a slightly-modified elliptic-PDE apparent horizon finding algorithm as a “subroutine”.

  • Flow Algorithms (Section 8.7)

    These algorithms start with a large 2-surface (larger than any possible apparent horizon in the slice) and shrink it inwards using an algorithm which ensures that the surface will stop shrinking when it coincides with the apparent horizon.

I describe the major algorithms and codes in these categories in detail in the following.

11.1 Zero-finding in spherical symmetry

In a spherically symmetric slice, any apparent horizon must also be spherically symmetric, so the apparent horizon equation (16) becomes a 1-dimensional nonlinear algebraic equation Θ(h) = 0 for the horizon radius h. For example, adopting the usual (symmetry-adapted) polar-spherical spatial coordinates xi = (r, θ, ϕ), we have [154, Equation (B7)]

$$\Theta \equiv {{{\partial _r}{g_{\theta \theta}}} \over {{g_{\theta \theta}}\sqrt {{g_{rr}}}}} - 2{{{K_{\theta \theta}}} \over {{g_{\theta \theta}}}} = 0.$$
(17)

Given the geometry variables grr, gθθ, rgθθ, and Kθθ, this equation may be easily and accurately solved using one of the zero-finding algorithms discussed in Appendix AFootnote 29.

Zero-finding has been used by many researchers, including [141, 142, 143, 144, 119, 47, 139, 9, 154, 155]Footnote 30. For example, the apparent horizons shown in Figure 2 were obtained using this algorithm. As another example, Figure 9 shows Θ(r) and h at various times in a (different) spherically symmetric collapse simulation.

Figure 9
figure 9

This figure shows the apparent horizons (actually MOTSs) for a spherically symmetric numerical evolution of a black hole accreting a narrow shell of scalar field, the 800.pqw1 evolution of Thornburg [155]. Part (a) of this figure shows Θ(r) (here labelled H) for a set of equally-spaced times between t=19 and t=20, while Part (b) shows the corresponding MOTS radius h(t) and the Misner-Sharp [111], [112, Box 23.1] mass m(h) internal to each MOTS. Notice how two new MOTSs appear when the local minimum in Θ(r) touches the 0=0 line, and two existing MOTS disappear when the local maximum in Θ(r) touches the Θ=0 line.

11.2 The shooting algorithm in axisymmetry

In an axisymmetric spacetime we can use symmetry-adapted coordinates (θ, ϕ), so (given the Strahlkörper assumption) without further loss of generality we can write the horizon shape function as h = h(θ). The apparent horizon equation (16) then becomes a nonlinear 2-point boundary-value ODE for the horizon shape function h [146, Equation (1.1)]

$$\Theta \equiv \Theta (h,{\partial _\theta}h,{\partial _{\theta \theta}}h;{g_{ij}},{\partial _k}{g_{ij}},{K_{ij}}) = 0,$$
(18)

where Θ(h) is a nonlinear 2nd order (ordinary) differential operator in h as shown.

Taking the angular coordinate θ to have the usual polar-spherical topology, local smoothness of the apparent horizon gives the boundary conditions

$${\partial _\theta}h = 0\quad {\rm{at}}\theta {\rm{= 0}}\;{\rm{and}}\;\theta ={\theta _{{{\max}}}},$$
(19)

where θmax is π/2 if there is “bitant” reflection symmetry across the z = 0 plane, or π otherwise.

As well as the more general algorithms described in the following, this may be solved by a shooting algorithmFootnote 31:

  1. 1.

    Guess the value of h at one endpoint, say h(θ=0) ≡ h*.

  2. 2.

    Use this guessed value of h(θ=0) together with the boundary condition (19) there as initial data to integrate (“shoot”) the ODE (18) from θ=0 to the other endpoint θ=θmax. This can be done easily and efficiently using one of the ODE codes described in Appendix B.

  3. 3.

    If the numerically computed solution satisfies the other boundary condition (19) at θ=θmax to within some tolerance, then the just-computed h(θ) describes the (an) apparent horizon, and the algorithm is finished.

  4. 4.

    Otherwise, adjust the guessed value h(θ=θ) ≡ h* and try again. Because there is only a single parameter (h*) to be adjusted, this can be done using one of the 1-dimensional zero-finding algorithms discussed in Appendix A.

This algorithm is fairly efficient and easy to program. By trying a sufficiently wide range of initial guesses h* this algorithm can give a high degree of confidence that all apparent horizons have been located, although this, of course, increases the cost.

Shooting algorithms of this type have been used by many researchers, for example [159, 66, 2, 29, 30, 145, 3, 4].

11.3 Minimization algorithms

This family of algorithms defines a scalar norm ‖ · ‖ on the expansion Θ over the space of possible trial surfaces, typically the mean-squared norm

$$\Vert\Theta \Vert \equiv \int {{\Theta ^2}d\Omega,}$$
(20)

where the integral is over all solid angles on a trial surface.

Assuming the horizon surface to be a Strahlkörper and adopting the spectral representation (5) for the horizon surface, we can view the norm (20) as being defined on the space of spectral coefficients {aℓm}.

This norm clearly has a global minimum ‖Θ‖ = 0 for each solution of the apparent horizon equation (16). To find the apparent horizon we numerically search the spectral-coefficient space for this (a) minimum, using a general-purpose “function-minimization” algorithm (code) such as Powell’s algorithmFootnote 32.

Evaluating the norm (20) requires a numerical integration over the horizon surface: We choose some grid of Nang points on the surface, interpolate the slice geometry fields (gij, kgij, and Kij) to this grid (see Section 7.5), and use numerical quadrature to approximate the integral. In practice this must be done for many different trial surface shapes (see Section 8.3.2), so it is important that it be as efficient as possible. Anninos et al. [8] and Baumgarte et al. [26] discuss various ways to optimize and/or parallelize this calculation.

Unfortunately, minimization algorithms have two serious disadvantages for apparent horizon finding: They can be susceptible to spurious local minima, and they’re very slow at high angular resolutions. However, for the (fairly common) case where we want to find a common apparent horizon as soon as it appears in a binary black-hole (or neutron-star) simulation, minimization algorithms do have a useful ability to “anticipate” the formation of the common apparent horizon, in a manner similar to the pretracking algorithms discussed in Section 8.6. I discuss the properties of minimization algorithms further in the following.

11.3.1 Spurious local minima

While the norm (20) clearly has a single global minimum ‖Θ‖ = 0 for each MOTS Θ = 0, it typically also has a large number of other local minima with Θ ≠ 0, which are “spurious” in the sense that they do not correspond (even approximately) to MOTSsFootnote 33. Unfortunately, generalpurpose “function-minimization” routines only locate local minima and, thus, may easily converge to one of the spurious Θ ≠ 0 minima.

What this problem means in practice is that a minimization algorithm needs quite a good (accurate) initial guess for the horizon shape in order to ensure that the algorithm converges to the true global minimum Θ = 0 rather than to one of the spurious Θ ≠ 0 local minima.

To view this problem from a different perspective, once the function-minimization algorithm does converge, we must somehow determine whether the “solution” found is the true one, Θ = 0, or a spurious one, Θ ≠ 0. Due to numerical errors in the geometry interpolation and the evaluation of the integral (20), ‖Θ‖ will almost never evaluate to exactly zero; rather, we must set a tolerance level for how large ‖Θ‖ may be. Unfortunately, in practice it is hard to choose this tolerance: If it is too small, the genuine solution may be falsely rejected, while if it is too large, we may accept a spurious solution (which may be very different from any of the true solutions).

Anninos et al. [8] and Baumgarte et al. [26] suggest screening out spurious solutions by repeating the algorithm with varying resolutions of the horizon-surface grid and checking that |O| shows the proper convergence towards zero. This seems like a good strategy, but it is tricky to automate and, again, it may be difficult to choose the necessary error tolerances in advance.

When the underlying simulation is a spectral one, Pfeiffer et al. [124, 121] report that in practice, spurious solutions can be avoided by a combination of two factors:

  • The underlying spectral solution can inherently be “interpolated” (evaluated at arbitrary positions) to very high accuracy.

  • Pfeiffer et al. use a large number of quadrature points (typically an order of magnitude larger than the number of coefficients in the expansion (5)) in numerically evaluating the integral (20).

11.3.2 Performance

For convenience of exposition, suppose the spectral representation (5) of the horizon-shape function h uses spherical harmonics Yℓm. (Symmetric trace-free tensors or other basis sets do not change the argument in any important way.) If we keep harmonics up to some maximum degree max, the number of coefficients is then Ncoeff = (max+1)2. max is set by the desired accuracy (angular resolution) of the algorithm and is typically on the order of 6 to 12.

To find a minimum in an Ncoeff-dimensional space (here the space of surface-shape coefficients {aℓm}), a general-purpose function-minimization algorithm typically needs on the order of \(5N_{{\rm{coeff}}}^2\) to \(10N_{{\rm{coeff}}}^2\) iterationsFootnote 34. Thus the number of iterations grows as \(\ell _{\max}^4\).

Each iteration requires an evaluation of the norm (20) for some trial set of surface-shape coefficients {aℓm}, which requires \({\mathcal O}({N_{{\rm{coeff}}}}) = {\mathcal O}(\ell _{\max}^2)\) work to compute the surface positions, together with \({\mathcal O}({N_{{\rm{ang}}}})\) work to interpolate the geometry fields to the surface points and compute the numerical quadrature of the integral (20).

Thus the total work for a single horizon finding is \({\mathcal O}(\ell _{\max}^6 + {N_{{\rm{ang}}}}\ell _{max}^4)\). Fortunately, the accuracy with which the horizon is found generally improves rapidly with max, sometimes even exponentiallyFootnote 35. Thus, relatively modest values of max (typically in the range 8–12) generally suffice for adequate accuracy. Even so, minimization horizon finders tend to be slower than other methods, particularly if high accuracy is required (large max and Nang). The one exception is in axisymmetry, where only spherical harmonics Yℓm with m=0 need be considered. In this case minimization algorithms are much faster, though probably still slower than shooting or elliptic-PDE algorithms.

11.3.3 Anticipating the formation of a common apparent horizon

Consider the case where we want to find a common apparent horizon as soon as it appears in a binary black-hole (or neutron-star) simulation. In Section 8.6 I discuss “horizon pretracking” algorithms which can determine — before it appears — approximately where (in space) and when (in time) the common apparent horizon will appear.

Minimization algorithms can provide a similar functionality: Before the common apparent horizon forms, trying to find it via a minimization algorithm will (hopefully) find the (a) surface which minimizes the error norm ∥Θ∥ (defined by Equation (20)). This surface can be viewed as the current slice’s closest approximation to a common apparent horizon, and as the evolution proceeds, it should converge to the actual common apparent horizon.

However, it is not clear whether minimization algorithms used in this way suffer from the problems discussed in Section 8.6.2. In particular, it is not clear whether, in a realistic binary-coalescence simulation, the minimum-∥Θ∥ surfaces would remain smooth enough to be represented accurately with a reasonable max.

11.3.4 Summary of minimization algorithms/codes

Minimization algorithms are fairly easy to program and have been used by many researchers, for example [43, 69, 102, 8, 26, 4]. However, at least when the underlying simulation uses finite differencing, minimization algorithms are susceptible to spurious local minima, have relatively poor accuracy, and tend to be quite slow. I believe that the other algorithms discussed in the following sections are generally preferable. If the underlying simulation uses spectral methods, then minimization algorithms may be (relatively) somewhat more efficient and robust.

Alcubierre’s apparent horizon finder AHFinder [4] includes a minimization algorithm based on the work of Anninos et al. [8]Footnote 36. It is implemented as a freely available module (“thorn”) in the Cactus computational toolkit (see Table 2).

11.4 Spectral integral-iteration algorithms

Nakamura, Kojima, and Oohara [113] developed a functional-iteration spectral algorithm for solving the apparent horizon equation (16).

This algorithm begins by choosing the usual polar-spherical topology for the angular coordinates (θ,ϕ), and rewriting the apparent horizon equation (16) in the form

$$\Delta h \equiv {\partial _{\theta \theta}}h + {{{\partial _\theta}h} \over {\tan \theta}} + {{{\partial _{\phi \phi}}h} \over {{{\sin}^2}\theta}} = G({\partial _{\theta \phi}}h,{\partial _{\phi \phi}}h,{\partial _\theta}h,{\partial _\phi}h;{g_{ij}},{K_{ij}},\Gamma _{ij}^k),$$
(21)

where Δ is the flat-space angular Laplacian operator, and G is a complicated nonlinear algebraic function of the arguments shown, which remains regular even at θ=0 and θ=π.

Next we expand h in spherical harmonics (5). Because the left hand side of Equation (21) is just the flat-space angular Laplacian of h, which has the Yℓm as orthogonal eigenfunctions, multiplying both sides of Equation (21) by Y*ℓm (the complex conjugate of Yℓm) and integrating over all solid angles gives

$${a_{\ell m}} = - {1 \over {\ell (\ell + 1)}}\int {Y_{\ell m}^{\ast}G\;d\Omega}$$
(22)

for each and m except = m = 0.

Based on this, Nakamura, Kojima, and Oohara [113] proposed the following functional-iteration algorithm for solving Equation (21):

  1. 1.

    Start with some (initial-guess) set of horizon-shape coefficients {aℓm}. These determine a surface shape via Equation (5).

  2. 2.

    Interpolate the geometry variables to this surface shape (see Section 7.5).

  3. 3.

    For each and m except = m = 0, evaluate the integral (22) by numerical quadrature to obtain a next-iteration coefficient am.

  4. 4.

    Determine a next-iteration coefficient a00 by numerically solving (finding a root of) the equation

    $$\int {Y_{00}^{\ast}G\;d\Omega = 0.}$$
    (23)

    This can be done using any of the 1-dimensional zero-finding algorithms discussed in Appendix A.

  5. 5.

    Iterate until all the coefficients {aℓm} converge.

Gundlach [80] observed that the subtraction and inversion of the flat-space angular Laplacian operator in this algorithm is actually a standard technique for solving nonlinear elliptic PDEs by spectral methods. I discuss this observation and its implications further in Section 8.7.4.

Nakamura, Kojima, and Oohara [113] report that their algorithm works well, but Nakao (cited as personal communication in [146]) has argued that it tends to become inefficient (and possibly inaccurate) for large (high angular resolution) because the Yℓm fail to be numerically orthogonal due to the finite resolution of the numerical grid. I know of no other published work addressing Nakao’s criticism.

11.4.1 Kemball and Bishop’s modifications of the Nakamura-Kojima-Oohara algorithm

Kemball and Bishop [93] investigated the behavior of the Nakamura-Kojima-Oohara algorithm and found that its (only) major weakness seems to be that the a00-update equation (23) “may have multiple roots or minima even in the presence of a single marginally outer trapped surface, and all should be tried for convergence”.

Kemball and Bishop [93] suggested and tested several modifications to improve the algorithm’s convergence behavior. They verified that (either in its original form or with their modifications) the algorithm’s rate of convergence (number of iterations to reach a given error level) is roughly independent of the degree max of spherical-harmonic expansion used. They also give an analysis that the algorithm’s cost is \({\mathcal O}(\ell _{\max}^4)\), and its accuracy \(\varepsilon = {\mathcal O}(1/{\ell _{\max}})\), i.e. the cost is \({\mathcal O}(1/{\varepsilon ^4})\). This accuracy is surprisingly low: Exponential convergence with max is typical of spectral algorithms and would be expected here. I do not know of any published work which addresses this discrepancy.

11.4.2 Lin and Novak’s variant of the Nakamura-Kojima-Oohara algorithm

Lin and Novak [104] have developed a variant of the Nakamura-Kojima-Oohara algorithm which avoids the need for a separate search for a00 at each iteration: Write the apparent horizon equation (16) in the form

$$\Delta h - 2h = \lambda \Theta + \Delta h - 2h,$$
(24)

where Δ is again the flat-space angular Laplacian operator and where λ is a nonzero scalar function on the horizon surface. Then choose λ as

$$\lambda = {\left({{{\det {g_{ij}}} \over {\det {f_{ij}}}}} \right)^{1/3}}{[{g^{mn}}({\bar \nabla _m}F)({\bar \nabla _n}F)]^{1/2}}{h^2},$$
(25)

where fij the flat metric of polar spherical coordinates, ∇ is the associated 3-covariant derivative operator, and F is the level set function (7).

Lin and Novak [104] showed that all the spherical-harmonic coefficients aℓm (including a00) can then be found by iteratively solving the equation

$${a_{\ell m}} = - {1 \over {\ell (\ell + 1) + 2}}\int {Y_{\ell m}^{\ast}\lambda (\Theta + \Delta h - 2h)d\Omega}.$$
(26)

Lin and Novak [104] find that this algorithm gives robust convergence and is quite fast, particularly at modest accuracy levels. For example, running on a 2 GHz processor, their implementation takes 3.1, 5.8, 17, 88, and 313 seconds to find the apparent horizon in a test slice to a relative error (measured in the horizon area) of 9 × 10−4, 5 × 10−5, 6 × 10−7, 9 × 10−10, and 3 × 10−12 respectivelyFootnote 37. This implementation is freely available as part of the Lorene toolkit for spectral computations in numerical relativity (see Table 2).

11.4.3 Summary of spectral integral-iteration algorithms

Despite what appears to be fairly good numerical behavior and reasonable ease of implementation, the original Nakamura-Kojima-Oohara algorithm has not been widely used apart from later work by its original developers (see, for example, [115, 114]). Kemball and Bishop [93] have proposed and tested several modifications to the basic Nakamura-Kojima-Oohara algorithm. Lin and Novak [104] have develped a variant of the Nakamura-Kojima-Oohara algorithm which avoids the need for a separate search for the a00 coefficient at each iteration. Their implementation of this variant is freely available as part of the Lorene toolkit for spectral computations in numerical relativity (see Table 2).

11.5 Elliptic-PDE algorithms

The basic concept of elliptic-PDE algorithms is simple: We view the apparent horizon equation (16) as a nonlinear elliptic PDE for the horizon shape function h on the angular-coordinate space and solve this equation by standard finite-differencing techniquesFootnote 38, generally using Newton’s method to solve the resulting set of nonlinear algebraic (finite-difference) equations. Algorithms of this type have been widely used both in axisymmetry and in fully generic slices.

11.5.1 Angular coordinates, grid, and boundary conditions

In more detail, elliptic-PDE algorithms assume that the horizon is a Strahlkörper about some local coordinate origin, and choose an angular coordinate system and a finite-difference grid of Nang points on S2 in the manner discussed in Section 2.2.

The most common choices are the usual polar-spherical coordinates (θ, ϕ) and a uniform “latitude/longitude” grid in these coordinates. Since these coordinates are “unwrapped” relative to the actual S2 trial-horizon-surface topology, the horizon shape function h satisfies periodic boundary conditions across the artificial grid boundary at ϕ = 0 and ϕ = 2π. The north and south poles θ = 0 and θ = π are trickier, but Huq et al. [89, 90], Shibata and Uryū [147], and Schnetter [132, 133]Footnote 39 “reaching across the pole” boundary conditions for these artificial grid boundaries.

Alternatively, Thornburg [156] avoids the z axis coordinate singularity of polar-spherical coordinates by using an “inflated-cube” system of six angular patches to cover S2. Here each patch’s nominal grid is surrounded by a “ghost zone” of additional grid points where h is determined by interpolation from the neighboring patches. The interpatch interpolation thus serves to tie the patches together, enforcing the continuity and differentiability of h across patch boundaries. Thornburg reports that this scheme works well but was quite complicated to program.

Overall, the latitude/longitude grid seems to be the superior choice: it works well, is simple to program, and eases interoperation with other software.

11.5.2 Evaluating the expansion Θ

The next step in the algorithm is to evaluate the expansion Θ given by Equation (16) on the angular grid given a trial horizon surface shape function h on this same grid (6).

Most researchers compute Θ via 2-dimensional angular finite differencing of Equation (16) on the trial horizon surface. 2nd order angular finite differencing is most common, but Thornburg [156] uses 4th order angular finite differencing for increased accuracy.

With a (θ, ϕ) latitude/longitude grid the Θ(h, uh, uvh) function in Equation (16) is singular on the z axis (at the north and south poles θ = 0 and θ = π) but can be regularized by applying L’Hôpital’s rule. Schnetter [132, 133] observes that using a Cartesian basis for all tensors greatly aids in this regularization.

Huq et al. [89, 90] choose, instead, to use a completely different computation technique for Θ, based on 3-dimensional Cartesian finite differencing:

  1. 1.

    They observe that the scalar field F defined by Equation (7) can be evaluated at any (3-dimensional) position in the slice by computing the corresponding (r, θ, ϕ) using the usual flat-space formulas, then interpolating h in the 2-dimensional (θ, ϕ) surface grid.

  2. 2.

    Rewrite the apparent horizon condition (15) in terms of F and its (3-dimensional) Cartesian derivatives,

    $$\Theta \equiv \Theta (F,{\partial _i}F,{\partial _{ij}}F;{g_{ij}},{\partial _k}{g_{ij}},{K_{ij}}) = 0.$$
    (27)
  3. 3.

    For each (latitude/longitude) grid point on the trial horizon surface, define a 3×3×3-point local Cartesian grid centered at that point. The spacing of this grid should be such as to allow accurate finite differencing, i.e. in practice it should probably be roughly comparable to that of the underlying numerical-relativity simulation’s grid.

  4. 4.

    Evaluate F on the local Cartesian grid as described in Step 1 above.

  5. 5.

    Evaluate the Cartesian derivatives in Equation (27) by centered 2nd order Cartesian finite differencing of the F values on the local Cartesian grid.

Comparing the different ways of evaluating Θ, 2-dimensional angular finite differencing of Equation (16) seems to me to be both simpler (easier to program) and likely more efficient than 3-dimensional Cartesian finite differencing of Equation (27).

11.5.3 Solving the nonlinear elliptic PDE

A variety of algorithms are possible for actually solving the nonlinear elliptic PDE (16) (or (27) for the Huq et al. [89, 90] horizon finder).

The most common choice is to use some variant of Newton’s method. That is, finite differencing Equation (16) or (27) (as appropriate) gives a system of Nang nonlinear algebraic equations for the horizon shape function h at the Nang angular grid points; these can be solved by Newton’s method in Nang dimensions. (As explained by Thornburg [153, Section VIII.C], this is usually equivalent to applying the Newton-Kantorovich algorithm [37, Appendix C] to the original nonlinear elliptic PDE (16) or (27).)

Newton’s method converges very quickly once the trial horizon surface is sufficiently close to a solution (a MOTS). However, for a less accurate initial guess, Newton’s method may converge very slowly or even fail to converge at all. There is no usable way of determining a priori just how large the radius of convergence of the iteration will be, but in practice \({1 \over 4} \ {\rm to} \ {1 \over 3}\) of the horizon radius is often a reasonable estimateFootnote 40.

Thornburg [153] described the use of various “line search” modifications to Newton’s method to improve its radius and robustness of convergence, and reported that even fairly simple modifications of this sort roughly doubled the radius of convergence.

Schnetter [132, 133] used the PETSc general-purpose elliptic-solver library [22, 23, 24] to solve the equations. This offers a wide variety of Newton-like algorithms already implemented in a highly optimized form.

Rather than Newton’s method or one of its variants, Shibata et al. [146, 147] use a functional-iteration algorithm directly on the nonlinear elliptic PDE (16). This seems likely to be less efficient than Newton’s method but avoids having to compute and manipulate the Jacobian matrix.

11.5.4 The Jacobian matrix

Newton’s method, and all its variants, require an explicit computation of the Jacobian matrix

$${{\bf{J}}_{{\rm{IJ}}}} = {{\partial {\Theta _{\rm{I}}}} \over {\partial {h_{\rm{J}}}}},$$
(28)

where the indices I and J label angular grid points on the horizon surface (or equivalently on S2).

Notice that J includes contributions both from the direct dependence of Θ on h, uh, and uvh, and also from the indirect dependence of Θ on h through the position-dependence of the geometry variables gij, kgij, and Kij (since Θ depends on the geometry variables at the horizon surface position, and this position is determined by h). Thornburg [153] discusses this indirect dependence in detail.

There are two basic ways to compute the Jacobian matrix.

  • Numerical Perturbation:

    The simplest way to determine the Jacobian matrix is by “numerical perturbation”, where for each horizon-surface grid point j, h is perturbed by some (small) amount ε at the J th grid point (that is, hIhI + εδIJ), and the expansion Θ is recomputedFootnote 41. The J th column of the Jacobian matrix (28) is then estimated as

    $${{\bf{J}}_{{\rm{IJ}}}} \approx {{{\Theta _{\rm{I}}}(h + \varepsilon {\delta _{{\rm{IJ}}}}){\Theta _{\rm{I}}}(h)} \over \varepsilon}.$$
    (29)

    Curtis and Reid [53], and Stoer and Bulirsch [150, Section 5.4.3] discuss the optimum choice of ε in this algorithmFootnote 42.

    This algorithm is easy to program but somewhat inefficient. It is used by a number of researchers including Schnetter [132, 133], and Huq et al. [89, 90].

  • Symbolic Differentiation:

    A more efficient, although somewhat more complicated, way to determine the Jacobian matrix is the “symbolic differentiation” algorithm described by Thornburg [153] and also used by Pasch [118], Shibata et al. [146, 147], and Thornburg [156]. Here the internal structure of the finite differenced Θ(h) function is used to directly determine the Jacobian matrix elements.

    This algorithm is best illustrated by an example which is simpler than the full apparent horizon equation: Consider the flat-space Laplacian in standard (θ, ϕ) polar-spherical coordinates,

    $$\Delta h \equiv {\partial _{\theta \theta}}h + {{{\partial _\theta}h} \over {\tan \theta}} + {{{\partial _{\phi \phi}}h} \over {{{\sin}^2}\theta}}.$$
    (30)

    Suppose we discretize this with centered 2nd order finite differences in θ and ϕ. Then neglecting finite-differencing truncation errors, and temporarily adopting the usual notation for 2-dimensional grid functions, hi, j = h(θ=θi, ϕ=ϕj), our discrete approximation to Δh is given by

    $${(\Delta h)_{i,j}} = {{{h_{i - 1,j}} - 2{h_{i,j}} + {h_{i + 1,j}}} \over {{{(\Delta \theta)}^2}}} + {1 \over {\tan \theta}}{{{h_{i + 1,j}} - {h_{i + 1,j}} - {h_{i - 1,j}}} \over {2\Delta \theta}} + {1 \over {{{\sin}^2}\theta}}{{{h_{i,j - 1}} - 2{h_{i,j}} + {h_{i,j + 1}}} \over {{{(\Delta \phi)}^2}}}.$$
    (31)

    The Jacobian of Δh is thus given by

    $${{\partial {{(\Delta h)}_{(i,j)}}} \over {\partial {h_{(k,\ell)}}}} = \left\{{\begin{array}{*{20}c} {{1 \over {{{(\Delta \theta)}^2}}} \pm {1 \over {2\tan \theta \Delta \theta}}} & {{\rm{if}}\;(k,\ell) = (i \pm 1,j),} \\ {{1 \over {{{\sin}^2}\theta {{(\Delta \phi)}^2}}}} & {{\rm{if}}\;(k,\ell) = (i \pm 1,j),} \\ {- {2 \over {{{(\Delta \theta)}^2}}} - {2 \over {{{\sin}^2}\theta {{(\Delta \phi)}^2}}}} & {{\rm{if}}\;(k,\ell) = (i,j),} \\ 0 & {{\rm{otherwise}}.} \\ \end{array}} \right.$$
    (32)

    Thornburg [153] describes how to generalize this to nonlinear differential operators without having to explicitly manipulate the nonlinear finite difference equations.

11.5.5 Solving the linear equations

All the algorithms described in Section 8.5.3 for treating nonlinear elliptic PDEs require solving a sequence of linear systems of Nang equations in Nang unknowns. Nang is typically on the order of a few thousand, and the Jacobian matrices in question are sparse due to the locality of the angular finite differencing (see Section 8.5.4). Thus, for reasonable efficiency, it is essential to use linear solvers that exploit this sparsity. Unfortunately, many such algorithms/codes only handle symmetric positive-definite matrices while, due to the angular boundary conditionsFootnote 43 (see Section 8.5.1), the Jacobian matrices that arise in apparent horizon finding are generally neither of these.

The numerical solution of large sparse linear systems is a whole subfield of numerical analysis. See, for example, Duff, Erisman, and Reid [65], and Saad [130] for extensive discussionsFootnote 44. In practice, a numerical relativist is unlikely to write her own linear solver but, rather, will use an existing subroutine (library).

Kershaw’s [94] ILUCG iterative solver is often used; this is only moderately efficient, but is quite easy to programFootnote 45. Schnetter [132, 133] reports good results with an ILU-preconditioned GMRES solver from the PETSc library. Thornburg [156] experimented with both an ILUCG solver and a direct sparse LU decomposition solver (Davis’ UMFPACK library [57, 58, 56, 55, 54]), and found each to be more efficient in some situations; overall, he found the UMFPACK solver to be the best choice.

11.5.6 Sample results

As an example of the results obtained with this type of apparent horizon finder, Figure 10 shows the numerically-computed apparent horizons (actually, MOTSs) at two times in a head-on binary black hole collision. (The physical system being simulated here is very similar to that simulated by Matzner et al. [108], a view of whose event horizon is shown in Figure 5.)

Figure 10
figure 10

This figure shows the numerically-computed apparent horizons (actually MOTSs) at two times in a head-on binary black hole collision. The black holes are colliding along the z axis. Figure reprinted with permission from [156]. © 2004 by IOP Publishing Ltd.

As another example, Figure 11 shows the time dependence of the irreducible masses of apparent horizons found in a (spiraling) binary black hole collision, simulated at several different grid resolutions, as found by both AHFinderDirect and another Cactus apparent horizon finder, AHFinderFootnote 46. For this evolution, the two apparent horizon finders give irreducible masses which agree to within about 2% for the individual horizons and 0.5% for the common horizon.

Figure 11
figure 11

This figure shows the irreducible masses \((\sqrt {{\rm{area}}}/(16\pi))\) of individual and common apparent horizons in a binary black hole collision, as calculated by two different apparent horizon finders in the Cactus toolkit, AHFinder and AHFinderDirect. (AHFinderDirect was also run in simulations at several different resolutions.) Notice that when both apparent horizon finders are run in the same simulation (resolution dx = 0.080), there are only small differences between their results. Figure reprinted with permission from [5]. © 2005 by the American Physical Society.

As a final 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 event horizons were computed using the EHFinder code described in Section 5.3.3.)

11.5.7 Summary of elliptic-PDE algorithms/codes

Elliptic-PDE apparent horizon finders have been developed by many researchers, including Eardley [67]Footnote 47, Cook [50, 52, 51], and Thornburg [153] in axisymmetry, and Shibata et al. [146, 147], Huq et al. [89, 90], Schnetter [132, 133], and Thornburg [156] in fully generic slices.

Elliptic-PDE algorithms are (or can be implemented to be) among the fastest horizon finding algorithms. For example, running on a 1.7 GHz processor, Thornburg’s AHFinderDirect [156] averaged 1.7 s per horizon finding, as compared with 61 s for an alternate “fast-flow” apparent horizon finder AHFinder (discussed in more detail in Section 8.7)Footnote 48. However, achieving maximum performance comes at some cost in implementation effort (e.g. the “symbolic differentiation” Jacobian computation discussed in Section 8.5.4).

Elliptic-PDE algorithms are probably somewhat more robust in their convergence (i.e. they have a slightly larger radius of convergence) than other types of local algorithms, particularly if the “line search” modifications of Newton’s method described by Thornburg [153] are implementedFootnote 49. Their typical radius of convergence is on the order of \({1 \over 3}\) of the horizon radius, but cases are known where it is much smaller. For example, Schnetter, Herrmann, and Pollney [135] report that (with no “line search” modifications) it is only about 10% for some slices in a binary black hole coalescence simulation.

Schnetter’s TGRapparentHorizon2D [132, 133] and Thornburg’s AHFinderDirect [156] are both elliptic-PDE apparent horizon finders implemented as freely available modules (“thorns”) in the Cactus computational toolkit (see Table 2). Both work with either the PUGH unigrid driver or the Carpet mesh-refinement driver for Cactus. TGRapparentHorizon2D is no longer maintained, but AHFinderDirect is actively supported and is now used by many different research groupsFootnote 50.

11.6 Horizon pretracking

Schnetter et al. [133, 135] introduced the important concept of “horizon pretracking”. They focus on the case where we want to find a common apparent horizon as soon as it appears in a binary black-hole (or neutron-star) simulation. While a global (flow) algorithm (Section 8.7) could be used to find this common apparent horizon, these algorithms tend to be very slow. They observe that the use of a local (elliptic-PDE) algorithm for this purpose is somewhat problematic:

The common [apparent] horizon […] appears instantaneously at some late time and without a previous good guess for its location. In practice, an estimate of the surface location and shape can be put in by hand. The quality of this guess will determine the rate of convergence of the finder and, more seriously, also determines whether a horizon is found at all. Gauge effects in the strong field region can induce distortions that have a large influence on the shape of the common horizon, making them difficult to predict, particularly after a long evolution using dynamical coordinate conditions. As such, it can be a matter of some expensive trial and error to find the common apparent horizon at the earliest possible time. Further, if a common apparent horizon is not found, it is not clear whether this is because there is none, or whether there exists one which has only been missed due to unsuitable initial guesses — for a fast apparent horizon finder, a good initial guess is crucial.

Pretracking tries (usually successfully) to eliminate these difficulties by determining — before it appears — approximately where (in space) and when (in time) the common apparent horizon will appear.

11.6.1 Constant-expansion surfaces

The basic idea of horizon pretracking is to consider surfaces of constant expansion (“CE surfaces”), i.e. smooth closed orientable 2-surfaces in a slice satisfying the condition

$$\Theta = E,$$
(33)

where the expansion E is a specified real number. Each marginally outer trapped surface (including the apparent horizon) is thus a CE surface with expansion E = 0; more generally Equation (33) defines a 1-parameter family of 2-surfaces in the slice. As discussed by Schnetter et al. [133, 135], for asymptotically flat slices containing a compact strong-field region, some of the E < 0 members of this family typically foliate the weak-field region.

In the binary-coalescence context, for each t = constant slice we define E* to be the smallest E ≥ 0 for which a CE surface (containing both strong-field regions) exists with expansion E. If E* = 0 this “minimum-expansion CE surface” is the common apparent horizon, while if E* > 0 this surface is an approximation to where the common apparent horizon will appear. We expect the minimum-expansion CE surface to change continuously during the evolution and its expansion E* to decrease towards 0. Essentially, horizon pretracking follows the time evolution of the minimum-expansion CE surface and uses it as an initial guess for (searching for) the common apparent horizon.

11.6.2 Generalized constant-expansion surfaces

Schnetter [133] implemented an early form of horizon pretracking, which followed the evolution of the minimum-expansion constant-expansion surface, and found that it worked well for simple test problems. However, Schnetter et al. [135] found that for more realistic binary-black-hole coalescence systems the algorithm needs to be extended:

  • While the expansion is zero for a common apparent horizon, it is also zero for a 2-sphere at spatial infinity. Figure 12 illustrates this for Schwarzschild spacetime. Notice that for small positive E* there will generally be two distinct CE surfaces with E = E*, an inner surface just outside the horizon and an outer one far out in the weak-field region. The inner CE surface converges to the common apparent horizon as E* decreases towards 0; this surface is the one we would like the pretracking algorithm to follow. Unfortunately, without measures such as those described below, there is nothing to prevent the algorithm from following the outer surface, which does not converge to the common apparent horizon as E* decreases towards 0.

  • In a realistic binary-coalescence simulation, the actual minimum-expansion CE surface may be highly distorted and, thus, hard to represent accurately with a finite-resolution angular grid.

Figure 12
figure 12

This figure shows the expansion Θ (left scale), and the “generalized expansions” r Θ (left scale) and r2Θ (right scale), for various r = constant surfaces in an Eddington-Finkelstein slice of Schwarzschild spacetime. Notice that all three functions have zeros at the horizon r = 2m, and that while Θ has a maximum at r ≈ 4.4 m, both r Θ and r2Θ increase monotonically with r.

Schnetter et al. [135] discuss these problems in more detail, arguing that to solve them, the expansion Θ should be generalized to a “shape function” H given by one of

$$\begin{array}{*{20}c} {{H_1} = \Theta,} \\ {\;{H_r} = h\Theta,} \\ {{H_{{r^2}}} = {h^2}\Theta,} \\ \end{array}$$
(34)

CE surfaces are then generalized to surfaces satisfying

$$H = E$$
(35)

for some specified E ≥ 0.

Note that unlike H1, both Hr and Hr2 are typically monotonic with radius. Neither Hr nor Hr2 are 3-covariantly defined, but they both still have the property that E = 0 in Equation (35) implies the surface is a MOTS, and in practice they work better for horizon pretracking.

11.6.3 Goal functions

To define the single “smallest” surface at each time, Schnetter et al. [135] introduce a second generalization, that of a “goal function” G, which maps surfaces to real numbers. The pretracking search then attempts, on each time slice, to find the surface (shape) satisfying H = E with the minimum value of G. They experimented with several different goal functions,

$$\begin{array}{*{20}c} {{G_H} = \overline{H},} \\ {{G_{rH}} = \overline{h} \overline{H},} \\ {\;\;{G_r} = \overline{h},} \\ \end{array}$$
(36)

where in each case the overbar (̅) denotes an average over the surfaceFootnote 51.

11.6.4 The pretracking search

Schnetter’s [133] original implementation of horizon pretracking (which followed the evolution of the minimum-expansion CE surface) used a binary search on the desired expansion E. Because E appears only on the right hand side of the generalized CE condition (35), it is trivial to modify any apparent horizon finder to search for a surface of specified expansion E. (Schnetter used his TGRapparentHorizon2D elliptic-PDE apparent horizon finder described in Section 8.5.7 for this.) A binary search on E can then be used to find the minimum value E*.Footnote 52

Implementing a horizon pretracking search on any of the generalized goal functions (36) is conceptually similar but somewhat more involved: As described by Schnetter et al. [135] for the case of an elliptic-PDE apparent horizon finderFootnote 53, we first write the equation defining a desired pretracking surface as

$$H - \overline{H} + G - p = 0,$$
(37)

where p is the desired value of the goal function G. Since H is the only term in Equation (37) which varies over the surface, it must be constant for the equation to be satisfied. In this case \(H - \bar H\) vanishes, so the equation just gives G = p, as desired.

Because \({\bar H}\) depends on H at all surface points, directly finite differencing Equation (37) would give a non-sparse Jacobian matrix, which would greatly slow the linear-solver phase of the elliptic-PDE apparent horizon finder (Section 8.5.5). Schnetter et al. [135, Section III.B] show how this problem can be solved by introducing a single extra unknown into the discrete system. This gives a Jacobian which has a single non-sparse row and column, but is otherwise sparse, so the linear equations can still be solved efficiently.

When doing the pretracking search, the cost of a single binary-search iteration is approximately the same as that of finding an apparent horizon. Schnetter et al. [135, Figure 5] report that their pretracking implementation (a modified version of Thornburg’s AHFinderDirect [156] elliptic-PDE apparent horizon finder described in Section 8.5.7) typically takes on the order of 5 to 10 binary-search iterationsFootnote 54. The cost of pretracking is thus on the order of 5 to 10 times that of finding a single apparent horizon. This is substantial, but not prohibitive, particularly if the pretracking algorithm is not run at every time step.

11.6.5 Sample results

As an example of the results obtained from horizon pretracking, Figure 13 shows the expansion Θ for various pretracking surfaces (i.e. various choices for the shape function H in a head-on binary black hole collision). Notice how all three of the shape functions (34) result in pretracking surfaces whose expansions converge smoothly to zero just when the apparent horizon appears (at about t = 1.1).

Figure 13
figure 13

This figure shows the expansion 0 for various pretracking surfaces, i.e. for various choices for the shape function H, in a head-on binary black hole collision. Notice how the three shape functions (34) (here labelled Θ, rΘ, and r2Θ) result in pretracking surfaces whose expansions converge smoothly to zero just when the apparent horizon appears (at about t = 1.1). Notice also that these three expansions have all converged to each other somewhat before the common apparent horizon appears. Figure reprinted with permission from [135]. © 2005 by the American Physical Society.

As a further example, Figure 14 shows the pretracking surfaces (more precisely, their cross sections projected into the black holes’ orbital plane) at various times in a spiraling binary black hole collision.

Figure 14
figure 14

This figure shows the pretracking surfaces at various times in a spiraling binary black hole collision, projected into the black holes’ orbital plane. (The apparent slow drift of the black holes in a clockwise direction is an artifact of the corotating coordinate system; the black holes are actually orbiting much faster, in a counterclockwise direction.) Notice how, even well before the common apparent horizon first appears (t = 16.44 mAMD, bottom right plot), the rΘ pretracking surface is already a reasonable approximation to the eventual common apparent horizon’s shape. Figure reprinted with permission from [135]. © 2005 by the American Physical Society.

11.6.6 Summary of horizon pretracking

Pretracking is a very valuable addition to the horizon finding repertoire: It essentially gives a local algorithm (in this case, an elliptic-PDE algorithm) most of the robustness of a global algorithm in terms of finding a common apparent horizon as soon as it appears. It is implemented as a higher-level algorithm which uses a slightly-modified elliptic-PDE apparent horizon finding algorithm as a “subroutine”.

The one significant disadvantage of pretracking is its cost: Each pretracking search typically takes 5 to 10 times as long as finding an apparent horizon. Further research to reduce the cost of pretracking would be useful.

Schnetter et al.’s pretracking implementation [135] is implemented as a set of modifications to Thornburg’s AHFinderDirect [156] apparent horizon finder. Like the original AHFinderDirect, the modified version is a freely available “thorn” in the Cactus toolkit (see Table 2).

11.7 Flow algorithms

Flow algorithms define a “flow” on 2-surfaces, i.e. they define an evolution of 2-surfaces in some pseudo-time λ, such that the apparent horizon is the λ → ∞ limit of a (any) suitable starting surface. Flow algorithms are different from other apparent horizon finding algorithms (except for zero-finding in spherical symmetry) in that their convergence does not depend on having a good initial guess. In other words, flow algorithms are global algorithms (Section 7.7).

To find the (an) apparent horizon, i.e. an outermost MOTS, the starting surface should be outside the largest possible MOTS in the slice. In practice, it generally suffices to start with a 2-sphere of areal radius substantially greater than 2 madm.

The global convergence property requires that a flow algorithm always flow from a large starting surface into the apparent horizon. This means that the algorithm gains no particular benefit from already knowing the approximate position of the apparent horizon. In particular, flow algorithms are no faster when “tracking” the apparent horizon (repeatedly finding it at frequent intervals) in a numerical time evolution. (In contrast, in this situation a local apparent horizon finding algorithm can use the most recent previously-found apparent horizon as an initial guess, greatly speeding the algorithm’s convergenceFootnote 55).

Flow algorithms were first proposed for apparent horizon finding by Tod [157]. He initially considered the case of a time-symmetric slice (one where Kij = 0). In this case, a MOTS (and thus an apparent horizon) is a surface of minimal area and may be found by a “mean curvature flow”

$${\partial _\lambda}{x^i} = - \kappa {s^i},$$
(38)

where xi are the spatial coordinates of a horizon-surface point, si is as before the outward-pointing unit 3-vector normal to the surface, and k ≡ ∇k sk is the mean curvature of the surface as embedded in the slice. This is a gradient flow for the surface area, and Grayson [79] has proven that if the slice contains a minimum-area surface, this will in fact be the stable λ → ∞ limit of this flow. Unfortunately, this proof is valid only for the time-symmetric case.

For non-time-symmetric slices, Tod [157] proposed generalizing the mean curvature flow to the “expansion flow”

$${\partial _\lambda}{x^i} = - \Theta {s^i}.$$
(39)

There is no theoretical proof that this flow will converge to the (an) apparent horizon, but several lines of argument make this convergence plausible:

  • The expansion flow is identical to the mean curvature flow (38) in the principal part.

  • The expansion flow’s velocity is clearly zero on an apparent horizon.

  • More generally, a simple argument due to Bartnik [25]Footnote 56 shows that the expansion flow can never move a (smooth) test surface through an apparent horizon. Suppose, to the contrary, that the test surface \({\mathcal T}\) is about to move through an apparent horizon \({\mathcal H}\), i.e. since both surfaces are by assumption smooth, that \({\mathcal T}\) and \({\mathcal H}\) touch at single (isolated) point P. At that point, \({\mathcal T}\) and \({\mathcal H}\) obviously have the same gij and Kij, and they also have the same si (because P is isolated). Hence the only term in Θ (as defined by Equation (15)) which differs between \({\mathcal T}\) and \({\mathcal H}\) is ∇isi. Clearly, if \({\mathcal T}\) is outside \({\mathcal H}\) and they touch at the single isolated point P, then relative to \({\mathcal H},{\mathcal T}\) must be concave outwards at P, so that \({\nabla _i}{s^i}({\mathcal T}) < {\nabla _i}{s^i}({\mathcal H})\). Thus the expansion flow (39) will move \({\mathcal T}\) outwards, away from the apparent horizon. (If \({\mathcal T}\) lies inside \({\mathcal H}\) the same argument holds with signs reversed appropriately.)

Numerical experiments by Bernstein [28], Shoemaker et al. [148, 149], and Pasch [118] show that in practice the expansion flow (39) does in fact converge robustly to the apparent horizon.

In the following I discuss a number of important implementation details for, and refinements of, this basic algorithm.

11.7.1 Implicit pseudo-time stepping

Assuming the Strahlkörper surface parameterization (4), the expansion flow (39) is a parabolic equation for the horizon shape function h.Footnote 57 This means that any fully explicit scheme to integrate it (in the pseudo-time λ) must severely restrict its pseudo-time step Δλ for stability, and this restriction grows (quadratically) worse at higher spatial resolutionsFootnote 58. This makes the horizon finding process very slow.

To avoid this restriction, practical implementations of flow algorithms use implicit pseudo-time integration schemes; these can have large pseudo-time steps and still be stable. Because we only care about the λ → ∞ limit, a highly accurate pseudo-time integration is not important; only the accuracy of approximating the spatial derivatives matters. Bernstein [28] used a modified Du Fort-Frankel scheme [64]Footnote 59 but found some problems with the surface shape gradually developing high-spatial-frequency noise. Pasch [118] reports that an “exponential” integrator (Hochbrucket al. [85]) works well, provided the flow’s Jacobian matrix is computed accuratelyFootnote 60. The most common choice is probably that of Shoemaker et al. [148, 149], who use the iterated Crank-Nicholson (“ICN”) schemeFootnote 61. They report that this works very well; in particular, they do not report any noise problems.

By refining his finite-element grid (Section 2.3) in a hierarchical manner, Metzger [109] is able to use standard conjugate-gradient elliptic solvers in a multigrid-like fashionFootnote 62, using each refinement level’s solution as an initial guess for the next higher refinement level’s iterative solution. This greatly speeds the flow integration: Metzger reports that the performance of the overall surface-finding algorithm is “of the same order of magnitude” as that of Thornburg’s AHFinderDirect [156] elliptic-PDE apparent horizon finder (described in Section 8.5.7).

In a more general context than numerical relativity, Osher and Sethian [116] have discussed a general class of numerical algorithms for integrating “fronts propagating with curvature-dependent speed”. These flow a level-set function (Section 2.1) which implicitly locates the actual “front”.

11.7.2 Varying the flow speed

Another important performance optimization of the standard expansion flow (39) is to replace Θ in the right-hand side by a suitable nonlinear function of Θ, chosen so the surface shrinks faster when it is far from the apparent horizon. For example, Shoemaker et al. [148, 149] use the flow

$${\partial _\lambda}{x^i} = - \left[ {(\Theta - c){{\arctan}^2}\left({{{\Theta - c} \over {{\Theta _0}}}} \right)} \right]{s^i}$$
(40)

for this purpose, where Θ0 is the value of Θ on the initial-guess surface, and c (which is gradually decreased towards 0 as the iteration proceeds) is a “goal” value for Θ.

11.7.3 Surface representation and the handling of bifurcations

Since a flow algorithm starts with (topologically) a single large 2-sphere, if there are multiple apparent horizons present the surface must change topology (bifurcate) at some point in the flow. Depending on how the surface is represented, this may be easy or difficult.

Pasch [118] and Shoemaker et al. [148, 149] use a level-set function approach (Section 2.1). This automatically handles any topology or topology change. However, it has the drawback of requiring the flow to be integrated throughout the entire volume of the slice (or at least in some neighborhood of each surface). This is likely to be much more expensive than only integrating the flow on the surface itself. Shoemaker et al. also generate an explicit Strahlkörper surface representation (Section 2.2), monitoring the surface shape to detect an imminent bifurcation and reparameterizing the shape into 2 separate surfaces if a bifurcation happens.

Metzger [109] uses a finite-element surface representation (Section 2.3), which can represent any topology. However, if the flow bifurcates, then to explicitly represent each apparent horizon the code must detect that the surface self-intersects, which may be expensive.

11.7.4 Gundlach’s “fast flow”

Gundlach [80] introduced the important concept of a “fast flow”. He observed that the subtraction and inversion of the flat-space Laplacian in the Nakamura-Kojima-Oohara spectral integral-iteration algorithm (Section 8.4) is an example of “a standard way of solving nonlinear elliptic problems numerically, namely subtracting a simple linear elliptic operator from the nonlinear one, inverting it by pseudo-spectral algorithms and iterating”. Gundlach then interpreted the Nakamura-Kojima-Oohara algorithm as a type of flow algorithm where each pseudo-time step of the flow corresponds to a single functional-iteration step of the Nakamura-Kojima-Oohara algorithm.

In this framework, Gundlach defines a 2-parameter family of flows interpolating between the Nakamura-Kojima-Oohara algorithm and Tod’s [157] expansion flow (39),

$${\partial _\lambda}h = - A{(1 - B\Delta)^{- 1}}\rho \Theta,$$
(41)

where A ≥ 0 and B ≥ 0 are parameters, ρ > 0 is a weight functional which depends on h through at most 1st derivatives, Δ is the flat-space Laplacian operator, and (1 — BΔ)−1 denotes inverting the operator (1 — BΔ). Gundlach then argues that intermediate “fast flow” members of this family should be useful compromises between the speed of the Nakamura-Kojima-Oohara algorithm and the robustness of Tod’s expansion flow. Based on numerical experiments, Gundlach suggests a particular choice for the weight functional ρ and the parameters A and B. The resulting algorithm updates high-spatial-frequency components of h essentially the same as the Nakamura-Kojima-Oohara algorithm but should reduce low-spatial-frequency error components faster. Gundlach’s algorithm also completely avoids the need for numerically solving Equation (23) for the a00 coefficient in the Nakamura-Kojima-Oohara algorithm.

Alcubierre’s AHFinder [4] horizon finder includes an implementation of Gundlach’s fast flow algorithm63. AHFinder is implemented as a freely available module (“thorn”) in the Cactus computational toolkit (see Table 2) and has been used by many research groups.

11.7.5 Summary of flow algorithms/codes

Flow algorithms are the only truly global apparent horizon finding algorithms and, as such, can be much more robust than local algorithms. In particular, flow algorithms can guarantee convergence to the outermost MOTS in a slice. Unfortunately, these convergence guarantees hold only for time-symmetric slices.

In the forms which have strong convergence guarantees, flow algorithms tend to be very slow. (Metzger’s algorithm [109] is a notable exception: It is very fast.) There are modifications which can make flow algorithms much faster, but then their convergence is no longer guaranteed. In particular, practical experience has shown that in some binary black hole coalescence simulations (Alcubierre et al. [5], Diener et al. [62]), “fast flow” algorithms (Section 8.7.4) can miss common apparent horizons which are found by other (local) algorithms.

Alcubierre’s apparent horizon finder AHFinder [4] includes a “fast flow” algorithm based on the work of Gundlach [80]Footnote 63. It is implemented as a freely available module (“thorn”) in the Cactus computational toolkit (see Table 2) and has been used by a number of research groups.

12 Summary of Algorithms/Codes for Finding Apparent Horizons

There are many apparent horizon finding algorithms, with differing trade-offs between speed, robustness of convergence, accuracy, and ease of programming.

In spherical symmetry, zero-finding (Section 8.1) is fast, robust, and easy to program. In axisymmetry, shooting algorithms (Section 8.2) work well and are fairly easy to program. Alternatively, any of the algorithms for generic slices (summarized below) can be used with implementations tailored to the axisymmetry.

Minimization algorithms (Section 8.3) are fairly easy to program, but when the underlying simulation uses finite differencing these algorithms are susceptible to spurious local minima, have relatively poor accuracy, and tend to be very slow unless axisymmetry is assumed. When the underlying simulation uses spectral methods, then minimization algorithms can be somewhat faster and more robust.

Spectral integral-iteration algorithms (Section 8.4) and elliptic-PDE algorithms (Section 8.5) are fast and accurate, but are moderately difficult to program. Their main disadvantage is the need for a fairly good initial guess for the horizon position/shape.

In many cases Schnetter’s “pretracking” algorithm (Section 8.6) can greatly improve an elliptic-PDE algorithm’s robustness, by determining — before it appears — approximately where (in space) and when (in time) a new outermost apparent horizon will appear. Pretracking is implemented as a modification of an existing elliptic-PDE algorithm and is moderately slow: It typically has a cost 5 to 10 times that of finding a single horizon with the elliptic-PDE algorithm.

Finally, flow algorithms (Section 8.7) are generally quite slow (Metzger’s algorithm [109] is a notable exception) but can be very robust in their convergence. They are moderately easy to program. Flow algorithms are global algorithms, in that their convergence does not depend on having a good initial guess.

Table 2 lists freely-available apparent horizon finding codes.