Curvature Invariants in a Binary Black Hole Merger

Black holes are solutions of general relativity and are most naturally characterized by their event horizon. The event horizon of a black hole (BH) is defined as the boundary of the causal past of future null infinity. Intuitively, this means that on the inner side of the event horizon, light cannot escape to null infinity. Notice that event horizons require knowledge of the global structure of spacetime [8, 9, 10]. However, for numerical relativity it is more convenient to use an initial value formulation of GR (a 3+1 approach), where initial data is given on a Cauchy hypersurface and is then evolved forward in time. This approach requires an alternative description of BH horizons which is not dependent of the BH’s future. [61, 60, 3, 4, 36, 29, 31, 32, 33].

Let $\mathcal{S}$ be a compact 2D surface without border and of spherical topology, and consider light rays leaving and entering $\mathcal{S}$, with directions $l$ and $n$, respectively. Let $q_{ab}$ be the induced metric on $\mathcal{S}$ and denote the respective expansions as $\Theta_{(l)}=q^{ab}\nabla_{a}l_{b}$ and $\Theta_{(n)}=q^{ab}\nabla_{a}n_{b}$ [55]. Then, $\Theta_{(l)}$ and $\Theta_{(n)}$ are quantities which are positive if the light rays locally diverge, and negative if the light rays locally converge, and are zero if the light rays are locally parallel. We say that $\mathcal{S}$ is a trapping surface if $\Theta_{(l)}<0$ and $\Theta_{(n)}<0$ [8, 48, 53, 50]. Define $\mathcal{S}$ to be a marginally outer trapped surface (MOTS) if it has zero expansion for the outgoing light rays, $\Theta_{(l)}=0$ [53, 50, 51, 57, 27, 30]. ($\mathcal{S}$ is a future MOTS if $\Theta_{(l)}=0$ and $\Theta_{(n)}<0$ and a past MOTS if $\Theta_{(l)}=0$ and $\Theta_{(n)}>0$ [51]). MOTSs turn out to be well-behaved numerically, and can be used to trace physical quantities of a BH as they evolve over time and through a BBH merger [57, 27, 30, 53, 54].

In practice, it is common to view MOTSs as contained in a given 3D Cauchy surface. Within such a surface, the outermost MOTS is called the apparent horizon (AH) [53, 50, 51, 57, 27, 30]. AHs have many applications to numerical relativity, since tracking an AH only requires knowledge of the intrinsic metric $q_{ab}$ restricted to the spacetime hypersurface and the extrinsic curvature of that hypersurface at a given time [28, 10, 27]. For example, AHs are useful to study gravitational waves, as gravitational fields at the AH are correlated with gravitational wave signals [28, 39, 38, 35, 30, 55, 34]. AHs are also used to numerically simulate binary black hole (BBH) mergers and the collapse of a star to form a BH [13]. As another example, AHs play a role in checking initial parameters and reading off final parameters of Kerr black holes in gravitational wave simulations at LIGO [13, 2, 1]. One possible disadvantage of AHs is that the definition of AHs as the ”outermost MOTS” relies on the given foliation of the spacetime into Cauchy surfaces [7, 57].

If one smoothly evolves a given MOTS forward in time, one obtains a world tube which is foliated by these MOTS. This world tube is known as a dynamical horizon (DH) [14, 8, 9, 10]. One application of DHs is that they could contribute to our understanding of BH formation [8, 9, 10, 13]. As is the case with AHs, DHs are dependent on the foliation of the spacetime into Cauchy surfaces, as this spacetime foliation corresponds uniquely to a DHs with a unique foliation into MOTSs [36, 7, 3, 4]. Furthermore, the above definitions of MOTSs, AHs and DHs serve as a quasi-local description of BHs [22, 21].

It has been conjectured that one can uniquely define a smooth, locally determined and foliation invariant horizon based on the algebraic (Petrov) classification of the Weyl tensor [22, 21]. The necessary conditions for the Weyl tensor to be of a certain Petrov type can be stated in terms of scalar polynomials in the Riemann tensor and its contractions which are called scalar polynomial (curvature) invariants (SPIs). The first aim of this work is to study certain SPIs numerically during a BBH merger.

The Petrov classification is an eigenvalue classification of the Weyl tensor, valid in 4 dimensions (D). Based on this classification, there are six different Petrov types for the Weyl tensor in 4D: types ${\bf I}$, ${\bf II}$, ${\bf D}$, ${\bf III}$, ${\bf N}$ and ${\bf O}$ (which is flat spacetime). One can also use the boost weight decomposition to classify the Weyl tensor, which is equivalent in 4D to the Petrov classification. One can also algebraically classify the symmetric trace free operator, $S_{ab}$, that is the trace free Ricci tensor, which is equivalent to the Segre classification [59].

The boost weight algebraic classification generalizes the Petrov classification to $N$-dimensional spacetimes [22, 21, 23, 24, 15, 45]. In $N$ D, and with Lorentzian signature $(+1,-1,\ldots,-1)$, we start with the frame of $N$–vectors, $\{{\bf l},{\bf n},\{{\bf m}_{i}\}_{i=2}^{N-1}\}$, where ${\bf l}$ and ${\bf n}$ are null and future pointing, ${\bf l}\cdot{\bf n}=1$, and the $\{{\bf m}_{i}\}$ are real, spacelike, mutually orthonormal, and span the orthogonal complement to the plane spanned by ${\bf l}$ and ${\bf n}$. The possible orthochronous Lorentz transformations are generated by null rotations about ${\bf l}$, null rotations about ${\bf n}$, spins (which involve rotations about ${\bf m}_{i}$), and boosts [45]. With respect to the given frame, boosts are given by the transformation:

	$\displaystyle{\bf l}$	$\displaystyle\rightarrow\lambda{\bf l}$
	$\displaystyle{\bf n}$	$\displaystyle\rightarrow\lambda^{-1}{\bf n}$
	$\displaystyle{\bf m_{i}}$	$\displaystyle\rightarrow{\bf m_{i}}$

for all $i\in\{2,\ldots,N-1\}$ and for some $\lambda\in\mathbb{R}\backslash\{0\}$. (The remaining transformations are given in [23, 24, 15, 45].) It is possible to decompose the Weyl tensor into components organized by boost weight [23, 24, 15].

It is of particular interest to know whether a given 4D spacetime is of special algebraic type II or D. We can state the necessary conditions as discriminant conditions in terms of simple $SPI$s [22, 21, 16, 18]. Just as an $SPI$ is a scalar obtained from a polynomial in the Riemann tensor and its contractions [22, 21], an $SPI$ of order $k$ is a scalar given as a polynomial in various contractions of the Riemann tensor and its covariant derivatives up to order $k$ [22, 21]. It turns out that BH spacetimes are completely characterized by their $SPI$s [17]. The necessary discriminant conditions on the 4D Weyl tensor for the spacetime to be of type II/D can be stated as two real conditions and are given in [17].

Contracting the 4D (complex) null tetrad, $({\bf l},{\bf n},{\bf m},{\overline{\bf m}})$ where ${\bf m}$ and ${\overline{\bf m}}$ are complex conjugates, with the Weyl tensor, $C_{abcd}$, one may form the complex scalars, ${\bf\Psi_{0}},{\bf\Psi_{1}},{\bf\Psi_{2}},{\bf\Psi_{3}},{\bf\Psi_{4}}$ and, in terms of these scalars, as in the Newman-Penrose (NP) formalism (discussed later), one may define the scalar invariants:

	$\displaystyle I$	$\displaystyle={\bf\Psi_{0}}{\bf\Psi_{4}}-4{\bf\Psi_{1}}{\bf\Psi_{3}}+3{\bf\Psi_{2}}^{2}$		(1)
	$\displaystyle J$	$\displaystyle=\begin{vmatrix}{\bf\Psi_{4}}&{\bf\Psi_{3}}&{\bf\Psi_{2}}\\ {\bf\Psi_{3}}&{\bf\Psi_{2}}&{\bf\Psi_{1}}\\ {\bf\Psi_{2}}&{\bf\Psi_{1}}&{\bf\Psi_{0}}\end{vmatrix}$		(2)

It can be shown that the two aforementioned real scalar conditions are equivalent to the real and imaginary parts of the following complex syzygy [59]:

$$\mathcal{D}\equiv I^{3}-27J^{2}=0$$

(3)

Thus for Petrov types ${\bf II}$ and ${\bf D}$, equation (3) holds everywhere. It also turns out that for Petrov types ${\bf III}$, ${\bf N}$, and for ${\bf O}$, we have $I=J=0$, so (3) is satisfied trivially.

Having discussed the Petrov and boost weight classifications, we now turn to the Geometric Horizon Conjecture (GHC) in which we define the geometric horizon (GH) as the set on which the $SPI$s, defined in (3), vanish [22, 21]. The level–$0$ sets of these $SPI$s might not form a horizon with nice properties, however, since these $SPI$s could vanish additionally on axes of symmetry or fixed points of isometries [22, 21]. We know from (3) that if the spacetime is algebraically special, then the given complex $SPI$ vanishes. More precisely, the GHC is given as follows [22, 21]:

There are many examples of spacetimes that support the plausibility of the GHC either by explicitly exhibiting GHs or by finding other established BH horizons on which the Weyl and Ricci tensor are algebraically special [44, 20, 8, 9, 10, 22, 21, 19, 59]. For example, in the Kerr spacetime, by invoking the notion of a non-expanding weakly isolated null horizon and an isolated horizon, it can be proven, using the induced metric and induced covariant derivatives on the submanifold and assuming the dominant energy condition, that the Weyl and Ricci tensors are both of type ${\bf II/D}$ on the event horizon. This means that one can extract a subset of the set of points where the Weyl and Ricci tensors are both of algebraic type ${\bf II/D}$, to define a smooth BH horizon, namely the event horizon [5, 41, 8, 9, 10, 22, 21]. It can also be shown that the covariant derivatives of the Riemann tensor are of type ${\bf II}$ on this horizon [22, 21]. The Kerr geometry approximates the spacetime outside the event horizon of a BH formed by a collapsing star. By continuity, the region just inside the event horizon must closely approximate the Kerr geometry and it is suggested in [22] that there is a surface inside this horizon which is smooth and uniquely identified by geometric constraints [8, 9, 10, 5, 41, 22, 21], thereby qualifying as a GH.

Another example to support the GHC comes from a family of exact closed universe solutions to the Einstein-Maxwell equations with a cosmological constant representing an arbitrary number of BHs, discovered by Kastor and Traschen (KT) [40]. Consider the merger of 2 BHs. At early times, there are two 3D disjoint GHs forming around each BH [22, 21]. However, at intermediate times, it turns out that the invariant, $\mathcal{D}=I^{3}-27J^{2}=0$ as in (3) only at the coordinate positions of each of the BHs, along certain segments of the symmetry axis, and along a 2D cylindrical surface, which expands to engulf the 2 BHs as they coalesce [22, 21]. During the intermediate process, there are 3D surfaces located at a finite distance from the axis of symmetry for which the traceless Ricci tensor (and hence the Ricci tensor, $R_{ab}$) is of algebraic type ${\bf II/D}$. There is also evidence of a minimal 3D dynamically evolving surface where a scalar invariant, ${\cal W}_{1}$, assumes a constant minimum value. This suggests that there is a GH during the dynamical regime between the spacetimes [22, 21], but further investigation is needed. At late times, the spacetime then settles down to a type ${\bf D}$ Reissner-Nordstrom-de-Sitter BH spacetime with mass $M=m_{1}+m_{2}$, which turns out to have a GH [44, 20]. So a GH forms at the beginning and end of the coalescence. For further information on the two-BH solution, see [40]. The KT solution for multiple BHs was studied and GHs around each BH were found in [43]. The results were compared with the previously mentioned 2-BH solution. Additionally, three black-hole solutions were studied and GHs were found around these BHs also [22, 21]. For information on more than two BHs, see [47].

There are additional examples of spacetimes that support the GHC by identifying GHs with the level–$0$ sets of $\rho$ and $\mu$. These examples include stationary spacetimes with stationary horizons (e.g. Kerr-Newman-NUT-AdS) [44]; spherically symmetric spacetimes such as vacuum solutions or known exact solutions (e.g. Vaidya or LTB dust models) [22]; quasi-spherical Szekeres spacetimes [19]; and the Kastor–Traschen solution for $N>1$ BHs, as mentioned previously [40]. The authors have also studied vacuum solutions in the case of axisymmetry (so $R_{ab}=0$) and where the Weyl tensor is of algebraic type ${\bf I}$ [25]. Based on the previous examples, it is also natural to study the covariant derivative of the Weyl tensor in this setting and the level–$0$ sets of $\rho$ and $\mu$.

We wish to study the behaviour of the complex $SPI$, $\mathcal{D}$, as defined in (3), through a BBH merger. Since the Kerr geometry is type ${\bf D}$ everywhere, it follows that $\mathcal{D}=0$ everywhere for a Kerr BH. Based on our understanding of the features of gravitational collapse and by a plethora of numerical simulations, it is believed that in a BBH merger the merged BHs at late times settle down to a solution well described by the Kerr metric [22, 21]. Thus, for a merger of 2 initially Kerr BHs, a plot of the real part and imaginary part of $\mathcal{D}$ should be roughly zero everywhere at early and late times. However, in the intermediate “dynamical” region (during the actual merger and coalescence at intermediate times), these same zero plots should yield important information. This is what we wish to study.

We highlight some known features of a binary black hole merger, as described in [53, 54, 28]. This also serves to set up our notation. In [53, 54], it was found that there is a connected sequence of MOTSs, which interpolate between the initial and final states of the merger (two separate BHs to one BH, respectively) [53, 54]. The dynamics are as follows: Initially, there are two BHs with disjoint MOTS (which are AHs at this point [28]), $\mathcal{S}_{1}$, and $\mathcal{S}_{2}$, one around each BH. Then, as the two BHs evolve, a common MOTS forms around the two separate BHs and bifurcates into an inner MOTS, $\mathcal{S}_{i}$, which surrounds the MOTS and an outer MOTS, $\mathcal{S}_{c}$. $\mathcal{S}_{c}$ increases in area and encloses $\mathcal{S}_{1}$, $\mathcal{S}_{2}$ and $\mathcal{S}_{i}$, and is the AH at the time of the merger [28, 53, 54]. The fate of $\mathcal{S}_{c}$ and the bifurcation at the time of the merger is well understood [28, 9, 57, 6, 30, 53, 54, 37, 46, 52].

Instead of studying a head-on collision, in this paper we shall study a quasi-circular orbit of two merging, equal mass and non-spinning BHs. Apart from [49, 25], the analysis of the quantity $\cal{D}$ employed in this simulation is new. In these simulations, the Einstein toolkit infrastructure was used [42] and the simulations are run using $4^{\text{th}}$ order finite differencing on an adaptive mesh grid, with adaptive refinement level of 6 [56, 12]. Brill-Lindquist initial data are used, with BH positions and momenta set up to satisfy the initial conditions necessary for a quasi-circular orbit which evolves for less than $1$ orbit before merging (QC0-initial condition). For more details, see Table I of [26]. Instead of analyzing a sequence of MOTS throughout the merger, we seek to define and study a GH as it evolves through the merger, in accordance with the GHC. Since (3) sets necessary conditions for the Weyl tensor to be of algebraic type II, we seek to analyze the constant contours of the difference $\mathcal{D}=I^{3}-27J^{2}$. In the simulations, the real and imaginary parts of $I$ and $J$ are calculated using the Weyl scalars, $\{{\bf\Psi}_{i}\}_{i=0}^{5}$, as given in equations (1) and (2), and the calculations are carried out using the orthonormal fiducial tetrad, as given by [11]. Note that in the rest of this paper, we will use the notation from [53, 54, 28] to describe the various MOTSs that appear in our simulations. To recapitulate, $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are the 1st and 2nd initial MOTS and $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$ describe the inner and common MOTS as they appear in our simulation, respectively.

Figures 1–4 provide plots of various level sets of $\mathcal{D}_{r}=\text{Re}\{D\}$, $\mathcal{D}_{i}=\text{Im}\{D\}$ and $|\mathcal{D}|$ as functions of $(x,y)\in\mathbb{R}^{2}$ at selected instances of the time parameter, $t$, where $t=0$ indicates the start of the numerical computation. The level–$\varepsilon$ sets of $|\mathcal{D}|,\mathcal{D}_{r},\mathcal{D}_{i}$ in Figures 1–4 are calculated with a fixed spatial coordinate value of $z=0.03125$, as are points displayed in Figure 2. We choose this value of $z=0.03125$ to illustrate most clearly the main features of the plots. However, the MOTSs $\mathcal{S}_{1}\cup\mathcal{S}_{2}$ in Figure 4 are calculated with $z$ lying in a range $z\in[0.02,0.04]$ and $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$ are calculated with $z\in[-0.1,+0.1]$ to accomodate the grid spacing in the spatial coordinates $(x,y,z)$ so that these MOTSs can be displayed fully. In Figure 4, these level sets of $|\mathcal{D}|,\mathcal{D}_{r},\mathcal{D}_{i}$ are also compared with $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$. The full compliment of pictures describing this BBH merger are displayed in [49]. We present a subset of those figures to illustrate the essential features. In each of Figures 1-4, the data corresponding to $x<0$ was obtained by rotating the data corresponding to $x>0$ by $180$ degrees about the $x=y=0$ axis. In Figures 1–4, we plot the centroids and outlines of $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ along with $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$, when they have formed in Figure 4. The centroids have been added as a marker to track the locations of the BHs through the merger and provide a reference against which we can compare our candidate GHs, namely the level–$\varepsilon$ sets of $|\mathcal{D}|$.

Figure 1 provides the contour plots of the magnitude of $\mathcal{D}=I^{3}-27J^{2}$, denoted $|\mathcal{D}|$, on a log scale (see (3)) for $t=8,\;12,\;16,\;18,\;20,\;24,\;26,\;27,\;28$ in the upper left, upper center, upper right, middle left, center, middle right, lower right, lower center and lower right panels, respectively, for fixed $z=0.03125$. Since $|\mathcal{D}|$ is positive definite, the level–$0$ sets of $|\mathcal{D}|$ are impossible to find precisely due to discrete resolution and numerical error. Instead, we highlight the evolution of the level–$\varepsilon$ sets, where $\varepsilon=3\times 10^{-4},\;5\times 10^{-4},\;1\times 10^{-3}$. The overlaid green, red and white contours of each frame are the level–$3\times 10^{-4}$, level–$5\times 10^{-4}$ and level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$, respectively.

At early times (i.e., at $t=8$ and $t=12$ in the upper left and upper center panels, respectively), each of the level–$\varepsilon$ sets are partitioned into pairs of simple closed curves. At $t=16$ (upper right panel), the red level–$5\times 10^{-4}$ set and the white level–$1\times 10^{-3}$ set of $|\mathcal{D}|$ each form a third simple closed curve between $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, which is centred at the origin. At times $t=18$ (middle left panel), $t=20$ (center panel) and $t=24$ (middle right panel), for each respective $\varepsilon=3\times 10^{-4},\;5\times 10^{-4},\;1\times 10^{-3}$, the multiple simple closed curves partitioning the level–$\varepsilon$ set of $|\mathcal{D}|$ have joined so that each level–$\varepsilon$ curve is now a single simple closed curve. At times $t=26,\;27,\;28$, the level–$\varepsilon$ sets of $|\mathcal{D}|$ continue to track the merged BHs and are displayed in the lower three panels of Figure 1. It follows that the level–$\varepsilon$ curves for each $\varepsilon=3\times 10^{-4},\;5\times 10^{-4},\;1\times 10^{-3}$ at each $t$ form an invariantly defined, foliation invariant horizon that contains each separate BH at early times, and contains the merged BH at late times.

The evolution of the level–$\varepsilon$ curves through the quasi-circular BBH merger in Figure 1 is reminiscent of the sequence of MOTS that take place during the head-on collision simulation in [53, 54]. In particular, the bifurcation into $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$ described in [53, 54, 28] can be compared to our present quasi-circular simulations. This comparison is most striking at $t=16$, in the upper right panel of Figure 1, when the white level–$1\times 10^{-3}$ and red level–$5\times 10^{-4}$ sets are partitioned into three simple closed curves. However, our numerical computations are not precise enough to study the details of the bifurcation in [53, 54, 28]. At times $t=24,\;26,\;27,\;28$, in the middle right panel and lower three panels of Figure 1 respectively, it also seems that $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ do not merge fully. Thus, it is possible that at at late times, the level–$\varepsilon$ sets of $|\mathcal{D}|$ for $\varepsilon=3\times 10^{-4},\;5\times 10^{-4},\;1\times 10^{-3}$ may track $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, which have been found in [28] to overlap but not intersect at late times. However, our simulations did not run to late enough times to make this clear.

In any case, Figure 1 provides strong evidence that for each $\varepsilon=3\times 10^{-4},\;5\times 10^{-4},\;1\times 10^{-3}$, the level–$\varepsilon$ sets of $|\mathcal{D}|$ track a unique GH, which can be identified by the level–$0$ set of $\mathcal{D}$. It is of interest to study the level–$\varepsilon$ curves as $\varepsilon\rightarrow 0$ and extrapolate from our results the appearance of level–$\varepsilon$ curves for arbitrary small $\varepsilon$. This could be aided with improved numerical resolution. Such an extrapolation scheme is beyond the scope of this paper, however, and in the meantime we study the features of level–$\varepsilon$ curves for an appropriate value of $\varepsilon$.

We observe that for $\varepsilon=3\times 10^{-4},\;5\times 10^{-4},\;1\times 10^{-3}$, the level–$\varepsilon$ contours are very close to each other, showing that the level–$\varepsilon$ sets vary continuously with $\varepsilon$. We also observe that if $\varepsilon_{1}\leq\varepsilon_{2}$, then the 2D area enclosed by the level–$\varepsilon_{1}$ curve encloses the 2D area enclosed by the level–$\varepsilon_{2}$ curve. Thus, each panel of Figure 1 indicates that $|\mathcal{D}|$ decreases on average away from $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, and the plots of $|\mathcal{D}|$ show no global minima. However, Figure 2 indicates that the plots of $|\mathcal{D}|$ do have have local minima which approximately coincide with the level–$0$ sets of $\mathcal{D}_{r}$ and with the level–$\varepsilon$ sets for $\varepsilon=3\times 10^{-4},\;1\times 10^{-3}$.

In order to investigate further the level–$0$ sets of $|\mathcal{D}|$ (or equivalently the level–$0$ sets of $\mathcal{D}$), we find out where $|\mathcal{D}|$ takes a local minimum value. If the value of $|\mathcal{D}|$ itself is small, then these locations of the local minima could possibly indicate positions of the actual zeros of $|\mathcal{D}|$, which would be caused by numerical errors. It could also be the case that the GHC should be modified so that the GH is defined as the set of points where $|\mathcal{D}|$ reaches the local minimum instead of being identically zero. If this were the case, then locating the local minima of $|\mathcal{D}|$ would locate the GH precisely instead of approximating it. However, further evidence from the analysis of $\mathcal{D}_{r}$ below perhaps suggests that this is not the case.

To this end, in Figure 2, we examine 1D plots of $|\mathcal{D}|$ as functions of $y$ for fixed $x=x_{0}$. At time $t=12$, (resp. $t=\;16,\;20$), we display the $|\mathcal{D}|$ vs $y$ plot in the upper-right panel (resp. middle-right panel, lower-right panel) of Figure 2 for $x_{0}=0.000$ (resp. $x=0.0625,\;0.125$). Along each $|\mathcal{D}|$ vs $y$ plot, we find and track the values of $y=y_{min}$, where $|\mathcal{D}|$ assumes a local minimum value and lies in the range $\left[1\times 10^{-4},1.2\times 10^{-3}\right]$. The resulting coordinates $(x_{0},y_{min})$ are then superimposed on the level–$1\times 10^{-3}$ and level–$3\times 10^{-4}$ sets of $|\mathcal{D}|$ at $t=12$ (resp. $t=16,\;t=20$) in the upper-left panel (resp. middle-left, lower-left) panel of Figure 2.

More specifically, in the upper-right panel of Figure 2, where $t=12$ and $x_{0}=0.000$, we see that local minima of $|\mathcal{D}|$ occur roughly at $y_{min}=-0.8,\;0,\;0.8$. The points $(x_{0},y_{min})=(0,-0.8),(0,0),(0,+0.8)$ are then marked with green dots on the level–$\varepsilon$ plots of $\mathcal{D}$ at $t=12$, as shown on the upper-left panel of Figure 2. The remaining green dots on this upper-left panel are found similarly, using the positions $(x_{0},y_{min})$ of the local minima of $|\mathcal{D}|$, with $|\mathcal{D}|\in\left[1\times 10^{-4},1.2\times 10^{-3}\right]$, but with varying $x_{0}$.

One can similarly inspect the plot in the middle-right panel of Figure 2, where $t=16$ and $x_{0}=0.0625$, to find that $|\mathcal{D}|$ is minimized roughly where $y_{min}=-0.2,\;0.1$. The corresponding points $(x_{0},y_{min})=(0.125,-0.2),(0.125,0.1)$ are then recorded with green dots on the level–$\varepsilon$ plots of $\mathcal{D}$ at $t=16$, in the middle-left panel of Figure 2. The remaining green dots on this panel are again obtained by varying $x_{0}$ and mark the positions $(x_{0},y_{min})$ of the local minima of $|\mathcal{D}|\in\left[1\times 10^{-4},1.2\times 10^{-3}\right]$.

Finally, in the lower-right panel of Figure 2, where $t=20$ and $x_{0}=0.125$, the quantity $|\mathcal{D}|$ is minimized at roughly $y_{min}=0.2$, so that the point $(x_{0},y_{min})=(0.125,0.2)$ is marked with a green dot on the lower-left panel of Figure 2 and the remaining green dots track the positions $(x_{0},y_{min})$ of the local min of $|\mathcal{D}|$, $|\mathcal{D}|\in\left[1\times 10^{-4},1.2\times 10^{-3}\right]$, as above.

By studying the left three panels of Figure 2, we observe that at times $t=12$ and $t=16$, these local minima appear to track the green level–$3\times 10^{-4}$ sets, while at time $t=20$, these local minima appear to track more closely the white level–$1\times 10^{-3}$ sets. This shows that the positions of the local minima of $|\mathcal{D}|$ align closely with the level-$\varepsilon$ sets of $|\mathcal{D}|$ for $\varepsilon=3\times 10^{-4},\;1\times 10^{-3}$. In [49], the positions of $(x_{0},y_{min})$ are compared with the positions of the minima $(x_{min},y_{0})$–obtained through the above procedure but with the roles of $x$ and $y$ reversed–and it is possible that $(x_{min},y_{0})$ are indeed the overall local minima of $|\mathcal{D}|$. In this case, Figure 2 provides supporting evidence that the local minima of $|\mathcal{D}|$ and hence the level–$\varepsilon$ sets of $|\mathcal{D}|$, accurately track the GH. It remains as future work, however, to study the local minima more closely and devise and implement an algorithm to track its $(x,y)$ position through a BBH merger.

When studying the zero set of $\mathcal{D}$, it is also helpful to study separately $\mathcal{D}_{r}=\text{Re}(\mathcal{D})$ and $\mathcal{D}_{i}=\text{Im}(\mathcal{D})$, as these quantities change sign through a zero. To wit, we plot in Figure 3 the contour plots of $\mathcal{D}_{r}$ (left panels) and $\mathcal{D}_{i}$ (right panels), with magnified resolution, along with their level–$-0.01$ sets in yellow and their level–$+0.01$ sets in lime green at times $t=12,\;16,\;20$ in the upper, middle and lower panels, respectively.

In each of the frames in Figure 3, the grey regions correspond to $-0.01<\mathcal{D}_{r}<0.01$ (resp. $-0.01<\mathcal{D}_{i}<0.01$), the black regions correspond to $\mathcal{D}_{r}\geq 1$ (resp. $\mathcal{D}_{i}\geq 1$), and the white regions correspond to $\mathcal{D}_{r}\leq-1$ (resp. $\mathcal{D}_{i}\leq-1$). Interpolating between the level–$+0.01$ sets and the level–$-0.01$ sets of $\mathcal{D}_{r}$ (resp. $\mathcal{D}_{i}$), we deduce that there must be a surface among the level–$\pm 0.01$ sets of $\mathcal{D}_{r}$ (resp. $\mathcal{D}_{i}$) across which $\mathcal{D}_{r}$ (resp. $\mathcal{D}_{i}$) changes sign. This surface is then the level–$0$ set of $\mathcal{D}_{r}$ (resp. $\mathcal{D}_{i}$).

Upon inspection of each frame in Figure 3, we see that the level–$\pm 0.01$ sets of $\mathcal{D}_{r}$ (resp. $\mathcal{D}_{i}$) occur in close proximity with, but are contained in, the interior of the level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$. Therefore, Figure 3 provides strong evidence that the level–$1\times 10^{-3}$ set of $|\mathcal{D}|$ well approximates the level–$0$ set of $\mathcal{D}_{r}$ and $\mathcal{D}_{i}$ and, hence, the level–$0$ set of $\mathcal{D}$.

We next explicitly compare the level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$ with the corresponding MOTSs $\mathcal{S}_{1,2,i,c}$ in Figure 4. We display the 2D contour plots of $|\mathcal{D}|$ with magnified resolution at times $t=12$ and $t=16$ in the upper left and upper right panels, respectively, and we display the 2D contour plots of $|\mathcal{D}|$ at $t=20$ and in the lower left and lower right panels. As in Figures 1–3, the white curves denote the white level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$ and the blue curves mark the $(x,y)$ coordinates of points on $\mathcal{S}_{1}\cup\mathcal{S}_{2}$ whose corresponding $z$ coordinate values lie in the range $\left[0.02,0.04\right]$, as mentioned previously.

In the present quasi-circular simulation, the bifurcation of the the third MOTS into $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$ occurs between times $t=18.5$ and $t=18.75$. Once this happens, $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are no longer AHs, as the MOTS, $\mathcal{S}_{c}$, now surrounds $\mathcal{S}_{1}$, $\mathcal{S}_{2}$ and $\mathcal{S}_{i}$. Thus, in order to compare our level–$\varepsilon$ sets of $|\mathcal{D}|$ with AHs (the outermost MOTSs), we have included plots of $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$ here. The purple dots on the bottom left (resp. bottom right) panel of Figure 4 label the $(x,y)$ coordinates of the points on $\mathcal{S}_{i}$ (resp. $\mathcal{S}_{c}$) whose corresponding $z$ value lies in the range, $\left[-0.1,+0.1\right]$. Note that in the bottom right–hand corner, the outer MOTS at late times is so big that the entire two dots and the scale for scale of the white level–$0.001$ sets of $|\mathcal{D}|$ and $\mathcal{S}_{i,1,2}$ are squashed at the origin.

We find that the white level–$1\times 10^{-3}$ set of $|\mathcal{D}|$ coincides closely with $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, especially at early times. Since the MOTSs $\mathcal{S}_{c}$ and $\mathcal{S}_{i}$ have not yet formed, $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are AHs at early times. Hence, Figure 4 shows that the white level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$ well approximate the AH at early times. This lends support to the choice of the white level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$ as a representative approximation to the level–$0$ set of $|\mathcal{D}|$ (and hence of $\mathcal{D}$).

At later times, however, it appears that the AH, $\mathcal{S}_{c}$, diverges from the white level–$1\times 10^{-3}$ sets, so that the white level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$ no longer approximate the AH in this régime. Thus, it appears in this current simulation that the level-sets of the curvature invariant $\mathcal{D}$ does not detect the common outer horizon forming at late times. Furthermore, in this simulation, no significant patterns were observed in the level–$\varepsilon$ sets of $|\mathcal{D}|$ immediately prior to the formation of $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$ [49]. However, the authors intend to further study the GHs at the formation of $\mathcal{S}_{i}$ and $\mathcal{S}_{c}$, and also at later times, which would require rerunning these simulations with improved numerical resolution.

Therefore, Figures 1–4 provide strong evidence that one can define a unique smooth GH, theoretically given by the level–$0$ set of the complex invariant $\mathcal{D}=I^{3}-27J^{2}$, which we have found is best approximated in the numerics by the level–$1\times 10^{-3}$ sets of $|\mathcal{D}|$.