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 points on 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 trial-horizon-surface topology, the horizon shape function satisfies periodic boundary conditions across the artificial grid boundary at and . The north and south poles and are trickier, but Huq et al. [89, 90], Shibata and Uryū , and Schnetter [132, 133]39 “reaching across the pole” boundary conditions for these artificial grid boundaries.
Alternatively, Thornburg  avoids the axis coordinate singularity of polar-spherical coordinates by using an “inflated-cube” system of six angular patches to cover . Here each patch’s nominal grid is surrounded by a “ghost zone” of additional grid points where is determined by interpolation from the neighboring patches. The interpatch interpolation thus serves to tie the patches together, enforcing the continuity and differentiability of 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.
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 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  uses 4th order angular finite differencing for increased accuracy.
With a latitude/longitude grid the function in Equation (16) is singular on the axis (at the north and south poles 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:
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).
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 nonlinear algebraic equations for the horizon shape function at the angular grid points; these can be solved by Newton’s method in 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 to of the horizon radius is often a reasonable estimate40.
Thornburg  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.
Newton’s method, and all its variants, require an explicit computation of the Jacobian matrix
Notice that includes contributions both from the direct dependence of on , , and , and also from the indirect dependence of on through the position-dependence of the geometry variables , , and (since depends on the geometry variables at the horizon surface position, and this position is determined by ). Thornburg  discusses this indirect dependence in detail.
There are two basic ways to compute the Jacobian matrix.
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].
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, describes how to generalize this to nonlinear differential operators without having to explicitly manipulate the nonlinear finite difference equations.
All the algorithms described in Section 8.5.3 for treating nonlinear elliptic PDEs require solving a sequence of linear systems of equations in unknowns. 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 conditions43 (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 , and Saad  for extensive discussions44. In practice, a numerical relativist is unlikely to write her own linear solver but, rather, will use an existing subroutine (library).
Kershaw’s  ILUCG iterative solver is often used; this is only moderately efficient, but is quite easy to program45. Schnetter [132, 133] reports good results with an ILU-preconditioned GMRES solver from the PETSc library. Thornburg  experimented with both an ILUCG solver and a direct sparse 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.
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. , a view of whose event horizon is shown in Figure 5.)
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, AHFinder46. 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.
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.)
Elliptic-PDE apparent horizon finders have been developed by many researchers, including Eardley 47, Cook [50, 52, 51], and Thornburg  in axisymmetry, and Shibata et al. [146, 147], Huq et al. [89, 90], Schnetter [132, 133], and Thornburg  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  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)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  are implemented49. Their typical radius of convergence is on the order of of the horizon radius, but cases are known where it is much smaller. For example, Schnetter, Herrmann, and Pollney  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  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 groups50.
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License.