Self-intersecting marginally outer trapped surfaces

TL;DR

This paper confirms the merger of marginally outer trapped surfaces (MOTSs) in binary black hole mergers, introduces improved numerical methods for locating MOTSs, and demonstrates the formation of self-intersecting MOTSs post-merger, aiding gravitational wave analysis.

Contribution

It provides numerical confirmation of MOTS mergers, details improved methods for locating MOTSs, and shows the formation of self-intersecting MOTSs after black hole mergers.

Findings

Confirmed MOTS merger scenario numerically

Demonstrated existence of self-intersecting MOTSs post-merger

Enhanced numerical methods for locating MOTSs

Abstract

We have shown previously that a merger of marginally outer trapped surfaces (MOTSs) occurs in a binary black hole merger and that there is a continuous sequence of MOTSs which connects the initial two black holes to the final one. In this paper, we confirm this scenario numerically and we detail further improvements in the numerical methods for locating MOTSs. With these improvements, we confirm the merger scenario and demonstrate the existence of self-intersecting MOTSs formed in the immediate aftermath of the merger. These results will allow us to track physical quantities across the non-linear merger process and to potentially infer properties of the merger from gravitational wave observations.

Equations

$${\Theta }_{(\ell )}=q^{ab}{\nabla }_a{\ell }_b,{\Theta }_{(n)}=q^{ab}{\nabla }_a{n}_b.$$

$${\delta }_{fr}{\Theta }_{(\ell )}=:Lf.$$

$$Lf=-\Delta f+{\gamma }^a{\partial }_{a}f+\beta f,$$

$${\ell }^a=\frac{1}{\sqrt{2}}\left({\tau }^a+r^a\right),n^a=\frac{1}{\sqrt{2}}\left({\tau }^a-r^a\right).$$

$$D_a{r}^a+K_{ab}{r}^a{r}^b-K=0.$$

$$F(r,\theta ,\varphi )=r-h(\theta ,\varphi ),$$

$$x=\frac{a\sin s}{\cosh t-\cos s},z=\frac{a\sinh t}{\cosh t-\cos s}+c.$$

$$h_{ij}$$

$$K_{ij}$$

$$\bar{x}=x\left(1-\frac{\beta }{\cosh ((z-z_0)/\gamma )}\right),\bar{z}=z.$$

$$\mathcal{H}=K_{ab}{K}^{ab}-K-R$$

Full text

I Introduction

Numerous binary black hole merger events have now been observed by gravitational wave detectors Abbott et al. (2016a, b, 2018); Green and Moffat (2018); Zackay et al. (2019); Nitz et al. (2019); Venumadhav et al. (2019). The general features of the gravitational wave signal from such events are now well known. The first is the inspiral regime where the signal is a chirp of increasing amplitude and frequency, and the system is effectively modeled as two point particles orbiting around each other and emitting gravitational waves as the orbit decays. As the two black holes approach each other and coalesce to form a final common black hole, the inspiral description is no longer valid, and non-perturbative aspects of general relativity become important; this is the merger regime. Eventually, as the final black hole reaches equilibrium, the gravitational wave signal can be well modeled as a superposition of damped sinusoids (and, in principle, much weaker power-law tails). Corresponding to this behavior of the gravitational wave signals, one visualizes the black holes themselves separately in the three different regimes. The inspiral regime consists of two disjoint black hole horizons slightly distorted by each other’s gravitational field. The merger is visualized as two horizons very close to each other and merging to form a single horizon which is initially very distorted. Finally, the ringdown is modeled as a perturbed Kerr horizon settling down to a final equilibrium Kerr black hole.

These features of the waveform must be correlated in some way with properties of the gravitational field in the strong field region. In particular, the three regimes must correspond in some way to properties of the black hole horizons. The details of the correlations between the different portions of the gravitational wave signal and the behavior of the horizons, and the precise demarcations between the three regimes are yet to be fully quantified. A full understanding of these correlations is obviously necessary to have a complete picture of a binary black hole merger (see e.g. Gupta et al. (2018); Jaramillo et al. (2012, 2011); Kamaretsos et al. (2012a, b); Bhagwat et al. (2018)). It is also of interest to understand further quantitative features of the merger, such as the evolution of physical quantities across the merger. This includes, among other things, the fluxes of energy and angular momentum, and the evolution of higher order multipoles during the merger. These might be correlated with interesting features of the radiative multipoles found in Borhanian et al. (2019). Numerical simulations are capable of solving the Einstein equations with high accuracy for binary black hole mergers (see e.g. Pretorius (2005a); Campanelli et al. (2006); Baker et al. (2006); Szilagyi et al. (2009)). Such simulations provide an obvious avenue for exploring such questions.

To understand the correlations between the gravitational wave signal and the black hole horizons, we need to first decide precisely what geometrical quantities on the horizon should be considered. In fact, we need to go a further step backwards and decide what kind of horizons should be considered. There are two different ways of visualizing horizons using either event horizons or marginally trapped surfaces. Both of these descriptions are in good agreement in the inspiral and ringdown regimes, but differ substantially during the merger where non-linear and non-perturbative effects of general relativity are especially important. Consider first the event horizon description. An event horizon is the boundary of the region which is causally disconnected from an asymptotically far away observer. It is clear that locating an event horizon requires knowledge of the global properties of the spacetime infinitely far into the future. It is possible, though not trivial, to locate event horizons in numerical binary black hole simulations Hughes et al. (1994); Diener (2003); Anninos et al. (1995); Thornburg (2007), and this yields the well known “pair of pants” picture Matzner et al. (1995). The cross-sections of the “pair of pants” corresponds with the expectations described above. At early times, the cross-section of the event horizon consist of two disjoint surfaces corresponding to the two separate black holes, and a single spherical surface at the end. There are several interesting features of the event horizon in the merger, including the existence of a toroidal phase early in the merger and the non-differentiability of the event horizon Husa and Winicour (1999); the non-differentiability is in fact a general feature of event horizons Chrusciel et al. (2002); Chrusciel and Galloway (1998).

The “pair of pants” picture is intuitively appealing and moreover it seems to provide a complete picture of the black hole merger in accordance with our physical expectations. In reality however, this picture is not so useful, both as a matter of principle and therefore also for any detailed quantitative studies. The problems can be traced back to the global and teleological nature of event horizons: to locate them, one needs to know what happens in the spacetime far in the future. In perturbative situations and when the end-state is known or assumed, it is indeed possible to obtain expressions for the fluxes of energy and angular momentum through the event horizon Hawking and Hartle (1972). In general dynamical situations however, this is not true. There are simple examples, even in spherical symmetry, when the area of the event horizon grows without any corresponding flux of energy Ashtekar and Krishnan (2004). Due to these teleological properties, there is no possible local expression of general validity for, say, the fluxes of energy and angular momentum through event horizons. It is thus not clear how to carry out the program of understanding the merger and relating it to gravitational wave observations outlined at the beginning of the previous paragraph. As a side remark, the teleological property also makes it difficult to locate event horizons in numerical simulations in real time, but in any case, it is certainly possible to locate them once the simulations are complete.

There is an alternate way of visualizing a binary black hole merger which, for both conceptual and practical reasons, is of much greater importance in numerical simulations. The starting point is an unusual property of certain surfaces in the black hole region, first pointed out by Penrose Penrose (1965). This requires the notion of the expansion $\Theta$ of a congruence of light rays; $\Theta$ is the logarithmic rate of change of an infinitesimal cross-section transverse to the null geodesics. A round sphere in flat space has $\Theta>0$ for the outgoing light rays and $\Theta<0$ for the ingoing ones. In the black hole region, there exist spheres (the trapped surfaces) for which both sets of light rays have negative expansion. The outermost such sphere at any given time has vanishing outgoing expansion; these are the marginally trapped surfaces. In stationary situations such as for a Schwarzschild or Kerr black hole, cross-sections of the event horizon are also marginally trapped surfaces, but this correspondence is not true in non-stationary situations. Thus, cross-sections of the event horizon are marginally trapped surfaces very early in the inspiral regime or at very late times. At intermediate times, especially near the merger, the two notions are very different. Furthermore, unlike event horizons, marginal surfaces are not teleological and can be located at any given time without reference to any future properties of spacetime. It is possible to define physical quantities such as mass, angular momentum, multipole moments, and fluxes of energy and angular momentum quasi-locally, i.e. on the marginal surfaces. For this reason, marginal surfaces are widely used in numerical simulations when referring to the properties of black holes. There is a large literature on these quasi-local definitions and their applications to various problems in classical and quantum black hole physics (see Ashtekar and Krishnan (2004); Booth (2005); Faraoni and Prain (2015); Krishnan (2014) for reviews).

Despite this progress, there is still a missing ingredient, namely a unified treatment of inspiral, merger and ringdown. Thus far, all studies of binary black hole coalescence using marginal surfaces have considered the pre- and post-merger regimes separately. The reason for this is that, until recently, it was not known how marginal surfaces behave across the merger; near the merger the marginal surfaces are extremely distorted and previous numerical methods were not successful in tracking such highly distorted surfaces. Using improved numerical methods Pook-Kolb et al. (2019a), we have recently shown the first evidence for the existence of a continuous sequence of marginal surfaces which interpolates between the two disjoint initial black holes and the single final remnant black hole Pook-Kolb et al. (2019b). This is the analog of the “pair of pants” picture for event horizons. In the present work, with further improvements in numerical methods for locating marginal surfaces, we shall provide further unambiguous evidence for this scenario. We shall also show the existence of marginal surfaces with self-intersections. In a companion paper we shall study physical characteristics of the world-tube of marginal surfaces, which is the other important ingredient for physical applications.

The scenario we obtain for the merger is summarized in Fig. 1. The details showing how these results are obtained will be explained in the next sections. The figure shows four snapshots of the MOTSs at various times111We define the factor $\mathcal{M}:=M_{\rm ADM}/1.3$ to be able to state our coordinate quantities in terms of the ADM mass, which in our simulations was chosen to be $1.3$.

in a head-on binary black hole merger starting with Brill-Lindquist initial data. We initially have only the two individual MOTSs without a common horizon. As the black holes get closer, a common MOTS is formed which immediately bifurcates into outer and inner portions visible in the second snapshot. The outer portion loses its distortions as it approaches its equilibrium state, while the inner MOTS becomes increasingly distorted. At some point, just shortly after the third snapshot, the two individual MOTSs touch each other exactly at the time when they merge with the inner common MOTS. After this merger, the two individual MOTSs go through each other. Surprisingly, it turns out that the inner common MOTS continues to exist after the merger and now has self-intersections as shown in the last snapshot. The remainder of this paper will be devoted to explaining how we arrive at this result. A detailed study of the physical aspects of this scenario will be presented elsewhere.

Sec. II summarizes the basic definitions and concepts that we shall need for this paper. The improved numerical algorithm for locating marginal surfaces is described in Sec. III and Sec. IV shows various numerical tests to validate the method. Sec. V discusses our modifications to the numerical methods used to evolve Cauchy data using the Einstein equations. These modifications allow us to reach the required numerical accuracy and convergence, and to carry out our simulations more efficiently. Sec. VI puts together all these ingredients and presents our main results. For a particular initial configuration (the head on collision of comparable mass non-spinning black holes), the merger of marginally trapped surfaces is demonstrated with high numerical accuracy. The merger involves the formation of a marginally trapped surface with self-intersections, showing topology change in a binary black hole merger.

II Marginally outer trapped surfaces

Let $\ell^{a}$ be a congruence of future directed null geodesics, and let $n^{a}$ be another such congruence satisfying $\ell^{a}n_{a}=-1$. Let $q_{ab}$ be the Riemannian metric in the 2-dimensional space transverse to both $\ell^{a}$ and $n^{a}$. The divergence of $\ell^{a}$ and $n^{a}$ are respectively

[TABLE]

Let $\mathcal{S}$ be a closed spacelike 2-surface with null normal fields $\ell^{a}$ and $n^{a}$ respectively. We assume that it is possible to assign outgoing and ingoing directions on $\mathcal{S}$, and by convention, $\ell^{a}$ and $n^{a}$ are the outgoing and ingoing null normals respectively. The classification of $\mathcal{S}$ based on conditions on the expansions are the following:

• •

Trapped: $\Theta_{(n)}<0$, $\Theta_{(\ell)}<0$

• •

Un-trapped: $\Theta_{(n)}<0$, $\Theta_{(\ell)}>0$

• •

Marginally trapped: $\Theta_{(n)}<0$, $\Theta_{(\ell)}=0$

• •

Marginally outer-trapped: $\Theta_{(\ell)}=0$ (no condition on $\Theta_{(n)}$)

All of these refer to future-directed $\ell^{a}$. Thus we should say future-trapped rather than just trapped, but we shall only consider future directed cases. The most important case for us is the marginally outer trapped surface (MOTS) lying within a spatial slice $\Sigma$.

As mentioned in the introduction, there is a large literature on the application of MOTSs to study black holes in various contexts (see e.g. Ashtekar and Krishnan (2004); Booth (2005); Gourgoulhon and Jaramillo (2006); Hayward (2004); Jaramillo (2011); Krishnan (2014, 2008)). They are regularly used in numerical relativity simulations to compute physical quantities Dreyer et al. (2003); Schnetter et al. (2006); Gupta et al. (2018), and this formalism leads naturally to various versions of quasi-local black hole horizons.

While we shall not delve into the mathematical and physical characteristics of MOTSs here, it shall be useful to understand the stability operator for a MOTS and its relevance for time evolution. For a given MOTS $\mathcal{S}$ consider a smooth one-parameter family of closed spherical surfaces $\mathcal{S}_{\lambda}$ which are variations of $\mathcal{S}$ in the normal direction Newman (1987) within the spatial hypersurface $\Sigma$.

On each $\mathcal{S}_{\lambda}$, just as for $\mathcal{S}$, we can define the null normals and calculate the expansion $\Theta_{(\ell)}(\lambda)$, which will of course generally not vanish. The differentiation of $\Theta_{(\ell)}(\lambda)$ leads to an operator $L$ on $\mathcal{S}$:

[TABLE]

Here $r^{a}$ refers to the unit outward pointing spacelike normal to $\mathcal{S}$ (within $\Sigma$) and $f$ is a scalar function on $\mathcal{S}$. Along the 1-parameter family $\mathcal{S}_{\lambda}$, every point on $\mathcal{S}$ traces out a curve with tangent vector $fr^{a}$. The variation of the expansion, i.e. the left hand side of the above equation, is the derivative of the expansion along these curves. This procedure defines an elliptic operator $L$ on a MOTS and the precise expression for $L$ can be worked out. Generically it is of the form

[TABLE]

Here $\Delta$ is the Laplace-Beltrami operator on $\mathcal{S}$ compatible with $q_{ab}$, $\gamma^{a}$ is a vector field on $\mathcal{S}$ related to black hole spin, and $\beta$ is a scalar related to the intrinsic (two-dimensional) Ricci scalar of $\mathcal{S}$. Thus, $L$ is not necessarily a self-adjoint operator due to the presence of $\gamma^{a}$, and its eigenvalues are not necessarily real. Nevertheless, its principal eigenvalue $\Lambda_{0}$, i.e. the eigenvalue with the smallest real part is indeed real. In this paper we shall restrict ourselves to non-spinning black holes with vanishing $\gamma^{a}$ so that all eigenvalues are real.

The primary utility of $L$ is that it determines the behavior of $\mathcal{S}$ under time evolution. It was shown that if the principal eigenvalue is positive, then the MOTS evolves smoothly in time Andersson et al. (2005, 2008, 2009). This stability condition is equivalent to saying that an outward deformation of $\mathcal{S}$ makes it untrapped which is what we expect to happen for the apparent horizon. While not emphasized in Andersson et al. (2005, 2008, 2009), the condition for the existence of $\mathcal{S}$ under time evolution is the invertibility of $L$. Thus, if 0 is not in the spectrum of $L$, then $\mathcal{S}$ continues to evolve smoothly. In the case when $\Lambda_{0}<0$ (which will happen in our case), we must consider the next eigenvalue $\Lambda_{1}$ and check that it does not vanish. See e.g. Booth et al. (2017); Sherif et al. (2018); Mach and Xie (2017) as examples of studies which consider this notion of stability in specific examples.

III Numerical methods for locating highly distorted MOTSs

Consider a Cauchy surface $\Sigma$ on which we wish to locate a MOTS $\mathcal{S}$. Let $\Sigma$ be equipped with a Riemannian metric $h_{ij}$ with the associated Levi-Civita connection $D_{a}$, and let the extrinsic curvature of $\Sigma$ be $K_{ij}$. Let $r^{a}$ be the unit-spacelike normal to $\mathcal{S}$ within $\Sigma$ and let $\tau^{a}$ be the unit-timelike normal to $\Sigma$. Then, a suitable choice of null-normals to $\mathcal{S}$ is

[TABLE]

The condition $\Theta_{(\ell)}=0$ is rewritten as

[TABLE]

This is the equation that we must solve to find $\mathcal{S}$. The conventional approach Thornburg (2007, 2004) assumes that the surface is defined by a level-set function

[TABLE]

where $(r,\theta,\phi)$ are spherical coordinates on $\Sigma$. This assumes that $\mathcal{S}$ is star-shaped with respect to the origin in the chosen coordinate system. In other words, any ray drawn from the origin must intersect the surface only once. This assumption will not hold for the surfaces of interest for us. A variant of this method was proposed in Pook-Kolb et al. (2019a) and shown to be capable of locating extremely distorted surfaces. This new method is based on using a reference surface $\sigma_{R}$, and representing $\mathcal{S}$ in terms of distances $h(\lambda,\mu)$ from $\sigma_{R}$, where $\lambda,\mu$ parameterize $\sigma_{R}$. As long as the reference surface is chosen appropriately, the method can be used to locate almost arbitrarily distorted surfaces. For example, in a numerical evolution, one could choose $\sigma_{R}$ to be the MOTS located in the previous time step. The problem of locating $\mathcal{S}$ then translates to solving a nonlinear partial differential equation for the horizon function $h$. This can be done e.g. via a pseudospectral method, which is what we chose.

For our present application, we have implemented two additional features compared to what was used in Pook-Kolb et al. (2019a). These features are meant to deal with two additional complications that we must necessarily deal with: i) surfaces which have a very narrow “neck” (almost like a figure-eight), and in some instances have features like cusps and self intersections. For this purpose, motivated by the methods used in Jaramillo et al. (2009), we employ bi-spherical coordinates Ansorg (2005). ii) Unlike in Pook-Kolb et al. (2019a) where the MOTS finder was applied to analytical initial data, we now have to deal with numerically generated data on a finite mesh. This requires the use of interpolation schemes some of which were already used in Pook-Kolb et al. (2019b). We now describe in turn both of these additional features. We shall still be restricted to axisymmetry in this work, reducing the task of finding the horizon function $h$ to a one-dimensional problem. However, no in-principle difficulties are foreseen for general non-axisymmetric cases.

III.1 Bi-spherical coordinates

For axisymmetric surfaces, choosing the symmetry axis to be the $z$ axis, we can restrict ourselves to the $(x,z)$ plane and it is often convenient to characterize any point using polar coordinates, i.e. using the distance from the origin and the angle of the position vector with the $z$ axis. However these coordinates are not optimal for describing surfaces with a very narrow neck connecting two spherical portions, i.e. close to a figure-eight in shape. We use instead the bipolar coordinates $(s,t)$ which are based on two foci located at $x=0$, $z=c\pm a$:

[TABLE]

The $(s,t)$ coordinates make the highly distorted inner common MOTS $\mathcal{S}_{\rm inner}$ much easier to parameterize.

Examples demonstrating the effect of this coordinate transformation for three different simulation times are shown in Fig. 2. The three snapshots are at times i) $T=3\,\mathcal{M}$ which is a bit after the top right panel of Fig. 1 and $\mathcal{S}_{\rm inner}$ does not have extreme distortions; ii) $T=5.5\,\mathcal{M}$, shortly before the bottom left panel in Fig. 1 where $\mathcal{S}_{\rm inner}$ has a very narrow neck, and finally iii) $T=6.5\,\mathcal{M}$, a little bit before the bottom right panel of Fig. 1, and $\mathcal{S}_{\rm inner}$ has self-intersections.

The bi-spherical coordinates are employed only for $\mathcal{S}_{\rm inner}$; none of the other horizons have the narrow neck and these coordinates are unnecessary to locate them. To determine the value of $c$ in (7), we first find the two individual MOTSs $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ and choose $c$ to lie in the coordinate center between the lowest point of $\mathcal{S}_{1}$ and the upper-most point of $\mathcal{S}_{2}$. As detailed below, we find the various MOTSs in a series of time slices produced by the numerical simulation. During this tracking of $\mathcal{S}_{\rm inner}$, we numerically approximate the optimal value for $a$ as a post-processing step once the MOTS is located. In practice, this is done by representing $\mathcal{S}_{\rm inner}$ in bi-spherical coordinates and expressing the coordinate functions $s(\lambda),t(\lambda)$ as a truncated series of sines and cosines, respectively, which have the correct symmetry for the problem. We use a slightly lower number of basis functions than necessary to obtain convergence and check the residual expansion of the now imperfect representation. Varying the parameter $a$, we repeat this process to find the value resulting in the lowest residual. The value for $a$ determined this way is then used for finding the MOTS in the next slice, assuming the optimal parameter varies slowly with simulation time.

A further optimization is to re-parameterize the reference surface $\sigma_{R}$ prior to finding the MOTS. A natural choice of parameterization would be the proper length or proper length in coordinate space, the latter obviously being better suited for our numerical representation of the surface. If the curve representing $\sigma_{R}$ in coordinate space is $\lambda\mapsto\gamma_{R}(\lambda)$, this would mean that $\|\gamma_{R}^{\prime}(\lambda)\|_{2}\equiv\text{const}$. However, we obtained faster convergence by taking a non-constant speed function such that $\|\gamma_{R}^{\prime}(\lambda)\|_{2}$ is roughly222We smoothen the speed function along the MOTS by exponentially damping the coefficients of a cosine series representation. This reduces higher frequencies in the density of collocation points along $\mathcal{S}$.

proportional to $1/k_{AB}k^{AB}$, where $k_{AB}$ is the second fundamental form of $\sigma_{R}$ embedded in coordinate space.

Utilization of bi-spherical coordinates together with the above re-parameterization has led to convergent solutions $\mathcal{S}_{\rm inner}$ with about one order of magnitude fewer collocation points compared to the previous method.

III.2 Interpolating numerical data

In each time step, our axisymmetric numerical simulations produce data on a 2-dimensional grid of points lying equidistant in the $(x,z)$ coordinate plane. However, for the nonlinear search for a MOTS $\mathcal{S}$, the expansion $\Theta_{(\ell)}$ and its derivatives have to be computed on a set of points $x_{n}\in\mathbb{R}^{2}$ along trial surfaces $\mathcal{S}^{i}$, c.f. Pook-Kolb et al. (2019a), Section III.B. This requires evaluating the components of the metric $h_{ij}$, its first and second spatial derivatives, the extrinsic curvature $K_{ij}$ and its first spatial derivatives at the points $x_{n}$ which generally do not coincide with any of the grid points of the simulation.

In Pook-Kolb et al. (2019b) we used $4$th order accurate $5$-point Lagrange interpolation. Derivatives were obtained by evaluating $4$th order accurate finite differencing derivatives using the data on the grid and then interpolating the results using $5$-point Lagrange interpolation. For the present paper, however, we switched to quintic Hermite interpolation, which allows us to control the values along with first and second derivatives of the interpolant at the grid points. These derivatives are evaluated using $6$th order accurate finite differencing. Derivatives between the grid points are then computed by analytically differentiating the interpolating polynomial. The advantage is that now first and second derivatives are continuous throughout, which is not the case with Lagrange interpolation.

Interpolation of discrete data will be more accurate with increased grid resolution. However, it will never be exact and even floating point accuracy cannot be neglected, especially near the punctures at computationally feasible resolutions. These additional inaccuracies may limit the numerical convergence as they move the plateau we see below in Fig. 7 up—for example when moving closer to the punctures or reducing the grid resolution—or down. To account for this effect while tracking a MOTS through simulation time, we compute the expansion between the collocation points each time the expansion drops below a pre-set tolerance at the collocation points. After this, we increase the spectral resolution and continue until the tolerance is met at the now larger set of collocation points. This is repeated until the expansion between the collocation points no longer improves, signaling that we have reached the plateau.

A second criterion for stopping to increase the spectral resolution is derived from the absolute values of the coefficients $a_{n}$ of the spectral representation of the horizon function $h$. In a pseudospectral method using a basis of cosines, one expects these coefficients to fall off exponentially for large $n$ if the solution exists. We hence stop increasing the resolution if sub-exponential fall-off of the $a_{n}$ is found following a region of exponential convergence. This prevents our code from overfitting $\mathcal{S}$ to features introduced by the interpolation method, which happens especially for lower resolution simulations.

IV Validating the MOTS finder

With the addition of numerical simulations, the task for our MOTS finder has become more general compared to the purely time-symmetric cases considered in Pook-Kolb et al. (2019a). Therefore, and in light of the surprising result of a self-intersecting MOTS, it is important to validate the method and test it for correctness in an analytic case where the result is known. We shall later present convergence results for further validation.

For this purpose we construct a non-time-symmetric slice with analytically known horizon shape. We choose a slice of the Schwarzschild spacetime in Kerr-Schild coordinates Matzner et al. (1999), i.e.

[TABLE]

where $\delta_{ij}$ is the flat metric, $x_{i}$ are the standard Cartesian coordinates for the flat metric, and we shall often use $(x,y,z)$ instead of $x_{i}$ when no confusion can arise. For Schwarzschild, the radial coordinate is just $r^{2}=x^{2}+y^{2}+z^{2}$. These data have nontrivial extrinsic curvature with the horizon being located at $r=2m$.

To make the horizon non-star-shaped and thus the task more difficult (but still axisymmetric), we transform the coordinates $(x,z)\rightarrow(\bar{x},\bar{z})$ via

[TABLE]

These equations are used to sample $h_{ij}$ and $K_{ij}$ on grids of various resolutions from $1/h=30$ to $1/h=1920$. We choose a reference shape that is close but not identical to the horizon. The MOTS $\mathcal{S}$ and the reference shape $\sigma_{R}$ are shown in the first panel of Fig. 3. For this test we compute the area $A$ of $\mathcal{S}$ and compare it to the exact area $A_{\rm exact}=16\pi m^{2}$, where $m=1$. We also compute the maximum coordinate distance $\|\mathcal{S}-\mathcal{S}_{\rm exact}\|_{\infty}$ of the numerical solutions to the exact horizon. The second panel demonstrates that our numerical solutions converge to the expected solutions as the resolution of the numerical grid is increased.

V The numerical evolutions

V.1 Formulations, Discretization, and Implementation

We set up initial conditions for the spacetime geometry as two puncture black hole using the method of Brill and Lindquist Brill and Lindquist (1963). To evolve the geometry, we use the BSSN formulation of the Einstein equations with a $1+\log$ slicing and a $\Gamma$-driver shift condition Alcubierre et al. (2000, 2003). We also impose axisymmetry throughout the calculation.

For our setup (see below), we choose a domain with $x\in[0;10]$, $z\in[-10;10]$, and $T\in[0;7]$. (Due to axisymmetry, we only consider the hyperplane $y=0$.) For simplicity, we use Dirichlet boundary conditions to set all time derivatives to zero at the outer boundary. We check that the errors introduced by the artificial boundary conditions do not affect the geometry near the MOTSs.

We choose a Cartesian basis for the tangent space, i.e. we represent vectors and tensors via their $x,y,z$ components. Although axisymmetry requires that certain components or linear combinations of components must vanish, we do not explicitly impose such conditions. Instead, we only impose axisymmetry on spatial derivatives: We require that the Lie derivatives of all quantities in the $\phi$ direction be zero, and we use this to remove all $y$ derivatives. ($y$ derivatives are then either [math], or are replaced by combinations of various $x$ derivatives.) We use l’Hôpital’s rule to regularize these expressions on the axis. This closely follows the approach described in Pretorius (2005b), extended to handle second derivatives as well. The set of expressions for handling first and second $y$ derivatives for all tensor ranks appearing in the BSSN formulation is lengthy, and is available in a Mathematica script as part of Kranc Husa et al. (2006); Kranc .

In our discretization, we also require a small region “on the other side” of the axis (where $x<0$), which we calculate by rotating the region with $x>0$ by $\pi$.

We also experimented with the Cartoon method Alcubierre et al. (2001) to impose axisymmetry. Cartoon uses a spatial rotation in the $\phi$ direction and then spatial interpolation to populate points away from the $y$ axis, so that $y$ derivatives can be calculated in the standard manner. We found that the Cartoon method does not work well with higher order (higher than $4$th) finite differencing: The result of a Lagrange interpolation is not continuous, which leads to large oscillations when derivatives are taken near the axis where the Cartoon rotation angle is large.

In our setup, the punctures are located on the $z$ axis and are initially at $z_{\pm}=\pm 0.65$. The puncture masses are $m_{+}=0.5$ and $m_{-}=0.8$ (i.e. the “upper” black hole is smaller). The punctures have no linear or angular momentum.

Details of initial and gauge conditions are described in Wardell et al. (2016). Our exact parameter settings are available in the parameter files in the repository Pook-Kolb et al. (2019c).

We use $6$th order finite differencing to discretize space. We also add a $6$th order Kreiss-Oliger artificial dissipation, which reduces our spatial accuracy to $5$th order. We use a $6$th order accurate Runge-Kutta time integrator. Our discretization is globally $5$th order accurate, as we demonstrate below in section V.2. We do not use mesh refinement nor multiple grid patches as these would not be beneficial for our calculations that span only a short time and a small region of space, compared to systems of orbiting binary black holes.

Compared to $4$th and $8$th order discretizations, $6$th order is most efficient for us. $4$th order calculations require significantly higher resolutions, and $8$th order calculations are significantly slower since they use larger stencils and require more integrator substeps. $8th$ order calculations also require higher resolutions before their error falls below that of $6$th order calculations.

We perform our calculation via the Einstein Toolkit Löffler et al. (2012); EinsteinToolkit . We use TwoPunctures Ansorg et al. (2004) to set up initial conditions and an axisymmetric version of McLachlan Brown et al. (2009) to solve the Einstein equations, which uses Kranc Husa et al. (2006); Kranc to generate efficient C++ code.

V.2 Accuracy, Convergence

To demonstrate the accuracy of our discretization, we plot in Fig. 4 the Hamiltonian constraint

[TABLE]

on grid points close to the inner common MOTS at two different times for different grid resolutions. Here, $R$ is the Ricci scalar of the slice $\Sigma$. There is no significant difference between the two times. Note that in coordinate space, $\mathcal{S}_{\rm inner}$ lies closer to the punctures in its upper than in its lower half, compare also Fig. 1. In terms of the curve’s proper length parameter $\bar{\lambda}$ (scaled to $\bar{\lambda}\in[0,\pi]$), this corresponds to $\bar{\lambda}\lesssim\pi/2$ and $\bar{\lambda}\gtrsim\pi/2$, respectively, where our representation only covers half of the plotted MOTS (say for positive $x$ values) due to axisymmetry.

The results have been scaled to account for $5$th order convergence. We indeed find $5$th order convergence for $1/h\geq 240$ closer to the punctures and for $1/h\geq 120$ further away from the punctures. In that latter region, the highest resolution results with $1/h=960$ show slightly larger errors than expected from $5$th order accuracy.

This is caused by round-off errors starting to dominate the finite difference derivatives, as is demonstrated in Fig. 5. Here, the different curves represent the results obtained using stencils of $3$ to $9$ points for the derivatives of the metric components, corresponding to $2$nd to $8$th order accuracy. We see the typical behavior of convergence up to the resolution at which the round-off error becomes dominant. This happens at lower resolutions for the higher order methods as these reach the round-off limit earlier. Note that the optimal resolution depends on the function being approximated and in our case becomes larger the closer we get to the puncture. This explains the different behavior in the first and second half of the plots in Fig. 4.

VI The existence of self-intersecting MOTSs

With the technical improvements at hand, we now turn to the main result of this paper, namely the merger of the inner MOTS with the two individual horizons, and the occurrence of self intersecting MOTSs just after this merger (see Fig. 1). We will study a single configuration with high resolution. We focus primarily on numerical accuracy and convergence to confirm the merger scenario and the existence of self-intersecting MOTSs. There are obviously numerous physical and geometrical properties of great interest. First however, we need to prove this scenario numerically beyond any reasonable doubt, which is what we shall do here. A detailed discussion of the interesting physical and geometrical properties of the world tube of MOTSs will be postponed to a forthcoming paper. Similarly, we shall not discuss here the various extensions to non-time symmetric and non-axisymmetric data. As mentioned previously, we start with Brill-Lindquist initial data with the bare masses $m_{+}=0.5$ and $m_{-}=0.8$. The initial coordinate separation between the punctures is $1.3\,\mathcal{M}$ (i.e. $1$ in units of the total ADM mass $M_{\rm ADM}=m_{+}+m_{-}$). Simulations are performed at various grid resolutions: $1/h=60$, $120$, $240$, $480$, $960$. We have already shown in the previous section that the numerical solution of the Einstein equations for the given initial data is sufficiently accurate and all constraint violations converge at the expected rate when $h$ is varied. Given this numerical spacetime, we can use our horizon finder to locate the various MOTSs. It remains to be shown now that the surfaces thus found are indeed MOTSs.

Before proceeding further, it might be useful to clarify the nature of the MOTS with self-intersections shown in the bottom right panel of Fig. 1. Viewed as a submanifold of the 3-dimensional Riemannian spatial slice $\Sigma$, this manifold might appear to be non-differentiable at the point of self-intersection and one might be concerned that there is no well defined normal to the manifold at that point (and hence no well defined expansion either). This is however incorrect, and formally the curve is simply understood as an immersion instead of an embedding. In the present case, because of axisymmetry, we can restrict ourselves to a two-dimensional section (say the $x$-$z$ plane as we have been using so far). Then the horizon is simply a parameterized curve, i.e. a mapping of the circle $S^{1}$ into $\Sigma$, $f:S^{1}\rightarrow\Sigma$ (this is precisely how this curve is defined numerically). Using the map $f$, we can push forward tangent vectors to $\Sigma$ and thus we have well defined normals depending on which direction one traverses the point of self-intersection (see Fig. 6). The relevant topological property of the curve is the winding number, i.e. the number of rotations that a tangent vector undergoes when we go all the way around the curve; each loop adds +1 to the winding number. Curves with different winding numbers cannot be smoothly deformed into each other Whitney, Hassler (1937). This is why in order to get the self-intersections, it is necessary to go through the cusp at the merger. In non-axisymmetric situations, we have to necessarily deal with mappings of $S^{2}$ into the 3-manifold $\Sigma$ (which are in fact simpler Smale (1959)), but we shall not discuss this here.

VI.1 Convergence

Except for the modifications introduced earlier in Sec. III, we employ the same basic Newton-Kantorovich search as in Pook-Kolb et al. (2019a) with each step being performed using a pseudospectral method. If the nonlinear search converges, we expect the exponential convergence of the individual pseudospectral steps to carry over to the solution of the full nonlinear problem.

This is indeed the case, as can be seen in Fig. 7. It shows the maximum residual expansion between the collocation points for $\mathcal{S}_{\rm inner}$ at two different times of the simulation: one at $T=5.35\,\mathcal{M}$, where the MOTS is already highly distorted, and one at $T\approx 5.7333\,\mathcal{M}$. This second case is after the individual MOTSs touch. At this stage, $\mathcal{S}_{\rm inner}$ lies in the inside of $\mathcal{S}_{1}\cup\mathcal{S}_{2}$ and intersects itself. There is no qualitative difference in convergence and the plateau is approximately at the same level for the same grid resolution.

We also see in Fig. 7 that the negative slope continues into the plateau region. This effect is more pronounced for lower grid resolutions and not noticeable for $1/h=960$. It is caused by fitting the horizon to features introduced by the interpolation. We avoid this unphysical effect in practice by limiting the pseudospectral resolution as described at the end of Sec. III.2.

Instead of varying the pseudospectral resolution, we can test convergence for different grid resolutions $1/h$ of the simulation. The quantity we use here is the convergence of the coordinate shapes of the curves representing the MOTSs. Fig. 8 shows that we indeed find convergence of the shapes.

We show in Fig. 9, as a function of time, the residual expansion of the various MOTSs for the highest resolution that we have considered, namely $1/h=960$. The residual expansion is one of the key ingredients which gives us confidence that the surfaces we find are indeed MOTSs. Note first that for all the “easy” cases, namely for the two individual MOTSs $\mathcal{S}_{1,2}$ and for the apparent horizon, the residual expansion is no more than $\mathcal{O}(10^{-11})$. These horizons do not have any portions with extreme curvatures and there is no difficulty in locating them. In fact, the residual expansion is largest for the smaller horizon, and is $\mathcal{O}(10^{-12})$ for the larger horizon and the apparent horizon. The difficult case is of course the inner common horizon, which required the various technical improvements detailed earlier. The most difficult cases are those which have the narrow neck and correspondingly highly curved portions. There is a small duration of time near $T_{\rm touch}$ where we are not able to locate $\mathcal{S}_{\rm inner}$. At all the other times shown in the plot, the residual expansion is no more than $\mathcal{O}(10^{-9})$. In fact, away from $T_{\rm touch}$, the residual expansion is as good as for the other MOTSs. In particular, this is true after $T\sim 5.7\,\mathcal{M}$. At these times $\mathcal{S}_{\rm inner}$ has developed self-intersections. Thus, our confidence in the existence of self-intersecting MOTSs is the same as our confidence in the existence of the other MOTSs, which of course are already well established.

VI.2 Area and stability

Some quantitative numbers for this evolution are:

• •

The common horizon forms at $T_{\rm common}\approx 1.37460222\,\mathcal{M}$.

• •

The two individual horizons touch at $T_{\rm touch}\approx 5.5378176\,\mathcal{M}$.

• •

The area of the inner horizon reaches a minimum at $T_{\rm min}\approx 5.50592\,\mathcal{M}$, i.e. just a little bit before $T_{\rm touch}$. This behavior of $\mathcal{S}_{\rm inner}$ was previously noted in Pook-Kolb et al. (2019b).

These values were computed at the various resolutions up to $1/h=960$ and converge up to the shown number of decimal places; compare also Fig. 10.

The areas of the various horizons are plotted as functions of time in Fig. 11. The bottom-right panel presents a useful picture of the merger process. It shows the areas of the apparent horizon, the inner common horizon and the sum of the areas of the individual horizons. It shows the formation and bifurcation of the apparent horizon and it also shows the merger, i.e. the crossing of the curves for the inner horizon and the individual horizons.

The principal eigenvalue of the stability operator for the various horizons is shown in Fig. 12. We see that $\Lambda_{0}$ is always positive for $\mathcal{S}_{1,2}$ and for the apparent horizon, and that it is not strongly varying. $\mathcal{S}_{\rm inner}$ is more interesting. When it is initially born, it coincides with the apparent horizon and has $\Lambda_{0}=0$. At all subsequent times, $\mathcal{S}_{\rm inner}$ has $\Lambda_{0}<0$; to understand its stability we need to consider the next eigenvalue $\Lambda_{1}$. But already from Fig. 12, we see interesting behavior of $\Lambda_{0}$ for $\mathcal{S}_{\rm inner}$, namely a cusp at $T_{\rm touch}$. Fig. 13 shows $\Lambda_{1}$ for the inner horizon, and it is seen to be positive thus demonstrating stability. Again, we see a cusp-like behavior near $T_{\rm touch}$.

VII Conclusions

In this paper we examined in detail the scenario for the merger of MOTSs outlined previously in Pook-Kolb et al. (2019b). We have done this by evolving a particular Brill-Lindquist setup and finding all MOTSs at various times. We have tracked the inner common horizon with high accuracy. In particular, we present strong numerical evidence that the inner horizon merges with the two individual horizons precisely at the time when they touch. Moreover, we find that the inner horizon develops self-intersections just after the merger. This provides then a connected sequence of MOTSs taking us from the two disjoint initial horizons to the final apparent horizon. We have also studied some basic properties of the MOTSs including their area and stability. There are numerous other interesting physical and geometric properties of the world tube of MOTSs which shall be studied in detail in forthcoming work.

Acknowledgements.

We thank Abhay Ashtekar, Bernd Brugmann, Luis Lehner, and Andrey Shoom for valuable discussions and comments. We are especially grateful to Jose-Luis Jaramillo for extensive discussions and for suggesting the use of bipolar coordinates. The MOTS finder Pook-Kolb et al. (2019c) used in this research is developed in Python with SimulationIO Schnetter and Miller (2019) being used for reading the numerical simulation data. The libraries SciPy Jones et al. (2001–), NumPy van der Walt et al. (2011), mpmath Johansson et al. (2018), SymPy Meurer et al. (2017) and Matplotlib Hunter (2007); Droettboom et al. (2018) were used for certain numerical, validation and plotting tasks. O.B. acknowledges the National Science Foundation (NSF) for financial support from Grant No. PHY-1607520. This research was also supported by the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade. This research was enabled in part by support provided by SciNet (www.scinethpc.ca) and Compute Canada (www.computecanada.ca). Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium Loken et al. (2010). SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund – Research Excellence; and the University of Toronto.