List of Footnotes

1 An algorithm’s actual “convergence region” (the set of all initial guesses for which the algorithm converges to the correct solution) may even be fractal in shape. For example, the Julia set is the convergence region of Newton’s method on a simple nonlinear algebraic equation.
2 For convenience of exposition I use spherical harmonics here, but there are no essential differences if other basis sets are used.
3 I discuss the choice of this angular grid in more detail in Section 8.5.1.
4 There has been some controversy over whether, and if so how quickly, Regge calculus converges to the continuum Einstein equations. (See, for example, the debate between Brewin [40] and Miller [110], and the explicit convergence demonstration of Gentle and Miller [73].) However, Brewin and Gentle [41] seem to have resolved this: Regge calculus does, in fact, converge to the continuum solution, and this convergence is generically 2nd order in the resolution.
5 See, for example, Choptuik [48], Pretorius [127], Schnetter et al. [134Jump To The Next Citation Point], and Pretorius and Choptuik [126] for general surveys of the uses of, and methods for, mesh refinement in numerical relativity.
6 Chruściel and Galloway [49] showed that if a “cloud of sand” falls into a large black hole, each “sand grain” generates a non-differentiable caustic in the event horizon.
7 This is a statement about the types of spacetimes usually studied by numerical relativists, not a statement about the mathematical properties of the event horizon itself.
8 I briefly review ODE integration algorithms and codes in Appendix B.
9 In practice the differentiation can usefully be combined with the interpolation; I outline how this can be done in Section 7.5.
10 This convergence is only true in a global sense: locally the event horizon has no special geometric properties, and the Riemann tensor components which govern geodesic convergence/divergence may have either sign.
11 Diener [60Jump To The Next Citation Point] describes how the algorithm can be enhanced to also determine (integrate) individual null generators of the event horizon. This requires interpolating the 4-metric to the generator positions but (still) not taking any derivatives of the 4-metric.
12 Walker [162Jump To The Next Citation Point] mentions an implementation for fully generic slices but only presents results for the axisymmetric case.
13 See [7, 103Jump To The Next Citation Point, 162].
14 Equivalently, Diener [60Jump To The Next Citation Point] observed that the locus of any given nonzero value of the level-set function F is itself a null surface and tends to move (exponentially) closer and closer to the event horizon as the backwards evolution proceeds.
15 They describe how Richardson extrapolation can be used to improve the accuracy of the solutions from 𝒪 (ɛ) to 𝒪(ɛ2), but it appears that this has not been done for their published results.
16 Note that the surface must be smooth everywhere. If this condition were not imposed, then MOTSs would lose many of their important properties. For example, even a standard t = constant slice of Minkowski spacetime contains many non-smooth “MOTSs”: The surface of any closed polyhedron in such a slice satisfies all the other conditions to be an MOTS.
17 Andersson and Metzger [6] have shown that MOTSs can only intersect if they are contained within an outer common MOTS. Szilágyi et al. [151] give a numerical example of such overlapping MOTSs found in a binary black hole collision.
18 As an indication of the importance of the “closed” requirement, Hawking [81] observed that if we consider two spacelike-separated events in Minkowski spacetime, the intersection of their backwards light cones satisfies all the conditions of the MOTS definition, except that it is not closed.
19 The proof is given by Hawking and Ellis [82, Proposition 9.2.8] and by Wald [160, Propositions 12.2.3 and 12.2.4].
20 Wald and Iyer [161] proved this by explicitly constructing a family of angularly anisotropic slices in Schwarzschild spacetime which approach arbitrarily close to r = 0 yet contain no apparent horizons. However, Schnetter and Krishnan [136] have recently studied the behavior of apparent horizons in various anisotropic slices in Schwarzschild and Vaidya spacetimes, finding that the Wald and Iyer behavior seems to be rare.
21 This world-tube is sometimes called “the apparent horizon”, but this is not standard terminology. In this review I always use the terminology that an MOTS or apparent horizon is a 2-surface contained in a (single) slice.
22 Ashtekar and Galloway [17] have recently proved “a number of physically interesting constraints” on this slicing-dependence.
23 The derivation of this condition is given by (for example) York [164], Gundlach [80Jump To The Next Citation Point, Section IIA], and Baumgarte and Shapiro [27, Section 6.1].
24 Notice that in order for the 3-divergence in Equation (15View Equation) to be meaningful, i s (defined only as a field on the MOTS) must be smoothly continued off the surface and extended to a field in some 3-dimensional neighborhood of the surface. The off-surface continuation is non-unique, but it is easy to see that this does not affect the value of Θ on the surface.
25 Or, in the Huq et al. [89Jump To The Next Citation Point, 90Jump To The Next Citation Point] algorithm described in Section 8.5.2, at the local Cartesian grid point positions.
26 If the underlying simulation uses spectral methods then the spectral series can be evaluated anywhere, so no actual interpolation need be done, although the term “spectral interpolation” is still often used. See Fornberg [70Jump To The Next Citation Point], Gottlieb and Orszag [75Jump To The Next Citation Point], and Boyd [37Jump To The Next Citation Point] for general discussions of spectral methods, and (for example) Ansorg et al. [12Jump To The Next Citation Point, 11Jump To The Next Citation Point, 10Jump To The Next Citation Point, 13Jump To The Next Citation Point], Bonazzola et al. [35Jump To The Next Citation Point, 33Jump To The Next Citation Point, 34Jump To The Next Citation Point], Grandclément et al. [77Jump To The Next Citation Point], Kidder et al. [95Jump To The Next Citation Point, 96Jump To The Next Citation Point, 97Jump To The Next Citation Point], and Pfeiffer et al. [120Jump To The Next Citation Point, 124Jump To The Next Citation Point, 123Jump To The Next Citation Point, 122Jump To The Next Citation Point] for applications to numerical relativity.
27 Conceptually, an interpolator generally works by locally fitting a fitting function (usually a low-degree polynomial) to the data points in a neighborhood of the interpolation point, then evaluating the fitting function at the interpolation point. By evaluating the derivative of the fitting function, the ∂kgij values can be obtained very cheaply at the same time as the gij values.
28 Thornburg [154Jump To The Next Citation Point, Appendix F] gives a more detailed discussion of the non-smoothness of Lagrange-polynomial interpolation errors.
29 Note that ∂rgθθ is a known coefficient field here, not an unknown; if necessary, it can be obtained by numerically differentiating gθθ. Therefore, despite the appearance of the derivative, Equation (17View Equation) is still an algebraic equation for the horizon radius h, not a differential equation.
30 See also the work of Bizoń, Malec, and Ó Murchadha [32] for an interesting analytical study giving necessary and sufficient conditions for apparent horizons to form in non-vacuum spherically symmetric spacetimes.
31 Ascher, Mattheij, and Russel [15, Chapter 4] give a more detailed discussion of shooting methods.
32 See, for example, Dennis and Schnabel [59], or Brent [39] for general surveys of general-purposes function-minimization algorithms and codes.
33 There is a simple heuristic argument (see, for example, Press et al. [125Jump To The Next Citation Point, Section 9.6]) that at least some spurious local minima should be expected. We are trying to solve a system of Nang nonlinear equations, Θi = 0 (one equation for each horizon-surface grid point). Equivalently, we are trying to find the intersection of the Nang codimension-one hypersurfaces Θi = 0 in surface-shape space. The problem is that anywhere two or more of these hypersurfaces approach close to each other, but do not actually intersect, there is a spurious local minimum in ∥Θ∥.
34 A simple counting argument suffices to show that any general-purpose function-minimization algorithm in n dimensions must involve at least 𝒪 (n2) function evaluations (see, for example, Press et al. [125Jump To The Next Citation Point, Section 10.6]): Suppose the function to be minimized is f : ℜn → ℜ, and suppose f has a local minimum near some point x0 ∈ ℜn. Taylor-expanding f in a neighborhood of x0 gives f(x) = f(x0)+ aT(x− x0)+ (x− x0)TB(x− x0)+ 𝒪 (∥x− x0∥3), where a ∈ ℜn, B ∈ ℜn×n is symmetric, and vT denotes the transpose of the column vector v ∈ ℜn. Neglecting the higher order terms (i.e. approximating f as a quadratic form in x in a neighborhood of x0), and ignoring f(x0) (which does not affect the position of the minimum), there are a total of N = n+ 12n (n + 1) coefficients in this expression. Changing any of these coefficients may change the position of the minimum, and at each function evaluation the algorithm “learns” only a single number (the value of f at the selected evaluation point), so the algorithm must make at least N = 𝒪(n2) function evaluations to (implicitly) determine all the coefficients. Actual functions are not exact quadratic forms, so in practice there are additional 𝒪 (1) multiplicative factors in the number of function evaluations. Minimization algorithms may also make additional performance and/or space-versus-time trade-offs to improve numerical robustness or to avoid explicitly manipulating n × n Jacobian matrices.
35 In the context of an underlying simulation with spectral accuracy, Pfeiffer [122Jump To The Next Citation Point, 121] reports exponential convergence of the horizon finding accuracy with ℓmax.
36 AHFinder also includes a “fast flow” algorithm (Section 8.7).
37 For comparison, the elliptic-PDE AHFinderDirect horizon finder (discussed in Section 8.5.6), running on a roughly similar processor, takes about 1.8 seconds to find the apparent horizon in a similar test slice to a relative error of −4 4× 10.
38 In theory this equation could also be solved by a spectral method on S2, using spectral differentiation to evaluate the angular derivatives. (See [70, 75, 37Jump To The Next Citation Point, 12, 11, 10, 13, 35, 33, 34, 77, 95, 96, 97, 120, 124, 123, 122] for further discussion of spectral methods.) This should yield a highly efficient apparent horizon finder. However, I know of no published work taking this approach.
39 See [133Jump To The Next Citation Point] for a substantially revised version of [132Jump To The Next Citation Point].
40 Thornburg [153Jump To The Next Citation Point] used a Monte-Carlo survey of horizon-shape perturbations to quantify the radius of convergence of Newton’s method for apparent horizon finding. He found that if strong high-spatial-frequency perturbations are present in the slice’s geometry then the radius of convergence may be very small. Fortunately, this problem rarely occurs in practice.
41 A very important optimization here is that Θ only needs to be recomputed within the finite-difference domain of dependence of the j th grid point.
42 Because of the one-sided finite differencing, the approximation (29View Equation) is only 𝒪 (ɛ) accurate. However, in practice this does not seriously impair the convergence of a horizon finder, and the extra cost of a centered-finite-differencing 2 𝒪 (ɛ ) approximation is not warranted.
43 Or the interpatch interpolation conditions in Thornburg’s multiple-grid-patch scheme [156Jump To The Next Citation Point].
44 Multigrid algorithms are also important here; these exploit the geometric structure of the underlying elliptic PDE. See Briggs, Henson, and McCormick [42Jump To The Next Citation Point] and Trottenberg, Oosterlee, and Schüller [158Jump To The Next Citation Point] for general introductions to multigrid algorithms.
45 Madderom’s Fortran subroutine DILUCG [107], which implements the method of [94], has been used by a number of numerical relativists for both this and other purposes.
46 AHFinder incorporates both a minimization algorithm (Section 8.3) and a fast-flow algorithm (Section 8.7.4); these tests used the fast-flow algorithm.
47 This paper does not say how the author finds apparent horizons, but [68, page 135] cites a preprint of this as treating the apparent-horizon equation as a 2-point (ODE) boundary value problem: Eardley uses a ‘beads on a string’ technique to solve the set of simultaneous equations, i.e., imagining the curve to be defined as a bead on each ray of constant angle. He solves for the positions on each ray at which the relation is satisfied everywhere.
48 As another comparison, the Lorene apparent horizon finder (discussed in more detail in Section 8.4.2), running on a roughly similar processor, takes between 3 and 6 seconds to find apparent horizons to comparable accuracy.
49 The convergence problems, which Thornburg [153] noted when high-spatial-frequency perturbations are present in the slice’s geometry, seem to be rare in practice.
50 In addition, at least two different research groups have now ported, or are in the process of porting, AHFinderDirect to their own (non-External LinkCactus) numerical relativity codes.
51 Schnetter et al. [135Jump To The Next Citation Point] use a simple arithmetic mean over all surface grid points. In theory this average could be defined 3-covariantly by taking the induced metric on the surface into account, but in practice they found that this was not worth the added complexity.
52 There is one complication here: Any local apparent horizon finding algorithm may fail if the initial guess is not good enough, even if the desired surface is actually present. The solution is to use the constant-expansion surface for a slightly larger expansion E as an initial guess, gradually “walking down” the value of E to find the minimum value E∗. Thornburg [156Jump To The Next Citation Point, Appendix C] describes such a “continuation-algorithm binary search” algorithm in detail.
53 As far as I know this is the only case that has been considered for horizon pretracking. Extension to other types of apparent horizon finders might be a fruitful area for further research.
54 This refers to the period before a common apparent horizon is found. Once a common apparent horizon is found, then pretracking can be disabled as the apparent horizon finder can easily “track” the apparent horizon’s motion from one time step to the next. With a binary search the number of iterations depends only weakly (logarithmically) on the pretracking algorithm’s accuracy tolerance. It might be possible to replace the binary search by a more sophisticated 1-dimensional search algorithm (I discuss such algorithms in Appendix A), potentially cutting the number of iterations substantially. This might be a fruitful area for further research.
55 Alternatively, a flow algorithm could use the most recent previously-found apparent horizon as an initial guess. In this case the algorithm would have only local convergence (in particular, it would probably fail to find a new outermost MOTS that appeared well outside the previously-found MOTS). However, the algorithm would only need to flow the surface a small distance, so the algorithm should be fairly fast.
56 Cited as Ref. [17] by [80Jump To The Next Citation Point].
57 Linearizing the Θ (h) function (16View Equation) gives a negative Laplacian in h as the principal part.
58 For a spatial resolution Δx, an explicit scheme is generally limited to a pseudo-time step Δλ ≲ (Δx)2.
59 Richtmyer and Morton [129, Section 7.5] give a very clear presentation and analysis of the Du Fort–Frankel scheme.
60 More precisely, Pasch [118Jump To The Next Citation Point] found that that an exponential integrator worked well when the flow’s Jacobian matrix was computed exactly (using the symbolic-differentiation technique described in Section 8.5.4). However, when the Jacobian matrix was approximated using the numerical-perturbation technique described in Section 8.5.4, Pasch found that the pseudo-time integration became unstable at high numerical resolutions. Pasch [118Jump To The Next Citation Point] also notes that the exponential integrator uses a very large amount of memory.
61 Teukolsky [152], and Leiler and Rezzolla [101] have analyzed ICN’s stability under various conditions.
62 See [42, 158] for general introductions to multigrid algorithms for elliptic PDEs.
63 AHFinder also includes a minimization algorithm (Section 8.3).
64 The parabola generically has two roots, but normally only one of them lies between x− and x+.
65 The numerical-analysis literature usually refers to this as the “initial value problem”. Unfortunately, in a relativity context this terminology often causes confusion with the “initial data problem” of solving the ADM constraint equations. I use the term “time-integration problem for ODEs” to (try to) avoid this confusion. In this appendix, sans-serif lower-case letters abc...z denote variables and functions in ℜn (for some fixed dimension n), and sans-serif upper-case letters ABC ...Z denote n× n real-valued matrices.
66 LSODA can also automatically detect stiff systems of ODEs and adjust its integration scheme so as to handle them efficiently.