Finding Apparent Horizons in Numerical Relativity

TL;DR

This work presents a thorough framework for finding apparent horizons in 3+1 numerical relativity using Newton's method applied to a discretized elliptic horizon equation for the horizon shape $h(\theta,\phi)$. A key contribution is the symbolic-differentiation approach to the horizon function Jacobian, notably the 3d.sd.1s method, which outperforms traditional numerical perturbation by a substantial margin and approaches the cost of evaluating the horizon function itself. The paper also analyzes convergence properties, including global convergence and the adverse impact of high-spatial-frequency initial errors, supported by 2-grid and 3-grid convergence tests and a Monte Carlo horizon-perturbation survey. It demonstrates high-precision horizon positions (often $\sim 10^{-5}$ to $10^{-6}$) with 4th-order finite differencing and discusses the open problem of reliably locating outermost horizons in nonspherical spacetimes. Overall, the results indicate that symbolic Jacobian computation combined with Newton’s method yields a fast, accurate horizon finder, while highlighting the need for robust global strategies in challenging geometries and for outer-horizon recognition.

Abstract

This paper presents a detailed discussion of the ``Newton's method'' algorithm for finding apparent horizons in 3+1 numerical relativity. We describe a method for computing the Jacobian matrix of the finite differenced $H(h)$ function by symbolically differentiating the finite difference equations, giving the Jacobian elements directly in terms of the finite difference molecule coefficients used in computing $H(h)$. Assuming the finite differencing scheme commutes with linearization, we show how the Jacobian elements may be computed by first linearizing the continuum $H(h)$ equations, then finite differencing the linearized (continuum) equations. We find this symbolic differentiation method of computing the $H(h)$ Jacobian to be {\em much} more efficient than the usual numerical perturbation method, and also much easier to implement than is commonly thought. When solving the discrete $H(h) = 0$ equations, we find that Newton's method generally converges very rapidly. However, if the initial guess for the horizon position contains significant high-spatial-frequency error components, Newton's method has a small (poor) radius of convergence. This is {\em not} an artifact of insufficient resolution in the finite difference grid; rather, it appears to be caused by a strong nonlinearity in the continuum $H(h)$ function for high-spatial-frequency error components in $h$. Robust variants of Newton's method can boost the radius of convergence by O(1) factors, but the underlying nonlinearity remains, and appears to worsen rapidly with increasing initial-guess-error spatial frequency. Using 4th~order finite differencing, we find typical accuracies for computed horizon positions in the $10^{-5}$ range for $Δθ= \frac{π/2}{50}$.

Figures (7)

• Figure 1: This figure illustrates the various 2-stage and 1-stage computation methods for the horizon function ${\sf H}({\sf h})$. The solid arrows denote finite differencing operations, the dotted arrow denotes an algebraic computation, and the dashed arrow denotes a radial interpolation to the horizon position $r = {\sf h}(\theta,\phi)$. Each path from ${\sf h}$ to ${\sf H}$ represents a separate computation method. Notice that there are three distinct 2-stage methods (using the upper arrows from ${}^{(2)} {\sf h}$ to ${}^{(3)} {\sf H}$ in the figure) and one 1-stage method (using the lower arrow from ${}^{(2)} {\sf h}$ to ${}^{(3)} {\sf H}$).
• Figure 2: This commutative diagram illustrates the two different ways a Jacobian matrix can be computed. Given a nonlinear continuum function $P(Q)$, the Jacobian matrix ${\bf J} \Bigl[ {\sf P}({\sf Q}) \Bigr]$ is logically defined in terms of the lower-left path in the diagram, i.e. it's defined as the Jacobian of a nonlinear discrete (finite difference) approximation ${\sf P}({\sf Q})$ to $P(Q)$. However, if the operations of discretization (finite differencing) and linearization commute, we can instead compute the Jacobian by the upper-right path in the diagram, i.e. by first linearizing the continuum $P(Q)$ function, then discretizing (finite differencing) this linearization $\delta P(\delta Q)$.
• Figure 3: This figure shows $H(r)$ for spherical trial horizon surfaces with coordinate radius $r$ in an Eddington-Finkelstein slice of a unit-mass Schwarzschild spacetime. Notice that for $r > r^{\max} \approx 4.372$, $H > 0$ and $dH / dr < 0$, so Newton's method diverges in this region.
• Figure 4: This figure illustrates how the convergence behavior of the basic and modified Newton iterations depends on the spatial-frequency spectrum of the initial guess's error ${\sf h}^{(0)} - {\sf h}^\ast$. In each part of the figure, the true continuum horizon $h^\ast$ is plotted as a solid line, while the horizon finder's first few iterates (trial horizon surfaces) ${\sf h}^{(k)}$ are plotted with dots at the grid points. Part (a) of the figure shows the behavior of Newton's method for an initial-guess-error containing only low spatial frequencies, part (b) shows the behavior of Newton's method for an initial-guess-error containing significant high spatial frequencies, and part (c) shows the behavior of the modified Newton iteration for the same initial guess as part (b). In parts (a) and (c), where the iteration is converging, the final iterates shown are indistinguishable from the true continuum horizon at the scale of the figure. In part (b), where the iteration is diverging, the computed values for the next iterate ${\sf h}^{(3)}$ (not shown) are almost all far outside the scale of the figure; many of them are in fact negative!
• Figure 5: This figure shows the results of a 3-grid covergence test for the 2nd-iteration Newton iterate (trial horizon surface) ${\sf h}^{(2)}$ plotted in figure \ref{['fig-Kerr-hp4-hp10']}(b). The line has slope ${ \frac{1}{16}}$, appropriate for 4th order convergence. (Recall that this line isn't fitted to the data, but is rather an a priori prediction with no adjustable parameters.) (The absolute magnitude of the errors shown here is much larger than is typical for our horizon finder, due to a combination of the compounding of smaller errors in the earlier Newton iterate ${\sf h}^{(1)}$, and the very strong angular variation in both iterates ${\sf h}^{(1)}$ and ${\sf h}^{(2)}$.)