Quasi-Periodicity of Supermassive Binary Black Hole Accretion Approaching Merger

TL;DR

This study investigates the quasi-periodic accretion behavior in supermassive binary black holes nearing merger, revealing resonant oscillations and modulation frequencies that could influence their electromagnetic signals.

Contribution

It extends previous simulations to longer durations, demonstrating how $m=1$ density features induce characteristic accretion modulations in merging supermassive black hole systems.

Findings

Accretion rates are modulated at beat frequencies between binary and lump orbitals.

Mini-disk masses and accretion rates are strongly modulated at specific frequencies.

System variability could produce distinctive electromagnetic fluctuations near merger.

Abstract

In this paper we continue the first ever study of magnetized mini-disks coupled to circumbinary accretion in a supermassive binary black hole (SMBBH) approaching merger reported in Bowen et al. 2018. We extend this simulation from 3 to 12 binary orbital periods. We find that relativistic SMBBH accretion acts as a resonant cavity, where quasi-periodic oscillations tied to the the frequency at which the black hole's orbital phase matches a non-linear $m=1$ density feature, or ``lump'', in the circumbinary accretion disk permeate the system. The rate of mass accretion onto each of the mini-disks around the black holes is modulated at the beat frequency between the binary frequency and the lump's mean orbital frequency, i.e., $\Omega_{\rm beat} = \Omega_{\rm bin} - \bar{\Omega}_{\rm lump}$, while the total mass accretion rate of this equal-mass binary is modulated at two different frequencies, $\gtrsim \bar{\Omega}_{\rm lump}$ and $\approx 2 \Omega_{\rm beat}$. The instantaneous rotation rate of the lump itself is also modulated at two frequencies close to the modulation frequencies of the total accretion rate, $\bar{\Omega}_{\rm lump}$ and $2 \Omega_{\rm beat}$. Because of the compact nature of the mini-disks in SMBBHs approaching merger, the inflow times within the mini-disks are comparable to the period on which their mass-supply varies, so that their masses---and the accretion rates they supply to their black holes---are strongly modulated at the same frequency. In essence, the azimuthal symmetry of the circumbinary disk is broken by the dynamics of orbits near a binary, and this $m=1$ asymmetry then drives quasi-periodic variation throughout the system, including both accretion and disk-feeding. In SMBBHs approaching merger, such time variability could introduce distinctive, increasingly rapid, fluctuations in their electromagnetic emission.

Tables

Table 1. Table 1 Warped Grid Parameters

Parameter	Value
${\delta }_{x1}$, ${\delta }_{x2}$, ${\delta }_{x3}$	0.2
${\delta }_{x4}$	0.1
${\delta }_z$	0.4
${\delta }_{y3}$, ${\delta }_{y4}$	0.15
$h_{x1}$, $h_{x2}$, $h_{x3}$, $h_{x4}$, $h_z$	20.
$h_{y3}$, $h_{y4}$	10.
$s_1$, $s_2$, $s_3$	0.01
$b_1$, $b_2$, $b_3$	15.
$a_{x1}$, $a_{x2}$	4.0
$a_z$	4.3
$R_{out}$	$260M$

Equations

$${\partial }_t\mathbf{U}\left(\mathbf{P}\right)=-{\partial }_i{\mathbf{F}}^i\left(\mathbf{P}\right)+\mathbf{S}\left(\mathbf{P}\right),$$

$$\mathbf{U}\left(\mathbf{P}\right)$$

$${\mathbf{F}}^i\left(\mathbf{P}\right)$$

$$\mathbf{S}\left(\mathbf{P}\right)$$

$$T_{\alpha \beta }=\left(\rho h+2p_m\right)u_{\alpha }{u}_{\beta }+\left(p+p_m\right)g_{\alpha \beta }-b_{\alpha }{b}_{\beta },$$

$${\mathcal{L}}_c=\frac{\rho ϵ}{t_{\mathrm{c}ool}}{\left(\frac{\Delta S}{S_0}+\left|\frac{\Delta S}{S_0}\right|\right)}^{1/2},$$

$$D_m={\int }_{r_{\min }}^{r_{\max }}\rho \sqrt{-g}{e}^{-im\varphi }{d}^{3}x.$$

$${\varphi }_m={\tan }^{-1}\left(-\frac{\Im \left(D_m\right)}{\Re \left(D_m\right)}\right),$$

$${\Omega }_{\mathrm{lump}}={\partial }_t{\varphi }_{m=1}.$$

$${\mathcal{M}}_i={\int }_{r=m_i}^{0.4a(t)}\rho u^0\sqrt{-g}{d}^{3}x,$$

Full text

††slugcomment: Submitted to ApJ

1. Introduction

In addition to the stellar-mass binary black holes (BBHs) already detected by LIGO (Abbott et al., 2016a, b, 2017, 2017, 2018b, 2018a), supermassive binary black holes (SMBBHs) are expected to form during galactic mergers (see(Khan et al., 2016; Kelley et al., 2017) for more details). While there remains uncertainty about the exact processes which extract sufficient angular momentum from the binary to create separations of $\mathcal{O}\left(10^{3}\right)$ gravitational radii (Begelman et al., 1980), interactions with both stars (Khan et al., 2011; Vasiliev et al., 2015; Gualandris et al., 2017) and gas (Shi et al., 2012; Dotti et al., 2012) may overcome this so-called “final parsec problem” (but see, e.g., Muñoz et al. (2019); Moody et al. (2019)). Beyond this point, gravitational wave emission can be expected to drive the binary to merger within a Hubble time(Milosavljević & Phinney, 2005). Unfortunately, the detection of gravitational waves from SMBBHs awaits the efforts of other observational campaigns, such as LISA (Amaro-Seoane et al., 2017) and Pulsar Timing Arrays (Shannon et al., 2015).

However, the galactic environments in which SMBBHs form should provide ample amounts of gas for them to be electromagnetically bright(Cuadra et al., 2009; Chapon et al., 2013; Colpi, 2014). With this consideration, and the plethora of current and upcoming electromagnetic survey campaigns (LSST Science Collaboration et al., 2009; Chambers et al., 2016), SMBBHs represent an excellent candidate for multimessenger astrophysics following the launch of LISA (Kara et al., 2019; Baker et al., 2019). However, the chief question remains: “What are the electromagnetic signals associated with SMBBH merger?”. To zeroth order, any electromagnetic counterpart to the SMBBH gravitational wave signal will be directly related to the amount and structure of gas in the immediate vicinity of the black holes (BHs) to be gravitationally heated at merger (Krolik, 2010).

The earliest estimates predicted that SMBBHs would exhibit a dry merger, in which there would be little to no gas available for electromagnetic emission. When the binary mass ratio is near unity, torques exerted by the binary on inflowing gas open up a cavity of radius $\approx 2a$, where $a$ is the binary separation (Pringle, 1991; Artymowicz & Lubow, 1994, 1996). It was thought that this process also prevented any mass from crossing the cavity. However, more recent numerical simulations have shown that in fact the BHs efficiently peel streams of gas off the inner edge of circumbinary disks, and these streams rapidly traverse the central cavity(MacFadyen & Milosavljević, 2008; Bode et al., 2010; Palenzuela et al., 2010; Farris et al., 2011; Bode et al., 2012; Farris et al., 2012; Giacomazzo et al., 2012; Noble et al., 2012; Shi et al., 2012; D’Orazio et al., 2013; Farris et al., 2014; Gold et al., 2014; Farris et al., 2015a, b; Shi & Krolik, 2015; D’Orazio et al., 2016; Bowen et al., 2018; Tang et al., 2018).

Moreover, these streams of gas are able to form individual accretion disks, or mini-disks, around each member of the binary(Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016; Tang et al., 2018; Moody et al., 2019). Because the mini-disks represent some of the most significant departures from standard active galactic nuclei, much effort has been undertaken to understand their structure (Ju et al., 2016; Bowen et al., 2017; Ryan & MacFadyen, 2017; Bowen et al., 2018) and mechanisms by which they could produce distinct electromagnetic signatures (Roedig et al., 2014; D’Orazio et al., 2015; Bowen et al., 2017; Ryan & MacFadyen, 2017; Bowen et al., 2018; d’Ascoli et al., 2018; Tang et al., 2018).

Finally, numerical studies have shown that the inner edge of the circumbinary disk itself differs significantly from a standard accretion disk around a single BH (MacFadyen & Milosavljević, 2008; Noble et al., 2012; Shi et al., 2012; D’Orazio et al., 2013; Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016; Tang et al., 2018). As streams of gas are peeled off the inner edge of the circumbinary disk, a portion of the stream is granted angular momentum from the binary and flung into the inner edge of the circumbinary disk. This serves as an initial seed for a feedback cycle resulting in an $m=1$ azimuthal asymmetry, or lump, at the inner edge of the circumbinary disk (Shi et al., 2012). The lump quasi-periodically modifies the accretion flux into the central cavity and therefore onto the mini-disks(Shi et al., 2012; Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016; Bowen et al., 2018).

This paper picks up from where the equal-mass binary simulation of Bowen et al. (2018) left off, and seeks to explore the time-dependent structure of relativistic mini-disks for binary separations of a\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}20M111We adopt geometrized units with $G=c=1$ where distance has units of $GM/c^{2}$ and $M$ is the total mass of the binary.. The lump is located $\approx 2.4a$ from the center-of-mass and orbits at nearly the local Keplerian orbital frequency around the center-of-mass, ${\bar{\Omega}}_{\rm lump}\approx 0.28\Omega_{\rm bin}$, where $\Omega_{\rm bin}$ is the binary orbital frequency. The quasi-periodic modulations of the accretion into the central cavity due to the lump occur at twice the frequency of a BH coming into phase with the lump: $2\Omega_{\rm beat}=2\left(\Omega_{\rm bin}-\bar{\Omega}_{\rm lump}\right)$ (Shi et al., 2012; Noble et al., 2012; D’Orazio et al., 2013; Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016; Bowen et al., 2018).

While Newtonian studies do note the modulation of the accretion streams, the inflow times in the mini-disks are sufficiently long that no significant mini-disk asymmetry was achieved (Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016). However, once the binary separation shrinks to a\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}20M, the tidal truncation radius of a mini-disk is only a factor of a few times the innermost stable circular orbit (ISCO). This scenario leads to radial pressure gradients accelerating mass inflow well beyond the rate associated with stresses arising from MHD turbulence (Beloborodov & Illarionov, 2001; Krolik et al., 2005) and in turn dynamically couples the mini-disk mass to the lump (Bowen et al., 2018).

By extending the original simulation, we now numerically extract the modulation frequencies of the mini-disk masses, combined mini-disk mass, and lump orbital frequency. We find that nearly every component of the lump-mini-disk system exhibits quasi-periodic modulations associated with the frequency at which the BH’s orbital phase matches the lump’s: $\Omega_{\rm beat}\approx 0.72\Omega_{\rm bin}$. For instance, we observe that the mini-disk masses have a modulation frequency of $\Omega_{\rm beat}$. Meanwhile, both the combined mini-disk mass and the instantaneous lump angular velocity $\Omega_{\rm lump}$ experience modulation at frequencies $\approx\left(0.2-0.4\right)\Omega_{\rm bin}$ and $2\Omega_{\rm beat}$. We therefore conclude that the internal time-dependent structure of the central cavity, streams, lump, and mini-disks are all intimately coupled together. Such rich time variability could introduce significant quasi-periodic signatures in the electromagnetic emission of SMBBHs (d’Ascoli et al., 2018).

The remainder of this paper is organized as follows. In Section 2 we present the details of our simulation. In Section 3 we present the results of our analysis including the extraction of characteristic frequencies of oscillations in the lump-mini-disk system. In Section 4 we discuss the implications of our findings. Finally, in Section 5 we present our concluding remarks.

2. Simulation Details

2.1. Overview

The primary objective of our simulation is to further explore the dynamic coupling of the mini-disks around each BH to a $m=1$ azimuthal Fourier mode in the circumbinary disk reported in (Bowen et al., 2018). We continue the simulation from three binary orbital periods, where Bowen et al. (2018) ended, out to 12 binary orbital periods. Our simulation starts at a binary separation of $20M$, where relativistic contributions to the spacetime and inspiral are significant (Zilhão et al., 2015), and the circumbinary disk is taken from a snapshot of $t=50,000M$ in runSE of Noble et al. (2012).

As the total mass in the accretion disks for astrophysical SMBBHs will be negligible compared to the combined mass of the binary, we neglect the self-gravity of the fluid and any feedback into the spacetime from the matter. We approximate the spacetime of the SMBBH with a fully analytic prescription which asymptotically matches BH perturbation theory near each individual BH to post-Newtonian (PN) theory at 2.5PN order (for full details see Noble et al. (2012); Mundim et al. (2014); Ireland et al. (2016); Bowen et al. (2017)). The BH trajectories are calculated to 3.5PN order accuracy. We evolve the gas using the equations of general relativistic magnetohydrodynamics (GRMHD) in flux-conservative form using the Harm3d code (see Noble et al. (2009, 2012); Bowen et al. (2017)).

Throughout the paper, unless otherwise noted, we use geometrized units in which $G=c=1$. When used as tensorial indices, we reserve Greek letters (e.g., $\alpha,\beta,\gamma,\ldots$) for spacetime indices and Roman letters (e.g., $i,j,k,\ldots$) as indices spanning spatial dimensions.

2.2. General Relativistic Magnetohydrodynamics

The equations of motion for GRMHD on a background spacetime are expressed through conservation laws for baryon number density and stress-energy coupled to the Maxwell induction equations and the divergence free constraint on the magnetic field. These may be written as the following system of conservation laws:

[TABLE]

where ${\bf P}$ are the “primitive” variables, ${\bf U}$ the “conserved” variables, ${\bf F}^{i}$ the fluxes, and $\mathbf{S}$ the source terms. In terms of the primitive variables and metric functions they can be expressed as

[TABLE]

where $g$ is the determinant of the metric, ${\Gamma^{\lambda}}_{\alpha\beta}$ are the Christoffel symbols, $b^{\alpha}=\left(1/u^{t}\right)\left({\delta^{\alpha}}_{\nu}+u^{\alpha}u_{\nu}\right)B^{\nu}$ is the magnetic 4-vector projected into the fluid’s comoving reference frame, and $u^{\alpha}$ are the components of the fluid’s 4-velocity. The stress-energy tensor is written as

[TABLE]

where $h=1+\epsilon+p/\rho$ is the specific enthalpy, $\epsilon$ is the specific internal energy, $p$ is the gas pressure, $p_{m}=\frac{1}{2}b^{2}$ is the magnetic pressure, and $\rho$ is the rest-mass density. The initial value of the divergence constraint is maintained to machine precision using FluxCT (Tóth, 2000).

We assume the accretion disks are radiatively efficient and cool away any local increases in entropy on a timescale $t_{\mathrm{c}ool}$, specified as the local fluid orbital period (for full details on how these are calculated see (Bowen et al., 2017; d’Ascoli et al., 2018)). We include cooling as a source term to the stress-energy conservation equations: $\nabla_{\lambda}{T^{\lambda}}_{\beta}=-{\cal L}_{c}u_{\beta}$. The fluid rest-frame cooling rate per unit volume ${\cal L}_{c}$ is determined via the prescription of Noble et al. (2012) in which the gas is cooled at a rate

[TABLE]

where $\Delta S\equiv S-S_{0}$ and $S=p/\rho^{\Gamma}$ is the local entropy. Our target entropy, $S_{0}=0.01$, is the initial entropy of each accretion disk. Finally, we close the system using a gamma-law equation of state of adiabatic index $\Gamma=5/3$.

2.3. Grid and Boundary Conditions

The simulation performed in Bowen et al. (2018), which we continue here, is performed in a time-dependent, double fish-eye (warped) spherical coordinate system whose origin is at the center-of-mass (Zilhão & Noble, 2014). We tabulate the full set of grid parameters in Table 1. This warped gridding scheme facilitates a focusing of cells in the immediate vicinity of the BHs while preserving spherical coordinates in the circumbinary disk. The focusing tracks the BHs along their shrinking orbits. We plot an equatorial and poloidal slice of the grid used in our simulation in Figure 1. For the poloidal grid, we employ the same polynomial focusing described in Noble et al. (2012).

We excise the coordinate singularity at the origin by placing a spherical cutout of radius $r=2M$ at the center-of-mass. The cutout size was chosen to sufficiently enlarge the timestep of our simulation to allow it to be feasible with available computational resources. Although this cutout removes the ability to study any transfer of material between the mini-disks, or sloshing (Bowen et al., 2017), it should not impact the physics of the mini-disk–stream–lump interaction presented here significantly. This is because at any given moment \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1\% of the combined mini-disk mass is in the cutout/sloshing region, and in a dynamical time only a fraction of it enters the cut-out (Bowen et al., 2017). The radial extent of the computational domain was selected as $13a(t=0)$ to fully encompass the circumbinary disk of Noble et al. (2012) used as initial data. We employ outflow boundary conditions on the radial boundaries, enforcing that $u^{r}$ be oriented out of the domain. If not, then we reset the radial velocity to zero and solve for the remaining velocity components. Our poloidal “$x^{2}$” coordinates have reflective, axisymmetric boundary conditions at the polar axis cutout and our azimuthal “$x^{3}$” boundary conditions are periodic.

The grid contains $600{\times}160{\times}640$ ($x^{1}{\times}x^{2}{\times}x^{3}$) cells. In the region of the circumbinary disk, it exactly matches the grid used in the simulations of Noble et al. (2012). Due to the polar grid spacing and off-grid-center location of the BHs, our grid does not include a full 32 cells per scale height in the mini-disks on the side of the BH farthest from the center-of-mass. Figure 1 illustrates this asymmetry. However, its cell density is not far below this number and is sufficient to fully resolve the magnetorotational instability in the circumbinary, so that MHD effects facilitate the accretion of material into the central cavity (Noble et al., 2012). Increasing the poloidal cell count would have further decreased our timestep and required an increase in the central radial cutout size to compensate. Finally, the azimuthal cell count in the equatorial plane was selected to be sufficiently large to be comparable to previous hydrodynamic studies (Bowen et al., 2017).

3. Results

3.1. Overview

Our simulation starts at a binary separation of $20M$, where PN corrections and inspiral are significant (Zilhão et al., 2015). Over the course of the 12 binary orbital periods simulated here, the binary orbital frequency increases by $\approx 24\%$ and the binary separation decays by $\approx 15\%$ to $\approx 17M$. The lump in the circumbinary disk quasi-periodically modulates the accretion into the central cavity at twice the BH-lump beat frequency, as each BH comes into orbital phase with the lump and drains mass from it: $2\Omega_{\rm beat}=2\left(\Omega_{\rm bin}-\bar{\Omega}_{\rm lump}\right)\approx 1.44\Omega_{\rm bin}$ (Shi et al., 2012; Noble et al., 2012; D’Orazio et al., 2013; Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016; Bowen et al., 2018). The association of accretion with the lump creates an asymmetry in the mass flux carried by the accretion streams, despite the nominal mirror symmetry of the equal-mass binary (see Figure 2) (Bowen et al., 2018).

As was the case in the beginning of the simulation, we observe that at any given time throughout the 12 orbits one mini-disk is the beneficiary of the primary accretion stream, while the other mini-disk is starved of mass supply. Unlike previous Newtonian studies of mini-disks (Farris et al., 2014, 2015b; D’Orazio et al., 2016), our binary separation is so small that the tidal truncation radii of the mini-disks are only a factor of a few times the ISCO radius (r_{\rm t}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}2.4r_{\rm ISCO}). The radially compact nature of the mini-disks in our simulation reduces the inflow times to values of order the BH-lump coupling period (Beloborodov & Illarionov, 2001; Krolik et al., 2005; Bowen et al., 2018). Therefore, during this coupling period one mini-disk is capable of depleting nearly all of its mass, while the other mini-disk is reformed from significant mass supply as the accreting stream taps into the over-density of the lump.

Beyond the previously presented three binary orbital periods, we find an increased clarity of the effects of the dynamic coupling between the lump and mini-disks. Furthermore, we observe that the state of the mini-disks throughout the inspiral more closely resembles the most extreme case presented in  Bowen et al. (2018), where nearly $75\%$ of the total mini-disk mass can be contained within a single disk at the peak of the cycle. In Section 3.3 we present the long term behavior of the mini-disk refilling and depletion cycle. To understand this cycle in the context of the mass supply asymmetry, we first study the azimuthal Fourier structure of the circumbinary disk in Section 3.2. Finally, we find that the modulation frequencies of mini-disk mass and the lump are all tied together through $\Omega_{\rm beat}$ (see Section 3.4).

3.2. Circumbinary Disk Azimuthal Structure

Azimuthal structure can be characterized in terms of its Fourier modes: $f\left(\phi\right)=\sum_{m}D_{m}e^{im\phi}$. We compute the amplitude of the $m$-th component of the circumbinary disk’s density as (see e.g. (Zurek & Benz, 1986; Heemskerk et al., 1992))

[TABLE]

We set our radial inner boundary of the integration as the location of the inner edge of the circumbinary disk, $r_{\rm min}=2a(t)$. At large radii, turbulence acts to disrupt the spiral density structure present within the circumbinary disk. This effectively removes amplitude in $m\neq 0$ while the larger integration volume increases the amplitude of $m=0$. As we are interested in the normalized amplitudes, $D_{m}/D_{0}$, we elect to impose an outer radial boundary on our integration to eliminate $m=0$ noise. The outer radial boundary on our integration is set to $r_{\rm max}=4a(t)$. This is roughly where the density in the circumbinary disk begins to fall off by eye in our density contours (see Figure 2 at $r\approx 80M$), and is well beyond the $2.4a(t)$ where we expect the lump $m=1$ contribution.

As shown in Figure 3, the spectrum of these modes is dominated by $m=1$, the primary contributor to the lump. Meanwhile, the gravitational quadrupole drives $m=2$ spiral density waves in the circumbinary disk. The normalized strength of the $m=1$ component can be as much as five times greater than the $m=2$ component at the same time. Next, although it is beyond the scope of this paper, we note that there appear to be high frequency oscillations within the dominant $m=1$ and $m=2$ components. Finally, there may also be a low frequency modulation in the overall strength of the lump itself, much like that seen in the simulation of Shi et al. (2012).

The primary advantage of our Fourier series expansion is that we can now use it to track the phase of the lump in time. We calculate the lump phase with the $m=1$ azimuthal pattern location as (Zurek & Benz, 1986; Heemskerk et al., 1992)

[TABLE]

where $\Im\left(D_{m}\right)$ and $\Re\left(D_{m}\right)$ are the imaginary and real parts of the $m$ component of our Fourier series expansion. In the top frame of Figure 4 we plot $\phi_{m}$ for the lump and the phase of one of the BHs as a function of time for our simulation. A point where the two lines cross would correspond to a BH coming into phase with the lump.

In Figure 4 we observe that the lump phase is not advancing at a constant rate. To highlight this, we define the instantaneous lump angular velocity as the $m=1$ pattern speed,

[TABLE]

We plot this angular velocity in ratio to the instantaneous binary orbital frequency in the bottom frame of Figure 4. First, we note that there does not appear to be any significant secular trend in the ratio $\Omega_{\rm lump}/\Omega_{\rm bin}$. This implies that the lump is readily able to track the binary inspiral throughout our simulation. Second, there are quasi-periodic oscillations in the ratio that perturb it by as much as 30–40% from its mean value $\approx 0.28$. Finally, we observe that the mean value of the ratio between the lump’s angular velocity and the binary orbital frequency is very close to the ratio that would obtain if the lump followed a Keplerian orbit with semi-major axis $2.4a(t)$, i.e., $0.27$. On this basis, as well as seeing that the lump’s density concentration moves very coherently, we argue that the lump is better thought of as a physical object following an orbit rather than a pattern through which orbiting fluid moves.

3.3. Mini-Disk Refilling and Depletion Cycle

We plot an equatorial slice of the gas density at $t=5960M$ in Figures 2 and 5. The latter figure shows clearly how large the contrast in mass between the two disks can be, a contrast stemming from the alternation in feeding between the two. In Figure 2, a dense band can be seen coming off the inner edge of the circumbinary disk on the right side and connecting to the BH on the left. This mini-disk has just past its time of maximum mass, so its mass supply, having exceeded its internal accretion rate for the past $\approx(\pi/2)/\Omega_{\rm beat}$ has now fallen a little below balance with its internal accretion rate. Conversely, the left side of the circumbinary edge produces a significantly less dense stream that deposits material onto the BH on the right. At the moment illustrated in the figures, the mass supply for this BH is beginning to overcome accretion mass-loss, so the mass of this disk is starting to grow.

The result of the asymmetric inflow into the central cavity is that the mini-disks exist in a persistent state of inflow disequilibrium, with their masses constantly increasing and decreasing. This can be most readily observed by looking at the fraction of total mini-disk mass that exists within a single mini-disk at a time. We define the mass contained within a mini-disk ($\mathcal{M}_{i}$) from the continuity equation as

[TABLE]

where $r=m_{i}$ is the radial location of the horizon for our non-spinning BH of mass $m_{i}$ in the PN Harmonic coordinates which we evolve in. We integrate over the full solid angle in angular coordinates. In Figure 6, we plot each mini-disk mass, the combined mini-disk mass, and $\mathcal{M}_{1}/\left(\mathcal{M}_{1}+\mathcal{M}_{2}\right)$ as a function of time for our simulation.

Observing the top panel of Figure 6, we note that the overall secular decay in the mini-disk masses appears to level out by the end of our simulation. Next, we note that there are quasi-periodic modulations in each mini-disk mass and the total mini-disk mass. However, neither of the oscillation periods of the combined mini-disk mass matches the oscillation period of the individual mini-disk masses (see Section 3.4). Finally, we observe that while the mass contained within one mini-disk increases, the mass contained within the other mini-disk decreases. This leads to quasi-periodic fluctuations in the fractional mass contained within a single mini-disk.

We now turn our attention to the mini-disk mass fractions in the bottom panel of Figure 6. First, at least $70\%$ of the total mini-disk mass is contained within a single mini-disk at six times throughout the simulation. Second, the quasi-periodic nature of the mini-disk filling/depletion cycle is fairly regular throughout the simulation and significantly more noticeable than in the first few orbits. Finally, the degree to which the mass fraction oscillates through our simulation does not appear to have significant dependence on binary separation.

3.4. Frequencies of Quasi-Periodic Behavior

Just as the azimuthal Fourier modes cleanly revealed azimuthal spatial structure, the power spectral density (PSD) picks out the principal timescales of variation. From Figure 7, we see that the time-dependence of four properties of the system—total mass of the mini-disks, the two individual mini-disk masses, and the orbital frequency of the lump—can all be described in terms of periodic behavior at three frequencies: the orbital frequency of the lump, the beat between the lump’s orbital frequency and the binary orbital frequency, and twice the beat frequency. In other words, the interaction of the lump with the binary orbit drives the mass of the mini-disks.

The total mini-disk mass is modulated at a frequency slightly greater than the lump orbital frequency and at twice the beat frequency, but not at $\Omega_{\rm beat}$, the oscillation frequency of the individual mini-disks. The reason is that the symmetry of the binary makes the filling/depletion cycles for the two BHs exactly $\pi$ out of phase and of roughly equal amplitude. However, the total mass supply into the central cavity is tied to phase-matching the lump and either partner in the binary. The total mass of the mini-disks can oscillate because the lump’s injection of mass each time it comes into phase with a BH is irregular, and the inflow time in each mini-disk is so short ($\sim 2\pi/\Omega_{\rm beat}$) that inflow equilibrium within the mini-disks is never achieved.

The two mini-disk masses vary primarily at the beat frequency, $\approx 0.72\bar{\Omega}_{\rm bin}$, where $\bar{\Omega}_{\rm bin}$ is the average binary orbital frequency from $N_{\rm orbits}=2$ to the end of the simulation (the same range over which the PSD is calculated). As already pointed out, this value corresponds to the mean beat frequency $\Omega_{\rm beat}=\bar{\Omega}_{\rm bin}-\bar{\Omega}_{\rm lump}$, when $\bar{\Omega}_{\rm lump}$ is taken from Figure 4. Although there is power outside of this peak in the PSD for $\mathcal{M}_{1}$, we do not believe this to be a robust signal because the location and structure of this peak depend on the starting point of data we include in the PSD calculation.

The lump orbital frequency (in units of the binary orbital frequency) is itself modulated at two frequencies. The low frequency modulation in the lump motion appears at $\approx 0.2\bar{\Omega}_{\rm bin}$. Due to the close proximity to the orbital frequency of the lump, we speculate that this modulation may be a result of eccentricity of the lump’s orbit in the circumbinary disk (Shi et al., 2012).

The lump motion’s high frequency modulations correspond to precisely twice the beat frequency, or the frequency at which the lump comes into phase with either BH. We speculate that this modulation is due to the time-dependent gravitational force exerted by a BH onto the lump. As an individual BH starts coming into phase with the lump, it is gravitationally pulled toward the BH in the opposite direction of the lump orbital velocity. Since the BH orbits more quickly than the lump, it then passes the lump and exerts a gravitational pull in the same direction of the lump orbital velocity. This could result in a gravitationally-driven slow down and then reacceleration of the lump with each BH passing at twice the beat frequency, $\approx 1.44\Omega_{\rm bin}$. Alternatively, this could be the result of shocks when the streams are flung back into the lump (which would also occur at twice the beat frequency).

4. Discussion

4.1. Quasi-Periodic Oscillations in Relativistic Binaries

Our simulation finds that the dynamics of the central cavity and circumbinary are tangled together, with one informing the evolution of the other in a feedback cycle. As the binary performs many orbits, the azimuthal structure of the circumbinary disk becomes increasingly more $m=1$ dominant (Noble et al., 2012; Shi et al., 2012; D’Orazio et al., 2013; Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016). However, previous simulational studies (Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016; Tang et al., 2018) have all reported that the individual mini-disks around each BH should be nearly symmetric for $q=1$ binaries.

Our results contrast with the previous work because our binary has such a small separation; there is then only a short distance between the outer and inner radii of the mini-disks (r_{\rm t}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}2.4r_{\rm ISCO}), and their inflow times become as short as the binary period. As a result, the quasi-periodic modulations in the mass supply into the central cavity directly imprint themselves onto the total mass within the central cavity, and therefore the rest-mass energy available for photon radiation. For example, the primary modulations in mini-disk mass occurs at the frequency of their respective BHs coming into phase with the lump. Meanwhile, the total combined cavity mass is modulated at the same frequencies as the lump orbital frequency. This seems to imply that the binary evolution significantly impacts the mini-disk evolution via binary–lump and lump–mini-disk coupling.

Remarkably, after a transient period lasting $\sim 3$ orbits, our simulation settles into a very regular behavior pattern in which all characteristic timescales maintain constant ratios even while the fundamental timescale, the binary orbital period, evolves. The longest timescale, the orbital period of the lump, is $\propto a^{3/2}$ exactly like the binary orbital period: this determines the slow modulation of the lump’s angular velocity and the slow modulation of the total mini-disk mass. The middle timescale is the inverse of the beat frequency, and it governs the modulation of the individual mini-disk masses. The beat frequency is $\Omega_{\rm bin}-\bar{\Omega}_{\rm lump}$; with both pieces $\propto a^{-3/2}$, the inverse is also $\propto a^{3/2}$. The shortest timescale, accounting for the rapid modulation of both $\Omega_{\rm lump}$ and the total mass, is half the inverse beat frequency, so it is likewise $\propto a^{3/2}$. Even the inflow time is related to $a^{3/2}$ because the tidal truncation radius $r_{t}\propto a$ in both Newtonian (Paczynski, 1977; Papaloizou & Pringle, 1977; Lin & Papaloizou, 1979; Artymowicz & Lubow, 1994; Mayama et al., 2010; de Val-Borro et al., 2011; Nelson & Marzari, 2016) and relativistic conditions (Bowen et al., 2017), and the inflow time is most closely tied to the orbital period at the mini-disk’s outer edge.

However, if the rapid acceleration in orbital evolution due to gravitational radiation causes the binary to shrink faster than the inner edge of the circumbinary disk, the ratio $\left(\Omega_{\rm lump}/\Omega_{\rm bin}\right)$ would exhibit a secular trend toward smaller values. This would increase $\left(\Omega_{\rm beat}/\Omega_{\rm bin}\right)$, and therefore drive more rapid cycles of mass supply and deprivation in the mini-disks.

During this evolution, the inflow rate in the mini-disks also accelerates because it is $\propto a^{-3/2}$. We expect, therefore, that the mini-disk masses would continue to follow the cycle of filling and depleting much as they have during our simulation. Eventually, however, tidal truncation of the mini-disks will come at such a small radius that they can no longer exist: $r_{t}=r_{ISCO}$ when $a\approx 8.3M$ for non-spinning BHs

As a final comment, we note that the precise frequencies, while relatively constant in units of $\Omega_{\rm bin}$, vary significantly throughout our simulation as a result of inspiral. Over the span of 12 binary orbital periods, the binary orbital frequency increases by as much as $\approx 24\%$ of its initial value. In addition to an increased binary orbital frequency, the binary separation (and therefore tidal truncation radii) have shrunk by $15\%$ of their original values. These significant changes in binary orbital parameters account for the broadened peak structure exhibited in Figure 7. Furthermore, our PSD results could be weighted more heavily towards the values at larger separations. This is because the binary spends significantly more time at larger separations than smaller due to the non-linear growth in inspiral rate from gravitational waves. It is possible that electromagnetic counterparts for merging SMBBHs, particularly in the high-energy bands emanating from the mini-disks, could be marked by a “chirp-like” growth in quasi-periodic fluctuations (Haiman, 2017).

4.2. Possible Parameter Dependence

The magnitude of mini-disk oscillations reported here depends on a number of parameters. The mini-disk cycle depends on binary separation, binary mass ratio, BH spin, and the accretion rate, which determines the run of temperature across the mini-disks and the circumbinary disk. It would therefore be presumptuous to claim that the oscillations reported here are completely generic to relativistic binaries.

Although the lump appears to be a generic result of circumbinary accretion simulations  (Noble et al., 2012; Shi et al., 2012; D’Orazio et al., 2013; Farris et al., 2014, 2015a, 2015b; D’Orazio et al., 2016; Tang et al., 2018), its quantitative structure remains subject to some uncertainty. In particular, 3D MHD simulations of circumbinary disks around non-evolving binaries have shown secular growth in its amplitude over timescales $\sim 10^{2}$ orbits (Shi et al., 2012; Noble et al., 2012). The amplitude of the lump in our simulation, drawn from a somewhat arbitrary moment in the simulation of Noble et al. (2012), may not necessarily match the amplitude of the lump found in a binary with separation $a=20M$ that evolved over a much longer period of time from an initial state of much greater separation. On the other hand, pseudo-Newtonian 2D hydro simulations with a phenomenological viscosity can follow circumbinary disks for much longer times: for example, all the way from $30M$ separation to merger (Tang et al., 2018), albeit with severe physics limitations. The fact that the lump amplitude they find (factor of several contrast in surface density) is similar to what we find suggests that the amplitude may be relatively insensitive to the arbitrary assumptions of different simulations.

The binary mass ratio could also play a significant role in the results presented here. Previous studies have shown that the strength of the lump is a function of the mass ratio of the binary. As the mass ratio deviates further from unity, the overall size of the lump in 2D hydro simulations diminishes (Farris et al., 2014; D’Orazio et al., 2016). This effect continues until a binary mass ratio of $q=0.04$, where the structure of the central cavity is altered dramatically (D’Orazio et al., 2016). More recent studies including MHD stresses find that the lump is only significant for mass ratios $q\geq 0.2$  (Noble et al., 2019).

In addition to altering the strength of the lump and mass supply fluctuations, the binary mass ratio will lead to preferential accretion onto the secondary (Farris et al., 2014; Shi & Krolik, 2015). In this scenario, the quasi-periodic oscillations in each mini-disk mass would not be the result of a flipping back and forth between each mini-disk coupling to the lump presented here, but rather a quasi-periodic modulation of mass supply to the secondary. In other words, the modulations in mass of the mini-disk around the more massive binary member would likely be diminished. Conversely, the modulations in the mini-disk around the secondary could exhibit even more pronounced fluctuations. In addition, the binary mass ratio would alter the tidal truncation radii of each mini-disk (potentially creating larger ratios of $r_{\rm t}/r_{\rm ISCO}$ for the primary’s mini-disk and smaller fractions for the secondary’s mini-disk).

There is, however, a possible self-regulation of the mass-ratio. Preferential accretion onto the secondary directly drives the mass ratio toward unity. In the relativistic regime, it is possible that there may be net transfer due to sloshing that acts in the same direction.

Another important caveat is that astrophysical BHs are expected to have spin (Lynden-Bell, 1969). The dominant effect of relevance is the modification to $r_{\rm ISCO}$ (Bardeen et al., 1972). If the angular momentum axis of the circumbinary disk, and therefore the angular momentum axis of the mini-disk, is aligned with the BH spin axis, the ratio of $r_{\rm t}/r_{\rm ISCO}$ will increase (Campanelli et al., 2006). This would serve to increase the mini-disk inflow times and prevent the mini-disk from fully depleting while starved of mass supply. It is therefore conceivable that spin effects could delay the onset of the mini-disk cycle reported here to smaller binary separations. Conversely, if the mini-disk angular momentum axis were counter-aligned with the BH spin axis, the ISCO would grow. This would either decrease the mini-disk inflow times, or prevent the formation of a mini-disk at larger binary separations. Finally, with certain misaligned spin configurations, the BH spin axes have been shown to change direction (Campanelli et al., 2007), or flip-flop, during the final stages of inspiral for mass ratios $q\geq 0.5$ (Lousto & Healy, 2015, 2016; Lousto et al., 2016; Gerosa et al., 2018). Such binary dynamics could have a profound effect on the mini-disk structure.

Furthermore, the internal stresses of the mini-disks and circumbinary disk depend on the disk temperature and resulting disk thickness. Our mini-disks may be hotter and thicker than mini-disks in nature, and therefore may have unphysically short inflow periods. However, astrophysical mini-disks will still be very hot. These high temperatures would serve to augment secondary accretion processes such as spiral shocks (Ju et al., 2016; Bowen et al., 2017; Ryan & MacFadyen, 2017; Bowen et al., 2018).

As we have outlined, the effects presented here depend strongly on the ratio $r_{\rm t}/r_{\rm ISCO}$, and therefore the binary separation. At larger separations, such as those in Newtonian studies, the inflow periods within the mini-disks will be too long for any significant depletion in the beat period. The binary separation where this mini-disk cycle takes hold remains unknown.

5. Conclusions

In this paper we have extended the first ever study of magnetized mini-disks coupled to circumbinary accretion in a relativistic SMBBH approaching merger. By extending the simulation to four times the original number of binary orbital periods, we were able to extract PSDs of the quasi-periodicity reported in Bowen et al. (2018) and explore longer term behavior through inspiral.

We found that nearly every aspect of our late-inspiral SMBBH can be characterized by quasi-periodic oscillations associated with the interaction of the individual BHs and the lump in the circumbinary disk. The individual BHs interact with the circumbinary disk by imparting angular momentum into the streams, flinging them into the circumbinary disk and forming the lump. The lump then interacts with the BHs by quasi-periodically modulating the mass accretion into the central cavity at frequencies of $2\Omega_{\rm beat}\approx 1.44\Omega_{\rm bin}$ and $\approx 0.2\Omega_{\rm bin}$.

Because of the small binary separation, the radial extent of our mini-disks is \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}2.4r_{\rm ISCO}. This leads to rapid depletion timescales in the mini-disks of order the period associated with an individual BH interacting with the lump: $\Omega_{\rm beat}\approx 0.72\Omega_{\rm bin}$. Furthermore, because the individual mini-disk masses oscillate at the beat frequency and deplete nearly all of their mass over this timescale before recoupling to the lump, we observe that the total mass within the mini-disks quasi-periodically oscillates at the same frequencies as the lump.

An astrophysical SMBBH resembling our simulation would exhibit quasi-periodicity in response to the strong fluctuations in the quantity and structure of gas in the immediate vicinity of the BHs. These quasi-periodic oscillations would be intimately tied to the structure and dynamics of the lump, and would augment Doppler variability. Finally, because our mini-disks live far outside the analytic regime of inflow equilibrium, our results emphasize the need for generating self-consistent light curves using non-linear MHD simulations.

Acknowledgments

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). For this work, D. B. B. also acknowledges partial support from NSF grants AST-1028087, AST-1516150, PHY-1305730, PHY-1707946, OAC-1550436 and OAC-1516125. M. A. is a Fellow of the RIT’s Frontier of Gravitational Wave Astronomy. M. A., M. C. and V. M. received support from NSF grants AST-1028087, AST-1516150, PHY-1305730, PHY-1707946, OAC-1550436 and OAC-1516125. 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 partially supported by NSF grants AST-1516299, PHYS-1707826 and OAC-1811287, as well as Simons Foundation grant 559794. V.M. also acknowledges partial support from AYA2015-66899-C2-1-P.

Computational resources were provided by the Blue Waters sustained-petascale computing NSF projects OAC-1811228, OAC-0832606, OAC-1238993, OAC-1516247 and OAC-1515969, OAC-0725070. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. Additional resources were provided by XSEDE allocation TG-PHY060027N and by the BlueSky Cluster at Rochester Institute of Technology. The BlueSky cluster was supported by NSF grants AST-1028087, PHY-0722703 and PHY-1229173.