Electromagnetic Signatures from Supermassive Binary Black Holes Approaching Merger

We use data from two GRMHD simulations: one with aligned spins, hereafter S06, performed in Combi et al. (2021b), and the other with zero spins, hereafter S0, performed in Bowen et al. (2018, 2019). In both simulations, the GRMHD equations for the plasma are evolved using the finite-volume code Harm3d with the same initial data and numerical setup, with the only difference being the metric spacetime.

In S0, the plasma evolves on top of a binary black hole spacetime represented by a semi-analytical approximate metric constructed using a matching technique (Mundim et al. 2014), and in which the black holes have zero spins. In S06, the spacetime is represented by a different approximate metric, built through a superposition, in which black holes have aligned spins of ${\chi }_{i}\equiv {J}_{i}/{m}_{i}^{2}=0.6$, where J_i and m_i are the angular momentum and mass of the ith black hole. Combi et al. (2021a) demonstrated that these two spacetime metrics—in both of which the binary separation is a = 20M in the initial state—are physically equivalent, and the only meaningful difference in the simulations is the spin of the black holes.

The initial data for the CBD are taken from a snapshot in Noble et al. (2012) in which the system evolved for 50,000M, reaching a relaxed state in which an m = 1 overdensity, known as the lump, orbits around the inner edge of the circumbinary cavity. Computing the time-average of the accretion rate, we found that the CBD is in excellent inflow equilibrium up to a distance r ∼ 50M, a little bit outside its inner edge. From there to ≃100M = 5a, the mean accretion rate increases slowly, rising by ≃50%. For this reason, the lowest-frequency portions of the CBD disk spectrum contribution are somewhat stronger than they would be if we had a genuine inflow equilibrium. However, because these low frequencies are not in a part of the spectrum we use (see Section 3), this inconsistency does not affect any of our conclusions.

The black holes are separated by a distance of r₁₂ = 20M, where M = m₁ + m₂ is the total mass of the system. This data is interpolated into a new grid that includes the black holes and two mini-disks on quasi-equilibrium and then cleaned of magnetic divergences; see Bowen et al. (2018) for more details. The simulations use outflow boundary conditions on the radial boundaries, reflective, axisymmetric boundary conditions at the polar axis cutout, and periodic boundary conditions on the azimuthal coordinates, and follows the black hole inspiral using post-Newtonian trajectories.

After a transient of ∼3 orbits, the mini-disks settle into a slowly-evolving limit cycle: they fill and drain, but with slowly-changing maximum and minimum masses. This cycle is quasi-steady after the transient, as can be seen in Figure 3 of Bowen et al. (2019) and Figure 4 of Combi et al. (2021b). In Figure 1 we show an equatorial view of the density and cooling function of the GRMHD simulation. Simulation S06 (S0) lasts for a total of 15 (12) orbits, ending with the black holes at a separation of ∼16.7M (17.3M).

Zoom In Zoom Out Reset image size

In Harm3d, the energy-momentum conservation equations have a loss-term that radiates away orbital energy dissipated into heat:

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{T}^{\mu \nu }={{ \mathcal L }}_{c}{u}^{\nu },\end{eqnarray} \tag{ 1 }$$

The cooling function is

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\rm{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),\end{eqnarray} \tag{ 2 }$$

where ρ is the rest-mass density, is the specific internal energy, ΔS ≡ S − S₀, S is the specific entropy, S₀ = 0.01 is the initial specific entropy, and the cooling timescale t_cool is defined as

$$\begin{eqnarray}{t}_{\mathrm{cool}}(r)=\left\{\begin{array}{ll}2\pi {\left(r+M\right)}^{3/2}/\sqrt{M}, & \mathrm{if}\ r\geqslant 1.5a,\\ 2\pi {r}_{\mathrm{BL},i}^{3/2}/\sqrt{{m}_{i}}, & \mathrm{if}\ {r}_{i}\leqslant 0.45a,\\ 2\pi {\left(1.5a+M\right)}^{3/2}/\sqrt{M}, & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

In the equation above, r_BL,i is the Boyer–Lindquist radial coordinate in the co-moving frame of the ith black hole, and m_i is the black hole's mass. Fixing the target entropy S₀, we control the temperature of the fluid and thus its height-radius aspect to h/r ≈ 0.1.

The cooling term assures that all of the energy dissipated is radiated, and the distribution in space and time of this radiation is very close to the dissipation distribution. In codes modeling ideal MHD, dissipation is a numerical grid-scale effect, but because physical dissipation tends to be enhanced by sharp gradients, whether of velocity or magnetic field, grid-scale dissipation mimics physical dissipation.

We use the code Bothros (Noble et al. 2007; d'Ascoli et al. 2018) to calculate the EM emission from our GRMHD simulations. Our methodology follows the one used in d'Ascoli et al. (2018). We use a camera-to-source approach, in which photons are launched from the camera going backward in time toward the source. We also perform the integration using the fast-light approximation, which assumes that photons move much faster than the plasma and the spacetime itself. ⁶

In Figure 2, we show the paths of a sample of geodesics launched from a camera located on the equatorial plane of the binary system at (x = 0, y = −1000 M). We distinguish the geodesics that arise from the event horizon of the holes with dark solid lines; at infinity, these define the black hole shadows. ⁷

Zoom In Zoom Out Reset image size

After calculating the full sample of geodesic paths for a given time (once per pixel), we integrate the radiative transfer equation along the geodesics to produce maps of specific intensity ${I}_{{\nu }_{\infty }}$ at the camera for a given observed frequency ν_∞. From these intensity maps, we then calculate flux spectra and light curves. To faithfully recover the flux spectrum for each snapshot, we use a resolution of 800 × 800 pixels (∼7 pixels/M) for the region r < 60M and 500 × 500 pixels (∼2 pixels/M) for the region 60M < r < 150M. This resolution is sufficient to compute the flux to a precision of ∼1%.

To obtain ${I}_{{\nu }_{\infty }}$ at each pixel, we first calculate the optical depth along the geodesic and determine whether the disk's photosphere, defined as the loci of points where τ = 1, is reached or not. If τ < 1 along the geodesic, then we start the radiative transfer integration from the last point in the geodesic (either at the end of the simulation domain or at the event horizon of one of the black holes) with the initial condition I₀ = 0. On the other hand, if the photosphere is reached, we start the integration from the photosphere with the initial condition ${I}_{0}={I}_{\mathrm{photosphere}}={B}_{\nu }({T}_{\mathrm{eff}})$. Here, B_ν (T) is the Planck function, and T_eff is the photosphere's effective temperature at the point where the given geodesic intersected it. Placing this initial condition is equivalent to assuming that the photosphere radiates a blackbody spectrum at the local value of the effective temperature.

The effective temperature is calculated using the Stefan–Boltzmann law from the radiative cooling flux ${ \mathcal F }$:

$$\begin{eqnarray}&&{T}_{\mathrm{eff}}={\left({ \mathcal F }/\sigma \right)}^{1/4},\end{eqnarray} \tag{ 4 }$$

where σ is the Stefan–Boltzmann constant. In turn, the radiative flux is calculated by integrating vertically the cooling function inside the photosphere:

$$\begin{eqnarray}&&{ \mathcal F }=\displaystyle \frac{1}{2}{\int }_{\tau \gt 1}{{ \mathcal L }}_{{\rm{c}}}(z){dz},\end{eqnarray} \tag{ 5 }$$

where the factor 1/2 comes from the fact that the disk has two surfaces through which to cool.

Regardless of the initial point, the radiative transfer equation is solved only in the flow's optically thin region, referred to as the corona. To perform the integration, we must assume an emission and absorption model. Electrons in AGN coronae reach relativistic temperatures and cool primarily through inverse Compton scattering of low-energy photons from the optically thick disk. This is a nonlocal process in which a photon may suffer many scatterings before escaping the system. It produces a nonthermal spectrum (power law) with a high energy cutoff, regardless of whether particles follow a thermal energy distribution, which is usually the case. A detailed treatment of this process is very complex and computationally costly (see, however, Kinch et al. 2020, 2021). For simplicity, we follow d'Ascoli et al. (2018) and model the multiple nonlocal Compton scatterings using a local effective emissivity coefficient of the form

$$\begin{eqnarray}&&{j}_{\nu }=\displaystyle \frac{{{ \mathcal L }}_{{\rm{c}}}}{4{\pi }^{3/2}}\displaystyle \frac{h}{{k}_{{\rm{B}}}T}{\left(\displaystyle \frac{h\nu }{{k}_{{\rm{B}}}T}\right)}^{-1/2}{e}^{-\tfrac{h\nu }{{k}_{{\rm{B}}}T}}\end{eqnarray} \tag{ 6 }$$

that mimics the actual output from the Comptonization. The normalization is such that $4\pi \int {j}_{\nu }d\nu ={{ \mathcal L }}_{{\rm{c}}}$, namely that the bolometric emissivity matches the cooling function at every point, and we fix the dimensionless temperature and spectral index of the corona to typical values observed in luminous single AGNs: Θ = k_B T/m_e c² = 0.2 and p = 0.5.

At the temperatures and densities of interest in our scenario, the dominant source of opacity is electron scattering. We assume a gray (frequency-independent) Thomson opacity,

$$\begin{eqnarray}&&{\alpha }_{\nu }=\mathrm{const}.\,=\,{\sigma }_{{\rm{T}}}{n}_{{\rm{e}}},\end{eqnarray} \tag{ 7 }$$

where n_e = ρ/m_H is the electron number density, and σ_T is the Thomson cross section.

Finally, after calculating the specific intensity at every pixel for a given time t, we integrate over all pixels to obtain the flux:

$$\begin{eqnarray}&&{F}_{{\nu }_{\infty }}=\int {I}_{{\nu }_{\infty }}(\xi ,\eta )d\xi d\eta ,\end{eqnarray} \tag{ 8 }$$

where (ξ, η) is the pair of angular coordinates at the camera, and we have assumed that this is located sufficiently far from the system so that the rays arrive approximately parallel to each other. ⁸ Then, the spectral luminosity is simply${L}_{{\nu }_{\infty }}=4\pi {r}_{\mathrm{cam}}^{2}{F}_{{\nu }_{\infty }}$.

Our GRMHD simulations ignore the fluid self-gravity, which is not important in this scenario. Thus, there is scale freedom for both the total mass of the binary system and the physical mass density scale of the gas. The latter can be set indirectly using a more easily interpretable quantity such as the accretion rate. This means that for each of the two sets of simulation data (spinning and nonspinning), we can investigate various scenarios in which either the total mass M, the accretion rate $\dot{M}$, or both change.

The accretion rate is calculated in the simulation in code units (CU) as

$$\begin{eqnarray}&&\dot{{ \mathcal M }}(r,t)=-{\int }_{{{ \mathcal S }}^{2}(r)}\rho {u}^{r}\sqrt{-g}d\theta d\phi ,\end{eqnarray} \tag{ 9 }$$

where ρ is the fluid mass density, u^r is the radial component of the fluid 4-velocity, g is the determinant of the metric, and ${{ \mathcal S }}^{2}(r)$ is a spherical surface of fixed radius r at time t in the harmonic (center of mass) coordinate system. We average $\dot{{ \mathcal M }}$ over the range 2a < r < 4a for the initial snapshot; this calculation gives a code-unit value of $\langle \dot{{ \mathcal M }}\rangle \approx 0.015$. Comparing this to a given assumed physical accretion rate $\dot{M}$, we obtain the mass density scale. ⁹

In both S06 and S0 simulations, the flow cools efficiently enough to maintain an average height-to-radius ratio of h/r ∼ 0.1. This aspect ratio is what is usually expected in quasars, where typical estimates of the accretion rate are ∼10^−1±1 × the Eddington rate. For rough consistency, we compute the radiative output assuming $\dot{M}=0.25{\dot{M}}_{\mathrm{Edd}}$, where ${\dot{M}}_{\mathrm{Edd}}\equiv 1.39\times {10}^{18}\left(M/{M}_{\odot }\right)\,{\rm{g}}\,{{\rm{s}}}^{-1}$ is the Eddington accretion rate with an assumed efficiency of 10%. We also set the total mass of the binary to M = 10⁶ M_⊙ as a fiducial value. In Section 4, we explore other values for the black hole mass.

We choose a face-on view (i = 0°) to investigate general features in the emission of these systems that are independent of inclination. ¹⁰ At higher viewing angles, the detected radiation would be strongly modulated by relativistic effects such as gravitational lensing and Doppler boosting. These modulations may be a useful tool for identifying SMBBHs through timing analysis, and they will be fully investigated in future work. Note, however, that the emission in the face-on case is still affected by gravitational redshift and transverse Doppler effects.

The spectrum of the system can be analyzed in terms of its three main contributors: the CBD, the streams, and the mini-disks. In each case, the emission is the sum of two components: a multitemperature blackbody, effectively emitted by the accretion flow's photosphere, and a Comptonized power law, extending up to hard X-rays, that reflects the energy dissipation in the corona. We assumed that this emission is produced by the Compton up-scattering of low-energy photons in the optically thin hot corona. Figure 3 shows an example of a spectrum for S06 at a specific time, t = 3290M. The CBD is defined as the region r > 2r₁₂, the streams as the region r₁₂ < r < 2r₁₂, and the mini-disks as the regions r_i < 0.45r₁₂, where r_i denotes the radial coordinate in the reference frame of the ith black hole. The remaining region is almost devoid of matter and has a negligible contribution; however, it is included in the total emission. The CBD spectrum is depicted by the blue solid line. Its thermal component peaks in the UV band because its inner boundary reaches an effective temperature of ∼0.9 × 10⁵K; this component dominates the bolometric luminosity of the system. The plasma in the mini-disks reaches higher effective temperatures (≳3.2 × 10⁵ K), and its spectrum peaks in the far-UV. A typical spectrum from S0 shows the same features described above.

Zoom In Zoom Out Reset image size

In accretion disks onto single black holes, one would expect the density and temperature to follow a continuous trend down to ≲1.5r_ISCO. On the contrary, the region between the CBD and the mini-disks in SMBBHs is usually rather empty and cold. This could result in a notch in the spectrum between the two thermal peaks (Roedig et al. 2014). However, this effect is missing in our spectrum because the mini-disks' luminosity is too low to create the expected rise at short frequencies (as we will show, this condition is specific to very close binaries). An additional minor effect is that the streams' emission lies in the middle of the two thermal peaks, though it always lies below the sum of the CBD's and mini-disk's emission.

In contrast to typical single black hole thin disks, the global peak in our spectrum does not correspond to the highest temperatures, reached here in the mini-disks. The mini-disk's low luminosity is caused by a combination of two effects. On one hand, the accretion rate is lower than in the CBD; on the other hand, a portion of the plasma in the mini-disks has low angular momentum; not requiring much angular momentum loss to fall into the black hole, less dissipation takes place within it. Indeed, as shown in Combi et al. (2021b), while some of the material in the mini-disk can orbit the black hole, another component plunges directly into the hole (see also Beloborodov & Illarionov 2001 for a similar situation potentially occurring in X-ray black hole binaries). In Section 4.2, we investigate the relative importance of these two effects in depth.

To compare the emission of S06 and S0 simulations, we first show in Figure 4 brightness maps at three different frequencies, time-averaged over the fifth binary orbit and displayed in the binary's co-rotating frame. The choice of the fifth orbit is arbitrary, but it secures that both simulations have passed the initial transient. The upper panel corresponds to S06 and the lower panel to S0. Both simulations share several emission characteristics. The UV map (first panel, ν = 6.6 × 10¹⁵ Hz or λ = 22.0 nm) shows a surface brightness that is nearly constant from the CBD down to the mini-disks, except for the low-density cavity. As a result, the CBD dominates the total flux due to its much larger area (see Figure 3). The far-UV emission (second panel, ν = 2.8 × 10¹⁶ Hz or λ = 99.3 nm) reveals higher temperature regions, and mini-disks and streams now dominate. The thermal emission from the CBD is lower, indicating that this frequency has passed its thermal maximum. In X-rays (third panel, ν = 10¹⁸ Hz or h ν = 4.1 keV), the emission is from the optically thin regions and is dominated by the mini-disks, showing that a significant fraction of their emission occurs in the corona rather than in the optically thick disk.

Zoom In Zoom Out Reset image size

The main difference between the S06 and S0 panels is in the mini-disks' thermal emission. While the X-ray map is fairly similar in both scenarios, the first two panels show that the thermal emission from the mini-disks and streams is ∼3 times higher in S06. This happens because the smaller value of l_crit when the black holes spin channels a larger fraction of the accretion rate into orbits within the mini-disks rather than plunging orbits. In order to move inward, the matter held within the mini-disks must lose angular momentum, and the processes that transfer angular momentum are accompanied by dissipation whether they are magnetic stresses associated with MHD turbulence or shocks. Furthermore, the streams in S06 are brighter than in S0 both where they strike the CBD and as they fall toward the mini-disks.

To investigate the radial distribution of luminosity, we integrate the intensity maps in the azimuthal coordinate of the image. Figure 5 depicts $d(\nu {L}_{\nu })/d\mathrm{log}r$ in the radial range 4–60 M for the three frequencies considered in Figure 4 (here r is the radius from the center of mass). The upper panel shows the luminosity in the UV band, where the CBD dominates the emission. The mini-disks and streams in S06 are brighter than in S0, but the difference fades and becomes insignificant in the CBD region. In the far-UV band (middle panel), the mini-disks surrounding the spinning black holes are a factor of two to five times brighter than in the nonspinning simulation. The relative variance in the emission is also higher than at lower frequencies, particularly at the boundary between the mini-disk and stream regions. This could be a result of our different prescriptions for the cooling function in the mini-disks and CBD. On the contrary, the optically thin emission (lower panel) from the mini-disks is nearly indistinguishable in the two cases. The optical depth to Compton scattering in the corona is always unity, and hence the coronal surface density is 1/κ_T, where κ_T is the Thomson opacity. Since the area of the mini-disks changes very little with spin, so does the total mass in the corona. Therefore, the near-constancy of the coronal luminosity with respect to spin suggests that the dissipation rate per unit mass in the corona is also spin-insensitive. For the same reason, the coronal luminosity should also not be overly affected by a change in accretion rate provided the mini-disks' surface density remains >1/κ_T over most of their area.

Zoom In Zoom Out Reset image size

Due to the complex dynamics associated with orbital motion, accretion flows onto SMBBHs are intrinsically variable. It is worthwhile to investigate how much of this variability is inherited by the EM emission. In Figure 6, we show the time-average of the spectral energy distributions (SEDs) over each binary orbit. The first orbits correspond to the transient phase, during which the mini-disk emission unphysically increases rapidly and then gradually decreases and stabilizes, becoming quasi-steady after the ∼4th orbit. After the ∼8th orbit, the mini-disks in the S0 simulation barely contribute to the thermal spectrum. In S06, on the other hand, their contribution in the far-UV band is nonnegligible until the last (15th) orbit, owing to the persistence of a disk-like structure (Combi et al. 2021b).

Zoom In Zoom Out Reset image size

In Figure 7, we compare the SEDs of the two simulations for the averaged fifth and 10th orbits. Though the bolometric luminosity, which is dominated by CBD emission, is comparable for both simulations, the mini-disk's thermal peak (shown in thicker lines) is ∼3–5 times brighter for the spinning black holes. Once again, this difference is explained by the fact that, in the spinning case, the Keplerian angular momentum at the ISCO is smaller, so that a greater amount of mass must lose more angular momentum, dissipating more energy in the process.

Zoom In Zoom Out Reset image size

We then use the spectral flux at each time to calculate light curves for S06 and S0 simulations. The light curves at different frequencies illustrate the variability in the properties of the various contributors to the emission. The lump modulates the accretion rate and mass in the mini-disks, which is inherited by the luminosity (see Figure 8). Nonetheless, because we would, in principle, observe the joint emission from the two mini-disks, we also look at the variability in the total emission from the system. Figure 9 displays light curves in the UV band (upper panel), where the CBD emission peaks, far-UV (middle panel), where the mini-disk emission peaks, and soft X-rays (lower panel), where the optically thin coronal emission dominates (see also Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To further investigate the putative periodicities, we calculate the power spectral density (PSD) of the three light curves for each simulation, using a time sampling of 10M. The spinning case is depicted in the upper panel of Figure 10, where the PSD has been normalized to its maximum value, and the frequency is shown in units of the mean orbital frequency. The mean period for S06 is 〈P〉 ≈ 505M, and for S0 it is 〈P〉 ≈ 530M. The three curves peak at $f\sim 0.2{f}_{{ \mathcal B }}$, which approximately corresponds to the frequency of the lump's radial oscillation, since the lump follows a slightly eccentric path around the cavity (Armengol et al. 2021; Noble et al. 2021); at the point of closest approach, the rate of matter falling into the cavity increases and so does the mini-disk luminosity. A second, smaller peak at $f\sim 1.4{f}_{{ \mathcal B }}$ is visible primarily in X-rays. This is twice the beat frequency and corresponds to the accretion event that happens when one of the two mini-disks passes near the lump. In the lower panel of Figure 10, we show the PSD for one mini-disk alone. The peak at the beat frequency is quite visible here, and it is even greater than the one associated with the lump's radial oscillation.

Zoom In Zoom Out Reset image size

In Figure 11, we show the PSD for S0. The three curves in the upper panel, which represent the total emission at different frequencies, have some differences with respect to the spinning case. The low-frequency variability is less predominant since S0 runs for a shorter period of time, not long enough to resolve this periodicity well. The peak at ∼2f_beat, on the other hand, is more intense than in S06 because the mini-disks are less massive and deplete completely during a beat period. Another distinct feature of this peak is that it is shifted to higher frequencies in the X-ray curve. This appears to be an effect of the binary frequency increasing during the inspiral. If we account for this time-dependent orbital frequency in the light-curve data by redefining the unit of time to be the instantaneous binary orbital frequency, the peak in the PSD shifts back to f ∼ 2f_beat (see the green dashed curve in the upper panel of Figure 11). This periodic signal is likely stronger at later times when the orbital period is shorter than 〈P〉. The lower panel, which depicts the PSD for a single mini-disk, shows no significant differences from the spinning case.

Zoom In Zoom Out Reset image size

Single supermassive black holes accreting at moderately high rates are thought to consist of a geometrically thin, optically thick, and radiatively efficient disk that emits a multitemperature blackbody spectrum peaking in the UV band, and an optically thin hot corona that emits a power-law spectrum at higher energies. The simplest and most used model for the thermal component is the Novikov–Thorne (NT) disk model (Novikov & Thorne 1973), defined by the assumptions of nearly-circular orbits, time-steadiness, azimuthal symmetry, local radiation of dissipated energy, and no stress inside the ISCO.

To determine whether NT models provide a good approximation to the emission from an accreting binary black hole system approaching merger, we compare the ray-traced SED from simulation S06, averaged during the fifth orbit, to the face-on SED of various NT models and show the results in the left panel of Figure 12. We calculate the NT flux as

$$\begin{eqnarray}&&{F}_{{\nu }_{\mathrm{obs}}}=\displaystyle \frac{2\pi }{{d}^{2}}{\int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}{dr}\,r\,{g}^{3}(r)\,{B}_{{\nu }_{\mathrm{em}}}({T}_{\mathrm{eff}}(r)),\end{eqnarray} \tag{ 10 }$$

where ν_obs = g ν_em and

$$\begin{eqnarray}&&g:=\displaystyle \frac{{({k}_{\mu })}_{\mathrm{obs}}{u}_{\mathrm{obs}}^{\mu }}{{({k}_{\nu })}_{\mathrm{em}}{u}_{\mathrm{em}}^{\nu }}=\sqrt{| -{g}_{{tt}}-2{\rm{\Omega }}{g}_{t\phi }-{{\rm{\Omega }}}^{2}{g}_{\phi \phi }| },\end{eqnarray} \tag{ 11 }$$

is the redshift factor for face-on emission. Here, g_{μ ν} are the components of the Kerr metric in Boyer–Lindquist coordinates and Ω is the orbital velocity of a Keplerian circular orbit in the Kerr spacetime. Equation (10) neglects light bending, but this phenomenon has little effect on the spectrum for face-on emission.

Zoom In Zoom Out Reset image size

SMBBH emission comes from both the CBD and the mini-disks. Not surprisingly, given the CBD's state of quasi-inflow equilibrium and its distance from the nominal ISCO (corresponding to the total mass of the binary), its emission averaged over the fifth orbit is well reproduced by an NT disk extending from R_in,circ = 2〈r₁₂〉 ∼ 38M to R_out,circ = 150M, accreting at a rate of $0.25{\dot{M}}_{\mathrm{Edd}}$, which are the same values used in our simulations. Quantities enclosed in brackets are averaged during the fifth orbit. Even though we integrate from the inner edge to the outer part of the disk, for this NT disk model, we have set the stress to zero at R_ISCO,circ = 6M, which corresponds to the system's fictitious ISCO.

To analyze the SED of the mini-disks, we compare them to three different NT models with varying accretion rates. None is a good match to the spectrum we calculate. In all cases, we set the outer radius equal to the truncation radius, R_out,md ∼0.4〈r₁₂〉 ∼ 8M = 16m_i , the inner radius equal to the individual ISCOs, R_ISCO,md(a = 0.6) = 3.8m_i = 1.9M, and the mass and spin equal to those of the black holes in the simulation: (m_i = 0.5 × 10⁶ M_⊙, χ_i = 0.6). At a fixed accretion rate in Eddington units, the frequencies of features are ∝ M^−1/4.

In Figure 12, MD Model A (red curve) represents the case in which the accretion rate onto both mini-disks is equal to that in the CBD region, $0.25{\dot{M}}_{\mathrm{Edd}}$. For simplicity, we assume that it divides evenly between the two mini-disks. ¹¹ The spectrum obtained is ∼3.7 times brighter than the one obtained from S06. One source of this large discrepancy is a breakdown in the NT model assumption of inflow equilibrium: the accretion rate in the mini-disks is approximately half that in the CBD region. More precisely, the averaged accretion rates onto the black holes during the fifth orbit are ∼$6.8\times {10}^{-2}{\dot{M}}_{\mathrm{Edd}}$ and $6\times {10}^{-2}{\dot{M}}_{\mathrm{Edd}}$ when the CBD accretion rate is $0.25{\dot{M}}_{\mathrm{Edd}}$.

However, this accretion rate contrast does not completely explain the shortfall. Model B shows the combined spectrum of two NT mini-disk models with the actual accretion rates. They are still ∼1.7 times brighter than the numerical spectrum, indicating that the mini-disks have a lower radiative efficiency than the NT disk. At least part of this diminished radiative efficiency is due to some of the accreting matter at each radius having less angular momentum than the value required for a circular orbit at that radius, i.e., l(r) < l_K(r). This material, which follows a decidedly noncircular orbit, is able to reach the event horizon with higher orbital energy (lower binding energy) than matter following stable circular orbits. To distinguish the luminosity from the fluid that follows quasi-circular orbits from that radiated by the fluid on noncircular orbits, we define the "circularized" accretion rate as the rate delivered by matter with l(r) ≥ l_K(r), and averaging from r_ISCO to r_trunc. The "circularized" accretion rates are $2.3\times {10}^{-2}{\dot{M}}_{\mathrm{Edd}}$ and $2.6\times {10}^{-2}{\dot{M}}_{\mathrm{Edd}}$ for the two black holes, respectively.

Model C shows the spectrum for two NT disks with the circularized accretion rates of the real mini-disks, but this model still departs from the simulation spectrum in significant ways; its luminosity is a factor ∼1.5 lower than the simulated SED, and it is significantly "softer" at frequencies above the peak.

In fact, the three analytical models all produce softer thermal spectra, which could be an effect of the higher temperatures achieved by the shocked plasma in the mini-disks. Another difference is that in the simulation, part of the cooling occurs above the thermalized photosphere producing a different spectrum: a power law up to hard X-rays. In fact, the absence of hard X-rays with significant luminosity is one of the principal failings of the NT model prediction for ordinary AGNs.

The analysis above shows that a typical mini-disk behaves differently than an NT disk onto a single black hole for the separations considered. Moreover, analytical models designed for larger separation binaries, in which the mini-disk accretion rate matches the CBD accretion rate and the angular momentum of all material delivered to the mini-disks is large enough that it cannot plunge, predict a clear "notch" in the thermal spectrum (see, e.g., Roedig et al. 2014) due to the absence of emission from the region 0.4r₁₂ < r < 2r₁₂. In contrast, the low radiative efficiency of the mini-disks at these separations completely overcomes this effect in our simulation. This is likely to change with increasing spins or at larger black hole separations.

It is also worthwhile to compare the predicted emission from an accreting binary black hole system with that of an accreting single black hole system, given the same global parameters: the total mass, outer radius, and accretion rate. Detailed GRMHD simulations of single black hole accretion disks show that the stresses in the disk do not completely vanish at the ISCO, and the dissipation profile deviates from NT. Schnittman et al. (2016) derived a very good analytical fit for the radial luminosity profile of single black hole disks in a steady state, using ray-tracing of GRMHD simulations with cooling similar to that used here (see Equation (3) in the Appendix of that paper). They found that ∼10% of the emission arises from an optically thin corona rather than from the optically thick disk. Based on these results, we show in the right panel of Figure 12 a comparison between the SED from our simulation data and the SED from a single black hole disk model that mimics the results from the GRMHD+postprocessing simulations by Schnittman et al. (2016). Here, we use their analytical fitting function for the radial luminosity profile and assume that 10% of the luminosity is in the form of a power law with an exponential cutoff, with the same shape as in our simulation. The mass of the black hole is M = 10⁶ M_⊙, the accretion rate is $0.25{\dot{M}}_{\mathrm{Edd}}$, the inner radius is R_in = 1.2r_H ≈ 2.16M, and the outer radius is R_out = 150M.

The total SED from the binary looks significantly different from that of the single black hole: the latter has a single broad peak whose luminosity is ∼3 times greater than the binary's. Moreover, the binary's spectrum peaks at a frequency ∼3.4 times smaller than the single black hole's spectrum and has a different slope between the frequency of the maximum and the frequency where the corona starts to dominate. In this region, the binary's SED is a broken power law with a break at the mini-disks' peak. Above this frequency, the spectrum softens but not so abruptly as in the case of a single black hole disk.

The coronal luminosity is also lower for the binary because the bulk of the coronal emission comes from the mini-disks, where the accretion rate is lower than in the CBD and the overall radiative efficiency is lower than that of the single black hole disk (See Section 4.2). These two facts also translate into the fact that the ratio of the luminosity from the inner region (r ≲ 20M) to the total luminosity (up to r = 150M) is ≲0.2 for the binary, whereas it is ≳0.5 for a single black hole accretion disk.

We have assumed an accretion rate of $0.25{\dot{M}}_{\mathrm{Edd}}$ in the CBD, which would correspond to a luminosity of 0.25L_Edd for a typical AGN with a standard radiative efficiency of 10%. This approximately coincides with the face-on value predicted by the NT model for a radiatively efficient disk onto a black hole with normalized spin ∼0.5. The bolometric luminosity, which is dominated by the CBD emission, is ≈0.1L_Edd in our simulations. Thus, the global radiative efficiency is ≈4%. According to the results discussed in Section 4.1, this contrast is the result of a combination of two effects: the lower accretion rate in the cavity and the mini-disks' low radiative efficiency.

To examine the mini-disks' radiative efficiency in greater detail, we choose one black hole, BH₁, and calculate the accretion rate onto it in CU as

$$\begin{eqnarray}&&{\dot{M}}_{{\mathrm{BH}}_{1}}(t)={\oint }_{{r}_{{\rm{H}}}}d\bar{A}\,{u}^{\bar{r}}\rho ,\end{eqnarray} \tag{ 12 }$$

where the overbars denote harmonic coordinates centered at the hole. We transform the accretion rate to physical units and calculate the accretion power as a function of time: ${L}_{\mathrm{acc},{\mathrm{BH}}_{1}}(t)\equiv {\dot{M}}_{{\mathrm{BH}}_{1}}(t){c}^{2}$. Then, we calculate the radiative efficiency of the mini-disk as ${\eta }_{\mathrm{eff}}\equiv {L}_{{\mathrm{MD}}_{1}}/{L}_{\mathrm{acc},{\mathrm{BH}}_{1}}$.

In Figure 13, we show the mini-disk radiative efficiency as a function of time for both simulations. The upper, middle, and lower panels depict the total (ν ∈ (0, ∞) Hz), thermal (ν ∈ (0, 2 × 10¹⁷) Hz), and coronal (ν ∈ (2 × 10¹⁷, ∞) Hz) efficiencies, respectively. The three plots show a quasiperiodic oscillatory behavior of the efficiency, with the effect being more pronounced for the S06 simulation. The spin of the black holes has little effect on the CBD in terms of modifying the bolometric emission, but it does have an effect on the emission from the mini-disk. The mean values for the total, thermal, and coronal radiative efficiency of the mini-disks are ≈1.9% (1.4%), ≈1.5% (0.8%), and ≈0.4% (0.7%), respectively, for the S06 (S0) simulation. For an NT disk extending from the ISCO to ∼0.4〈r₁₂〉, where 〈r₁₂〉 ≈ 18M = 36m₁, the radiative efficiency calculated in the same way as in our simulation (face-on emission) is η_NT(χ = 0.6) ≈ 4% and η_NT(χ = 0) ≈1.7% for a spinning and a nonspinning black hole, respectively. The thermal radiative efficiency of the mini-disks is ∼40%–60% of that predicted by the NT model for the same spin and disk extension in both S06 and S0. On the other hand, the total radiative efficiency of the mini-disks relative to the NT prediction is a fraction ∼55% for S06 and ∼90% for S0.

Zoom In Zoom Out Reset image size

The greater fraction of NT efficiency generated in the S0 mini-disks is due to the persistence of a coronal region even when the surface density in the mini-disk is relatively small. In other words, for equal accretion rate, more matter is below the thermalized Thomson photosphere in S06, resulting in a lower fraction in the corona, and vice versa for S0. These results might drastically change during other stages of the inspiral. In particular, at earlier stages and beyond some black hole separation, the single disk's predictions probably would be recovered.

Previous attempts to model EM observables from SMBBHs have been based upon 2D-viscous Newtonian simulations (Farris et al. 2015a, 2015b; Tang et al. 2018; Westernacher-Schneider et al. 2021), 2D-inviscid GRHD simulations of isolated mini-disks in a Kerr spacetime, i.e., without any tidal gravity from the companion (Ryan & MacFadyen 2017), and fully relativistic 3D-MHD simulations (d'Ascoli et al. 2018).

Our work agrees with the Newtonian 2D-viscous calculations in finding that the shock created when stream material is propelled back out to the CBD can create strong enough heating to influence the output spectrum. However, we disagree in several respects. The thermal peaks of the spectra in Farris et al. (2015a) and Tang et al. (2018) are shifted to higher frequencies compared to our results (see, e.g., Figure 2 in Tang et al. 2018) because the effective temperature in these works is much higher than our estimates. The reason for this is that they derive the effective temperature assuming the disk pressure is dominated by gas pressure. If radiation pressure is significant, as it likely is for near-Eddington accretion rates and separations out to several hundred gravitational radii, this approach overestimates the effective temperature. In our simulations, the effective temperature is linked directly to the actual dissipation rate (from shocks and turbulence). Note that Westernacher-Schneider et al. (2021) made efforts to account for this effect, but also assumed a very high accretion rate, $10{\dot{M}}_{\mathrm{Edd}}$.

In some of this work, a spectral notch appears (Tang et al. 2018) but not in others (Farris et al. 2015a). When it is absent, it is because radiation from the shocked streams fills the gap between the CBD and mini-disk peaks. In our simulations, although the stream emission does fall in this gap, it is not strong enough to fill it; the notch disappears because the mini-disks are too faint. In the Newtonian calculations, the mini-disks are much more luminous relative to the CBD because they are able to maintain inflow equilibrium and do not permit rapid inflow because their truncation radii are very large compared to their ISCOs. A final point of contrast in predicted spectra is that these papers considered only thermal emission, whereas we included also coronal X-rays.

There are also points of both agreement and disagreement in the predicted light curves. The equal-mass circular cases considered by Westernacher-Schneider et al. (2021) exhibited a modulation at the lump orbital frequency and at twice the beat frequency between the lump orbital frequency and the binary orbital frequency, just as in our work. However, when they weakened the coupling between the mini-disks' outer edges and their sink regions by refining the grid and shrinking the sink region, they found that the beat frequency feature shifted to the binary orbital frequency. Thus, they implicitly confirm the presence of this feature in the mass-delivery rate, a feature found in many previous calculations (Noble et al. 2012; Shi et al. 2012; Bowen et al. 2019), but the dynamics producing the binary orbital frequency modulation in the light curve remain undetermined.

We extended the calculations made by d'Ascoli et al. (2018) in two ways: adding a case with spinning black holes to the Schwarzschild case examined in that paper and lengthening the run from three orbits to more than 11. The greater duration enabled us to identify truly periodic features in the emission after the system had reached a quasi-steady state.

We have explored various features of the spectrum and its time variability for a binary black hole system with a total mass of M = 10⁶ M_⊙. We chose this value considering that mergers of these systems are future LISA targets and thus candidates for future multimessenger observations. However, because an AGN's luminosity is expected to scale with the black hole mass, heavier SMBBHs may be detected via their EM emission at greater distances. Additionally, as mass increases, the peak of the thermal spectrum shifts to lower energies, allowing optical/UV telescopes to observe them. As an illustration, Figure 14 shows the SED from S06 at t = 3290M for five different the total masses: M = 10⁵, 10⁶, 10⁷, 10⁸, and 10⁹ M_⊙. The total luminosity scales as L ∝ M, whereas the thermal peak scales as ν_peak ∝ M^−0.25. The orbital period is also ∝ M and, at the separations considered, ranges from ∼few minutes for M = 10⁵ to >30 days for M = 10⁹. This large range of timescales implies the need for different observational strategies to detect SMBBHs with different masses. On the other hand, we must not only detect SMBBHs but also distinguish them from the many more single AGNs in the universe. The best way to accomplish this is to identify distinct periodic modulations in the emission.

Zoom In Zoom Out Reset image size

Using the results of the S06 simulation, we investigate two scenarios: one with a total mass of 10⁹ M_⊙ and another with a mass of 10⁶ M_⊙. In the first case, the two periodic signals we identified at ∼$0.2{f}_{{ \mathcal B }}$ and ∼$1.4{f}_{{ \mathcal B }}$ would have periods of ∼150 days and ∼20 days, respectively. If such a source is identified as an SMBBH candidate, a follow-up observational campaign with observations every <10 days may detect these variabilities and help to confirm the nature of the source. This could be accomplished using NICER ¹² ; since the NICER X-ray instrument is on board the International Space Station in a low Earth orbit, depending on the visibility, it can monitor X-ray sources up to 16 times per day. In Figure 15, we show a simulated light curve obtained by convolving our theoretical prediction with response matrices from NICER using the software XSPEC (Arnaud 1996). We simulate a scenario in which recurrent observations of 2 ks are made every ∼5 days. We assume a distance of 500 Mpc from the source and a typical Hydrogen column density of N_H = 10²² cm⁻². To detect a source of this type at a farther distance, longer exposure times would be required.

Zoom In Zoom Out Reset image size

On the other hand, if the total mass of the system is 10⁶ M_⊙, the total duration of our simulation corresponds to 10 hr. This period can be covered by a single observation, which should be divided appropriately to account for the light curve's variability. Athena, ¹³ which will be launched in the next decade, is the best future X-ray facility for studying this type of source. Athena's effective area will be orders of magnitude larger than those of today's X-ray observatories. Assuming once more a hydrogen column density of N_H = 10²² cm⁻², the variabilities could be discerned up to a distance of d ∼ 50 Mpc.

With regard to this timing analysis, the nonspinning situation would be very similar, with the only difference being that the signal at ∼$1.4{f}_{{ \mathcal B }}$ would be stronger than in the spinning case.

A further question concerns the prospects for observing EM emission during the merger itself. The luminosity of this event is proportional to the amount of gas that may be heated during the merger. The total mass contained within the cavity at the end of our simulations provides an upper bound on this quantity. We estimate this by taking an average of the mini-disks' mass during the final 2000M of each simulation. For a total black hole mass of 10⁶ M_⊙ and an accretion rate of $0.25{\dot{M}}_{\mathrm{Edd}}$, the average mass in the cavity is ∼5.7 × 10²⁵ g and 3.6 × 10²⁵ g for S06 and S0 simulations, respectively. These two values are well-fitted by the following phenomenological expression:

$$\begin{eqnarray}&&{M}_{\mathrm{cav}}\lesssim 1.42\times {10}^{26}{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{10}^{6}{M}_{\odot }}\right)}^{2}\left(\displaystyle \frac{\dot{M}}{{\dot{M}}_{\mathrm{Edd}}}\right)\left(1+\chi \right)\,{\rm{g}}.\end{eqnarray} \tag{ 13 }$$

The mass and accretion rate dependency in Equation (13) are directly derived from the unit conversion between the numerical and the physical mass units, whereas the spin dependence is the simplest linear fit for the two values considered (0 and 0.6) and should not be extrapolated directly to other spin values. For an optimistic scenario involving a 10⁹ M_⊙ SMBBH system accreting at a rate of ∼${\dot{M}}_{\mathrm{Edd}}$, Equation (13) yields a mass upper bound of M_cav ∼ 1.5 × 10³² g, which corresponds to a maximum possible energy ∼10⁵³ erg, ∼10⁴ × the total radiated energy in a supernova.