$$Date$$

Abstract

This thorn locates the Event Horizon (EH) in an analytic or numerical spacetime by evolving a null surface backwards in time. The null surface is described at each time step as the 0-level iso-surface of a 3D scalar function $f$. This level set description of the surface allows, trivially, changes of the topology of the surface so it is possible to follow the merger of two (or more) black holes into a ﬁnal black hole.

This thorn attempts to locate the Event Horizon (EH) in an analytic or numerical spacetime by evolving a null surface backwards in time. This method depends on the fact that, except in cases where the coordinates are adapted to outgoing null geodesics, an outgoing null surface started close to the EH, when evolved forward in time, is diverging exponentially from the EH. Reversing the time evolution then means that an outgoing null surface will converge exponentially to the EH. The level set function, $f$, is evolved according to

$$\begin{array}{rcll}{\partial}_{t}f& =& \frac{-{g}^{ti}{\partial}_{i}f+\sqrt{{\left({g}^{ti}{\partial}_{i}f\right)}^{2}-{g}^{tt}{g}^{ij}{\partial}_{i}f{\partial}_{j}f}}{{g}^{tt}}& \text{}\\ & =& {\beta}^{i}{\partial}_{i}f-\sqrt{{\alpha}^{2}{\gamma}^{ij}{\partial}_{i}f{\partial}_{j}f},& \text{(1)}\text{}\text{}\end{array}$$

where in the second equation the lapse, shift and 3-metric has been substituted for the 4-metric. For more details on the theory and implementation see [1].

This thorn uses a level set description of the null surface, i.e.the surface is the 0-level isosurface of a 3D scalar function, $f$, that is negative inside and positive outside the surface. With this choice of surface description one level-set function can describe multiple surfaces at the same time, simply by having several, disconnected regions with negative values. The biggest advantage, however, is that any change of topology of the surface is handled naturally and simply by the level-set function changing sign. Therefore this EHFinder can follow the EH during the merger of two (or more) black holes into one black hole.

To ﬁnd the EH in a numerical spacetime several points have to be taken into consideration. Since the null surface has to be evolved backwards in time, the EHFinder has to be seen as a pre-processing analysis thorn. Therefore it is necessary to evolve the initial data forward in time while outputting enough 3D data, that the full 4-metric can be recovered at each timestep. The 3D data is then read back in, in reverse order, while the level-set function is evolved backwards in time. The thorn can evolve more than one level set function at a time using diﬀerent initial guesses for the surfaces (NOTE: this is a quite recent feature and has not yet been tested extensively). More details about the actual use of the thorn in section 7

The evolution equation for $f$, equation (1), causes steepening of the gradient of $f$, which is undesireble numerically. For that reason, $f$ is periodically re-initialized to a distance function. That is, the values of $f$ are changed so that the the value of $f$ in a grid point is equal to the (signed) distance from the grid point to the surface $f=0$. This is done by evolving $f$ according to the following evolution equation (in the parameter $\lambda $)

$$\frac{df}{d\lambda}=-\frac{f}{\sqrt{{f}^{2}+1}}\left(|\nabla f|-1\right)$$ | (2) |

until a steady state is achieved. This method is called the pde-method since it is basically evolving a pde. Sometimes the $f=0$ surface can be moved slightly during the re-initialization procedure. This happens when the surface develops a narrow neck just before a topology change. For that reason, there is an option to detect when this is about to happen and undo the re-initialization.

There used to be another method doing the re-initialization, but it proved inferior to the pde-method and was removed. Other methods may be implemented in the future.

Currently three diﬀerent choices for the initial shape of the surface are implemented. The simplest choice is a sphere in which case $f$ is set according to

$$f=\sqrt{{\left(x-{x}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}+{\left(z-{z}_{0}\right)}^{2}}-{r}_{0},$$ | (3) |

where ${r}_{0}$ is the radius of the sphere and ${x}_{0}$, ${y}_{0}$ and ${z}_{0}$ are the coordinates of the center of the sphere. The second choice is a rotated and translated ellipsoid. The basic equation is here

$$f=\sqrt{\frac{{x}^{2}}{{a}^{2}}+\frac{{y}^{2}}{{b}^{2}}+\frac{{z}^{2}}{{c}^{2}}}-1$$ | (4) |

This ellipsoid is ﬁrst rotated an angle $\alpha $ around the $z$-axis, then rotated an angle $\beta $ around the $y$-axis, then rotated an angle $\gamma $ around the $x$-axis and ﬁnally the rotated ellipsoid is translated to have its “center” at the point $\left({x}_{0},{y}_{0},{z}_{0}\right)$. The ﬁnal possible shape of the initial surface is an ovaloid of Cassini. This was implented initially to test changing the topology in ﬂat space. it is most likely not useful for numerical data. In this case $f$ is set according to

$$f={\left({x}^{2}+{y}^{2}+{z}^{2}\right)}^{2}+{a}^{4}-2{a}^{2}\left({x}^{2}-{y}^{2}-{z}^{2}\right)-{b}^{4}.$$ | (5) |

There are no translation or rotations implemented for the ovaloid of Cassini.

Diﬀerent initial shapes can be used for the diﬀerent level set functions used in the same run.

Even though the level set function, $f$, in principle can be deﬁned everywhere it is often advantageous to only evolve it in a certain region around the surface $f=0$. Since $f$ is re-initialized regularly to become a distance function, $f$ itself can be used as a measure of the distance from a certain grid point to the surface $f=0$. The parameter ehfinder::shell_width speciﬁes the size of the active region around $f=0$. However the interior and exterior excisions are handled diﬀerently. The interior excision is most simple, since here all grid points with $f<-\text{ehfinder::shell\_width}$ are marked as inactive. This should work in all cases when the excised region is everywhere concave, since then all points on the boundary of the excised region will have enough neighbouring active points to be able to calculate second order accurate one sided derivatives. If the interior excised region happens to have a convex region, this might fail. To avoid a similiar problem at the outer excised boundary, this boundary is shaped as a rectangular box. The box is chosen so that all points with $f<\text{ehfinder::shell\_width}$ are in the active region. This is illustrated in Figure 1, for the case

ehfinder::shell_width = 4, where the excised regions are hashed.

Changes to the excision regions are only done after re-initialization, since it is only at this time that $f$ is a distance function. The excision regions can move across the grid, tracking the surface $f=0$.

Also if the numerical run was done with excision using SpaceMask, the EHFinder excision region is guaranteed to completely cover the numerical excision region. The EHFinder currently only supports the old style excision mask but the support of the new style excision mask should be trivial and fast to implement.

All ﬁnite diﬀerences used in the evolution of the null surface are second order one sided diﬀerences. For that reason a ghost_size larger or equal to 2 should always be used. It is possible to choose between three diﬀerent upwinding schemes. This is done by setting the parameter ehfinder::upwind_type to either intrinsic, shift or characteristic.

The intrinsic scheme, looks at the values of $f$ itself, to determine the direction of the stencil. This is basically in order to be able to handle situations like the one illustrated in 1D in Figure 2.

If the stencil for calculating derivatives in the point labeled 1 is taken to consist of the points 1, 2 and 3’, the non diﬀerentiablility of $f$ may cause problems. The algorithm detects this and instead uses the points 1, 2, 3 as the stencil. This ensures that a non diﬀerentiable feature can be maintained in the evolution.

However this scheme only works in a few simple cases. For more general cases it is important to use characteristic information in order to ensure that the numerical stencil contains the domain of dependence. Therefore the characteristic scheme was introduced. Fortunately, by using the information contained in the level set functions the characteristic of the level set function can be calculated using

$$\frac{d{x}^{i}}{dt}=-{\beta}^{i}+\frac{{\alpha}^{2}{\gamma}^{ij}{\partial}_{j}f}{\sqrt{{\alpha}^{2}{\gamma}^{kl}{\partial}_{k}f{\partial}_{l}f}}.$$ | (6) |

The method estimates the characteristic direction using centered ﬁnite diﬀerences in equation (6) and then recalculates the one sided ﬁnite diﬀerences in the appropriate direction.

It might happen that the upwinding direction based on the characteristic direction results in the stencil to consist of the points 1, 2, 3’ in Figure 2. However, if the re-initialization is done often enough, this turns out not to cause any problems.

The shift scheme was implemented as an improvement to the intrinsic scheme but it turned out that the characteristic scheme was superior. Therefore shift upwinding should not be used. It might be completely removed in future revisions.

For the re-initialization the default is to use the intrinsic second order scheme (the re-initialization doesn’t depend on where the surface is moving), so characteristic upwinding is not applicable. It is possible to use a ﬁrst order intrinsic scheme, but this is, in my experience, not accurate enough. A centered diﬀerencing scheme is also available, but is only there for testing purposes and should never be used. These alternative schemes will probably be removed in the future.

Here the most important parameters are described.

- ehfinder::mode

The mode can either be set to normal (normal event horizon ﬁnder mode), analysis (compare a previously calculated level set function to a small number of analytic spacetimes) and generator (only evolve the generators while keeping the level set function ﬁxed). The default is normal and should normally not be changed. The other modes are only useful for debugging and testing purposes. - ehfinder::eh_number_level_sets

An integer parameter specifying how many individual level set functions to evolve at a time. Currently it has to be between 1 and 10. - ehfinder::eh_metric_type

The metric type can either be set to numerical or analytic. If it is set to numerical the EHFinder will attempt to read in the metric from ﬁles in the directory speciﬁed by the io::recover_dir parameter. At present all the timesteps has to be saved in the same ﬁle. Note that if the numerical data was produced with admbase::metric_type = "static conformal" this parameter has to be speciﬁed again. In this case the EHFinder will attempt to also read in the conformal factor from a ﬁle. If metric type is set to analytic another thorn needs to set up the metric. It is possible to only set the metric on the initial slice, but it is also possible to have a thorn (like thorn Exact) set the metric at CCTK_PRESTEP if the analytic metric is time dependent. The default is numerical. - ehfinder::eh_lapse_type

The same for the lapse. - ehfinder::eh_shift_type

The same for the shift. - initial_f[i]

A vector parameter specifying the initial shape of the null surface for the individual level set functions. The initial shape can currently be chosen from sphere, ellipsoid and cassini as described in section 3. The default is sphere. - initial_rad[i]

A vector parameter specifying the radius of the initial sphere (${r}_{0}$ in equation 3). The deafault is 1. - translate_x[i]

A vector parameter specifying how much to translate the initial surface in the $x$-direction (${x}_{0}$ in equation 3). Also used for the initial ellipsoid. The default is 0. - translate_y[i]

A vector parameter specifying how much to translate the initial surface in the $y$-direction (${y}_{0}$ in equation 3). Also used for the initial ellipsoid. The default is 0. - translate_z[i]

A vector parameter specifying how much to translate the initial surface in the $z$-direction (${z}_{0}$ in equation 3). Also used for the initial ellipsoid. The default is 0. - initial_a[i]

A vector parameter specifying $a$ in equation 4. The default is 1. - initial_b[i]

A vector parameter specifying $b$ in equation 4. The default is 1. - initial_c[i]

A vector parameter specifying $c$ in equation 4. The default is 1. - rotation_alpha[i]

A vector parameter specifying the rotation angle $\alpha $ for the ellipsoid around the $z$-axis. The default is 0. - rotation_beta[i]

A vector parameter specifying the rotation angle $\beta $ for the ellipsoid around the $y$-axis. The default is 0. - rotation_gamma[i]

A vector parameter specifying the rotation angle $\gamma $ for the ellipsoid around the $x$-axis. The default is 0. - shell_width

The width of the active evolution region. Grid points more than shell_width gridspacings away from the $f=0$ surface are marked as inactive and are not evolved as described in section 4. The default is $7$ gridspacings. - use_inner_excision

A boolean parameter specifying whether the interior excision should be used or not. - use_outer_excision

A boolean parameter specifying whether the exterior excision should be used or not. - upwind_type

The type of upwinding to be used (either intrinsic, shift or characteristic). See the detailed description of the upwinding types in section 5. The default is characteristic. - surface_direction

The code can track both outgoing and ingoing null surfaces. Choose the direction by using outward or inward. The default is outward. Note that the code only works as an event horizon ﬁnder when evolving outward going null surfaces backwards in time. - re_init_undo

Should the re-initialization be undone just before pinch-oﬀ or not as described in section 2. The default is "no". - re_init_int_method

Choose the integration method in the pde-re-initialization method. Choose either a simple Euler (euler) integration scheme or a second order Runge-Kutta (rk2) scheme. Since a pde is evolved to steady state, it seems that the Euler scheme works just ﬁne and is faster than the Runge-Kutta scheme. The default is euler. - re_init_max_iter

The maximum number of iterations in the pde-re-initialization scheme, before giving up. Unless you are working at high resolution the default should be enough. The default is 800. - pde_differences

Choose the type of ﬁnite diﬀerencing used in the pde re-initialization. Don’t ever use anything else than second order upwinding (upwind2). The other choices (centered and upwind) are there only for testing purposes and might be removed. - re_initialize_every

How often to re-initialize using the pde re-initialization. This depends on the problem. For some problems it is necessary to do it more often than for other problems. You’ll have to experiment to ﬁgure out what works best. The default is 10. - last_iteration_number

The last iteration number of the numerical simulation that produced the metric data. Active when eh_metric_type is equal to numerical. - saved_iteration_every

How often was the numerical data saved? This and the above parameter is used in the code to ﬁgure out which data set iteration number to read in from ﬁle

(last_iteration_number-saved_iteration_every*cctk_iteration). - ntheta

How many points in the $\mathit{\theta}$-direction should be used for integrations over the surface. - nphi

How many points in the $\varphi $-direction should be used for integrations over the surface. - maximum_surface_number

The maximal number of surfaces expected at any given time during the whole evolution. - surface_interpolator

Which interpolator to use for the location of points on the surface. It should be the name of a valid interpolator provided by one of the available interpolators. The interpolator should be able to return both the interpolated value and derivatives. At present this is provided by AEILocalInterp. The default is Hermite polynomial interpolation. - surface_interpolation_order

The interpolation order used for ﬁnding points on the surface. Higher orders require larger number of ghost zones for parallel runs. The default value is 2 (consistent with a ghostsize of 2). - area_interpolator

Which interpolator should be used to interpolate metric information onto the surface once the surface points have been found. The default is Lagrange polynomial interpolation. - area_interpolation_order

The interpolation order used for calculating the area of the surface. The default value is 3 (consistent with a ghostsize of 2).

EHFinder also extends the following parameter from admbase in order to be able to read in initial data.

- admbase::initial_data is extended with "read from file".
- admbase::initial_lapse is extended with "read from file".
- admbase::initial_shift is extended with "read from file".

In this section I will try to describe in little more detail how EHFinder can be used to ﬁnd the EH in a numerical spacetime.

The ﬁrst thing to make sure, is that enough data is output during the numerical run, to be able to reconstruct
the full 4-metric. The required output therefore consists of the ADMBase::metric, ADMBase::lapse and
ADMBase::shift. However if the evolution was done with zero shift and/or lapse equal to one, it is not necessary
to output these grid functions, as long as storage is turned on and they are set to the correct value initially when
EHFinder is run. If the evolution was done with a conformal factor then StaticConformal::psi has to be
output as well, since it is necessary in order to reconstruct the 4-metric. However it is only necessary to output
this once since it is constant during the evolution. It is not necessary to output the derivatives of
the conformal factor. If excision was used in the numerical run then it is also necessary to output
SpaceMask::emask^{1} ,
since EHFinder has to make sure that its internal mask covers the space mask.

At present it is necessary to output all timesteps into the same ﬁle (use IO::out_timesteps_per_file = -1, which is the default). In principle both FLEXIO and HDF5 output should be supported, but only HDF5 output has been tested. Since EHFinder can be run on a lot less processors compared to the spacetime evolution, it is often advantageous to either do unchuncked output or to recombine the output ﬁles, since it is then possible to read the data onto a smaller number of processors (use IO::out_unchunked = "yes" to write unchunked data). If the numerical run is larger than the EH containing region (hopefully that is the case; otherwise the boundaries are deﬁnitely to close in), it is possible to use hyperslabbing to just output the EH containing region (see for example CactusPUGHIO/IOHDF5 for details on this). If hyperslabbing is used it is deﬁnitely necessary to do the output unchunked or recombine afterwards. An example parameter ﬁle can be seen in the par/Misner_2.2_80_3D.par.

In principle EHFinder should also work for downsampled (in both space and time) data, but no experiments have been done so far to estimate the loss of accuracy (I have always used the full resolution and done output at every timestep).

If hyperslabbing and/or downsampling is used, it is the users responsibility (by specifying the right parameters in the parameter ﬁle) to ensure that EHFinder is run with the correct grid spacing and time step.

EHFinder is still under development and testing and can as yet not be used as a black box. But still I can give some guidelines and advice on how to proceed.

The ﬁrst concern is to setup the initial guess for the surface. Ideally one would like to use at least two diﬀerent initial surfaces. One surface completely inside the EH and one surface completely outside the EH. In practice it is often a good idea to use more than two diﬀerent initial surfaces, since then the initial surfaces closest to the EH can be identiﬁed. Using the feature of evolving multiple level set functions at a time, this can be done while reading in the numerical data only once. The easiest way to get an initial surface inside the EH, is to set up the initial guess to be completely inside the apparent horizon (AH).