A binary black hole metric approximation from inspiral to merger

We start with the Kerr-Schild metric in Cartesian coordinates $\{X^{a}\}$ with arbitrary spin given by

with the null covector given by

the Boyer-Lindquist radius is defined as

the Cartesian radius is $R^{2}=X^{2}+Y^{2}+Z^{2}$, and the function $H$ is defined as

where $\epsilon^{i}_{jk}$ is the Levi-Civita symbol, $\delta_{ij}$ is the Kronecker delta, $a^{i}$ is the spin vector, and $a=\sqrt{a^{i}a^{j}\delta_{ij}}$ This is a compact form of writing a rotated Kerr-Schild metric [54]. We will use these coordinates to build a moving superposition of black holes.

In the Kerr-Schild form (1), the metric is represented by a superposition of a background flat metric, $\eta_{ab}$, and a black hole term, $2H(Xi)l^{a}l^{b}$. This provides a natural way of adding a second black hole to the spacetime. Because Einstein’s equations are highly non-linear, this superposition is, of course, not a vacuum solution; however, it naturally approaches an exact solution as the separation of the black holes increases. We will show in the next section that it is a good approximation for our purposes, i.e. investigating the behavior of matter around dynamic spacetime.

In the post-Newtonian sense, black holes in a binary are moving on accelerated orbits. To include the orbital motion of the BHs in the superposition, we build a time-dependent boost transformation, that can be considered as a generalization of the linear Lorentz boost. A time-dependent boost in general can be defined as a transformation from an inertial frame to an accelerated frame that is moving in a time-like trajectory. We now show how to build this transformation for an arbitrary orbit moving according to the post-Newtonian equations of motion.

Let us associate each moving black hole with an accelerated frame and the center of mass with the global inertial frame. We denote $\{X^{b}\}$ as the coordinates of the accelerated frame and $\{x^{a}\}$ as the coordinates of the inertial frame. In coordinates of the inertial frame, the trajectory of one of the BHs is given by $s^{a}(\tau)=\{t(\tau),s^{i}[t_{\tau}(\tau)]\}$ where $\tau$ is the proper time in the BH frame related to the global time as $\tau=\int^{t_{\tau}}dt\gamma^{-1}$. The trajectories we are using are parameterized in terms of the global time, so we will consider $s^{a}[t(\tau)]\equiv s^{a}(t_{\tau})$. We now construct the transformation between these coordinate systems for each BH, i.e., $X^{b}(x^{a})$, and use them to boost the BH terms, similar to an active Lorentz transformation.

The local coordinates of a moving frame in a general spacetime are the so-called normal or Fermi coordinates [55, 56, 57, 58], in which the metric is locally flat. We review how to build these coordinates, generalizing the methods in Refs. [59, 42]. In Fermi coordinates, the time component at the origin of the moving frame is defined as the proper time of the trajectory, $T=\tau$, while the space coordinates of a given point $X^{i}$ are constructed by finding the unique space-like geodesics passing through that point and orthogonal to the four-velocity $u^{a}=ds^{a}/d\tau$ [58]. In a Minkowski background spacetime, the full transformation can be written as:

where $e_{i}=e^{a}_{i}\partial_{a}$ is an orthonormal vector basis, Fermi-Walker transported by the accelerated frame and adapted to the local coordinates, i.e., the spatial axis of the frame [58]. Because the tetrad field lives on the tangential space of a moving frame, for each time $t$ we can find the components of the basis in global coordinates $\{x^{a}\}$ by performing a Lorentz boost. This assumes that, for a fixed time, there is an inertial frame that is equivalent instantaneously to an accelerated frame [59].

At $X^{i}=0$, the orthonormal frame in Fermi normal coordinates is simply $e^{b}_{i}=\delta^{b}_{i}$ and thus applying a local boost, we have, in the global coordinate system, $e^{a}_{i}=(\Lambda^{-1})^{a}_{i}$, where $(\Lambda^{-1})^{a}_{i}$ is the inverse of the Lorentz boost matrix, given by:

where $\xi=\tan^{-1}{(\beta)}$ is the rapidity, $n^{i}=\beta^{i}/\beta=(n_{x},n_{y},n_{z})$ is the unit direction of the velocity, and ${K}^{i}$ are the generators of the Lorentz group. Notice that $\Lambda^{-1}({\beta}^{i})=\Lambda(-\beta^{i})$.

Because for binary black holes $\gamma(t)$ changes slowly in time for most of the evolution, we approximate $t_{\tau}=\gamma\tau\equiv\gamma T$; The transformation from the inertial frame to the accelerated frame is then:

The transformation reduces to a standard linear boost when the acceleration is zero; indeed, for a frame moving at uniform velocity, we have $s^{i}(t_{\tau})=\beta n^{i}t_{\tau}=\beta n^{i}\gamma T$, and the Jacobian of the transformation in Eqs. (7) is simply given by $\Lambda^{a}_{b}$. Let us note two issues with this time-dependent boost. First, the transformation to (and from) the accelerated frame is valid only in a neighborhood of the worldline. The transformation diverges when $|X|\rightarrow c^{2}/\beta$; for a linearly accelerating frame, this is known as the Rindler horizon [60]. Beyond this horizon, the space-like geodesics used to construct the coordinates become time-like. A second issue is that to obtain $X^{b}(x^{a})$ from $x^{a}(X^{b})$ we need to invert Eqs. (7), which is a set of non-linear equations for arbitrary trajectories because $s^{i}(\tau)$ is, in general, non-linear unless the velocity is constant.

To avoid a non-linear inversion, we can expand the trajectory around the global time as $s^{i}(t_{\tau})=s^{i}(t)+\beta^{i}(t)(t_{\tau}-t)$. As we show in Appendix A, this approximation is good for a binary black hole and allows us to write the inverse as

where we define

To avoid coordinate pathologies at large distances, we discard the acceleration term when calculating the Jacobian of Eqs. (8). In this way, we obtain that the transformation is simply given by:

i.e. a simple Lorentz boost as defined before (see Appendix A for more details). The transformation is, in this way, well-defined at infinity.

The full binary black hole metric, which we name Superposed Kerr-Schild (SKS), is then given by the boosted superposition:

where subscript $A=1,2$ refers to each BH; the null covectors, $l_{c}\,(X^{d})$, and $H(X^{d})$ function in the frame of the BH are given by Eqs. (3) and (4), the coordinate transformation to the inertial frame is given by Eq. (8) and the Lorentz boost matrix $\Lambda^{d}_{a\,(A)}$ is given by Eq. (6).

At large separations, the SKS metric tends naturally to an exact solution of Einstein’s equation. We notice that even when the velocities are small (e.g. at very large separations), we still need to transform the superposition terms with the Lorentz boosts at first order in $c\beta$, i.e. we need to do the proper Galilean transformation of the metric. Otherwise, this solution does not have the right transformation properties in the Newtonian limit. In other words, even if $(1-\gamma)\sim 0$, there are still terms proportional to $c\beta$ in the transformation, so discarding them amounts to a first-order error in the solution.

The free parameters of the metric are given by the spin ${a}^{i}_{(A)}$, the mass, $M_{(A)}$, the spatial trajectories, ${s}^{i}_{(A)}$, and three-velocities, ${v}^{i}_{(A)}$, for each black hole.

Although the metric represents a relativistic, strong-field spacetime, the motion of the black holes can be computed using a post-Newtonian approximation where the black holes are treated as point masses on a flat background and Einstein’s equations are expanded in terms of $\epsilon\approx v^{2}\approx m/r$ to compute corrections to Newtonian theory. Because we are interested in modeling the late stages of the inspiral where velocities reach tens of percent of the speed of light, we use a high-order expansion up to 4PN order.

Given the position of each black hole in the (perturbed) flat background spacetime, ${s}^{i}_{(A)}$, we can eliminate the center of mass in the orbital equation, reducing the problem to a one-body equation of motion for ${x}^{i}={s}^{i}_{(1)}-{s}^{i}_{(2)}$. This is given in Harmonic coordinates by:

where ${a}_{\rm N}$, ${a}_{\rm 1PN}$, ${a}_{\rm 2PN}$, ${a}_{\rm 3PN}$, and ${a}_{\rm 4PN}$, are the contributions to the acceleration due to different PN orders; ${a}_{\rm SO}$, ${a}_{\rm PNSO}$, and ${a}_{\rm 3.5PNSO}$, are the spin-orbit coupling contribution and its 1PN, and 3.5PN order corrections respectively; ${a}_{\rm SS}$ is the spin-spin coupling term; and finally ${a}^{\rm BT}_{\rm RR}$, ${a}_{\rm RR1PN}$, ${a}_{\rm RRSO}$, ${a}_{\rm RRSS}$ are the radiation-reaction contribution with its 1PN correction, and contributions from spin-orbit and spin-spin interactions, respectively. Explicit expressions for all these terms can be found in Ref. [61], in the Appendix of Ref. [62], and explicitly in the CBWaves source code with corresponding references ¹¹1cbwave.

The spin-orbit coupling term, which can induce orbital precession in the system, is not unique and depends on additional spin supplementary conditions. We fix these with the covariant expression $S^{ab}u_{a}=0$, where $u^{a}$ is the four-velocity of each black hole and $S^{ab}$ is the spin tensor [64]. The spin vector is defined as $S^{i}=(1/2)\epsilon^{i}_{jk}S^{jk}$, related to the spin parameter of Kerr-Schild as $S^{i}=Ma^{i}$ and the dimensionless spin $S^{i}=\chi^{i}M^{2}$. Notice that the spin vector components are defined in the coordinate system associated with the moving BH.

Spins can precess due to spin-spin and spin-orbit coupling as the binary evolves. The equations of motions for the spins of one of the BHs are given by:

where $\epsilon^{i}_{jk}$ is the Levi-Civita symbol, and the precession vector ${\Omega^{j}}$ depends on the angular momentum of the binary and the spins of both black holes. We use the expression from Equation (3) in Ref. [65] containing terms up to 3.5PN order (see also references therein for explicit formulae).

We solve the coupled equation of motions for ${x}^{i}$, ${S}^{i}_{(1)}$, and ${S}^{i}_{(2)}$ with a Runge-Kutta scheme of fourth-order using a version of the public code CBWaves, a modulated C code that solves the PN equations efficiently. Initial conditions are set by choosing the semi-major axis, $r_{12}$, spin directions, and eccentricity, defined as $e=(r_{\rm max}-r_{\rm min})/(r_{\rm max}+r_{\rm min})$, where $r_{\rm max}$ and $r_{\rm min}$ are the apocenter and pericenter of the orbit, respectively. Once we obtain ${x}^{i}(t)$, we compute the individual trajectories ${s}^{i}_{(A)}$ and then the velocities ${v}^{i}_{(A)}$. Because our metric construction was formulated in Kerr-Schild coordinates, we must in principle transform the trajectories from Harmonic to Kerr-Schild using, e.g., the equations in Appendix A of Ref. [54]. However, we have verified that this transformation usually involves a small offset for the trajectories and barely changes the constraints.

When the binary approaches merger, the dynamics of the spacetime near the black hole horizon become highly non-linear; accurate evolution during merger and gravitational wave emission can only be predicted by solving Einstein’s equation numerically. Nevertheless, when separations are small $r_{12}\lesssim 5\,M$, the characteristic timescales from inspiral to merger [6], $t_{\rm insp}\sim[r_{12}/(2M)]^{4}\,M$, are comparable or smaller than the dynamical timescale of the gas surrounding the binary, $t_{\rm dyn}\sim r/v_{\rm dyn}\sim(r/M)^{3/2}\,M$, with $r\gtrsim r_{12}$. Moreover, we know that shortly after the merger and ringdown phase, the remnant will settle into a Kerr black hole, with properties that we can estimate well following analytical arguments [66], together with fitting formulae from NR [67]. This motivates us to approximate the merger process interpolating the superposed binary metric to a single black hole metric. We argue that this is a reasonable approximation for the binary spacetime if we only care about the behavior of matter surrounding the merger.

We interpolate between inspiral and post-merger through the mass and spin of the black holes in the following way. During binary evolution, the mass, spin, and velocity of each black hole are defined as:

where $M_{0\,(A)}$, ${a}^{i}_{0\,(A)}$, ${v}^{i}_{0\,(A)}$ are the mass, spin, and velocity during inspiral and $M_{\rm f}$, ${a}^{i}_{\rm f}=\chi^{i}_{\rm f}M_{\rm f}$, ${v}^{i}_{\rm f}$ are the corresponding properties of the remnant black hole after merger.

We use an interpolating function that is $0$ for $t<t_{\rm merger}-dt_{\rm buffer}$, $1$ for $t>t_{\rm merger}+dt_{\rm buffer}$ and smoothly interpolates from 0 to 1 as

where $t_{\rm merger}$ is the time of merger, defined as the time when the light rings of each black hole touch, and $dt_{\rm buffer}$ a transition buffer time to implement the interpolation. After $t>t_{\rm merger}+dt_{\rm buffer}$, the interpolation of BH masses and spin brings each superposed term to a single boosted Kerr black hole.

The remnant properties are computed using fitting formulas from numerical relativity. Following Ref. [68], we use Eq. (18) from Ref. [69] to compute the final remnant mass as

where $E_{\rm rad}$ is the radiated energy. To compute ${a}^{i}_{\rm f}$ we calculate the dimensionless spin $|{\chi}_{\rm f}|={a}_{\rm f}/M_{\rm f}$ and the direction of the spin with respect to the total angular momentum using Eqs. (12) and (19) from Ref.[70]. To compute the BH kick velocity ${v}^{i}_{\rm f}$ we use the compiled methods in Ref. [68], which in turn, makes use of Ref. [71] formulas.

The BH remnant properties depend on the mass-ratio $q=M_{(1)}/M_{(2)}$, the dimensionless spins $\chi^{i}_{(A)}$, the angle between spins and orbital angular momentum, $\cos{\theta_{(A)}}=\hat{S}^{(A)}\cdot\hat{L}$, and the azimuthal angle between the projections of spins onto the orbital plane, defined as $\cos{\Delta\Phi}=(\hat{S}_{(1)}\times\hat{L})/|\hat{S}_{(1)}\times\hat{L}|\cdot(\hat{S}_{(2)}\times\hat{L})/|\hat{S}_{(2)}\times\hat{L}|$. Here, ${L}$ is the orbital angular momentum and hats indicate unit vectors. These quantities are calculated at $t_{\rm merger}-dt_{\rm buffer}$.

The full expression of the metric given in Eq. (11) is written first in symbolic form using Mathematica notebook. This allows us to perform controlled checks and avoid typos in writing large expressions. The covariant metric, $g_{ab}$, as a function returning a 4x4 matrix, is then exported to C language in an optimized way as a function of coordinates and trajectories; in particular, we optimize syntactic mathematical operations before exporting to avoid multiple evaluations of expressions when possible. The black hole trajectories are calculated using a modified version of the public software CBwave, which is written in C and executed via a Python wrapper. The transition from inspiral to merger is calculated after this step in post-processing. The trajectories are then exported to an HDF5 table containing 1D arrays of each property as a function of time. We use linear interpolation to obtain the orbital parameters at any given time.

The implementation of the metric is now available in a public repository [72] which includes Mathematica notebooks, the PN solver (our version of CBWave), and explicit C functions to use it for numerical computations. The metric has been implemented, tested, and used in production in the public GRMHD codes Athena++ and within the framework of the EinsteinToolkit using either GRHydro and IllinoisGRMHD codes. These implementations will become publicly available soon.

For Run-BSSN we use puncture initial data to set an equal-mass binary at a coordinate separation of $10M$. This method solves the constraint equations using the so-called conformal-transverse-traceless decomposition [92, 93] where we choose a conformally flat metric $\gamma_{ij}=\psi^{4}\delta_{ij}$ with $\psi=1+2m_{(1)}/r_{1}+2m_{(2)}/r_{2}+u$; here, $m_{(A)}$ is the puncture parameter (or bare mass) of the A-th black hole and $u$ is a function that must be solved through the constraint equations. The extrinsic curvature is chosen to incorporate the initial momentum of the BH. We use $m_{(1)}=m_{(2)}=0.48595$ and an initial momentum in the $y$ direction of $P_{(1)\,y}=0.095433M$ and $P_{(2)\,y}=-0.095433M$ for each BH, which sets the binary on quasi-circular orbit and a total mass of $M\sim 1$.

For Run-SKS, we first solve the PN equations of two equal-mass black holes on a quasi-circular orbit (set via eccentricity reduction), choosing a coordinate separation of $10\>M$, and a total mass of $M=1$. We transition from inspiral to merger when the separation of the point masses is $r_{12}=3M$; we interpolate the spin, velocity, and mass of both black holes to the remnant properties given in this case by $M_{\rm f}=0.95$, zero kick velocity, and spin given by $\chi^{3}=0.66$, using a buffer time of $dt_{\rm buffer}=0.01\>M$ (see Section II.4).

For the fluid part, we set up a static uniform ideal gas of $\rho_{0}=10^{-2}$ (in code units) all over the domain assuming an adiabatic equation of state (EOS) with an adiabatic index $\Gamma=4/3$ to calculate the pressure and internal energy. We notice that the total density is decoupled from spacetime so the choice of $\rho_{0}$ is arbitrary. For this test, we set the magnetic field to a very low value but we effectively solve all MHD equations; we have also tested that the SKS metric works well with magnetic fields.

For Run-BSSN we evolve Einstein’s equation using the BSSN formulation and puncture gauge as explained in Section III. We neglect the energy-momentum tensor of matter assuming that the gas does not affect the spacetime

For Run-SKS, we implement the metric in the framework of EinsteinToolkit first decomposing the 4D metric (11) in 3+1 form as in Eq. (19). We select the natural slicing adapted to the coordinates given by $\alpha^{2}=-1/g^{tt}$, $\beta_{i}=g_{ti}$, and $\gamma_{ij}=g_{ij}$. We also need to compute the extrinsic curvature, $K_{ab}$, for finding the apparent horizon and its quasilocal properties; we compute derivatives using a finite-difference stencil of order two for the spatial and time derivatives.

On top of the spacetime, we evolve the finite-volume GRMHD equations in the Valencia formulation [94] using our version of the public code GRHydro. We solve the Riemann problem using an HLLE solver [95], fifth-order WENO-Z for reconstruction [96], and a scheme for transforming conservative to primitive variables following Ref. [97]. We use an ideal EOS with an adiabatic index of $\Gamma=4/3$.

The grid structure consists of 2 moving blocks with 7 (6) refinement levels each following the black holes, with a domain $-260\,M<x,y,z<260\,M$ for Run-BSSN (Run-SKS). The finer resolution for Run-BSSN is $dx=2/2^{6}\,M=M/32$ and for Run-BSSN is $dx=2/2^{5}\>M=M/16$. In both cases, the horizon is covered by 16 points, which is on the lower end of what is recommended for having a stable evolution. For Run-BSSN we evolve in time using a method of lines and a fourth-order Runge-Kutta with a Courant factor of $0.3$. For Run-SKS, because the spacetime is fixed and the convergence properties of the plasma are less strict than the spacetime evolution [98], we use a lower-order time integrator given by a second-order Runge-Kutta with a Courant factor of $0.3$. Since the velocity of the BHs is $v\sim 0.1-0.2c$, the metric changes on a cell crossing time given by $dt\sim 10dx/c$ where $dx$ is the size of the smallest cell in the simulation grid. The Courant condition for the plasma, however, requires a much shorter time-step of the order of $dt\sim 0.3dx/c$ (set by the Courant number times the cell crossing time for speeds $\sim$ $c$). In this way, we can safely update the metric every $\sim 10$ time-steps, saving additional computational time.

Using a prescribed metric offers several computational advantages compared to a full numerical evolution. Most importantly, evolving MHD equations on a semi-analytical metric background is much cheaper than performing a full NR simulation for several reasons. First of all, an evolution scheme for Einstein’s equations such as BSSN requires solving many extra variables and lengthy source terms while communicating information among possibly many CPU cores. The cost of our analytical metric comes from calculating all terms in Eq. (11), calculating the inverse and derivatives, and doing a binary search in a table to interpolate and obtain the position and velocity of BHs, both being cheap and local operations. Furthermore, a substantial boost in efficiency can be obtained by updating the metric every $N$ iteration and using a lower-order Runge-Kutta integrator. Finally, the resolution requirements of the SKS metric are less restrictive than for puncture coordinates since the horizon can be resolved with fewer cells, which allows us to use less resolution.

Using all these advantages, when we run our simulation on a total number of $567$ CPU cores we obtain that Run-BSSN evolves with an update rate of $15\,M/$hr while Run-SKS has an update rate of $77\,M/$hr. We remark again that these runs have the same effective resolution on the BH horizon but Run-BSSN uses 7 levels of refinement while Run-SKS uses 6. Using the same grid and time evolution, the performance boost is closer to $\sim 2.5$.

Because the gauge ($\alpha$ and $\beta^{i}$) evolve in time in a full numerical relativity evolution, it is hard to make a one-to-one comparison with our SKS metric (see Ref. [99]). Moreover, the puncture gauge used in the BSSN evolution is initially different from the Kerr-Schild superposition, even at large separations. For a single non-spinning BH in puncture coordinates, the radial coordinate of the horizon is $r\sim 0.75M$ [100]. In Kerr-Schild coordinates, the radial coordinate of the horizon is instead $R=2M$. While both coordinate systems agree far away from the BH, the horizon can be more easily resolved in Kerr-Schild coordinates; this is especially useful for avoiding numerical artifacts and the necessity of using more resolution. Nevertheless, this coordinate difference, noticeable for instance in the lapse in Fig. 1 (plotted for both simulations at $t\sim 600M$), could be problematic for comparing diagnostics. As we will show, these differences amount to a fixed constant offset in some quantities.

To show that the SKS metric is a good approximation for accretion physics even at very small separations, we evolve an initially static uniform gas onto the spacetime. This resembles a Bondi flow problem, which has been investigated for BH binaries [106, 47, 107]. For this system, the Bondi radius is given by $R_{\rm b}=2M/c^{2}_{\rm s,\infty}=2M/(\Gamma P_{\infty}/\rho_{\infty})\approx 25\,M$.

After a short radially infalling phase, the flow reaches a quasi-steady state as can be seen in the accretion rate evolution onto the BHs in Fig. 8. Close to the black holes, the flow resembles a Bondi-Hoyle-Lyttleton fluid solution [108], with a shock front that propagates in the direction of motion of the black holes, creating spiral waves due to the orbital motion of the binary. In between the black holes, a tidal bridge is formed creating a stationary high-density region.

All these features are remarkably similar for Run-SKS and Run-BSSN even one orbit before the merger, as shown, e.g., in the equatorial density in Fig. 7. We find that Run-SKS has higher density regions near the horizon, by a factor of a few, while the flow velocity in Run-BSSN is slightly higher (see Fig.8). The accretion rate measured at the event horizon $\dot{M}=\int\sqrt{-g}d\Omega u^{r}\rho$, where $u^{r}$ is the radial four-velocity in the BH frame, however, follows the same behavior in both simulations: it rises in the first $\sim 100M$, reaching a steady state that lasts until merger, when it settles into a new steady state very quickly. We observe only tens of percent differences between the accretion rates in each simulation. The specific energy of the flow averaged over the horizon differs by a factor of two, although this quantity is strongly gauge-dependent, especially near the horizon. The average velocity of the flow onto the horizon follows the same behavior in each run but differs throughout the beginning by a fixed constant offset that could be attributed to gauge differences in how we measure the velocity.