Finding apparent horizons in numerical relativity

We review various algorithms for finding apparent horizons in 3+1 numerical relativity. We then focus on one particular algorithm, in which we pose the apparent horizon equation H=del(i)n(i)+K(ij)n(i)n(j)-K=0 as a nonlinear elliptic (boundary-value) PDE on angular-coordinate space for the horizon shape function r=h(theta,phi), finite difference this PDE, and use Newton's method or a variant to solve the finite difference equations. We describe a method for computing the Jacobian matrix of the finite differenced H(h) function H(h) 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(iz) equations, then finite differencing the linearized continuum equations. (This is essentially just the ''Jacobian part'' of the Newton-Kantorovich method for solving nonlinear PDEs). We tabulate the resulting Jacobian coefficients for a number of different H(h) and Jacobian computation schemes. We find this symbolic differentiation method of computing the H(h) Jacobian to be 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 shows robust convergence. However, we find that it has a small (poor) radius of convergence if the initial guess for the horizon position contains significant high-spatial-frequency error components, i.e., angular; Fourier components varying as (say) cosm theta with m greater than or similar to 8. (Such components occur naturally if spacetime contains significant amounts of high-frequency gravitational radiation.) We show that this poor convergence behavior is 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. We find that a simple ''line search'' modification of Newton's method roughly doubles the horizon finder's radius of convergence, but both the unmodified and modified methods' radia of convergence still fall rapidly with increasing spatial frequency, approximately as 1/m(3/2). Further research is needed to explore more robust numerical algorithms for solving the H(h)=0 equations. Provided it converges, the Newton's-method algorithm for horizon finding is potentially very accurate, in practice limited only by the accuracy of the H(h) finite differencing scheme. Using fourth order finite differencing, we demonstrate that the error in the numerically computed horizon position shows the expected O((Delta theta)(4)) scaling with grid resolution a Delta theta, and is typically similar to 10(-5)(10(-6)) for a grid resolution of Delta theta=pi/2/50(pi/2/100). Finally, we briefly discuss the global problem of finding or recognizing the outermost apparent horizon in a slice. We argue that this is an important problem, and that no reliable algorithms currently exist for it except in spherical symmetry.