# Mini-disk accretion onto spinning black hole binaries:
quasi-periodicities and outflows

Luciano Combi 1,2, Federico G. López Armengol 2, Manuela Campanelli 2, Scott
C. Noble 3, Mark Avara 4, Julian H. Krolik 5 Dennis Bowen 6,7 1 Instituto
Argentino de Radioastronomía (IAR, CCT La Plata, CONICET/CIC), C.C.5, (1984)
Villa Elisa, Buenos Aires, Argentina
2 Center for Computational Relativity and Gravitation, Rochester Institute of
Technology, Rochester, NY 14623
3 Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center,
Greenbelt, MD 20771.
4 Institute of Astronomy, University of Cambridge Madingley Road, CB3 0HA
Cambridge, UK
5 Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218
6 Center for Theoretical Astrophysics, Los Alamos National Laboratory, P.O.
Box 1663, Los Alamos, NM 87545.
7 X Computational Physics, Los Alamos National Laboratory, P.O. Box 1663, Los
Alamos, NM 87545<EMAIL_ADDRESS>
###### Abstract

We perform a full 3D general relativistic magnetohydrodynamical (GRMHD)
simulation of an equal-mass, spinning, binary black hole approaching merger,
surrounded by a circumbinary disk and with mini-disks around each black hole.
For this purpose, we evolve the ideal GRMHD equations on top of an
approximated spacetime for the binary that is valid in every position of
space, including the black hole horizons, during the inspiral regime. We use
relaxed initial data for the circumbinary disk from a previous long-term
simulation, where the accretion is dominated by an $m=1$ overdensity called
the lump. We compare our new spinning simulation with a previous non-spinning
run, studying how spin influences the mini-disk properties. We analyze the
accretion from the inner edge of the lump to the black hole, focusing on the
angular momentum budget of the fluid around the mini-disks. We find that mini-
disks in the spinning case have more mass over a cycle than the non-spinning
case. However, in both cases we find most of the mass received by the black
holes is delivered by direct plunging of material from the lump. We also
analyze the morphology and variability of the electromagnetic fluxes and we
find they share the same periodicities of the accretion rate. In the spinning
case, we find that the outflows are stronger than the non-spinning case. Our
results will be useful to understand and produce realistic synthetic light
curves and spectra, which can be used in future observations.

###### Subject headings:

Black hole physics - magnetohydrodynamics - accretion, accretion disks - spins
- jets

††slugcomment: Accepted for publication in ApJ

##  1\. Introduction

When two galaxies merge, a supermassive binary black hole (SMBBH) is expected
to form (Merritt & Milosavljević, 2005). The potential interaction of the new
system with the surrounding gas and the dynamical friction of stars might
shrink the binary separation to sub-parsec scales (Begelman et al., 1980;
Mayer et al., 2007; Escala et al., 2004, 2005; Merritt, 2004, 2006; Dotti et
al., 2007, 2009b; Shi et al., 2012; Sesana & Khan, 2015; Mirza et al., 2017;
Khan et al., 2019; Tiede et al., 2020a; Dittmann & Ryan, 2022). At those
separations, energy and angular momentum is extracted from the system by
gravitational radiation until the BHs merge (Pretorius, 2005; Campanelli et
al., 2006a; Baker et al., 2006). In the near future, gravitational waves from
SMBBH mergers might be observable in the mHZ frequency band by the Laser
Interferometer Space Antenna (LISA, Amaro-Seoane et al., 2017) and by pulsar
timing techniques in the nHz range (Alam et al., 2020; Babak et al., 2016;
Reardon et al., 2016).

Since the environment of these systems is most likely gas-rich, a SMBBH could
also emit electromagnetic radiation through accretion (Barnes & Hernquist,
1992, 1996; Mihos & Hernquist, 1996; Mayer et al., 2007; Dotti et al., 2012;
Mayer, 2013; Derdzinski et al., 2019). At subparsec scales, SMBBHs cannot be
spatially resolved and they might be hard to distinguish from ordinary active
galactic nuclei (AGN). There are many proposed signatures to identify the
presence of a SMBBH using electromagnetic waves, for instance: Doppler
variations due to the orbital motion (D’Orazio et al., 2015), binary
periodicities in the light curves (Valtonen et al., 2006; Graham et al.,
2015b, a; Liu et al., 2019; Saade et al., 2020), interruption of jet emission
(Schoenmakers et al., 2000; Liu et al., 2003) and “spin-flips” of the BH after
merger (Merritt & Ekers, 2002), dual-radio cores (Rodriguez et al., 2006),
profile shifts of broad emission lines (Dotti et al., 2009a; Bogdanović et
al., 2009), a “notch” in the optical/IR spectrum (Roedig et al., 2014; Sesana
et al., 2012), periodicities in the thermal spectrum due to a short residence
time for gas in the disk (Bowen et al., 2019) and X-ray periodicities (Sesana
et al., 2012; Roedig et al., 2014). The feasibility of detecting some of these
signatures depends strongly on binary properties such as mass-ratio and
orbital separation (see also Krolik et al. (2019) for likely source counts).

Because the interstellar gas of the merged galaxy would have a considerable
amount of angular momentum, a circumbinary disk should form around the binary
(Springel et al., 2005; Chapon et al., 2013). For mass-ratios close to one,
the system would present a gap of radius $\sim 2a$ between the binary’s
semimajor axis $a$ and the edge of the circumbinary disk. Accretion then
occurs through two streams from the circumbinary to each black hole. Depending
on its angular momentum, the material from the streams can eventually start
orbiting the BHs, forming “mini-disks”. On the other hand, the time-dependent
quadrupole potential of the binary system can induce a density concentration
in the edge of the circumbinary disk on a small azimuthal range, breaking the
axisymmetric accretion (MacFadyen & Milosavljević, 2008). This feature is
usually referred as the lump (Shi et al., 2012; Noble et al., 2012, 2021).
When the lump is formed, it behaves as a coherent $m=1$ density mode orbiting
the system at a frequency $\sim 0.25\>\Omega_{\rm bin}$, where $\Omega_{\rm
bin}$ is the binary frequency (Noble et al., 2021; Lopez Armengol et al.,
2021). The time-dependence of the gas available for accretion onto the black
holes is then dominated by the lump. If the residence time of the gas in the
mini-disk is shorter (or comparable to) the modulation period of the lump,
then the accretion rate and luminosity of the mini-disks are modulated by the
lump (Bowen et al., 2018; d’Ascoli et al., 2018; Bowen et al., 2019). Whether
this modulation happens depends mostly on the orbital separation (see Section
3)

Numerical simulations are key tools to make accurate models from these highly
non-linear systems. Circumbinary accretion has been largely investigated using
viscous hydrodynamical models in 2D with Newtonian gravity (MacFadyen &
Milosavljević, 2008; D’Orazio et al., 2013; Farris et al., 2014a, b; D’Orazio
et al., 2016; Muñoz & Lai, 2016; Miranda et al., 2017; Derdzinski et al.,
2019; Muñoz et al., 2019; Moody et al., 2019; Mösta et al., 2019; Duffell et
al., 2020; Zrake et al., 2020; Muñoz et al., 2020; Muñoz & Lithwick, 2020;
Tiede et al., 2020b; Derdzinski et al., 2021). In these simulations, the
system can be evolved for $\mathcal{O}(100-1000)$ binary orbits, when a
relaxed stage is reached, and the effects of the initial conditions are
suppressed. However, these simulations rely on artificial sink conditions for
the black hole and their stress model is not self-consistent, i.e. they adopt
an ad-hoc description for the viscosity.

More realistic simulations using MHD have been performed in Newtonian gravity
(Shi et al., 2012; Shi & Krolik, 2015), in approximated GR (Noble et al.,
2012; Zilhão et al., 2015; Lopez Armengol et al., 2021; Noble et al., 2021;
Bowen et al., 2018, 2019), and in full numerical relativity (Gold et al.,
2014, 2014; Farris et al., 2011; Paschalidis et al., 2021; Cattorini et al.,
2021; Giacomazzo et al., 2012). MHD simulations are more computationally
expensive than two-dimensional $\alpha-$viscous simulations because they
demand three-dimensional domains and fine resolutions. Moreover, GRMHD
simulations including both the mini-disks and the circumbinary disk need to
resolve different dynamical timescales, evolve the spacetime metric, and thus
are usually too expensive to evolve for many orbits.

Since the properties of the mini-disks are necessarily tied to the
circumbinary accretion, it is important to perform simulations linking these
regimes in order to obtain a proper global description of the system. A viable
method to accomplish this was presented in Bowen et al. (2018), where a
snapshot of an evolved circumbinary disk simulation (Noble et al., 2012) was
used as initial data for studying mini-disk accretion onto non-spinning black
holes. In this way, Bowen et al. (2019) showed that mini-disks at close binary
separations exhibit a quasi-periodic filling and depletion cycle determined by
the lump and the short inflow time of the mini-disks.

An important property of supermassive binary black holes that many of these
studies miss is the spin (Lin & Papaloizou, 1979). When binary BHs are
approaching merger, spins can have an important impact on the spacetime
evolution: they can alter the orbital motion of the system (Campanelli et al.,
2006c; Hemberger et al., 2013; Healy & Lousto, 2018), induce precession and
nutation (Campanelli et al., 2007), repeatedly flip their sign (Lousto &
Healy, 2015; Lousto et al., 2016a, b; Lousto & Healy, 2019) and even tilt the
orbital orientation (Kesden et al., 2014). Spins also have a key role in the
accretion of matter into BHs (Krolik et al., 2005) and binary BHs. For
instance, accretion rate per unit mass near the circumbinary disk’s inner edge
depends on the spin, altering the mass profile in the inner part of the disk
(Lopez Armengol et al., 2021). Moreover, the character of the flow within the
mini-disks depends on the ratio between the radius at which they are fed and
the radius of the innermost stable circular orbit (ISCO) (Bowen et al., 2018;
Gold et al., 2014; Paschalidis et al., 2021), which is strongly dependent on
the spin.

On the other hand, spinning black holes are expected to launch electromagnetic
outflows (De Villiers & Hawley, 2003; Hirose et al., 2004; McKinney & Gammie,
2004; De Villiers et al., 2005; Krolik et al., 2005; Tchekhovskoy, 2015). The
generation of a net Poynting flux in a binary BH system approaching merger has
been modeled within General Relativistic Force-Free Electrodynamics (GRFFE)
(Palenzuela et al., 2010a, b; Neilsen et al., 2011; Moesta et al., 2012) and
ideal GRMHD (Giacomazzo et al., 2012; Kelly et al., 2017, 2021; Farris et al.,
2012; Gold et al., 2013, 2014; Paschalidis et al., 2021). In the force-free
regime, given a homogeneous plasma threaded by a constant magnetic field, the
Blandford-Znajek (BZ) mechanism (Blandford & Znajek, 1977) proved to operate
efficiently around each BH, leading to a pair of collimated jets in a double-
helical structure that coalesces after merger (Palenzuela et al., 2010a). The
dynamics of binary BHs in a homogeneous medium have also been investigated in
the ideal GRMHD regime, where the inertia of the plasma is taken into account;
in this case, accretion of the plasma onto the BHs leads to a steeper growth
of the initial magnetic field during the inspiral, via compression and
magnetic winding, resulting in higher luminosities (Giacomazzo et al., 2012).
At the same time, the inertia of the (homogenous) plasma falling onto the BHs
interferes with the propagation of the electromagnetic flux, causing less
collimated magnetic structures compared with GRFFE (see also Kelly et al.,
2017, 2021). Further, Farris et al. (2012); Gold et al. (2013, 2014);
Paschalidis et al. (2021) performed GRMHD simulations that include the
circumbinary disk in the domain and found the development of net Poynting
fluxes emitted from the polar regions of each BH that coalesce at larger
distances.

There are only a few GRMHD simulations of circumbinary accretion into spinning
SMBBHs. Very recently, Lopez Armengol et al. (2021) presented the first long-
term circumbinary GRMHD simulation of spinning black holes (with the inner
cavity excised), and Paschalidis et al. (2021) presented the first numerical
relativity simulation of mini-disks around spinning BHs, confirming previous
expectations that the ISCO plays a key role in the mini-disk mass (Bowen et
al., 2018, 2019; Gold et al., 2014). A natural next step is to consider more
realistic simulations where both mini-disks and a properly relaxed
circumbinary disk around spinning BBHs are taken into account. This is
important to make accurate predictions of the light curves and spectra of
these systems.

In this work, we present the results of a GRMHD simulation of mini-disk
accretion around spinning black holes of spins $a=0.6~{}M$ aligned with the
orbital angular momentum. We evolve the ideal GRMHD equations on top of a
binary BH spacetime that is moving on a quasi-circular orbit starting at $20M$
separation. We use an approximate BBH spacetime that uses Post-Newtonian (PN)
trajectories for the BHs in the inspiral regime, but is valid at every space
point, including the BH horizons. As initial data for the plasma, we use a
steady-state snapshot of a circumbinary disk simulation performed in Noble et
al. (2012). We compare our new simulation with a previous non-spinning
simulation (Bowen et al., 2018, 2019) that uses the same initial data.

The paper is structured as follows: In Section 2, we present the simulations
set-up, the numerical methods, the spacetime approximation that we use, and
the initial data. In Section 3, we present our results. First, we give an
overview of the main features of the system and the relation with previous
simulations. In Section 3.2 we investigate the accretion rate, inflow time,
and mass evolution of the mini-disks, and in Section 3.3 we analyze the mini-
disks structure, the specific angular momentum distribution, and the azimuthal
density modes. Then, in Section 3.4 we analyze the outflows, the magnetized
structure of the system, and the variability of the Poynting flux, comparing
spinning and non-spinning results. Finally, in Section 4 we discuss some of
the implications of our results and in Section 5 we summarize the main points
and results of the paper.

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.

##  2\. Simulation set-up

We evolve the equations of ideal GRMHD in the spacetime of a binary black hole
system using the finite-volume code Harm3d. Our goal is to analyze the effects
of the black hole spin in the minidisks. For that purpose, we compare two
simulations of an equal-mass binary black hole, with and without spins. The
non-spinning simulation, denoted as S0, was performed in Bowen et al. (2019)
using an analytical metric built by matching different spacetimes (Mundim et
al., 2014). We perform a new simulation, denoted S06, with BHs having spins of
$\chi=0.6$ aligned with the orbital angular momentum of the binary, using an
approximate analytical metric described in Combi et al. (2021) (a concise
summary is given in Sec. 2.3). For both simulations, we use the same grid and
initial data in order to have a faithful comparison (see Section 2.3 for
details).

###  2.1. GRMHD equations of motion

Assuming that the plasma does not influence the spacetime, we evolve the ideal
GRMHD equations in the dynamical binary BH metric described in Section 2.3.
The equations of motion are given by the conservation of the baryon number,
the conservation of the energy-momentum tensor, and Maxwell’s equations with
the ideal MHD condition:

$\nabla_{a}(\rho u^{a})=0,\quad\nabla_{a}{T^{a}}_{b}=\mathcal{F}_{b},$
$\nabla_{a}\>^{*}F^{ab}=0,\quad u_{a}F^{ab}=0,$ (1)

where $\rho$ is the rest-mass density, $F^{ab}/\sqrt{4\pi}$ is the Faraday
tensor111Following Noble et al. (2009), we absorb the factor $1/\sqrt{4\pi}$
in the definition of the tensor $F^{ab}$., $u^{a}$ is the four-velocity of the
fluid, $\mathcal{F}_{b}$ is the radiated energy-momentum per 4-volume unit,
and the MHD energy-momentum tensor is

$T^{ab}:=(\rho h+b^{2})u^{a}u^{b}+(P+\frac{1}{2}b^{2})g^{ab}-b^{a}b^{b},$ (2)

where $h:=(1+\epsilon+P/\rho)$ is the specfic enthalpy, $\epsilon$ is the
specific internal energy, $P$ is the pressure, $b^{a}:=\>^{*}F^{ab}u_{b}:$ is
the four-vector magnetic field, and $b^{2}:=b^{a}b_{a}$ is proportional to the
magnetic pressure, $p_{\rm m}:=b^{2}/2$. We follow Noble et al. (2012, 2009)
and write these coupled equations of motion in manifest conservative form as

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

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

$\displaystyle\mathbf{P}:$ $\displaystyle=$
$\displaystyle[\rho,u,\tilde{u}^{j},B^{j}]\ ,$ (4) $\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^{j}\right]\ ,$ (5) $\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^{j}-b^{j}u^{i}\right)\right],$ (6)
$\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],$
(7)

where $g$ is the determinant of the metric, ${\Gamma^{\lambda}}_{\alpha\beta}$
are the Christoffel symbols, $u:=\rho\epsilon$ is the internal energy,
$\tilde{u}^{j}:=u^{j}-g^{tj}/g^{tt}$ is the velocity relative to the normal
spacelike hypersurface, and $B^{j}:=\>^{*}F^{it}$ is the magnetic field, which
is both a conserved and a primitive variable 222We denote the magnetic field
in the frame of normal observers (proportional to the constrained-transported
field) as $B^{i}$, while we denote the magnetic field in the frame of the
fluid as
$b^{a}=\left(1/u^{t}\right)\left({\delta^{a}}_{\nu}+u^{a}u_{b}\right)B^{b}$..
We close the system with a $\Gamma$-law equation of state,
$P=(\Gamma-1)\rho\epsilon$, where we set $\Gamma=5/3$.

The source term in the energy-momentum conservation ensures that part of the
dissipated energy caused by MHD turbulence is converted to radiation that
escapes from system. We assume radiation is removed from each cell
independently of all the others, isotropically in the fluid frame. In this
way, we set $\mathcal{F}_{a}={\cal L}_{c}u_{\beta}$, where $\mathcal{L}_{c}$
is the cooling function. We use the prescription used in Noble et al. (2012)
for the rest-frame cooling rate per unit volume:

$\mathcal{L}_{c}=\frac{\rho\epsilon}{t_{\mathrm{cool}}}\left(\frac{\Delta
S}{S_{0}}+\left|\frac{\Delta S}{S_{0}}\right|\right)^{1/2}\ ,$ (8)

where $t_{\mathrm{cool}}$ is the cooling timescale, where the disk radiates
away any local increase in entropy, $\Delta S:=S-S_{0}$, where
$S:=p/\rho^{\Gamma}$ is the local entropy. Our target entropy, $S_{0}=0.01$,
is the initial entropy of each accretion disk in the simulation. The timescale
$t_{\mathrm{cool}}$ is determined by the local fluid orbital period, following
the prescriptions in Bowen et al. (2017); d’Ascoli et al. (2018).

Table 1

Simulation | S06 | S0
---|---|---
Spin parameter [$\chi$] | 0.6 | 0.0
BH1,2 mass [$M_{\rm BH}$] | 0.5 |
Mass-ratio [$q$] | 1 |
Final time [$t_{f}$] | $8000\>M$ | $6000\>M$
Final separation [$r_{12}(t_{f})$] | $16.6\>M$ | $17.8\>M$
# Orbits | $15$ | $12.5$
Init. separation [$r_{12}(0)$] | $20\>M$ |
Init. total mini-disk mass [$M_{0}$] | $20$ |
Average orbital period [$T_{\mathcal{B}}$] |  | $530\>M$
Lump orbital frequency [$\Omega_{\rm lump}$] | $0.28\>\Omega_{\mathcal{B}}$ |
ISCO radius [$r_{\rm ISCO}$] | $2.82\>M_{\rm BH}$ | $5.0\>M_{\rm BH}$
Truncation radius [$r_{\rm trunc}$] | $0.4\>r_{12}$ |
Grid [$(x^{1}\times x^{2}\times x^{3})$] |  | $(600,160,640)$
Physical Size [($r_{\rm min},r_{\rm max}$)] |  | $(2M,260M)$

Note. — Physical and grid parameters of both non-spinning and spinning
simulations.

###  2.2. Code details, grid, and boundary conditions

We solve equations (3) using the high-resolution, shock-capturing methods
implemented in Harm3d. Following Noble et al. (2012), we use a piecewise
parabolic reconstruction of the primitive variables for the local Lax-
Friedrichs flux at each cell interface and the Flux CT method to maintain the
solenoidal constraint (Tóth, 2000). Once the numerical fluxes are found, the
equations are evolved in time using the method of lines with a second-order
Runge-Kutta method. The primitive variables are recovered from the evolved
conserved variables using the 2D method developed in Noble et al. (2006).

Figure 1.— Initial data used in the simulation with quasi-equilibrated mini-
toris around the black holes. In white thin lines we plot the warped spherical
grid every each 50 cells.

Figure 2.—  Left Panel: Rest-mass density snapshot of the fluid in the
equatorial plane for S06. Right panel: Rest-mass density snapshot of the fluid
in the meridional plane corrotating with the (second) black hole. White stream
lines represent the comoving magnetic field projected on the meridional plane.

The grid and boundary conditions used in this simulation are the same as in
Bowen et al. (2018) and Bowen et al. (2019). We use a time-dependent, double
fish-eye, warped spherical grid, centered in the center of mass of the binary
system, developed in Zilhão & Noble (2014) (see full details of the grid used
in Bowen et al. (2019)). The maximum physical size of the grid is set to
$r_{\rm max}=13\>r_{12}(0)$, containing the circumbinary disk of RunSE in
Noble et al. (2012), that we use as initial data. We use outflow boundary
conditions on the radial ($x^{1}$) boundaries, demanding the physical radial
velocity $u^{r}$ to be oriented out of the domain; if not, we reset the radial
velocity to zero and solve for the remaining velocity components. Poloidal
coordinates ($x^{2}$) have reflective, axisymmetric boundary conditions at the
polar axis cutout and the azimuthal coordinates ($x^{3}$) have periodic
boundary conditions.

The resolution is given by $N_{r}\times N_{\theta}\times
N_{\phi}=600{\times}160{\times}640$ cells. The shape of the grid in the
circumbinary region matches the grid used in Noble et al. (2012) and is
sufficient to resolve the magnetorotational instability (MRI) in the
circumbinary disk. Because of our polar grid resolution and off-grid-center
location of the BHs in the spherical grid, our configuration does not include
a full 32 cells per scale height in the mini-disks on the side farthest from
the center-of-mass.

###  2.3. Binary BH spacetime and initial data

The spacetime of the binary black hole is approximated by superposing two Kerr
spacetimes on a Minkowski background. We describe it in terms of harmonic
coordinates. The metric can be written schematically as

$g_{ab}=\eta_{ab}+M_{1}\mathcal{H}^{1}_{ab}(\vec{x}_{1},\vec{v}_{1})+M_{2}\mathcal{H}^{2}_{ab}(\vec{x}_{2},\vec{v}_{2}),$
(9)

where $\eta_{ab}$ is the Cartesian Minkowski metric, $\mathcal{H}^{A}_{ab}$ is
the boosted black hole term for $A=1,2$, $M_{A}$ is the mass, and
$\\{\vec{x}_{A}(t),\vec{v}_{A}(t)\\}$ are the position and velocity of the
black hole. The trajectories are obtained by solving the Post-Newtonian
equations of motion for a spinning binary BH in quasi-circular motion at 3.5
PN order. In Combi et al. (2021), we showed that this analytical metric
constitutes a good approximation to a vacuum solution of Einstein’s field
equations for a binary approaching merger , see also East et al. (2012), Varma
et al. (2018) and Ma et al. (2021) for similar approaches in the context of
numerical relavity . The metric (9) is computationally efficient, compared
with previous approaches, and easy to handle for different parameters. In the
non-spinning simulation, the spacetime was represented by an semi-analytical
metric built by stitching different approximate solutions of Einstein’s
equation (Mundim et al., 2014). In Combi et al. (2021), we compared the
matching and superposed metrics evolving two GRMHD simulations for the non-
spinning case and we found that they are completely equivalent in this
context. Moreover, we analyze the spacetime scalars and integrated Hamiltonian
constraints for each one and we found that (a) they remain small and well-
behaved up to a separation of $10~{}M$, and (b) the constraints remain
invariant when we change the BH spin, and thus no pathologies are introduced
by the spin .

The level at which these constraints are violated is comparable to the very
low level achieved in the numerical simulations performed in Zlochower et al.
(2016), who used the matching metric as initial data for evolving Einstein’s
equations. Constraint violations in this evolution, which can be damped using
the CCZ4 scheme, mainly introduce deviations to the trajectories
(eccentricity) and errors in the masses and spins of the BHs. In our
analytical metric, however, there are no such dynamical effects since we solve
the trajectories using the PN approximation and the BH masses/spins are fixed.
The small constraint violations, relative to the mass of the BHs, might
produce small errors in the gas dynamics, but these are washed out by MHD
turbulence in our simulation (Zilhão et al., 2015).

Although the metric uses PN trajectories and thus is only valid in the
inspiral regime, it is mathematically well-defined at every point in space,
including the horizons of the black holes; hence, no artificial sink terms or
large excisions are needed in the evolution. We apply, however, a mask inside
the horizon to avoid the singularity of each black hole. In particular, we do
not evolve the hydro fluxes inside the masked region and we set to zero the
magnetic fluxes. This allows us to evolve the induction equations in the whole
domain and preserve the solenoidal constraints.

For this simulation, we use an equal-mass black hole binary,
$M_{1}=M_{2}=M/2$, with an initial separation of $r_{12}(0)=20\>M$, where
relativistic effects are important (Bowen et al., 2017) and the orbit is
shrinking due to gravitational radiation. In S06 we set the spins of the black
holes perpendicular to the orbital plane (i.e. no precession) with a moderate
value of $\chi_{i}\equiv a_{i}/M=0.6$. Since spin couples with the orbital
motion, the trajectories of the holes change with respect to a non-spinning
system. In particular, the inspiral is delayed because of the hang-up effect
(Campanelli et al., 2006b), and the orbital frequency increases (see Figure
3).

Figure 3.— Properties a binary black hole system with aligned spins of
$\chi=0.0$ (blue dashed line), $\chi=0.6$ (solid orange line), and $\chi=0.9$
(green dashed line). In the top panel, we plot the radial velocity,
$\dot{r}_{12}$, normalized with its initial value, in the middle panel the
orbital separation $r_{12}$, and in the bottom panel the period of the orbit,
$P$.

In order to compare with the previous simulation, we take the same initial
data for the matter fields as in Bowen et al. (2018) . We start with a
snapshot of the circumbinary disk from Noble et al. (2012), previously evolved
for $5\times 10^{4}M$ ($\sim 80$ orbits). At this time, the disk is in a
turbulent state and accretion into the cavity is dominated by a $m=1$ density
mode, the so-called lump, that is orbiting at the inner edge of the
circumbinary disk. In Noble et al. (2012), a zero spin PN metric in harmonic
coordinates was used to evolve the system. As shown in Lopez Armengol et al.
(2021), the bulk properties of the circumbinary disk are not sensitive to the
spin, even for high values, so it is a good initial state for our spinning
simulation. We interpolate this data onto our grid and we initialize two mini-
disks inside the cavity, see Figure 1. We then clean magnetic divergences
introduced by the interpolation to the new grid using a projection method as
explained in Bowen et al. (2018).

###  2.4. Diagnostics

Properties of the circumbinary disk and other global properties of the system
are better analyzed in the center of mass coordinates. On the other hand, to
analyze mini-disk properties such as fluxes, we compute quantities on the
(moving) BH frame. We define the BH frame at a given time slice with a boosted
coordinate system centered at each BH (see Combi et al. (2021)); we denote
these BH coordinates with a bar,
$\\{\bar{t},\bar{r},\bar{\theta},\bar{\phi}\\}$. We notices that all our
diagnostic are written in the harmonic coordinate gauge . Fluxes and other
local properties in this frame are computed in post-process, interpolating the
global grid into a spherical grid centered in the BH with the Python package
naturalneighbor 333https://github.com/innolitics/natural-neighbor-
interpolation that implements a fast Discrete Sibson interpolation (Park et
al., 2006).

Weighted surface averages of an MHD quantity $\mathcal{Q}$ with respect to a
quantity $\sigma$ are defined as

$\langle\mathcal{Q}\rangle_{\sigma}:=\frac{\int dA\>\sigma\>\mathcal{Q}}{\int
dA\>\sigma},$ (10)

where $dA:=d\theta d\phi\sqrt{-g}$. A time-average of a surface-average is
defined as

$\langle\langle\mathcal{Q}\rangle\rangle:=\frac{1}{\Delta
t}\int^{t_{i}}_{t_{f}}dt\langle\mathcal{Q}\rangle,$ (11)

where we always sum over a given time interval after the initial transient.

##  3\. Results

###  3.1. Overview of the system and previous studies

In the steady-state of an equal-mass SMBBH, a lump orbits the edge of the
circumbinary disk at an average frequency $\langle\Omega_{\rm
lump}\rangle=0.28~{}\Omega_{\mathcal{B}}$ (Noble et al., 2012; Shi et al.,
2012; D’Orazio et al., 2013, 2016; Farris et al., 2015b, a; Lopez Armengol et
al., 2021; Noble et al., 2021) modulating the accretion into the cavity. When
one of the BHs passes near the lump, it peels off part of lump’s inner edge,
forming a stream that feeds the black hole with a beat frequency of
$\Omega_{\rm beat}:=\Omega_{\mathcal{B}}-\langle\Omega_{\rm lump}\rangle\sim
0.72\>\Omega_{\mathcal{B}}$. This stream is almost ballistic, formed by fluid
particles with relatively low angular momentum (Shi & Krolik, 2015). This
material can start orbiting the black hole as it approaches, forming a mini-
disk. The maximum size of the mini-disk is determined by the tidal truncation
radius of the binary, or Hill’s sphere. The residence time of matter in the
mini-disks is determined by the ratio of the truncation radius and the radius
of the ISCO. At close relativistic separations, such as the ones here, the
mini-disks will be out of inflow equilibrium with the circumbinary lump
accretion, and thus the masses oscillate quasi periodically in a filling and
depletion cycle (Bowen et al., 2018, 2019).

For relativistic binaries, the tidal truncation radius is approximately at
$r_{t}\sim 0.4\>r_{12}(t)$ (Bowen et al., 2017), similar to the Newtonian
value, estimated to be $\sim 0.3\>r_{12}$ (Paczynski, 1977; Papaloizou &
Pringle, 1977; Artymowicz & Lubow, 1994). Since spin is a second order effect
in the effective potential (Lopez Armengol et al., 2021), mild values of spin
do not change the truncation radius significantly for binary separations
greater or close to $20M$. The most relevant difference between our two
simulations, S06 and S0, is the location of the ISCO: $r_{\rm
ISCO}(\chi=0.6)=2.82\>M_{\rm BH_{i}}$ for S06, and $r_{\rm
ISCO}(\chi=0.0)=5.0\>M_{\rm BH_{i}}$ for S0 (both in given here in harmonic
coordinates). The smaller ISCOs of the spinning black holes allow material
with lower angular momentum to maintain circular orbits closer to the BH
instead of plunging in directly (see Section 3.3).

In the following sections, we analyze how the size of the ISCO plays a role in
the accretion rate, inflow time, and periodicities, as well as in the
structure of the mini-disks. We will also examine how the presence of an
ergosphere in the spinning case helps the black holes to produce more Poynting
flux. We will also analyze how the variability of the fluxes is connected with
the variability of the accretion rate, dominated by the lump.

###  3.2. Mass evolution, accretion rate, and inflow time

We start the analysis by calculating the integrated rest-mass of each mini-
disk, defined as

$M:=\int^{r_{t}(t)}_{r_{\mathcal{H}}}dV\>\rho u^{t},$ (12)

where $r_{t}(t):=0.4\>r_{12}(t)$ is the truncation radius, $r_{\mathcal{H}}$
is the BH horizon, and $dV:=\sqrt{-g}d^{3}x$ is the volume element. In both
simulations there is an initial transient due to the initial conditions that
lasts approximately $\sim 3$ orbits for S0 and $\sim 4$ orbits for S06 (see
Figure 4). Both simulations start with two quasi-equilibrated mini-tori around
the holes, with a specific angular momentum distribution adapted specifically
for non-spinning black holes, following the prescription described in Bowen et
al. (2017). Since we are also using this initial data for the spinning
simulation, the initial tori have an excess of angular momentum, making the
transient slightly longer in S06. We analyze each simulation after this
transient, marked in the plots as a vertical line. As a time unit, we use
$T_{\mathcal{B}}=530\>M$, the average binary period of the non-spinning
simulation.

Figure 4.— Upper panel: the mass fraction evolution $M_{i}$ of each minidisk
for S06 (solid lines) and S0 (dashed lines), where we define $M_{\rm
tot}:=M_{1}(t)+M_{2}(t)$. Lower panel: total mass evolution for S06 (solid
line) and S0 (dashed line), where $M_{0}:=M_{1}(0)+M_{2}(0)$. The black line
indicates the end of the transient phase.

We find that although the mini-disks of S06, like those in S0, go through a
filling-depletion cycle, the mini-disks around spinning BHs in S06 are more
massive than in S0 by a factor of $2$ through most of the evolution (Figure
4), although they both follow the same decay. When the mass fraction of a
mini-disk is more than $50\%$ of the total mass, we say that the disk is in
its high state; otherwise, it is in its low state. The cycle of the mass
fraction is similar in both simulations, although marginally smaller in
amplitude for S06. The frequency of the cycle is associated with the orbital
frequency and thus is higher in the spinning case, as can be seen plainly
after $\sim 8$ orbits (upper panel colorblue of Fig. 4). On the other hand, at
$t=11\>T_{\mathcal{B}}$, we observe a slight increase of mass in the system.
Because the lump grows and oscillates radially around the cavity (Lopez
Armengol et al., 2021), it generates stronger accretion events onto the black
holes with a lower frequency rate.

We also analyze the accretion rate evolution at the horizon in the black hole
rest-frame, which is given by

$\dot{M}:=\oint_{r_{\mathcal{H}}}d\bar{A}\>u^{\bar{r}}\rho.$ (13)

Here, and in the remainder of this paper, overbars indicate (harmonic)
coordinates whose origin is the center of one of the black holes. We plot the
sum of the accretion rate in each BH, $\dot{M}_{\rm Tot}$, for each simulation
in Figure 5. We find that the accretion rate evolution is overall similar
after the transient in both S06 and S0. Moreover, the time-dependence of the
mini-disk mass is, to a first approximation, a smoothed version of the
accretion rate’s time-dependence, e.g. see Figure 6 for BH1 in S06.

Figure 5.— Total accretion rate evolution $\dot{M}_{\rm Tot}\equiv\dot{M}_{\rm
BH_{2}}+\dot{M}_{\rm BH_{2}}$ in S06 (solid lines) and S0(dashed lines)

The differences in the masses per cycle are closely related to the inflow time
of particles in the mini-disk. It is useful then to define an (Eulerian)
inflow time as the characteristic time for a fluid element to move past a
fixed radius $r$ (Krolik et al., 2005). On average, this can be defined as

$t_{\rm inflow}^{-1}:=\frac{1}{\bar{r}}\Big{\langle}\langle
V^{\bar{r}}\rangle_{\rho}\Big{\rangle},$ (14)

where $V^{\bar{r}}:=u^{\bar{r}}/u^{\bar{t}}$ is the transport velocity. In
Figure 7, we show the inflow time as a function of coordinate radius for BH1,
averaged in time for the first part (solid lines) and second part (dashed
lines) of both simulations. Inside the truncation radius, the inflow time is
consistently longer in S06 than in S0, generally by tens of percent.

Figure 6.— Accretion rate ($\dot{M}_{\rm BH_{1}}$) in red solid lines and mass
($M_{\rm BH_{1}}$) in black dot-dashed lines for a mini-disk around BH1 in
S06. Both quantities were rescaled for plotting.

As the binary shrinks, the inflow time at the truncation radius diminishes
from $\sim 0.41\>T_{\mathcal{B}}$ to $\sim 0.31\>T_{\mathcal{B}}$ for the
spinning simulation. The average inflow time at the truncation radius is well
below the beat period, $T_{\rm beat}=1.3T_{\mathcal{B}}$, on which the mini-
disk refills. Our mean inflow time is also significantly shorter than typical
inflow timescales of Keplerian orbits around single black holes (Krolik et
al., 2005). This suggests that accretion in our mini-disks is driven by
different mechanisms than single BH disks. Moreover, the similar measures of
accretion rates at the horizon for S06 and S0 (see again Figure 5), seem to
imply a shared accretion mechanism between spinning and non-spinning systems
(see Sec. 3.3).

Figure 7.— Inflow time as a function of radius in harmonic coordinates for BH1
in S06 (maroon) and S0 (dark blue). Solid lines are time averages over the
first half of the simulation, while dot-dashed lines are time averages of the
last half of the simulation. Thin lines show the instantaneous inflow time
every $300M$ for the spinning case. The maroon (dark blue) vertical line
represents the ISCO of the spinning (non-spinning) black hole. The vertical
black line is the initial truncation radius of the binary. Figure 8.— Power
spectral density of the mini-disk’s masses for S06 (upper panel) and S0 (lower
pannel) using a Welch algorithm with a Hamming window size and a frequency of
$10M$. The confidence intervals at $3\sigma$ are shown as shadowed areas.

The quasi-periodic behavior of the system can be described through a power
density spectrum (PSD) of the mini-disk’s masses. In Figure 8 we show the PSD,
normalized with the power of the highest peak within the simulation. Most
features of the system’s quasi-periodicities are shared in both simulations
and were described in detail in Bowen et al. (2019). We find, still, some
interesting differences. The PSDs for $M_{1}$ and $M_{2}$ taken individually
peak at the beat frequency in both simulations. The PSD of the total mass of
the mini-disks, $M_{1}+M_{2}$, has a peak at $2\Omega_{\rm beat}$ in S0, while
the latter is severely damped in S06. Indeed, if the individual masses vary
with a characteristic frequency $\Omega_{\rm beat}$, and these are out of
phase with the same amplitude, we expect their sum to vary with $2\Omega_{\rm
beat}$. In S06, however, the inflow time of the mini-disks is larger and the
depletion period of a mini-disk briefly coexists with the filling period of
the other mini-disk, reducing the variability of the total mass. On the other
hand, the beat frequency is slightly higher for S06 ($\Omega_{\rm
beat}=0.71\Omega_{\mathcal{B}}$) than S0 ($\Omega_{\rm
beat}=0.68\Omega_{\mathcal{B}}$) as the orbital frequency of the spinning BHs
is higher. Further, since we evolved the binary for longer, in S06 we find a
more prominent amplitude at $\sim 0.20\Omega_{\mathcal{B}}$. We can associate
this low-frequency power to the radial oscillations of the lump around the
cavity that produces additional accretion events (Lopez Armengol et al.,
2021). We observe this frequency is closely related but different from the
orbital frequency of the lump at $0.28\Omega_{\mathcal{B}}$, which could be
related to the orbital frequency increasing during the inspiral (notice we
measure the PSD at fixed frequencies).

Finally, because we use a spherical grid with a central cutout, we cannot
analyze the effects of the sloshing of matter between mini-disks (Bowen et
al., 2017). To estimate how much mass we lose through the cutout, we compute
the accretion rate at the inner boundary of the grid. This mass loss
constitutes only $5\%$ of the total mass accreted by the BHs throughout the
simulation, although the instantaneous accretion can be close to $20\%$ of the
accretion onto a single BH. We do not expect this small mass loss to alter the
main conclusions of this work, namely, the differences between mini-disks in
spinning and non-spinning BBH.

###  3.3. Structure and orbital motion in mini-disks

In this section, we analyze in detail how the structure of the mini-disks in
S06 compares with mini-disks in S0. We first focus on how the spin changes the
surface density distribution and the azimuthal density modes in the mini-
disks. We then investigate the angular momentum of the fluid and how it
compares with the angular momentum at the ISCO.

Figure 9.— Surface density snapshot for S06 (upper row) and S0 (lower row) at
$t=4000M$ and $t=4060M$ respectively, where the phase of the binary is the
same in both simulations. White dashed lines indicate the truncation radius
and solid white lines indicate the ISCO. The sense of rotation of the binary
is counter-clockwise. Figure 10.— Surface density average in the azimuthal
ranges $\Delta\phi_{1}=(\pi/4,3\pi/4)$ (positive $y_{\rm BH}$-axis) and
$\Delta\phi_{2}=(5\pi/4,7\pi/4)$ (negative $y_{\rm BH}$-axis) for BH1 in
S06(upper panel) and S0(lower panel). Solid lines represent a time average on
the high state of the cycle, while dot-dashed lines represent a time average
over the low state. For reference, we indicate the direction of the orbital BH
velocity. Dashed vertical lines indicate the position of the ISCO.

The surface density is defined as $\Sigma(t,r,\phi):=\int
d\theta\rho\sqrt{-g}/\sqrt{-\sigma}$, where we use
$\sqrt{-\sigma}=\sqrt{g_{\phi\phi}g_{rr}}$ as the surface metric of the
equatorial plane. In Figure 9, we plot the surface density in both S06 and S0,
for the same orbital phase at the $7$th orbit. In this plot, the mini-disk
around BH1 (right side) is in the peak of the mass cycle. In both simulations
we can clearly notice the lump stream plunging directly into the hole. There
is also circularized gas orbiting BH1 in both simulations, but there is much
more of it in S06. On the other hand, we observe that BH2 (left side), in its
low state, has a noticeable disk structure in S06, while the material is
already depleted for S0.

We can quantify these differences in structure by computing the average
surface density over two ranges of $\bar{\phi}$, as measured in the BH frame,
representing the front and back of the mini-disk with respect to the orbital
motion. We define

$\langle\Sigma(r,t)\rangle:=\frac{\int_{\Delta\bar{\phi}}d\bar{l}\>\Sigma}{\int_{\Delta\bar{\phi}}d\bar{l}},$
(15)

where
$d\bar{l}:=d\bar{\phi}\sqrt{g_{\bar{\phi}\bar{\phi}}(\bar{\theta}=\pi/2)}$. In
Figure 10, we plot $\langle\langle\Sigma(r)\rangle\rangle$ for
$\Delta\bar{\phi}_{1}=(\pi/4,3\pi/4)$ and
$\Delta\bar{\phi}_{2}=(5\pi/4,7\pi/4)$, averaging in time over the high and
low state of the mini-disk separately. In the high-state, the mini-disk
accumulates more material at the front while the back of the mini-disk is
flatter. In the low-state, both simulations show a flatter profile, with a
slightly higher density at the front. The asymmetry between front and back
arises because of the orbital motion of the black holes, capturing and
accumulating the stream material as they orbit. In S06, the density profile is
steeper near the ISCO for high and low states. In the outer part of the mini-
disk, the slope of the surface density for both S0 and S06 have a similar
profile, indicating a common truncation radius. The density is higher in S06
by a factor of $\sim 2$ in both states.

Figure 11.— Azimuthal density modes for BH1 in S06 (upper panel) and S0 (lower
panel) normalized with the zero mode $D_{0}$. The black vertical line
represents the end of the transient in S06.

Another important property of mini-disks in relativistic binaries is the
presence of non-trivial azimuthal density modes. When the accreting stream of
the lump impacts the mini-disk, it generates a pressure wave, forming strong
spiral shocks (Bowen et al., 2018). This induces an $m=1$ density mode in the
mini-disk that competes with the $m=2$ mode excited by the tidal interaction
of the companion black hole. The spiral wave patterns can be analyzed
decomposing the mini-disk rest-mass density in azimuthal Fourier modes (Zurek
& Benz, 1986), $\rho(\bar{\phi})\equiv\sum_{m}D_{m}\exp{(-im\bar{\phi})}$,
where

$D_{m}:=\int^{r_{t}(t)}_{r_{\mathcal{H}}}d\bar{V}\rho\exp{(-im\bar{\phi})}.$
(16)

Let us compare these modes in S0 and S06 for BH1. From Figure 11, we observe
that both simulations share common features. In both simulations, the mini-
disks are mainly dominated by $m=1$ modes, followed closely by $m=2$ modes. We
also observe important $m=3,4$ contributions. The modes are excited in the
high state of the cycle, where the mini-disks increase their mass and the
stream is accreted onto the black holes. In S0, the amplitudes of the modes
are noticeably larger than S06, while in the latter the amplitudes grow as the
system evolves. The growth of modes in S06 is correlated with the mass
decrease of the mini-disks. This behavior could indicate that, as the mini-
disks become less massive, the density modes are really representing the
single-arm stream of the lump that plunges directly into the hole. Azimuthal
modes grow to large amplitudes in the high phase of a mini-disk only when the
material orbiting the BH is less or equally massive than the material with low
angular momentum that plunges from the lump, occurring around
$10T_{\mathcal{B}}$ in S06 (see Figure 13 and discussion below). This is
another consequence of the disk-like structure of the mini-disks surviving for
longer time in the spinning case.

Figure 12.— Specific angular momentum as a function of radius for S06 (upper
panel) and S0 (lower panel) for both BHs. The time-averages are in solid
lines, and the individual values are the very thin lines. The Keplerian value
is plotted in dashed green lines.

Our analysis so far indicates that a considerable amount of the matter in the
mini-disk region is plunging directly from the lump to the black hole. In
order to further analyze the orbital motion of the fluid in the mini-disk, we
compute the density-weighted specific angular momentum in the BH frame:

$\langle\bar{\ell}\rangle_{\rho}:=\langle-
u_{\bar{\phi}}/u_{\bar{t}}\rangle_{\rho}.$ (17)

In Figure 12, we show the time-average of $\langle\bar{\ell}\rangle_{\rho}$
for both BHs and both simulations. In absolute terms,
$\langle\bar{\ell}\rangle_{\rho}$ is nearly the same for both the spinning and
non-spinning cases, with the spinning case only slightly greater. This is
because the specific angular momentum of the material that falls into the
cavity is essentially determined by the stresses at the inner edge of the
circumbinary disk. These stresses are determined by binary torques and the
plasma Reynolds and magnetic stresses (Shi et al., 2012; Noble et al., 2012).
Indeed, in Lopez Armengol et al. (2021) we found that these quantities depend
weakly on spin outside the cavity. On the other hand, their relation to their
respective Keplerian (circular orbit) values, $\bar{\ell}_{K}(\bar{r},\chi)$,
is quite different because they depend strongly on the spin. In S06, the
distribution of the angular momentum tracks closely the Keplerian value. For
S0, the behavior is always sub-Keplerian on average.

Let us assume that the angular momentum distribution of the circumbinary
streams is independent of spin for a fixed binary separation and mass-ratio.
In that case, our simulation data indicate that the angular momentum with
which the streams arrive at the mini-disk would be greater than the ISCO
angular momentum when the BH spin is $\chi>0.45$. This estimate could serve as
a crude criterion for determining whether mini-disks form in relativistic
binaries.

Figure 13.— Sub-Keplerian and (super-)Keplerian components of the mass for BH1
in S06 (upper panel) and S0 (lower panel)

We can also use the specific angular momentum to distinguish the material in
the mini-disk with high angular momentum that manages to orbit the black hole
from the low angular momentum part that plunges in. To do so, we recompute the
mass as in equation (12), taking fluid elements with
$\bar{\ell}<\bar{\ell}_{K}$ and $\bar{\ell}\geq\bar{\ell}_{K}$ separately. In
Figure 13 we plot the evolution of the sub-Keplerian and super-Keplerian mass
components for BH1 in S06 and S0. In S06 after the initial transient, a little
more than half of the mass comes from relatively high angular momentum fluid.
As the system inspirals, however, the truncation radius decreases, and the
masses of these two components become nearly equal. In S0, on the other hand,
most of the fluid has relatively low angular momentum. This sub-Keplerian
component has roughly the same mass in S06 and S0, while the mass of the high
angular momentum component of the fluid is much greater in S06, as expected.

Although a fair amount of the mass in the mini-disk has relatively high
angular momentum and manages to orbit the black hole in S06, the accreted mass
onto the BH, in both simulations, is always dominated by the low angular
momentum part that plunges directly. To demonstrate this, we compute the
average accretion rates for low and high angular momentum particles as we did
with the mass. Figure 14 shows that the total accretion rate onto the BH has a
flat radial profile in both S06 and S0, with very similar average values.
Accretion by low angular momentum particles dominates at all radii, although
the high angular momentum contribution becomes comparable to the low angular
momentum one near the ISCO for S06.

Figure 14.— Time averaged accretion rate in BH1 for S06 (solid lines) and S0
(dashed lines) considering particles with low (blue) and high (red) angular
momentum. The vertical dashed blue and dot-dashed red lines mark the ISCO for
S0 and S06, respectively.

We can also compute the density-weighted specific energy, $E:=\langle-
u_{\bar{t}}\rangle_{\rho}$, the mass-weighted sum of rest-mass, kinetic, and
binding energy for individual fluid elements. As can be seen in Figure 15, on
average, fluid in the mini-disks around the spinning black holes is more bound
than in the non-spinning case. On the other hand, fluid in both S06 and S0 is
more bound than particles on circular orbits. Near the ISCO, the specific
energy drops sharply inward in both cases, as is often found when accretion
physics is treated in MHD: stress does not cease at the ISCO when magnetic
fields are present.

Figure 15.— Density-weighted specific energy $E=-u_{\bar{t}}$, averaged in
time as a function of radius for BH1 in S06 and S0. The dashed line represents
the specific energy for a geodesic particle in circular motion in each case.
Vertical lines indicate the location of the ISCO.

When the mini-disk is in its high state, the spiral shocks heat up the mini-
disks and increase their aspect ratio, $h/r$. Given our cooling prescription,
the entropy is kept close to its original value, regulating the aspect ratio
to $h/r=0.1$. However, the gas scale height increases dramatically where the
lump stream impacts the mini-disk. At the peak of the accretion cycle, the
aspect ratio rises to $h/r\sim 2$, but the gas cools before the next accretion
event.

###  3.4. Electromagnetic and hydrodynamical fluxes

In this section we analyze the extraction of energy from the system in two
forms: outward electromagnetic luminosities, arising from a Poynting flux, and
unbound material. In particular, we analyze how these energy fluxes change
with spin and how their variability is characterized by the same periodicity
as the accretion.

The electromagnetic luminosity from each mini-disk evaluated in the BH frame
is:

$L_{\rm EM}(t,\bar{r})=\oint_{\bar{r}}d\bar{A}\>\mathcal{S}^{\bar{r}},$ (18)

where the Poynting flux is $\mathcal{S}^{\bar{i}}:=(T_{\rm
EM})_{\>\>\bar{t}}^{\bar{i}}$. In Figure 16, we plot the EM luminosity as a
function of time for S06 and S0, evaluated on spheres of radius
$\bar{r}=10\>M$ that follow each black hole. These luminosities are normalized
to the average accretion rate $\langle\dot{M}\rangle=0.002$, so they are
equivalent to the rest-mass efficiency of the jets. The most noteworthy
element in Figure 16 is that the EM luminosity is an order of magnitude larger
in S06 than S0. In S0, there are no clear long-term trends while in S06, there
is a secular growth of $L_{\rm EM}$ until $7.5T_{\mathcal{B}}$, when it starts
declining.

Figure 16.— EM luminosity evolution in the BH frame from a sphere at $r=10M$
for S06 (upper panel) and S0 (lower panel).

The instantaneous efficiency $\eta=L_{\rm EM}(t)/\dot{M}(t)$ increases during
the first several orbits, saturating at $\approx 0.05$ in the case of S06, but
an order of magnitude lower for S0, as can be seen in Figure 17. The
efficiency in S06 with spin parameter 0.6 is rather larger than $0.013$, the
value found by De Villiers et al. (2005) for the efficiency of a single BH
with specific angular momentum of 0.5 surrounded by a statistically time-
steady disk. Both S0 and S06 present a similar secular growth of $\eta(t)$
until $10T_{\mathcal{B}}$, where they slightly drop and plateau. Measuring
fluxes in the comoving frame, we found that non-spinning black holes produce
negligible EM luminosity, which is consistent with the fundamental idea of the
Blandford-Znajek mechanism (Blandford & Znajek, 1977; Komissarov, 2001) and
has been confirmed in many simulations (McKinney & Gammie, 2004; De Villiers
et al., 2005; Hawley & Krolik, 2006; Tchekhovskoy, 2015). However, because the
BHs are orbiting around the center of mass this could power additional EM
fluxes (Neilsen et al., 2011; Palenzuela et al., 2010b). In order to capture
the electromagnetic fluxes from the entire binary, we move to the center of
mass frame and calculate the Poynting flux through a sphere of radius
$r=100\>M$ surrounding the binary system.

In Figure 19, we plot this quantity for both S06 and S0 as a function of
retarded time $t-r/\langle v\rangle$, where $\langle v\rangle$ is the mean
velocity of the outflow. We also plot the sum of the Poynting fluxes around
each mini-disk, measured in the BH comoving frame. We notice that the fluxes
in S0 are in average five times larger in the center of mass frame compared
with the BH frame, suggesting that there is a kinematic contribution to
$L_{\rm EM}$ from the orbital motion of the black holes and possibly from the
circumbinary disk. The luminosities in S06, on the other hand, differ between
frames only by tens of percent, suggesting that the rotation of the black hole
dominates the extraction of EM energy here 444We have also checked the
Poynting luminosity as a function of radius and we observed that, in average,
it has varations of around $30\%$ far from the source, which is expectable as
we lose resolution and the outflow interacts with the atmosphere floor of the
simulation.. During our simulations, the speed of the black holes remains
fairly constant at $v\sim 0.1$ but, closer to the merger, the orbital speed
might enhance significantly the electromagnetic luminosity, as seen e.g. by
(Palenzuela et al., 2010b; Farris et al., 2012; Kelly et al., 2020).

In both S06 and S0 the EM luminosities are variable, and in both a Fourier
power spectrum reveals a periodic modulation at the frequency of the
circumbinary disk’s inner edge, i.e., the “lump” frequency. However, in S06
there is an additional modulation of similar amplitude at twice the beat
frequency, proving it is tied closely to accretion, see Figure 22.

Figure 17.— Instantaneous efficiency $\eta$ of the Poynting flux for S06
(solid lines) and S0 (dashed lines) in the BH frame measured at $\bar{r}=10M$.

Figure 18.— Meridional plot of a time average Poynting scalar for BH1 in S06
(left) and in S0 (right). The black hole is at $x\sim 10M$ and the center of
mass is at $x=0M$. The red lines represent the division between bound and
unbound material, while the dot-dashed white lines represent the magnetically
dominated material. Figure 19.— Evolution of the total Poynting flux measured
in the BH frame (dashed lines) and in the (inertial) center of mass frame at
$100~{}M$(solid lines) for both S06 and S0. For the center of mass fluxes, we
use the retarded time $t-r/\langle v\rangle$ to account for the delay.

To further study the spatial distribution of the electromagnetic flux, in
Figure 18, we plot a meridional slice of the Poynting scalar
$\mathcal{P}:=\mathcal{S}^{\bar{i}}\mathcal{S}_{\bar{i}}$ for both S06 and S0,
averaged in time over half an orbit at $4000~{}M$. The dotted-dashed white
lines are defined by $b^{2}/\rho=1$, containing the regions that better
approximate the magnetically dominated region, while the red solid lines are
defined by the region where $hu_{t}=-1$, where the fluid becomes unbound. As
expected, S06 has a much more prominent Poynting jet structure than S0, which
has almost zero jet power. The poloidal distribution of $\mathcal{P}$ in S06
around the black hole has a parabolic shape, with most of the flux being
emitted at mid-latitudes. Each individual jet shape is similar to those found
around single BHs (Nakamura et al., 2018). In both cases, we can also notice
the strong (bound) Poynting fluxes generated by magnetic stresses in the disk
at the equatorial plane.

Figure 20.— Time average of Poynting scalar $\mathcal{P}$ projected on a
sphere of radius 60 M for spinning (left sphere) and non-spinning (right
sphere) for unbound elements of fluid. Figure 21.— Hydro luminosity as a
function of time measured in the center of mass frame for S06 (green lines)
and S0 (blue lines) at different radii.

In Figure 20, we plot $\mathcal{P}$, averaged over half an orbit in the
corotating frame of the binary, for a sphere of radius $r=60~{}M$ at the
center of mass. In S06 the Poynting flux has a double cone structure that
extends to larger polar angles than does the Poynting flux in S0, likely
because of interaction between the two jets. On the other hand, in S0, the
Poynting flux is distributed on a uniform ring around the axis of the binary.
Unfortunately, our simulation coordinates require a cutout in the grid
covering the polar axis running through the center of mass. This prevents an
accurate study of the interaction of the jets. Nonetheless, because this
cutout is small ($2^{\circ}$) compared to the angular size of the jets, we are
still able to pick out important features.

Besides the electromagnetic fluxes, there are also hydrodynamical fluxes from
the system. We define the hydro luminosity as the integral of the energy flux
component of the hydrodynamic stress tensor minus the contribution of the
rest-mass energy flux at that radius in order to get the ‘usable’ energy flux
(Hawley & Krolik, 2006):

$L_{\rm H}(t,r):=\Big{(}\oint_{r}dA\>(T_{\rm H})^{r}_{t}\Big{)}-\dot{M}_{\rm
jet}$ (19)

In Figure 21, we plot the hydro luminosity in the center of mass frame as a
function of time for both S0 and S06, at different radii, normalized by the
averaged value of the accretion rate, so that these luminosities, too, can be
described in the language of rest-mass efficiency. In the integral over flux
(Eqn. (19)), we include only fluid elements that are both unbound according to
the Bernoulli criterion ($-hu_{t}>1$) and moving outward ($u^{r}>0$) (cf. De
Villiers et al. (2005)). Like the EM luminosity in S06, but not S0, the hydro
luminosities are modulated for both spin cases at the “lump” frequency and at
twice the beat frequency. Figure 22 shows the Fourier power spectrum for the
spinning case. We observe in Figure 21 a secular growth of the hydro fluxes in
S06 while in S0 they remain rather constant, with an average efficiency of
$\sim 2\%$. Also like the EM case, the hydro energy flux is considerably
greater in S06 than S0, but by a factor $\sim 5$. Such a contrast resembles
the differences between the hydro efficiencies measured in spinning and non-
spinning single BH simulations (Hawley & Krolik, 2006).

We caution, however, that the luminosities measured at $100~{}M$ may not be
the luminosities received at infinity. Energy can be easily converted from EM
to hydro or vice versa. Here we quote the values at $100~{}M$ because they are
the largest we can measure within our grid. Nonetheless the comparison between
S06 and S0 demonstrates clearly that spinning BBH are much more efficient at
creating coherent outflows carrying energy both electromagnetically and
hydrodynamically than non-spinning BBH.

Figure 22.— Power spectral density of the hydro ($L_{\rm H}$) and EM
luminosities ($L_{\rm EM}$) for S06 (thick lines) and S0 (dashed lines) at
$100~{}M$ using a Welch algorithm with a Hamming window size and a frequency
of $10~{}M$. The confidence intervals at $3\sigma$ are shown as shadowed areas
for S06. The two main peaks are given by twice the beat frequency,
$2\Omega_{\rm beat}=1.4\Omega_{\rm bin}$, and the lump accretion periodicity
$\sim 0.22\Omega_{\rm bin}$

##  4\. Discussion

Although a few candidates have been identified, the existence of supermassive
black hole binaries has not been confirmed. The direct detection of their
gravitational waves by LISA or pulsar timing arrays (PTA) remains at least a
decade into the future. Nevertheless, upcoming wide-field surveys such as the
Vera C. Rubin Observatory, SDSS-V, and DESI, may discover many SMBBH
candidates through their electromagnetic emission.

In order to confirm the presence of a SMBBH, we need to build accurate models
and predictions of their electromagnetic signatures. Our GRMHD simulations
will be useful for this purpose: as a next step, in Gutiérrez et al. (2021),
we use these simulations to extract light curves and spectra using ray-tracing
techniques (Noble et al., 2007; d’Ascoli et al., 2018) with different
radiation models and different masses. The results in this paper constitute
the foundations to interpret the underlying physics of those predictions.

Circumbinary and mini-disk accretion onto an equal-mass binary system has been
largely studied in the past in the context of 2D $\alpha-$viscous simulations.
These simulations are particularly good for analyzing the very long-term
behavior of the system, evolving sometimes for $1000$ orbits. Close to the
black holes and at close separations, however, the inclusion of 3D MHD and
accurate spacetime dynamics becomes necessary in order to describe the proper
mechanisms of accretion and outflow. 2D $\alpha-$disk simulations are not able
to include spin effects and most of them do not include GR effects (see,
however, Ryan & MacFadyen (2017)). On the other hand, in this work we analyze
the balance of hydro accretion from the circumbinary streams and conventional
accretion from the internal stresses of the mini-disk; to properly model the
latter, we need MHD. Moreover, the presence of a proper black hole, and its
horizon, makes the accretion processes entirely self-consistent without adding
adhoc sink conditions as used in Newtonian simulations (see, however, Dittmann
& Ryan (2021)). Finally, 3D MHD simulations are necessary to model
magnetically-dominated regions and jets. The connection of the accretion and
the production of electromagnetic luminosity was one of the main motivations
of this work, and impossible to analyze in 2D hydro simulations.

Recently, Paschalidis et al. (2021) presented GRMHD simulations of a system
similar to the one analyzed in this paper: equal-mass, spinning binary black
holes approaching merger. It is then interesting to compare our results and
highlight the differences with their model and analysis. In their paper, they
use a slightly higher spin value ($\chi=0.75$) and explore different spin
configurations, including antialigned and up-down directions with respect to
the orbital angular momentum. Their system has different thermodynamics than
ours, using an ideal-gas state equation with $\Gamma=4/3$ and no cooling.
Their focus is on the mass budget of the mini-disk (as in Bowen et al. (2019))
and the electromagnetic luminosity when spin is included. They report that
spinning black holes have more massive mini-disks and the electromagnetic
luminosity is higher, with quantitative measures similar to what we find in
this paper.

In our work, we analyze in great detail, for the first time, the accretion
mechanisms onto the mini-disk and their connection to the circumbinary disk.
We show that the BHs accretes in two different ways: through direct plunging
of the stream from the lump’s inner edge (that dominates the accretion), and
through ‘conventional’ stresses of the circular component orbiting the mini-
disk. This is qualitatively different than single BHs disks and a direct
consequenece of the short inflow time determined by $r_{\rm ISCO}/r_{\rm
trunc}$; for larger separations and higher spins, we expect mini-disks to
behave closer to conventional single BH disks . Our simulations also differ
significantly in the grid setup and initial data. We start our simulations
with an evolved circumbinary disk snapshot, taken from Noble et al. (2012),
which is already turbulent and presents a lump (starting the simulation from a
quasi-stationary torus, the lump appears after $\sim 50$ orbits at these
seperations, once the inner edge has settled). This is very important to
accurately describe the periodicities of the system given by the beat
frequency, which is set by the orbital motion of the lump. These quasi-
periodicities might be different if the thermodynamics change, e.g. if there
is no cooling, although currently there are no sufficiently long 3D GRMHD
simulations of circumbinary disks exploring this. Interestingly, we found that
the Poynting flux is also modulated by the beat frequency. For BBH approaching
merger, this constitutes a possible independent observable if this periodicity
is translated to jet emission. As expected, for spinning BHs, we also found
more powerful Poynting fluxes, in agreement with Paschalidis et al. (2021).

With our careful analysis of the accretion onto the mini-disks, we show that a
disk-like structure survives for longer as the binary shrinks when the black
holes have spin. Further explorations with higher spins will show how far
these structures survive very close to merger.

##  5\. Conclusions

We have performed a GRMHD accretion simulation of an equal-mass binary black
hole with aligned spins of $a=0.6~{}M_{\rm BH}$ approaching merger. We have
compared this simulation with a previous non-spinning simulation of the same
system, analyzing the main differences in mini-disk accretion and the
variabilities induced by the circumbinary disk accretion. Our main findings
can be summarized as follows:

  * •

Mini-disks in S06, where BHs have aligned spins $\chi=0.6$, are more massive
than in S0, where BHs have zero spins, by a factor of two. The mass and
accretion rate of mini-disks have quasi-periodicities determined by the beat
frequency in both simulations (see Section 3.2).

  * •

The material in the mini-disk region can be separated into two components of
relatively high and low angular momentum. The low angular momentum component
mostly plunges directly from the lump edge, forming a strong single-arm
stream. The (supra) Keplerian angular momentum component of the fluid is
determined by the size of the ISCO and the truncation radius. We have shown
that most of the mini-disk mass in S0 is sub-Keplerian while in S06 most of
the material follows closely the Keplerian value up to the end of the
evolution when it becomes comparable to the low angular momentum component.
For binary parameters as in this simulation, we have also predicted a critical
value of $\chi\sim 0.45$ for which most of the mass will have low angular
momentum relative to the ISCO.

  * •

The accretion rates at the horizon in S06 and S0 are very similar through the
evolution (see Fig. 5), as they are dominated by the plunging material from
the lump stream (see Figure 14).

  * •

In S06 a jet-like structure is formed self-consistently around each BH (see
Figures 2 and 18). Due to the black hole spin, there is a well-defined
Poynting flux (see Figure 18, left panel) with an efficiency of $\eta\sim 8\%$
(see Figure 17). On the other hand, in S0 the efficiency is closer to
$\eta\sim 4\%$, and its Poynting flux is more homogeneous in space (see Figure
18, right panel).

  * •

The time evolution of the Poynting flux is modulated by the quasi-periodicity
of the accretion, determined by the beat frequency. In the spinning case, the
fluxes in the comoving frame grow and start decreasing at
$7.5T_{\mathcal{B}}$, while in the non-spinning these remain fairly constant.

## Acknowledgments

We thank Eduardo Gutiérrez, Vassilios Mewes, Carlos Lousto, and Alex Dittman,
for useful discussions. L. C., F. L. A, M. C., acknowledge support from
AST-2009330, AST- 754 1028087, AST-1516150, PHY-1707946. 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. J.H.K. was
supported by AST-1028111, PHY-1707826, and AST-2009260 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).

Computational resources were provided by the Blue Waters sustained-petascale
computing NSF Projects No. OAC-1811228 and No. OAC-1516125”, and replace with
”Computational resources were provided by the TACC’s Frontera supercomputer
allocation No. PHY-20010 and AST-20021 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.

## References

  * Alam et al. (2020) Alam, M. F., et al. 2020, arXiv e-prints, arXiv:2005.06495
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., et al. 2017, arXiv e-prints, arXiv:1702.00786
  * Artymowicz & Lubow (1994) Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651
  * Babak et al. (2016) Babak, S., et al. 2016, MNRAS, 455, 1665
  * Baker et al. (2006) Baker, J. G., Centrella, J., Choi, D.-I., Koppitz, M., & van Meter, J. 2006, Phys. Rev. Lett., 96, 111102
  * Barnes & Hernquist (1992) Barnes, J. E., & Hernquist, L. 1992, ARA&A, 30, 705
  * Barnes & Hernquist (1996) Barnes, J. E., & Hernquist, L. 1996, The Astrophysical Journal, 471, 115
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, Monthly Notices of the Royal Astronomical Society, 179, 433
  * Bogdanović et al. (2009) Bogdanović, T., Eracleous, M., & Sigurdsson, S. 2009, ApJ, 697, 288
  * Bowen et al. (2017) Bowen, D. B., Campanelli, M., Krolik, J. H., Mewes, V., & Noble, S. C. 2017, ApJ, 838, 42
  * Bowen et al. (2018) Bowen, D. B., Mewes, V., Campanelli, M., Noble, S. C., Krolik, J. H., & Zilhão, M. 2018, The Astrophysical Journal Letters, 853, L17
  * Bowen et al. (2019) Bowen, D. B., Mewes, V., Noble, S. C., Avara, M., Campanelli, M., & Krolik, J. H. 2019, ApJ, 879, 76
  * Campanelli et al. (2006a) Campanelli, M., Lousto, C. O., Marronetti, P., & Zlochower, Y. 2006a, Phys. Rev. Lett., 96, 111101
  * Campanelli et al. (2006b) Campanelli, M., Lousto, C. O., & Zlochower, Y. 2006b, Phys. Rev. D, 74, 084023
  * Campanelli et al. (2006c) —. 2006c, Phys. Rev., D74, 041501
  * Campanelli et al. (2007) Campanelli, M., Lousto, C. O., Zlochower, Y., Krishnan, B., & Merritt, D. 2007, Phys. Rev., D75, 064030
  * Cattorini et al. (2021) Cattorini, F., Giacomazzo, B., Haardt, F., & Colpi, M. 2021, Physical Review D, 103, 103022
  * Chapon et al. (2013) Chapon, D., Mayer, L., & Teyssier, R. 2013, Monthly Notices of the Royal Astronomical Society, 429, 3114
  * Combi et al. (2021) Combi, L., Armengol, F. G. L., Campanelli, M., Ireland, B., Noble, S. C., Nakano, H., & Bowen, D. 2021, Phys. Rev. D, 104, 044041
  * d’Ascoli et al. (2018) d’Ascoli, S., Noble, S. C., Bowen, D. B., Campanelli, M., Krolik, J. H., & Mewes, V. 2018, ApJ, 865, 140
  * De Villiers & Hawley (2003) De Villiers, J.-P., & Hawley, J. F. 2003, ApJ, 592, 1060
  * De Villiers et al. (2005) De Villiers, J.-P., Hawley, J. F., Krolik, J. H., & Hirose, S. 2005, ApJ, 620, 878
  * Derdzinski et al. (2021) Derdzinski, A., D’Orazio, D., Duffell, P., Haiman, Z., & MacFadyen, A. 2021, MNRAS, 501, 3540
  * Derdzinski et al. (2019) Derdzinski, A. M., D’Orazio, D., Duffell, P., Haiman, Z., & MacFadyen, A. 2019, MNRAS, 486, 2754
  * Dittmann & Ryan (2021) Dittmann, A., & Ryan, G. 2021, arXiv preprint arXiv:2102.05684
  * Dittmann & Ryan (2022) Dittmann, A. J., & Ryan, G. 2022
  * D’Orazio et al. (2016) D’Orazio, D. J., Haiman, Z., Duffell, P., MacFadyen, A., & Farris, B. 2016, MNRAS, 459, 2379
  * D’Orazio et al. (2013) D’Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2997
  * D’Orazio et al. (2013) D’Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2997
  * D’Orazio et al. (2015) D’Orazio, D. J., Haiman, Z., & Schiminovich, D. 2015, Nature, 525, 351
  * Dotti et al. (2007) Dotti, M., Colpi, M., Haardt, F., & Mayer, L. 2007, MNRAS, 379, 956
  * Dotti et al. (2009a) Dotti, M., Montuori, C., Decarli, R., Volonteri, M., Colpi, M., & Haardt, F. 2009a, Monthly Notices of the Royal Astronomical Society: Letters, 398, L73
  * Dotti et al. (2009b) Dotti, M., Ruszkowski, M., Paredi, L., Colpi, M., Volonteri, M., & Haardt, F. 2009b, Monthly Notices of the Royal Astronomical Society, 396, 1640
  * Dotti et al. (2012) Dotti, M., Sesana, A., & Decarli, R. 2012, Advances in Astronomy, 2012, 940568
  * Duffell et al. (2020) Duffell, P. C., D’Orazio, D., Derdzinski, A., Haiman, Z., MacFadyen, A., Rosen, A. L., & Zrake, J. 2020, ApJ, 901, 25
  * East et al. (2012) East, W. E., Ramazanoğlu, F. M., & Pretorius, F. 2012, Physical Review D, 86, 104053
  * Escala et al. (2004) Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2004, ApJ, 607, 765
  * Escala et al. (2005) —. 2005, ApJ, 630, 152
  * Farris et al. (2014a) Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2014a, The Astrophysical Journal, 783, 134
  * Farris et al. (2014b) —. 2014b, Monthly Notices of the Royal Astronomical Society: Letters, 446, L36
  * Farris et al. (2015a) Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2015a, MNRAS, 447, L80
  * Farris et al. (2015b) —. 2015b, MNRAS, 446, L36
  * Farris et al. (2012) Farris, B. D., Gold, R., Paschalidis, V., Etienne, Z. B., & Shapiro, S. L. 2012, Phys. Rev. Lett., 109, 221102
  * Farris et al. (2012) Farris, B. D., Gold, R., Paschalidis, V., Etienne, Z. B., & Shapiro, S. L. 2012, Phys. Rev. Lett., 109, 221102
  * Farris et al. (2011) Farris, B. D., Liu, Y. T., & Shapiro, S. L. 2011, Phys. Rev. D, 84, 024024
  * Giacomazzo et al. (2012) Giacomazzo, B., Baker, J. G., Miller, M. C., Reynolds, C. S., & van Meter, J. R. 2012, The Astrophysical Journal Letters, 752, L15
  * Giacomazzo et al. (2012) Giacomazzo, B., Baker, J. G., Miller, M. C., Reynolds, C. S., & van Meter, J. R. 2012, ArXiv e-prints
  * Giacomazzo et al. (2012) Giacomazzo, B., Baker, J. G., Miller, M. C., Reynolds, C. S., & van Meter, J. R. 2012, Astrophys. J., 752, L15
  * Gold et al. (2013) Gold, R., Farris, B., Paschalidis, V., Etienne, Z., & Shapiro, S. 2013, in APS April Meeting Abstracts, H8002
  * Gold et al. (2014) Gold, R., Paschalidis, V., Etienne, Z. B., Shapiro, S. L., & Pfeiffer, H. P. 2014, Phys. Rev. D, 89, 064060
  * Gold et al. (2014) Gold, R., Paschalidis, V., Ruiz, M., Shapiro, S. L., Etienne, Z. B., & Pfeiffer, H. P. 2014, Phys. Rev. D, 90, 104030
  * Graham et al. (2015a) Graham, M. J., et al. 2015a, Monthly Notices of the Royal Astronomical Society, 453, 1562
  * Graham et al. (2015b) —. 2015b, Nature, 518, 74
  * Gutiérrez et al. (2021) Gutiérrez, E. M., Combi, L., Noble, S. C., Campanelli, M., Krolik, J. H., Armengol, F. G. L., & García, F. 2021, arXiv e-prints, arXiv: 2112.09773
  * Hawley & Krolik (2006) Hawley, J. F., & Krolik, J. H. 2006, ApJ, 641, 103
  * Healy & Lousto (2018) Healy, J., & Lousto, C. O. 2018, Physical Review D, 97, 084002
  * Hemberger et al. (2013) Hemberger, D. A., Lovelace, G., Loredo, T. J., Kidder, L. E., Scheel, M. A., Szilágyi, B., Taylor, N. W., & Teukolsky, S. A. 2013, Physical Review D, 88, 064014
  * Hirose et al. (2004) Hirose, S., Krolik, J. H., De Villiers, J.-P., & Hawley, J. F. 2004, ApJ, 606, 1083
  * Kelly et al. (2017) Kelly, B. J., Baker, J. G., Etienne, Z. B., Giacomazzo, B., & Schnittman, J. 2017, Phys. Rev. D, 96, 123003
  * Kelly et al. (2020) Kelly, B. J., Eteinne, Z. B., Golomb, J., Schnittman, J. D., Baker, J. G., Noble, S. C., & Ryan, G. 2020, Electromagnetic Emission from a Binary Black Hole Merger Remnant in Plasma: Field Alignment and Plasma Temperature
  * Kelly et al. (2021) Kelly, B. J., Etienne, Z. B., Golomb, J., Schnittman, J. D., Baker, J. G., Noble, S. C., & Ryan, G. 2021, Phys. Rev. D, 103, 063039
  * Kesden et al. (2014) Kesden, M., Gerosa, D., Berti, E., O’Shaughnessy, R., & Sperhake, U. 2014, in APS April Meeting Abstracts, Vol. 2014, C15–004
  * Khan et al. (2019) Khan, F. M., Mirza, M. A., & Holley-Bockelmann, K. 2019, Monthly Notices of the Royal Astronomical Society, 492, 256
  * Komissarov (2001) Komissarov, S. S. 2001, MNRAS, 326, L41
  * Krolik et al. (2005) Krolik, J. H., Hawley, J. F., & Hirose, S. 2005, ApJ, 622, 1008
  * Krolik et al. (2019) Krolik, J. H., Volonteri, M., Dubois, Y., & Devriendt, J. 2019, The Astrophysical Journal, 879, 110
  * Lin & Papaloizou (1979) Lin, D. N. C., & Papaloizou, J. 1979, Monthly Notices of the Royal Astronomical Society, 188, 191
  * Liu et al. (2003) Liu, F. K., Wu, X.-B., & Cao, S. L. 2003, Monthly Notices of the Royal Astronomical Society, 340, 411
  * Liu et al. (2019) Liu, T., et al. 2019, ApJ, 884, 36
  * Lopez Armengol et al. (2021) Lopez Armengol, F. G., et al. 2021, The Astrophysical Journal, 913, 16
  * Lousto & Healy (2015) Lousto, C. O., & Healy, J. 2015, Phys. Rev. Lett., 114, 141101
  * Lousto & Healy (2019) —. 2019, Phys. Rev. D, 99, 064023
  * Lousto et al. (2016a) Lousto, C. O., Healy, J., & Nakano, H. 2016a, Phys. Rev. D, 93, 044031
  * Lousto et al. (2016b) —. 2016b, Physical Review D, 93, 044031
  * Ma et al. (2021) Ma, S., Giesler, M., Scheel, M. A., & Varma, V. 2021, Physical Review D, 103, 084029
  * MacFadyen & Milosavljević (2008) MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83
  * Mayer (2013) Mayer, L. 2013, Classical and Quantum Gravity, 30, 244008
  * Mayer et al. (2007) Mayer, L., Kazantzidis, S., Madau, P., Colpi, M., Quinn, T., & Wadsley, J. 2007, Science, 316, 1874
  * McKinney & Gammie (2004) McKinney, J. C., & Gammie, C. F. 2004, The astrophysical journal, 611, 977
  * McKinney & Gammie (2004) McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977
  * Merritt (2004) Merritt, D. 2004, in Coevolution of Black Holes and Galaxies, ed. L. C. Ho, 263
  * Merritt (2006) Merritt, D. 2006, Reports on Progress in Physics, 69, 2513
  * Merritt & Ekers (2002) Merritt, D., & Ekers, R. D. 2002, Science, 297, 1310
  * Merritt & Milosavljević (2005) Merritt, D., & Milosavljević, M. 2005, Living Reviews in Relativity, 8, 8
  * Mihos & Hernquist (1996) Mihos, J. C., & Hernquist, L. 1996, ApJ, 464, 641
  * Miranda et al. (2017) Miranda, R., Muñoz, D. J., & Lai, D. 2017, MNRAS, 466, 1170
  * Mirza et al. (2017) Mirza, M. A., Tahir, A., Khan, F. M., Holley-Bockelmann, H., Baig, A. M., Berczik, P., & Chishtie, F. 2017, Monthly Notices of the Royal Astronomical Society, 470, 940
  * Moesta et al. (2012) Moesta, P., Alic, D., Rezzolla, L., Zanotti, O., & Palenzuela, C. 2012, The Astrophysical Journal Letters, 749, L32
  * Moody et al. (2019) Moody, M. S. L., Shi, J.-M., & Stone, J. M. 2019, ApJ, 875, 66
  * Mösta et al. (2019) Mösta, P., Taam, R. E., & Duffell, P. C. 2019, ApJ, 875, L21
  * Muñoz & Lai (2016) Muñoz, D. J., & Lai, D. 2016, ApJ, 827, 43
  * Muñoz et al. (2020) Muñoz, D. J., Lai, D., Kratter, K., & Miranda, R. 2020, ApJ, 889, 114
  * Muñoz & Lithwick (2020) Muñoz, D. J., & Lithwick, Y. 2020, ApJ, 905, 106
  * Muñoz et al. (2019) Muñoz, D. J., Miranda, R., & Lai, D. 2019, ApJ, 871, 84
  * Mundim et al. (2014) Mundim, B. C., Nakano, H., Yunes, N., Campanelli, M., Noble, S. C., & Zlochower, Y. 2014, Phys. Rev. D, 89, 084008
  * Nakamura et al. (2018) Nakamura, M., et al. 2018, The Astrophysical Journal, 868, 146
  * Neilsen et al. (2011) Neilsen, D., Lehner, L., Palenzuela, C., Hirschmann, E. W., Liebling, S. L., Motl, P. M., & Garrett, T. 2011, Proceedings of the National Academy of Sciences, 108, 12641
  * Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626
  * Noble et al. (2021) Noble, S. C., Krolik, J. H., Campanelli, M., Zlochower, Y., Mundim, B. C., Nakano, H., & Zilhão, M. 2021, The Astrophysical Journal, 922, 175
  * Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411
  * Noble et al. (2007) Noble, S. C., Leung, P. K., Gammie, C. F., & Book, L. G. 2007, Classical and Quantum Gravity, 24, 259
  * Noble et al. (2012) Noble, S. C., Mundim, B. C., Nakano, H., Krolik, J. H., Campanelli, M., Zlochower, Y., & Yunes, N. 2012, The Astrophysical Journal, 755, 51
  * Paczynski (1977) Paczynski, B. 1977, ApJ, 216, 822
  * Palenzuela et al. (2010a) Palenzuela, C., Garrett, T., Lehner, L., & Liebling, S. L. 2010a, Physical Review D, 82, 044045
  * Palenzuela et al. (2010b) Palenzuela, C., Lehner, L., & Liebling, S. L. 2010b, Science, 329, 927
  * Papaloizou & Pringle (1977) Papaloizou, J., & Pringle, J. E. 1977, MNRAS, 181, 441
  * Park et al. (2006) Park, S., Linsen, L., Kreylos, O., Owens, J., & Hamann, B. 2006, IEEE Transactions on Visualization and Computer Graphics, 12, 243
  * Paschalidis et al. (2021) Paschalidis, V., Bright, J., Ruiz, M., & Gold, R. 2021
  * Pretorius (2005) Pretorius, F. 2005, Phys. Rev. Lett., 95, 121101
  * Reardon et al. (2016) Reardon, D. J., et al. 2016, MNRAS, 455, 1751
  * Rodriguez et al. (2006) Rodriguez, C., Taylor, G. B., Zavala, R. T., Peck, A. B., Pollack, L. K., & Romani, R. W. 2006, The Astrophysical Journal, 646, 49
  * Roedig et al. (2014) Roedig, C., Krolik, J. H., & Miller, M. C. 2014, The Astrophysical Journal, 785, 115
  * Ryan & MacFadyen (2017) Ryan, G., & MacFadyen, A. 2017, ApJ, 835, 199
  * Saade et al. (2020) Saade, M. L., et al. 2020, arXiv e-prints, arXiv:2001.08870
  * Schoenmakers et al. (2000) Schoenmakers, A. P., de Bruyn, A. G., Rottgering, H. J. A., van der Laan, H., & Kaiser, C. R. 2000, Monthly Notices of the Royal Astronomical Society, 315, 371
  * Sesana & Khan (2015) Sesana, A., & Khan, F. M. 2015, MNRAS, 454, L66
  * Sesana et al. (2012) Sesana, A., Roedig, C., Reynolds, M. T., & Dotti, M. 2012, MNRAS, 420, 860
  * Shi & Krolik (2015) Shi, J.-M., & Krolik, J. H. 2015, ApJ, 807, 131
  * Shi et al. (2012) Shi, J.-M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, The Astrophysical Journal, 749, 118
  * Shi et al. (2012) Shi, J.-M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118
  * Springel et al. (2005) Springel, V., Matteo, T. D., & Hernquist, L. 2005, The Astrophysical Journal, 620, L79
  * Tchekhovskoy (2015) Tchekhovskoy, A. 2015, Launching of Active Galactic Nuclei Jets, ed. I. Contopoulos, D. Gabuzda, & N. Kylafis, Vol. 414, 45
  * Tiede et al. (2020a) Tiede, C., Zrake, J., MacFadyen, A., & Haiman, Z. 2020a, ApJ, 900, 43
  * Tiede et al. (2020b) —. 2020b, ApJ, 900, 43
  * Tóth (2000) Tóth, G. 2000, Journal of Computational Physics, 161, 605
  * Valtonen et al. (2006) Valtonen, M. J., et al. 2006, ApJ, 643, L9
  * Varma et al. (2018) Varma, V., Scheel, M. A., & Pfeiffer, H. P. 2018, Physical Review D, 98, 104011
  * Zilhão et al. (2015) Zilhão, M., Noble, S. C., Campanelli, M., & Zlochower, Y. 2015, Phys. Rev. D, 91, 024034
  * Zilhão & Noble (2014) Zilhão, M., & Noble, S. C. 2014, Classical and Quantum Gravity, 31, 065013
  * Zlochower et al. (2016) Zlochower, Y., Nakano, H., Mundim, B. C., Campanelli, M., Noble, S., & Zilhao, M. 2016, Physical Review D, 93, 124072
  * Zrake et al. (2020) Zrake, J., Tiede, C., MacFadyen, A., & Haiman, Z. 2020, arXiv e-prints, arXiv:2010.09707
  * Zurek & Benz (1986) Zurek, W. H., & Benz, W. 1986, ApJ, 308, 123