Nearly extremal apparent horizons in simulations of merging black holes

In the decade following Pretorious's breakthrough simulation [1], and the development of the moving puncture technique [2, 3], several research groups have made great strides toward simulating merging binary black holes (BBHs) with a variety of mass ratios and spins (for recent reviews, see, e.g., [4–7]). These BBH simulations reveal not only the emitted gravitational waves but also the behavior of the strongly warped, highly dynamical spacetime near the holes' horizons. For example, a number of research groups have explored the properties (mass, spin, and recoil velocity) of the final, remnant hole in BBH mergers [8–19]. These studies typically follow a 'scattering matrix' approach, understanding the nonlinear dynamics by exploring how different initial BBH configurations produce different remnant Kerr black holes. Some studies have explored the highly nonlinear dynamics of the strongly curved spacetime during the merger itself; these include recent work using tendex and vortex lines (analogous to electric and magnetic field lines) to visualize the curvature of simulated spacetimes [20] and recent work exploring how the remnant properties are imprinted on the emitted gravitational waves [21].

In this paper, we consider the extremality of dynamical black holes in numerical spacetimes. A single Kerr black hole with spin magnitude S, horizon area A, and mass M obeys the inequalities

$$\begin{eqnarray}&&8\pi S\leqslant A,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&S\leqslant {{M}^{2}},\end{eqnarray} \tag{ 2 }$$

or

$$\begin{eqnarray}&&\zeta \equiv \frac{8\pi S}{A}\leqslant 1,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\chi \equiv \frac{S}{{{M}^{2}}}\leqslant 1.\end{eqnarray} \tag{ 4 }$$

For a Kerr black hole, both χ and ζ can be interpreted as measures of the hole's extremality; a Kerr hole is nearly extremal if $\chi \approx 1$ and $\zeta \approx 1$.

The actual spins and extremalities of astrophysical black holes are uncertain, but there is observational evidence that nearly extremal black holes could exist. For instance, recent measurements using both continuum fitting and x-ray reflection fitting suggest that Cygnus X-1 (the first black hole discovered) is nearly extremal [22–24]; there are also measurements of nearly extremal spin in other stellar-mass black holes (such as in x-ray binaries GRS 1915+105 [25] and GX 339-4 [26]) and in supermassive black holes (e.g. Swift J0501.9-3239 [27]). (For reviews of black hole spin measurements, see, e.g., [28, 29].) Therefore, BBHs targeted by gravitational-wave detectors could contain nearly extremal black holes; this has motivated previous and ongoing efforts to simulate BBH mergers with nearly extremal spins [10, 30–37].

Can such numerical simulations of merging, nearly extremal black holes contain superextremal horizons? Answering this question requires a generalization of extremality from the Kerr solution, preferably one that can be measured quasilocally (e.g. on apparent horizons). For exact axisymmetry, an unambiguous spin S (the Komar spin angular momentum [38]) can be defined on the apparent horizon; for this spin measure, inequality (1) has been proven to hold for spacetimes satisfying the Einstein equations with non-negative cosmological constant and with non-electromagnetic matter fields obeying the dominant energy condition [39]. (See [40] and the references therein for a review of geometric inequalities that have been proven for axisymmetric black holes.)

BBHs, in contrast, do not have exact symmetries to facilitate definitions of conserved mass and angular momentum. In [41], Booth and Fairhurst argue that extremality bounds such as inequalities (1) or (2) are not well posed, since there are no symmetries to define mass and angular momentum. Booth and Fairhurst propose an alternative extremality parameter e that does not require approximate axial symmetry. In practice, numerical relativists typically measure black hole spin angular momentum χ (and thus the extremality ζ) in terms of approximate symmetries characterized by approximate Killing vectors.

In this paper, we explore two alternative measures of extremality in numerical simulations of merging black holes. First, we explore whether inequality (1) is violated on apparent horizons of area A, evaluating the spin S of the apparent horizon using approximate Killing vectors, which identify the best approximate symmetries available. We confine our attention to a set of BBH simulations performed using the Spectral Einstein Code (SpEC) [42]. Each simulation has equal masses and equal, aligned spins. Even with spins as high as $S/{{M}^{2}}=0.994$, the highest yet simulated, we observe no violation⁸ of inequality (1). When a common horizon first forms, however, A and $8\pi S$ can be very close to each other: for a BBH with an initial spin of $\chi =0.994$, we find ζ as high as 0.997.

Second, we revisit the Booth–Fairhurst extremality e, whose value depends on a choice of scaling for the ingoing and outgoing null normals to the horizon. We define a gauge-invariant lower bound of the extremality, e₀, as the smallest value e can take for any scaling of the null vectors. Whatever gauge-fixing prescription might be used, the extremality e will always be larger than or equal to e₀. The lower bound e₀ can be moderately large: for merging holes with initial spins $\chi =0.994$, e₀ of the common apparent horizon, just after it appears, is comparable to e₀ of a Kerr black hole with a spin of $\chi \approx 0.97$. This suggests that these newly formed common apparent horizons must be at least moderately close to extremality.

Third, we construct binary black hole initial data with superextremal marginally outer trapped surfaces that violate equation (1). As in [30], a larger marginally outer trapped surface (the apparent horizon) which satisfies (1) always encloses these overspun surfaces. Our lower bound of extremality e₀ exceeds unity on the inner, overspun surfaces, verifying that they are in fact superextremal. On the apparent horizons, we find that e₀ is less than unity but can exceed the value of e₀ for an extremal Kerr black hole.

The rest of this paper is organized as follows. In section 2, we summarize our numerical methods and our methods for measuring extremality, including our lower bound for the extremality of an apparent horizon. In section 3, we present results for the extremalities of apparent horizons in BBH simulations with merging, rapidly spinning black holes. We conclude in section 4.

The standard method for computing spin in numerical relativity is the following integral, carried out on the apparent-horizon 2-surface [44–46]:

$$\begin{eqnarray}&&S=\frac{1}{8\pi }\oint {{\omega }_{B}}{{\phi }^{B}}{\rm d}A.\end{eqnarray} \tag{ 5 }$$

Here capital latin indices index the tangent bundle to the 2-surface within spacetime. The one-form ${{\omega }_{B}}$ physically represents a surface angular momentum density and mathematically represents a connection on the normal bundle, defined as

$$\begin{eqnarray}&&{{\omega }_{A}}=e_{A}^{\mu }{{n}_{\rho }}{{\nabla }_{\mu }}{{\ell }^{\rho }},\end{eqnarray} \tag{ 6 }$$

where $\vec{n}$ and $\vec{\ell }$ are (respectively) the ingoing and outgoing null normals to the 2-surface in spacetime, arrows denote spacetime 4-vectors, $\vec{\nabla }$ is the spacetime covariant derivative, and $e_{A}^{\mu }$ is the spacetime representation of a basis for the tangent bundle (a projector to the 2-surface). In practice, the form ${{\omega }_{A}}$ is usually computed from the extrinsic curvature of the spatial slice (e.g. equation (A1) of [30]).

The vector field $\vec{\phi }$ on the 2-surface is a generator of rotations, and this vector field encodes the directional nature of the angular momentum. Here, we choose this vector field as in previous papers, using methods of approximate Killing vectors [30, 47–49]. An arbitrary 2-surface will generally not have symmetries, but this method finds the vector field that comes closest to a symmetry on an arbitrary 2-surface, in the sense of minimizing a residual of Killing's equation with respect to variations in the space of smooth vector fields.

First, because Killing's equation implies that $\vec{\phi }$ is divergence-free, we start with the condition

$$\begin{eqnarray}&&{{\phi }^{A}}={{\epsilon }^{AB}}{{D}_{B}}z,\end{eqnarray} \tag{ 7 }$$

for some smooth function z on the 2-surface, where D and are, respectively, the covariant derivative and Levi-Civita tensor intrinsic to the 2-surface. This condition renders the computed spin invariant under the boost-gauge ambiguity, which is an ambiguity of the scaling of the null normals. The null normals $\vec{\ell }$ and $\vec{n}$ are normalized only relative to one another, through the standard condition

$$\begin{eqnarray}&&\vec{\ell }\cdot \vec{n}=-1.\end{eqnarray} \tag{ 8 }$$

This condition is preserved by the rescaling

$$\begin{eqnarray}\begin{array}{rcl} {\vec{\ell }} & \mapsto & {\rm exp} (a)\vec{\ell }, \\ {\vec{n}} & \mapsto & {\rm exp} (-a)\vec{n}, \\ \end{array}\end{eqnarray} \tag{ 9 }$$

where a is an arbitrary function on the 2-surface, which can be interpreted as a rapidity. The surface angular momentum density transforms as

$$\begin{eqnarray}&&{{\omega }_{A}}\mapsto {{\omega }_{A}}+{{D}_{A}}a.\end{eqnarray} \tag{ 10 }$$

Thus, under a rescaling of the null normals, the quasilocal angular momentum transforms as:

$$\begin{eqnarray}&&S\mapsto S+\oint \left( {{D}_{B}}a \right)\left( {{\epsilon }^{BC}}{{D}_{C}}z \right){\rm d}A.\end{eqnarray} \tag{ 11 }$$

An integration by parts and the condition of zero torsion then show that S is unchanged by the boost transformation (9).

A minimization problem for the integral of the square of the shear of $\vec{\phi }$ (the remainder of Killing's equation) then implies that z must satisfy a generalized eigenproblem

$$\begin{eqnarray}&&{{D}^{4}}z+\vec{D}\cdot \left( R\vec{D}z \right)=\lambda {{D}^{2}}z,\end{eqnarray} \tag{ 12 }$$

where R is the intrinsic scalar curvature on the horizon and the eigenvalue λ is related to the overall shear of the vector field $\vec{\phi }$. We solve a spectral representation of this generalized eigenproblem, find the eigenfunction z with minimum eigenvalue⁹ , normalize it according to a prescription described in [30], and compute the quasilocal spin angular momentum as

$$\begin{eqnarray}&&S=\frac{1}{8\pi }\oint {{\omega }_{B}}{{\epsilon }^{BC}}{{D}_{C}}z{\rm d}A.\end{eqnarray} \tag{ 13 }$$

Through an integration by parts, this equation can equivalently be written as

$$\begin{eqnarray}&&S=\frac{1}{8\pi }\oint z\Omega {\rm d}A,\end{eqnarray} \tag{ 14 }$$

where $\Omega \equiv {{\epsilon }^{AB}}{{D}_{A}}{{\omega }_{B}}$ is a scalar curvature of the normal bundle of the 2-surface in spacetime. This quantity Ω is the same (up to a constant factor) as the imaginary part of the complex curvature defined in [50]. Also, apart from correction terms that happen to vanish on isolated horizons, it is also equal to the imaginary part of the Weyl scalar ${{\Psi }_{2}}$ in a tetrad adapted to the 2-surface, or equivalently, to the normal–normal component of the magnetic component of the Weyl tensor, referred to as the horizon vorticity in [20].

This approach to defining black hole spin springs from an assumption that the horizon is nearly axisymmetric, a condition that is blatantly violated in some of the cases we consider here. Thus, it is surprising and rather mysterious that we get any useful physics at all—such as the apparent satisfaction of the extremality bound $8\pi S\leqslant A$ (section 3). Moreover, while the above construction provides an apparently reasonable implementation of the intuitive idea of 'approximate symmetry,' it can be generalized rather significantly, providing alternative measures of spin that agree with this measure on axisymmetric surfaces, but which differ from it when the axisymmetry is broken. One such modified form is that in [48], and a broader family is presented in [51].

To make a clearer case for extremality bounds in highly-deformed horizons, we would like a measure of black hole extremality that is independent of any assumptions of symmetry and any particular definition of black hole spin. We address this in the next section.

On highly deformed horizons, standard measures of black hole spin, such as the one given in equation (14) and the surrounding discussion, have a questionable physical motivation. The root of the difficulty is the need to define the vector field $\vec{\phi }$ on the surface. In the absence of even approximate axisymmetry, there is no obviously preferred vector field to choose. We continue to define spin using approximate Killing vectors, even though their use is no longer geometrically well-motivated.

In [41], Booth and Fairhurst introduced a measure that avoids this issue. Their extremality parameter is

$$\begin{eqnarray}&&e=\frac{1}{4\pi }\oint {{\omega }_{B}}{{\omega }^{B}}{\rm d}A.\end{eqnarray} \tag{ 15 }$$

Intuitively, this might be understood as a squared 'quasilocal spin magnitude,' though as was noted in [41], because we are integrating over the continuous horizon 2-surface, it might better be understood as including information from all of the current multipole moments. We will specify this relationship more precisely below.

The practical benefit of working with e rather than S is clear: it removes the need for solving the eigenproblem given in equation (12), and it removes questions of interpreting the results when axisymmetry is strongly broken. However, e is not simply an ad hoc quantity chosen to avoid these practical issues. Reference [41] shows that the value of e is closely related to the question of whether fully-trapped surfaces exist within the marginally-trapped surface on which it might be computed, and that subextremal isolated horizons (subextremal in the sense of having positive surface gravity) satisfy $e\lt 1$.

In spite of these appealing features, there is one major practical drawback to working with the Booth–Fairhurst extremality, related to the issue of boost gauge described in section 2.1. Under the rescaling of $\vec{\ell }$ and $\vec{n}$ given in equation (9), e transforms as

$$\begin{eqnarray}&&e\mapsto e+2\oint {{\omega }^{B}}{{D}_{B}}a{\rm d}A+\oint |\vec{D}a{{|}^{2}}{\rm d}A.\end{eqnarray} \tag{ 16 }$$

To talk about the extremality of the horizon in terms of the quantity e, one must fix a boost gauge.

The suggestion made in [41] for fixing the boost gauge is a standard one in the dynamical horizon literature. Because the dynamical horizon is a spacelike object, it is natural to scale $\vec{\ell }$ and $\vec{n}$ such that the unit spacelike normal to the 2-surface, $\hat{s}=\left( \vec{\ell }-\vec{n} \right)/\sqrt{2}$, is tangent to the dynamical horizon. This condition has many attractive mathematical features [46], but it has one particularly worrisome drawback in the numerical context: as a horizon settles down to that of a Kerr black hole, it approaches a null surface in spacetime. Boosting the spatial vector $\hat{s}$ to this asymptotically null surface requires an arbitrarily large boost rapidity a. Reference [52] notes that the process of calculating certain quantities in the boost gauge adapted to the dynamical horizon becomes numerically ill-behaved as the horizon settles down to Kerr. Even more troubling, [52] finds inner horizons that are spacelike in some regions and timelike in others.

For these reasons, we have opted to take a different approach to fixing the boost gauge. Our approach begins with a decomposition of the angular momentum surface density into two scalar potentials ϖ and π,

$$\begin{eqnarray}&&{{\omega }_{A}}={{\epsilon }_{A}}^{B}{{D}_{B}}\varpi +{{D}_{A}}\pi .\end{eqnarray} \tag{ 17 }$$

Given any ${{\omega }_{A}}$, the two potentials can be computed (up to an irrelevant constant) by solving Poisson equations on the horizon 2-surface:

$$\begin{eqnarray}&&{{D}^{2}}\pi ={{D}^{A}}{{\omega }_{A}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{D}^{2}}\varpi ={{\epsilon }^{AB}}{{D}_{A}}{{\omega }_{B}}=\Omega .\end{eqnarray} \tag{ 19 }$$

Note that, by construction, the right-hand side of both of these equations averages to zero over the 2-surface, as is required for a solution to exist.

The two potentials ϖ and π are distinguished not only by parity considerations, but more importantly by how they behave under changes in boost gauge. Equation (9) implies that the sources of the above Poisson equations transform as

$$\begin{eqnarray}&&{{D}^{A}}{{\omega }_{A}}\mapsto {{D}^{A}}{{\omega }_{A}}+{{D}^{2}}a,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\Omega \mapsto \Omega ,\end{eqnarray} \tag{ 21 }$$

so that the potentials transform as

$$\begin{eqnarray}&&\pi \mapsto \pi -a,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\varpi \mapsto \varpi .\end{eqnarray} \tag{ 23 }$$

In other words, given any ${{\omega }_{A}}$, one can always find the transformation into a preferred family of boost gauges in which π is constant, and therefore

$$\begin{eqnarray}&&{{\omega }_{A}}={{\epsilon }_{A}}^{B}{{D}_{B}}\varpi ,\end{eqnarray} \tag{ 24 }$$

where $\varpi ={{D}^{-2}}\Omega$. This technique for fixing boost gauge is used, for example, in [53]. In this special family of boost gauges, the extremality is:

$$\begin{eqnarray}&&4\pi {{e}_{0}}=\oint \left( {{\epsilon }_{B}}^{C}{{D}_{C}}\varpi \right)\left( {{\epsilon }^{BE}}{{D}_{E}}\varpi \right){\rm d}A\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&=-\oint \varpi {{D}^{2}}\varpi {\rm d}A\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&=-\oint \varpi \Omega {\rm d}A\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&=-\oint \Omega {{D}^{-2}}\Omega {\rm d}A.\end{eqnarray} \tag{ 28 }$$

By fixing to the gauge where $\vec{D}\cdot \vec{\omega }=0$, we have derived an expression for extremality entirely computable in terms of boost-invariant data on the 2-surface. For this reason, the quantity above can be considered a boost-invariant version of the Booth–Fairhurst extremality. It should be noted, though, that in general this quantity is not the same as the extremality computed in a boost gauge adapted to a dynamical horizon.

Despite not being adapted to a dynamical horizon, the extremality defined here has some appealing properties. The one most relevant to our work here is that it provides a lower bound on the extremality, in the space of all choices of boost gauge. In our preferred family of boost gauges in which equation (24) holds, we have:

$$\begin{eqnarray}&&{{\omega }_{A}}={{\epsilon }_{A}}^{B}{{D}_{B}}\varpi ,\end{eqnarray} \tag{ 29 }$$

where, again, ϖ is a boost-invariant quantity. In any other boost gauge, ${{\omega }_{A}}$ must differ from this by a gradient:

$$\begin{eqnarray}&&{{\omega }_{A}}={{\epsilon }_{A}}^{B}{{D}_{B}}\varpi +{{D}_{A}}a.\end{eqnarray} \tag{ 30 }$$

The extremality computed in this other boost gauge is:

$$\begin{eqnarray}&&4\pi e=\oint \left( {{\epsilon }_{A}}^{B}{{D}_{B}}\varpi +{{D}_{A}}a \right)\left( {{\epsilon }^{BE}}{{D}_{E}}\varpi +{{D}^{B}}a \right){\rm d}A\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&=\;4\pi {{e}_{0}}+2\oint {{\epsilon }_{A}}^{B}{{D}_{B}}\varpi {{D}^{A}}a{\rm d}A+\oint |\vec{D}a{{|}^{2}}{\rm d}A.\end{eqnarray} \tag{ 32 }$$

The second term on the right-hand side vanishes after an integration by parts, and the final term is manifestly non-negative. Thus

$$\begin{eqnarray}&&e\geqslant {{e}_{0}},\end{eqnarray} \tag{ 33 }$$

with equality only in our preferred family of gauges. Thus, this measure of extremality can exceed 1 only if the Booth–Fairhurst extremality measure exceeds 1 in all boost gauges, including the one adapted to the dynamical horizon.

As noted above, this measure of extremality is closely related to the multipolar structure of the horizon. Closely related formalisms for defining such multipolar structure can be found in [54–56]. In particular, in [56], the eigenfunctions in equation (12) are taken to define the current multipoles of a horizon 2-surface. Owing to the self-adjointness of the operators on both sides of the generalized eigenproblem, the eigenfunctions satisfy the condition

$$\begin{eqnarray}&&\oint {{z}_{i}}{{D}^{2}}{{z}_{j}}{\rm d}A=0\quad {\rm if}\;{{\lambda }_{i}}\ne {{\lambda }_{j}}.\end{eqnarray} \tag{ 34 }$$

Hence the eigenfunctions can be normalized in such a way that

$$\begin{eqnarray}&&\oint {{z}_{i}}{{D}^{2}}{{z}_{j}}{\rm d}A=-{{\delta }_{ij}}.\end{eqnarray} \tag{ 35 }$$

The negative sign is necessary because D² is a negative-definite operator in the space of non-constant functions. Note also that this is not the normalization condition used when computing spin by approximate Killing vectors.

Using this orthogonality relation, the potential ϖ can be expanded as

$$\begin{eqnarray}&&\varpi =\mathop{\sum }\limits_{i}{{L}_{i}}{{z}_{i}},\end{eqnarray} \tag{ 36 }$$

where the L_i are current multipole moments, defined in [56] as

$$\begin{eqnarray}&&{{L}_{i}}=\oint \varpi {{D}^{2}}{{z}_{i}}{\rm d}A=\oint {{z}_{i}}{{D}^{2}}\varpi {\rm d}A=\oint {{z}_{i}}\Omega {\rm d}A.\end{eqnarray} \tag{ 37 }$$

Inserting this expansion into equation (26), we have

$$\begin{eqnarray}&&4\pi {{e}_{0}}=-\mathop{\sum }\limits_{i,j}{{L}_{i}}{{L}_{j}}\oint {{z}_{i}}{{D}^{2}}{{z}_{j}}{\rm d}A\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&=\mathop{\sum }\limits_{i,j}{{L}_{i}}{{L}_{j}}{{\delta }_{ij}}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&=\mathop{\sum }\limits_{i}L_{i}^{2}.\end{eqnarray} \tag{ 40 }$$

Finally, note that the relationship between dimensionless spin χ and our gauge-invariant extremality parameter e₀ can be calculated explicitly on a Kerr horizon. A straightforward calculation gives:

$$\begin{eqnarray}\begin{array}{rcl} e_{0}^{{\rm Kerr}}(\zeta ) & = & \frac{1}{8\zeta (1+{{\zeta }^{2}})}\left[ 3\zeta +8{{\zeta }^{3}}-3{{\zeta }^{5}} \right. \\ {} & {} & +\left. \left( -3+3{{\zeta }^{2}}+7{{\zeta }^{4}}+{{\zeta }^{6}} \right){{{\rm tan} }^{-1}}(\zeta ) \right], \\ \end{array}\end{eqnarray} \tag{ 41 }$$

where

$$\begin{eqnarray}&&\zeta \equiv \frac{8\pi S}{A}=\frac{\chi }{1\pm \sqrt{1-{{\chi }^{2}}}}=\frac{a}{M\pm \sqrt{{{M}^{2}}-{{a}^{2}}}},\end{eqnarray} \tag{ 42 }$$

where the + sign is chosen for subextremal holes. Figure 1 compares the extremality measures e₀ and ζ for Kerr horizons.

Zoom In Zoom Out Reset image size

Note one interesting fact about the lower bound e₀: while the $\chi =0$ case gives ${{e}_{0}}=0$, as one would expect, in the extremal case, $\chi =1$, $\zeta =1$, we have $e_{0}^{{\rm Kerr}}(\zeta =1)=(4+\pi )/8\approx 0.893$. Thus, in this measure, even extremal Kerr does not have an extremality of 1. This could be seen as a weakness of our choice of boost-gauge, and a reason to prefer the family of boost gauges favored in [41], in which extremal Kerr does indeed have an extremality of 1. On the other hand, in section 3.3, we will construct BBH initial data containing highly distorted, highly dynamical horizons on which ${{e}_{0}}\gt e_{0}^{{\rm Kerr}}(\zeta =1)$ while remaining subextremal by all of our measures (cf figure 7). Thus, one might intuitively conclude that some amount of dynamics is necessary for the extremality parameter e₀ to approach unity.

We performed the numerical simulations used in this paper with SpEC [42]. We construct [57] quasi-equilibrium [30, 58] initial data to solve the Einstein constraint equations [59] for binaries with low ($\sim {{10}^{-4}}$) eccentricity [60–62]. In particular, following [30], we base our initial data on a weighted superposition of two boosted, spinning Kerr–Schild black holes.

We evolve the initial data using a generalized harmonic formulation [63–66] of Einstein's equations and damped harmonic gauge [67–69]. The adaptively-refined [34, 70] grid extends from pure-outflow excision boundaries conforming to the shapes of the apparent horizons [69, 71–73] to an artificial outer boundary, where we enforce constraint-preserving boundary conditions [66, 74, 75]. After the holes merge, the grid has only one excision boundary [71, 72]. We use a pseudospectral fast-flow algorithm [76] to find apparent horizons.

Building on the methods of [16, 35], we evolve this initial data through inspiral, merger, and ringdown. Evolutions with very high black hole spins are particularly challenging; for instance, maintaining pure-outflow excision boundaries near the time of merger is especially challenging in this case. A companion paper [77] discusses our methods for handling nearly extremal black hole spins.

We consider a family of BBH simulations with equal masses and equal spins aligned with the orbital angular momentum (table 1). Results from most of these have been previously published [16, 35, 36]. The simulations $S_{0.99}^{++}$ and $S_{0.994}^{++}$ [77] are new, as are the simulations $S_{0.75}^{++}$ and $S_{0.96}^{++}$. We measure the approximate-Killing-vector spin S and area A on each apparent horizon. The mass M is given by the Christodoulou formula¹⁰

$$\begin{eqnarray}&&{{M}^{2}}\equiv M_{{\rm irr}}^{2}+\frac{{{S}^{2}}}{4M_{{\rm irr}}^{2}},\end{eqnarray} \tag{ 43 }$$

where ${{M}_{{\rm irr}}}\equiv \sqrt{A/16\pi }$ is the irreducible mass. For most of the simulations, we perform calculations at three different numerical resolutions; however, for practical reasons¹¹ we consider only the highest resolution for simulations $S_{0.85}^{++}$ and $S_{0.97}^{++}$. This does not significantly affect our results or conclusions.

Table 1.  The numerical binary black hole simulations examined in this paper. All simulations have equal masses and equal spins aligned with the orbital angular momentum. We name each simulation $S_{\chi }^{++}$, where χ is the spin of the holes after the initial relaxation. Each simulation available in the SXS catalog [78] has the label SXS:BBH:ID. (Note that the publicly released simulation SXS:BBH:0178 used a newer version of SpEC than the version used for S⁺⁺_0.994 in this paper. The differences between the two versions are on the order of the errors we quote here). The 'Ref.' column lists references describing those simulations that are presented elsewhere. The quantity ${{\zeta }_{{\rm relax}}}$ represents ζ and χ respectively, on the individual horizons after the initial relaxation; ${{\zeta }_{{\rm final}}}$ and χ_final show ζ and χ, respectively, on the common horizon at the final time. Uncertainties are estimated by comparing at three different numerical resolutions. There is no uncertainty listed for χ and ${{\zeta }_{{\rm relax}}}$ in case $S_{0.994}^{++}$, because in this case lower resolutions only differ from the high resolution after the first 3.5 orbits (see [77] for details).

Name	ID	Ref.	$\chi ={{\chi }_{{\rm relax}}}$	${{\zeta }_{{\rm relax}}}$	${{\zeta }_{{\rm final}}}$	χfinal
$S_{0.75}^{++}$	0175	—	0.74994(7)	0.45136(6)	0.61982(6)	0.89558(4)
$S_{0.8}^{++}$	0155	[16]	0.79987(1)	0.499868(5)	0.63910(2)	0.90753(1)
$S_{0.85}^{++}$	0153	[16]	0.84983	0.55651	0.659292	0.91909
$S_{0.90}^{++}$	0160	[16]	0.899737(3)	0.626370(5)	0.68047(5)	0.93021(2)
$S_{0.95}^{++}$	0157	[16]	0.949586(7)	0.722940(2)	0.70275(6)	0.94085(3)
$S_{0.96}^{++}$	0176	—	0.95956(9)	0.7488(2)	0.70733(4)	0.94291(2)
$S_{0.97}^{++}$	0158	[35]	0.96950	0.778672	0.712011	0.94496
$S_{0.98}^{++}$	0172	[36]	0.97941(6)	0.8149(3)	0.71666(3)	0.94696(1)
$S_{0.99}^{++}$	0177	[77]	0.9893(2)	0.86306(6)	0.72135(1)	0.948930(5)
$S_{0.994}^{++}$	(0178)	[77]	0.9942	0.89805	0.723761(5)	0.949924(2)

In this section, we evaluate the extremality measure $\zeta \equiv 8\pi S/A$ for the numerical simulations described in section 2.

In each simulation, after the initial data relax and emit spurious gravitational radiation, the extremality ζ remains nearly constant until near the time of merger. Just before merger, ζ of the individual horizons decreases as A increases and S decreases; this is caused by tidal heating. The common horizon initially has a very large ζ that quickly decreases as the hole relaxes to the Kerr geometry.

Figure 2 illustrates this for $S_{0.994}^{++}$, and figure 3 zooms in on times near the formation of the common apparent horizon. During the inspiral and late in the ringdown, when the holes are nearly Kerr, χ and ζ are consistent with equation (42). For our higher-spin simulations, including $S_{0.994}^{++}$, we find that ζ is closest to unity (but rapidly decreasing) on the common horizon just after merger. This motivates us to determine with high precision the time when the common horizon first appears and the value of ζ at that time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

During a simulation, we typically search for the common horizon at regular time intervals, starting at a time when the holes' separation becomes sufficiently small. This gives us a rough measure of when the common horizon forms. For the simulations described here, we reran the portion of each simulation near common horizon formation, searching for a common horizon at each time step of the simulation. In each simulation, we find a common horizon at some earliest time step t₁ but not at previous time steps.

We expect, by analogy with figure 5 of [79], that the common horizon should first appear as a single surface that immediately bifurcates into an inner and outer surface. We find that this expectation holds in the simulations we have examined; as a concrete example, we plot ζ as a function of time $(t-{{t}_{1}})/{{M}_{{\rm tot}}}$ for simulation $S_{0.994}^{++}$ (figure 4), where ${{M}_{{\rm tot}}}$ is the sum of the initial Chrsitodoulou masses of the individual apparent horizons. We compute ζ at the ten earliest time steps where we can find a common horizon, and we fit these ten points to a parabola. We find that the local minimum of this parabola lies at some time t₀ between the earliest time step t₁ where we can find the horizon and the latest step ${{t}_{-1}}$ where we cannot find it. This gives us confidence that the common horizon does appear at the local minimum of this parabola.

Zoom In Zoom Out Reset image size

Thus we conclude that the common apparent horizon has its largest value of ζ at time t₀, the moment when it first appears. For our simulations with nearly extremal spins, this is the largest value of ζ on any apparent horizon during the simulation. In practice, t₀ lies between simulation time steps; because ζ is changing so quickly as a function of time, accurate estimates of ζ at $t={{t}_{0}}$ (which we will refer to as ${{\zeta }_{Ci}}$) require extrapolation.

We compare four different estimates of t₀ and ${{\zeta }_{Ci}}$:

• (1)  
  cubic Hermite extrapolation of $\zeta (t)$ to find t₀ and ${{\zeta }_{Ci}}$,
• (2)  
  cubic Hermite extrapolation of the irreducible mass¹² ${{M}_{{\rm irr}}}(t)$ to find t₀, followed by cubic Hermite extrapolation of $\zeta (t)$ to find ζ at t₀,
• (3)  
  fit a parabola to the 10 earliest measurements of time t as a function of the common horizon's extremality ζ, and extrapolate using the fitted parabola to find ${{\zeta }_{Ci}}$, as in figure 4,
• (4)  
  estimate ${{\zeta }_{Ci}}$ as the value of ζ at time t₁, the earliest simulation time step where we find the common horizon.

In figure 5, we compare these four estimates of ${{\zeta }_{Ci}}$ as functions of the dimensionless spin $\chi \equiv S/{{M}^{2}}$ measured on one of the individual horizons after the initial relaxation. We find an approximately linear relationship. The first three estimates, each of which involve extrapolation, give consistent estimates of ${{\zeta }_{Ci}}$ that are typically larger than $\zeta (t={{t}_{1}})$.

Zoom In Zoom Out Reset image size

We conclude from figure 5 that the common horizon satisfies the inequality $\zeta \equiv 8\pi S/A\leqslant 1$ at all times, even when it first appears. However, for nearly extremal initial black hole spins, the inequality is almost violated at time t₀, when the common horizon first appears. The approximately linear dependence on the initial, relaxed spin χ suggests that this inequality could be slightly violated with even higher initial black hole spins than we have simulated. We speculate that this does not happen in practice, but that instead ${{\zeta }_{Ci}}$ depends on χ such that ${{\zeta }_{Ci}}$ never exceeds unity for any $\chi \lt 1$. Verifying this speculation would require additional simulations beyond SpEC's current capabilities.

In the previous section, we found that the inequality $\zeta \equiv 8\pi S/A\leqslant 1$ was satisfied for all of our simulations, but that it was nearly violated on the common apparent horizons when they first appeared. This demonstrates that our measure of the horizon spins, based on approximate Killing vectors, behaves at least somewhat sensibly, even when the horizons are not very axisymmetric. But does this necessarily mean that these early common horizons are very close to extremality? Or, might it be, because the common horizons are far from axisymmetry, that the approximate-Killing-vector spin measure S does not represent spin in a physically meaningful way?

To address this question, in this section we consider an independent measure of the extremality of the common horizons, the lower bound e₀ introduced in section 2.2. This quantity is gauge and slicing invariant and does not rely on approximate axisymmetry. However, it is only a lower bound, and even an extremal Kerr hole has $e_{0}^{{\rm Kerr}}(\zeta =1)\approx 0.89$.

Figure 6 compares $e_{0}^{{\rm Kerr}}(\zeta )$ (dashed curve) to e₀ measured on the common apparent horizons when they first appear (data points, shown for all resolutions), using cubic Hermite interpolation and extrapolation to determine ${{e}_{0,Ci}}$ (method 1 in the previous section). The newborn common apparent horizons have ${{e}_{0,Ci}}$ as large as $e_{0}^{{\rm Kerr}}(\zeta \approx 0.8)$. While this is significantly below $e_{0}^{{\rm Kerr}}(\zeta =1)$, it is still moderately extremal; note that a Kerr black hole with $\zeta \approx 0.8$ has a spin $\chi \approx 0.97$. This suggests that the common horizons are at least moderately close to extremality when they first appear, independent of whether the approximate-Killing-vector method is a good measure of spin on apparent horizons that are far from axisymmetric.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Is it ever possible for apparent horizons in numerical simulations of merging black holes to violate the inequality ${{e}_{0}}\leqslant 1$? We observe no such violations in the numerical simulations we have considered so far. In this section, we attempt to construct initial data for merging black holes with apparent horizons with ${{e}_{0}}\gt 1$.

We follow the methods of [30] and the references therein to construct superposed-Kerr–Schild initial data for two equal-mass black holes with spins of equal magnitude aligned with the orbital angular momentum. In this initial data method, the initial spatial metric g_ij is proportional to a conformal metric that is a weighted superposition of two boosted, spinning Kerr–Schild black holes $\tilde{g}_{ij}^{{\rm SKS}}$:

$$\begin{eqnarray}&&{{g}_{ij}}={{\psi }^{4}}\tilde{g}_{ij}^{{\rm SKS}}.\end{eqnarray} \tag{ 44 }$$

Two regions are excised from the computational domain, and boundary conditions are imposed on the excision surfaces. One boundary condition enforces that the excision surfaces have zero expansion, i.e., that they are marginally outer trapped surfaces. Another controls the holes' spins via a parameter ${{\Omega }_{r}}$, which adjusts the tangential part of the shift on the boundary.

Following [30], we construct a family of initial data sets with the same conformal metric $\tilde{g}_{ij}^{{\rm SKS}}$ whose Kerr–Schild black holes have dimensionless spins $\tilde{\chi }\equiv \tilde{S}/{{\tilde{M}}^{2}}=0.99$. Each initial data set has a different choice of ${{\Omega }_{r}}$. Figure 7 shows ζ and e₀ as a function of ${{\Omega }_{r}}$. As in [30], we find that for sufficiently large values of ${{\Omega }_{r}}$, we can construct initial data with zero-expansion 'inner horizons' that have $\zeta \gt 1$. These surfaces are enclosed by 'outer horizons' (the apparent horizons) satisfying $\zeta \lt 1$.

Figure 7 also shows e₀ for these initial data sets. We are able to construct initial data with inner horizons with extremality lower bounds ${{e}_{0}}\gt 1$. These surfaces are superextremal, but they are always enclosed by larger apparent horizons that are subextremal ($\zeta \lt 1$ and ${{e}_{0}}\lt 1$). However, note that some of the apparent horizons do have ${{e}_{0}}\gt e_{0}^{{\rm Kerr}}(\zeta =1)$.