Shock-driven periodic variability in a low-mass-ratio supermassive black hole binary

ABSTRACT

We investigate the time-varying electromagnetic emission of a low-mass-ratio supermassive black hole binary (SMBHB) embedded in a circumprimary disc, with a particular interest in variability of shocks driven by the binary. We perform a 2D, locally isothermal hydrodynamics simulation of an SMBHB with mass ratio q = 0.01 and separation a = 100 R_g, using a physically self-consistent steady disc model. We estimate the electromagnetic variability from the system by monitoring accretion on to the secondary and using an artificial viscosity scheme to capture shocks and monitor the energy dissipated. The SMBHB produces a wide, eccentric gap in the disc, previously only observed for larger mass ratios, which we attribute to our disc model being much thinner (H/R ≈ 0.01 near the secondary) than is typical of previous works. The eccentric gap drives periodic accretion on to the secondary SMBH on a time-scale matching the orbital period of the binary, |$t_{\rm {bin}}\approx 0.1\,\,\rm {yr}$|⁠, implying that the variable accretion regime of the SMBHB parameter space extends to lower mass ratios than previously established. Shocks driven by the binary are periodic, with a period matching the orbital period, and the shocks are correlated with the accretion rate, with peaks in the shock luminosity lagging peaks in the accretion rate by 0.43 t_bin. We propose that the correlation of these quantities represents a useful identifier of SMBHB candidates, via observations of correlated variability in X-ray and ultraviolet monitoring of active galactic nuclei, rather than single-waveband periodicity alone.

accretion, accretion discs, black hole physics, hydrodynamics, shock waves

1 INTRODUCTION

Since most galaxies are host to a supermassive black hole (SMBH; Kormendy & Ho 2013) and galaxy mergers are a regular occurrence (Lotz et al. 2011), it is expected that there should exist a population of post-merger galaxies containing a supermassive black hole binary (SMBHB; Begelman, Blandford & Rees 1980). These binaries are a necessary precursor to SMBH mergers, and both binaries and mergers are sources of gravitational waves (GWs) which may be detected by pulsar timing array (PTA) collaborations and the upcoming Laser Interferometer Space Antenna (LISA) mission (Amaro-Seoane et al. 2017), respectively. While there is no conclusive detection yet for a continuous-wave GW signal of an individual binary (Arzoumanian et al. 2023; Antoniadis et al. 2023a; Agazie et al. 2023c), several PTA collaborations have recently reported detection of a nanohertz gravitational wave background (GWB), consistent with what would be produced by a population of SMBHBs (EPTA Collaboration 2023; Xu et al. 2023; Agazie et al. 2023a, 2023b; Reardon et al. 2023a; Antoniadis et al. 2023b; Reardon et al. 2023b).

These discoveries and future continuous-wave GW detections of individual binaries can best be used for astrophysical inference when combined with electromagnetic (EM) observations, necessitating the development of methods for identifying and studying candidate SMBHBs with EM data alone. At wider separations, candidate multiple-SMBH systems can be identified through spatially resolving them (down to ∼1 kpc; Komossa et al. 2003; Comerford et al. 2012; Blecha, Loeb & Narayan 2013; Foord et al. 2019), through H i absorption lines with different Doppler shifts implying the presence of a second relativistic jet (down to approximately a few pc; Rodriguez et al. 2009; Tremblay et al. 2016; Bansal et al. 2017), or through identifying spectrally distinct broad line emission regions (down to ∼0.1 pc; Eracleous et al. 2012; Runnoe et al. 2017). However, at the very close (≪1 pc) separations most relevant to multimessenger studies with GWs, we rely on searches for active galactic nuclei (AGNs) with periodically modulated light curves (Graham et al. 2015; Liu et al. 2019; Chen et al. 2020; Liu et al. 2020, among others). Unfortunately, even single-SMBH AGNs are intrinsically variable sources, such that distinguishing periodicity from stochastic variability requires observing many cycles of the hypothetical period (Vaughan et al. 2016). It is important, then, to understand the properties of binary-induced EM periodicity which distinguish it from other sources of AGNs variability.

A great deal of numerical simulations have been performed to understand the nature of SMBHB variability and evolution, investigating trends with binary mass ratio, separation, and eccentricity, as well as with disc temperature and viscosity (MacFadyen & Milosavljević 2008; D’Orazio, Haiman & MacFadyen 2013; Farris et al. 2014, 2015; D’Orazio et al. 2016; Tang, MacFadyen & Haiman 2017; Tang, Haiman & MacFadyen 2018; Moody, Shi & Stone 2019; Duffell et al. 2020; Tiede et al. 2020; D’Orazio & Duffell 2021; Derdzinski et al. 2021; Westernacher-Schneider et al. 2022).

The vast majority of the existing literature has used a locally isothermal equation of state, with accretion rate on to the binary components being used as the proxy for EM variability. However, Tang, Haiman & MacFadyen (2018) have shown in a non-isothermal simulation that shock-heating of the cavity wall produces significant periodicity at high energies. This shock-induced periodicity can still be monitored in locally isothermal simulations via shock capturing, as binary-induced shocks are still present and dissipating energy, offering a second proxy for EM variability to analyse alongside accretion rates.

A particularly interesting case to examine is low-q SMBHBs. The majority of galaxy mergers are unequal-mass, and thus are expected to produce unequal-mass SMBHBs, with the caveat that smaller BHs are likely to form binaries less efficiently. Additionally, there are other proposed formation channels for low-q SMBHBs, such as the presence of low-mass (10²–10⁴ M_⊙) SMBH seeds (Bellovary et al. 2011), in situ growth of stellar mass BHs to intermediate mass black hole (IMBH) masses in the disc (McKernan et al. 2012), and disruption of IMBH-hosting globular clusters (Fragione & Capuzzo-Dolcetta 2016; Arca-Sedda & Gualandris 2018; Fragione, Ginsburg & Kocsis 2018a; Fragione et al. 2018b; Rasskazov et al. 2019). However, it has been found previously that accretion on to the binary components becomes steady for q < 0.05 (D’Orazio et al. 2016). It is of interest, then, to determine whether periodic shocks also vanish for low-q SMBHBs, or whether the parameter space for EM-variable SMBHBs extends to lower mass ratios than previously established.

In this work, we perform a 2D hydrodynamic simulation of a q = 0.01 SMBHB and examine its implied variable electromagnetic properties via monitoring of both accretion and shocks. Shock capturing is performed using an artificial viscosity to both identify shocks and calculate the energy dissipated by them (Von Neumann & Richtmyer 1950). Further, we relax the common assumption of constant disc aspect ratio in favour of a physically self-consistent disc model in order to better reproduce the real environments of these systems.

This paper is organized as follows. In Section 2, we describe the simulation set-up, including the models used for the disc, gravity, accretion, and shock capturing. Section 3 presents our analysis of the disc dynamics and morphological evolution and our monitored accretion rate and shock outputs. We discuss the implications of these results in Section 4, in particular the dynamical and morphological differences from past works at this mass ratio (Section 4.1) and EM observables (Section 4.2). Section 5 summarizes our findings. Appendix A provides a derivation for the disc model used in this work.

2 METHODS

In this section, we describe the set-up of the numerical simulation performed for this work. We review the initial conditions used for the system, the physics included in the simulation, and methodology for monitoring accretion and shock dissipation, which will serve as proxies for electromagnetic variability in our analysis.

2.1 FARGO3D

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

$$\begin{eqnarray} \frac{\partial \Sigma }{\partial t} + \nabla \cdot \left(\Sigma {\bf v}\right) = 0, \end{eqnarray}$$

(1)

$$\begin{eqnarray} \Sigma \left(\frac{\partial {\bf v}}{\partial t} + {\bf v} \cdot \nabla {\bf v}\right) = -\nabla P-\Sigma \nabla \Phi + \nabla \cdot \Pi , \end{eqnarray}$$

(2)

where Σ is the surface density, v is the fluid velocity, |$P=\Sigma c_\mathrm{ s}^2$| is the isothermal pressure with c_s the isothermal sound speed, Φ is the gravitational potential, and Π is the viscous stress tensor. We use fargo3d’s built-in α-prescription to set the viscosity as ν = αc_sH (Shakura & Sunyaev 1973), where H is the local scale height of the disc.

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 time-step for systems dominated by rotation, such as accretion discs.

2.2 Disc set-up

Our simulated domain is that of a disc centred on a primary SMBH with mass M₁ = 10⁸ M_⊙ and extends from R_min = 10 R_g to R_max = 1000 R_g, where R_g = GM₁/c² is the gravitational radius of the primary SMBH. The secondary SMBH has mass M₂ = 10⁶ M_⊙ and is placed at a separation a = 100 R_g from the primary, resulting in a binary orbital time-scale |$t_{\rm {bin}}\approx 0.1\,\rm {yr}$|⁠. We use the torque-free boundary conditions described by Dempsey, Lee & Lithwick (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 disc are derived from the steady disc equations (Frank, King & Raine 2002),

$$\begin{eqnarray} \left. \begin{aligned} &\rho = \Sigma /H;\\ &H = c_\mathrm{ s} R^{3/2}/\left(GM\right)^{1/2};\\ &c_\mathrm{ s}^2 = P/\rho ;\\ &P = \frac{\rho k_\mathrm{ B} T_\mathrm{ c}}{\mu m_\mathrm{ p}} + \frac{4\sigma }{3c} T_\mathrm{ c}^4;\\ &\frac{4\sigma T_\mathrm{ c}^4}{3\tau } = \frac{3GM\dot{M}}{8\pi R^3}\left[1-\left(\frac{2R_\mathrm{ g}}{R}\right)^{1/2}\right];\\ &\tau = \Sigma \kappa _{\rm R}\left(\rho ,T_\mathrm{ c}\right)=\tau \left(\Sigma ,\rho ,T_\mathrm{ c}\right);\\ & \nu \Sigma = \frac{\dot{M}}{3\pi } \left[1-\left(\frac{2R_\mathrm{ g}}{R}\right)^{1/2}\right];\\ &\nu = \nu \left(\rho ,T_\mathrm{ c},\Sigma ,\alpha ,...\right), \end{aligned} \right\rbrace , \end{eqnarray}$$

(3)

where ρ is the gas density, Σ 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 disc, P is the total pressure, T_c is the temperature at the disc midplane, ν is the viscosity, τ is the optical depth, and κ_R is the Rosseland mean opacity, which throughout this work is assumed to be dominated by electron scattering, κ_R ≈ κ_e = 0.4 cm² g⁻¹. We choose an α-viscosity, ν = αc_sH (Shakura & Sunyaev 1973), and introduce a prescription for the scale height in order to control where the disc transitions from radiation-dominated to gas-dominated and the proportionality of H to r in the gas-dominated region:

$$\begin{eqnarray} H(r) = h(r) d(r)^{1/2}; \end{eqnarray}$$

(4)

$$\begin{eqnarray} h(r) \equiv \frac{3}{2} \frac{GM_\odot }{\eta c^2} m\dot{m}\left[1-\left(\frac{2R_\mathrm{ g}}{R}\right)^{1/2}\right], \end{eqnarray}$$

(5)

$$\begin{eqnarray} d(r) \equiv \left[\left(\frac{r}{r_\mathrm{ c}}\right)^{2k} + \frac{1}{2}\left(1+\sqrt{1+4\left(\frac{r}{r_\mathrm{ c}}\right)^{2k}}\right)\right]. \end{eqnarray}$$

(6)

Here, m ≡ M/M_⊙ is the mass of the central object in units of solar masses, r ≡ R/R_g is the radial coordinate in units of gravitational radii (R_g ≡ GM/c²), 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_\mathrm{ p}/(\eta \sigma _\mathrm{ T} c)$| with m_p being the mass of a proton, σ_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 disc pressure transitions from radiation to gas-dominated and the power-law slope of H in the gas-dominated region (i.e. H ∝ r^k far from the central object). With these prescriptions for ν and H in hand, we are able to solve the system in equation (3) for the surface density and sound speed distributions,

$$\begin{eqnarray} \Sigma (r) = \frac{16}{27}\frac{m_\mathrm{ p}}{\sigma _\mathrm{ T}} \eta \alpha ^{-1} \dot{m}^{-1} r^{{3}/{2}} \left(1-\sqrt{\frac{2}{r}}\right)^{-1} d\left(r\right)^{-1}, \end{eqnarray}$$

(7)

$$\begin{eqnarray} c_\mathrm{ s}(r) = \frac{3}{2}c\eta ^{-1}\dot{m}r^{-{3}/{2}}\left(1-\sqrt{\frac{2}{r}}\right)^2 d\left(r\right)^{{1}/{2}}. \end{eqnarray}$$

(8)

A full derivation of these expressions [equations (4)–(8)] can be found in Appendix A. For this work, we use α = 10⁻², |$\dot{m}=10^{-1}$|⁠, η = 10⁻¹, 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_gas + P_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 initial conditions used for the simulation, calculated from the disc model described in Section 2.2. The vertical dotted line on each panel shows the position of the secondary SMBH. For H/R (solid purple), Hr/R (dashed orange) and Hg/R (dotted magenta), where $H_\mathrm{ r}^2\equiv P_\mathrm{ r}/(\Sigma \Omega _\mathrm{ K}^2)$ and $H_\mathrm{ g}^2\equiv P_\mathrm{ g}/(\Sigma \Omega _\mathrm{ K}^2)$, are also plotted to illustrate the relative contributions of radiation and gas pressure, respectively, throughout the disc.$

Figure 1.

The initial conditions used for the simulation, calculated from the disc model described in Section 2.2. The vertical dotted line on each panel shows the position of the secondary SMBH. For H/R (solid purple), H_r/R (dashed orange) and H_g/R (dotted magenta), where |$H_\mathrm{ r}^2\equiv P_\mathrm{ r}/(\Sigma \Omega _\mathrm{ K}^2)$| and |$H_\mathrm{ g}^2\equiv P_\mathrm{ g}/(\Sigma \Omega _\mathrm{ K}^2)$|⁠, are also plotted to illustrate the relative contributions of radiation and gas pressure, respectively, throughout the disc.

Open in new tab Download slide

The simulation is run on a cylindrical grid of N_r = 2560, N_θ = 3493 cells. The radial spacing of the grid is logarithmic, with Δ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 disc (Teyssandier & Ogilvie 2017).

2.3 Gravity

The regions of interest to this work, namely the secondary’s minidisc 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 disc scale height, giving each black hole a potential of the form

$$\begin{eqnarray} \Phi _\mathrm{ i} = \frac{GM_\mathrm{ i}}{\sqrt{r_\mathrm{ i}^2 + \xi ^2 H^2}}, \end{eqnarray}$$

(9)

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/Ω_K is the disc scale height at that point, with Ω_K the local Keplerian angular velocity, and ξ is a constant parameter controlling the magnitude of the softening. We choose ξ = 0.6, following from considerations made by Duffell & MacFadyen (2013). The calculations do not include the self-gravity of the disc, as the disc is much less massive than the secondary SMBH (M_disc/M₂ ∼ 0.04). The mass of the secondary is increased from 0 to M₂ 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 disc 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 Section 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),

$$\begin{eqnarray} \frac{t_{\mathrm{sim}}}{t_{\mathrm{GW}}} = \frac{512 \pi }{5} N_{\mathrm{orb}} \tilde{r}^{-5/2} q (1 + q)^{1/2} (1-f_\mathrm{sep}^4)^{-1}, \end{eqnarray}$$

(10)

where N_orb = 10³ is the number of orbits in the simulation, |$\tilde{r}=10^2$| is the initial separation in units of R_g, q = 10⁻² is the mass ratio, and f_sep is the fraction of the initial separation to which the orbit decays. For a small decay to f_sep = 0.99 and our binary parameters, this ratio is below unity, indicating that the orbit will decay by |$\lt 1~{{\rm per\ cent}}$| of the initial separation over the course of our simulation.

2.4 Accretion model

The purpose of including accretion of gas on to the secondary SMBH is twofold: first, the removal of gas from the minidisc is a natural part of the disc accretion flow, preventing spurious buildup of material at the location of the secondary; second, the mass accretion rate of an α-disc is directly linked to its radiative emission, so monitoring this property allows an estimation of the minidisc’s luminosity over time. The secondary removes gas from within an accretion radius r_accr = 0.5r_hill, where r_hill = a(q/3)^1/3 is the Hill radius of the secondary, on an accretion time-scale, t_accr. t_accr varies as a function of distance r from the secondary SMBH, based on the viscous time-scale for our disc model,

$$\begin{eqnarray} t_{\rm {accr}} = \frac{4}{9} \frac{GM_\odot }{c^3} \eta ^2 \alpha ^{-1} m \dot{m}^{-2} r^{7/2} \left(1-\sqrt{\frac{2}{r}}\right)^2 d\left(r\right), \end{eqnarray}$$

(11)

where we use the same values as in our initial conditions for all parameters except that m = 10⁶ and r is now measured in gravitational radii of the secondary. A disc-like accretion time-scale such as this steeply increases with distance from the central body, causing accretion to be dominated by the inner parts of the minidisc and relatively insensitive to the specific value of r_accr chosen.

The fraction of mass removed from cells within r_accr is computed as

$$\begin{eqnarray} f_{\rm {accr}} = \frac{\Delta t}{t_{\rm {accr}}}, \end{eqnarray}$$

(12)

where Δt is the simulation time-step. We monitor mass, momentum, and angular momentum accreted on to the secondary, binned over every 10⁻³ t_bin.

2.5 Shock capturing

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

$$\begin{eqnarray} {\bf Q} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}l^2\Delta x^2 \Sigma \nabla \cdot {\bf v}\left[\nabla {\bf v}-\frac{1}{3}\left(\nabla \cdot {\bf v}\right){\bf e}\right],& \text{if } \nabla \cdot {\bf v}\lt 0\\ 0,& \text{otherwise}, \end{array}\right. \end{eqnarray}$$

(13)

where Δx is the zone size, ∇v is the symmetrized velocity gradient tensor, 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 (∇ · 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 disc and does not significantly contribute to heating the gas on long enough time-scales to affect the disc evolution.

For the range of typical shock Mach numbers in the simulation (⁠|$\mathcal {M}\gtrsim 10$|⁠), the thermal time-scale due to radiative diffusion t_diff = H²(3C_p/16σ)(ρ²κ_R/T³) (Kippenhahn & Weigert 1990), with C_p the specific heat at constant pressure, is less than t_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⁻³ t_bin. The Hill sphere of the secondary is excluded from the analysis of this monitoring.

3 RESULTS

In this section, we present outputs from the simulation as well as our analysis of the dynamics and of the accretion rate and shock dissipation, which stand in as predictors of the electromagnetic variability of the system. When indicated, we restrict our analysis to simulation outputs after the system has reached a quasi-steady state, when t > 700 t_bin, the angular momentum deficit (AMD) has saturated, and the r > a shock luminosity and average accretion rate on to the secondary have plateaued.

3.1 Disc morphology and dynamics

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 disc. In Fig. 2, we show snapshots of the disc surface density at various times and overplot ellipses showing the orbits taken by gas in the disc. 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 semimajor axis a, eccentricity e, and argument of pericentre |$\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 δa. Throughout this analysis, we use a bin width of δa = 1.0 R_g.

Figure 2.

Snapshots of the gas surface density over time. The white circle indicates the extent of the secondary SMBH’s Hill sphere, and the black ellipses show fits to the orbits taken by gas in the disc, with the circles on each ellipse indicating pericentre for that orbit. The snapshot at t = 1000 t_bin includes a 50 R_g × 50 R_g inset showing the minidisc and accretion streams, with the colourbar rescaled to span 3–3 × 10³ g cm⁻² to better highlight the structure in this region. The snapshot at t = 450 t_bin has dashed lines showing the boundaries used in the torque analysis in Section 3.1. Early on, the gap opened by the secondary is circular, but interactions with the binary subsequently drive rapid eccentricity growth, which saturates between 400–500 t_bin. The orientation of this eccentric gap precesses relative to the binary over long time-scales.

Open in new tab Download slide

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

$$\begin{eqnarray} \mathrm{AMD} = \int ^{R_\mathrm{out}}_{R_\mathrm{in}} \int _0^{2\pi } \Sigma R^2 \Omega \left(1-\sqrt{1 - e^2}\right)R\,\,\mathrm{ d}R\,\,\mathrm{ d}\phi , \end{eqnarray}$$

(14)

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 behaviour in other monitored quantities (see Section 3.2), motivating our use of it as an indicator of reaching a quasi-steady state.

Figure 3.

Total angular momentum deficit between r = 200 R_g and r = 400 R_g. The orange box highlights the time during which the AMD grows exponentially. Linear fits to the growth rate are plotted for two separate times, one during the exponential phase and one immediately after as the AMD begins to saturate. These fits are labelled with the corresponding eccentricity growth rate, calculated by multiplying the slope of each fit by log (10)/(2π × 2), which corrects for the log scaling, the time unit, and the proportionality AMD ∝ e², as was done in Teyssandier & Ogilvie (2017).

Open in new tab Download slide

Fig. 4 shows the time-averaged eccentricity profile of our disc at the end of the simulation, both in terms of the radial coordinate of the simulation, r, and the semimajor axes, a, of the gas orbits. The gap edge, at approximately r, a = 200 R_g, reaches an eccentricity of 0.25, larger even than the e ≈ 0.1–0.2 measured for higher mass ratios in past works with hotter discs (MacFadyen & Milosavljević 2008; Farris et al. 2014).

Figure 4.

Eccentricity profile of the disc as a function of radial coordinate r (solid purple) and semimajor axis a (dashed orange). These profiles were time-averaged over 50 t_bin from t = 950 t_bin to t = 1000 t_bin. The peak in eccentricity occurs at the approximate location of the gap edge.

Open in new tab Download slide

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_g ellipse, corresponding to the gap edge, over time. From this, we obtain a precession time-scale of ∼1.41 × 10³ t_bin, comparable to the ∼3.85 × 10³ t_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

$$\begin{eqnarray} \frac{\mathrm{ d}T}{\mathrm{ d}r} = -2 \pi R \left\langle \Sigma \frac{\partial \Phi }{\partial \phi }\right\rangle , \end{eqnarray}$$

(15)

and integrated torque, |$T(r)=\int ^{r}_{0}\left({\mathrm{ d}T}/{\mathrm{ d}r}\right)\mathrm{ d}r$|⁠, of the secondary SMBH acting on the disc, averaged over 50 t_bin. The strongest peak in the torque density occurs near r = 2.5a, coinciding with the peak in the surface density distribution of the disc. Beyond this, the torque density quickly decays to oscillations about zero. The total torque experienced by the disc is positive, with T(∞) ≈ 7.89 × 10⁵⁰ g cm²s⁻². The reciprocal torque of the disc on the secondary is thus negative, serving to shrink the binary separation on a time-scale t_T = M₂a²Ω_BH2/T(∞) = 3.62 × 10⁶ t_bin. This is substantially longer than the runtime of the simulation, but much shorter than the time-scales 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 disc at this time.

Figure 5.

Torque density (thin purple) and integrated torque (thick orange) exerted by the binary on the disc. Quantities were time-averaged over 50 t_bin from t = 950 t_bin to t = 1000 t_bin. The total torque on the disc is positive, and thus the reciprocal torque on the binary is negative, shrinking the binary separation over time. The disc can be separated into an inner, gap, and outer region based on the distribution of the torque density, indicated with the vertical grey lines at 0.7a, 1.4a, and 3a. The total torque is plainly dominated by the contributions of the material near the outer disc edge.

Open in new tab Download slide

We can further explore where this negative torque comes from by following a similar procedure to Tiede et al. (2020) and Muñoz, Miranda & Lai (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 disc, (2) an outer region R > 1.4a, which includes the entire outer disc, 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 disc, with |T_out| ≈ 9|T_in|, and, further, that these contributions become negligible past R ∼ 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 discs drives the total torque on the binary to be negative.

3.2 Electromagnetic emission

The primary aim of this work is to investigate the time-varying electromagnetic properties of a low-mass-ratio SMBHB embedded in a realistic disc. 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 (Section 2.5), and accretion rate on to the secondary (Section 2.4), standing in for the radiative emission of the minidisc. 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 ∼1.2.

3.2.1 Shocks

There are two major modes in which shocks occur in our disc: spiral shocks in the inner disc and shocks at the outer edge of the eccentric gap, driven by a portion of the gas streams feeding the secondary being flung back into the gap edge. The behaviour of shocks in these regions can be observed in Fig. 6, which shows the surface density and shock dissipation flux as a functions of time and radius through the disc. We see that, in the quasi-steady state, the shock dissipation in the inner disc is steady, while the outer gap edge displays clear repeating shock waves, which map on to similar behaviour in the surface density. It is these shock waves in the outer gap edge that drive variability in the total shock emission.

$Azimuthally averaged surface density $\bar{\Sigma }$ (top) and azimuthally summed shock flux Fshock (bottom) as functions of time. The left column shows the quantities over the full simulation, while the right zooms in on 10 orbits during the quasi-steady state to highlight the parallel periodic structure in both quantities. The white bar on the left column plots indicates the times covered by the right column plots. The inner disc behaviour is steady on short time-scales, while a regular pattern is seen in both quantities at the outer edge of the gap (∼150–250 Rg), indicating periodic shock waves excited by the binary. Shock dissipation in the inner disc gradually decreases over the course of the simulation, caused by depletion of surface density in the inner disc as gas is accreted across the inner boundary of the domain.$

Figure 6.

Azimuthally averaged surface density |$\bar{\Sigma }$| (top) and azimuthally summed shock flux F_shock (bottom) as functions of time. The left column shows the quantities over the full simulation, while the right zooms in on 10 orbits during the quasi-steady state to highlight the parallel periodic structure in both quantities. The white bar on the left column plots indicates the times covered by the right column plots. The inner disc behaviour is steady on short time-scales, while a regular pattern is seen in both quantities at the outer edge of the gap (∼150–250 R_g), indicating periodic shock waves excited by the binary. Shock dissipation in the inner disc gradually decreases over the course of the simulation, caused by depletion of surface density in the inner disc as gas is accreted across the inner boundary of the domain.

Open in new tab Download slide

We produce a shock ‘light curve’ (Fig. 7, top left panel) by monitoring the total energy dissipated by artificial viscosity across the disc over our output time-scale DT = 10⁻³ t_bin. First, we note that the shocks are bright, producing an average luminosity |$L_{\rm {shock}}\approx 0.13\,\,\rm {L}_{\rm {Edd}}$| late in the simulation, with |$\rm {L}_{\rm {Edd}}$| being the Eddington luminosity of the primary. For comparison, the unperturbed disc of our initial conditions has a luminosity |$L_{\rm {disc}}\approx 0.25\,\,\rm {L}_{\rm {Edd}}$|⁠, assuming an opacity dominated by electron scattering.

Figure 7.

Top: Shock dissipation rate summed over the disc (left) and accretion rate onto the secondary SMBH (right), shown over the full runtime of the simulation. The insets zoom in on 10 orbits late in the run to highlight the periodic structure of each quantity. The thin orange curves are moving averages of each quantity, taken over 10 orbits. For the shock light curve, the horizontal dotted line shows, for comparison, the thermal luminosity of the unperturbed disc model used for our initial conditions, though we note that the shocks will primarily emit in X-rays, while the thermal luminosity will peak at optical wavelengths. The pink and grey curves are, respectively, the shock luminosity from outside and inside the secondary’s orbit, revealing that the variability is dominated by regions outside the secondary’s orbit (see also Fig. 6) and that the behaviour of shocks in that region have reached a steady state. The thin lilac curve is a moving-average of the r > a shock luminosity, taken over 10 orbits. The large, long-duration flare in the shock luminosity between 300 t_bin and 400 t_bin coincides with the exponential growth in gap eccentricity illustrated in Fig. 3. Bottom: Fourier transforms of the shock light curve (left) and accretion rate (right). Only data from after the system reached a quasi-steady state around ∼700 t_bin were included in these calculations. Both quantities exhibit clear periodicity on the orbital time-scale of the binary.

Open in new tab Download slide

The bottom-left panel of Fig. 7 shows the Fourier transform of the light curve. There is a clear peak corresponding to the orbital frequency of the binary, with the first few harmonics of this frequency also present, but weak. This periodicity is plainly visible in the zoomed-in inset of the light curve itself, where we also note that the peak-to-trough ratio of the shock luminosity is only 1.17, though this rises to ∼1.4 when considering emission from r > a only, as may be relevant if the r < a emission continues decaying to zero on longer time-scales.

3.2.2 Accretion

The accretion rate on to the secondary, shown in the top-right panel of Fig. 7, is highly super-Eddington, hovering at around |$\sim 130\,\,\rm {\dot{M}}_{\rm {Edd}}$| in the steady state. Enhancement of accretion through the circumbinary disc and comparable accretion rates on to unequal-mass binary components has also been seen in preceding works such as Farris et al. (2014). Enhancement of accretion through the disc is also, in general, expected given the enhancement of Reynolds stresses due to the secondary’s perturbation of the flow.

The bottom-right panel of Fig. 7 shows the Fourier transform of the secondary’s accretion rate. Like the shock light curve, the accretion rate varies on the orbital frequency of the binary, with clear but weaker harmonics also apparent. The accretion rate is less steady from orbit to orbit than the shock light curve, but has a similar average peak-to-trough ratio of ∼1.2. The once-per-orbit peaks in accretion rate occur when the secondary makes its closest approach to the gap edge, stripping off material, some of which is accreted on to the minidisc, and the rest of which is flung back into the gap edge. Similar behaviour is discussed in e.g. D’Orazio, Haiman & MacFadyen (2013).

We also monitored the linear and angular momentum accreted by the secondary to consider the resulting evolution of the secondary’s orbit and spin, respectively. In the steady state, the secondary accreted momentum at a rate Δp/t_bin = 1.01 × 10⁻¹⁰p_BH2, where p_BH2 is the initial momentum of the secondary on its orbit. This net accretion torque is positive and therefore acts to increase the binary separation, but its magnitude is very low, validating our choice not to evolve the binary orbit based on the accreted momentum during the simulation. The accreted angular momentum is likewise inconsequential. Expressed as a rate of change of the dimensionless spin parameter, the secondary SMBH accretes angular momentum at a rate Δa/t_bin = 1.30 × 10⁻¹⁰.

4 DISCUSSION

In this section, we discuss the implications of our results in greater detail. First, we examine the disc dynamics and why the dynamics observed in this work differ from previous low-mass-ratio SMBHB simulations. Then, we discuss our proxies for electromagnetic variability and their implications for observational identification of SMBHB candidates. Finally, we look at the various torques acting on the binary and their implications for the binary orbital evolution.

4.1 Disc dynamics

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 disc model for the surface density and scale height, we show that such behaviour can be present in mass ratios as low as q ∼ 0.01.

At a similar mass-ratio regime, simulations of super-Jupiters embedded in protoplanetary discs also produce eccentric gap systems with comparable eccentricity distributions to those in Fig. 4 (Dempsey, Muñoz & Lithwick 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 disc being much thinner than the thinnest H/R tested in Dempsey, Muñoz & Lithwick (2021), our H/R being non-uniform throughout our domain, or some other difference in the two set-ups is an interesting question, but beyond the scope of the present investigation.

The greatest difference between the disc in this study and those of preceding works is the disc model used for our initial conditions. In particular, our disc aspect ratio H/R varies with radius and is, in general, much thinner than the constant H/R = 0.1 discs common in the literature, with the disc having H/R ≈ 10⁻² 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_tid ∝ (H/R)⁻³, and the viscous torque T_visc ∝ ν ∝ (H/R)², which work to open the gap and to fill it in, respectively. In a colder, thinner disc, 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, Morbidelli & Masset 2006; Duffell 2015, 2020). This phenomenon of wider gaps in colder discs has been observed in simulations of equal-mass binaries which test different values of H/R (Ragusa, Lodato & Price 2016; Tiede et al. 2020), and here we demonstrate that this behaviour continues to lower mass ratios.

So a colder, thinner disc makes gap-clearing more effective. In clearing a deeper gap, the secondary removes material near its orbit, weakening the damping of eccentricity occurring at ECRs, enhancing the ability for the binary to drive eccentricity growth in the disc (Teyssandier & Ogilvie 2017). In summary, thinner discs are conducive to opening deeper gaps at lower mass ratios, which in turn drives eccentricity evolution in the outer disc and the consequent variable accretion on to the binary. Given that AGNs have discs with H/R ≈ 10⁻²–10⁻³ (Krolik 1999), consistent with the disc 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.

4.2 Variability

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 ≲ 0.05 (Farris et al. 2014; D’Orazio et al. 2016; Duffell et al. 2020). The time-scale of accretion variability depends on the mass ratio, with lower mass ratios 0.05 ≲ q ≲ 0.25 being modulated only on the binary orbital time-scale, while at larger mass ratios an overdense lump forms on the gap wall, causing periodic accretion on a time-scale linked to its own orbital time (D’Orazio, Haiman & MacFadyen 2013; Farris et al. 2014).

In contrast, we find periodic accretion onto our secondary even with q = 0.01, with a time-scale matching the orbital time-scale of the binary. This matches the properties seen for binaries in the 0.05 ≲ q ≲ 0.25 regime in these past works, and implies that this regime may simply extend to lower mass ratios given a thinner disc. 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 time-scale used determines how quickly the minidisc around the secondary is drained. If this time-scale is very short, the minidisc 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 minidisc persists over many orbits, the minidisc acts as a ‘buffer,’ smoothing out the ‘spiky’ events of material being added to the minidisc (see Duffell et al. (2020) for one investigation of this relationship between accretion variability and sink accretion time-scale). Here, we have chosen to model the minidisc as a steady α-disc using the same model as for our initial conditions, which resulted in a stable minidisc with moderate accretion variability. Other disc models, such as an eccentric or shock-dominated disc, may be appropriate to consider for modelling the sink accretion time-scale, but will likely impact the variability seen here primarily by being either faster or slower than the time-scale 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, Haiman & MacFadyen (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, Haiman & MacFadyen (2018), shocks occur periodically from rejected accretion streams striking the gap wall. This periodic shocking, like our variable accretion, occurs on a time-scale equal to the orbital time-scale 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~{{\rm per\ cent}}$| above their mean. However, the time-scale 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 ultraviolet (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 light curve 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 Δt_lag ≈ 0.43 t_bin between peaks in the accretion rate and peaks in the shock luminosity. This parallels the relationship found between the optical light curve 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, Haiman & MacFadyen 2018), while the minidisc 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 AGNs could be a valuable observational pursuit for identifying candidate SMBHBs.

Figure 8.

Examination of the timing relationship between shocks and accretion. Left: Shock light curve (thick purple) and accretion rate (thin orange), each normalized by their mean value. Right: Cross-correlation of the shock and accretion ‘light curves.’ The cross-correlation is itself periodic due to both signals being periodic, and peaks at a lag of 0.43 t_bin. Since the two signals have a period of t_bin, this affirms that they are almost fully out of phase with one another, as is apparent by eye in the left panel.

Open in new tab Download slide

4.3 Binary orbital evolution

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 disc 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 ∼ 0.05, albeit for a much warmer disc 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 discs 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 disc 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 disc 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 disc, the disc should excite eccentricity in the binary orbit. There are several works which have explored the case of eccentric binaries (Miranda, Muñoz & Lai 2017; D’Orazio & Duffell 2021) finding that the eccentricity of the binary affects the variability of accretion, the morphology and dynamics of the disc 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 disc model and how these compare to the existing warmer disc works is an interesting path for future investigation.

5 CONCLUSIONS

We performed a 2D, locally isothermal hydrodynamic simulation of a low-mass-ratio (q = 10⁻²) SMBHB embedded in a disc with initial surface density and sound speed profiles derived from a physically self-consistent disc model. We monitored the dynamical evolution of the system as well as two proxies for electromagnetic variability, the accretion rate on to the secondary SMBH, and the energy dissipated by shocks.

Our main findings can be summarized as follows.

The binary opens a wide, eccentric gap which precesses on long time-scales. This behaviour is seen throughout the literature at larger mass ratios, but has not generally been observed for q ≲ 0.05. We attribute this change to our disc model, which is much thinner than the H/R = 0.1 discs common in the existing literature, and in line with the H/R ∼ 0.01–0.001 expected for real AGNs discs. A thinner disc is both more susceptible to gravitational torques from the binary, which open the gap, and experiences weaker viscous torques, which serve to close it. Wider cavities in thinner discs have been seen in a few works which explored changing H/R with q = 1 binaries, and here we have shown that this behaviour extends to q = 0.01 binaries.
We find that accretion onto the secondary SMBH is clearly variable, with a period matching the orbital period of the binary and a peak-to-trough ratio of ∼1.2. As with the eccentric gap, this behaviour was previously not seen for q ≲ 0.05 in works using thicker discs. We find, as has been the case in previous works, that the peaks in the accretion rate occur due to the passage of the secondary near the overdensity at the edge of the gap, stripping off material to feed the minidisc. This process is necessarily linked to the gap becoming eccentric. Since we find that a thinner disc leads to eccentricity at smaller mass ratios than previously seen, it is unsurprising that accretion variability follows. Importantly, since real AGN discs are expected to be thin (H/R ∼ 0.01–0.001), we expect that the parameter space for which real binaries are variable is wider than has been previously established by works with thicker discs.
We find that shocks excited by the binary are also periodic, with a period matching the orbital period of the binary and a peak-to-trough ratio of 1.17. Shocks have been found to be an important source of periodic X-ray emission in non-isothermal works, but have not previously been monitored in isothermal simulations.
We find that there is a correlated lag between the accretion and shock light curves. The two quantities are nearly fully out of phase, with the shocks lagging behind accretion by 0.43 t_bin. This presents a potential smoking gun for binary candidacy in observations. Since accretion tracks the minidisc luminosity, which will typically be brightest at UV wavelengths, and shocks are seen to be bright in X-rays, the timing correlation between shock dissipation and accretion rate implies correlated variability in separate wavebands. Observations can then be made to search for such correlated variability as a sign of a possible SMBHB rather than more ambiguous single-waveband variability.

ACKNOWLEDGEMENTS

This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. This work was funded in part by Michigan Space Grant Consortium, National Aeronautics and Space Administration (NASA) grant #NNX15AJ20H.

AK would like to acknowledge support provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51463.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

MR acknowledges support from NASA grants 80NSSC20K1541 and 80NSSC20K1583 and National Science Foundation grants AST-1715140 and AST-2009227.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author

References

Agazie

et al. ,

2023a

ApJ

951

Month:	Total Views:
November 2023	8
December 2023	64
January 2024	41
February 2024	24
March 2024	17
April 2024	15
May 2024	26
June 2024	26
July 2024	8
August 2024	10
September 2024	26
October 2024	26
November 2024	5
December 2024	15
January 2025	17
February 2025	24
March 2025	13
April 2025	13
May 2025	4

Article Contents

Shock-driven periodic variability in a low-mass-ratio supermassive black hole binary

ABSTRACT

1 INTRODUCTION

2 METHODS

2.1 FARGO3D

2.2 Disc set-up

2.3 Gravity

2.4 Accretion model

2.5 Shock capturing

3 RESULTS

3.1 Disc morphology and dynamics

3.2 Electromagnetic emission

3.2.1 Shocks

3.2.2 Accretion

4 DISCUSSION

4.1 Disc dynamics

4.2 Variability

4.3 Binary orbital evolution

5 CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

References

APPENDIX: DISC MODEL

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only