Go to previous page Go up Go to next page

5.3 Integrating null surfaces backwards in time

Anninos et al. [7Jump To The Next Citation Point], Libson et al. [103Jump To The Next Citation Point], and Walker [162Jump To The Next Citation Point] introduced the important concept of explicitly (numerically) finding the event horizon as a null surface in spacetime. They observed that if we parameterize the event horizon with any level-set function F satisfying the basic level-set definition (3View Equation), then the condition for the surface F = 0 to be null is just
ab g ∂aF ∂bF = 0. (9 )
Applying a 3 + 1 decomposition to this then gives a quadratic equation which can be solved to find the time evolution of the level-set function,
ti ∘ ------------------------ −-g-∂iF-+----(gti∂iF)2-−-gttgij ∂if-∂jf ∂tF = gtt . (10 )

Alternatively, assuming the event horizon in each slice to be a Strahlkörper in the manner of Section 2.2, we can define a suitable level-set function F by Equation (7View Equation). Substituting this definition into Equation (10View Equation) then gives an explicit evolution equation for the horizon shape function,

tr ru ∘ --tr---tu----2----tt--rr-----ru--------uv--------- ∂th = −-g--+-g--∂uh-+----(g--−-g--∂uh-)-−-g--(g--−-2g--∂uh-+--g--∂uh-∂vh). (11 ) gtt

Surfaces near the event horizon share the same “attraction” property discussed in Section 5.2 for geodesics near the event horizon. Thus by integrating either surface representation (10View Equation) or (11View Equation) backwards in time, we can refine an initial guess into a very accurate approximation to the event horizon.

In contrast to the null geodesic equation (8View Equation), neither Equation (10View Equation) nor Equation (11View Equation) contain any derivatives of the 4-metric (or equivalently the 3 + 1 geometry variables). This makes it much easier to integrate these latter equations accurately11. This formulation of the event-horizon finding problem also completely eliminates the tangential-drifting problem discussed in Section 5.2, since the level-set function only parameterizes motion normal to the surface.

5.3.1 Error bounds: Integrating a pair of surfaces

For a practical algorithm, it is useful to integrate a pair of trial null surfaces backwards: an “inner-bound” one which starts (and thus always remains) inside the event horizon and an “outer-bound” one which starts (and thus always remains) outside the event horizon. If the final slice contains an apparent horizon then any 2-surface inside this can serve as our inner-bound surface. However, choosing an outer-bound surface is more difficult.

It is this desire for a reliable outer bound on the event horizon position that motivates our requirement (Section 4) for the final slice (or at least its strong-field region) to be approximately stationary: In the absence of time-dependent equations of state or external perturbations entering the system, this requirement ensures that, for example, any surface substantially outside the apparent horizon can serve as an outer-bound surface.

Assuming we have an inner- and an outer-bound surface on the final slice, the spacing between these two surfaces after some period of backwards integration then gives an error bound for the computed event horizon position. Equivalently, a necessary (and, if there are no other numerical problems, sufficient) condition for the event-horizon finding algorithm to be accurate is that the backwards integration must have proceeded far enough for the spacing between the two trial surfaces to be “small”. For a reasonable definition of “small”, this typically takes at least 15 m ADM of backwards integration, with 20 m ADM or more providing much higher accuracy.

In some cases it is difficult to obtain a long enough span of numerical data for this backwards integration. For example, in some simulations of binary black hole collisions, the evolution becomes unstable and crashes soon after a common apparent horizon forms. This means that we cannot compute an accurate event horizon for the most interesting region of the spacetime, that which is close to the black-hole merger. There is no good solution to this problem except for the obvious one of developing a stable (or less-unstable) simulation that can be continued for a longer time.

5.3.2 Explicit Strahlkörper surface representation

The initial implementations of the “integrate null surface backwards” algorithm by Anninos et al. [7Jump To The Next Citation Point], Libson et al. [103Jump To The Next Citation Point], and Walker [162Jump To The Next Citation Point] were based on the explicit Strahlkörper surface integration formula (11View Equation), further restricted to axisymmetry12.

For a single black hole the coordinate choice is straightforward. For the two-black-hole case, the authors used topologically cylindrical coordinates (ρ,z, φ), where the two black holes collide along the axisymmetry (z) axis. Based on the symmetry of the problem, they then assumed that the event horizon shape could be written in the form

ρ = h(z) (12 )
in each t = constant slice.

This spacetime’s event horizon has the now-classic “pair of pants” shape, with a non-differentiable cusp along the “inseam” (the z axis ρ = 0) where new generators join the surface. The authors tried two ways of treating this cusp numerically:

View Image

Figure 4: This figure shows a view of the numerically-computed event horizon in a single slice, together with the locus of the event horizon’s generators that have not yet joined the event horizon in this slice, for a head-on binary black hole collision. Notice how the event horizon is non-differentiable at the cusp where the new generators join it. Figure reprinted with permission from [103]. © 1996 by the American Physical Society.
View Image

Figure 5: This figure shows a perspective view of the numerically-computed event horizon, together with some of its generators, for the head-on binary black hole collision discussed in detail by Matzner et al. [108Jump To The Next Citation Point]. Figure courtesy of Edward Seidel.

Caveny et al. [44Jump To The Next Citation Point46Jump To The Next Citation Point] implemented the “integrate null surfaces backwards” algorithm for fully generic numerically-computed spacetimes using the explicit Strahlkörper surface integration formula (11View Equation). To handle moving black holes, they recentered each black hole’s Strahlkörper parameterization (4View Equation) on the black hole’s coordinate centroid at each time step.

For single-black-hole test cases (Kerr spacetime in various coordinates), they report typical accuracies of a few percent in determining the event horizon position and area. For binary-black-hole test cases (Kastor–Traschen extremal-charge black hole coalescence with a cosmological constant), they detect black hole coalescence (which appears as a bifurcation in the backwards time integration) by the “necking off” of the surface. Figure 6View Image shows an example of their results.

View Image

Figure 6: This figure shows the cross-section of the numerically-computed event horizon in each of five different slices, for the head-on collision of two extremal Kastor–Traschen black holes. Figure reprinted with permission from [46]. © 2003 by the American Physical Society.

5.3.3 Level-set parameterization

Caveny et al. [44Jump To The Next Citation Point45Jump To The Next Citation Point] and Diener [60Jump To The Next Citation Point] (independently) implemented the “integrate null surfaces backwards” algorithm for fully generic numerically-computed spacetimes, using the level-set function integration formula (10View Equation). Here the level-set function F is initialized on the final slice of the evolution and evolved backwards in time using Equation (10View Equation) on (conceptually) the entire numerical grid. (In practice, only a smaller box containing the event horizon need be evolved.)

This surface parameterization has the advantage that the event-horizon topology and (non-)smoothness are completely unconstrained, allowing the numerical study of configurations such as toroidal event horizons (discussed in Section 4). It is also convenient that the level-set function F is defined on the same numerical grid as the spacetime geometry, so that no interpolation is needed for the evolution.

The major problem with this algorithm is that during the backwards evolution, spatial gradients in F tend to steepen into a jump discontinuity at the event horizon14, eventually causing numerical difficulty.

Caveny et al. [4445] deal with this problem by adding an artificial viscosity (i.e. diffusion) term to the level-set function evolution equation, smoothing out the jump discontinuity in F. That is, instead of Equation (10View Equation), they actually evolve F via

∂tF = rhsEq.(10) + ɛ2∇2F, (13 )
where rhsEq.(10) is the right hand side of Equation (10View Equation) and ∇2 is a generic 2nd order linear (elliptic) spatial differential operator, and ɛ > 0 is a (small) dissipation constant. This scheme works, but the numerical viscosity does seem to lead to significant errors (several percent) in their computed event-horizon positions and areas15, and even failure to converge to the correct solution for some test cases (e.g. rapidly-spinning Kerr black holes).

Alternatively, Diener [60Jump To The Next Citation Point] developed a technique of periodically reinitializing the level-set function to approximately the signed distance from the event horizon. To do this, he periodically evolves

∂λF = − √--F-----(|∇F | − 1) (14 ) F 2 + 1
in an unphysical “pseudo-time” λ until an approximate steady state has been achieved. He reports that this works well in most circumstances but can significantly distort the computed event horizon if the F = 0 isosurface (the current approximation to the event horizon) is only a few grid points thick in any direction, as typically occurs just around the time of a topology change in the isosurface. He avoids this problem by estimating the minimum thickness of this isosurface and, if it is below a threshold, deferring the reinitialization.

In various tests on analytical data, Diener [60Jump To The Next Citation Point] found this event-horizon finder, EHFinder, to be robust and highly accurate, typically locating the event horizon to much less than 1% of the 3-dimensional grid spacing. As an example of results obtained with EHFinder, Figure 7View Image shows two views of the numerically-computed event horizon for a spiraling binary black hole collision. As another example, Figure 8View Image shows the numerically-computed event and apparent horizons in the collapse of a rapidly rotating neutron star to a Kerr black hole. (The apparent horizons were computed using the AHFinderDirect code described in Section 8.5.7.)

EHFinder is implemented as a freely available module (“thorn”) in the External LinkCactus computational toolkit (see Table 2). It originally worked only with the PUGH unigrid driver, but work is ongoing [61] to enhance it to work with the External LinkCarpet mesh-refinement driver [134131].

View Image

Figure 7: This figure shows two views of the numerically-computed event horizon’s cross-section in the orbital plane for a spiraling binary black hole collision. The two orbital-plane dimensions are shown horizontally; time runs upwards. The initial data was constructed to have an approximate helical Killing vector, corresponding to black holes in approximately circular orbits (the D = 18 case of Grandclément et al. [78]), with a proper separation of the apparent horizons of 6.9 m. (The growth of the individual event horizons by roughly a factor of 3 in the first part of the evolution is an artifact of the coordinate choice – the black holes are actually in a quasi-equilibrium state.) Figure courtesy of Peter Diener.
View Image

Figure 8: This figure shows the polar and equatorial radii of the event horizon (solid lines) and apparent horizon (dashed lines) in a numerical simulation of the collapse of a rapidly rotating neutron star to form a Kerr black hole. The dotted line shows the equatorial radius of the stellar surface. These results are from the D4 simulation of Baiotti et al. [21Jump To The Next Citation Point]. Notice how the event horizon grows from zero size while the apparent horizon first appears at a finite size and grows in a spacelike manner. Notice also that both surfaces are flattened due to the rotation. Figure reprinted with permission from [21]. © 2005 by the American Physical Society.

  Go to previous page Go up Go to next page