Shock-Driven Periodic Variability in a Low-Mass-Ratio Supermassive Black Hole Binary

The simulations run for this work were performed on a two-dimensional, Eulerian, cylindrical grid using the FARGO3D code (Benítez-Llambay & Masset, 2016). FARGO3D solves the hydrodynamics equations

where $\Sigma$ is the surface density, $\bf{v}$ is the fluid velocity, $P=\Sigma c_{s}^{2}$ is the isothermal pressure with $c_{s}$ the isothermal sound speed, $\Phi$ is the gravitational potential, and $\Pi$ is the viscous stress tensor. We use FARGO3D’s built-in $\alpha$-prescription to set the viscosity as $\nu=\alpha c_{s}H$ (Shakura & Sunyaev, 1973), where $H$ is the local scale height of the disk.

For the transport step in the azimuthal direction, FARGO3D employs the FARGO orbital advection algorithm (Masset, 2000), which greatly increases the solution accuracy and allowable simulation timestep for systems dominated by rotation, such as accretion disks.

Our simulated domain is that of a disk centered on a primary SMBH with mass $M_{1}=10^{8}\ M_{\odot}$ and extends from $R_{\mathrm{min}}=10\;R_{g}$ to $R_{\mathrm{max}}=1000\;R_{g}$, where $R_{g}=GM_{1}/c^{2}$ is the gravitational radius of the primary SMBH. The secondary SMBH has mass $M_{2}=10^{6}\ M_{\odot}$ and is placed at a separation $a=100\;R_{g}$ from the primary, resulting in a binary orbital timescale $t_{\rm{bin}}\approx 0.1\;\rm{yr}$. We use the torque-free boundary conditions described by Dempsey et al. (2020) for our inner boundary, including the use of the de Val-Borro et al. (2006) wave-killing prescription, but allow outflow at the outer boundary with no wave-killing.

The initial conditions of the disk are derived from the steady disk equations (Frank et al., 2002),

where $\rho$ is the gas density, $\Sigma$ is the surface density, $H$ is the scale height, $c_{s}$ is the sound speed, $R$ is the radial distance from the central object, $M$ is the mass of the central object, $\dot{M}$ is the accretion rate through the disk, $P$ is the total pressure, $T_{c}$ is the temperature at the disk midplane, $\nu$ is the viscosity, $\tau$ is the optical depth, and $\kappa_{\rm R}$ is the Rosseland mean opacity, which throughout this work is assumed to be dominated by electron scattering, $\kappa_{\rm R}\approx\kappa_{\rm e}=0.4\;\mathrm{cm^{2}/g}$. We choose an $\alpha$-viscosity, $\nu=\alpha c_{s}H$ (Shakura & Sunyaev, 1973), and introduce a prescription for the scale height in order to control where the disk transitions from radiation-dominated to gas-dominated and the proportionality of $H$ to $r$ in the gas-dominated region:

Here, $m\equiv M/M_{\odot}$ is the mass of the central object in units of solar masses, $r\equiv R/R_{g}$ is the radial coordinate in units of gravitational radii ($R_{g}\equiv GM/c^{2}$), and $\dot{m}\equiv\dot{M}/\dot{M}_{\mathrm{Edd}}$ is the accretion rate in units of Eddington, where $\dot{M}_{\mathrm{Edd}}\equiv 4\pi GMm_{p}/(\eta\sigma_{T}c)$ with $m_{p}$ being the mass of a proton, $\sigma_{T}$ the Thomson cross-section, and $\eta\equiv L/(\dot{M}c^{2})$ the accretion efficiency. The parameters $r_{c}$ and $k$ set, respectively, the radius at which the disk pressure transitions from radiation to gas-dominated and the powerlaw slope of $H$ in the gas-dominated region (i.e., $H\propto r^{k}$ far from the central object). With these prescriptions for $\nu$ and $H$ in hand, we are able to solve the system in equation (3) for the surface density and sound speed distributions,

A full derivation of these expressions (equations (4)–(8)) can be found in Appendix A. For this work, we use $\alpha=10^{-2}$, $\dot{m}=10^{-1}$, $\eta=10^{-1}$, $r_{c}=300$, and $k=1$. Fig. 1 shows the surface density, temperature, and scale height distributions calculated from equations (7), (8), and (4) for our initial conditions. Note that throughout this derivation we consider $P=P_{\rm{gas}}+P_{\rm{rad}}$, the sum of the gas and radiation pressures. The resulting $c_{s}(r)$ from equation (8) used in our simulation is then an effective isothermal sound speed which includes the contributions from both the gas and radiation pressure.

The simulation is run on a cylindrical grid of $N_{r}=2560,\;N_{\theta}=3493$ cells. The radial spacing of the grid is logarithmic, with $\Delta R/R=0.0018$. This is chosen in order to resolve the 2:1 eccentric corotation resonance (ECR), which is necessary to correctly capture the evolution of eccentricity in the disk (Teyssandier & Ogilvie, 2017).

The regions of interest to this work, namely the secondary’s minidisk and the gas comprising and surrounding the gap, are non-relativistic, and so we use a purely Newtonian potential for each black hole. We apply a softening length based on the local disk scale height, giving each black hole a potential of the form

where $M_{i}$ is the mass of black hole $i$, $r_{i}$ is the distance between black hole $i$ and the point of interest, $H=c_{s}/\Omega_{K}$ is the disk scale height at that point, with $\Omega_{K}$ the local Keplerian angular velocity, and $\xi$ is a constant parameter controlling the magnitude of the softening. We choose $\xi=0.6$, following from considerations made by Duffell & MacFadyen (2013). The calculations do not include the self-gravity of the disk, as the disk is much less massive than the secondary SMBH ($M_{\mathrm{disk}}/M_{2}\sim 0.04$). The mass of the secondary is increased from $0$ to $M_{2}$ over the first 20 orbits of the simulation in order to soften spurious waves generated by the introduction of this asymmetry to the initially symmetric system.

We set the binary on a fixed, circular orbit. The disk is less massive than the secondary SMBH and so torques from the surrounding gas are not expected to significantly alter the orbit over the duration of the simulation (verified in §4.3). Likewise, we do not expect significant orbital evolution due to the emission of gravitational waves. We can demonstrate this by calculating the ratio of the simulation time to the gravitational inspiral time (Peters & Mathews, 1963),

where $N_{\mathrm{orb}}=10^{3}$ is the number of orbits in the simulation, $\tilde{r}=10^{2}$ is the initial separation in units of $R_{g}$, $q=10^{-2}$ is the mass ratio, and $f_{\mathrm{sep}}$ is the fraction of the initial separation to which the orbit decays. For a small decay to $f_{\mathrm{sep}}=0.99$ and our binary parameters, this ratio is below unity, indicating that the orbit will decay by $<1\%$ of the initial separation over the course of our simulation.

The purpose of including accretion of gas onto the secondary SMBH is twofold: first, the removal of gas from the minidisk is a natural part of the disk accretion flow, preventing spurious buildup of material at the location of the secondary; second, the mass accretion rate of an $\alpha$-disk is directly linked to its radiative emission, so monitoring this property allows an estimation of the minidisk’s luminosity over time. The secondary removes gas from within an accretion radius $r_{\rm{accr}}=0.5r_{\rm{hill}}$, where $r_{\rm{hill}}=a(q/3)^{1/3}$ is the Hill radius of the secondary, on an accretion timescale, $t_{\rm{accr}}$. $t_{\rm{accr}}$ varies as a function of distance $r$ from the secondary SMBH, based on the viscous timescale for our disk model,

where we use the same values as in our initial conditions for all parameters except that $m=10^{6}$ and $r$ is now measured in gravitational radii of the secondary. A disk-like accretion timescale such as this steeply increases with distance from the central body, causing accretion to be dominated by the inner parts of the minidisk and relatively insensitive to the specific value of $r_{\rm{accr}}$ chosen.

The fraction of mass removed from cells within $r_{\rm{accr}}$ is computed as

where $\Delta t$ is the simulation timestep. We monitor mass, momentum, and angular momentum accreted onto the secondary, binned over every $10^{-3}\;t_{\rm{bin}}$.

For shock capturing and characterization, we adapt FARGO3D’s von Neumann-Richtmyer artificial viscosity (Von Neumann & Richtmyer, 1950), using the tensor implementation described in Appendix B of Stone & Norman (1992). This artificial viscosity has the form

where $\Delta x$ is the zone size, $\nabla\bf{v}$ is the symmetrized velocity gradient tensor, $\bf{e}$ is the unit tensor, and $l$ is a parameter controlling the number of zone widths a shock is spread over. In this work, we use $l=5$. The key properties of the artificial viscosity are that it broadens shocks such that they are resolved on the grid, is sensitive only to compressive flows ($\nabla\cdot\bf{v}<0$), and is large inside shocks, but small elsewhere.

In this work, we estimate the EM radiation emitted by shocks from the energy dissipated by the shock capturing scheme, $\partial e/\partial t=-\bf{Q}:\nabla\bf{v}$. This estimate assumes that shock energy is radiated efficiently, such that the thermal energy escapes the disk and does not significantly contribute to heating the gas on long enough timescales to affect the disk evolution.

For the range of typical shock Mach numbers in the simulation ($\mathcal{M}\gtrsim 10$), the thermal timescale due to radiative diffusion $t_{\rm{diff}}=H^{2}(3C_{p}/16\sigma)(\rho^{2}\kappa_{R}/T^{3})$ (Kippenhahn & Weigert, 1990), with $C_{p}$ the specific heat at constant pressure, is less than $t_{\rm{bin}}$ in the region of interest $R<300\;R_{g}$, supporting this assumption.

Like the accretion monitoring, we bin the shock emission over every $10^{-3}\;t_{\rm{bin}}$. The Hill sphere of the secondary is excluded from the analysis of this monitoring.

In contrast to previous works at this mass ratio, we find that a $q=0.01$ binary is able to sustain a lopsided gap in the disk. In Fig. 2, we show snapshots of the disk surface density at various times and overplot ellipses showing the orbits taken by gas in the disk. The ellipse-fitting procedure is based on that of Teyssandier & Ogilvie (2017, see their §4). First, we calculate orbital elements for each cell in the simulation: the semi-major axis $a$, eccentricity $e$, and argument of pericenter $\overline{\omega}$. We then bin cells by semimajor axis and average the orbital properties of cells in the bin to obtain an overall fit to the orbit for a given value of $a$ and bin width $\delta a$. Throughout this analysis, we use a bin width of $\delta a=1.0\;R_{\rm g}$.

The eccentricity of the disk can be seen to grow over time and, as was done by Teyssandier & Ogilvie (2017), we measure this growth rate using the disk’s angular momentum deficit,

which is shown in Fig. 3. Similar to Kley & Dirksen (2006) and Teyssandier & Ogilvie (2017), we observe an early relaxing phase where the gap is opened, followed by an exponential growth phase of the eccentric mode, and, finally, saturation. The saturation of the AMD indicates the settling of the large-scale dynamics of the system and coincides with the settling of behavior in other monitored quantities (see Section 3.2), motivating our use of it as an indicator of reaching a quasi-steady state.

Fig. 4 shows the time-averaged eccentricity profile of our disk at the end of the simulation, both in terms of the radial coordinate of the simulation, $r$, and the semi-major axes, $a$, of the gas orbits. The gap edge, at approximately $r,a=200\;R_{\mathrm{g}}$, reaches an eccentricity of 0.25, larger even than the $e\approx 0.1-0.2$ measured for higher mass ratios in past works with hotter disks (Farris et al., 2014; MacFadyen & Milosavljević, 2008).

It can also be seen from Fig. 2 that the eccentric gap is precessing. Using the same orbit-fitting as was used to produce the orbital contours shown in Fig. 2, we calculate the rotation angle of the $a=200\;R_{\rm g}$ ellipse, corresponding to the gap edge, over time. From this, we obtain a precession timescale of $\sim 1.41\times 10^{3}\;t_{\mathrm{bin}}$, comparable to the $\sim 3.85\times 10^{3}\;t_{\mathrm{bin}}$ calculated from the linear theory developed by Teyssandier & Ogilvie (2016).

Fig. 5 shows the time-averaged torque density, defined in MacFadyen & Milosavljević (2008) as

and integrated torque, $T(r)=\int^{r}_{0}\left({dT}/{dr}\right)dr$, of the secondary SMBH acting on the disk, averaged over 50 $t_{\mathrm{bin}}$. The strongest peak in the torque density occurs near $r=2.5a$, coinciding with the peak in the surface density distribution of the disk. Beyond this, the torque density quickly decays to oscillations about zero. The total torque experienced by the disk is positive, with $T(\infty)\approx 7.89\times 10^{50}\;{\mathrm{g\;cm^{2}}}{\mathrm{s^{-2}}}$. The reciprocal torque of the disk on the secondary is thus negative, serving to shrink the binary separation on a timescale $t_{T}={M_{2}a^{2}\Omega_{\mathrm{BH2}}}/{T\left(\infty\right)}=3.62\times 10^{6}\;t_{\mathrm{bin}}$. This is substantially longer than the runtime of the simulation, but much shorter than the timescales computed for the accretion of momentum and the emission of gravitational waves, implying that the binary’s evolution is dominated by the gravitational torque of the disk at this time.

We can further explore where this negative torque comes from by following a similar procedure to Tiede et al. (2020) and Muñoz et al. (2019). We divide our domain into three regions based on the distribution of the torque density shown in Fig. 5: (1) an inner region $R<0.7a$, consisting of the inner disk, (2) an outer region $R>1.4a$, which includes the entire outer disk, and (3) a gap region $0.7a<R<1.4a$. It can be seen from the integrated torque in Fig. 5 that the torque is dominated by the contributions from the outer disk, with $|T_{\mathrm{out}}|\approx 9|T_{\mathrm{in}}|$, and, further, that these contributions becomes negligible past $R\sim 3a$. We plot guides showing each of these regional boundaries in the bottom-left panel of Fig. 2, which reveals that the outer region $1.4a<R<3a$ dominating the torque corresponds to the approximate extent of the overdense region near the eccentric gap edge. This is in line with the results of Tiede et al. (2020), who found that the enhanced buildup of material near the gap edge in cold disks drives the total torque on the binary to be negative.

The primary aim of this work is to investigate the time-varying electromagnetic properties of a low-mass-ratio SMBHB embedded in a realistic disk. To this end, we have employed two key proxies for variable electromagnetic emission generated from the system: shocks, tracked via dissipation in the artificial viscosity (§2.5), and accretion rate onto the secondary (§2.4), standing in for the radiative emission of the minidisk. We find clear evidence for periodicity in both of these markers, each with a period matching the binary orbital period and a peak-to-trough ratio of $\sim 1.2$.

Previous SMBHB works have found that eccentric cavities as drivers of accretion variability are not present at low mass ratios, $q<0.05$ (D’Orazio et al., 2016). Conversely, in this work, with the use of a self-consistent disk model for the surface density and scale height, we show that such behavior can be present in mass ratios as low as $q\sim 0.01$.

At a similar mass-ratio regime, simulations of super-Jupiters embedded in protoplanetary disks also produce eccentric gap systems with comparable eccentricity distributions to those in Fig.  4 (Dempsey et al., 2021). They also find that the mass ratio for which eccentricity is excited coincides with the transition from inward to outward migration due to gravitational torques. Contrary to this, we are able to develop an eccentric gap while maintaining an inward migration of the secondary. Whether this decoupling of the two phenomena is due to our disk being much thinner than the thinnest $H/R$ tested in Dempsey et al. (2021), our $H/R$ being non-uniform throughout our domain, or some other difference in the two setups is an interesting question, but beyond the scope of the present investigation.

The greatest difference between the disk in this study and those of preceding works is the disk model used for our initial conditions. In particular, our disk aspect ratio $H/R$ varies with radius and is, in general, much thinner than the constant $H/R=0.1$ disks common in the literature, with the disk having $H/R\approx 10^{-2}$ at the location of the secondary (see Fig. 1). The efficiency of gap-opening depends on the relationship between opposing torques: the tidal torque, $T_{\rm{tid}}\propto\left(H/R\right)^{-3}$, and the viscous torque $T_{\rm{visc}}\propto\nu\propto\left(H/R\right)^{2}$, which work to open the gap and to fill it in, respectively. In a colder, thinner disk, the tidal torque gets stronger and the viscous torque gets weaker, leading to wider, deeper gaps and allowing gap-opening to extend to SMBHBs with lower secondary masses (Crida et al., 2006; Duffell, 2015, 2020). This phenomenon of wider gaps in colder disks has been observed in simulations of equal-mass binaries which test different values of $H/R$ (Ragusa et al., 2016; Tiede et al., 2020), and here we demonstrate that this behavior continues to lower mass ratios.

So a colder, thinner disk makes gap-clearing more effective. In clearing a deeper gap, the secondary removes material near its orbit, weakening the damping of eccentricity occurring at eccentric corotation resonances (ECRs), enhancing the ability for the binary to drive eccentricity growth in the disk (Teyssandier & Ogilvie, 2017). In summary, thinner disks are conducive to opening deeper gaps at lower mass ratios, which in turn drives eccentricity evolution in the outer disk and the consequent variable accretion onto the binary. Given that AGN have disks with $H/R\approx 10^{-2}$–$10^{-3}$ (Krolik, 1999), consistent with the disk model used in this study, this suggests that the electromagnetic variability signatures seen in this and previous works may exist for lower mass ratios than previously indicated, increasing the population of SMBHB candidates identifiable in observations.

The underlying motivation of this work and many other simulations of SMBHBs is to better characterize the periodic electromagnetic emission which can differentiate binaries from single-SMBH AGN. Previous works generally find that high-mass-ratio binaries exhibit periodic accretion rates, but that accretion becomes steady for $q\lesssim 0.05$ (Farris et al., 2014; D’Orazio et al., 2016; Duffell et al., 2020). The timescale of accretion variability depends on the mass ratio, with lower mass ratios $0.05\lesssim q\lesssim 0.25$ being modulated only on the binary orbital timescale, while at larger mass ratios an overdense lump forms on the gap wall, causing periodic accretion on a timescale linked to its own orbital time (Farris et al., 2014; D’Orazio et al., 2013).

In contrast, we find periodic accretion onto our secondary even with $q=0.01$, with a timescale matching the orbital timescale of the binary. This matches the properties seen for binaries in the $0.05\lesssim q\lesssim 0.25$ regime in these past works, and implies that this regime may simply extend to lower mass ratios given a thinner disk. We note that the exact structure of this periodicity is dependent on the sink prescription used for the accretion. In the broadest sense, the accretion timescale used determines how quickly the minidisk around the secondary is drained. If this timescale is very short, the minidisk vanishes, and the amplitude of variability is magnified, being dependent only on the rate at which material enters the accretion zone. For slower sink accretion, where the minidisk persists over many orbits, the minidisk acts as a “buffer," smoothing out the “spiky" events of material being added to the minidisk (see Duffell et al. (2020) for one investigation of this relationship between accretion variability and sink accretion timescale). Here, we have chosen to model the minidisk as a steady $\alpha$-disk using the same model as for our initial conditions, which resulted in a stable minidisk with moderate accretion variability. Other disk models, such as an eccentric or shock-dominated disk, may be appropriate to consider for modeling the sink accretion timescale, but will likely impact the variability seen here primarily by being either faster or slower than the timescale used in this work, resulting in variability which is more or less pronounced, respectively.

Shocks have not been used as a signifier of variable EM emission in preceding isothermal works, though they have been shown to be a source of periodic X-rays in non-isothermal simulations such as Tang et al. (2018), implying that shock capturing and characterization is necessary to understand a potentially crucial source of periodic emission driven by the binary. We find from our shock monitoring scheme that, like Tang et al. (2018), shocks occur periodically from rejected accretion streams striking the gap wall. This periodic shocking, like our variable accretion, occurs on a timescale equal to the orbital timescale of the binary.

One important caveat to the variability seen in our accretion and shocks is that the amplitude is modest, with both proxies peaking $\lesssim 10\%$ above their mean. However, the timescale is well-constrained and short ($t_{\rm{bin}}\approx 0.1\;\rm{yr}$), allowing many cycles of such a system to be obtained from continued monitoring with UV or X-ray instruments, alleviating the difficulties discussed by Vaughan et al. (2016) in distinguishing genuine periodicity from insufficiently-sampled red noise.

Further still, in Fig. 8 we examine the timing relationship between the shock lightcurve and the secondary’s accretion rate via cross-correlation and find that they are nearly fully out-of-phase with one another, with a timing offset of $\Delta t_{\rm{lag}}\approx 0.43\;t_{\rm{bin}}$ between peaks in the accretion rate and peaks in the shock luminosity. This parallels the relationship found between the optical lightcurve and accretion rate for a circular binary in the fully adiabatic simulations performed by Westernacher-Schneider et al. (2022). This timing offset is potentially valuable as an identifier of SMBHB candidates, as the shocks can produce X-ray variability (Tang et al., 2018), while the minidisk emission will typically peak at UV or longer wavelengths, depending on the mass of the secondary (Shakura & Sunyaev, 1973). This type of correlated variability across wavebands is not expected for red noise, and thus cross-correlation analysis between UV and X-ray monitoring of AGN could be a valuable observational pursuit for identifying candidate SMBHBs.

A recurring point of discussion around SMBHBs is the sign of the gravitational torque acting on the binary, as this tends to dominate the evolution of the separation prior to the gravitational wave-dominated phase. We find that our disk exerts a negative torque on the binary, reducing the binary separation over time, which matches Duffell et al. (2020), which finds that the torque is negative for mass ratios below $q\sim 0.05$, albeit for a much warmer disk with constant $H/R=0.1$. The majority of the negative torque can be attributed to the build up of material at the outer edge of an eccentric gap – as seen in the case of equal-mass binaries in Tiede et al. (2020) in which this effect is greater in magnitude for colder disks that accumulate more material at the gap edge. Conversely, Derdzinski et al. (2021) find, for lower mass ratios than ours, that the torque on the binary flips from negative to positive as the disk becomes thinner, though they note that the density asymmetry key to determining the gravitational torque is unresolved, falling within the smoothing length of the secondary in some of their simulations. In general, these works and ours suggest that the magnitude of this effect is dependent on the specific disk model as well as the mass ratio of the embedded binary, likely due to whether the system leads to the opening of an eccentric gap. A systematic study of the gravitational torque over a wide range of gap opening scenarios is an interesting avenue for future investigation.

While not explored in this work, it is expected that, just as the binary excites eccentricity in the disk, the disk should excite eccentricity in the binary orbit. There are several works which have explored the case of eccentric binaries (D’Orazio & Duffell, 2021; Miranda et al., 2017) finding that the eccentricity of the binary affects the variability of accretion, the morphology and dynamics of the disk gap, and the evolution of the binary orbit. Franchini et al. (2023) have also found that fixing the binary orbit leads to overestimation of the gravitational torque and underestimation of accretion torques for equal-mass binaries. Exploring the effects of eccentric and live binaries in our cold disk model and how these compare to the existing warmer disk works is an interesting path for future investigation.