Minidisk Accretion onto Spinning Black Hole Binaries: Quasi-periodicities and Outflows

Assuming that the plasma does not influence the spacetime, we evolve the ideal GRMHD equations in the dynamical BBH 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:

$$\begin{eqnarray*}&&{{\rm{\nabla }}}_{a}(\rho {u}^{a})=0,\quad {{\rm{\nabla }}}_{a}{{T}^{a}}_{b}={{ \mathcal F }}_{b},\end{eqnarray*}$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{a}{\ }^{* }{F}^{{ab}}=0,\quad {u}_{a}{F}^{{ab}}=0,\end{eqnarray} \tag{ 1 }$$

where ρ is the rest-mass density, ${F}^{{ab}}/\sqrt{4\pi }$ is the Faraday tensor, ⁸ 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

$$\begin{eqnarray}&&{T}^{{ab}}:= (\rho h+{b}^{2}){u}^{a}{u}^{b}+\left(P,+,\displaystyle \frac{1}{2},{b}^{2}\right){g}^{{ab}}-{b}^{a}{b}^{b}\end{eqnarray} \tag{ 2 }$$

where h $:=$ (1 + + P/ρ) is the specific enthalpy, is the specific internal energy, P is the pressure, b^a $:=$ ^* F^ab u_b : is the four-vector magnetic field, and b² $:=$ b^a b_a is proportional to the magnetic pressure, p_m $:=$ b²/2. We follow Noble et al. (2012, 2009) and write these coupled equations of motion in manifest conservative form as

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{U}}\left({\boldsymbol{P}}\right)=-{\partial }_{i}{{\boldsymbol{F}}}^{i}\left({\boldsymbol{P}}\right)+{\boldsymbol{S}}\left({\boldsymbol{P}}\right),\end{eqnarray} \tag{ 3 }$$

where P are the primitive variables, U the conserved variables, F ⁱ the fluxes, and S the source terms. These are given explicitly as:

$$\begin{eqnarray}&&{\boldsymbol{P}}:=\,[\rho ,u,{\tilde{u}}^{j},{B}^{j}],\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\boldsymbol{U}}\left({\boldsymbol{P}}\right)=\sqrt{-g}\left[\rho {u}^{t},{{T}^{t}}_{t}+\rho {u}^{t},{{T}^{t}}_{j},{B}^{j}\right],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}^{i}\left({\boldsymbol{P}}\right)=\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],\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\boldsymbol{S}}\left({\boldsymbol{P}}\right)=\sqrt{-g}\left[0,{{T}^{\kappa }}_{\lambda }{{{\rm{\Gamma }}}^{\lambda }}_{t\kappa }-{{ \mathcal F }}_{t},{{T}^{\kappa }}_{\lambda }{{{\rm{\Gamma }}}^{\lambda }}_{j\kappa }-{{ \mathcal F }}_{j},0\right],\end{eqnarray} \tag{ 7 }$$

where g is the determinant of the metric, ${{{\rm{\Gamma }}}^{\lambda }}_{\alpha \beta }$ are the Christoffel symbols, u $:=$ ρ 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. ⁹ We close the system with a Γ-law equation of state, P = (Γ − 1)ρ , where we set Γ = 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 the 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}={{ \mathcal L }}_{c}{u}_{\beta }$, where ${ \mathcal L }$ 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}$:

$$\begin{eqnarray}&&{{ \mathcal L }}_{c}=\displaystyle \frac{\rho \epsilon }{{t}_{\mathrm{cool}}}{\left(\displaystyle \frac{{\rm{\Delta }}S}{{S}_{0}}+\left|\displaystyle \frac{{\rm{\Delta }}S}{{S}_{0}}\right|\right)}^{1/2},\end{eqnarray} \tag{ 8 }$$

where t_cool is the cooling timescale where the disk radiates away any local increase in entropy, ΔS $:=$ S − S₀, where S $:=$ p/ρ^Γ is the local entropy. Our target entropy, S₀ = 0.01, is the initial entropy of each accretion disk in the simulation. The timescale t_cool is determined by the local fluid orbital period, following the prescriptions in Bowen et al. (2017) and d'Ascoli et al. (2018).

We solve Equation (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).

The grid and boundary conditions used in this simulation are the same as in Bowen et al. (2018, 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}_{\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¹) 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²) have reflective, axisymmetric boundary conditions at the polar axis cutout and the azimuthal coordinates (x³) have periodic boundary conditions.

The resolution is given by 600 × 160 × 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 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 minidisks on the side farthest from the center of mass, see Figure 1.

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

$$\begin{eqnarray}&&{g}_{{ab}}={\eta }_{{ab}}+{M}_{1}{{ \mathcal H }}_{{ab}}^{1}({{\boldsymbol{x}}}_{1},{{\boldsymbol{v}}}_{1})+{M}_{2}{{ \mathcal H }}_{{ab}}^{2}({{\boldsymbol{x}}}_{2},{{\boldsymbol{v}}}_{2}),\end{eqnarray} \tag{ 9 }$$

where η_ab is the Cartesian Minkowski metric, ${{ \mathcal H }}_{{ab}}^{A}$ is the boosted black hole term for A = 1, 2, M_A is the mass, and { x _A (t), 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 BBH 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 BBH 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 relativity. 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 analyzed 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 10M, and that (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 because 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 valid only 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₁ = M₂ = M/2, with an initial separation of r₁₂(0) = 20M, 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 χ ≡ a/M = 0.6. Because 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. 2006c), and the orbital frequency increases (see Figure 2).

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 × 10⁴ M (∼80 orbits). At this time, the disk is in a turbulent state and the accretion into the cavity is dominated by the 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 of spin, so it is a good initial state for our spinning simulation. We interpolate these data onto our grid and we initialize two minitori 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). Further details of the simulations setup are outlined in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 minidisk properties such as fluxes, we shall 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 notice 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 ¹⁰ , which implements a fast Discrete Sibson interpolation (Park et al. 2006).

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

$$\begin{eqnarray}&&\langle { \mathcal Q }{\rangle }_{\sigma }:= \displaystyle \frac{\int {dA}\ \sigma \ { \mathcal Q }}{\int {dA}\ \sigma },\end{eqnarray} \tag{ 10 }$$

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

$$\begin{eqnarray} &&\langle\langle { \mathcal Q } \rangle\rangle := \displaystyle \frac{1}{{\rm{\Delta }}t}{\int }_{{t}_{f}}^{{t}_{i}}{dt}\langle { \mathcal Q }\rangle ,\end{eqnarray} \tag{ 11 }$$

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

In the steady state of an equal-mass SMBBH, a lump orbits the edge of the circumbinary disk at an average frequency $\langle {{\rm{\Omega }}}_{\mathrm{lump}}\rangle =0.28{{\rm{\Omega }}}_{{ \mathcal B }}$ (Noble et al. 2012, 2021; Shi et al. 2012; D'Orazio et al. 2013, 2016; Farris et al. 2015a, 2015b; Lopez Armengol et al. 2021) modulating the accretion into the cavity. In Figure 3, we show a global look of the system. When one of the BHs passes near the lump, it peels off part of the lump's inner edge, forming a stream that feeds the black hole with a beat frequency of ${{\rm{\Omega }}}_{\mathrm{beat}}:= {{\rm{\Omega }}}_{{ \mathcal B }}-\langle {{\rm{\Omega }}}_{\mathrm{lump}}\rangle \sim 0.72\ {{\rm{\Omega }}}_{{ \mathcal B }}$. This stream is almost ballistic, formed by fluid particles with relatively low angular momentum (Shi & Krolik 2015). As it approaches the black hole, this material can start orbiting the black hole, forming a minidisk. The maximum size of the minidisk is determined by the tidal truncation radius of the binary, or Hill's sphere. The residence time of matter in the minidisks 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 minidisks 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 ∼ 0.4 r₁₂(t) (Bowen et al. 2017), similar to the Newtonian value, estimated to be ∼0.3 r₁₂ (Paczynski 1977; Papaloizou & Pringle 1977; Artymowicz & Lubow 1994). Because 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 than 20M. The most relevant difference between our two simulations S06 and S0 is the location of the ISCO: ${r}_{\mathrm{ISCO}}(\chi =0.6)=2.82\ {M}_{{\mathrm{BH}}_{i}}$ for S06, and ${r}_{\mathrm{ISCO}}(\chi =0.0)=5.0\ {M}_{{\mathrm{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 minidisks. 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.

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

$$\begin{eqnarray}&&M:= {\int }_{{r}_{{ \mathcal H }}}^{{r}_{t}(t)}{dV}\ \rho {u}^{t},\end{eqnarray} \tag{ 12 }$$

where r_t (t) $:=$ 0.4 r₁₂(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 ∼3 orbits for S0 and ∼4 orbits for S06 (see Figure 4). Both simulations start with two quasi-equilibrated minitori 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). Because we are also using these 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 spinning simulation.

Zoom In Zoom Out Reset image size

We find that although the minidisks of S06, like those in S0, go through a filling-depletion cycle, the minidisks 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 minidisk 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 ∼8 orbits (upper panel, color blue, of Figure 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 (see below).

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

$$\begin{eqnarray}&&\dot{M}:= {\oint }_{{r}_{{ \mathcal H }}}d\bar{A}\ {u}^{\bar{r}}\rho .\end{eqnarray} \tag{ 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}}_{\mathrm{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 minidisk mass is, to a first approximation, a smoothed version of the accretion rate's time dependence, e.g., see Figure 6 for BH₁ in S06.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The differences in the masses per cycle are closely related to the inflow time of particles in the minidisk. 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

$$\begin{eqnarray}&&{t}_{\mathrm{inflow}}^{-1}:= \displaystyle \frac{1}{\bar{r}}\langle \langle {V}^{\bar{r}}{\rangle }_{\rho }\rangle ,\end{eqnarray} \tag{ 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 BH₁, 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.

Zoom In Zoom Out Reset image size

As the binary shrinks, the inflow time at the truncation radius diminishes from ∼$0.41\ {T}_{{ \mathcal B }}$ to ∼$0.31\ {T}_{{ \mathcal B }}$ for the spinning simulation. The average inflow time at the truncation radius is well below the beat period, ${T}_{\mathrm{beat}}=1.3{T}_{{ \mathcal B }}$, on which the minidisk 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 minidisks 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 Section 3.3).

The quasi-periodic behavior of the system can be described through a power density spectrum (PSD) of the minidisk'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₁ and M₂ taken individually peak at the beat frequency in both simulations. The PSD of the total mass of the minidisks, M₁ + M₂, has a peak at 2Ω_beat in S0, while the latter is severely damped in S06. Indeed, if the individual masses vary with a characteristic frequency Ω_beat, and these are out of phase with the same amplitude, we expect their sum to vary with 2Ω_beat. In S06, however, the inflow time of the minidisks is larger and the depletion period of a minidisk briefly coexists with the filling period of the other minidisk, reducing the variability of the total mass. On the other hand, the beat frequency is slightly higher for S06 (${{\rm{\Omega }}}_{\mathrm{beat}}=0.71{{\rm{\Omega }}}_{{ \mathcal B }}$) than S0 (${{\rm{\Omega }}}_{\mathrm{beat}}=0.68{{\rm{\Omega }}}_{{ \mathcal B }}$) as the orbital frequency of the spinning BHs is higher. Further, because we evolved the binary for longer, in S06 we find a more prominent amplitude at ∼$0.20{{\rm{\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 that this frequency is closely related but different from the orbital frequency of the lump at $0.28{{\rm{\Omega }}}_{{ \mathcal B }}$, which could be related to the orbital frequency increasing during the inspiral (notice we measure the PSD at fixed frequencies).

Zoom In Zoom Out Reset image size

Finally, because we use a spherical grid with a central cutout, we cannot analyze the effects of the sloshing of matter between minidisks (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 minidisks in spinning and non-spinning BBH.

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

The surface density is defined as ${\rm{\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 seventh orbit. In this plot, the minidisk around BH₁ (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 BH₁ in both simulations, but there is much more of it in S06. On the other hand, we observe that BH₂ (left side), in its low state, has a noticeable disk structure in S06, while the material is already depleted for S0.

Zoom In Zoom Out Reset image size

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 minidisk with respect to the orbital motion. We define

$$\begin{eqnarray}&&\langle {\rm{\Sigma }}(r,t)\rangle := \displaystyle \frac{{\int }_{{\rm{\Delta }}\bar{\phi }}d\bar{l}\ {\rm{\Sigma }}}{{\int }_{{\rm{\Delta }}\bar{\phi }}d\bar{l}},\end{eqnarray} \tag{ 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$Σ(r)$\rangle\rangle$ for ${\rm{\Delta }}{\bar{\phi }}_{1}=(\pi /4,3\pi /4)$ and ${\rm{\Delta }}{\bar{\phi }}_{2}=(5\pi /4,7\pi /4)$, averaging in time over the high and low state of the minidisk separately. In the high state, the minidisk accumulates more material at the front, while the back of the minidisk 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 minidisk, 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 ∼2 in both states.

Zoom In Zoom Out Reset image size

Another important property of minidisks in relativistic binaries is the presence of nontrivial azimuthal density modes. When the accreting stream of the lump impacts the minidisk, it generates a pressure wave, forming strong spiral shocks (Bowen et al. 2018). This induces an m = 1 density mode in the minidisk 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 minidisk rest-mass density in azimuthal Fourier modes (Zurek & Benz 1986), $\rho (\bar{\phi })\equiv {\sum }_{m}{D}_{m}\exp (-{im}\bar{\phi })$, where

$$\begin{eqnarray}&&{D}_{m}:= {\int }_{{r}_{{ \mathcal H }}}^{{r}_{t}(t)}d\bar{V}\rho \exp (-{im}\bar{\phi }).\end{eqnarray} \tag{ 16 }$$

Let us compare these modes in S0 and S06 for BH₁. From Figure 11, we observe that both simulations share common features. In both simulations, the minidisks 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 minidisks increase their mass and the stream is accreted into 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 minidisks. This behavior could indicate that, as the minidisks 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 minidisk 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 $10{T}_{{ \mathcal B }}$ in S06 (see Figure 13 and discussion below). This is another consequence of the disk-like structure of the minidisks surviving for longer time in the spinning case.

Zoom In Zoom Out Reset image size

Our analysis so far indicates that a considerable amount of the matter in the minidisk region is plunging directly from the lump to the black hole. To further analyze the orbital motion of the fluid in the minidisk, we compute the density-weighted specific angular momentum in the BH frame $\langle \bar{{\ell }}{\rangle }_{\rho }$, where $\bar{{\ell }}:= -{u}_{\bar{\phi }}/{u}_{\bar{t}}$. 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 (Noble et al. 2012; Shi 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 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 to the Keplerian value. For S0, the behavior is always sub-Keplerian on average. If the angular momentum distribution of the circumbinary streams is independent of spin for a fixed binary separation and mass-ratio, the angular momentum with which the streams arrive at the minidisk is greater than the ISCO angular momentum for spin χ > 0.45. This estimate could serve as a crude criterion for determining whether minidisks form in relativistic binaries.

Zoom In Zoom Out Reset image size

We can also use the specific angular momentum to distinguish the material in the minidisk 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 }}\lt {\bar{{\ell }}}_{K}$ and $\bar{{\ell }}\geqslant {\bar{{\ell }}}_{K}$ separately. In Figure 13 we plot the evolution of the sub-Keplerian and super-Keplerian mass components for BH₁ 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.

Zoom In Zoom Out Reset image size

Although a fair amount of the mass in the minidisk 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.

Zoom In Zoom Out Reset image size

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 minidisks 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.

Zoom In Zoom Out Reset image size

When the minidisk is in its high state, the spiral shocks heat up the minidisks 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 minidisk. At the peak of the accretion cycle, the aspect ratio rises to h/r ∼ 2, but the gas cools before the next accretion event.

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 minidisk evaluated in the BH frame is:

$$\begin{eqnarray}&&{L}_{\mathrm{EM}}(t,\bar{r})={\oint }_{\bar{r}}d\bar{A}\ {{ \mathcal S }}^{\bar{r}},\end{eqnarray} \tag{ 17 }$$

where the Poynting flux is ${{ \mathcal S }}^{\bar{i}}:= ({T}_{\mathrm{EM}}{)}_{\ \ \bar{t}}^{\bar{i}}$. In Figure 16, we plot the emission measure (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 both S06 and S0 the EM luminosities are variable, and in both a Fourier power spectrum reveals a periodic modulation at the orbital frequency of the circumbinary disk's inner edge, i.e., the "lump" frequency, see Figure 17. However, in S06 there is an additional modulation of similar amplitude at twice the beat frequency, proving it is tied closely to accretion. In S0, there are no clear long-term trends, whereas in S06, there is a secular growth of L_EM until $7.5{T}_{{ \mathcal B }}$, when it starts declining.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Physical and Grid Parameters of both Non-spinning and Spinning Simulations

Simulation	S06	S0
Spin parameter [χ]	0.6	0.0
BH1,2 mass [MBH]	0.5
Mass-ratio [q]	1
Final time [tf ]	8000M	6000M
Final separation [r12(tf )]	16.6M	17.8M
# Orbits	15	12.5
Init. separation [r12(0)]	20M
Init. total minidisk mass [M0]	20
Average orbital period $[{T}_{{ \mathcal B }}]$	530M
Lump orbital frequency [Ωlump]	$0.28\ {{\rm{\Omega }}}_{{ \mathcal B }}$
ISCO radius [rISCO]	2.82 MBH	5.0 MBH
Truncation radius [rtrunc]	0.4 r12
Grid [(x1 × x2 × x3)]		(600, 160, 640)
Physical Size [(${r}_{\min },{r}_{\max }$)]		(2M, 260M)

Download table as:  ASCIITypeset image

The instantaneous efficiency $\eta ={L}_{\mathrm{EM}}(t)/\dot{M}(t)$ increases during the first several orbits, saturating at ≈0.05 in the case of S06, but an order of magnitude lower for S0, as can be seen in Figure 18. 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 η(t) until $10{T}_{{ \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 (Palenzuela et al. 2010b; Neilsen et al. 2011). 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 = 100M surrounding the binary system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 19, we plot this quantity for both S06 and S0 as a function of retarded time t − r/〈v〉, where 〈v〉 is the mean velocity of the outflow. We also plot the sum of the Poynting fluxes around each minidisk, measured in the BH comoving frame. We notice that the fluxes in S0 are on average five times larger in the center of mass frame compared with the BH frame, suggesting that there is a kinematic contribution to L_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. ¹¹ During our simulations, the speed of the black holes remains fairly constant at v ∼ 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), and Kelly et al. (2021).

To further study the spatial distribution of the electromagnetic flux, in Figure 20, 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 4000M. The dotted-dashed white lines are defined by b²/ρ = 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.

In Figure 21, we plot ${ \mathcal P }$, averaged over half an orbit in the corotating frame of the binary, for a sphere of radius r = 60M 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 compared with the angular size of the jets, we are still able to pick out important features.

Zoom In Zoom Out Reset image size

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):

$$\begin{eqnarray}&&{L}_{{\rm{H}}}(t,r):= \left({\oint }_{r}{dA}\ ({T}_{{\rm{H}}}{)}_{t}^{r}\right)-{\dot{M}}_{\mathrm{jet}}.\end{eqnarray} \tag{ 18 }$$

In Figure 22, 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 (Equation (18)), we include only fluid elements that are both unbound according to the Bernoulli criterion ( − hu_t > 1) and moving outward (u^r > 0) (see 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 17 shows the Fourier power spectrum for spinning and non-spinning cases. We observe in Figure 22 a secular growth of the hydro fluxes in S06, while in S0 they remain rather constant, with an average efficiency of ∼2%. Also like the EM case, the hydro energy flux is considerably greater in S06 than S0, but by a factor ∼5. Such a contrast resembles the differences between the hydro efficiencies measured in spinning and non-spinning single BH simulations (Hawley & Krolik 2006).

Zoom In Zoom Out Reset image size

We caution, however, that the luminosities measured at 100M 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 100M 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.