A Fast Apparent-Horizon Finder for 3-Dimensional Cartesian Grids in Numerical Relativity

TL;DR

AHFinderDirect solves the apparent-horizon equation by treating the horizon as a Strahlkörper with shape $r = h(\rho,\sigma)$ on $S^2$ and solving a nonlinear elliptic PDE for $h$ on the sphere. It discretizes $S^2$ with a multi-patch angular grid, interpolates 3D geometry to the horizon using cubic Hermite interpolation, and constructs a symbolic Jacobian to drive Newton iterations efficiently, with parallelization for multiple horizons. The method achieves typical horizons at $\sim 10^{-5} m$ positional accuracy in seconds to a few seconds per horizon and outperforms prior fast-flow approaches by an order of magnitude in speed and by 1–2 orders in accuracy across diverse test spacetimes, including boosted Kerr, Misner/Brill-Lindquist initial data, and binary black-hole mergers. This enables horizon finding at every time step of a 3D evolution, facilitating improved coordinate control and diagnostics in numerical relativity simulations. The code AHFinderDirect is implemented as a Cactus thorn and released under the GNU GPL, promoting easy porting and broad adoption.

Abstract

In 3+1 numerical simulations of dynamic black hole spacetimes, it's useful to be able to find the apparent horizon(s) (AH) in each slice of a time evolution. A number of AH finders are available, but they often take many minutes to run, so they're too slow to be practically usable at each time step. Here I present a new AH finder,_AHFinderDirect_, which is very fast and accurate: at typical resolutions it takes only a few seconds to find an AH to $\sim 10^{-5} m$ accuracy on a GHz-class processor. I assume that an AH to be searched for is a Strahlkörper (star-shaped region) with respect to some local origin, and so parameterize the AH shape by $r = h(angle)$ for some single-valued function $h: S^2 \to \Re^+$. The AH equation then becomes a nonlinear elliptic PDE in $h$ on $S^2$, whose coefficients are algebraic functions of $g_{ij}$, $K_{ij}$, and the Cartesian-coordinate spatial derivatives of $g_{ij}$. I discretize $S^2$ using 6 angular patches (one each in the neighborhood of the $\pm x$, $\pm y$, and $\pm z$ axes) to avoid coordinate singularities, and finite difference the AH equation in the angular coordinates using 4th order finite differencing. I solve the resulting system of nonlinear algebraic equations (for $h$ at the angular grid points) by Newton's method, using a "symbolic differentiation" technique to compute the Jacobian matrix._AHFinderDirect_ is implemented as a thorn in the_Cactus_ computational toolkit, and is freely available by anonymous CVS checkout.

Figures (6)

• Figure 1: This figure shows the errors for cubic Lagrange and Hermite interpolation of the function $f(x) = \exp[\sin(2 \pi x)]$ with grid spacing $\Delta x = 0.1$. Notice that the Lagrange error (and hence the Langrange interpolant itself) is non-differentiable at the grid points, whereas the Hermite error (and interpolant) is differentiable everywhere.
• Figure 2: This figure shows a multiple-grid-patch system covering the $(+,+,+)$ octant of $S^2$ with 3 patches, at an angular resolution of $5^\circ$. The $+z$, $+x$, and $+y$ patches are shown in red, green, and blue respectively. Each patch's nominal grid is shown in thick lines; the ghost zones are shown in thin lines.
• Figure 3: This figure shows the areas of the various AHs in the Misner $\mu = 2.0$ collision described in table \ref{['tab-misner12-066-parameters']}. The black points are the areas found by the fast flow AH finder; the other curves are all from AHFinderDirect. The gradual rise in the area of the outer common AH after $t \approx 9$, and in the area of the individual AH after $t \approx 15$, is due to outer boundary reflections making the overall evolution inaccurate.
• Figure 4: This figure shows the three AHs in the Misner $\mu = 2.0$ collision described in table \ref{['tab-misner12-066-parameters']}. Part (a) shows the horizons at $t = 5.00$ ($3.93m$); part (b) shows them at $t = 8.00$ ($6.28m$). In both parts the color coding matches that of figure \ref{['fig-misner12-066/AH-areas']}.
• Figure 5: This figure shows the combinations of symmetry operations and interpatch interpolations used to fill in grid function values in ghost-zone corners. Part (a) (where there are both $x \leftrightarrow -x$ and $y \leftrightarrow -y$ reflection symmetries) shows a corner between two symmetry ghost zones. Part (b) (where there is a $y \leftrightarrow -y$ reflection symmetry) shows a corner between a symmetry and an interpatch ghost zone; sample output points for each phase of the 3-phase algorithm described in the text are labelled as 1, 2, and 3. Part (c) (where there are no symmetries) shows a corner between two interpatch ghost zones (this only happens when 3 patches meet at a corner). In each part, arrows show the symmetry operations, and for parts (b) and (c) the boxed and circled points show the inputs and outputs for the interpatch interpolations.