In more detail, ellipticPDE algorithms assume that the horizon is a Strahlkörper about some local coordinate origin, and choose an angular coordinate system and a finitedifference grid of points on in the manner discussed in Section 2.2.
The most common choices are the usual polarspherical coordinates and a uniform “latitude/longitude” grid in these coordinates. Since these coordinates are “unwrapped” relative to the actual trialhorizonsurface 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Ć« [147], and Schnetter [132, 133]^{39} “reaching across the pole” boundary conditions for these artificial grid boundaries.
Alternatively, Thornburg [156] avoids the axis coordinate singularity of polarspherical coordinates by using an “inflatedcube” 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 2dimensional 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 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 3dimensional Cartesian finite differencing:
Comparing the different ways of evaluating , 2dimensional angular finite differencing of Equation (16) seems to me to be both simpler (easier to program) and likely more efficient than 3dimensional 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 estimate^{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 generalpurpose ellipticsolver library [22, 23, 24] to solve the equations. This offers a wide variety of Newtonlike 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 functionaliteration 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
where the indices and label angular grid points on the horizon surface (or equivalently on ).Notice that includes contributions both from the direct dependence of on , , and , and also from the indirect dependence of on through the positiondependence of the geometry variables , , and (since depends on the geometry variables at the horizon surface position, and this position is determined by ). Thornburg [153] 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 flatspace Laplacian in standard polarspherical coordinates,
Suppose we discretize this with centered 2nd order finite differences in and . Then neglecting finitedifferencing truncation errors, and temporarily adopting the usual notation for 2dimensional grid functions, , our discrete approximation to is given by The Jacobian of is thus given by Thornburg [153] 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 positivedefinite matrices while, due to the angular boundary conditions^{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 discussions^{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 program^{45}. Schnetter [132, 133] reports good results with an ILUpreconditioned GMRES solver from the PETSc library. Thornburg [156] 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 numericallycomputed apparent horizons (actually, MOTSs) at two times in a headon 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.)

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, AHFinder^{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.

As a final example, Figure 8 shows the numericallycomputed 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.)
EllipticPDE apparent horizon finders have been developed by many researchers, including Eardley [67]^{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.
EllipticPDE 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 “fastflow” 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).
EllipticPDE 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 implemented^{49}. 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 [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 ellipticPDE 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 meshrefinement driver for Cactus. TGRapparentHorizon2D is no longer maintained, but AHFinderDirect is actively supported and is now used by many different research groups^{50}.
http://www.livingreviews.org/lrr20073 
This work is licensed under a Creative Commons AttributionNoncommercialNo Derivative Works 2.0 Germany License. Problems/comments to 