# Superposed metric for spinning black hole binaries approaching merger

Luciano Combi<EMAIL_ADDRESS>Instituto Argentino de Radioastronomía
(IAR, CCT La Plata, CONICET/CIC), C.C.5, (1984) Villa Elisa, Buenos Aires,
Argentina  Center for Computational Relativity and Gravitation, Rochester
Institute of Technology, Rochester, NY 14623.
Federico G. Lopez Armengol  Manuela Campanelli  Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, Rochester, NY
14623.
Brennan Ireland  Center for Computational Relativity and Gravitation,
Rochester Institute of Technology, Rochester, NY 14623.
The US Agency for International Development - 1300 Pennsylvania ave NW,
Washington, DC, 20004
Scott C. Noble  Gravitational Astrophysics Laboratory, Goddard Space Flight
Center, Greenbelt, MD 20771.
Hiroyuki Nakano  Faculty of Law, Ryukoku University, Kyoto 612-8577, Japan
Dennis Bowen  Center for Theoretical Astrophysics, Los Alamos National
Laboratory, P.O. Box 1663, Los Alamos, NM 87545.
X Computational Physics, Los Alamos National Laboratory, P.O. Box 1663, Los
Alamos, NM 87545.

###### Abstract

We construct an approximate metric that represents the spacetime of spinning
binary black holes (BBH) approaching merger. We build the metric as an
analytical superposition of two Kerr metrics in harmonic coordinates, where we
transform each black hole term with time-dependent boosts describing an
inspiral trajectory. The velocities and trajectories of the boost are obtained
by solving the post-Newtonian (PN) equations of motion at 3.5 PN order. We
analyze the spacetime scalars of the new metric and we show that it is an
accurate approximation of Einstein’s field equations in vacuum for a BBH
system in the inspiral regime. Furthermore, to prove the effectiveness of our
approach, we test the metric in the context of a 3D general relativistic
magneto-hydrodynamical (GRMHD) simulation of accreting mini-disks around the
black holes. We compare our results with a previous well-tested spacetime
construction based on the asymptotic matching method. We conclude that our new
spacetime is well-suited for long-term GRMHD simulations of spinning binary
black holes on their way to the merger.

Usage

Secondary publications and information retrieval purposes.

Structure

You may use the description environment to structure your abstract; use the
optional argument of the `\item` command to give the category of each item.

††preprint: APS/123-QED

##  I  Introduction

There is abundant evidence that most galaxies harbor supermassive black holes
(SMBHs) at their centers [1, 2]. A non-negligible fraction of these galaxies
undergoes one or more mergers within a Hubble time [3, 4]. After two galaxies
merge, various processes, such as dynamical friction, might drive their SMBHs
to close separations of sub-parsec scales [5, 6, 7].

In this situation, the binary system starts emitting gravitational radiation
efficiently, losing energy, and eventually merging [8, 9, 10]. The frequencies
of these gravitational waves (GWs) span from nano-Hertz for the inspiral phase
[11], up to milli-Hertz for the merger. Pulsar timing array consortiums and
LISA are actively working towards detecting these GW for the first time in the
next decade [12, 13, 14, 15].

Because galaxy mergers can be very efficient at driving interstellar gas
toward the galactic center [16, 17, 3], SMBH binary mergers should accrete
enough gas and emit observable electromagnetic (EM) radiation [18]. The total
energy radiated is proportional to the gas mass present during the SMBH binary
merger [19], and its characteristic form of this EM emission is expected to be
distinct relative to ordinary accreting supermassive black hole holes, such as
a single active galactic nucleus (AGN) [20, 21, 22]. Identifying those
signatures, however, requires complex calculations of the SMBH binary mergers
and their associated luminosity, spectrum, and time-dependence. Some numerical
calculations have started to reveal interesting properties associated with
these EM signals [23, 24, 25, 26, 27]. However, because they depend strongly
on the system’s properties such as the total binary mass, mass-ratio, spins
(magnitude and direction), and accretion rate, much work remains to be done.

Since the interstellar gas of the merged galaxies is expected to have a
considerable amount of angular momentum, a circumbinary disk will form around
the BHs [28]. Semi-analytical models of these systems predicted that the
binary would decouple from the fluid and coalesce in a dry merger; in other
words, the inflow time of the accreted matter would become larger than the
inspiral time at short orbital separations [29, 30].

Since the equations of magneto-hydrodynamics (MHD) in dynamical spacetimes are
highly non-linear, we need numerical simulations to make accurate predictions.
In the past decade, $\alpha$-viscous simulations and 3D MHD simulations
demonstrated that accretion onto the binary occurs even in the late inspiral
phase [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 25,
27], with appreciable accretion sustained right up to the time of merger when
using relativistic inspiral rates [38]. These simulations also showed that the
circumbinary disk is truncated at a distance approximately twice the binary
separation from the system’s center-of-mass. Outside this truncation radius,
mass piles up, forming a local peak in the surface density profile; inside
this radius, the accretion flow onto the binary is confined within two narrow
streams traversing a low-density gap. Each of these streams is associated with
one of the BHs, forming mini-disks around each hole.

Moreover, most of these simulations revealed the formation of a characteristic
$m=1$ mode overdensity, or lump, in the circumbinary disk for mass-ratios
close to unity [31, 38, 39, 40, 41, 43, 44, 46, 27]. In that case, the lump
modulates the accretion of the system and feeds the BHs with a single-arm
stream. If the mass ratio is small, the lump is weak [47], and the lighter BH
receives most of the mass, carving a path near the inner edge of the
circumbinary disk.

Solving Einstein’s field equations for the metric of SMBH binaries coupled to
MHD fields is computationally challenging. This has been done in the past in
the force-free regime [48, 49, 50, 51], or for close binary separations [37,
36, 52, 53, 54]. Recently, authors in Ref. [55] have evolved a spinning BBH in
full GRMHD for a few orbits in the inspiral regime, focusing on mini-disks
dynamics.

In order to model realistic scenarios in these systems, however, we need first
to evolve the circumbinary disk for many binary orbits until reaching a
steady-state [56]. In particular, the presence of the $m=1$ lump mode in the
circumbinary disk is very important to determine the dynamics of what happens
in the inner cavity [57]. An alternative approach to achieve this is to use
approximate, semi-analytical, solutions of Einstein’s equations for spacetime,
and evolve the MHD equations on it. This allows one to explore the parameter
space of the spacetime more efficiently and to focus computational resources
(such as the configuration of the grid) on the MHD fields. As an example of
this approach, previous work evolved relativistic circumbinary disks with a
post-Newtonian (PN) metric during the inspiral regime [58, 59, 47]. Since this
metric is only valid far from the BHs, these simulations must excise the
binary region from the computational domain.

To analyze the strong-field behavior of the plasma near the BHs, we need a
background metric that is valid at these scales. Such metric can be built, for
instance, through the so-called asymptotic matching approach [60, 61, 62] that
stitches different known analytical approximations for a binary black hole
(BBH) metric. In this approach, a perturbed Schwarzschild or Kerr solution is
used for the inner-zone, a PN expansion for the near-zone, and a post-
Minkowskian expansion for the far-zone are glued together via the transition
techniques developed in Refs. [63, 64]. This metric has been used to perform
GRMHD simulations of accretion flows with mini-disks around non-spinning BHs
for the first time [65, 66, 57].

This approach for the spacetime construction can be generalized to spinning
BHs [67, 61, 68], but the analytical matching metric becomes too complex and
computationally expensive for long-term GRMHD simulations. This motivates the
search for more efficient approaches for building an analytical spinning BBH
metric. From numerical relativity simulations, we know that spins play a key
role in BBH inspirals and mergers. Spins aligned with the orbit, for example,
can significantly alter the pace of orbital evolution by gravitational
radiation [69, 70, 71]. Oblique spins can drive complex precession and
nutation whose amplitude increases rapidly at smaller separations [72]; spins
with partial orbital alignment can repeatedly flip sign [73, 74, 75, 76];
spin-orbit PN resonances can also tilt the orbital orientation [77].

In Ref. [56], we introduced a new approach for spinning BBHs, building the
approximate metric as a linear superposition of two boosted Kerr-Schild BHs.
We used this new approximate spacetime to analyze the accretion of a
circumbinary disk around the BBH in a Keplerian orbit and unveil the influence
of the spin in the circumbinary flow. We found that streams falling into the
binary cavity as well as the accretion rate are affected by the magnitude and
direction of the BH spins, while other properties in the bulk of the
circumbinary disk remain unaffected. In particular, due to frame dragging
effects, accretion decreases (increases) when the spins are (anti-)aligned,
with important effects on the overall luminosity of the system. We also find
that the circumbinary disk is stabilized after more than 100 orbits, which
implies that long simulations are required for making realistic predictions.
In the present work, we formalize and extend this previous approach by
superposing two boosted Kerr BHs in harmonic coordinates, solving the PN
equations of motion (EOM) for the BH trajectories. This allows us to have a
more accurate approximation for the spacetime and to analyze the influence of
the inspiral on matter orbiting the BBHs. We test this metric by analyzing its
spacetime scalars and using it in a full 3D GRMHD simulation of accreting
mini-disks. We compare our results with the more expensive and complex
matching metric for non-spinning BHs. Having tested the viability and accuracy
of our spacetime, in an upcoming work we will use our new metric for analyzing
the influence of spins in the mini-disk dynamics and their outflows.

We organize the paper as follows: in Section II, we build the approximate
metric in the harmonic gauge by boosting and superposing two BH terms. In
Section III, we test the metric by analyzing its spacetime scalars and
comparing them with the asymptotic matching approach. In Section IV, we test
the metric as a background spacetime for a GRMHD simulation, comparing again
with previous results from the asymptotic matching approach. We conclude that
the metric is accurate and robust to be used in accretion disk simulations of
BBH in the inspiral regime.

Notation and conventions. We use the signature $(-,+,+,+)$ and we follow the
Misner-Thorne-Wheeler convention for tensor signs. We use geometrized units,
$G=c=1$. We use Latin letters $a,b,c,...=0,1,2,3$ for four dimensional
components of tensors, and $i,j,k,...=1,2,3$ for space components. An
orthonormal space basis is written as $\vec{e}_{(i)}=e^{a}_{(i)}\partial_{a}$,
where its components are denoted as $(i),(j),(k),...=1,2,3$.

##  II Construction of superposed binary black holes metric

To solve Einstein’s field equations with numerical methods, one usually starts
from a three-dimensional slice of the spacetime metric and matter fields as
initial data for the problem [78, 79, 80]. This initial data cannot be
arbitrary since it has to satisfy the constraints of Einstein’s equations;
there is, however, a significant amount of freedom to choose it because the
equations are invariant under diffeomorphisms. In General Relativistic
simulations, some of these choices are preferred over others by their
numerical robustness and accuracy. For instance, in the Extended Conformal
Thin Sandwich formalism [81, 82] used by the Spectral Einstein Code (SpEC)
[83], we can freely specify the conformal metric, the trace of the extrinsic
curvature, and their time derivatives. For BBH simulations, a widely used and
well-tested choice for the conformal metric is a superposition of two BHs in
Kerr-Schild coordinates [84, 85, 86, 87].

Motivated by this approach, we shall test an ansatz for a 4-dimensional BBH
metric constructed as a superposition of the form
$g_{ab}(t)\sim\eta_{ab}+H^{(1)}_{ab}(t)+H^{(2)}_{ab}(t)$, where $\eta_{ab}$ is
the Minkowski background metric and the terms $H^{(n)}_{ab}(t)$, $n=1,2$,
correspond to each BH. Each BH term is boosted with a time-dependent
transformation following the trajectories of the holes. These trajectories can
be accurately described by solving the PN equations of motion [88] for orbital
separations larger than $\sim 10\>M$. The final BBH metric is a simple, time-
dependent, analytical function that we can use as a background spacetime for
MHD simulations. In this section, we show how to build this ansatz.

###  II.1 Kerr black hole in harmonic coordinates

The metric of a spinning BH in Kerr-Schild (KS) coordinates, $\\{t_{\rm
KS},\,x_{\rm KS},\,y_{\rm KS},\,z_{\rm KS}\\}$, is the natural choice for
building a superposed metric since it has the form of a background term plus a
BH term (see Ref. [89, 90]):

$g_{ab}=\eta_{ab}+2Hl_{a}l_{b},$ (1)

where $\eta_{ab}$ is the Minkowski metric in Cartesian coordinates, and the
null covector $l_{a}$ is defined as

$\displaystyle-l_{a}dx^{a}_{\rm KS}:=$ $\displaystyle dt_{\rm
KS}+\frac{r\>x_{\rm KS}+a\>y_{\rm KS}}{r^{2}+a^{2}}dx_{\rm KS}$ (2)
$\displaystyle+\frac{r\>y_{\rm KS}-a\>x_{\rm KS}}{r^{2}+a^{2}}dy_{\rm
KS}+\frac{z_{\rm KS}}{r}dz_{\rm KS},$ (3)

with the spin parameter $a$, the function:

$H:=\frac{2Mr^{3}}{r^{4}+a^{2}z_{\rm KS}^{2}},$ (4)

where $M$ is the mass of the black hole, and the Boyer-Lindquist radius $r$ is
given by

$r^{2}:=\frac{1}{2}(r_{\rm KS}^{2}-a^{2})\left(1+\sqrt{1+\frac{4a^{2}z_{\rm
KS}^{2}}{(r_{\rm KS}^{2}-a^{2})^{2}}}\right),$ (5)

where

$r_{\rm KS}^{2}:=x_{\rm KS}^{2}+y_{\rm KS}^{2}+z_{\rm KS}^{2}.$ (6)

This coordinate system was used in our previous work [56] to build an
approximate spacetime metric for a BBH system in a Keplerian orbit, using a
simple prescription to boost the BHs. Since we want to describe an inspiraling
BBH using PN trajectories, we must use a coordinate system compatible with the
PN gauge. As we explain in the next subsection, we use the PN trajectories in
the standard harmonic coordinate system [88] and, for this reason, we build
our superposition directly in harmonic coordinates 111It was shown recently
[87] that building superposed free initial data in harmonic coordinates,
rather than in the standard Kerr-Schild coordinates, is more accurate for
numerical simulations of BBH mergers. We also found that our superposed metric
in harmonic coordinates is more accurate than the KS gauge..

The Kerr metric has a well-known harmonic coordinate system introduced by Cook
and Scheel [92], which is also horizon penetrating. This is an important
feature of the coordinates for doing GRMHD simulations because the excision
can be placed inside the horizon (we analyze more features of this harmonic
coordinate system in Appendix A) . Using the known transformation from in-
going Kerr to this harmonic coordinate system (c.f. Appendix B in Ref. [68]),
we express the transformation from Kerr-Schild coordinates, $\\{x^{a}_{\rm
KS}\\}$, to Cook-Scheel harmonic coordinates, $\\{x^{a}_{\rm H}\\}$, as:

$\displaystyle t_{\rm KS}=$ $\displaystyle t_{\rm
H}+2M\log{(r-r_{-})}-2M\log(2M),$ (7) $\displaystyle x_{\rm KS}=$
$\displaystyle x_{\rm H}+M\left[\frac{(r-M)y_{\rm H}-ax_{\rm
H}}{(r-M)^{2}+a^{2}}\right],$ (8) $\displaystyle y_{\rm KS}=$ $\displaystyle
y_{\rm H}+M\left[\frac{(r-M)x_{\rm H}+ay_{\rm H}}{(r-M)^{2}+a^{2}}\right],$
(9) $\displaystyle z_{\rm KS}=$ $\displaystyle z_{\rm H}+M\left(\frac{z_{\rm
H}}{r-M}\right),$ (10)

where we have:

$r-M=\sqrt{(Q+W)/2},\quad W:=\sqrt{Q^{2}+4a^{2}z_{\rm H}^{2}},$ (11)

and

$Q:=r_{\rm H}^{2}-a^{2},\quad r_{\rm H}^{2}:=x_{\rm H}^{2}+y_{\rm
H}^{2}+z_{\rm H}^{2}.$ (12)

The space components are thus related by the elegant relation:

$(x^{i}_{\rm KS}-x^{i}_{\rm H})\delta_{ij}(x^{j}_{\rm KS}-x^{j}_{\rm
H})=M^{2}.$ (13)

If we apply this transformation to the Cartesian Minkowski part of the Kerr-
Schild metric (1):

$\eta^{\rm H}_{ab}(a,M)=\frac{\partial x_{\rm KS}^{a^{\prime}}}{\partial
x_{\rm H}^{a}}\frac{\partial x_{\rm KS}^{b^{\prime}}}{\partial x_{\rm
H}^{b}}\eta_{a^{\prime}b^{\prime}},$ (14)

we note that the transformed quantity $\eta^{H}_{ab}(a,M)$ now depends on the
spin and mass of the BHs. However, we can still write this as a flat Cartesian
metric plus a source term:

$\eta^{\rm H}_{ab}(a,M)=\eta_{ab}+M\mathcal{A}_{ab}(a,M),$ (15)

where $\eta_{ab}$ is again the Cartesian Minkowski metric. It can be shown
that $\mathcal{A}_{ab}(a,M)$ is well-behaved at spatial infinity:

$\mathcal{A}\sim 1/r,\quad\text{for }r\rightarrow\infty.$ (16)

The second term of the Kerr-Schild metric can be transformed in the same
manner and is also well behaved at infinity. We conclude that the Kerr metric
in harmonic coordinates can be written as a background plus a BH term,
suitable for superposition, as:

$g_{ab}=\eta_{ab}+M\mathcal{H}_{ab},$ (17)

where $\mathcal{H}_{ab}:=2Hl^{\rm H}_{a}l^{\rm H}_{b}+\mathcal{A}_{ab}$. In
these coordinates, we shall build an effective metric (superposed harmonic PN,
or SHPN) of the form:

$g_{ab}=\eta_{ab}+\phi^{*}_{(1)}\Big{(}M_{(1)}\mathcal{H}^{(1)}_{ab}\Big{)}+\phi^{*}_{(2)}\Big{(}M_{(2)}\mathcal{H}^{(2)}_{ab}\Big{)},$
(18)

where $M_{(1)}$ and $M_{(2)}$ are the masses of each BH,
$\mathcal{H}^{(1)}_{ab}$ and $\mathcal{H}^{(2)}_{ab}$ their corresponding
tensorial functions $\mathcal{H}_{ab}$, and $\phi^{*}_{(1)}$ and
$\phi^{*}_{(2)}$ are transformations that boost the BH terms to describe the
global metric of a BBH system, using a PN approximation for the trajectories.
We show how to build this transformation in the next sub-section.

###  II.2 Moving superposed black holes with PN trajectories

The Kerr metric in harmonic coordinates (17) represents a BH at rest with
respect to an asymptotically inertial frame. To describe a uniformly moving
BH, we can apply a Lorentz boost transformation and convert our coordinates to
boosted coordinates (c.f. Ref. [93]). Physical quantities at spatial infinity
transform as four-vectors in Minkowski space-time. For instance, the
asymptotic observer will measure that a boosted BH has a mass $M_{B}=\gamma
M$, where $M$ is its rest ADM mass. This is simply a frame transformation, and
it does not change gauge-invariant quantities such as the Ricci or Kretschmann
scalars.

Let us suppose that we have a binary system, where the BHs move in an
inspiraling orbit with respect to the origin of a (Cartesian) coordinate
system $\mathcal{O}$, with their trajectories given by:

$s^{a}_{A}(t)=(t,\vec{s}_{A}(t))\equiv(\gamma_{A}\tau_{A},\vec{s}_{A}(\tau_{A})),$
(19)

where $A\in\\{$ BH1, BH${}_{2}\\}$, $\vec{s}_{A}(t)$ is the spatial trajectory
of a BH, $\tau_{A}$ the proper time, and $\gamma_{A}$ the Lorentz factor.
Throughout this work, we assume that the BH spins are (anti-)aligned, which
implies that there is no precession and the orbit lies on the $xy$ plane.

For our approximate BBH metric, we are going to boost two BH terms
$\mathcal{H}_{ab}$ from Eq. (17) and superpose them, as sketched in Eq. (18).
We build the time-dependent boosts as coordinate transformations from the BH
frames $\mathcal{O}^{\prime}\\{X^{a}\\}$, to the (global) center of mass frame
$\mathcal{O}\\{x^{a}\\}$. In $\\{X^{a}\\}$, the BH is at rest and its metric
is locally given by Eq. (17); in the global coordinates $\\{x^{a}\\}$, on the
other hand, the BH is moving according to the worldline (19). This
transformation constitutes a generalized boost since the BHs are not in
uniform motion, i.e., the BH coordinates $\\{X^{a}\\}$ are non-inertial
coordinates.

The natural (pseudo-Cartesian) coordinate system associated with the frame of
an accelerated worldline is called a Fermi Normal coordinate system [94, 95].
These widely used coordinates generalize the boost transformation for time-
dependent velocities (see Ref. [96] for details). Let us say we want to build
this coordinate transformation for a given event $e$ in spacetime (see Figure
1). First, find the proper time for which $s^{a}(\tau)$ is simultaneous to $e$
in the non-inertial BH frame. Then define the time coordinate of the system
$\mathcal{O^{\prime}}$ to be the proper time of the worldline $X^{0}=T=\tau$.
Finally, assume that the hypersurface orthogonal to the worldline is
approximately Euclidean, so the event $e$ described in the global coordinate
system $\mathcal{O}$ is connected with $X^{a}$ as:

$x^{a}=s^{a}(\tau)+X^{i}e^{a}_{(i)}(\tau),$ (20)

where $e^{a}_{(i)}$ are the components of the orthonormal basis carried by the
BH, in global coordinates.

To find the coordinate transformation we need to find the components of the
orthonormal basis of the BH in the global coordinates. Let us assume that the
frame carried by the BH is parallel to the axes of the inertial system
$\mathcal{O}$. Since at each point the BH has a time-dependent velocity with
respect to $\mathcal{O}$, locally we have to boost at each point in time the
spatial frame to compare this with the global frame. The general Lorentz
transformation in the $xy$ plane, given the spatial velocity
$\vec{v}=\beta\vec{n}(t)$ of the BH, can be obtained with the boost
generators, $\vec{K}$, and rapidity, $\xi=\tanh^{-1}(\beta)$, as:

$\displaystyle\Lambda(t)$ $\displaystyle=\exp{(\xi\vec{n}(t)\cdot\vec{K})}$
$\displaystyle=\begin{pmatrix}\gamma&\gamma\beta n_{x}&\gamma\beta n_{y}&0\\\
\gamma\beta n_{x}&1+(\gamma-1)n_{x}^{2}&(\gamma-1)n_{x}n_{y}&0\\\ \gamma\beta
n_{y}&(\gamma-1)n_{x}n_{y}&1+(\gamma-1)n_{y}^{2}&0\\\ 0&0&0&1\end{pmatrix}.$

We use this transformation to boost the spatial basis of the BH. In the BH
coordinates, this basis is simply given by
$e^{a^{\prime}}_{(i)}=\delta^{a^{\prime}}_{i}$ (i.e. the Cartesian spatial
basis). Then, in global coordinates this is given by
$e^{a}_{(i)}=\Lambda^{a}_{i}(t)$.

Figure 1: Diagram of the BH rest frame, $\mathcal{O}^{\prime}$, which is non-
inertial, and the global frame, $\mathcal{O}$. The coordinates $\\{X^{a}\\}$
and $\\{x^{a}\\}$ both describe the event $e$.

The coordinate transformation, using (20), is given by:

$\displaystyle t$ $\displaystyle=\gamma\left[T-\beta(n_{x}X+n_{y}Y)\right],$
(21) $\displaystyle x$
$\displaystyle=s_{x}(t)+X\left[1+(\gamma-1)n_{x}^{2}\right]+Y\left[-(\gamma-1)n_{x}n_{y})\right],$
(22) $\displaystyle y$
$\displaystyle=s_{y}(t)+X\left[-(\gamma-1)n_{x}n_{y}\right]+Y\left[1+(\gamma-1)n_{y}^{2}\right],$
(23) $\displaystyle z$ $\displaystyle=Z.$ (24)

The non-inertial coordinates $X^{a}$ in terms of the global coordinates
$x^{a}$ are easily obtained inverting these equations. Note that this
transformation reduces to a standard Lorentz boost if the trajectory of the BH
is a straight line with uniform rapidity. superposing the two BH terms and
performing this transformation for each term, with worldlines $s^{a}_{1}(t)$
and $s^{a}_{2}(t)$, we have explicitly our SHPN metric:

$\boxed{g_{ab}=\eta_{ab}+M_{1}\Big{(}\frac{\partial
X_{1}^{\overline{a}}}{\partial x^{a}}\frac{\partial
X_{1}^{\overline{b}}}{\partial
x^{b}}\mathcal{H}_{\overline{a}\overline{b}}\Big{)}+M_{2}\Big{(}\frac{\partial
X_{2}^{\overline{a}}}{\partial x^{a}}\frac{\partial
X_{2}^{\overline{b}}}{\partial
x^{b}}\mathcal{H}_{\overline{a}\overline{b}}\Big{)},}$ (25)

where the tensors are transformed through the Jacobian of the coordinates
$X^{a}_{A}(x)$. We still have to supplement the metric with the position,
velocity, and acceleration of the BHs. In the case of a BBH, we can obtain
those solving the PN equations of the system in harmonic coordinates, as we
show in the next section.

###  II.3 Post-Newtonian trajectories for spinning BH binaries

We assume for now that the orbit of the binary has circularized and the system
is well described by the so-called adiabatic approximation [88]. We also
assume that the spins of the holes are (anti-)aligned with the orbit and thus
we ignore orbital precession. In this case, the inspiral is driven by the loss
of binding energy of the orbit, $E$, balanced by the gravitational wave flux
of energy, $\mathcal{F}$, and change in mass $\dot{M}$:

$\dot{E}=-\mathcal{F}-\dot{M},$ (26)

where a dot represents a derivative with respect to the global time $t$. From
this equation, we can obtain the orbital phase, $\Phi(t)$, and separation,
$r_{12}(t)$, of the system. In the case of quasi-circular orbits, the gauge-
dependent separation $r_{12}$ is linked to the orbital frequency through the
relativistic generalization of Kepler’s law [88].

First, we solve for the orbital phase. We rewrite Eq. (26) as two equations in
terms of the (gauge invariant) variable $v:=(Md\Phi/dt)^{1/3}$:

$\displaystyle\frac{dv}{dt}$
$\displaystyle=-\frac{\mathcal{F}(v)+\dot{M}(v)}{dE(v)/dv},$ (27)
$\displaystyle\frac{d\Phi}{dt}$ $\displaystyle=\frac{v^{3}}{M},$ (28)

and replace $E(v)$, $\mathcal{F}(v)$, and $\dot{M}(v)$, with their explicit
expressions at 3.5 PN order for the case of non-precessing binaries, as
presented in Ref. [97].

Following the TaylorT4 scheme [98], we expand the right hand side of Eq. (27)
in a Taylor series to the proper PN order, and integrate it to obtain $t(v)$
(see, e.g., Refs. [99, 100, 101] for Taylor PN approximants). Then, we invert
this quantity, and solve Eq. (28) for $\Phi(t)$. Next, we solve for the
orbital separation $r_{12}(t)$ within numerical accuracy although it is
directly derived from the orbital frequency in the PN approximation. We write
its time derivative in terms of the orbital energy:

$\dot{r}_{12}=\frac{dE/dt}{dE/dr_{12}}\equiv-\frac{\mathcal{F}(r_{12})+\dot{M}(r_{12})}{dE/dr_{12}},$
(29)

and integrate to find:

$t(r_{12})=t_{c}-\int^{r_{12}}_{0}d\tilde{r}_{12}\frac{dE/d\tilde{r}_{12}}{\mathcal{F}(\tilde{r}_{12})+\dot{M}(\tilde{r}_{12})}\,,$
(30)

where $t_{c}$ is the time until coalescence. We replace $E(r_{12})$,
$\mathcal{F}(r_{12})$, and $\dot{M}(r_{12})$ with their explicit expressions,
found in Ref. [88] and references therein, and solve for $t(r_{12})$. Finally,
we invert using a Newton-Rawson method to obtain $r_{12}(t)$. With $\Phi(t)$
and $r_{12}(t)$ we can reconstruct the worldlines (19) of each spinning hole
in harmonic coordinates, as required by our metric (25).

Note that, even though we use the PN approximation to obtain the BH
trajectories, our metric is valid in the inner zone because we are using the
full relativistic BH terms that include the ergosphere and horizon (see
Appendix A). In other words, we are restricted to binaries with separations
larger than $\sim 10M$, but the metric, as we will show now, is accurate in
both inner and near zones.

##  III Analysis of the spacetime metric

In this section, we test the global validity of our SHPN metric (25). The
metric of a BBH system must satisfy Einstein’s field equation in a vacuum, and
thus the Ricci tensor must be zero. In numerical relativity, violations of
Einstein’s equation are tracked using the Hamiltonian and momentum
constraints. Since we intend to use the four-dimensional form of the metric in
our applications, here we focus first on four-dimensional quantities to
quantify deviations from the exact solution. In particular, following Ref.
[60, 68], we investigate the Ricci scalar, $R:=g^{ab}R_{ab}$, where $R_{ab}$
is the Ricci tensor. Violations of the Ricci scalar $R$ are not absolute and,
thus, they are only meaningful when compared with other quantities. For
instance, if we have that $|R(t_{1})|>|R(t_{0})|$ for $t_{1}>t_{0}$, we can
state that the approximate metric has deteriorated or deviated from a vacuum
solution over time. Similarly, we can compare the Ricci scalar of different
systems or at different points in space to assess locally where the metric is
a better approximation to a vacuum solution.

Figure 2: Determinant of superposed metric for different values of the spin
and separation of $r_{12}=20M$. Note that the curves for different spins are
very similar. For comparison, we include the determinant of the matching
metric for $\chi=0$. Figure 3: Absolute value of several metric components for
the superposed and matching metric in the fiducial configuration. Note that
the superposed metric components are much smoother than the matching metric
because there are no transition regions.

We also compare the validity of our solution with the alternative approach
presented in Refs. [60, 68], where an analytical metric is built by stitching
different approximated solutions of Einstein’s equations at three
characteristics zones of a binary compact system, namely, the Inner-Zone (IZ),
the Near-Zone (NZ), and the Far-Zone (FZ) [102, 103]. This so-called
asymptotic matching procedure brings all these different parts into the same
harmonic coordinate system and the global metric can be written as:

$\displaystyle g_{ab}$ $\displaystyle=$ $\displaystyle(1-f_{\rm
FZ})\Bigl{\\{}f_{\rm NZ}\bigl{[}f_{{\rm IZ},1}\,g_{ab}^{\rm(NZ)}+(1-f_{{\rm
IZ},1})\,g_{ab}^{\rm(IZ1)}\bigr{]}$ (31) $\displaystyle+(1-f_{\rm
NZ})\bigl{[}f_{{\rm IZ},2}\,g_{ab}^{\rm(NZ)}+(1-f_{{\rm
IZ},2})\,g_{ab}^{\rm(IZ2)}\bigr{]}\Bigr{\\}}$ $\displaystyle+f_{\rm
FZ}\,g_{ab}^{\rm(FZ)}\,,$

where transition functions $f_{i}$ are used to go from one zone to the other.
This analytical metric, however, is computationally expensive and complex to
handle. The Jacobians required to stitch the different parts of the metric
into the same coordinate system are very long, and many operations are
required to compute them at each timestep. Moreover, for the spinning case,
the matching procedure renders the metric prohibitively expensive for MHD
simulations. In our new approach, we lose some accuracy in comparison with the
matching metric but we gain much more efficiency.

###  III.1 Spacetime scalars

Although the metric is analytical, we compute its spacetime scalars
numerically as it is faster and more practical to incorporate the PN
trajectories. The functional form of the metric is built using Mathematica
[104] and exported to C language in an optimized form. We use then a C-based
code that implements fourth-order finite differences in a Cartesian grid for
the derivatives of every metric function. We analyze and plot the outputs
using Numpy and Matplotlib [105, 106] The convergence analysis of these
methods is presented in Appendix B.

We are interested in using the metric in the inspiraling regime, where the PN
approximation holds, and the system is emitting a significant amount of
gravitational radiation. We explore the characteristics of the system for a
fiducial configuration, with a separation of $r_{12}(t_{0})=20M$, equal BH
masses, and the adimensional spin parameter, $\chi:=a/M$, in the interval
$0<\chi<0.9$.

As a first check of consistency, we analyze the metric determinant
$\sqrt{-g}$. In Figure 2, we plot the determinant for a separation of
$r_{12}=20M$ and several values of the spin parameter $\chi$, along with the
determinant for the matching metric. We see that for all these values, the
determinant for the superposed metric is globally well-behaved, free of
pathologies, and similar to the matching space-time.

In Figure 3, we plot some components of both metrics. It is interesting to
note that the $g_{tt}$ component of the SHPN metric is globally similar to the
matching one, meaning that the effective PN potential of both spacetimes is
much akin [65]. The differences between the two metrics are important in the
transition regions and the Far-Zone. In the latter, the matching metric
incorporates the post-Minkowski background of gravitational waves, while our
new SHPN is asymptotically flat; however, we do take into account the
gravitational radiation losses in the trajectories of the BHs.

Figure 4: Ricci scalar for the matching and superposed metric with
$r_{12}=20M$, equal mass-ratio, and different values of spin. The dashed green
lines denote the location of the horizon and the solid green lines the
location of the ISCO for a non spinning BH. Figure 5: Equatorial logarithmic
plot of the Ricci scalar for the SHPN metric with $r_{12}=20M$, equal masses,
and $\chi=0.9$

In Figure 4, we plot the Ricci scalar of the SHPN metric over the positive $x$
axis at $z=y=0$, for different values of spin, and we compare it with the
Ricci scalar of the matching metric for a binary of the same characteristics.
First, we see that $R$ varies very little under different spin parameters,
consistent with Ref. [68]. Note that the matching metric is better in the IZ
but the violations are worse at the transition regions outside the ISCO, where
the SHPN is smoother and performs better. A good metric accuracy in this
region is an important feature for determining the correct gas dynamics of an
accreting disk near the hole. In Figure 5, we show an equatorial plot of the
Ricci scalar. As expected, the higher violations are concentrated in the
middle region between the BHs and drop sharply with distance. In Figure 6, we
plot the Ricci violations for different mass-ratios $q:=M_{1}/M_{2}$. We find
the values of $R$ depend smoothly on $q$, improving in the middle region for
smaller $q$.

Figure 6: Ricci violations of the SHPN metric, for $r_{12}=20M$, $\chi=0.5$,
and different mass-ratio values.

Besides $R$, we can explore other curvature scalars to assess the global
behavior of the metric. In particular, considering the ADM equations for a
general spacetime, we can define the Hamiltonian constraint $\mathcal{H}$ as:

$\mathcal{H}:={}^{3}R+K-K_{ab}K^{ab}=16\pi\tilde{\rho},$ (32)

where ${}^{3}R$ is the spatial Ricci curvature, $K_{ab}$ the extrinsic
curvature, and $\tilde{\rho}$ the energy density of matter. For our BBH vacuum
metric, a non-zero $\mathcal{H}$ means that the spacetime has “fake mass” due
to the approximation. This will introduce errors in the true gravitational
potential and thus in the geodesic motion of matter. Since we are interested
in using this spacetime as a background scenario for evolving an MHD fluid, it
is important to analyze this quantity and its evolution. We consider the
volume-integrated value of $\mathcal{H}$ as a measure of the total fake mass
introduced by the approximated metric:

$M_{\rm fake}=\frac{1}{16\pi}\int_{\mathcal{V}}\mathcal{H}\>d\mathcal{V}.$
(33)

Considering a cube of radius $r=50M$ around the center of mass, we can track
the evolution of $M_{\rm fake}$ for different orbital separations. As we show
in Figure 7, this fake mass is overall small with respect to the total mass of
the BBH in both SHPN and matching metric but starts increasing exponentially
at $\sim 8M$, where the PN approximation breaks.

Figure 7: Fake mass introduced in the spacetime by the SHPN and matching
metric approximations for different orbital distances.  Figure 8: The Ricci
scalar ($R$), the Hamiltonian ($\mathcal{H}$) and momentum ($\mathcal{M}$)
constraint, and the Kretschmann scalar $\mathcal{K}$ of the SHPN metric, for a
BBH with $r_{12}=20M$, equal mass, and $\chi=0.9$. We can notice that the
Kretschmann scalar follows the decay $\sim 1/r^{6}$ typical for a single BH.
The solid green lines denotes the location of the ISCO for a non-spinning BH.

Finally, in Figure 8, we plot the Ricci scalar, the Hamiltonian constraint,
the square root of the Momentum constraints, $\mathcal{M}$, and the
Kretschmann scalar, $\mathcal{K}:=R_{abcd}R^{abcd}$, for the SHPN metric. We
observe here that the Hamiltonian constraint and the Ricci scalar have similar
behaviors, indicating that the errors of the approximation come essentially
from the fake mass component.

##  IV Superposed metric in GRMHD simulations

Figure 9: Equatorial plot of the rest-mass density $\rho$ in logarithmic scale
for (right) the superposed metric and (left) the matching metric, at
$T=1810M$, which represents approximately $\sim 3$ orbits.

We built our new superposed metric in harmonic coordinates, which allows us to
use accurate PN trajectories and directly compare our simulations with
previous results that use the same gauge. Moreover, the metric is accurate
enough near the BH, allowing us to analyze what happens with the plasma
physics around each BH.

The matching metric approach has been tested in simulations, both as a
background spacetime for MHD [66, 57], and in numerical relativity, [107]. MHD
simulations of BBHs accretion disks with the spinning matching metric, as we
said, are prohibitively expensive. In contrast, the SHPN metric is
computationally cheaper, allowing us to simulate this type of system for the
first time. In this section, we show the results of a full 3D GRMHD simulation
with the SHPN metric using the code Harm3d. In particular, we focus on
comparing an accretion disk simulation using the SHPN metric for non-spinning
BHs and the analog simulation presented in Refs. [23, 25], which uses the non-
spinning matching metric. We will present full details of the spinning BBH
simulation with the new metric in an upcoming work [108].

###  IV.1 GRMHD evolution

Assuming the surrounding gas does not influence the spacetime dynamics, we can
use our superposed metric to simulate the MHD evolution of accretion disks in
a BBH system. For that purpose, we implement the new metric in the GRMHD code
Harm3d [109, 110, 111, 38], which evolves the ideal GRMHD equations in flux-
conservative form, for an arbitrary metric and coordinate system. The
equations of motion are the continuity equation, the local conservation of
energy and momentum, and Maxwell’s equations, which can be written as:

$\partial_{t}{\bf U}\left({{\mathbf{P}}}\right)=-\partial_{i}{\bf
F}^{i}\left({{\mathbf{P}}}\right)+\mathbf{S}\left({{\mathbf{P}}}\right)\ ,$
(34)

where ${\bf P}$ are the primitive variables, ${\bf U}$ are the conserved
variables, ${\bf F}^{i}$ the fluxes, and $\mathbf{S}$ are the source terms.
These can be expressed as

$\displaystyle{\bf P}$ $\displaystyle:=$
$\displaystyle\left[\rho,u,\tilde{u}^{i},B^{i}\right],$ $\displaystyle{\bf
U}\left({{\mathbf{P}}}\right)$ $\displaystyle:=$
$\displaystyle\sqrt{-g}\left[\rho u^{t},\,{T^{t}}_{t}+\rho
u^{t},\,{T^{t}}_{j},\,B^{k}\right]\ ,$ $\displaystyle{\bf
F}^{i}\left({{\mathbf{P}}}\right)$ $\displaystyle:=$
$\displaystyle\sqrt{-g}\left[\rho u^{i},\,{T^{i}}_{t}+\rho
u^{i},\,{T^{i}}_{j},\,\left(b^{i}u^{k}-b^{k}u^{i}\right)\right]\ ,$
$\displaystyle\mathbf{S}\left({{\mathbf{P}}}\right)$ $\displaystyle:=$
$\displaystyle\sqrt{-g}\left[0,\,{T^{\kappa}}_{\lambda}{\Gamma^{\lambda}}_{t\kappa}-\mathcal{F}_{t},\,{T^{\kappa}}_{\lambda}{\Gamma^{\lambda}}_{j\kappa}-\mathcal{F}_{j},\,0\right],$

where ${\Gamma^{c}}_{ab}$ are the Christoffel symbols,
$b^{a}=\left(1/u^{t}\right)\left({\delta^{a}}_{b}+u^{a}u_{b}\right)B^{b}$ is
the magnetic 4-vector projected into the fluid’s comoving reference frame,
$B^{i}$ is the magnetic field in the reference frame of the space normal
hypersurface, $u$ is the internal energy density, $u^{\alpha}$ are the
components of the fluid’s 4-velocity, and $\tilde{u}^{i}$ is the fluid
velocity in the zero-angular-momentum observer (ZAMO) frame. The stress-energy
tensor is written as

$T_{ab}=\left(\rho
h+2p_{m}\right)u_{a}u_{b}+\left(p+p_{m}\right)g_{ab}-b_{a}b_{b}\ ,$ (35)

where $h=1+\epsilon+p/\rho$ is the specific enthalpy, $\epsilon:=u/\rho$ is
the specific internal energy, $p$ is the gas pressure,
$p_{m}=\frac{1}{2}b^{2}$ is the magnetic pressure, and $\rho$ is the rest-mass
density. We include a source term $\mathcal{F}_{\nu}$ into the local energy
conservation equation in order to approximate effects from radiative cooling,
designed to preserve the height ratio of the disk [23]. We assume and ideal
$\Gamma$-law equation of state: $p=(\Gamma-1)\rho\epsilon$, with $\Gamma=5/3$.

Harm3d uses high-resolution shock-capturing methods to integrate the
conservation equations (34). In particular, we use the piece-wise parabolic
reconstruction of primitive variables for the local Lax-Friedrichs fluxes, a
Flux-CT scheme to maintain the solenoidal constraint [112], and a robust
recovery procedure from conserved to primitive variables [113]. The code uses
fourth-order finite differences for spatial derivatives of the metric to find
the Christoffel symbols, and the method of lines for time integration with a
Runge-Kutta method of second-order (more details of the algorithm in Refs.
[111] and [38]).

###  IV.2 Mini-disk dynamics in a binary black hole system: comparison with
previous simulations

In Ref. [25], the matching metric (31) was used to simulate the MHD dynamics
of a circumbinary disk around a BBH and the formation of mini-disks. Ref. [66]
showed, for the first time, that mini-disks in tight binary systems are out of
inflow equilibrium, filling and depleting their mass in less than an orbital
period, showing interesting modulations [57, 26].

We perform a simulation with the same configuration and initial data as Refs.
[66, 57] but switching the matching metric for our superposed metric (31) with
zero spin. The simulation uses a double-fisheye spherical grid [114] that
focuses more cells in the vicinity of the BHs and maintains a spherical
topology at the circumbinary region. The initial separation of the BHs is
$r_{12}=20\>M$ and the initial data for the matter fields are constructed from
a stabilized snapshot of the circumbinary simulation performed in Ref. [38],
with additional quasi-equilibrated mini-disks around each BH (see Ref. [57]
for more details). Since we are starting with the same initial data and grid,
we re-normalize the primitive $B^{i}$ field by $\sqrt{-g_{\rm
match}}/\sqrt{-g_{\rm sup}}$ to maintain the solenoidal constraint of the
field after switching to the new metric.

As shown in Ref. [57], the simulation has an initial transient that lasts for
$\sim 2$ orbits. We thus evolve the system with the new metric for $\sim 3.5$
orbits as was done in [66]. After equilibration, both matching and superposed
metric simulations are in good agreement, as it can be seen, for instance,
from the equatorial density snapshots in Figure 9. For a more quantitative
assessment of both simulations we analyze the evolution of the mass contained
in each mini-disk (Figure 10), defined as the integrated rest-mass density:

$M_{i}=\int_{\mathcal{V}_{i}}\rho u^{0}\>\sqrt{-g}d^{3}x,$ (36)

where we take $\mathcal{V}_{i}$ as a spherical volume between the BH horizon,
$r_{i}=r_{\mathcal{H}}$, and $r_{f}=0.4a(t)$, which is close to the Newtonian
truncation radius [65]. We also investigate the volume integral of the
magnetic energy $b^{2}$ of each mini-disk (see Figure 11).

As we mentioned, for an equal-mass binary, the mini-disks are subject to a
filling and depletion cycle. While the circumbinary lump feeds material to one
of the mini-disks, the plasma in the other BH is completely accreted, and the
BH starves. After the initial transient, we see a remarkable overlap of each
mini-disk mass $M_{i}$ for both simulations. The cycle is evident earlier in
the superposed metric simulation. Since we are using the same initial data for
both simulations (not just the same prescription), the equilibration of the
mini-disks changes in the transient phase for the superposed, as the
equilibrated mini-tori were set up using the matching metric. We also observe
a good agreement in the behavior of the magnetic energy contained in each
mini-disk after the transient.

Figure 10: Mass in each BH’s mini-disk region for the superposed metric
simulation (thick lines) and the matching metric simulation (dashed lines).
Figure 11: Magnetic energy in each BH’s mini-disk region for the superposed
metric simulation (thick lines) and the matching metric simulation (dashed
lines).

##  V Conclusions

We have presented a new approximate solution of Einstein’s field equations for
a spinning BBH in the inspiral regime. We built this solution as a linear
superposition of boosted Kerr BHs in harmonic coordinates, supplemented with
PN trajectories at 3.5 PN order. We compared our new metric with the well-
tested asymptotic matching approach through an analysis of spacetime scalars.
Although the matching approach has better accuracy in some specific regions,
we found that the superposed metric has comparable accuracy, is smoother, and
much cheaper. We also compared the performance of the metric in a GRMHD
simulation using the same setup as previous simulations with the matching
metric. We found that the superposed metric reproduces the same physical
features of the matching metric simulation. We conclude that a superposed
metric is a robust approach for exploring the MHD plasma on BBH systems.

Based on these results, in an upcoming paper, we will analyze the effects of
the BH spins and orbital evolution on the mini-disk dynamics and outflows of a
BBH system embedded in a circumbinary disk [108]. When the black holes rotate,
jets of entirely new characteristics may emerge in a BBH system. Any jet
launched should have a helical structure with a diameter equal to the major
axis of the responsible black hole’s orbit. If each black hole produces a jet,
the two jets may collide or interact. Because jets are intrinsically unsteady,
and the minidisks’ mass accretion rates vary with a phase difference of
$\simeq\pi$ [57], intersection dynamics are expected to be asymmetric in
general even if the black holes have the same mass. This would generate a
whole range of unexplored phenomena, such as periodic non-thermal flares
produced by the collision region [115].

One of the main advantages of our approach is that it assumes very little of
the BBH properties. Although we have restricted to quasi-circular orbits in
this work, implementing a general orbital motion is straightforward because
the superposition does not assume any symmetries in the trajectories. As long
as the trajectories are solutions of the Post-Newtonian equations, the
constraints should remain low. On the other hand, even though we assume that
the spins of the BHs are aligned to the orbital plane, we can easily
generalize this by applying a rotation to the BH term in (25) before the boost
transformation. Moreover, this rotation can be time-dependent to take into
account the Post-Newtonian evolution of spins. In this way, we will be able to
describe the metric of precessing binaries approaching merger in all
generality. The last two points, however, must be tested in the same way we
did here. We leave them for future work. Finally, using the multi-grid
Patchwork MHD code [116], we would have the possibility to handle the entire
parameter space of the system for inspiralling binaries in a computationally
efficient way.

Figure 12: Left: Ergosphere region of Kerr BH with spin $\chi=0.9$ for Kerr-
Schild coordinates (blue) and harmonic coordinates (orange). Note that the
surfaces in harmonic coordinates are more oblique compared with the Kerr-
Schild coordinates. The radius of the singularity (green) is the same for both
coordinate systems. Right: Ergosphere regions for a x-boosted harmonic Kerr BH
with spin $\chi=0.9$ for different velocities ($v/c=0,\,0.1,\,0.5,\,0.9$). The
horizon (red) and singularity (green) are the same in each case, but the
ergosphere region increases with increasing velocity.

###### Acknowledgements.

We thank Carlos Lousto, Eduardo Gutiérrez, and Mark Avara for useful
discussions. L. C., F. L. A, M. C., and B. I acknowledge support from
AST-1028087, AST-1028087, AST-1516150, PHY-1707946 and from NASA TCAN grant
No. 80NSSC18K1488. L. C also acknowledges support from a CONICET (Argentina)
fellowship. S. C. N. was supported by AST-1028087, AST-1515982 and
OAC-1515969, and by an appointment to the NASA Postdoctoral Program at the
Goddard Space Flight Center administrated by USRA through a contract with
NASA. D. B. B. is supported by the US Department of Energy through the Los
Alamos National Laboratory. Los Alamos National Laboratory is operated by
Triad National Security, LLC, for the National Nuclear Security Administration
of U.S. Department of Energy (Contract No. 89233218CNA000001). H. N.
acknowledges support from JSPS KAKENHI Grant Nos. JP16K05347 and JP17H06358.
Computational resources were provided by the Blue Waters sustained-petascale
computing NSF projects OAC-1811228 and OAC-1516125. Blue Waters is a joint
effort of the University of Illinois at Urbana-Champaign and its National
Center for Supercomputing Applications. Additional resources were provided by
the RIT’s BlueSky and Green Pairie Clusters acquired with NSF grants
AST-1028087, PHY-0722703, PHY-1229173 and PHY-1726215. The views and opinions
expressed in this paper are those of the authors and not the views of the
agencies or US government.

##  Appendix A Horizons in harmonic coordinates

To perform GRMHD simulations in a BBH spacetime, we need a clear picture of
how the BH singularities and horizons behave in the chosen coordinates. In
this appendix, we will show how the harmonic coordinates compare with usual
Kerr-Schild coordinates defined in Eq. (1).

We are interested in the ergosphere $r_{\mathcal{E}}$ and the outer horizon
$r_{+}$ of our spacetime. In Boyer-Lindquist (BL) coordinates, these are given
by [89]:

$r^{\rm BL}_{\mathcal{E}}=M+\sqrt{M^{2}-a^{2}\cos^{2}(\theta_{\rm BL})},$ (37)
$r^{\rm BL}_{+}=M+\sqrt{M^{2}-a^{2}}.$ (38)

The radius in harmonic coordinates can be related with the BL radius, $r_{\rm
BL}$, as:

$r^{2}_{\rm H}=(r_{\rm BL}-M)^{2}+a^{2}(1-\cos^{2}(\theta_{\rm BL})),$ (39)

where

$\cos(\theta_{\rm BL})\equiv\frac{z_{\rm H}}{r_{\rm BL}-M},$ (40)

and $r_{\rm BL}=r(x_{\rm H},\,y_{\rm H},\,z_{\rm H})$ is calculated from
equation (5). From these expressions, we can derive parametric equations for
the horizons and ergosphere in harmonic coordinates:

$r^{\rm H}_{\mathcal{E}}=M\sqrt{1-\chi^{2}\Big{(}2\cos^{2}(\theta_{\rm
BL})-1\Big{)}},$ (41) $r^{\rm H}_{+}=M\sqrt{1-\chi^{2}\cos^{2}(\theta_{\rm
BL})}.$ (42)

From Eq. (40), we observe that the harmonic coordinates become singular at
$r^{\rm BL}_{\mathcal{S}}=M$, which means that there is a disk singularity at
$z^{\rm H}=0$ with radius given by:

$r^{\rm H}_{\mathcal{S}}=M\chi.$ (43)

The radius of the horizon in harmonic coordinates shrinks at the poles as
$\chi$ increases, while the radius of the horizon at $z=0$ is fixed at $r^{\rm
H}_{+}(z=0)=M$ for any value of the spin. Then, the distance between the
singularity and the horizon at $z=0$ shrinks with increasing spin as
$\delta_{\rm H}:=M(1-\chi)$. In contrast, in Kerr-Schild coordinates, the
horizon is further away from the singularity, with a separation of
$\delta_{\rm KS}:=\delta_{\rm H}+Mg(\chi)$, where
$g(\chi):=\sqrt{2+2\sqrt{1-\chi^{2}}}-1>0$ (see Figure 12).

To avoid any spurious effect of the inner boundary of the domain, it is
usually placed inside the horizon to mask the singularity of the spacetime. In
this way, the coordinates must be horizon penetrating and the singularity
should be sufficiently far from the horizon. In this regard, the Kerr-Schild
coordinate system is more convenient than harmonic coordinates since the
distance $\delta_{\rm KS}$ is bigger than $\delta_{\rm H}$. For performing
high-spin simulations with the harmonic coordinates, one could artificially
remove the singularity by implementing a modification of the metric inside the
horizon, e.g., modifying the function $r_{\rm BL}(x_{\rm H})$.

Figure 13: Convergence factor of our numerical scheme for different regions.
In the top panel we use $h/M=0.0125$, in the middle panel $h/M=0.1$, and in
the bottom panel $h/M=0.8$.

Note that the Cartesian Kerr-Schild coordinates used here are not the usual
coordinates that accretion disk theorists call ‘Kerr-Schild’ [117, 118]. The
‘accretion-disk Kerr-Schild’ coordinates are a modification of the BL
coordinates that renders the metric horizon-penetrating but maintains the
singularity at $r_{\rm AKS}=r_{\rm BL}=0$. The ‘Cartesian Kerr-Schild’
coordinates that we use here are more common in numerical relativity and
appear in the original work of Kerr [90, 119].

Finally, let us note that our spacetime contains moving BHs, boosted with
respect to the asymptotically flat region. This means that the morphology of
the ergosphere would be different from a static BH and these differences can
be significant for high velocities. As discussed in Ref. [93], even a non-
spinning BH acquires an ergosphere when the BH is boosted. In the case of a
spinning BH, we can see from Figure 12 that the ergosphere increases when the
BH has higher velocities.

##  Appendix B Convergence tests

Since we are using a finite difference scheme for computing the metric and
connection derivatives, we show here the convergence to the analytical
solution of the fourth-order discretization. As explained in Ref. [68], since
the metric spans several length scales, we need different mesh spacing to
resolve the solution. Given a numerical quantity $U$, we explore the
convergence factor, $\mathcal{Q}^{h}(U)$, defined as:

$\mathcal{Q}^{h}(U):=\frac{U^{(4h)}-U^{(2h)}}{U^{(2h)}-U^{(h)}},$ (44)

where $h$ is the size of the mesh spacing and $U^{(k)}$ is the numerical
approximation with spacing $k$. If the method is well behaved, we have [120]:

$\mathcal{Q}^{h}(U)=2^{p}+\mathcal{O}(h).$ (45)

We explore the convergence of the Ricci scalar along the $x$ axis in the
equatorial plane, which is the most relevant region. We calculate the
derivatives of the metric on a uniform Cartesian grid using a fourth-order
finite difference method, so we use $p\equiv 4$ in (45). In Figure 13, we show
convergence for different regions. In the top panel, for the vicinity of the
BH at $[0M,\,20M]$, we use a mesh spacing given by $h/M=0.0125$. Away from the
BH, a high-resolution mesh drops the convergence because of the limited
machine precision to represent numbers (double precision in our case) 222We
can confirm this using a quadrupole precision code, although it is
prohibitively expensive for our current simulations [68], so we change the
mesh to $h/M=0.1$ at $[20M,60M]$ and $h/M=0.8$ at $[60M,200M]$. From Figure 13
we can check that convergence is achieved in the different regions, where we
expect $\mathcal{Q}\sim 16$.

## References

  * [1] K. Gültekin, D. O. Richstone, K. Gebhardt, T. R. Lauer, S. Tremaine, M. C. Aller, R. Bender, A. Dressler, S. M. Faber, A. V. Filippenko, R. Green, L. C. Ho, J. Kormendy, J. Magorrian, J. Pinkney, and C. Siopis, “The M-$\sigma$ and M-L Relations in Galactic Bulges, and Determinations of Their Intrinsic Scatter,” Astrophys. J. , vol. 698, pp. 198–221, June 2009.
  * [2] M. Volonteri, “Formation of supermassive black holes,” The Astronomy and Astrophysics Review, vol. 18, no. 3, pp. 279–315, 2010.
  * [3] P. F. Hopkins and E. Quataert, “How do massive black holes get their gas?,” Monthly Notices of the Royal Astronomical Society, vol. 407, pp. 1529–1564, Jul 2010.
  * [4] T. Naab and J. P. Ostriker, “Theoretical challenges in galaxy formation,” Annual Review of Astronomy and Astrophysics, vol. 55, pp. 59–109, Aug 2017\.
  * [5] F. M. Khan, D. Fiacconi, L. Mayer, P. Berczik, and A. Just, “Swift Coalescence of Supermassive Black Holes in Cosmological Mergers of Massive Galaxies,” Astrophys. J. , vol. 828, p. 73, Sept. 2016.
  * [6] L. Z. Kelley, L. Blecha, and L. Hernquist, “Massive black hole binary mergers in dynamical galactic environments,” , vol. 464, pp. 3131–3157, Jan. 2017.
  * [7] H. Pfister, A. Lupi, P. R. Capelo, M. Volonteri, J. M. Bellovary, and M. Dotti, “The birth of a supermassive black hole binary,” Mon. Not. Roy. Astron. Soc., vol. 471, no. 3, pp. 3646–3656, 2017.
  * [8] F. Pretorius, “Evolution of binary black-hole spacetimes,” Phys. Rev. Lett., vol. 95, p. 121101, Sep 2005.
  * [9] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, “Accurate evolutions of orbiting black-hole binaries without excision,” Phys. Rev. Lett., vol. 96, p. 111101, Mar 2006.
  * [10] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter, “Gravitational-wave extraction from an inspiraling configuration of merging black holes,” Phys. Rev. Lett., vol. 96, p. 111102, Mar 2006.
  * [11] S. Burke-Spolaor, S. R. Taylor, M. Charisi, T. Dolch, J. S. Hazboun, A. M. Holgado, L. Z. Kelley, T. J. W. Lazio, D. R. Madison, N. McMann, et al., “The astrophysics of nanohertz gravitational waves,” The Astronomy and Astrophysics Review, vol. 27, no. 1, p. 5, 2019.
  * [12] W. Del Pozzo, A. Sesana, and A. Klein, “Stellar binary black holes in the lisa band: a new class of standard sirens,” Monthly Notices of the Royal Astronomical Society, vol. 475, pp. 3485–3492, Jan 2018.
  * [13] B. Perera, M. DeCesar, P. Demorest, M. Kerr, L. Lentati, D. Nice, S. Osłowski, S. Ransom, M. Keith, Z. Arzoumanian, et al., “The international pulsar timing array: second data release,” Monthly Notices of the Royal Astronomical Society, vol. 490, no. 4, pp. 4666–4687, 2019\.
  * [14] J. Baker, J. Bellovary, P. L. Bender, E. Berti, R. Caldwell, J. Camp, J. W. Conklin, N. Cornish, C. Cutler, R. DeRosa, et al., “The laser interferometer space antenna: unveiling the millihertz gravitational wave sky,” arXiv preprint arXiv:1907.06482, 2019.
  * [15] S. T. McWilliams, R. Caldwell, K. Holley-Bockelmann, S. L. Larson, and M. Vallisneri, “Astro2020 decadal science white paper: The state of gravitational-wave astrophysics in 2020,” arXiv preprint arXiv:1903.04592, 2019.
  * [16] D. Chapon, L. Mayer, and R. Teyssier, “Hydrodynamics of galaxy mergers with supermassive black holes: is there a last parsec problem?,” Monthly Notices of the Royal Astronomical Society, vol. 429, pp. 3114–3122, Jan 2013\.
  * [17] M. Tremmel, M. Karcher, F. Governato, M. Volonteri, T. R. Quinn, A. Pontzen, L. Anderson, and J. Bellovary, “The romulus cosmological simulations: a physical approach to the formation, dynamics and accretion models of smbhs,” Monthly Notices of the Royal Astronomical Society, vol. 470, pp. 1121–1139, May 2017.
  * [18] M. Colpi, “Massive Binary Black Holes in Galactic Nuclei and Their Path to Coalescence,” , vol. 183, pp. 189–221, Sept. 2014.
  * [19] J. H. Krolik, “Estimating the prompt electromagnetic luminosity of a black hole merger,” The Astrophysical Journal, vol. 709, pp. 774–779, Jan 2010\.
  * [20] C. Roedig, J. H. Krolik, and M. C. Miller, “Observational Signatures of Binary Supermassive Black Holes,” Astrophys. J. , vol. 785, p. 115, Apr. 2014.
  * [21] D. J. D’Orazio, Z. Haiman, and D. Schiminovich, “Relativistic boost as the cause of periodicity in a massive black-hole binary candidate,” Nature (London), vol. 525, pp. 351–353, Sept. 2015.
  * [22] J. H. Krolik, M. Volonteri, Y. Dubois, and J. Devriendt, “Population estimates for electromagnetically distinguishable supermassive binary black holes,” The Astrophysical Journal, vol. 879, p. 110, jul 2019.
  * [23] D. B. Bowen, M. Campanelli, J. H. Krolik, V. Mewes, and S. C. Noble, “Relativistic Dynamics and Mass Exchange in Binary Black Hole Mini-disks,” Astrophys. J. , vol. 838, p. 42, Mar. 2017.
  * [24] G. Ryan and A. MacFadyen, “Minidisks in Binary Black Hole Accretion,” Astrophys. J. , vol. 835, p. 199, Feb. 2017.
  * [25] D. B. Bowen, V. Mewes, M. Campanelli, S. C. Noble, J. H. Krolik, and M. Zilhão, “Quasi-periodic Behavior of Mini-disks in Binary Black Holes Approaching Merger,” , vol. 853, p. L17, Jan. 2018.
  * [26] S. d’Ascoli, S. C. Noble, D. B. Bowen, M. Campanelli, J. H. Krolik, and V. Mewes, “Electromagnetic Emission from Supermassive Binary Black Holes Approaching Merger,” Astrophys. J., vol. 865, no. 2, p. 140, 2018.
  * [27] Y. Tang, Z. Haiman, and A. MacFadyen, “The late inspiral of supermassive black hole binaries with circumbinary gas discs in the LISA band,” , vol. 476, pp. 2249–2257, May 2018.
  * [28] V. Springel, T. D. Matteo, and L. Hernquist, “Black holes in galaxy mergers: The formation of red elliptical galaxies,” The Astrophysical Journal, vol. 620, pp. L79–L82, jan 2005.
  * [29] M. Milosavljević and E. S. Phinney, “The Afterglow of Massive Black Hole Coalescence,” , vol. 622, pp. L93–L96, Apr. 2005.
  * [30] J. E. Pringle, “The properties of external accretion discs,” Monthly Notices of the Royal Astronomical Society, vol. 248, pp. 754–759, 02 1991.
  * [31] A. I. MacFadyen and M. Milosavljević, “An Eccentric Circumbinary Accretion Disk and the Detection of Binary Massive Black Holes,” Astrophys. J. , vol. 672, pp. 83–93, Jan. 2008.
  * [32] T. Bode, R. Haas, T. Bogdanović, P. Laguna, and D. Shoemaker, “Relativistic Mergers of Supermassive Black Holes and Their Electromagnetic Signatures,” Astrophys. J. , vol. 715, pp. 1117–1131, June 2010.
  * [33] C. Palenzuela, L. Lehner, and S. Yoshida, “Understanding possible electromagnetic counterparts to loud gravitational wave events: Binary black hole effects on electromagnetic fields,” Physical Review D, vol. 81, p. 084007, Apr. 2010.
  * [34] B. D. Farris, Y. T. Liu, and S. L. Shapiro, “Binary black hole mergers in gaseous disks: Simulations in general relativity,” Phys. Rev. D, vol. 84, p. 024024, July 2011.
  * [35] T. Bode, T. Bogdanović, R. Haas, J. Healy, P. Laguna, and D. Shoemaker, “Mergers of Supermassive Black Holes in Astrophysical Environments,” Astrophys. J. , vol. 744, p. 45, Jan. 2012.
  * [36] B. D. Farris, R. Gold, V. Paschalidis, Z. B. Etienne, and S. L. Shapiro, “Binary black-hole mergers in magnetized disks: Simulations in full general relativity,” Phys. Rev. Lett., vol. 109, p. 221102, Nov 2012.
  * [37] B. Giacomazzo, J. G. Baker, M. C. Miller, C. S. Reynolds, and J. R. van Meter, “General Relativistic Simulations of Magnetized Plasmas around Merging Supermassive Black Holes,” , vol. 752, p. L15, June 2012\.
  * [38] S. C. Noble, B. C. Mundim, H. Nakano, J. H. Krolik, M. Campanelli, Y. Zlochower, and N. Yunes, “Circumbinary magnetohydrodynamic accretion into inspiraling binary black holes,” The Astrophysical Journal, vol. 755, no. 1, p. 51, 2012.
  * [39] J.-M. Shi, J. H. Krolik, S. H. Lubow, and J. F. Hawley, “Three-dimensional magnetohydrodynamic simulations of circumbinary accretion disks: Disk structures and angular momentum transport,” The Astrophysical Journal, vol. 749, no. 2, p. 118, 2012.
  * [40] D. J. D’Orazio, Z. Haiman, and A. MacFadyen, “Accretion into the central cavity of a circumbinary disc,” Monthly Notices of the Royal Astronomical Society, vol. 436, pp. 2997–3020, Dec. 2013.
  * [41] B. D. Farris, P. Duffell, A. I. MacFadyen, and Z. Haiman, “Binary Black Hole Accretion from a Circumbinary Disk: Gas Dynamics inside the Central Cavity,” The Astrophysical Journal, vol. 783, p. 134, Mar. 2014\.
  * [42] R. Gold, V. Paschalidis, Z. B. Etienne, S. L. Shapiro, and H. P. Pfeiffer, “Accretion disks around binary black holes of unequal mass: General relativistic magnetohydrodynamic simulations near decoupling,” Phys. Rev. D, vol. 89, p. 064060, Mar 2014.
  * [43] B. D. Farris, P. Duffell, A. I. MacFadyen, and Z. Haiman, “Binary black hole accretion during inspiral and merger,” , vol. 447, pp. L80–L84, Feb. 2015.
  * [44] B. D. Farris, P. Duffell, A. I. MacFadyen, and Z. Haiman, “Characteristic signatures in the thermal emission from accreting binary black holes,” , vol. 446, pp. L36–L40, Jan. 2015.
  * [45] J.-M. Shi and J. H. Krolik, “Three-dimensional MHD Simulation of Circumbinary Accretion Disks. II. Net Accretion Rate,” Astrophys. J. , vol. 807, p. 131, July 2015.
  * [46] D. J. D’Orazio, Z. Haiman, P. Duffell, A. MacFadyen, and B. Farris, “A transition in circumbinary accretion discs at a binary mass ratio of 1:25,” , vol. 459, pp. 2379–2393, July 2016.
  * [47] S. C. Noble, J. H. Krolik, M. Campanelli, Y. Zlochower, B. C. Mundim, H. Nakano, and M. Zilhão, “Mass-ratio and magnetic flux-dependence of modulated accretion from circumbinary disks,” arXiv pre-Print 2103.12100, 2021.
  * [48] D. Alic, P. Moesta, L. Rezzolla, O. Zanotti, and J. L. Jaramillo, “Accurate simulations of binary black hole mergers in force-free electrodynamics,” The Astrophysical Journal, vol. 754, no. 1, p. 36, 2012.
  * [49] P. Moesta, D. Alic, L. Rezzolla, O. Zanotti, and C. Palenzuela, “On the detectability of dual jets from binary black holes,” The Astrophysical Journal Letters, vol. 749, no. 2, p. L32, 2012.
  * [50] C. Palenzuela, L. Lehner, and S. L. Liebling, “Dual jets from binary black holes,” Science, vol. 329, no. 5994, pp. 927–930, 2010.
  * [51] C. Palenzuela, T. Garrett, L. Lehner, and S. L. Liebling, “Magnetospheres of black hole systems in force-free plasma,” Physical Review D, vol. 82, no. 4, p. 044045, 2010.
  * [52] B. D. Farris, P. Duffell, A. I. MacFadyen, and Z. Haiman, “Binary Black Hole Accretion from a Circumbinary Disk: Gas Dynamics inside the Central Cavity,” Astrophys. J. , vol. 783, p. 134, Mar. 2014.
  * [53] R. Gold, V. Paschalidis, Z. B. Etienne, S. L. Shapiro, and H. P. Pfeiffer, “Accretion disks around binary black holes of unequal mass: General relativistic magnetohydrodynamic simulations near decoupling,” Physical Review D, vol. 89, p. 064060, Mar 2014.
  * [54] R. Gold, V. Paschalidis, M. Ruiz, S. L. Shapiro, Z. B. Etienne, and H. P. Pfeiffer, “Accretion disks around binary black holes of unequal mass: General relativistic mhd simulations of postdecoupling and merger,” Physical Review D, vol. 90, p. 104030, Nov 2014.
  * [55] V. Paschalidis, J. Bright, M. Ruiz, and R. Gold, “Minidisk dynamics in accreting, spinning black hole binaries: Simulations in full general relativity,” 2021.
  * [56] F. G. L. Armengol, L. Combi, M. Campanelli, S. C. Noble, J. H. Krolik, D. B. Bowen, M. J. Avara, V. Mewes, and H. Nakano, “Circumbinary disk accretion into spinning black hole binaries,” The Astrophysical Journal, vol. 913, no. 1, p. 16, 2021.
  * [57] D. B. Bowen, V. Mewes, S. C. Noble, M. Avara, M. Campanelli, and J. H. Krolik, “Quasi-periodicity of Supermassive Binary Black Hole Accretion Approaching Merger,” Astrophys. J. , vol. 879, p. 76, July 2019.
  * [58] S. C. Noble, B. C. Mundim, H. Nakano, J. H. Krolik, M. Campanelli, Y. Zlochower, and N. Yunes, “Circumbinary Magnetohydrodynamic Accretion into Inspiraling Binary Black Holes,” Astrophys. J. , vol. 755, p. 51, Aug. 2012\.
  * [59] M. Zilhao, S. C. Noble, M. Campanelli, and Y. Zlochower, “Resolving the relative influence of strong field spacetime dynamics and MHD on circumbinary disk physics,” Phys. Rev., vol. D91, no. 2, p. 024034, 2015.
  * [60] B. C. Mundim, H. Nakano, N. Yunes, M. Campanelli, S. C. Noble, and Y. Zlochower, “Approximate black hole binary spacetime via asymptotic matching,” Phys. Rev., vol. D89, p. 084008, 2014.
  * [61] H. Nakano, B. Ireland, M. Campanelli, and E. J. West, “Spinning, Precessing, Black Hole Binary Spacetime via Asymptotic Matching,” Class. Quant. Grav., vol. 33, no. 24, p. 247001, 2016.
  * [62] B. Ireland, B. C. Mundim, H. Nakano, and M. Campanelli, “Inspiralling, nonprecessing, spinning black hole binary spacetime via asymptotic matching,” Phys. Rev. D, vol. 93, p. 104057, May 2016.
  * [63] N. Yunes, W. Tichy, B. J. Owen, and B. Bruegmann, “Binary black hole initial data from matched asymptotic expansions,” Phys. Rev., vol. D74, p. 104011, 2006.
  * [64] N. Yunes and W. Tichy, “Improved initial data for black hole binaries by asymptotic matching of post-Newtonian and perturbed black hole solutions,” Phys. Rev., vol. D74, p. 064013, 2006.
  * [65] D. B. Bowen, M. Campanelli, J. H. Krolik, V. Mewes, and S. C. Noble, “Relativistic Dynamics and Mass Exchange in Binary Black Hole Mini-disks,” Astrophys. J. , vol. 838, p. 42, Mar. 2017.
  * [66] D. B. Bowen, V. Mewes, M. Campanelli, S. C. Noble, J. H. Krolik, and M. Zilhão, “Quasi-periodic Behavior of Mini-disks in Binary Black Holes Approaching Merger,” , vol. 853, p. L17, Jan. 2018.
  * [67] L. Gallouin, H. Nakano, N. Yunes, and M. Campanelli, “Asymptotically Matched Spacetime Metric for Non-Precessing, Spinning Black Hole Binaries,” Class. Quant. Grav., vol. 29, p. 235013, 2012.
  * [68] B. Ireland, B. C. Mundim, H. Nakano, and M. Campanelli, “Inspiralling, nonprecessing, spinning black hole binary spacetime via asymptotic matching,” Physical Review D, vol. 93, no. 10, p. 104057, 2016.
  * [69] M. Campanelli, C. O. Lousto, and Y. Zlochower, “Spinning-black-hole binaries: The orbital hang up,” Phys. Rev., vol. D74, p. 041501, 2006.
  * [70] D. A. Hemberger, G. Lovelace, T. J. Loredo, L. E. Kidder, M. A. Scheel, B. Szilágyi, N. W. Taylor, and S. A. Teukolsky, “Final spin and radiated energy in numerical simulations of binary black holes with equal masses and equal, aligned or antialigned spins,” Physical Review D, vol. 88, no. 6, p. 064014, 2013.
  * [71] J. Healy and C. O. Lousto, “Hangup effect in unequal mass binary black hole mergers and further studies of their gravitational radiation and remnant properties,” Physical Review D, vol. 97, no. 8, p. 084002, 2018.
  * [72] M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krishnan, and D. Merritt, “Spin flips and precession in black-hole-binary mergers,” Phys. Rev., vol. D75, p. 064030, 2007.
  * [73] C. O. Lousto and J. Healy, “Flip-flopping binary black holes,” Phys. Rev. Lett., vol. 114, no. 14, p. 141101, 2015.
  * [74] C. O. Lousto, J. Healy, and H. Nakano, “Spin flips in generic black hole binaries,” Phys. Rev. D, vol. 93, no. 4, p. 044031, 2016.
  * [75] C. O. Lousto, J. Healy, and H. Nakano, “Spin flips in generic black hole binaries,” Physical Review D, vol. 93, no. 4, p. 044031, 2016.
  * [76] C. O. Lousto and J. Healy, “Gravitational Wave Beacons,” Phys. Rev. D, vol. 99, no. 6, p. 064023, 2019.
  * [77] M. Kesden, D. Gerosa, E. Berti, R. O’Shaughnessy, and U. Sperhake, “The gravitational-wave signature of binary black holes in spin-orbit resonances,” in APS April Meeting Abstracts, vol. 2014, pp. C15–004, 2014\.
  * [78] T. W. Baumgarte and S. L. Shapiro, Numerical relativity: solving Einstein’s equations on the computer.  Cambridge University Press, 2010.
  * [79] J. W. York, “Gravitational degrees of freedom and the initial-value problem,” Phys. Rev. Lett., vol. 26, pp. 1656–1658, Jun 1971.
  * [80] J. M. Bowen and J. W. York, “Time-asymmetric initial data for black holes and black-hole collisions,” Phys. Rev. D, vol. 21, pp. 2047–2056, Apr 1980\.
  * [81] J. W. York Jr, “Conformal “thin-sandwich” data for the initial-value problem of general relativity,” Physical review letters, vol. 82, no. 7, p. 1350, 1999.
  * [82] H. P. Pfeiffer and J. W. York Jr, “Extrinsic curvature and the einstein constraints,” Physical Review D, vol. 67, no. 4, p. 044022, 2003.
  * [83] H. P. Lawrence Kidder and M. Scheel, “Spectral Einstein Code.” https://www.black-holes.org/code/SpEC.html, 2020.  [Online; accessed 20-July-2020].
  * [84] R. A. Matzner, M. F. Huq, and D. Shoemaker, “Initial data and coordinates for multiple black hole systems,” Phys. Rev. D, vol. 59, p. 024015, Dec 1998\.
  * [85] G. Lovelace, R. Owen, H. P. Pfeiffer, and T. Chu, “Binary-black-hole initial data with nearly extremal spins,” Physical Review D, vol. 78, no. 8, p. 084017, 2008.
  * [86] G. B. Cook and H. P. Pfeiffer, “Excision boundary conditions for black-hole initial data,” Physical Review D, vol. 70, no. 10, p. 104016, 2004.
  * [87] V. Varma, M. A. Scheel, and H. P. Pfeiffer, “Comparison of binary black hole initial data sets,” Physical Review D, vol. 98, no. 10, p. 104011, 2018\.
  * [88] L. Blanchet, “Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries,” Living Rev. Rel., vol. 17, p. 2, 2014\.
  * [89] E. Poisson, A relativist’s toolkit: the mathematics of black-hole mechanics.  Cambridge university press, 2004.
  * [90] M. Visser, “The kerr spacetime: A brief introduction,” arXiv preprint arXiv:0706.0622, 2007.
  * [91] It was shown recently [87] that building superposed free initial data in harmonic coordinates, rather than in the standard Kerr-Schild coordinates, is more accurate for numerical simulations of BBH mergers. We also found that our superposed metric in harmonic coordinates is more accurate than the KS gauge.
  * [92] G. B. Cook and M. A. Scheel, “Well-behaved harmonic time slices of a charged, rotating, boosted black hole,” Physical Review D, vol. 56, no. 8, p. 4775, 1997.
  * [93] R. F. Penna, “Energy extraction from boosted black holes: Penrose process, jets, and the membrane at infinity,” Physical Review D, vol. 91, no. 8, p. 084044, 2015.
  * [94] F. Manasse and C. W. Misner, “Fermi normal coordinates and some basic concepts in differential geometry,” Journal of mathematical physics, vol. 4, no. 6, pp. 735–745, 1963.
  * [95] E. Poisson, A. Pound, and I. Vega, “The motion of point particles in curved spacetime,” Living Reviews in Relativity, vol. 14, no. 1, p. 7, 2011.
  * [96] B. Mashhoon and U. Muench, “Length measurement in accelerated systems,” Annalen der Physik, vol. 11, no. 7, pp. 532–547, 2002.
  * [97] P. Ajith, M. Boyle, D. A. Brown, B. Brugmann, L. T. Buchman, et al., “The NINJA-2 catalog of hybrid post-Newtonian/numerical-relativity waveforms for non-precessing black-hole binaries,” Class. Quant. Grav., vol. 29, p. 124001, 2012.
  * [98] A. Buonanno, Y.-b. Chen, and M. Vallisneri, “Detecting gravitational waves from precessing binaries of spinning compact objects: Adiabatic limit,” Phys. Rev., vol. D67, p. 104025, 2003.  [Erratum: Phys. Rev.D74,029904(2006)].
  * [99] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, “A Comparison of search templates for gravitational waves from binary inspiral,” Phys. Rev., vol. D63, p. 044023, 2001.  [Erratum: Phys. Rev.D72,029902(2005)].
  * [100] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, “A Comparison of search templates for gravitational waves from binary inspiral - 3.5PN update,” Phys. Rev., vol. D66, p. 027502, 2002.
  * [101] A. Buonanno, B. Iyer, E. Ochsner, Y. Pan, and B. Sathyaprakash, “Comparison of post-Newtonian templates for compact binary inspiral signals in gravitational-wave detectors,” Phys. Rev. D, vol. 80, p. 084043, 2009\.
  * [102] N. Yunes, W. Tichy, B. J. Owen, and B. Brügmann, “Binary black hole initial data from matched asymptotic expansions,” Phys. Rev. D, vol. 74, p. 104011, Nov 2006.
  * [103] N. Yunes and W. Tichy, “Improved initial data for black hole binaries by asymptotic matching of post-Newtonian and perturbed black hole solutions,” Phys. Rev., vol. D74, p. 064013, 2006.
  * [104] W. R. Inc., “Mathematica, Version 12.2.”  Champaign, IL, 2020.
  * [105] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del R’ıo, M. Wiebe, P. Peterson, P. G’erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature, vol. 585, pp. 357–362, Sept. 2020.
  * [106] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering, vol. 9, no. 3, pp. 90–95, 2007.
  * [107] J. Sadiq, Y. Zlochower, and H. Nakano, “Comparing an analytical spacetime metric for a merging binary to a fully nonlinear numerical evolution using curvature scalars,” Physical Review D, vol. 97, no. 8, p. 084007, 2018\.
  * [108] L. Combi, F. G. L. Armengol, M. Campanelli, S. C. Noble, M. J. Avara, J. H. Krolik, and D. B. Bowen, “Accretion onto spinning black hole binaries: mini-disk and outflows,” In prep., 2021.
  * [109] C. F. Gammie, J. C. McKinney, and G. Tóth, “Harm: A numerical scheme for general relativistic magnetohydrodynamics,” Astrophys. J. , vol. 589, pp. 444–457, May 2003.
  * [110] S. C. Noble, C. F. Gammie, J. C. McKinney, and L. Del Zanna, “Primitive variable solvers for conservative general relativistic magnetohydrodynamics,” Astrophys. J. , vol. 641, pp. 626–637, Apr. 2006.
  * [111] S. C. Noble, J. H. Krolik, and J. F. Hawley, “Direct Calculation of the Radiative Efficiency of an Accretion Disk Around a Black Hole,” Astrophys. J. , vol. 692, pp. 411–421, Feb. 2009.
  * [112] G. Tóth, “The $\nabla$ · b = 0 constraint in shock-capturing magnetohydrodynamics codes,” Journal of Computational Physics, vol. 161, no. 2, pp. 605–652, 2000.
  * [113] S. C. Noble, C. F. Gammie, J. C. McKinney, and L. D. Zanna, “Primitive variable solvers for conservative general relativistic magnetohydrodynamics,” The Astrophysical Journal, vol. 641, pp. 626–637, apr 2006.
  * [114] M. Zilhao and S. C. Noble, “Dynamic fisheye grids for binary black hole simulations,” Classical and Quantum Gravity, vol. 31, no. 6, p. 065013, 2014.
  * [115] E. Gutierrez, L. Combi, G. E. Romero, and M. Campanelli, “Repeating flares from jets in supermassive black hole binaries,” In prep., 2021.
  * [116] M. Avara and et al., “Global 3d grmhd simulations of accreting supermassive black hole binaries approaching merger,” In prep., 2021.
  * [117] J. C. McKinney and C. F. Gammie, “A measurement of the electromagnetic luminosity of a kerr black hole,” The astrophysical journal, vol. 611, no. 2, p. 977, 2004.
  * [118] B. J. Kelly, Z. B. Etienne, J. Golomb, J. D. Schnittman, J. G. Baker, S. C. Noble, and G. Ryan, “Electromagnetic emission from a binary black hole merger remnant in plasma: Field alignment and plasma temperature,” Physical Review D, vol. 103, no. 6, p. 063039, 2021.
  * [119] R. P. Kerr, “Gravitational field of a spinning mass as an example of algebraically special metrics,” Physical review letters, vol. 11, no. 5, p. 237, 1963.
  * [120] L. Rezzolla and O. Zanotti, Relativistic hydrodynamics.  Oxford University Press, 2013.
  * [121] We can confirm this using a quadrupole precision code, although it is prohibitively expensive for our current simulations [68].