### 8.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 on the angular-coordinate space and solve this equation by standard finite-differencing techniques, 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.

#### 8.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 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. [8990], Shibata and UryĆ« [147], and Schnetter [132133] “reaching across the pole” boundary conditions for these artificial grid boundaries.

Alternatively, Thornburg [156] 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.

#### 8.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 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 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 [132133] observes that using a Cartesian basis for all tensors greatly aids in this regularization.

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

1. They observe that the scalar field defined by Equation (7) can be evaluated at any (3-dimensional) position in the slice by computing the corresponding using the usual flat-space formulas, then interpolating in the 2-dimensional surface grid.
2. Rewrite the apparent horizon condition (15) in terms of and its (3-dimensional) Cartesian derivatives,
Huq et al. [8990] give the function explicitly.
3. For each (latitude/longitude) grid point on the trial horizon surface, define a -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. Evaluate on the local Cartesian grid as described in Step 1 above.
5. Evaluate the Cartesian derivatives in Equation (27) by centered 2nd order Cartesian finite differencing of the 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).

#### 8.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. [8990] 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.

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 [132133] used the PETSc general-purpose elliptic-solver library [222324] 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. [146147] 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.

#### 8.5.4 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 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 [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 , is perturbed by some (small) amount at the th grid point (that is, ), and the expansion is recomputed. The th column of the Jacobian matrix (28) is then estimated as
Curtis and Reid [53], and Stoer and Bulirsch [150, Section 5.4.3] discuss the optimum choice of in this algorithm.

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

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. [146147], and Thornburg [156]. Here the internal structure of the finite differenced 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,

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, , 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.

#### 8.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 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 conditions (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. 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. Schnetter [132133] 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 decomposition solver (Davis’ UMFPACK library [5758565554]), and found each to be more efficient in some situations; overall, he found the UMFPACK solver to be the best choice.

#### 8.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.)

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. 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.)

#### 8.5.7 Summary of elliptic-PDE algorithms/codes

Elliptic-PDE apparent horizon finders have been developed by many researchers, including Eardley [67], Cook [505251], and Thornburg [153] in axisymmetry, and Shibata et al. [146147], Huq et al. [8990], Schnetter [132133], 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). 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 implemented. 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 [132133] 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 groups.