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. [134], 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 [60] 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 [162] mentions an implementation for fully generic slices but only presents results for the axisymmetric case. | |

13 | See [7, 103, 162]. | |

14 | Equivalently, Diener [60] observed that the locus of any given nonzero value of the level-set function 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 , 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 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 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 [80, Section IIA], and Baumgarte and Shapiro [27, Section 6.1]. | |

24 | Notice that in order for the 3-divergence in Equation (15) to be meaningful, (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. [89, 90] 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 [70], Gottlieb and Orszag [75], and Boyd [37] for general discussions of spectral methods, and (for example) Ansorg et al. [12, 11, 10, 13], Bonazzola et al. [35, 33, 34], Grandclément et al. [77], Kidder et al. [95, 96, 97], and Pfeiffer et al. [120, 124, 123, 122] 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 values can be obtained very cheaply at the same time as the values. | |

28 | Thornburg [154, Appendix F] gives a more detailed discussion of the non-smoothness of Lagrange-polynomial interpolation errors. | |

29 | Note that is a known coefficient field here, not an unknown; if necessary, it can be obtained by numerically differentiating . Therefore, despite the appearance of the derivative, Equation (17) is still an algebraic equation for the horizon radius , 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. [125, Section 9.6]) that at least some spurious local minima should be expected. We are trying to solve a system of nonlinear equations, (one equation for each horizon-surface grid point). Equivalently, we are trying to find the intersection of the codimension-one hypersurfaces 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 dimensions must involve at least function evaluations (see, for example, Press et al. [125, Section 10.6]): Suppose the function to be minimized is , and suppose has a local minimum near some point . Taylor-expanding in a neighborhood of gives , where , is symmetric, and denotes the transpose of the column vector . Neglecting the higher order terms (i.e. approximating as a quadratic form in in a neighborhood of ), and ignoring (which does not affect the position of the minimum), there are a total of 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 at the selected evaluation point), so the algorithm must make at least function evaluations to (implicitly) determine all the coefficients. Actual functions are not exact quadratic forms, so in practice there are additional 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 Jacobian matrices. | |

35 | In the context of an underlying simulation with spectral accuracy, Pfeiffer [122, 121] reports exponential convergence of the horizon finding accuracy with . | |

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 seconds to find the apparent horizon in a similar test slice to a relative error of . | |

38 | In theory this equation could also be solved by a spectral method on , using spectral differentiation to evaluate the angular derivatives. (See [70, 75, 37, 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 [133] for a substantially revised version of [132]. | |

40 | Thornburg [153] 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 th grid point. | |

42 | Because of the one-sided finite differencing, the approximation (29) 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 approximation is not warranted. | |

43 | Or the interpatch interpolation conditions in Thornburg’s multiple-grid-patch scheme [156]. | |

44 | Multigrid algorithms are also important here; these exploit the geometric structure of the underlying elliptic PDE. See Briggs, Henson, and McCormick [42] and Trottenberg, Oosterlee, and Schüller [158] 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-Cactus) numerical relativity codes. | |

51 | Schnetter et al. [135] 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 as an initial guess, gradually “walking down” the value of to find the minimum value . Thornburg [156, 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 [80]. | |

57 | Linearizing the function (16) gives a negative Laplacian in as the principal part. | |

58 | For a spatial resolution , an explicit scheme is generally limited to a pseudo-time step . | |

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 [118] 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 [118] 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 and . | |

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 denote variables and functions in (for some fixed dimension ), and sans-serif upper-case letters denote real-valued matrices. | |

66 | LSODA can also automatically detect stiff systems of ODEs and adjust its integration scheme so as to handle them efficiently. |

http://www.livingreviews.org/lrr-2007-3 |
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License. Problems/comments to |