### 8.7 Flow algorithms

Flow algorithms define a “flow” on 2-surfaces, i.e. they define an evolution of 2-surfaces in some pseudo-time , such that the apparent horizon is the limit of a (any) suitable starting surface. Flow algorithms are different from other apparent horizon finding algorithms (except for zero-finding in spherical symmetry) in that their convergence does not depend on having a good initial guess. In other words, flow algorithms are global algorithms (Section 7.7).

To find the (an) apparent horizon, i.e. an outermost MOTS, the starting surface should be outside the largest possible MOTS in the slice. In practice, it generally suffices to start with a 2-sphere of areal radius substantially greater than .

The global convergence property requires that a flow algorithm always flow from a large starting surface into the apparent horizon. This means that the algorithm gains no particular benefit from already knowing the approximate position of the apparent horizon. In particular, flow algorithms are no faster when “tracking” the apparent horizon (repeatedly finding it at frequent intervals) in a numerical time evolution. (In contrast, in this situation a local apparent horizon finding algorithm can use the most recent previously-found apparent horizon as an initial guess, greatly speeding the algorithm’s convergence).

Flow algorithms were first proposed for apparent horizon finding by Tod [157]. He initially considered the case of a time-symmetric slice (one where ). In this case, a MOTS (and thus an apparent horizon) is a surface of minimal area and may be found by a “mean curvature flow”

where are the spatial coordinates of a horizon-surface point, is as before the outward-pointing unit 3-vector normal to the surface, and is the mean curvature of the surface as embedded in the slice. This is a gradient flow for the surface area, and Grayson [79] has proven that if the slice contains a minimum-area surface, this will in fact be the stable limit of this flow. Unfortunately, this proof is valid only for the time-symmetric case.

For non-time-symmetric slices, Tod [157] proposed generalizing the mean curvature flow to the “expansion flow”

There is no theoretical proof that this flow will converge to the (an) apparent horizon, but several lines of argument make this convergence plausible:
• The expansion flow is identical to the mean curvature flow (38) in the principal part.
• The expansion flow’s velocity is clearly zero on an apparent horizon.
• More generally, a simple argument due to Bartnik [25] shows that the expansion flow can never move a (smooth) test surface through an apparent horizon. Suppose, to the contrary, that the test surface is about to move through an apparent horizon , i.e. since both surfaces are by assumption smooth, that and touch at single (isolated) point . At that point, and obviously have the same and , and they also have the same (because is isolated). Hence the only term in (as defined by Equation (15)) which differs between and is . Clearly, if is outside and they touch at the single isolated point , then relative to , must be concave outwards at , so that . Thus the expansion flow (39) will move outwards, away from the apparent horizon. (If lies inside the same argument holds with signs reversed appropriately.)

Numerical experiments by Bernstein [28], Shoemaker et al. [148149], and Pasch [118] show that in practice the expansion flow (39) does in fact converge robustly to the apparent horizon.

In the following I discuss a number of important implementation details for, and refinements of, this basic algorithm.

#### 8.7.1 Implicit pseudo-time stepping

Assuming the Strahlkörper surface parameterization (4), the expansion flow (39) is a parabolic equation for the horizon shape function . This means that any fully explicit scheme to integrate it (in the pseudo-time ) must severely restrict its pseudo-time step for stability, and this restriction grows (quadratically) worse at higher spatial resolutions. This makes the horizon finding process very slow.

To avoid this restriction, practical implementations of flow algorithms use implicit pseudo-time integration schemes; these can have large pseudo-time steps and still be stable. Because we only care about the limit, a highly accurate pseudo-time integration is not important; only the accuracy of approximating the spatial derivatives matters. Bernstein [28] used a modified Du Fort–Frankel scheme [64] but found some problems with the surface shape gradually developing high-spatial-frequency noise. Pasch [118] reports that an “exponential” integrator (Hochbruck et al. [85]) works well, provided the flow’s Jacobian matrix is computed accurately. The most common choice is probably that of Shoemaker et al. [148149], who use the iterated Crank–Nicholson (“ICN”) scheme. They report that this works very well; in particular, they do not report any noise problems.

By refining his finite-element grid (Section 2.3) in a hierarchical manner, Metzger [109] is able to use standard conjugate-gradient elliptic solvers in a multigrid-like fashion, using each refinement level’s solution as an initial guess for the next higher refinement level’s iterative solution. This greatly speeds the flow integration: Metzger reports that the performance of the overall surface-finding algorithm is “of the same order of magnitude” as that of Thornburg’s AHFinderDirect [156] elliptic-PDE apparent horizon finder (described in Section 8.5.7).

In a more general context than numerical relativity, Osher and Sethian [116] have discussed a general class of numerical algorithms for integrating “fronts propagating with curvature-dependent speed”. These flow a level-set function (Section 2.1) which implicitly locates the actual “front”.

#### 8.7.2 Varying the flow speed

Another important performance optimization of the standard expansion flow (39) is to replace in the right-hand side by a suitable nonlinear function of , chosen so the surface shrinks faster when it is far from the apparent horizon. For example, Shoemaker et al. [148149] use the flow

for this purpose, where is the value of on the initial-guess surface, and (which is gradually decreased towards as the iteration proceeds) is a “goal” value for .

#### 8.7.3 Surface representation and the handling of bifurcations

Since a flow algorithm starts with (topologically) a single large 2-sphere, if there are multiple apparent horizons present the surface must change topology (bifurcate) at some point in the flow. Depending on how the surface is represented, this may be easy or difficult.

Pasch [118] and Shoemaker et al. [148149] use a level-set function approach (Section 2.1). This automatically handles any topology or topology change. However, it has the drawback of requiring the flow to be integrated throughout the entire volume of the slice (or at least in some neighborhood of each surface). This is likely to be much more expensive than only integrating the flow on the surface itself. Shoemaker et al. also generate an explicit Strahlkörper surface representation (Section 2.2), monitoring the surface shape to detect an imminent bifurcation and reparameterizing the shape into 2 separate surfaces if a bifurcation happens.

Metzger [109] uses a finite-element surface representation (Section 2.3), which can represent any topology. However, if the flow bifurcates, then to explicitly represent each apparent horizon the code must detect that the surface self-intersects, which may be expensive.

#### 8.7.4 Gundlach’s “fast flow”

Gundlach [80] introduced the important concept of a “fast flow”. He observed that the subtraction and inversion of the flat-space Laplacian in the Nakamura–Kojima–Oohara spectral integral-iteration algorithm (Section 8.4) is an example of “a standard way of solving nonlinear elliptic problems numerically, namely subtracting a simple linear elliptic operator from the nonlinear one, inverting it by pseudo-spectral algorithms and iterating”. Gundlach then interpreted the Nakamura–Kojima–Oohara algorithm as a type of flow algorithm where each pseudo-time step of the flow corresponds to a single functional-iteration step of the Nakamura–Kojima–Oohara algorithm.

In this framework, Gundlach defines a 2-parameter family of flows interpolating between the Nakamura–Kojima–Oohara algorithm and Tod’s [157] expansion flow (39),

where and are parameters, is a weight functional which depends on through at most 1st derivatives, is the flat-space Laplacian operator, and denotes inverting the operator . Gundlach then argues that intermediate “fast flow” members of this family should be useful compromises between the speed of the Nakamura–Kojima–Oohara algorithm and the robustness of Tod’s expansion flow. Based on numerical experiments, Gundlach suggests a particular choice for the weight functional and the parameters and . The resulting algorithm updates high-spatial-frequency components of essentially the same as the Nakamura–Kojima–Oohara algorithm but should reduce low-spatial-frequency error components faster. Gundlach’s algorithm also completely avoids the need for numerically solving Equation (23) for the coefficient in the Nakamura–Kojima–Oohara algorithm.

Alcubierre’s AHFinder [4] horizon finder includes an implementation of Gundlach’s fast flow algorithm. AHFinder is implemented as a freely available module (“thorn”) in the Cactus computational toolkit (see Table 2) and has been used by many research groups.

#### 8.7.5 Summary of flow algorithms/codes

Flow algorithms are the only truly global apparent horizon finding algorithms and, as such, can be much more robust than local algorithms. In particular, flow algorithms can guarantee convergence to the outermost MOTS in a slice. Unfortunately, these convergence guarantees hold only for time-symmetric slices.

In the forms which have strong convergence guarantees, flow algorithms tend to be very slow. (Metzger’s algorithm [109] is a notable exception: It is very fast.) There are modifications which can make flow algorithms much faster, but then their convergence is no longer guaranteed. In particular, practical experience has shown that in some binary black hole coalescence simulations (Alcubierre et al. [5], Diener et al. [62]), “fast flow” algorithms (Section 8.7.4) can miss common apparent horizons which are found by other (local) algorithms.

Alcubierre’s apparent horizon finder AHFinder [4] includes a “fast flow” algorithm based on the work of Gundlach [80]. It is implemented as a freely available module (“thorn”) in the Cactus computational toolkit (see Table 2) and has been used by a number of research groups.