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

We investigate the time-varying electromagnetic emission of a low-mass-ratio supermassive black hole binary (SMBHB) embedded in a circumprimary disk, with a particular interest in variability of shocks driven by the binary. We perform a 2D, locally isothermal hydrodynamics simulation of a SMBHB with mass ratio 𝑞 = 0 . 01 and separation 𝑎 = 100 𝑅 𝑔 , using a physically self-consistent steady disk model. We estimate the electromagnetic variability from the system by monitoring accretion onto 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 disk, previously only observed for larger mass ratios, which we attribute to our disk model being much thinner ( 𝐻 / 𝑅 ≈ 0 . 01 near the secondary) than is typical of previous works. The eccentric gap drives periodic accretion onto the secondary SMBH on a timescale matching the orbital period of the binary, 𝑡 bin ≈ 0 . 1 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 𝑡 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 UV monitoring of AGN, rather than single-waveband periodicity alone.


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 et al. 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 detection yet for a continuouswave GW signal of an individual binary (Arzoumanian et al. 2023), the NANOGrav collaboration has reported evidence for a stochastic common-process signal, which is expected for a gravitational wave background (GWB) generated by a population of SMBHBs (Middleton et al. 2016), though the current evidence remains insufficient to claim detection of a GWB (Arzoumanian et al. 2020).
These hints and future GW detections 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 ★ E-mail: kwhitle@umich.eduspatially resolving them (down to ∼1 kpc; Komossa et al. 2003;Comerford et al. 2012;Blecha et al. 2013;Foord et al. 2019), through HI absorption lines with different Doppler shifts implying the presence of a second relativistic jet (down to ∼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 (AGN) with periodically-modulated lightcurves (Graham et al. 2015;Liu et al. 2019Liu et al. , 2020;;Chen et al. 2020, among others).Unfortunately, even single-SMBH AGN 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 AGN variability.
The vast majority of the existing literature has used a locally isothermal equation of state, with accretion rate onto the binary components being used as the proxy for EM variability.However, Tang et al. (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 binaryinduced shocks are still present and dissipating energy, offering a second proxy for EM variability to analyze alongside accretion rates.
A particularly interesting case to examine is low- 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- SMBHBs, such as the presence of low-mass (10 2 -10 4 M ⊙ ) SMBH seeds (Bellovary et al. 2011), in-situ growth of stellar mass BHs to intermediate mass black hole (IMBH) masses in the disk (McKernan et al. 2012), and disruption of IMBH-hosting globular clusters (Fragione & Capuzzo-Dolcetta 2016;Fragione et al. 2018a;Arca-Sedda & Gualandris 2018;Fragione et al. 2018b;Rasskazov et al. 2019).However, it has been found previously that accretion onto the binary components becomes steady for  < 0.05 (D'Orazio et al. 2016).It is of interest, then, to determine whether periodic shocks also vanish for low- 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  = 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 disk aspect ratio in favor of a physically self-consistent disk 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 setup, including the models used for the disk, gravity, accretion, and shock capturing.Section 3 presents our analysis of the disk 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 ( §4.1) and EM observables ( §4.2).Section 5 summarizes our findings.Appendix A provides a derivation for the disk model used in this work.

METHODS
In this section, we describe the setup 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.

FARGO3D
The simulations run for this work were performed on a twodimensional, Eulerian, cylindrical grid using the FARGO3D code (Benítez-Llambay & Masset 2016).FARGO3D solves the hydrodynamics equations where Σ is the surface density, v is the fluid velocity,  = Σ 2  is the isothermal pressure with   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  =    (Shakura & Sunyaev 1973), where  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.

Disk Setup
Our simulated domain is that of a disk centered on a primary SMBH with mass  1 = 10 8  ⊙ and extends from  min = 10   to  max = 1000   , where   =   1 / 2 is the gravitational radius of the primary SMBH.The secondary SMBH has mass  2 = 10 6  ⊙ and is placed at a separation  = 100   from the primary, resulting in a binary orbital timescale  bin ≈ 0.1 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  is the gas density, Σ is the surface density,  is the scale height,   is the sound speed,  is the radial distance from the central object,  is the mass of the central object,  is the accretion rate through the disk,  is the total pressure,   is the temperature at the disk 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 2 /g.We choose an -viscosity,  =    (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  to  in the gas-dominated region: ℎ() ≡ 3 2 Here,  ≡ / ⊙ is the mass of the central object in units of solar masses,  ≡ /  is the radial coordinate in units of gravitational radii (  ≡  / 2 ), and  ≡ /  Edd is the accretion rate in units of Eddington, where  Edd ≡ 4   /(  ) with   being the mass of a proton,   the Thomson cross-section, and  ≡ /(  2 ) the accretion efficiency.The parameters   and  set, respectively, the radius at which the disk pressure transitions from radiation to gasdominated and the powerlaw slope of  in the gas-dominated region (i.e.,  ∝   far from the central object).With these prescriptions for  and  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  = 10 −2 ,  = 10 −1 ,  = 10 −1 ,   = 300, and  = 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  =  gas +  rad , the sum of the gas and radiation pressures.The resulting   () 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   = 2560,   = 3493 cells.The radial spacing of the grid is logarithmic, with Δ/ = 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).

Gravity
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   is the mass of black hole ,   is the distance between black hole  and the point of interest,  =   /Ω  is the disk scale height at that point, with Ω  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 disk, as the disk is much less massive than the secondary SMBH ( disk / 2 ∼ 0.04).The mass of the secondary is increased from 0 to  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  orb = 10 3 is the number of orbits in the simulation, r = 10 2 is the initial separation in units of   ,  = 10 −2 is the mass ratio, and  sep is the fraction of the initial separation to which the orbit decays.For a small decay to  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.

Accretion Model
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 -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  accr = 0.5 hill , where  hill = (/3) 1/3 is the Hill radius of the secondary, on an accretion timescale,  accr . accr varies as a function of distance  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  = 10 6 and  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  accr chosen.
The fraction of mass removed from cells within  accr is computed as where Δ is the simulation timestep.We monitor mass, momentum, and angular momentum accreted onto the secondary, binned over every 10 −3  bin .

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 where Δ is the zone size, ∇v is the symmetrized velocity gradient tensor, e is the unit tensor, and  is a parameter controlling the number of zone widths a shock is spread over.In this work, we use  = 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, / = −Q : ∇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 (M ≳ 10), the thermal timescale due to radiative diffusion the specific heat at constant pressure, is less than  bin in the region of interest  < 300   , supporting this assumption.
Like the accretion monitoring, we bin the shock emission over every 10 −3  bin .The Hill sphere of the secondary is excluded from the analysis of this monitoring.

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  > 700  bin , the angular momentum deficit (AMD) has saturated, and the  >  shock luminosity and average accretion rate onto the secondary have plateaued.

Disk Morphology and Dynamics
In contrast to previous works at this mass ratio, we find that a  = 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 ellipsefitting 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 , eccentricity , and argument of pericenter .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  and bin width .Throughout this analysis, we use a bin width of  = 1.0  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, , and the semi-major axes, , of the gas orbits.The gap edge, at approximately ,  = 200  g , reaches an eccentricity of 0.25, larger even than the  ≈ 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  = 200  g ellipse, corresponding to the gap edge, over time.From this, we obtain a precession timescale of ∼ 1.41 × 10 3  bin , comparable to the ∼ 3.85 × 10 3  bin calculated from the linear theory developed by Teyssandier & Ogilvie (2016).show fits to the orbits taken by gas in the disk, with the circles on each ellipse indicating pericenter for that orbit.The snapshot at  = 1000  bin includes a 50   × 50   inset showing the minidisk and accretion streams, with the colorbar rescaled to span 3-3 × 10 3 g/cm 2 to better highlight the structure in this region.The snapshot at  = 450  bin has dashed lines showing the boundaries used in the torque analysis in §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  bin .The orientation of this eccentric gap precesses relative to the binary over long timescales.Fig. 5 shows the time-averaged torque density, defined in Mac-Fadyen & Milosavljević (2008) as and integrated torque,  () = ∫  0 (/) , of the secondary SMBH acting on the disk, averaged over 50  bin .The strongest peak in the torque density occurs near  = 2.5, 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  (∞) ≈ 7.89 × 10 50 g cm 2 s −2 .The reciprocal torque of the disk on the secondary is thus negative, serving to shrink the binary separation on a timescale   =  2  2 Ω BH2 / (∞) = 3.62 × 10 6  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  < 0.7, consisting of the inner disk, (2) an outer region  > 1.4, which includes the entire outer disk, and (3) a gap region 0.7 <  < 1.4.It can be seen from the integrated torque in Fig. 5 that the torque is dominated by the contributions from the outer disk, with | out | ≈ 9| in |, and, further, that these contributions becomes negligible past  ∼ 3.We plot guides showing each of these regional boundaries in the bottom-left panel of Fig. 2, which reveals that the outer region 1.4 <  < 3 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.

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 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 ∼ 1.2.The total torque on the disk is positive, and thus the reciprocal torque on the binary is negative, shrinking the binary separation over time.The disk 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.7, 1.4, and 3.The total torque is plainly dominated by the contributions of the material near the outer disk edge.

Shocks
There are two major modes in which shocks occur in our disk: spiral shocks in the inner disk 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 behavior 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 disk.We see that, in the quasi-steady state, the shock dissipation in the inner disk is steady, while the outer gap edge displays clear repeating shock waves, which map onto similar behavior in the surface density.
It is these shock waves in the outer gap edge that drive variability in the total shock emission.We produce a shock "lightcurve" (Fig. 7, top left panel) by monitoring the total energy dissipated by artificial viscosity across the disk over our output timescale DT = 10 −3  bin .First, we note that the shocks are bright, producing an average luminosity  shock ≈ 0.13 L Edd late in the simulation, with L Edd being the Eddington luminosity of the primary.For comparison, the unperturbed disk of our initial conditions has a luminosity  disk ≈ 0.25 L Edd , assuming an opacity dominated by electron scattering.
The bottom-left panel of Fig. 7 shows the Fourier transform of the lightcurve.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 lightcurve 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  >  only, as may be relevant if the  <  emission continues decaying to zero on longer timescales.

Accretion
The accretion rate onto the secondary, shown in the top-right panel of Fig. 7, is highly super-Eddington, hovering at around ∼ 130 M Edd in the steady state.Enhancement of accretion through the circumbi- nary disk and comparable accretion rates onto unequal-mass binary components has also been seen in preceding works such as Farris et al. (2014).Enhancement of accretion through the disk 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 lightcurve, 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 lightcurve, 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 onto the minidisk, and the rest of which is flung back into the gap edge.Similar behavior is discussed in, e.g., D'Orazio et al. (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 Δ/ bin = 1.01 × 10 −10  BH2 , where  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 Δ/ bin = 1.30 × 10 −10 .

DISCUSSION
In this section, we discuss the implications of our results in greater detail.First, we examine the disk 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.The insets zoom in on 10 orbits late in the run to highlight the periodic structure of each quantity.The orange curves are moving averages of each quantity, taken over 10 orbits.For the shock lightcurve, the horizontal dotted line shows, for comparison, the thermal luminosity of the unperturbed disk 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 behavior of shocks in that region have reached a steady state.The lilac curve is a moving-average of the  >  shock luminosity, taken over 10 orbits.The large, long-duration flare in the shock luminosity between 300  bin and 400  bin coincides with the exponential growth in gap eccentricity illustrated in Fig. 3. Bottom: Fourier transforms of the shock lightcurve (left) and accretion rate (right).Only data from after the system reached a quasi-steady state around ∼ 700  bin were included in these calculations.Both quantities exhibit clear periodicity on the orbital timescale of the binary.

Disk Dynamics
Previous SMBHB works have found that eccentric cavities as drivers of accretion variability are not present at low mass ratios,  < 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  ∼ 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 / tested in Dempsey et al. (2021), our / 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 / varies with radius and is, in general, much thinner than the constant / = 0.1 disks common in the literature, with the disk having / ≈ 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,  tid ∝ (/) −3 , and the viscous torque  visc ∝  ∝ (/) 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 2015Duffell , 2020)).This phenomenon of wider gaps in colder disks has been observed in simulations of equal-mass binaries which test different values of / (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 / ≈ 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.

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  ≲ 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 ≲  ≲ 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  = 0.01, with a timescale matching the orbital timescale of the binary.This matches the properties seen for binaries in the 0.05 ≲  ≲ 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 -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 ≲ 10% above their mean.However, the timescale is well-constrained and short ( bin ≈ 0.1 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 crosscorrelation and find that they are nearly fully out-of-phase with one another, with a timing offset of Δ lag ≈ 0.43  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.

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 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  ∼ 0.05, albeit for a much warmer disk with constant / = 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.

CONCLUSIONS
We performed a 2D, locally isothermal hydrodynamic simulation of a low-mass-ratio ( = 10 −2 ) SMBHB embedded in a disk with initial surface density and sound speed profiles derived from a physically self-consistent disk model.We monitored the dynamical evolution of the system as well as two proxies for electromagnetic variability, the accretion rate onto the secondary SMBH, and the energy dissipated by shocks.
Our main findings can be summarized as follows.
(i) The binary opens a wide, eccentric gap which precesses on long timescales.This behavior is seen throughout the literature at larger mass ratios, but has not generally been observed for  ≲ 0.05.We attribute this change to our disk model, which is much thinner than the / = 0.1 disks common in the existing literature, and in line with the / ∼ 0.01-0.001expected for real AGN disks.A thinner disk 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 disks have been seen in a few works which explored changing / with  = 1 binaries, and here we have shown that this behavior extends to  = 0.01 binaries.
(ii) 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 behavior was previously not seen for  ≲ 0.05 in works using thicker disks.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 minidisk.This process is necessarily linked to the gap becoming eccentric.Since we find that a thinner disk leads to eccentricity at smaller mass ratios than previously seen, it is unsurprising that accretion variability follows.Importantly, since real AGN disks are expected to be thin (/ ∼ 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 disks.
(iii) We find that shocks excited by the binary are also periodic, with a period matching the orbital period of the binary and a peakto-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.
(iv) We find that there is a correlated lag between the accretion and shock lightcurves.The two quantities are nearly fully out of phase, with the shocks lagging behind accretion by 0.43  bin .This presents a potential smoking gun for binary candidacy in observations.Since accretion tracks the minidisk luminosity, which will typically be brightest at ultraviolet 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.

Figure 2 .
Figure2.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 disk, with the circles on each ellipse indicating pericenter for that orbit.The snapshot at  = 1000  bin includes a 50   × 50   inset showing the minidisk and accretion streams, with the colorbar rescaled to span 3-3 × 10 3 g/cm 2 to better highlight the structure in this region.The snapshot at  = 450  bin has dashed lines showing the boundaries used in the torque analysis in §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  bin .The orientation of this eccentric gap precesses relative to the binary over long timescales.

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

Figure 5 .
Figure5.Torque density (purple) and integrated torque (orange) exerted by the binary on the disk.Quantities were time-averaged over 50  bin from  = 950  bin to  = 1000  bin .The total torque on the disk is positive, and thus the reciprocal torque on the binary is negative, shrinking the binary separation over time.The disk 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.7, 1.4, and 3.The total torque is plainly dominated by the contributions of the material near the outer disk edge.

Figure 6 .
Figure 6.Azimuthally-averaged surface density Σ (top) and azimuthally-summed shock flux  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 disk behavior is steady on short timescales, while a regular pattern is seen in both quantities at the outer edge of the gap (∼ 150-250  g ), indicating periodic shock waves excited by the binary.Shock dissipation in the inner disk gradually decreases over the course of the simulation, caused by depletion of surface density in the inner disk as gas is accreted across the inner boundary of the domain.

Figure 7 .
Figure7.Top: Shock dissipation rate summed over the disk (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 orange curves are moving averages of each quantity, taken over 10 orbits.For the shock lightcurve, the horizontal dotted line shows, for comparison, the thermal luminosity of the unperturbed disk 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 behavior of shocks in that region have reached a steady state.The lilac curve is a moving-average of the  >  shock luminosity, taken over 10 orbits.The large, long-duration flare in the shock luminosity between 300  bin and 400  bin coincides with the exponential growth in gap eccentricity illustrated in Fig.3.Bottom: Fourier transforms of the shock lightcurve (left) and accretion rate (right).Only data from after the system reached a quasi-steady state around ∼ 700  bin were included in these calculations.Both quantities exhibit clear periodicity on the orbital timescale of the binary.

Figure 8 .
Figure 8. Examination of the timing relationship between shocks and accretion.Left: Shock lightcurve (purple) and accretion rate (orange), each normalized by their mean value.Right: Cross-correlation of the shock and accretion "lightcurves."The cross-correlation is itself periodic due to both signals being periodic, and peaks at a lag of 0.43  bin .Since the two signals have a period of  bin , this affirms that they are almost fully out of phase with one another, as is apparent by eye in the left panel.
The initial conditions used for the simulation, calculated from the disk model described in Section 2.2.The vertical dotted line on each panel shows the position of the secondary SMBH.For / (solid purple),   / (dashed orange) and   / (dotted magenta), where  2  ≡   /(ΣΩ 2  ) and  2  ≡   /(ΣΩ 2  ), are also plotted to illustrate the relative contributions of radiation and gas pressure, respectively, throughout the disk.
Teyssandier & Ogilvie (2017)um deficit between  = 200  g and  = 400  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 labeled 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 ∝  2 , as was done inTeyssandier & Ogilvie (2017).