The dynamics and electromagnetic signatures of accretion in unequal-mass binary black hole inspirals

ABSTRACT

We present a theoretical study of the gravitational wave (GW) driven inspirals of accreting black hole binaries with mass |$M = 10^7 M_\odot$| and mass ratios between |$10^{-3}$| and |$10^{-1}$|⁠. Our results are based on analytic estimates, and grid-based hydrodynamics simulations run for many thousands of binary orbits before the merger. We show that the GW inspiral is evident in the light curves and colour evolution of a binary-hosting quasar over years to decades before a merger. The long-term electromagnetic (EM) signature is characterized by a gradual UV brightening and X-ray dimming, followed by an X-ray disappearance hours to days before the GW burst, and finally, a years-like re-brightening as the disc relaxes and refuels the remnant black hole. These time-scales are surprisingly insensitive to the normalization of the kinematic viscosity in the disc. The spectrum of quasi-thermal disc emission shows two peaks: one in the UV and another in the X-ray, associated with the outer and circum-secondary discs, respectively; emission from the inner disc is suppressed because the secondary consumes most of the inflowing gas. We discuss implications for real-time and archival EM follow-up of GW bursts detected by LISA.

1 INTRODUCTION

The gravitational wave (GW) driven inspirals of massive black hole binaries (MBHBs) are high-priority targets for the future space-based GW detectors such as LISA, Taiji, and TianQin (Ruan et al. 2018; Gong, Luo & Wang 2021; Amaro-Seoane et al. 2023). Within the next 15 yr, LISA will observe the GW’s emitted by MBHBs in the final hours to months of their lives (e.g. Sesana 2021), including the GW burst produced as they merge, and the ring-down of the final black hole remnant (Baibhav, Berti & Cardoso 2020; Pitte et al. 2023). These observations will directly constrain the binary component masses and spins (e.g. Trias & Sintes 2008), and could reveal the presence of a circumbinary gas disc (Yunes et al. 2011; Derdzinski et al. 2019, 2021; Garg et al. 2024).

The science potential of a LISA detection is greatly increased if the host galaxy can be identified. The host morphology could reveal evidence of recent galaxy mergers (e.g. Hughes 2002; Izquierdo-Villalba et al. 2022), potentially constraining the time-scale for MBHB formation (Begelman, Blandford & Rees 1980). Also, the host redshift, combined with the GW luminosity distance, yields a measurement of the Hubble constant (Schutz 1986). However, LISA observations of MBHB inspirals could be localized to a region of |$1\!-\!10 \deg ^2$|⁠, generally too large to identify the host galaxies of individual events prior to merger (Mangiagli et al. 2020). It should be noted that post-merger, the LISA localization region can be potentially as small as |$0.01 \deg ^2$|⁠. Nevertheless, host galaxy identification in the pre-merger phase will almost certainly require an electromagnetic (EM) counterpart, and pre-merger galaxy selection can uniquely allow for multiwavelength monitoring of the source during and after the merger.

The MBHBs that will merge as LISA sources likely exist in the gas-rich nuclear regions of post-merger galaxies (Barnes & Hernquist 1996) and could thus be surrounded by an optically thick accretion flow. The flow would consist of an outer (circumbinary) disc, an inner disc around the primary (more massive) black hole, and a smaller ‘minidisc’ around the secondary (less massive) black hole. The quasi-thermal radiation released from the accretion flow could exhibit distinctive temporal or spectral characteristics, which may aid in the identification of a host galaxy, either in archival data from time-domain EM surveys or with contemporaneous tiling of the LISA error box.

Predicting such EM signatures requires a detailed understanding of how the accretion flow evolves in time throughout the binary inspiral, and many studies have now been carried out based on magnetohydrodynamics simulations in both 2D (e.g. Noble et al. 2012; Farris et al. 2015; Tang, Haiman & MacFadyen 2018; Derdzinski et al. 2019, 2021; Dittmann, Ryan & Miller 2023; Krauth et al. 2023b) and 3D (e.g. Farris, Liu & Shapiro 2010, 2011; Farris et al. 2012; Bowen et al. 2018, 2019; Pereira et al. 2019; Combi et al. 2021; Gutiérrez et al. 2022; Ruiz, Tsokaros & Shapiro 2023; Avara et al. 2024; Cocchiararo et al. 2024; Franchini et al. 2024). Each of these is based on equal-mass binary inspirals, except for Ruiz et al. (2023), which includes a model with a mass ratio as low as |$q \equiv M_2/M_1 = 0.25$|⁠, and Pereira et al. (2019) and Derdzinski et al. (2019, 2021), which include intermediate mass ratio inspirals (IMRIs; Arca Sedda, Amaro Seoane & Chen 2021), in the range of |$q = 10^{-2}$| to |$10^{-4}$|⁠. The study led by Pereira reports the time-evolving mass accretion rates, but not EM signatures, and the studies led by Derdzinski analyse the magnitude of gas-induced deviations from the vacuum post-Newtonian inspiral, but also do not characterize the EM signatures. Very recently, Cocchiararo et al. (2024) reported EM signatures for mass ratios as small as |$q=0.1$|⁠.

In this paper, we use multidimensional hydrodynamics simulations to characterize the EM signatures of unequal-mass inspirals, with q in the range of |$10^{-3}$| to |$10^{-1}$|⁠. Our interest in unequal-mass inspirals is motivated by the possibility that a large fraction of the MBHB population could have significantly unequal masses (Sesana et al. 2012; Bellovary et al. 2019), and that accretion in low mass ratio binaries seems to operate with reduced stochastic variability, especially below |$q \sim 0.01\!-\! 0.04$| (e.g. Farris et al. 2014; Shi & Krolik 2015; D’Orazio et al. 2016; Duffell et al. 2020; Dittmann & Ryan 2024), which could allow for periodic or secular variability signatures to stand out more clearly. IMRIs are also of particular interest because the lower mass component black holes will be in the range of |$10^{3-5} \ \rm {M}_\odot$|⁠, meaning that LISA observations of them would directly probe a population of intermediate-mass black holes, which has so far been elusive (e.g. Miller 2009).

Accreting IMRIs might also present distinctive spectral signatures to distinguish them from MBHBs with nearly equal component masses. This is because the less massive black hole is generally expected to consume a majority of the inflowing gas as seen in numerical simulations with thin discs (e.g. Lubow, Seibert & Artymowicz 1999; Duffell et al. 2020; Muñoz et al. 2020; Siwek et al. 2023a), an effect referred to in literature as preferential accretion (Muñoz & Lai 2016). The effect implies that emission from the secondary minidisc could dominate the system’s overall EM luminosity or that the spectral energy distribution of quasi-thermal disc emission could exhibit multiple peaks. Also, if the degree of preferential accretion changes throughout the late inspiral, that might imprint a distinctive colour evolution that could be used to photometrically select host galaxy candidates from time-domain surveys. Thus, the main goal of this study is to measure the evolution of the component accretion rates |$\dot{M}_1$| and |$\dot{M}_2$| throughout a GW-driven inspiral, and use those measurements to approximate the colour evolution of accreting IMRIs.

To the best of our knowledge, this study is the first to use multidimensional hydrodynamics simulations to model the EM emission of accreting IMRIs. However, the earliest papers on EM emission from MBHB inspirals, which were based on one-dimensional disc models, focused specifically on systems with mass ratios |$q \lesssim 0.1$|⁠. Armitage & Natarajan (2002) used the ‘standard disc’ equations (e.g. Shakura & Sunyaev 1973; Pringle 1981) to model how the inner and outer discs evolve throughout the inspiral. They used an empirical term derived from the theory of disc–satellite interactions (Goldreich & Tremaine 1980; Lin & Papaloizou 1986) to model the gravitational torque applied to the disc by the secondary. Such treatments can lead to the formation of a vacuum gap surrounding the secondary’s orbit, prohibiting mass from moving between the inner and outer discs or on to the secondary black hole. As the GW inspiral speeds up, the vacuum gap implies ‘tidal squeezing’ of the inner disc towards the primary black hole. Armitage & Natarajan (2002) speculated that this process (now sometimes referred to as a ‘snowplow’ mechanism) force-feeds the larger black hole and could drive a massive outflow and EM flare in the very late inspiral. The studies of Lodato et al. (2009) and Chang et al. (2010) are similarly based on time-dependent one-dimensional disc models and also support the possibility of a snowplow mechanism.

However, the vacuum gap around the secondary’s orbit is now understood to be an artefact of one-dimensional disc treatments. Two-dimensional calculations show that gas crosses the secondary orbit on so-called ‘horseshoe’ orbits (Lubow et al. 1999; Duffell et al. 2014; Fung, Shi & Chiang 2014). It means that the inner disc gas could be expelled to the outer disc or consumed by the secondary and that a late-inspiral enhancement of the accretion rate may not take place (Baruteau, Ramirez-Ruiz & Masset 2012). In fact, the opposite effect has been seen in simulations of equal-mass binaries (e.g. Farris et al. 2015; Dittmann et al. 2023; Krauth et al. 2023b; Franchini et al. 2024); the binary accretion rate drops throughout the late inspiral. We refer to such mass starvation as ‘viscous decoupling’, an effect first described in Liu, Wu & Cao (2003). It happens because the binary eventually contracts faster than the inward viscous spreading of the outer disc, and this leads to a diminishing rate of mass supply to the binary components. It was also pointed out in Milosavljević & Phinney (2005) that viscous decoupling of the binary should lead to a post-merger re-brightening as the disc spreads inwards following the merger and eventually begins fuelling the black hole remnant. The refuelling has been simulated for the case of post-merger equal-mass binaries in Farris et al. (2015) and Krauth et al. (2023b), and we examine it in this study for the case of intermediate mass ratios |$q < 0.1$|⁠.

A cautionary note for readers of this paper and similar works: previous works have not always made the distinction between two different time-scales in the binary–disc interaction. One time-scale, the viscous decoupling time-scale, as mentioned in the preceding paragraph, describes the time at which the binary merger rate begins to contract faster than the viscous spreading at the inner edge of the outer disc. Previous works (e.g. Sesana et al. 2012) have referred to this time-scale as the ‘disc freezing time-scale’. The second time-scale is the time at which the binary switches from contraction/expansion due to the exchange of angular momentum with the gas disc to contraction due to GW emission. In this paper, we will refer to this time-scale as the ‘transition time-scale’, but others have referred to this time-scale as a ‘decoupling time-scale’ (e.g. Armitage & Natarajan 2002; Sesana et al. 2012). These time-scales are distinct and should be calculated separately; see Sections 2.3 and 2.5 for derivations of the viscous decoupling time-scale and transition time-scale, respectively.

Our study is based on time-dependent numerical solutions of the two-dimensional (vertically averaged) Navier–Stokes equations, where gas is subject to the time-dependent gravitational potential of an unequal-mass black hole binary, undergoing a GW-driven post-Newtonian inspiral (Peters 1964). We adopt a simplified thermodynamic treatment, in which the speed of sound is set to a small fraction, |$1/\mathcal {M}$|⁠, of the nominal orbital speed (⁠|$\mathcal {M}$| is a global Mach number, selected to be in the range of |$10\!-\!20$|⁠). The bolometric EM luminosity is reasonably well predicted by measurements of the mass accretion rate, provided the flow is radiatively efficient. Owing to the use of simplified thermodynamic treatment, our simulations do not yield a self-consistent disc photosphere temperature. To model the colour evolution, we have thus developed a ‘three-disc’ toy model for the disc surface temperature, based on approximating the inner, outer, and secondary discs as stationary |$\alpha$|-discs (Shakura & Sunyaev 1973), with simulation-calibrated mass inflow rates. We believe this approach yields a reasonable first approximation of the colour evolution of accreting IMRIs. Models with self-consistent thermodynamics (e.g. Westernacher-Schneider et al. 2022, 2024) for unequal-mass binaries will be reported in a follow-up study.

Our paper is organized as follows. In Section 2, we describe the physical picture of the MBHB inspiral and provide analytical estimates for the EM emission from these systems. In Section 4.7, we describe the equations of motion, initial conditions, and the numerical set-up. Section 6.1 contains the key dynamical measurements obtained from the numerical simulations. Section 5 includes the predictions of spectral evolution and light curves based on the three-disc spectral toy model and discusses prospects for EM detections of LISA sources. Finally, in Section 6, we summarize our results, discuss some caveats, and describe plans for future work.

2 ACCRETING UNEQUAL-MASS BLACK HOLE BINARIES

2.1 Fiducial system

We adopt a fiducial black hole binary that will be detectable by LISA out to a high redshift. The total mass of the fiducial system is |$M = 10^7 \ \rm {M}_\odot$|⁠, the mass ratio is |$q = 0.01$|⁠, and the orbit is circular. The assumption of a circular orbit is justified on the basis that in the late-inspiral stage, the GW radiative back-reaction exceeds the gravitational forces of the gas, so the orbital eccentricity, driven up at wider binary separations (D’Orazio & Duffell 2021; Zrake et al. 2021; Siwek, Weinberger & Hernquist 2023b), has already been dissipated. At the leading post-Newtonian order, the rate of orbital contraction due to GW’s is given in Peters (1964) as

where |$\eta \equiv q/(1 + q)^2 = M_1 M_2/M^2$| and a is the size of the orbit. If the fiducial binary were on track to merge in the next 50 yr, then it would presently have a semimajor axis |$a_0 \simeq 6.3\,\mathrm{au}$| and an orbital period |$P_0 \simeq 1.8\,\mathrm{d}$|⁠. The time-to-coalescence is denoted as |$t_{\rm minus} \equiv \tau /4$|⁠, where |$\tau \equiv -a/\dot{a}_{\rm gw}$| is the nominal orbital contraction time-scale. The fiducial system is very mildly relativistic, with the orbital velocity roughly equivalent to |$0.17c$| and gravitational radius |$r_g \equiv GM/c^2$| about 64 times smaller than the initial semimajor axis. For later reference, the three (small) relativistic parameters |$\upsilon /c$|⁠, |$r_g/a$|⁠, and |$P/\tau$| are mutually related via equation (1),

The Eddington luminosity of a |$10^7 M_\odot$| black hole is |$L_{\rm edd} \simeq 1.15 \times 10^{45}\,\mathrm{erg\,s^{-1}}$|⁠, corresponding to a bolometric absolute magnitude of |$M_{\rm bol} \simeq -24$|⁠. For reference, such modest-luminosity active galactic nuclei (AGNs) will be detectable by LSST out to redshifts of |$z \simeq 5.5$|⁠. (Brandt et al. 2009). The LISA detection horizon for these events will be at redshifts of |$z \simeq 20$| (Sesana 2021).

2.2 Physical picture and spectral appearance

The binary is surrounded by a thin, centrifugally supported, optically thick gas flow. In binaries with a small mass ratio, the overall flow is similar to that of planet-forming discs. It comprises an outer circumbinary disc, a low-density annular gap surrounding the secondary’s orbit, and an inner disc around the primary.

The accretion flow appears as a multitemperature blackbody; the spectrum is determined by the distribution of the local disc photosphere temperatures. To get a sense for the composite spectrum of the accretion flow, including the outer disc, the inner disc, and the secondary minidisc, it is instructive to consider each of these discs as a steady-state |$\alpha$|-disc (Shakura & Sunyaev 1973), and examine the composite quasi-thermal emission spectrum. The characteristic temperatures are given by

where |$i=1$| for the primary and |$i=2$| for the secondary, and |$R_i = 6 G M_i/c^2$| are the inner cut-off radii of the component discs. For the fiducial system discussed in Section 2.1, accreting at the Eddington limit of the primary,¹ with |$\dot{M}_2/\dot{M}_1 = 10$|⁠, the characteristic temperatures of the three discs are

The secondary black hole has a more compact disc and typically accretes faster, so the secondary disc emission is bluer and brighter than the primary and outer discs.

Fig. 1 shows a toy model spectrum based on this picture and the temperature profiles in equation (3). The inner disc extends to |$r=a$|⁠, and the secondary disc extends to the Hill radius |$r_{\rm hill} \simeq a (q/3)^{1/3}$|⁠. The outer disc starts at |$r=a$| and extends to infinity. The preferential accretion parameter is set to |$\dot{M}_2/\dot{M}_1 = 10$|⁠. The spectrum is double-peaked, with the UV and X-ray bumps formed by the outer and secondary discs, respectively. The inner disc is outshined by the other components, a consequence of the secondary black hole dominating the accretion. Note that a double-peaked spectrum similar to that shown in Fig. 1 is also presented in Chang et al. (2010), but the reason is different. In that study, the secondary does not accrete, and the two spectral bumps are produced by the inner and outer discs.

2.3 Viscous decoupling time-scale

At some point during the GW inspiral, the GW contraction rate |$\tau ^{-1}$| exceeds the viscous relaxation rate (e.g. Liu et al. 2003; Milosavljević & Phinney 2005). After that time, the disc is said to have ‘viscously decoupled’ from the binary, as it can no longer spread viscously inwards as fast as the speed |$\dot{a}_{\rm gw}$| of the binary contraction. The secular evolution of the brightness and colour of the accretion flow around the binary could be significantly altered by viscous decoupling. Note that previous works (e.g. Sesana et al. 2012) have referred to this time-scale as the ‘disc freezing time-scale’. Also, the term ‘decoupling time-scale’ has been used in other works to refer to the time after which the GW power exceeds the rate of energy loss due to tidal interactions with the disc. In this paper, we refer to the transition from gas to GW driving as the ‘transition time-scale’. The transition time-scale is estimated later in Section 2.5.

The nominal viscous relaxation time-scale of the accretion disc is given by

where |$\nu (r)$| is the kinematic viscosity at radius r. In a standard |$\alpha$|-disc with constant aspect ratio |$h/r$|⁠, the viscosity coefficient can be written as |$\nu (r) = \bar{\nu }\sqrt{G M r}$|⁠. Here, |$\bar{\nu }\equiv \alpha /\mathcal {M}^2$|⁠, where |$\alpha \lesssim 0.1$| represents the Reynolds stress associated with small-scale (numerically unresolved) turbulence in the disc (Shakura & Sunyaev 1973), and |$\mathcal {M}$| is the orbital Mach number. In addition to the fiducial binary parameters introduced in Section 2.1, we also define a set of fiducial disc parameters |$\alpha = 0.1$| and |$\mathcal {M} = 10$|⁠.

The ratio of the viscous time near the binary to the GW inspiral time is

Note that we have chosen to equate the inspiral time-scale to the viscous time-scale at |$r = a$| for the viscous decoupling condition. This choice seems more appropriate when the mass ratio is small (see Figs 2 and 3) rather than at |$r = 2.5a$|⁠, which is typically used as the radius for the outer disc cavity edge in equal-mass binary studies. This ratio, as shown in equation (5), starts off small and increases as the inspiral proceeds. It exceeds unity at the viscous decoupling radius |$a_{\rm dec}$|⁠,

which we have written here, first in terms of |$r_g$|⁠, and then for a typical LISA binary with the fiducial parameters described in Section 2.1. The method we have used to determine the decoupling separation in equation (6) has been questioned in works such as Dittmann et al. (2023), where the authors show that a more accurate method to derive the decoupling separation is to equate the binary contraction speed to the viscous drift speed rather than equating time-scales. Nevertheless, our result in equation (6) is the same regardless of equating time-scales or speeds. This is due to our assumption that |$t_{\rm visc} = \tau _{\rm gw}$|⁠, while Dittmann et al. (2023) assume |$t_{\rm visc} = \tau _{\rm gw}/4$|⁠. Thus, the viscous decoupling scales we derived here are equivalent to those in equations (5) and (6) of Dittmann et al. (2023), for their case of |$\xi =1$| and |$\lambda =2$|⁠.

The residence time for the fiducial system at the viscous decoupling stage is

which corresponds to a time-to-coalescence of |$t_{\rm minus} \simeq 1.4 \ \text{d}$|⁠. The orbital period at the viscous decoupling radius is

Equation (7) shows the time-scale at which the binary is generally expected to become starved of mass supply. Note the apparent high sensitivity to the viscosity suggests that thinner |$\alpha$|-discs might viscously decouple much earlier than warmer discs; if |$\bar{\nu }$| were instead |$10^{-5}$| (⁠|$\mathcal {M} = 100$|⁠), equation (7) indicates the viscous decoupling time-scale would then be of the order of years rather than days. In Section 4.2, we show that this expectation is not supported by our simulations; we find instead that viscous decoupling extends over thousands of |$\tau _{\rm dec}$|⁠, implying that colder and warmer discs alike could show signs of the inspiral long before the nominal decoupling time.

Measurements of the preferential accretion parameter |$\dot{M}_2/\dot{M}_1$| are reported in a number of studies (e.g. Lubow et al. 1999; Young & Clarke 2015; Duffell et al. 2020; Muñoz et al. 2020; Siwek et al. 2023a), and these generally find values of |$\dot{M}_2/\dot{M}_1$| in the range of |$1\!-\!10$|⁠, when the mass ratio is larger than roughly |$10^{-3}$|⁠. All of those studies were based on simulations with a viscously relaxed disc. In Section 4.2, we present calculations revealing how the preferential accretion parameter evolves throughout the long-term viscous decoupling.

2.4 An argument against tidal squeezing

As the disc loses its viscous coupling to the binary, the feeding rates |$\dot{M}_1$| and |$\dot{M}_2$| evolve away from their steady-state values. These long-term trends will largely determine the IMRI brightness and colour evolution before coalescence. However, the feeding rates are not easy to predict without detailed simulations. The main question here is whether the GW burst is preceded by a period of AGN brightening or dimming. Armitage & Natarajan (2002) computed the evolution of a one-dimensional Keplerian disc with an embedded low-mass secondary black hole (⁠|$q=0.01$|⁠) on a GW-driven inspiral. They accounted for the gravitational torque of the secondary on the inner and outer discs, using a prescription for Lindblad torques due to Lin & Papaloizou (1986). That prescription leads to the formation of an evacuated gap around the secondary orbit, which fully insulates the mass flow between the inner and outer discs and prohibits any movement of gas on to the secondary black hole. The result, which can be seen, for example, in fig. 3 in Armitage & Natarajan (2002), is that the inner disc is squeezed towards the primary as the secondary spirals in.

Tidal squeezing would require the secondary to exert sufficient gravitational torque to the inner disc to lower the orbits of gas parcels faster than |$\dot{a}_{\rm gw}$|⁠. It implies the positive torque, |$\dot{J}_{\rm in}$|⁠, exerted by the inner disc on the secondary, must exceed

where the second equality is precise in a steady state with |$\dot{M} = 3 \pi \nu (r) \Sigma (r)$|⁠. The dependence of this torque on |$a^{-2}$| indicates that to squeeze the inner disc would require the torque to diverge leading up to merger. The gravitational torque coming from the inner disc can be written as |$\dot{J}_{\rm in} = \ell _{\rm in} J_2 \dot{M}/M$|⁠, where |$J_2 \simeq M_2 \Omega a^2$| is the secondary angular momentum, and |$\ell _{\rm in}$| is a dimensionless ‘eigenvalue’, which is typically not greater in magnitude than about 10, which is confirmed for binary mass ratio as low as |$q=0.01$| (e.g. Duffell et al. 2020). It is then straightforward to see that |$\dot{J}_{\rm plow}$| (the torque required to squeeze the inner disc) exceeds |$\dot{J}_2$| some |$1/(6 \pi \eta \bar{\nu }\ell _{\rm in})$| orbits before merger. In the fiducial binary and disc system with |$q \approx \eta \simeq 0.01$|⁠, |$\bar{\nu } = 10^{-3}$|⁠, and torque eigenvalue of order unity, this corresponds to |$\sim \!5000$| orbits or |$\sim \!25$| yr before merger. Thus, on dynamical grounds, we expect that in the final thousands of orbits, the inner disc gas is either funnelled to the outer disc or falls on to the secondary. It would require an unrealistically large gravitational torque to confine the inner disc inside the secondary’s orbit in the late inspiral. In Section 4.6, we present numerical calculations confirming these estimates.

2.5 Transition from gas to GW driving

For much of the binary lifetime, the orbital contraction is thought to be governed by coupling to the circumbinary disc (e.g. Begelman et al. 1980). As the binary contracts and the GW power increases, there is a transition from gas to GW driving (Escala et al. 2005; Cuadra et al. 2009; Lodato et al. 2009; Muñoz et al. 2020; Bortolas et al. 2021). Once the binary reaches that separation, which we denote here as |$a_{\rm trans}$|⁠, an inspiral towards merger is inevitable as GW’s remove orbital energy from the binary at a runaway rate. Note that some authors use different terminology for this separation and time-scale, for example, Sesana et al. (2012) simply calls this ‘decoupling’. We adopted the term ‘transition separation’ for the switch to GW driving, in order to reserve the term decoupling for the process of viscous decoupling.

The transition separation occurs where the rate of orbital contraction by GW radiation (equation 1) is equivalent to the rate of orbital contraction by gas i.e. |$\dot{a}_{\rm gas} = \dot{a}_{\rm gw}$|⁠. Gas-induced changes of the semimajor axis are prescribed here in terms of another dimensionless constant |$\epsilon$|⁠, as |$\dot{a}_{\rm gas} = -\epsilon a \dot{M}/M$|⁠. Equating |$\dot{a}_{\rm gw}$| and |$\dot{a}_{\rm gas}$|⁠, we find the separation at which GW’s begin to dominate the inspiral,

where |$\dot{m} \equiv \dot{M}/\dot{M}_{\rm edd}$| is the mass accretion rate normalized by the Eddington accretion rate and |$\epsilon \approx 5.0$| is consistent with results from simulations (e.g. Tiede et al. 2020; Dittmann & Ryan 2024; Tiede et al. 2024a). The time-scale for the transition from gas to GW-driven inspiral is |$\tau _{\rm trans} \equiv \tau (a_{\rm trans}) = t_{\rm sal}/\epsilon \dot{m}$|⁠, where |$t_{\rm sal} = 5 \times 10^7\,\mathrm{yr}$| is the Salpeter time. This nominal time-scale for the end of the gas-driven evolutionary phase holds when the orbital decay follows |$\dot{a}/a \sim -\dot{M}/M$| and the mass accretion rate is a fraction of the Eddington accretion rate.

Therefore, LISA inspirals will be generally deep in the GW-driven regime, and it justifies the choice of circular binary orbits in this study. Even though a detectable level of orbital eccentricity could be left over from a gas-driven evolutionary phase in LISA sources (Armitage & Natarajan 2005; Roedig et al. 2011; Zrake et al. 2021; Siwek et al. 2023b; Garg et al. 2024; Valli et al. 2024), from the point of view of the circumbinary disc interaction years before a merger, the orbit can be considered as nearly circular.

2.6 Post-merger viscous closing

Viscous decoupling implies the binary merges in a relatively low-density ‘hole’ in the disc. Here, we estimate the time-scale for the hole to be refilled following the merger, noting that such effect has been surmised to lead to a late-time brightening, or ‘merger afterglow’ (e.g. Milosavljević & Phinney 2005). The afterglow is predicted to begin a time span |$t_{\rm refill}$| following the merger, given by the viscous relaxation time-scale at the radius where the binary became viscously decoupled,

This estimate is analogous to one given in Milosavljević & Phinney (2005). Note that it suggests the remnant black hole begins refuelling on a time-scale that is quite sensitive to the disc viscosity; in an |$\alpha$|-viscosity model, the refuelling time-scale would appear to increase from days for |$\mathcal {M}=10$| to decades for |$\mathcal {M}=100$|⁠. We show in Section 4.3 that such sensitivity is not supported by the simulations, instead the remnant begins fuelling over year-like time-scales, with surprisingly little sensitivity to the normalization of the kinematic viscosity.

2.7 Importance of self-gravity

Our simulations do not include the self-gravity of the disc gas, so it is important to check when that assumption could be valid. The importance of self-gravity can be estimated from the Toomre parameter, |$Q = c_s \Omega /\pi G \Sigma$| for a Keplerian disc; discs are generally unstable to fragmentation instabilities if |$Q \lesssim 1$|⁠. When |$\alpha =0.1$| and |$\mathcal {M} = 10$| the Toomre parameter is found to be

Colder discs are more unstable to gravitational fragmentation; however equation (12) indicates that even when the Mach number is of the order of hundreds, the disc still has |$Q \gg 1$| out to |$r \sim 10a$|⁠. The effect of gravitational instabilities on long-term binary orbital evolution has been examined in, e.g. Bortolas et al. (2021) and Franchini, Sesana & Dotti (2021). Although Franchini et al. (2021) only simulated equal-mass binaries and Bortolas et al. (2021) only studied |$q \gtrsim 0.1$|⁠, self-gravity is not likely to influence the EM signatures of the late inspiral regardless of binary mass ratio.

3 NUMERICAL METHODS

Our hydrodynamic simulations are performed using the publicly available sailfish code. sailfish is a GPU-accelerated, grid-based hydrodynamics code with physics and post-processing features specifically targeting problems related to binary–disc interactions. sailfish uses a fixed-mesh, second-order Godunov solver and has been used in many published studies of binary accretion (e.g. Westernacher-Schneider et al. 2022; Krauth et al. 2023a,b; Duffell et al. 2024; Tiede & D’Orazio 2024; Westernacher-Schneider et al. 2024; DeLaurentiis et al. 2025).

3.1 Equations of motion

Our simulations are based on solutions to the vertically averaged, time-dependent mass continuity and Navier–Stokes equations,

Here, |$\Sigma$| is the vertically integrated mass density, |$\mathbf {v}$| is the gas velocity, |$P \equiv \Sigma c_s^2$| is the vertically integrated gas pressure, and |$\mathbf {I}$| is the identity tensor (pressure is isotropic). We adopt a local isothermal equation of state, |$c_s^2 = -\mathbb {V}/\mathcal {M}^2$|⁠, where |$\mathcal {M}$| is the orbital Mach number, and the gravitational potential is given by

Here, |$r_1$| and |$r_2$| are the distances to each black hole, and |$r_s$| is the gravitational softening length to ensure the potential is finite at the component positions. Use of the local isothermal equation of state means that we do not solve an energy transport equation nor adopt a phenomenological cooling prescription. Certain thermodynamic processes, such as shock-heating, are not accounted for as a result of the isothermal equation of state.

The ‘sink’ term |$\dot{\Sigma }_{\rm sink}$|⁠, appearing in equations (13) and 14, is responsible for the direct exchange of mass and momentum between the gas and the binary components,

We have checked that our choice of sink rate, |$\tau _{\rm sink}^{-1}$|⁠, does not significantly affect our measurements of the preferential accretion rate, provided the sink rate is sufficiently large. This seems to be different from Dittmann & Ryan (2024), who reported modest sensitivity to the sink rate. We suspect the difference could be in the use of |$\alpha$| versus constant-|$\nu$| viscosity. In this study, we have adopted a sink radius of |$r_{\rm sink} = 0.05 a_0$| and a sink rate of |$\tau _{\rm sink}^{-1} = 10 \Omega _0$|⁠.

Our chosen sink radius is always smaller than the Hill radius of the secondary black hole, |$r_{\rm hill} = a_0 (q/3)^{1/3} \ge 0.07a_0$|⁠, at the start of the simulation, even for the lowest mass ratio considered (⁠|$q = 10^{-3}$|⁠). However, as the binary separation decreases, the sink size of |$0.05a_0$| will become too large to resolve the Hill radius for lower mass ratios. We acknowledge this as a caveat, as our results may be sensitive to sink size, particularly for the smallest mass ratios simulated. To address this issue, we distinguish our accretion rate results based on the sink radius of the secondary black hole. When the sink radius satisfies |$r_{\rm sink} \ge 0.5 \ r_{\rm hill}$|⁠, we represent the results with dashed, less opaque lines to indicate that gas flows in the circum-secondary disc are not well resolved. For a discussion of the sensitivity of the accretion rates to the sink size and sink rate in our simulations, see Appendix  A.

The vertically integrated gravitational force density, associated with the softened gravitational potential, is |$\mathbf {F}_g = -\Sigma \, \nabla \mathbb {V}$|⁠. We assume the binary mass |$M = M_1 + M_2$| is much larger than the disc mass so that self-gravity can be safely ignored (see Section 2.7). The components of the viscous stress tensor |$\mathbf {T}_{\text{vis}}$| in equation (14) are given by |$T^{ij}_{\text{vis}} = \nu \Sigma \left(\partial ^j v^i + \partial ^i v^j-\partial _k v^k \delta ^{ij} \right)$|⁠. sailfish supports both ‘constant-|$\nu$|’ and ‘constant-|$\alpha$|’ viscosity prescriptions. For this study, we have adopted the constant-|$\alpha$| viscosity prescription with kinematic viscosity coefficient

where the modified Keplerian orbital frequency |$\tilde{\Omega }$| is

At large radii |$r \gg a$|⁠, equation (17) reduces to |$\nu \simeq \bar{\nu } \sqrt{GMr}$|⁠, where |$\bar{\nu } = \alpha /\mathcal {M}^2$|⁠.

The black holes are moved according to the Kepler two-body problem, modified to account for the GW-driven inspiral up to the leading post-Newtonian order (Peters 1964). The semimajor axis as a function of time is the solution of equation (1),

where |$\tau _0$| is the contraction time-scale at the start of the simulation. Positioning the binary components on the grid requires an expression for the orbital phase |$\phi (t)$| as a function of the simulation time t. The orbital phase is obtained by integrating in time the instantaneous orbital frequency,

The gravitational force of the gas is neglected in the binary equation of motion, which is justified because we are only considering binaries deep in the GW-driven regime.

3.2 Initial data, outer boundary condition, and diagnostics

We start our simulations of accretion on to the binary with a circular gas disc: the azimuthal velocity profile is |$\upsilon _\phi = \upsilon _{\rm kep} = \sqrt{G M/r}$|⁠, the radial velocity is inwards at the viscous drift speed, |$\upsilon _r = -3 \nu (r)/2r$|⁠, and the surface density is written as

generally as a steady-state viscous disc in the middle and on the right-hand side approximated for an |$\alpha$|-disc with viscosity as in equation (17) and with uniform orbital Mach number. Note that equation (21) is the surface density of a steady disc with a zero net angular momentum current; this fact is revisited in Section 3.3. After the start of the simulation, it takes several viscous time-scales at |$r \sim a_0$| for the flow to relax to a statistically steady state. For reference, when |$\alpha = 0.05$| and |$\mathcal {M} = 10$|⁠, the viscous time-scale at |$r=a_0$| is |$t_{\rm visc} \simeq 200 \ {\rm orbits}$|⁠. For both |$q=0.01$| and |$q=0.1$| binaries, Fig. 2 shows a snapshot of the surface density once the system has viscously relaxed to a quasi-steady state, or in other words, the surface density snapshot is taken after a thousand orbits or five viscous time-scales at |$r=a_0$|⁠.

Because our code uses Cartesian coordinates with a square-shaped domain, we impose a ‘soft’ outer boundary condition within a circular buffer zone starting at radius |$r_{\rm buf} \simeq 7.5 a_0$| (we have confirmed that our results are not changed by an increase of |$r_{\rm buf}$|⁠, see Section 3.3). Inside the buffer region, the surface density and gas velocity are driven smoothly towards a zero-torque |$\alpha$|-disc (which is also the initial condition). The domain size of our simulations is |$20 a_0$|⁠, and the number of grid zones is |$2400 \times 2400$| unless mentioned otherwise; the simulations considered in this paper have a grid spacing |$\Delta x \simeq 0.008 a_0$|⁠.

The accretion rate for each black hole is measured by separately integrating the two terms in equation (16) corresponding to the primary and secondary black hole accretion rates. Smoothing is applied to the time series of the accretion rates unless specifically mentioned. The time series of the gravitational torque is computed by integrating the gravitational torque density over the whole surface of the disc,

The gravitational torque is decomposed into the torque on the primary or the secondary and also into contributions from the inner and outer discs. We ignore the accreted spin angular momentum as we have calculated that the accretion torques are negligible.

3.3 A note about accretion rates in the steady state

Ideally, our accretion rates |$\dot{M}_1$| and |$\dot{M}_2$| would be reported relative to the theoretical inflow rate |$\dot{M}_0$| in equation (21). However, we have observed that in our simulations, the mean binary accretion rate |$\langle \dot{M} \rangle$| generally exceeds |$\dot{M}_0$| by a significant margin. We understand this ‘overshoot’ to be an artefact of our treatment of the domain outer boundary (Section 3.2).

Near the outer boundary (or the buffer region), the surface density is driven smoothly towards |$\Sigma _0(r)$| in equation (21), which corresponds to a zero-torque disc. We make this choice mainly for convenience because the torque on the binary |$\dot{J}$| cannot be controlled. Rather, it needs to be measured from the simulation. When |$\dot{J}$| is positive, a better approximation of the surface density far from the binary is (see e.g. Pringle 1981)

Thus, the target surface density |$\Sigma _0(r)$| would be generally smaller than |$\Sigma (r)$| if we had |$\langle \dot{M} \rangle = \dot{M}_0$|⁠. However, the buffer zone forces the equality |$\Sigma (r_{\rm buf}) = \Sigma _0(r_{\rm buf})$|⁠, so in general we have |$\langle \dot{M} \rangle \gtrsim \dot{M}_0$|⁠.

The amount of this ‘overshoot’ can be estimated, noting that when the binary mass ratio is small, the gravitational torque is negligible, and the net angular momentum current is then |$\dot{J} \simeq \dot{M}_2 \Omega a^2$|⁠. We check if the gravitational torque is negligible in Section 4.5 and find that for a mass ratio of |$q = 0.01$|⁠, the gravitational torque is non-negligible. However, as an approximation, we find that the ratio of the measured inflow rate |$\langle \dot{M} \rangle$| to the nominal inflow rate |$\dot{M}_0$| is then

In equation (24) we have also assumed the secondary dominates a majority of the accretion, i.e. |$\dot{M}_2 \simeq \langle \dot{M} \rangle$|⁠. In simulations where |$r_{\rm buf} = 7.5 a_0$|⁠, equation (24) predicts |$\langle \dot{M} \rangle /\dot{M}_0 \simeq 1.57$|⁠, which is very close to what we observe.

We have checked that the overshoot does not contaminate any of our derived science results, such as the shape of the accretion rate time series. Overshoots of the kind we observed here could be mitigated by guessing the value of |$\dot{J}/\dot{M}$| and then driving the solution in the buffer region towards |$\Sigma (r)$| in equation (23). Note that overshoots of the binary accretion rate have also been reported in Farris et al. (2014) and Shi & Krolik (2015) and were also discussed in Rafikov (2016). We suggest here that the finite domain size could have been responsible for those observations.

Note, on the other hand, the accretion overshoot due to the finite outer boundary introduced in this section is not the same as the accretion overshoot described in Miranda, Muñoz & Lai (2017), Dittmann & Ryan (2022), and Duffell et al. (2024). The accretion overshoot described in those works will, over many orbits, slowly approach the theoretical inflow rate |$\dot{M}_0$| as the disc settles to larger radii. The lasting overshoot we observe here is due to the finite radius of the domain outer boundary, and the gradual decline of the accretion rate towards |$\dot{M}_0$| is halted when the disc would need to be viscously relaxed at radii larger than |$r_{\rm buf}$|⁠.

4 SIMULATION RESULTS

We simulated the GW-driven inspirals of black hole binaries with mass ratios between |$q=10^{-3}$| and |$q = 10^{-1}$|⁠, surrounded by discs with orbital Mach numbers |$\mathcal {M} = 10, 15, 20$|⁠. The binary mass is |$M = 10^7 M_\odot$|⁠; however, the accretion rate time series could be rescaled by an appropriate remapping of the time coordinate to describe a binary of any mass. Results related to observer time-scales, EM light curves, and emission spectra depend on M, but they also scale with it in obvious ways. Simulations generally begin of the order of |$\tau /P \sim 10^4$| binary orbits prior to the coalescence. In Section 4.7, we show that this provides sufficient time for the disc to relax and ‘forget’ the simulation’s initial condition well before the viscous decoupling stage. Times are measured in the binary rest-frame, not in the observer time, i.e. we do not adjust the time-scales for cosmological redshift.

We have summarized the initial binary parameters used in our results in Table 1 for a binary with a total mass of |$M = 10^7 \ \rm {M}_{\odot }$|⁠, a mass ratio of |$q = 0.01$|⁠, and |$10^4$| orbits until merger. As shown in the table, these values can be scaled for different total binary masses, mass ratios, and numbers of orbits before merger using equations (1) and 2. Consequently, not all results in this work use the exact initial values given in Table 1, as our study explores a broad range of mass ratios and varying initial orbits before merger.

4.1 Evolution of disc surface density

Fig. 3 shows the surface density radial profiles of the fiducial disc, initially at |$t_{\rm minus}=1$| yr and then at four other representative times leading up to merger. At radii |$r \gtrsim a_0$|⁠, the surface density profile is approximately |$\Sigma \propto r^{-1/2}$|⁠, as for a steady state |$\alpha$|-disc in equation (21). The surface density evolution is similar in certain ways to the one-dimensional disc model from Chang et al. (2010). An annular gap is formed around the secondary orbit, and the outer disc follows the secondary inwards up until the viscous decoupling time. However, the gap in our simulation is not fully evacuated; gas flows across the gap and accretes to the secondary. Also, the surface density of the inner disc in our simulation is seen to be modestly decreasing over time, whereas in Chang et al. (2010), the presence of the vacuum gap around the secondary orbit requires the surface density of the inner disc to increase. We show in the next subsection that the decrease of the inner disc mass is driven mainly by accretion to the secondary. In other words, negligible mass passes the secondary’s orbit, and instead, most gas is captured by the secondary, gradually starving the primary.

The surface density radial profiles for varying Mach number discs are shown in Fig. 4 at 4000 orbits before merger (dashed) and at merger (solid). As the Mach number increases, the viscosity decreases, and while the inflow rate |$\dot{M}_0 = 3\pi \nu \Sigma$| is held constant across all runs, the surface density magnitude must then increase. This explains the shift in surface density for higher Mach number discs.

4.2 Long-term viscous decoupling

Fig. 5 shows the binary accretion rate |$\dot{M}$| (top), the accretion rates to the binary components |$\dot{M}_1$| and |$\dot{M}_2$| (middle), and the preferential accretion rate |$\dot{M}_2/\dot{M}_1$| (bottom). The curves show the time series data for binaries with various mass ratios. Our main results here can be summarized: (1) the secondary generally starves more slowly than the primary, and (2) the binary accretion rate |$\dot{M}$| shows a long-term downward trend, being reduced from |$\langle \dot{M} \rangle$| by 10 per cent as early as |$\sim \! 500 \tau _{\rm dec}$| (decade-like time span) before the merger.

The second result is surprising, as one could reasonably guess that |$\dot{M}$| would be significantly reduced from |$\dot{M}_0$| only within a few |$\tau _{\rm dec}$| of the merger, perhaps following a functional form like |$\dot{M} = \dot{M}_0 e^{-\tau _{\rm dec}/t_{\rm minus}}$|⁠. Instead, our results suggest the quantity |$\delta \equiv 1-\dot{M}/\dot{M}_0$| might grow as a power law in |$t_{\rm minus}$|⁠, at least while |$\delta$| is small. We summarize next a model of the viscous decoupling process, which accounts for the time evolution of the binary-disc torque and predicts such scaling of |$\delta$|⁠. A full derivation is given in Zrake, Clyburn & Feyan (2025).

Our model characterizes the disc as having two distinct zones: a close zone that is tightly coupled to the binary and a distant zone that is relatively frozen over the remaining time before the merger. The two zones meet at a viscous radius, |$r_\nu$|⁠, satisfying |$t_{\rm visc}(r_\nu) = -a/\dot{a}_{\rm gw}$|⁠.² In the close zone |$r < r_\nu$|⁠, the disc evolves through a sequence of steady states (equation 23), each conducting angular momentum inwards at rate |$\dot{J} = \dot{M} \Omega a^2$|⁠. At each step n of the sequence, a boundary condition is imposed, |$\Sigma _n(r_n) = \Sigma _{n-1}(r_n)$|⁠, where |$r_n$| is the viscous radius at step n. Such recursion leads to a first-order ordinary differential equation for |$\delta$|⁠, and the solutions go asymptotically as |$\delta \propto t_{\rm minus}^{-5/24}$| (Zrake et al. 2025). Such weak dependence of |$\delta$| on |$t_{\rm minus}$| could account for the long-term viscous decoupling seen in our simulations. Fig. 6 shows a plot of |$\delta (t)$| for the fiducial binary compared to the asymptotic prediction |$\delta \propto t_{\rm minus}^{-5/24}$|⁠. With many decoupling time-scales before the merger, the simulations agree with the predicted power-law dependence.

In Section 2.3, we estimated that the viscous decoupling time-scale would be very sensitive to the Mach number of the disc, |$\tau _{\rm dec} \propto \mathcal {M}^{16/5}$|⁠. However, given our argument above, where |$\delta$| scales as a power law in |$t_{\rm minus}$|⁠, it seems reasonable to expect that the late-inspiral accretion rate curves might have only modest sensitivity to the Mach number. To test this, we ran a suite of simulations with a range of Mach numbers, |$\mathcal {M} = 10, 15, 20$|⁠. The results appear in Fig. 7, and do confirm this expectation. The run with Mach 20 has a nominal viscous decoupling time roughly 10 times larger than the run with Mach 10, and yet the curves differ mainly by a vertical offset. The long-term decoupling effect seen here also seems to be present in simulations of equal-mass inspirals in, e.g. fig. 6 in Farris et al. (2014), or fig. 3 of Dittmann et al. (2023).

4.3 Remnant fuelling

Equally important to the pre-merger accretion evolution is the post-merger remnant refuelling, where the disc viscously refills the cavity and accretes on to the merger remnant. In Section 2.6, we derived the refilling time-scale |$t_{\rm refill}$| (equation 11) and found that for the fiducial binary, the disc should refill the cavity, and begin accreting to the black hole merger remnant roughly 4 d after the merger. Note, however, that our estimate predicts |$t_{\rm refill} \propto \mathcal {M}^{16/5}$|⁠, analogous to the viscous decoupling time-scale |$\tau _{\rm dec}$|⁠. It would imply that the onset of refuelling could be delayed from days (Mach 10) to decades if the disc had a Mach number of 100.

In Fig. 8, we plot the total binary accretion rate post-merger out to 5 yr after coalescence for a binary with mass ratio |$q=0.01$| and for simulations with different Mach numbers. The result indicates that the onset of refuelling happens over a time span of months to years, with relatively little sensitivity to our derived estimate for |$t_{\rm refill}$|⁠. We thus suggest that the refilling time-scale post-merger is not well represented by a time constant such as that given in equation (11), as the recovery of the accretion rate to |$\dot{M}_0$| does not show any exponential characteristics. The exact functional form of the refuelling curve can be predicted from the model, which we summarized in Section 4.2, and comparisons are presented in Zrake et al. (2025).

4.4 Variability of mass accretion rates

Low mass ratio binaries seem to accrete with lower stochastic variability than binaries of near-equal mass (e.g. Farris et al. 2014; Shi & Krolik 2015; D’Orazio et al. 2016; Duffell et al. 2020; Dittmann & Ryan 2024). This suggests that accreting IMRIs could produce ‘cleaner’ light curves, possibly making it easier to identify them based on EM periodicity, or colour evolution, in the years before a merger.³ In this section, we show that the reduced stochasticity seen by Duffell et al. (2020) and Dittmann & Ryan (2024) applies throughout a GW-driven inspiral.

The time series of the secondary accretion rate |$\dot{M}_2$| is shown in Fig. 9 for a |$q=0.01$| binary, and in Fig. 10 for a |$q=0.1$| binary. The amplitude of the stochastic variability in |$\dot{M}_2$| is less than roughly 10 per cent for the case of |$q=0.01$| and greater than roughly 50 per cent for the case of |$q=0.1$|⁠, which confirms there is a significant increase in the variability amplitude somewhere between |$q=0.01$| and |$q=0.1$|⁠, and that this change persists also in the late-inspiral phase, at least before the viscous decoupling.

The periodic variability of unequal-mass systems can be seen in the insets of Figs 9 and 10. The secondary accretion rate modulates at the binary orbital period, and we believe this could indicate a slight eccentricity of the outer disc (e.g. Kley & Dirksen 2006). For lower mass ratios, a second peak is more predominant in the secondary’s accretion rate. This half-orbital periodicity has been seen in other works (e.g. Farris et al. 2014; Cocchiararo et al. 2024; Dittmann & Ryan 2024), and is likely the signature of an eccentric outer disc.⁴ We have measured the eccentricity of the outer disc in our runs by summing the local gas eccentricity vectors and found the eccentricity to be approximately |$2~{{\ \rm per\ cent}}$| for the |$q=0.01$| binary and around |$14~{{\ \rm per\ cent}}$| for the |$q=0.1$| binary.

4.5 Evolution of binary torques

Fig. 11 shows the time series of the gravitational torque on the binary for three mass ratios (⁠|$q = 0.002, 0.01, 0.02$|⁠). The top panel shows the non-dimensional torques |$\ell _{\rm tot} = (\dot{J}_{\rm tot}/J)/(\dot{M}/M)$| and |$\ell _{\rm in} = (\dot{J}_{\rm in}/J)/(\dot{M}/M)$|⁠, and the bottom panel shows the inner disc gravitational torque |$\dot{J}_{\rm in}$|⁠, normalized by |$\dot{M}_0 \sqrt{G M a_0}$|⁠. Fig. 12 shows the gravitational torque on the binary for simulations of varying Mach numbers. The top panel shows the total gravitational torque exerted on the binary, and the middle and bottom panels show, respectively, the inner and outer disc contributions.

In the top panel of Fig. 12, we find that the net gravitational torque on the binary is negative for all of our runs and that it becomes increasingly negative with increasing Mach number. This trend is consistent with other studies, including Tiede et al. (2020) and Penzlin et al. (2022), which also find that the torque on the binary becomes increasingly negative with increasing Mach number, albeit these studies focused on equal-mass or high mass ratio binaries. Dittmann & Ryan (2024), on the other hand, simulated unequal-mass binaries down to |$q=0.01$| and measured the gravitational torque of a |$q=0.01$| binary with |$\mathcal {M}=10$| to be marginally positive with |$\dot{J}_{\rm grav} \simeq 0.1$| and no clear trend for varying Mach number. Our results in Fig. 12, for the same binary system, show a marginally negative total gravitational torque of |$\dot{J}_{\rm grav} \simeq 0.1$| with a clear trend of monotonically decreasing torque with increasing Mach number. We believe the (small) difference might be attributable either to the difference in the viscosity prescription (constant-|$\alpha$| versus constant-|$\nu$|⁠) or differences in the sink rate, as Dittmann & Ryan (2024) found some sensitivity of their torques to the value selected for the sink rate.

The outer disc is responsible for the increasingly negative torque at a higher Mach number; the inner disc actually exerts a slightly positive torque on the binary, which is rather insensitive to the Mach number. The increasing magnitude of torque from the outer disc can be attributed to the accumulation of mass outside the secondary orbit (Fig. 4).

These results could be used to predict how the gravitational waveform of IMRIs surrounded by gas should differ from those evolving in vacuum (e.g. Derdzinski et al. 2019, 2021; Tiede et al. 2024b). Such predictions will be reported in a future study focusing specifically on GW signatures.

4.6 Tidal squeezing of the inner disc

As discussed in Section 2.4, the theory of disc–satellite interactions implies that the orbiting secondary generally exerts a negative gravitational torque to the inner disc. When prescriptions for that torque are utilized in a one-dimensional disc model, they can lead to the formation of a zero-density, insulating gap in the disc surrounding the secondary. If the secondary is on a GW-driven inspiral, the insulating gap (which prohibits the flow of gas across the secondary’s orbit) then implies tidal squeezing of the inner disc (so-called ‘snowplow’ mechanism), and a potentially dramatic enhancement of the accretion rate |$\dot{M}_1$| to the primary in the very late inspiral. In Section 2.4, we referenced other studies that have observed that the insulating gap does not form in two-dimensional hydrodynamics calculations, and we also pointed out that the magnitude of the gravitational torque, denoted as |$\dot{J}_{\rm plow}$|⁠, needed to reduce the inner disc radius at the rate |$\dot{a}_{\rm gw}$|⁠, diverges in the late inspiral. We report here the time series of gravitational torque |$\dot{J}_{\rm in}$|⁠, applied to the binary by the inner disc (as well as |$\dot{J}_{\rm tot}$|⁠, the torque from the whole disc), and confirm our estimate from Section 2.4 that in the final thousands of orbits, |$\dot{J}_{\rm in}$| is exceeded by |$\dot{J}_{\rm plow}$|⁠, meaning that tidal squeezing becomes ineffective well before the merger.

Fig. 11 shows time series data from the gravitational torques on the binary for three mass ratios (⁠|$q = 0.002, 0.01, 0.02$|⁠). The bottom panel shows |$\dot{J}_{\rm in}$| (solid) and |$\dot{J}_{\rm plow}$| (dashed) as derived in equation (9). The filled circles are where |$\dot{J}_{\rm in} = \dot{J}_{\rm plow}$|⁠. Notice that |$\dot{J}_{\rm plow}$| diverges close to merger; the secondary black hole cannot provide enough torque to squeeze the inner disc in the final |$\sim \!4000$| orbits, consistent with the estimate given in Section 2.4. The top panel of Fig. 11 confirms an assumption we made in that estimate: that the non-dimensional torques do not change radically throughout the inspiral.

4.7 Numerical convergence

We ran many simulations to check the sensitivity of our results to numerical parameters, including the grid resolution and the start time of the simulation relative to the binary merger; see Appendix  A for additional discussions on the sensitivity of our results to the sink size and sink rate. Fig. 13 presents the accretion rates of the fiducial binary for four simulations with different start times. The startup transients settle after roughly a hundred to a thousand orbits, and henceforth the accretion rates show very little sensitivity to the simulation start time. Also, in Fig. 13, we have plotted the total accretion rate for a fiducial binary that is not inspiralling. This also helps assure us that the long-term viscous decoupling discussed in Section 4.2 is uniquely a consequence of the inspiralling binary.

We also checked that our accretion rate and gravitational torque time series data are well converged with respect to the grid resolution. In Fig. 14, we have plotted the total accretion rate, preferential accretion rate, and gravitational torque for simulations performed with different resolutions and for two different mass ratios. Resolution is quantified by the grid spacing |$\Delta x/a_0$| (better resolution means smaller |$\Delta x$|⁠). Some disagreement is seen in both |$\dot{M}_1$| and |$\dot{M}_2$| and the gravitational torques when the grid spacing is greater than |$\Delta x \simeq 0.016a_0$|⁠, especially in the very late-inspiral stage when the distance between the binary components becomes poorly resolved. However, at early times, the shapes of the accretion rate curves, as well as the torques, are reasonably consistent with one another when the grid spacing is smaller than |$\Delta x \simeq 0.016 \ a_0$|⁠. Additionally, we see convergence for both the fiducial mass ratio |$q=0.01$| and for the lowest mass ratio tested |$q=0.001$| in both the accretion rates and torques in Fig. 14.

5 ELECTROMAGNETIC SIGNATURES

The inspiral of unequal-mass MBHBs might produce tell-tale EM signatures. The detections of such EM signatures by ground or space-based observatories, especially alongside a GW detection by LISA, would provide an abundance of knowledge on binary evolution, accretion, and environments of MBHBs within galaxies. To maximize the likelihood of an EM detection, we must understand how the hydrodynamical interactions between the binary and disc affect the thermal radiation released from the accretion flow. In this section, we use the spectral toy model from Section 2.2, combined with the measured mass accretion rates, to provide a reasonable first approximation of the brightness and colour evolution of an accretion flow surrounding an inspiralling black hole binary of significantly unequal mass. Note that a full model of the variable AGN emission around the time of a merger will require a more detailed treatment of the radiative transfer, both within the disc and also throughout the larger environment. Nevertheless, we do believe the analysis presented here qualifies as a first approximation of the variable AGN disc emission before and after an MBHB merger.

5.1 Spectral evolution and light curves

In Fig. 15 we plot the composite spectral energy distribution of the accretion flow in the fiducial binary simulation, at different times throughout the inspiral. We used the three-disc toy model from Section 2.2, in which the spectrum is computed by approximating each disc (inner, secondary, and outer) to emit as a blackbody. The inner cut-off radii for the inner and secondary disc are |$R_i = 6GM_i/c^2$|⁠, the innermost stable circular orbits of the respective black hole components. The inner radius of the outer disc and the outer radius of the inner disc are set to the instantaneous separation of the binary, while the secondary disc is extended to the instantaneous Hill radius. The result is not sensitive to the outer radius of the outer disc.

In Fig. 16, we plot the soft UV (top), hard UV (second), soft X-ray (third), and hard X-ray (bottom) light curves derived from our three-disc toy model of the disc temperatures. We define the soft UV energy range as |$1.5\!-\!4$| eV, the hard UV energy range as |$4\!-\!100$| eV, the soft X-ray energy range as |$0.2\!-\!2$| keV, and the hard X-ray energy range as |$2\!-\!12$| keV. Note that once the binary reaches the decoupling separation, we have manually imposed in our model to evolve the outer disc cavity wall at the viscous drift rate rather than at the binary merger rate to highlight the difference in the EM signatures pre- and post-viscous decoupling. This is the reason for the abrupt change in the light curves seen around a couple of days prior to the merger, at the vertical dashed line corresponding to the decoupling time-scale in Fig. 16.

Years before merger, the majority of the overall accretion power is released as X-rays from the small secondary disc. However, as the inspiral proceeds, the outer disc reaches deeper into the gravitational potential of the primary, glowing brighter in the UV, while the secondary is gradually starved and gets dimmer in the X-rays. The outer disc follows the inspiralling secondary black hole, and as it does so, the inner edge of the outer disc increases in speed and heats up the gas, leading to a brightening effect. The UV brightening peaks at the time when the disc viscously decouples from the binary as seen in Fig. 15 at |$t_{\rm minus} \simeq 1$| d. For the fiducial disc–binary model, the X-ray emission is surpassed by the UV emission of the order of 1 month before the merger. At the same time as the peak of the UV emission, the X-ray emission from the secondary disc decays because of the fast reduction in the rate of mass supplied to the secondary very late in the inspiral (Krauth et al. 2023b). This is also evident in the measured accretion rates in Fig. 5.

We do not show the post-merger EM signatures associated with the remnant fuelling curves presented in Section 4.3, because the disc is far from being viscously relaxed while the accretion rate recovers to |$\dot{M}_0$|⁠. Analytic estimates for this phase are presented in Milosavljević & Phinney (2005), Schnittman & Krolik (2008), and Shapiro (2010). We also do not show the orbital period variations in the light curves or spectrum, as these cannot generally be captured in the three-disc model. Spectral evolution on the time-scale of individual binary orbits requires a self-consistent thermodynamic treatment and will be presented in future work (see Section 6.2).

5.2 Prospects for EM detections

Here, we consider how our results might influence proposed strategies to discover the EM counterparts of LISA GW detections. As discussed in the introduction, EM counterparts from LISA sources will likely be needed to confidently identify the host galaxies of MBHB inspirals.

Our results indicate that in the years before a merger, accreting unequal-mass black hole binaries get brighter in the UV and dimmer in the X-ray. They also exhibit an abrupt disappearance in X-rays in the final hours to days before the merger. We further predict that variations at the orbital period should be discernible at the tens-of-per cent level, resulting from orbital-period modulations of the accretion on to the secondary. Observatories carrying out time domain surveys, or space observatories capable of slewing in response to a GW trigger, could thus be instrumental in EM counterpart identification. We discuss below the relevance of these three predictions (periodicity, UV brightening, X-ray dimming) to current and upcoming observatories.

With regard to periodicity in the EM signatures, the Nancy Grace Roman Space Telescope, a wide-field optical to near-infrared space telescope, could potentially detect these EM signatures in an accreting MBHB, corresponding to the quasi-periodic variability of the mass accretion rates shown in the insets of Figs 9 and 10. Such periodicity detections by Roman could become extremely valuable if the system is |$10\!-\!20$| yr away from merging because the merger might then take place during the LISA period of operation (Haiman et al. 2023). As discussed in Section 2.1, a |$10^7$| M|$_\odot$| binary that is 50 yr away from merging has a rest-frame orbital period of |$1.8\,\mathrm{d}$| and in Fig. 9, the period is seen to decrease to |$0.9\,\mathrm{d}$|⁠, around 10 yr before merging. At |$z \sim 1\!-\!2$|⁠, this corresponds to multiday periodicity, around a factor of 2 shorter than the anticipated |$\sim \!5$| d cadence planned for the Roman high latitude time domain survey. Our calculations in Section 4.4 only predict the periodic variations of the mass accretion rates, so follow-up calculations will be needed to predict the level of periodicity expected where Roman is sensitive. This is an intended area of future work; see Section 6.2.

The primary EM signature predicted in this work is a gradual increase in UV brightness from the circumbinary disc, aligning with the contraction rate of the binary driven by gravitational radiation emission. The Legacy Survey of Space and Time (LSST) of the Vera Rubin Observatory will observe nearly half the sky, |$18\,000\,\mathrm{deg^2}$|⁠, and will detect up to |$10^6$| galaxies per square degree, observing quasars with a cadence of once per day (Xin & Haiman 2021, 2024). Vera Rubin’s LSST camera will be sensitive in the soft UV energy range up to 4 eV. In the top panel of Fig. 16, we have plotted the light curve for the approximate LSST sensitivity. The steady UV brightening we are predicting to accompany a LISA inspiral could be evident in LSST archival data over the 10 yr prior to a GW detection by LISA; however, the UV-brightening is more pronounced and thus easier to detect in the hard UV energy ranges (second panel of Fig. 16). Therefore, telescopes in the hard-UV energy range, such as the recently selected NASA Explorer mission, UVEX, could be better suited to detect the predicted UV brightening (Kulkarni et al. 2021). Alternatively, since the characteristic temperature of the disc’s inner edge decreases with the black hole mass, it is possible that LSST will be sensitive to the variable emission from higher mass AGNs.

The X-ray dimming after viscous decoupling predicted in equal-mass binaries by Krauth et al. (2023b) and confirmed for unequal-mass binaries here could be detectable by the Athena X-ray telescope and/or eROSITA. Athena could potentially operate during the lifetime of LISA to have simultaneous multimessenger detections, while eROSITA is a current X-ray detector that could provide a reference sky template to aid in host galaxy characterization (Merloni et al. 2012). In the bottom two panels of Fig. 16, we have plotted the light curves for both the soft and hard X-ray energy ranges. The soft X-rays will be brighter initially, making them easier to observe; however, both energy ranges exhibit a significant drop in brightness due to viscous decoupling.

6 SUMMARY

6.1 Main results

Using analytic estimates and high-resolution grid-based hydrodynamics simulations, we have explored the gas dynamics and EM signatures of accreting binary black holes with significantly unequal masses, |$q = 10^{-3}\!-\!10^{-1}$| and |$M \sim 10^7 M_\odot$|⁠, over decade-like time-scales before a merger. Our aim is to help enable EM counterpart detection and host galaxy identification of binary black hole inspirals observed by LISA. Our main results are summarized here:

• The viscous decoupling time is generally of the order of |$\tau _{\rm dec}\sim$| days for |$\sim 10^7 M_\odot$| binaries. One might reasonably expect the binary accretion rate |$\dot{M}$| to remain steady until |$\sim \tau _{\rm dec}$| before merger; however, our simulations show significant reductions in |$\dot{M}$| over decade-like time-scales, |$10^3 \tau _{\rm dec}$|⁠. We attribute this effect to the decrease of angular momentum current through the disc, accompanying the binary contraction. A model of this process is the subject of a forthcoming work.

• Low mass ratio binaries show accretion rate variability at the orbital period, at the level of around 10 per cent. Stochastic variability is reduced throughout the inspiral when the mass ratio is |$q=0.01$| compared to |$q=0.1$|⁠, consistent with the trend reported in Dittmann & Ryan (2024).

  The secondary generally consumes most of the mass coming to the binary (preferential accretion effect). Both components are gradually starved of mass supply as the inspiral progresses, but the secondary starves more slowly than the primary.

• EM flares produced by tidal squeezing of the inner disc (snowplow effect) are unlikely. This is based on dynamical considerations and confirmed in the simulations; the gravitational torque becomes too weak years before merger to reduce the orbits of inner disc gas parcels faster than |$\dot{a}_{\rm gw}$|⁠. Instead, the inner disc gas accretes to the secondary during the inspiral, leading to the observed slower starvation of the secondary.

• We developed a toy model for the quasi-thermal emission spectrum of the accretion flow surrounding a low mass ratio binary black hole based on the assumption that the outer, inner, and secondary discs are viscously relaxed. The simple model predicts an emission spectrum with peaks in the UV and X-ray, reflecting the characteristic temperatures of the outer and circum-secondary discs, respectively. Those discs outshine the inner (primary) disc because of the preferential accretion effect.

• Together with the mass accretion rates measured in the simulations, the toy model predicts a decade-like gradual brightening in the UV, and dimming in the X-ray, followed by an abrupt disappearance of the X-rays around time |$\tau _{\rm dec} \sim \mathrm{days}$| before the merger (analogous to Krauth et al. 2023b), and finally year-like gradual X-ray re-brightening associated with fuelling of the remnant black hole. The long-term pre-merger UV brightening is the result of the increasing peak temperature of the outer disc, as it extends deeper into the potential well of the primary throughout the inspiral.

6.2 Caveats and future work

Our simulations are based on two-dimensional (vertically averaged) solutions of the Navier–Stokes equations in the local isothermal approximation; the effects of shock heating, radiation, magnetic fields, general relativity, and out-of-plane gas flows, including winds and jets, are all neglected. Other groups are actively exploring these, especially in the context of near-equal mass ratio binaries in the final tens of orbits before a merger (see references given in the introduction). The simplified simulations presented here were necessary to study the long-term dynamics and EM signatures and may be accurate enough to help identify host galaxies of LISA events.

In future work, we will present simulations of unequal-mass inspirals with a more realistic thermodynamics treatment to properly account for shock heating, as done with the sailfish code in Westernacher-Schneider et al. (2022), Krauth et al. (2023a,b), DeLaurentiis et al. (2025), and Westernacher-Schneider et al. (2024). In the context of near-equal-mass binaries, shock-heating does significantly alter the profile of disc surface temperature (e.g. Farris et al. 2015). However, when the mass ratio is small, we expect less of the gas orbital energy to be dissipated in strong shocks, so the result might not be very different from the spectrum toy model we presented here.

We do not expect general relativistic effects to significantly alter the long-term brightness and colour changes of a binary-host quasar in the years before a merger because the long-term changes are driven by large-scale viscous relaxation of the disc far from the gravitational radii of either black hole. Accounting for magnetic fields and radiation could change some of our predictions, especially in a regime where the mass flow to the binary is comparable to the Eddington rate of the primary. In such cases, preferential accretion implies the secondary tries to accrete much faster than its own Eddington rate and that radiation-driven outflows could significantly change the global mass currents in the system.

ACKNOWLEDGEMENTS

MC acknowledges support from the NASA Future Investigators Program (FINESST) through Award No. 80-NSSC-23K1443. JZ acknowledges support from the LISA Preparatory Science Program (LPS) through NASA Award No. 80-NSSC-24K0440. JZ also acknowledges valuable discussions with Alex Dittmann, Sterl Phinney, Julian Krolik, Elena Rossi, and Zoltan Haiman, some of which took place at the 2022 binary accretion workshop hosted by the Kavli Institute for Theoretical Physics (KITP); this research was supported in part by grant NSF PHY-1748958 to KITP. The authors would also like to thank the anonymous referee for their detailed reading of our manuscript and constructive input. All simulations were performed on Clemson University’s Palmetto cluster.

DATA AVAILABILITY

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

Footnotes

REFERENCES

APPENDIX A: SENSITIVITY TO NUMERICAL PARAMETERS

As mentioned in Section 3.1, the accretion rates of low mass ratio binaries are sensitive to the sink size of the black holes and should be interpreted with caution. In Figs A1 and A2, we plot the primary, secondary, and preferential accretion rates for different sink rates and sink radii, respectively. All simulations discussed in this section use a binary mass ratio of |$q=0.01$| and an orbital Mach number of |$\mathcal {M}=10$|⁠.

In Fig. A1, we present simulations with sink rates of |$\tau _{\rm sink}^{-1} = 5, 10, 20 \ \Omega _0$|⁠, each using a sink size of |$r_{\rm sink} = 0.05 \ a_0$|⁠. For the results presented in Section 6.1, we use a sink rate of |$\tau _{\rm sink}^{-1} = 10 \ \Omega _0$|⁠, which, according to Fig. A1, is high enough that the accretion rates are relatively insensitive to the sink rate. However, the accretion rate to the primary BH becomes spuriously large, and the accretion rate to the secondary BH becomes spuriously small when the sink rate is too low, i.e. |$\tau _{\rm sink}^{-1} < 10 \ \Omega _0$|⁠. Note that the preferential accretion rate increases with increasing sink rate. Although the components’ accretion rates are nearly identical at high sink rates, the preferential accretion depends on the ratio |$\dot{M}_2/\dot{M}_1$|⁠. Consequently, even small variations in |$\dot{M}_1$| between runs with different sink rates can lead to appreciable changes in the preferential accretion.

In Fig. A2, we present simulations with sink radii of |$r_{\rm sink} = 0.05, 0.025, 0.0125 \ a_0$|⁠, each with a sink rate of |$\tau _{\rm sink}^{-1} = 10 \ \Omega _0$|⁠. As mentioned in Section 3.1, all simulations in this manuscript begin with |$r_{\rm sink} < r_{\rm hill}$|⁠, but during inspiral, the Hill radius decreases, making the flow of gas around the secondary increasingly difficult to capture. For low mass ratios, the accretion rate to the primary BH becomes spuriously small, and the accretion rate to the secondary BH becomes spuriously large when the sink size is too large. Thus, the preferential accretion in the bottom panel of Fig. A2 is enhanced for smaller BH sinks. Nevertheless, the sensitivity of the accretion rates to the sink size does not affect our key results of a long-term viscous decoupling and a gradual UV brightening during inspiral.