# The late inspiral of supermassive black hole binaries with circumbinary gas discs in the LISA band

The late inspiral of supermassive black hole binaries with circumbinary gas discs in the LISA band Abstract We present the results of 2D, moving-mesh, viscous hydrodynamical simulations of an accretion disc around a merging supermassive black hole binary (SMBHB). The simulation is pseudo-Newtonian, with the BHs modelled as point masses with a Paczynski–Wiita potential, and includes viscous heating, shock heating, and radiative cooling. We follow the gravitational inspiral of an equal-mass binary with a component mass Mbh = 106 M⊙ from an initial separation of 60rg (where rg ≡ GMbh/c2 is the gravitational radius) to the merger. We find that a central, low-density cavity forms around the binary, as in previous work, but that the BHs capture gas from the circumbinary disc and accrete efficiently via their own minidiscs, well after their inspiral outpaces the viscous evolution of the disc. The system remains luminous, displaying strong periodicity at twice the binary orbital frequency throughout the entire inspiral process, all the way to the merger. In the soft X-ray band, the thermal emission is dominated by the inner edge of the circumbinary disc with especially clear periodicity in the early inspiral. By comparison, harder X-ray emission is dominated by the minidiscs, and the light curve is initially more noisy but develops a clear periodicity in the late inspiral stage. This variability pattern should help identify the electromagnetic counterparts of SMBHBs detected by the space-based gravitational-wave detector LISA. accretion, accretion discs, black hole physics, hydrodynamics 1 INTRODUCTION Supermassive black holes (SMBH) are currently believed to reside in most, if not all galactic nuclei, and SMBH binaries (SMBHBs) are expected to be produced frequently in mergers of galaxies (e.g. Kormendy & Ho 2013). The resulting compact SMBHBs are likely embedded in a gaseous environment. Electromagnetic (EM) signatures of such compact SMBHBs can arise from their interaction with gas. The space-based gravitational-wave (GW) detector LISA will be sensitive to SMBHBs with masses in the range 104–107 M⊙ (Amaro-Seoane et al. 2017). Identifying a GW source in the EM bands has considerable importance for astrophysics and cosmology (e.g. Phinney 2009). It has been proposed that periodic emission in the EM bands may track the orbital motion of SMBHBs throughout the late inspiral, and that this would then allow a unique identification of a LISA source (Kocsis, Haiman & Menou 2008). Comparing the EM and GW chirp signals would then also help probe the difference in the propagation speed of photons and gravitons (Kocsis et al. 2008; Haiman 2017), complementing measurements of the graviton speed from phasing of the GWs alone (see e.g. the review by Will 2006 and references therein). A similar constraint has recently been derived from the time delay between the GWs and the accompanying gamma-ray burst arriving from the neutron-star merger GW170817 (Abbott et al. 2017). An SMBHB offers an improved and more robust measurement of this time delay; this is especially the case if Doppler modulations of the EM chirp are detected from the orbital motion of the BHs themselves, fixing the relative phase of the EM and GW signals (Haiman 2017). In order to predict the EM signatures of SMBHB mergers in the LISA band, the first step is to resolve the gas dynamics during this late inspiral stage. The morphology of the gas, even on scales as large as that of the circumbinary disc, is sensitive to the behaviour of gas streams in the vicinity of the individual BHs (Tang, MacFadyen & Haiman 2017; Bowen et al. 2017a). One necessary requirement for reliably simulating the LISA stage is to resolve the accretion flows near the individual BHs, on scales down to that of the innermost stable circular orbit (ISCO; located at 6GMbh/c2 ≡ 6rg for a non-spinning BH). Furthermore, the gas dynamics cannot be in a steady state, since the binary orbit is evolving rapidly. Depending on the component masses, for a circular, comparable-mass binary, the LISA band corresponds to the inspiral from an initial separation of (50–200)GM/c2 all the way to the merger, and covers the last ≈100 to a few thousand orbits of the binary. The orbital evolution of the binary is likely overwhelmingly dominated by its gravitational radiation at this late stage (Haiman, Kocsis & Menou 2009; Dotti, Sesana & Decarli 2012; Tang et al. 2017). Because of the dual requirements of high spatial resolution and large number of orbits, simulating the LISA stage in its entirety has remained beyond the capability of fully three-dimensional (GR)MHD simulations. They have generally focused instead on the last ≈10 orbits preceding the merger (Bode et al. 2010; Farris et al. 2012; Giacomazzo et al. 2012; Noble et al. 2012; Shi et al. 2012; Bowen et al. 2017b). However, following the binary throughout the entire LISA band is feasible in 2D hydrodynamics. Previous high-resolution Newtonian simulations have indeed followed binaries for thousands of orbits. However, these works have focused on large separations, well before the LISA stage (Cuadra et al. 2009; Krolik 2010; Roedig et al. 2012; Farris et al. 2015a). As such, they did not directly resolve accretion inside the ISCO (typically either ignoring accretion, or implementing accretion via a sink prescription). In our earlier work (Farris et al. 2015b), we performed a 2D simulation of an inspiraling and merging BH binary, embedded in a viscous circumbinary disc. However, this earlier study was purely Newtonian, and employed an artificial sink that depleted the gas according to the fluid viscous time-scale. Moreover, the focus was on the overall accretion rate, and we did not present an analysis of the gas dynamics and did not compute the emerging EM emission and light curves. In this paper, we focus on modelling the entire inspiral stage observable with LISA, extending our earlier work. We consider an equal-mass, non-spinning BH binary with a component mass of Mbh = 106 M⊙ and initial separation 60rg. This corresponds to the last ≈100 orbits of the merger, or the last ≈(1 + z)/2 week observed for a binary at redshift z. The binary is assumed to be on a circular orbit, and we allow the orbit to shrink due to GW emission all the way to merger. The binary is embedded in a viscous, optical thick α-disc with radiative cooling. We employ a pseudo-Newtonian potential for the black holes. Our spatial resolution is sufficient to resolve the individual BH's ISCO, allowing us to remove gas inside the event horizon physically. We analyse the dynamics and the time-dependent spectral signatures of the emerging thermal emission. This paper is organized as follows. In Section 2, we summarize our numerical methods. In Section 3, we present our main findings, including disc spectra and light curves. We also examine the periodicities of the light curves separately in the soft and hard X-ray bands. In Section 4, we discuss the implications of our results, along with some caveats. Finally, in Section 5, we summarize our conclusions and indicate topics for future work. 2 NUMERICAL SET-UP The model and numerical set-up of this work closely resembles that in our previous work (Farris et al. 2015b). We briefly outline the key aspects here and highlight the modifications we made for this study. We refer the reader to the above paper for further details. All simulations are performed using the publicly available moving-mesh grid code DISCO in 2D (Duffell & MacFadyen 2012; Duffell 2016). We simulate an equal-mass SMBHB with Mbin = 2Mbh = 2 × 106 M⊙ and initial separation a0 = 60GMbh/c2. At this separation, a typical 106 M⊙ binary at z < 1 can be localized by LISA on the sky to a few square degrees, giving the chance to a wide-field telescope to locate and identify the source (Kocsis et al. 2008; Lang & Hughes 2008). We assume the BHs are non-spinning so that a0 = 10risco = 30rS, where rS is the Schwarzschild radius of an individual BH. The BHs are placed on the grid, and the vertically integrated fluid equations are solved assuming an α-viscosity with α = 0.1. Inside the cavity that forms around the binary (with a size roughly twice the initial binary separation), the grid spacing is roughly linear in radius and risco is resolved with ≈8 radial cells; there are ≈ 100 cells in total inside the ISCO. Outside the cavity, the radial grid spacing increases approximately logarithmically up to 30a0, allowing us to concentrate computational resources on the flow in the inner regions of the disc. The total number of radial cells is 512. At all radii, we chose our azimuthal grid spacing such that the aspect ratio of each cell is kept approximately unity, unless this requires a number of azimuthal grid cells Nϕ > 512, in which case we cap the number of ϕ-zones and set Nϕ = 512. To mimic the effect of black hole accretion, gas inside the event horizon (roughly the innermost 2 or 3 radial cells around each BH) is depleted effectively instantaneously at each time-step using a sink prescription similar to Tang et al. (2017). We use a sink radius rsink = rS and a sink time-scale ≈10−4 orbital period. In this paper, we assume the mass and orbit of the binary black hole are not affected by the accretion on the time-scale of our simulations. Our simulation is pseudo-Newtonian, i.e. we model the black holes as point masses with the modified potential (Paczyńsky & Wiita 1980)   $$\psi (r)=-GM/(r-r_{\rm S}).$$ (1)Since the Paczynski–Wiita potential diverges at r = rS, we employ an upper limit for the gravitational force: |f(r)| = |f(1.5rS)| for r < 1.5rS. We choose 1.5rS because our cell size is ≈0.5rS. We use a γ-law equation of state for the gas in the form P = (γ − 1)e, where P and e are the vertically integrated pressure and internal energy, respectively, and γ = 5/3 is the adiabatic index. Radiative cooling is incorporated naturally through the energy equation, as described in Farris et al. (2015a). We assume the disc is optically thick and geometrically thin (the optical depth due to electron scattering at r ∼ a near the time of decoupling is ∼105; Haiman et al. 2009). These assumptions are consistent with our simulation results. The cooling rate is $$q_{\rm cool}=\frac{4\sigma }{3\tau _{{\rm e}}} T^{4}$$, where T is the mid-plane temperature and τe is the optical depth for electron scattering (τe = Σσt/mp). For a Shakura-Sunyaev disc around a single BH, the cooling is in balance with the local viscous heating rate, qvis = 9/8αPΩ. Following Farris et al. (2015a) and Farris et al. (2015b), we parametrize our simulations by specifying the Mach number at the fixed initial radius r = a0, $$\mathcal {M}_{a}=v_{a}/\sqrt{P(a)/\Sigma (a)}$$, where va ≡ (GMbin/a)1/2 and Σ is the surface density, We then obtain the local cooling rate from the scaling   $$q_{\rm cool}=\frac{9}{8}\alpha \Omega (a)P(a)\mathcal {M}^{8}_{a}\left({\frac{P}{\Sigma v^{2}_a}}\right)^{4}\left(\frac{\Sigma }{\Sigma (a)}\right)^{-1},$$ (2)where Ω is the Keplerian angular frequency.1 Changing $$\mathcal {M}_{a}$$ is equivalent to changing the thickness h/r of the disc, and in this work we fix $$\mathcal {M}_{a}=(h/r)^{-1}=10$$. A standard steady-state quasar disc (Shakura & Sunyaev 1973) with a near-Eddington accretion rate $${\sim } \dot{M}_{\rm Edd}$$ has a scale height of h/r ∼ 10−3, and is much thinner than the disc in our simulations (h/r ∼ 0.1). Simulating an accretion disc with h/r ∼ 10−3 is prohibitively challenging, and our choice is typical of values employed in previous numerical works (e.g. D'Orazio, Haiman & MacFadyen 2013; Farris et al. 2014; Farris et al. 2015b; Ragusa, Lodato & Price 2016; Miranda, Muñoz & Lai 2017 have included slightly smaller h/r values ranging from 0.02 to 0.13) and allows for a direct comparison with those works. The optical thickness to electron scattering (τes > 10) and free–free absorption (τff > 104) are both very high in our simulation domain (see e.g. fig. 10 in Corrales, Haiman & MacFadyen 2010). Therefore, throughout this paper, we compute the spectrum assuming multicolour thermal body radiation   $$L_{\nu }=\int \frac{2h\nu ^{3}}{c^{2}(\exp \left[\frac{h\nu }{kT_{\rm eff}(r, \phi )}\right]-1)}{\rm d}A,$$ (3)where Teff(r, ϕ) is the effective temperature obtained from the local cooling rate qcool(r, ϕ) in the cell at (r, ϕ) via Teff = (qcool/σ)0.25 with σ the Stefan–Boltzmann constant. While computing the observed spectrum and light curve, we additionally consider the relativistic Doppler effect and the gravitational redshift. For simplicity, and to assess the maximal impact of the Doppler effect, we further assume the circumbinary disc is observed edge-on. Photons emitted by gas suffer a Doppler shift in frequency D = [Γ(1 − β∥)]−1, where Γ = (1 − β2)−1/2 is the Lorentz factor, and β∥ the component of the velocity along the line of sight. The apparent flux Fν at a fixed observed frequency ν is modified from the flux of a stationary source $$F_\nu ^{\rm 0}$$ to $$F_\nu =D^{3}F^{\rm 0}_{\nu /D}$$. For the gravitational redshift, we employ D = (1 − rS/Re)−0.5, where rS is the Schwarzschild radius and Re is the distance between the source of emission and the BH. The above is an approximation to a fully relativistic ray-tracing in the binary's dynamical metric (Schnittman et al. 2017). In particular, we ignore gravitational self-lensing (D'Orazio & Di Stefano 2017; Haiman 2017) and other relativistic effects, such as light-travel time modulations (i.e. Shapiro delay), which appear at order v/c. Our disc has an initial surface density profile Σ(r) = Σ0(r/a)−0.5 and contains an artificial cavity with the density depleted inside the radius 2.5a0, as in Tang et al. (2017). Viscosity in the circumbinary disc transports angular momentum outwards, causing the cavity to refill. We first hold the binary on the initial circular orbit and evolve the system for several hundred orbits to reach a quasi-steady state. We then allow the binary separation to shrink due to quadrupole formula for GW emission (Peters 1964),   $$a(t)=a_0(1-t/\tau )^{0.25}.$$ (4)Here, τ = 7.22 d is the GW inspiral time from the initial separation of 30 GMbin/c2. Finally, we employ a set of convenient code units, in which the initial binary period is T(0) = 2π and a0 = 1; as a result, τ = 385.17 in code units. 3 RESULTS As stated above, before allowing the black holes to inspiral and merge, we first hold the binary separation fixed, and evolve the system for several hundred orbits to reach a quasi-steady state. This quasi-steady state is similar to that in previous works (D'Orazio et al. 2013; Farris et al. 2014, 2015a; Tang et al. 2017) in that the cavity becomes lopsided, and two narrow streams periodically connect the minidiscs around the individual BHs with the circumbinary disc. Similar to Farris et al. (2015b), we find that these features persist until right before the merger, with only a gradual decline in the size and surface density of the minidiscs and the streams. In Fig. 1, we show the 2D distributions of surface density and effective temperature at three different snapshots, corresponding to t = 0 (the beginning of the inspiral), 0.5τ (half-way to merger), and 0.99τ (just before the merger). Figure 1. View largeDownload slide The top row shows snapshots of the surface density Σ/Σ0, where Σ0 is the initial surface density at r = a, and the bottom row shows snapshots of the effective temperature Teff in [eV]. The horizontal and vertical axes are in units of the initial binary separation a0. The snapshots are taken, from left to right, at the times t = 0.0τ, 0.5τ, and 0.99τ, where τ = 7.22 d corresponds to the merger. The cavity size decreases and the BHs accrete efficiently during the entire inspiral. Figure 1. View largeDownload slide The top row shows snapshots of the surface density Σ/Σ0, where Σ0 is the initial surface density at r = a, and the bottom row shows snapshots of the effective temperature Teff in [eV]. The horizontal and vertical axes are in units of the initial binary separation a0. The snapshots are taken, from left to right, at the times t = 0.0τ, 0.5τ, and 0.99τ, where τ = 7.22 d corresponds to the merger. The cavity size decreases and the BHs accrete efficiently during the entire inspiral. A spectrum computed at the beginning of the inspiral is also shown in Fig. 2. In this figure, we additionally show three distinct components of the spectrum, arising from outer disc (r > 3a0), the shocked streams inside the cavity and shocked hot gas near the cavity wall (a0 < r < 3a0) and minidiscs (r < a0). At frequencies below ∼1keV, the emission from the viscously heated outer disc dominates, while at higher frequencies, the emission mostly comes from the shock-heated gas in the minidiscs, streams, and near the cavity wall. At frequencies between ∼(1–20) keV, the total spectrum shows a plateau with a small depression near ∼3 keV, resembling the ‘notch’ discussed in Roedig, Krolik & Miller (2014). The spectrum computed in Farris et al. (2015a) has a larger magnitude at the highest frequencies, and does not show a plateau. This is because in our simulation we resolve the event horizon, and gas near the black hole is efficiently cleared out, as it falls rapidly into the BH. In Farris et al. (2015a), a much slower sink prescription was implemented, causing hot, shock-heated gas to accumulate around each black hole, and brightening the high-frequency emission. Figure 2. View largeDownload slide Thermal spectrum computed from a simulation snapshot at the beginning of the inspiral (t = 0). The full spectrum is shown by the blue curve. The distinct components arising from gas in the minidiscs (r < a), streams and cavity wall (a < r < 3a) and from the outer regions of the circumbinary disc (r > 3a) are represented by the orange dashed line, green dotted line, and red dash–dotted line, respectively. Figure 2. View largeDownload slide Thermal spectrum computed from a simulation snapshot at the beginning of the inspiral (t = 0). The full spectrum is shown by the blue curve. The distinct components arising from gas in the minidiscs (r < a), streams and cavity wall (a < r < 3a) and from the outer regions of the circumbinary disc (r > 3a) are represented by the orange dashed line, green dotted line, and red dash–dotted line, respectively. In Fig. 3, we show the total spectrum at three different snapshots, t = 0 (blue; same as in Fig. 2), 0.8τ (orange), and 0.99τ (green). The low-frequency part of the spectrum, which is dominated by the outer regions of the circumbinary disc, shows little change throughout the simulation. By comparison, the higher frequencies (≳ 1 keV), which arises from the minidiscs or streams in the cavity, declines noticeably by 0.8τ; the decline in the hard X-ray band becomes significant during the end of inspiral (0.99τ). Note that in this paper we chose Mbin = 2 × 106 M⊙, so that the GW signal is in the LISA band, and 0.99τ corresponds to the last ≈2 h of the inspiral. Figure 3. View largeDownload slide Total composite thermal spectra computed from simulation snapshots at t = 0 (blue curve), t = 0.8τ (orange dashed), and t = 0.99τ (green dot–dashed). Figure 3. View largeDownload slide Total composite thermal spectra computed from simulation snapshots at t = 0 (blue curve), t = 0.8τ (orange dashed), and t = 0.99τ (green dot–dashed). The scaling relationships in our simulation can be calculated by equating the viscous heating qvis( ∝ αΩP) to the radiative cooling qcool( ∝ Σ−1T4). Note that we fix a0 = 60GMbh/c2, so that the orbital velocity is a constant, independent of Mbin, and the angular velocity $$\Omega \propto M_{\rm bh}^{-1}$$. Since we have further fixed $$\mathcal {M}_{\rm a}=10$$, the mid-plane temperature T and ratio P/Σ are also fixed, i.e. P ∝ Σ. Finally, with α fixed as well, we have qvis ∝ Σ/Mbh and qcool ∝ Σ−1. This then yields the scaling $$\Sigma _{0}\propto M_{\rm bh}^{1/2}$$ and the effective temperature Teff ∝ Σ−1/4 scales as $$M_{\rm bh}^{-1/8}$$. Note that this implies that the photon frequency scales with total binary mass as $$\nu \propto T_{\rm eff}\propto M_{\rm bin}^{-1/8}$$, as well. For a steady disc around a solitary BH, it would further scale with Mach number as $$h\nu \propto \mathcal {M}_{\rm a}^{-5/4}$$ ($$\Sigma \propto \mathcal {M}^{-3}$$ and $$T\propto \mathcal {M}^{-2}$$). This suggests that for a thinner disc, the frequencies could be lower than shown for our fiducial $$\mathcal {M}_{a}=10$$ in Figs 2 and 3. However, for gas near the cavity wall and minidiscs, the variation of the disc thickness caused by its interaction with the binary BHs requires further study (e.g. Ragusa et al. 2016). In Fig. 4, we show h/r of our simulated disc. In the shock-heated region, h/r is increased to ∼0.2, while in the outer region h/r stays at ∼0.1. Figure 4. View largeDownload slide Local disc thickness h/r throughout the disc at t = 0. Figure 4. View largeDownload slide Local disc thickness h/r throughout the disc at t = 0. In Fig. 5, we illustrate the effects of the relativistic Doppler shift and gravitational redshift on the observed spectra. As mentioned in Section 2, we assume the disc is observed edge-on, so that the Doppler effect is maximized. Note the spectra in Fig. 5 do not include the inclination angle. The observed spectrum will be dimmer by a factor of cos (θ), where θ is the angle of inclination to the line of sight. We include four viewing angles in the plane of the disc, corresponding to the $$\pm \hat{x}$$ and $$\pm \hat{y}$$ directions. As the figure shows, at the low-frequency end (hν ≲ 1 keV), the spectrum is not modified by either the Doppler or gravitational shifts. At higher frequencies (1 keV ≲ hν ≲ 20 keV), the Doppler shift causes an overall dimming (i.e. red and blue curves), while at the highest frequencies (hν ≳ 20 keV), the spectrum is brightened significantly. Interestingly, we find that the dimming and the brightening retain their sign in all four different viewing directions; we offer a simple interpretation of this result in Section 4 below. For reference, the figure also shows a case with only the GR redshift (green curve), which is independent of viewing angle, and always reduces the observed flux. Figure 5. View largeDownload slide Spectra illustrating the impact of the special relativistic Doppler shift and the gravitational red/blue shift. For reference, the black curve shows the spectrum without considering either of these effects. The green curve shows the spectrum including only the gravitational redshift, while the other four spectra, shown in blue and red, include both relativistic effects (for the four different viewing directions $$\pm \hat{x}$$ and $$\pm \hat{y}$$, as labelled in the inset). All spectra are calculated from the snapshot at t = 0, when the orbital velocity of the BHs is 0.091c. Figure 5. View largeDownload slide Spectra illustrating the impact of the special relativistic Doppler shift and the gravitational red/blue shift. For reference, the black curve shows the spectrum without considering either of these effects. The green curve shows the spectrum including only the gravitational redshift, while the other four spectra, shown in blue and red, include both relativistic effects (for the four different viewing directions $$\pm \hat{x}$$ and $$\pm \hat{y}$$, as labelled in the inset). All spectra are calculated from the snapshot at t = 0, when the orbital velocity of the BHs is 0.091c. In order to study the variability of the spectrum at different frequencies, we calculate the light curves in the soft (at 2 keV) and hard (at 10 keV) X-ray bands. In Fig. 6, we show the light curve at the beginning of the simulation and right before merger, with and without the Doppler effect. In the 2 keV band, we find clear periodicity in the light curve throughout the inspiral. In the 10 keV band, the inspiral light curve at the beginning is much noisier. We can still see some periodicity by eye, but it is harder to distinguish than in the 2keV case. However, in the late stages of the inspiral, the 10 keV band also shows clear periodicity at the orbital period. The light curves also show an overall see-saw modulation, with a period approximately eight times longer than torb. Furthermore, the light curves show a clear a behaviour analogous to the GW chirp signal, in that the periodicity tracks the binary's increasing orbital frequency. Figure 6. View largeDownload slide Gravitational wave strain (top panels) and X-ray light curves (middle and bottom) calculated at the beginning (left) and end (right) of the simulated inspiral in the LISA band. The top panel shows the GW strain observed edge-on from a binary at redshift z = 1. The middle panel corresponds to the 2 keV band and the bottom panel to 10 keV band. The blue curves show the light curves neglecting the relativistic Doppler shift or gravitational redshift. The orange curves show the light curves including these effects. While the Doppler effect imposes a small additional periodic modulation due to the binary's orbital motion, its dominant effect is a nearly uniform dimming of the light curve (this results from the asymmetry between the large blue and redshifts due to the internal motion of gas in the minidiscs; see the text for discussion). The x-axis shows the look-back time in the source's rest frame, and the number of binary orbits (in parentheses). Figure 6. View largeDownload slide Gravitational wave strain (top panels) and X-ray light curves (middle and bottom) calculated at the beginning (left) and end (right) of the simulated inspiral in the LISA band. The top panel shows the GW strain observed edge-on from a binary at redshift z = 1. The middle panel corresponds to the 2 keV band and the bottom panel to 10 keV band. The blue curves show the light curves neglecting the relativistic Doppler shift or gravitational redshift. The orange curves show the light curves including these effects. While the Doppler effect imposes a small additional periodic modulation due to the binary's orbital motion, its dominant effect is a nearly uniform dimming of the light curve (this results from the asymmetry between the large blue and redshifts due to the internal motion of gas in the minidiscs; see the text for discussion). The x-axis shows the look-back time in the source's rest frame, and the number of binary orbits (in parentheses). In order to identify the origin of the periodicity in the light curves, in Fig. 7 we decompose them into emission arising from r < a (dominated by gas in/near the minidiscs) and from r > a (dominated by the circumbinary disc). This figure reveals that the emission from the circumbinary disc is clearly periodic from the beginning of the simulation all the way to the merger in both bands. This is because each BH periodically interacts with the cavity wall and creates a stream which hits the cavity wall, shock-heats the gas, and leads to enhanced X-ray emission. By comparison, the emission from the minidiscs is noisy in the beginning, and develops clear periodicity only closer to the merger. One explanation for this phenomenon is that in the earlier stage the minidiscs act as buffers of the accretion streams. Unlike steady-state α discs, the internal shocks in the minidiscs create significant noise; their overall shape often differs significantly from a circular disc. In the late stage of the inspiral (in the last 60 h before the merger), as the truncation radius decreases, this buffer effect declines, and the narrow accretion streams eventually connect the cavity wall and event horizon almost directly (see the right-hand panel of Fig. 1). Therefore, the noise produced by the minidiscs gradually disappears. Figure 7. View largeDownload slide Decomposition of the light curves. The orange and blue curves correspond to emission arising from r < a (dominated by gas in/near the minidiscs) and from r > a (dominated by the circumbinary disc). The green curves correspond to the total emission. The upper two panels show the light curves in the 2 keV band, and the bottom two panels in the 10 keV band. For clarity of display, the orange curves, which show sub-dominant emission from the minidiscs in the soft band in the upper two panels, use the y-axis labels on the right side. Figure 7. View largeDownload slide Decomposition of the light curves. The orange and blue curves correspond to emission arising from r < a (dominated by gas in/near the minidiscs) and from r > a (dominated by the circumbinary disc). The green curves correspond to the total emission. The upper two panels show the light curves in the 2 keV band, and the bottom two panels in the 10 keV band. For clarity of display, the orange curves, which show sub-dominant emission from the minidiscs in the soft band in the upper two panels, use the y-axis labels on the right side. In the 2keV band, emission from the circumbinary disc dominates, and the light curve therefore also shows very clear periodicity, during the entire inspiral process. The opposite situation occurs in the 10keV band, in which the total emission is dominated by the minidiscs. The light curve is noisy in the beginning, and shows much clearer periodicity in the end stage. From the bottom-most panel of Fig. 7, we see that the emission from the minidiscs and circumbinary discs has the same period, but with different phases. Therefore, in the later stages, when the strength of the emission from these two regions is comparable, the total light curve sometimes shows two maximas in a single period. Haiman (2017) and D'Orazio & Haiman (2017) proposed a toy model for the minidiscs as two ‘light bulbs’, both of which are steady in their own rest frame. Therefore, unless observed perfectly face-on, the Doppler modulation will lead to periodicity at the orbital frequency. Note that the light curve would still be modulated, even in the case of an equal-mass binary, because the red and blueshifts of the thermal emission from the two BHs, with equal but opposite line-of-sight velocities, will not exactly cancel. However, in our simulation, we find that the internal structure of the minidiscs is far from stationary, with shocks forming and streams of gas entering and leaving the fiducial Hill sphere. The emission from the minidisc regions is therefore quite noisy. As a result, the Doppler modulation may be difficult to extract from the observed light curve. We expect that this may still be feasible, given a GW template from LISA, with which the X-ray light curve can be cross-correlated, but this requires further analysis. Fortunately, we find that the shock-heated cavity wall provides a clear periodicity throughout the entire inspiral process. Moreover, since this periodicity does not originate from the Doppler effect (see Fig. 6), it should exist independently of the inclination angle. In order to assess how well the periodicity in the light  curves track the binary's orbital period, we rescale the light-curve time series with the evolving binary orbital period. In Fig. 8, we plot the Lomb–Scargle periodograms with and without this time rescaling. The rescaled periodograms show peaks at ≈2 × Ωorb, both at the beginning and the end of the inspiral, implying that this periodicity is able to follow the accelerating binary orbital frequency (i.e. the ‘chirp’) throughout the entire inspiral process. The first peak in the periodogram without rescaling arises from a dense lump in the lopsided cavity wall, as reported in a series of previous works (e.g. MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Krolik 2010; Roedig et al. 2011; Noble et al. 2012; Shi et al. 2012; D'Orazio et al. 2013; Farris et al. 2014; Bowen et al. 2017a). This feature is smeared out and disappears after time rescaling, because it does not follow the binary's frequency chirp. More specifically, the frequency of this low-frequency peak increases from ≈0.25Ωorb(0) to ≈0.4Ωorb(0), but this increase is slower than the evolving binary orbital frequency Ωorb(t). Just prior to the merger, the period of this low-frequency modulation is around 15 orbits. Similarly, in Fig. 1, we see that the overall size of the cavity is shrinking, but cannot keep up with the binary separation, which decreases more rapidly. Nevertheless, this low-frequency modulation persists throughout the inspiral process, as we can clearly see in Fig. 6. In principle, the evolution of this low-frequency modulation could be measured in the X-ray light curves, and would probe the gas morphology near the central cavity. Figure 8. View largeDownload slide Lomb–Scargle periodograms of the light curves. The four panels correspond to those in Fig. 6, namely to the 2 and 10 keV bands (top and bottom, respectively) and to the beginning and end of the simulated inspiral (left and right, respectively). In each panel, the red curve corresponds to the periodogram of the original light curves shown in Fig. 6. The blue curves correspond to scaled light curves, in which time is measured in units of the instantaneous binary orbital period. Figure 8. View largeDownload slide Lomb–Scargle periodograms of the light curves. The four panels correspond to those in Fig. 6, namely to the 2 and 10 keV bands (top and bottom, respectively) and to the beginning and end of the simulated inspiral (left and right, respectively). In each panel, the red curve corresponds to the periodogram of the original light curves shown in Fig. 6. The blue curves correspond to scaled light curves, in which time is measured in units of the instantaneous binary orbital period. 4 DISCUSSION One of the motivations of this work was to compute the modulation of the light curves due to relativistic Doppler effect. Naively, such a modulation inevitably arises from the same orbital motion that produces GWs, and therefore the GW and X-ray chirp signals should track one another with a known phase (Haiman 2017). This would be of particular interest for LISA, allowing a robust measurement of the speed difference of photons and gravitons, and more generally, a secure identification of the EM counterpart of a LISA binary. Our results suggest that it may be difficult to extract this Doppler modulation from the X-ray light curve. We find efficient fuelling of both BHs through the inspiral process, with gas funnelled inside each BH's Hill radius, forming tidally truncated ‘minidiscs’. These minidiscs persist nearly all the way to the merger, until the tidal truncation radius shrinks to a size comparable to that of the ISCOs (roughly 60 h prior to merger). However, the effective temperature, density, and velocity fields of the gas inside these minidiscs are not stationary, and the net Doppler shift is dominated by these large ‘internal’ gas motions, rather than the orbital velocity of the BHs. Interestingly, as Figs 5 and 6 show, the most prominent effect of the Doppler shift in the 2–10 keV range is an overall dimming of the light curve, rather than a sinusoidal modulation. This is because the hot gas dominating this emission is located very close to each BH (near or inside its Hill radius), and the velocity of this gas is larger than the orbital speed of the BHs. Regardless of the viewing direction, there is always gas with both large blue and redshifts along the line of sight, and the total composite spectrum includes both of these shifts. The net result can be either an overall dimming or brightening, depending on the local slope and curvature of the spectrum. In other words, for the steeply curved blackbody spectral shape, the Doppler effect is not symmetric under switching the sign of the line-of-sight velocity. Fig. 5 shows that in the 2–10 keV range, the redshift and the corresponding dimming from receding gas dominate over the blueshift and brightening from approaching gas. We emphasize that this result is sensitive to the spectral shape. One caveat is that, as mentioned above, due to numerical constraints our simulated disc is thicker than a typical active galactic nucleus (AGN) disc of a solitary BH. If we were to scale the simulation to correspond to a higher Mach number, this would shift the normalization of the photon energy. For example, the soft versus hard X-ray bands could become FUV and EUV bands for a disc with a Mach number of $$\mathcal {M}_{a} \sim 1000$$, instead of the $$\mathcal {M}_{a}=10$$ we adopted (recall that $$h\nu \propto \mathcal {M}_{\rm a}^{-5/4}$$ for single black hole accretion disc). On the other hand, the gas in the inner regions is being heated primarily by shocks, rather than viscosity, and we find that h/r in the inner regions of the flow remains ∼ 0.2 throughout the simulation, justifying the scaling to lower Mach numbers (see Fig. 4). However, it requires further study to understand how this depends on the choice of $$\mathcal {M}_{a}$$. A higher $$\mathcal {M}_{a}$$ could lead to a larger cavity and reduce the gas fuelling from the cavity wall (Ragusa et al. 2016). This may reduce the luminosity of the minidiscs. We also note that in equation (2), we assumed that the fluid is gas-pressure dominated; this assumption is likely violated in the shock-heated streams and minidiscs. Another caveat in our conclusions is that we focused on thermal emission. We have verified that the effective absorptive opacity, defined as τeff = [τabs(τabs + τscat)]1/2, where τabs is the optical depth to free–free absorption and τscat is the optical depth to electron scattering, remains τeff ∼ 104 ≫ 1 throughout our simulation. This justifies considering a thermalized spectrum. However, empirically, the X-ray emission from an AGN consists of a combination of thermal disc emission and a hot corona, and in this paper, we have not included the latter component. If the emission from a corona turns out to be more stable (i.e. if the corona lacks large internal bulk motions), the Doppler modulation in the X-ray band may be dominated by the BH's orbital motion, and much more prominent. Moreover, in this paper, we have focused on an equal-mass binary. We expect that the Doppler modulation would be more prominent in the case of low-mass ratio binaries, because the circum-primary and -secondary minidiscs would be more stable, and there would be less cancellation between their emission. The overall shape of the light curves we found show a clear ‘X-ray chirp’, caused by the interaction between the binary and the gas discs. The periodicity of the chirp tracks the shrinking binary orbital period torb, in tandem with the GW signal. Since the chirp is produced by hydrodynamical modulation, its absolute phase is not tied directly to the GW emission. Nevertheless, the speed difference between photons and gravitons should still be measurable, in principle, given a precise measurement of the phase-evolution of both the GW and the X-ray chirp signals. In the future, it would be interesting to compute the accuracy to which the phase difference can be measured in practice, if the absolute phase is unknown ab-initio, given realistic observational S/N from LISA, and from an X-ray instrument such as Athena,2 Lynx,3 ULTRASAT,4 or EXTP (Zhang et al. 2016). In this paper, we neglect spin of the BHs, but the spin would change the GW signal significantly, and would further complicate a direct comparison of the X-ray light curve with the LISA signal. The ISCO for a spinning black hole is smaller than a non-spinning one, which may enhance the luminosity at high frequencies. We have made several simplifying assumptions throughout this paper, which we intend to relax in our future work. We have focused here on an equal-mass binary, but we intend to study the dependence of the spectra and light curves on the binary mass ratio. Similarly, we intend to determine the dependence of these features on the disc thickness. We also used an α-prescription for the viscosity. We intend to perform magnetohydrodynamic (MHD) simulations in the future (e.g. Krolik 2010). In this work, we have also assumed that the binary BHs are on circular orbits, with zero eccentricity. On the other hand, eccentricity is likely to develop as a result of strong gas torques during the earlier stages of the inspiral (e.g. Cuadra et al. 2009; Roedig et al. 2011; Sesana, Gualandris & Dotti 2011; Roedig et al. 2012; Miranda, Muñoz & Lai 2017). In the future, we intend to study live binaries, and follow the development of eccentricity, and to assess the impact of eccentricity on the accretion rates and the predicted spectra and light curves. In the present simulation, we followed a pseudo-Newtonian approach; a full relativistic simulation of the gas in the vicinity of the event horizon could modify our results. Finally, we started our current simulations somewhat past the nominal ‘decoupling’ stage. In the cartoon picture of the binary inspiral (e.g. Milosavljevic & Phinney 2005), once the binary is sufficiently compact that the GW inspiral time-scale is shorter than the viscous time-scale in the nearby disc, the BHs outpace the disc and leave the disc gas behind. The decoupling separation adec can be calculated with tgw = tvis, where $$t_{\rm gw}=\frac{5c^{5}a^{4}}{64G^{3}m_{1}m_{2}(m_{1}+m_{2})}$$ is the GW inspiral time-scale (Peters 1964) and $$t_{\rm vis}=\frac{2(2a)^2}{3\nu }$$ is the viscous time in an unperturbed disc at the radius r = 2a. With the parameters (α = 0.1; $$\mathcal {M}_{a}=10$$) adopted in our simulation, we find the decoupling separation adec = 128GMbh/c2, which is somewhat larger than our initial separation of 60GMbh/c2. In this picture, our procedure of holding the binary's orbit initially fixed, and allowing it to reach a steady state would not be justified. However, we first note that even in a simple 1D toy model, the decoupling radius can only be taken as an order-of-magnitude estimate. For example, torques modify the density profile outside the binary's orbit, causing a pile-up, increasing the pressure, steepening the radial gradients, shortening the viscous time, and decreasing adec (e.g. Milosavljevic & Phinney 2005; Kocsis, Haiman & Loeb 2012). More importantly, from the results we find here, as well as from our earlier work (Farris et al. 2015b), we see that the cartoon picture is inaccurate, and the cavity can, in fact, shrink and follow the inspiraling binary, even well past the nominal decoupling. This is because the gas dynamics near the black holes is dominated by gravity and shocks, and not by viscosity, and the gas configuration is neither steady nor axisymmetric. More specifically, in Farris et al. (2015b), a simulation was started at the separation corresponding to tgw = 5tvis, i.e. safely before the nominal decoupling. Nevertheless, we found that the gas flows inwards, following the binary, and the accretion rate is only halved, even well past the decoupling, 0.1tvis before the merger. Similarly, in our simulation here, the spectrum is only slightly modified even at 0.8τ (see Fig. 3). We therefore expect our way of setting up the initial condition to be reasonably accurate. However, one caveat to note here is that the nominal decoupling separation is larger for a disc with smaller aspect ratios h/r (i.e. large Mach numbers). This means that the reduction of accretion rate and luminosity in more realistic, thinner discs may be more significant than the discs in our simulations. We intend to study this in future work, by beginning simulation near or prior to the decoupling stage and with larger Mach numbers. 5 CONCLUSIONS We have performed 2D simulations of accretion on to an inspiraling and merging binary BH system, and studied the corresponding EM signatures. We have computed the evolving multicolour blackbody emission from snapshots of our simulations, including the effects of the relativistic Doppler shift and gravitational redshift. The most important conclusions of this work can be summarized as follows: The light curves show clear periodicity. Accretion on to the individual BHs has long been known to be periodic (Artymowicz & Lubow 1994), and it is not surprising that any luminosity produced by gas near the individual BHs is periodic. However, we have found that the emission arising from gas farther out, near the cavity wall, which dominates the softer X-ray bands, shows a similar (and even clearer) periodicity. This is because streams of gas, flung outwards by the BHs, periodically hit and shock the cavity wall. The frequency is twice the binary's instantaneous orbital frequency. Distinct behaviour in different bands. In the beginning of the simulation, the emission from the minidiscs is noisy, but develops clear periodicity during the later stages. In the soft X-ray band, the emission from the circumbinary disc dominates, and the light curve shows a clear periodicity from the beginning of the simulation (a = 60GM/c2). By comparison, in the hard X-ray band, the emission in the beginning is dominated by the minidiscs, and periodicity in the early-inspiral light curve is obscured by noise, with clearer periodicity developing only in the late stages. Doppler modulation is sub-dominant. We find that the Doppler modulation of the light curve does not follow a simple toy model of two moving light bulbs in orbit. This is because the effective temperature, density, and velocity field of the gas in the minidiscs are not stationary, and the line-of-sight velocity of this gas is dominated by the large internal gas motions inside the minidiscs, rather than the clean orbital motion of the BHs. We find that the most conspicuous effect of the Doppler shift is an overall dimming (1keV≲ hν ≲ 20keV) or brightening (hν ≳ 20keV) of the light curves. We caution that these results are sensitive to the spectral shape and emission mechanism. In particular the Doppler effect could more closely track the binary's orbital motion for X-ray emission from more stable coronae around the individual BHs which do not have large internal motions. The EM chirp follows the GW inspiral. We have found that the minidiscs persist until the very end of the inspiral process. The discs are truncated at the gradually shrinking size of the Hill radius. The thermal emission from the minidiscs declines gradually, but remains overall significant all the way to the merger. The overall periodicities in both the soft and hard X-ray bands track the shrinking orbital period of the binary. No ‘second decoupling’. Fontecilla, Chen & Cuadra (2017) suggested that the inner (circum-primary) disc will become geometrically thick (h/r ∼ 1) in the late stages of the inspiral, as a result of tidal heating. They find, using 1D analytical models and numerical simulation, that the viscous time in the inner disc decreases to tvis ∼ tGW, and it effectively decouples from the binary. By comparison, in our 2D simulation, the size and surface density of the two minidiscs are continuously adjusted, and eventually, in the late stages, the minidiscs disappear, and narrow accretion streams instead fall inside the event horizon directly. The inspiral process does thicken the minidiscs, but the aspect ratios remain below h/r ∼ 0.2. The simulations presented here are based on many simplifications. Nevertheless, our results suggest that an EM chirp signal may accompany the GW emission, with the EM period tracking that of the GW signal right up to the merger of an SMBH binary. A detection of this EM chirp could provide a potential ‘smoking gun’ evidence to uniquely identify the EM counterpart of an SMBH binary detected by LISA. Acknowledgements We thank Daniel D'Orazio, Paul Duffell, Julian Krolik, and Geoffrey Ryan for useful discussions. Financial support was provided by NASA through ATP grant NNX15AB19G and ADAP grant NNX17AL82G and by NSF grants 1715661 and 1715356. ZH also gratefully acknowledges support from a Simons Fellowship in Theoretical Physics and hospitality by NYU during his sabbatical leave when this work began. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. Footnotes 1 We note there is a typographical error in equation (2) of Farris et al. (2015a); the correct exponent of the Mach number is $$\mathcal {M}_a^{8}$$, as in equation (2) above. 2 See http://www.cosmos.esa.int/web/athena. 3 See http://www.astro.msfc.nasa.gov/lynx. 4 See http://www.weizmann.ac.il/ultrasat. REFERENCES Abbott B. P. et al.  , 2017, ApJ , 848, L13 https://doi.org/10.3847/2041-8213/aa920c CrossRef Search ADS   Amaro-Seoane P. et al.  , 2017, preprint (arXiv:1702.00786) Artymowicz P., Lubow S. H., 1994, ApJ , 421, 651 https://doi.org/10.1086/173679 CrossRef Search ADS   Bode T., Haas R., Bogdanović T., Laguna P., Shoemaker D., 2010, ApJ , 715, 1117 https://doi.org/10.1088/0004-637X/715/2/1117 CrossRef Search ADS   Bowen D. B., Mewes V., Campanelli M., Noble S. C., Krolik J. H., Zilhao M., 2017a, ApJ , 853, L17 https://doi.org/10.3847/2041-8213/aaa756 CrossRef Search ADS   Bowen D. B., Campanelli M., Krolik J. H., Mewes V., Noble S. C., 2017b, ApJ , 838, 42 https://doi.org/10.3847/1538-4357/aa63f3 CrossRef Search ADS   Corrales L. R., Haiman Z., MacFadyen A., 2010, MNRAS , 404, 947 https://doi.org/10.1111/j.1365-2966.2010.16324.x CrossRef Search ADS   Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, MNRAS , 393, 1423 https://doi.org/10.1111/j.1365-2966.2008.14147.x CrossRef Search ADS   D'Orazio D. J., Di Stefano R., 2017, MNRAS , 474, 2975 https://doi.org/10.1093/mnras/stx2936 CrossRef Search ADS   D'Orazio D. J., Haiman Z., 2017, MNRAS , 470, 1198 https://doi.org/10.1093/mnras/stx1269 CrossRef Search ADS   D'Orazio D. J., Haiman Z., MacFadyen A., 2013, MNRAS , 436, 2997 https://doi.org/10.1093/mnras/stt1787 CrossRef Search ADS   Dotti M., Sesana A., Decarli R., 2012, Adv. Astron. , 2012 Duffell P. C., 2016, ApJS , 226, 2 https://doi.org/10.3847/0067-0049/226/1/2 CrossRef Search ADS   Duffell P. C., MacFadyen A. I., 2012, ApJ , 755, 7 https://doi.org/10.1088/0004-637X/755/1/7 CrossRef Search ADS   Farris B. D., Gold R., Paschalidis V., Etienne Z. B., Shapiro S. L., 2012, Phys. Rev. Lett. , 109, 221102 https://doi.org/10.1103/PhysRevLett.109.221102 CrossRef Search ADS PubMed  Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ , 783, 134 https://doi.org/10.1088/0004-637X/783/2/134 CrossRef Search ADS   Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015a, MNRAS , 446, L36 https://doi.org/10.1093/mnrasl/slu160 CrossRef Search ADS   Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015b, MNRAS , 447, L80 https://doi.org/10.1093/mnrasl/slu184 CrossRef Search ADS   Fontecilla C., Chen X., Cuadra J., 2017, MNRAS , 468, L50 Giacomazzo B., Baker J. G., Miller M. C., Reynolds C. S., van Meter J. R., 2012, ApJ , 752, L15 https://doi.org/10.1088/2041-8205/752/1/L15 CrossRef Search ADS   Haiman Z., 2017, Phys. Rev. D , 96, 023004 https://doi.org/10.1103/PhysRevD.96.023004 CrossRef Search ADS   Haiman Z., Kocsis B., Menou K., 2009, ApJ , 700, 1952 https://doi.org/10.1088/0004-637X/700/2/1952 CrossRef Search ADS   Kocsis B., Haiman Z., Menou K., 2008, ApJ , 684, 870 https://doi.org/10.1086/590230 CrossRef Search ADS   Kocsis B., Haiman Z., Loeb A., 2012, MNRAS , 427, 2680 https://doi.org/10.1111/j.1365-2966.2012.22118.x CrossRef Search ADS   Kormendy J., Ho L. C., 2013, ARA&A , 51, 511 https://doi.org/10.1146/annurev-astro-082708-101811 CrossRef Search ADS   Krolik J. H., 2010, ApJ , 709, 774 https://doi.org/10.1088/0004-637X/709/2/774 CrossRef Search ADS   Lang R. N., Hughes S. A., 2008, ApJ , 677, 1184 https://doi.org/10.1086/528953 CrossRef Search ADS   MacFadyen A. I., Milosavljević M., 2008, ApJ , 672, 83 https://doi.org/10.1086/523869 CrossRef Search ADS   Milosavljevic M., Phinney E. S., 2005, ApJ , 622, L93 https://doi.org/10.1086/429618 CrossRef Search ADS   Miranda R., Muñoz D. J., Lai D., 2017, MNRAS , 466, 1170 https://doi.org/10.1093/mnras/stw3189 CrossRef Search ADS   Noble S. C., Mundim B. C., Nakano H., Krolik J. H., Campanelli M., Zlochower Y., Yunes N., 2012, ApJ , 755, 51 https://doi.org/10.1088/0004-637X/755/1/51 CrossRef Search ADS   Paczyńsky B., Wiita P. J., 1980, A&A , 88, 23 Peters P. C., 1964, Phys. Rev. , 136, B1224 https://doi.org/10.1103/PhysRev.136.B1224 CrossRef Search ADS   Phinney E. S., 2009, Astro2010: The Astronomy and Astrophysics Decadal Survey, Science White Papers , National Research Council (2010) Washington, 235 Ragusa E., Lodato G., Price D. J., 2016, MNRAS , 460, 1243 https://doi.org/10.1093/mnras/stw1081 CrossRef Search ADS   Roedig C., Dotti M., Sesana A., Cuadra J., Colpi M., 2011, MNRAS , 415, 3033 https://doi.org/10.1111/j.1365-2966.2011.18927.x CrossRef Search ADS   Roedig C., Sesana A., Dotti M., Cuadra J., Amaro-Seoane P., Haardt F., 2012, A&A , 545, A127 https://doi.org/10.1051/0004-6361/201219986 CrossRef Search ADS   Roedig C., Krolik J. H., Miller M. C., 2014, ApJ , 785, 115 https://doi.org/10.1088/0004-637X/785/2/115 CrossRef Search ADS   Schnittman J. D., Dal Canton T., Camp J., Tsang D., Kelly B. J., 2017, ApJ , 853, 123 https://doi.org/10.3847/1538-4357/aaa08b CrossRef Search ADS   Sesana A., Gualandris A., Dotti M., 2011, MNRAS , 415, L35 https://doi.org/10.1111/j.1745-3933.2011.01073.x CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Shi J.-M., Krolik J. H., Lubow S. H., Hawley J. F., 2012, ApJ , 749, 118 https://doi.org/10.1088/0004-637X/749/2/118 CrossRef Search ADS   Tang Y., MacFadyen A., Haiman Z., 2017, MNRAS , 469, 4258 https://doi.org/10.1093/mnras/stx1130 CrossRef Search ADS   Will C. M., 2006, Living Rev. Relativ. , 9, 3 https://doi.org/10.12942/lrr-2006-3 CrossRef Search ADS PubMed  Zhang S. N. et al.  , 2016, in Jan-Willem A., den H., Tadayuki T., Marshall B., eds, Proc. SPIE Conf. Ser. Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray . SPIE, Bellingham, p. 99051Q © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# The late inspiral of supermassive black hole binaries with circumbinary gas discs in the LISA band

, Volume 476 (2) – May 1, 2018
9 pages

/lp/ou_press/the-late-inspiral-of-supermassive-black-hole-binaries-with-ZxKTGgOf4f
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty423
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present the results of 2D, moving-mesh, viscous hydrodynamical simulations of an accretion disc around a merging supermassive black hole binary (SMBHB). The simulation is pseudo-Newtonian, with the BHs modelled as point masses with a Paczynski–Wiita potential, and includes viscous heating, shock heating, and radiative cooling. We follow the gravitational inspiral of an equal-mass binary with a component mass Mbh = 106 M⊙ from an initial separation of 60rg (where rg ≡ GMbh/c2 is the gravitational radius) to the merger. We find that a central, low-density cavity forms around the binary, as in previous work, but that the BHs capture gas from the circumbinary disc and accrete efficiently via their own minidiscs, well after their inspiral outpaces the viscous evolution of the disc. The system remains luminous, displaying strong periodicity at twice the binary orbital frequency throughout the entire inspiral process, all the way to the merger. In the soft X-ray band, the thermal emission is dominated by the inner edge of the circumbinary disc with especially clear periodicity in the early inspiral. By comparison, harder X-ray emission is dominated by the minidiscs, and the light curve is initially more noisy but develops a clear periodicity in the late inspiral stage. This variability pattern should help identify the electromagnetic counterparts of SMBHBs detected by the space-based gravitational-wave detector LISA. accretion, accretion discs, black hole physics, hydrodynamics 1 INTRODUCTION Supermassive black holes (SMBH) are currently believed to reside in most, if not all galactic nuclei, and SMBH binaries (SMBHBs) are expected to be produced frequently in mergers of galaxies (e.g. Kormendy & Ho 2013). The resulting compact SMBHBs are likely embedded in a gaseous environment. Electromagnetic (EM) signatures of such compact SMBHBs can arise from their interaction with gas. The space-based gravitational-wave (GW) detector LISA will be sensitive to SMBHBs with masses in the range 104–107 M⊙ (Amaro-Seoane et al. 2017). Identifying a GW source in the EM bands has considerable importance for astrophysics and cosmology (e.g. Phinney 2009). It has been proposed that periodic emission in the EM bands may track the orbital motion of SMBHBs throughout the late inspiral, and that this would then allow a unique identification of a LISA source (Kocsis, Haiman & Menou 2008). Comparing the EM and GW chirp signals would then also help probe the difference in the propagation speed of photons and gravitons (Kocsis et al. 2008; Haiman 2017), complementing measurements of the graviton speed from phasing of the GWs alone (see e.g. the review by Will 2006 and references therein). A similar constraint has recently been derived from the time delay between the GWs and the accompanying gamma-ray burst arriving from the neutron-star merger GW170817 (Abbott et al. 2017). An SMBHB offers an improved and more robust measurement of this time delay; this is especially the case if Doppler modulations of the EM chirp are detected from the orbital motion of the BHs themselves, fixing the relative phase of the EM and GW signals (Haiman 2017). In order to predict the EM signatures of SMBHB mergers in the LISA band, the first step is to resolve the gas dynamics during this late inspiral stage. The morphology of the gas, even on scales as large as that of the circumbinary disc, is sensitive to the behaviour of gas streams in the vicinity of the individual BHs (Tang, MacFadyen & Haiman 2017; Bowen et al. 2017a). One necessary requirement for reliably simulating the LISA stage is to resolve the accretion flows near the individual BHs, on scales down to that of the innermost stable circular orbit (ISCO; located at 6GMbh/c2 ≡ 6rg for a non-spinning BH). Furthermore, the gas dynamics cannot be in a steady state, since the binary orbit is evolving rapidly. Depending on the component masses, for a circular, comparable-mass binary, the LISA band corresponds to the inspiral from an initial separation of (50–200)GM/c2 all the way to the merger, and covers the last ≈100 to a few thousand orbits of the binary. The orbital evolution of the binary is likely overwhelmingly dominated by its gravitational radiation at this late stage (Haiman, Kocsis & Menou 2009; Dotti, Sesana & Decarli 2012; Tang et al. 2017). Because of the dual requirements of high spatial resolution and large number of orbits, simulating the LISA stage in its entirety has remained beyond the capability of fully three-dimensional (GR)MHD simulations. They have generally focused instead on the last ≈10 orbits preceding the merger (Bode et al. 2010; Farris et al. 2012; Giacomazzo et al. 2012; Noble et al. 2012; Shi et al. 2012; Bowen et al. 2017b). However, following the binary throughout the entire LISA band is feasible in 2D hydrodynamics. Previous high-resolution Newtonian simulations have indeed followed binaries for thousands of orbits. However, these works have focused on large separations, well before the LISA stage (Cuadra et al. 2009; Krolik 2010; Roedig et al. 2012; Farris et al. 2015a). As such, they did not directly resolve accretion inside the ISCO (typically either ignoring accretion, or implementing accretion via a sink prescription). In our earlier work (Farris et al. 2015b), we performed a 2D simulation of an inspiraling and merging BH binary, embedded in a viscous circumbinary disc. However, this earlier study was purely Newtonian, and employed an artificial sink that depleted the gas according to the fluid viscous time-scale. Moreover, the focus was on the overall accretion rate, and we did not present an analysis of the gas dynamics and did not compute the emerging EM emission and light curves. In this paper, we focus on modelling the entire inspiral stage observable with LISA, extending our earlier work. We consider an equal-mass, non-spinning BH binary with a component mass of Mbh = 106 M⊙ and initial separation 60rg. This corresponds to the last ≈100 orbits of the merger, or the last ≈(1 + z)/2 week observed for a binary at redshift z. The binary is assumed to be on a circular orbit, and we allow the orbit to shrink due to GW emission all the way to merger. The binary is embedded in a viscous, optical thick α-disc with radiative cooling. We employ a pseudo-Newtonian potential for the black holes. Our spatial resolution is sufficient to resolve the individual BH's ISCO, allowing us to remove gas inside the event horizon physically. We analyse the dynamics and the time-dependent spectral signatures of the emerging thermal emission. This paper is organized as follows. In Section 2, we summarize our numerical methods. In Section 3, we present our main findings, including disc spectra and light curves. We also examine the periodicities of the light curves separately in the soft and hard X-ray bands. In Section 4, we discuss the implications of our results, along with some caveats. Finally, in Section 5, we summarize our conclusions and indicate topics for future work. 2 NUMERICAL SET-UP The model and numerical set-up of this work closely resembles that in our previous work (Farris et al. 2015b). We briefly outline the key aspects here and highlight the modifications we made for this study. We refer the reader to the above paper for further details. All simulations are performed using the publicly available moving-mesh grid code DISCO in 2D (Duffell & MacFadyen 2012; Duffell 2016). We simulate an equal-mass SMBHB with Mbin = 2Mbh = 2 × 106 M⊙ and initial separation a0 = 60GMbh/c2. At this separation, a typical 106 M⊙ binary at z < 1 can be localized by LISA on the sky to a few square degrees, giving the chance to a wide-field telescope to locate and identify the source (Kocsis et al. 2008; Lang & Hughes 2008). We assume the BHs are non-spinning so that a0 = 10risco = 30rS, where rS is the Schwarzschild radius of an individual BH. The BHs are placed on the grid, and the vertically integrated fluid equations are solved assuming an α-viscosity with α = 0.1. Inside the cavity that forms around the binary (with a size roughly twice the initial binary separation), the grid spacing is roughly linear in radius and risco is resolved with ≈8 radial cells; there are ≈ 100 cells in total inside the ISCO. Outside the cavity, the radial grid spacing increases approximately logarithmically up to 30a0, allowing us to concentrate computational resources on the flow in the inner regions of the disc. The total number of radial cells is 512. At all radii, we chose our azimuthal grid spacing such that the aspect ratio of each cell is kept approximately unity, unless this requires a number of azimuthal grid cells Nϕ > 512, in which case we cap the number of ϕ-zones and set Nϕ = 512. To mimic the effect of black hole accretion, gas inside the event horizon (roughly the innermost 2 or 3 radial cells around each BH) is depleted effectively instantaneously at each time-step using a sink prescription similar to Tang et al. (2017). We use a sink radius rsink = rS and a sink time-scale ≈10−4 orbital period. In this paper, we assume the mass and orbit of the binary black hole are not affected by the accretion on the time-scale of our simulations. Our simulation is pseudo-Newtonian, i.e. we model the black holes as point masses with the modified potential (Paczyńsky & Wiita 1980)   $$\psi (r)=-GM/(r-r_{\rm S}).$$ (1)Since the Paczynski–Wiita potential diverges at r = rS, we employ an upper limit for the gravitational force: |f(r)| = |f(1.5rS)| for r < 1.5rS. We choose 1.5rS because our cell size is ≈0.5rS. We use a γ-law equation of state for the gas in the form P = (γ − 1)e, where P and e are the vertically integrated pressure and internal energy, respectively, and γ = 5/3 is the adiabatic index. Radiative cooling is incorporated naturally through the energy equation, as described in Farris et al. (2015a). We assume the disc is optically thick and geometrically thin (the optical depth due to electron scattering at r ∼ a near the time of decoupling is ∼105; Haiman et al. 2009). These assumptions are consistent with our simulation results. The cooling rate is $$q_{\rm cool}=\frac{4\sigma }{3\tau _{{\rm e}}} T^{4}$$, where T is the mid-plane temperature and τe is the optical depth for electron scattering (τe = Σσt/mp). For a Shakura-Sunyaev disc around a single BH, the cooling is in balance with the local viscous heating rate, qvis = 9/8αPΩ. Following Farris et al. (2015a) and Farris et al. (2015b), we parametrize our simulations by specifying the Mach number at the fixed initial radius r = a0, $$\mathcal {M}_{a}=v_{a}/\sqrt{P(a)/\Sigma (a)}$$, where va ≡ (GMbin/a)1/2 and Σ is the surface density, We then obtain the local cooling rate from the scaling   $$q_{\rm cool}=\frac{9}{8}\alpha \Omega (a)P(a)\mathcal {M}^{8}_{a}\left({\frac{P}{\Sigma v^{2}_a}}\right)^{4}\left(\frac{\Sigma }{\Sigma (a)}\right)^{-1},$$ (2)where Ω is the Keplerian angular frequency.1 Changing $$\mathcal {M}_{a}$$ is equivalent to changing the thickness h/r of the disc, and in this work we fix $$\mathcal {M}_{a}=(h/r)^{-1}=10$$. A standard steady-state quasar disc (Shakura & Sunyaev 1973) with a near-Eddington accretion rate $${\sim } \dot{M}_{\rm Edd}$$ has a scale height of h/r ∼ 10−3, and is much thinner than the disc in our simulations (h/r ∼ 0.1). Simulating an accretion disc with h/r ∼ 10−3 is prohibitively challenging, and our choice is typical of values employed in previous numerical works (e.g. D'Orazio, Haiman & MacFadyen 2013; Farris et al. 2014; Farris et al. 2015b; Ragusa, Lodato & Price 2016; Miranda, Muñoz & Lai 2017 have included slightly smaller h/r values ranging from 0.02 to 0.13) and allows for a direct comparison with those works. The optical thickness to electron scattering (τes > 10) and free–free absorption (τff > 104) are both very high in our simulation domain (see e.g. fig. 10 in Corrales, Haiman & MacFadyen 2010). Therefore, throughout this paper, we compute the spectrum assuming multicolour thermal body radiation   $$L_{\nu }=\int \frac{2h\nu ^{3}}{c^{2}(\exp \left[\frac{h\nu }{kT_{\rm eff}(r, \phi )}\right]-1)}{\rm d}A,$$ (3)where Teff(r, ϕ) is the effective temperature obtained from the local cooling rate qcool(r, ϕ) in the cell at (r, ϕ) via Teff = (qcool/σ)0.25 with σ the Stefan–Boltzmann constant. While computing the observed spectrum and light curve, we additionally consider the relativistic Doppler effect and the gravitational redshift. For simplicity, and to assess the maximal impact of the Doppler effect, we further assume the circumbinary disc is observed edge-on. Photons emitted by gas suffer a Doppler shift in frequency D = [Γ(1 − β∥)]−1, where Γ = (1 − β2)−1/2 is the Lorentz factor, and β∥ the component of the velocity along the line of sight. The apparent flux Fν at a fixed observed frequency ν is modified from the flux of a stationary source $$F_\nu ^{\rm 0}$$ to $$F_\nu =D^{3}F^{\rm 0}_{\nu /D}$$. For the gravitational redshift, we employ D = (1 − rS/Re)−0.5, where rS is the Schwarzschild radius and Re is the distance between the source of emission and the BH. The above is an approximation to a fully relativistic ray-tracing in the binary's dynamical metric (Schnittman et al. 2017). In particular, we ignore gravitational self-lensing (D'Orazio & Di Stefano 2017; Haiman 2017) and other relativistic effects, such as light-travel time modulations (i.e. Shapiro delay), which appear at order v/c. Our disc has an initial surface density profile Σ(r) = Σ0(r/a)−0.5 and contains an artificial cavity with the density depleted inside the radius 2.5a0, as in Tang et al. (2017). Viscosity in the circumbinary disc transports angular momentum outwards, causing the cavity to refill. We first hold the binary on the initial circular orbit and evolve the system for several hundred orbits to reach a quasi-steady state. We then allow the binary separation to shrink due to quadrupole formula for GW emission (Peters 1964),   $$a(t)=a_0(1-t/\tau )^{0.25}.$$ (4)Here, τ = 7.22 d is the GW inspiral time from the initial separation of 30 GMbin/c2. Finally, we employ a set of convenient code units, in which the initial binary period is T(0) = 2π and a0 = 1; as a result, τ = 385.17 in code units. 3 RESULTS As stated above, before allowing the black holes to inspiral and merge, we first hold the binary separation fixed, and evolve the system for several hundred orbits to reach a quasi-steady state. This quasi-steady state is similar to that in previous works (D'Orazio et al. 2013; Farris et al. 2014, 2015a; Tang et al. 2017) in that the cavity becomes lopsided, and two narrow streams periodically connect the minidiscs around the individual BHs with the circumbinary disc. Similar to Farris et al. (2015b), we find that these features persist until right before the merger, with only a gradual decline in the size and surface density of the minidiscs and the streams. In Fig. 1, we show the 2D distributions of surface density and effective temperature at three different snapshots, corresponding to t = 0 (the beginning of the inspiral), 0.5τ (half-way to merger), and 0.99τ (just before the merger). Figure 1. View largeDownload slide The top row shows snapshots of the surface density Σ/Σ0, where Σ0 is the initial surface density at r = a, and the bottom row shows snapshots of the effective temperature Teff in [eV]. The horizontal and vertical axes are in units of the initial binary separation a0. The snapshots are taken, from left to right, at the times t = 0.0τ, 0.5τ, and 0.99τ, where τ = 7.22 d corresponds to the merger. The cavity size decreases and the BHs accrete efficiently during the entire inspiral. Figure 1. View largeDownload slide The top row shows snapshots of the surface density Σ/Σ0, where Σ0 is the initial surface density at r = a, and the bottom row shows snapshots of the effective temperature Teff in [eV]. The horizontal and vertical axes are in units of the initial binary separation a0. The snapshots are taken, from left to right, at the times t = 0.0τ, 0.5τ, and 0.99τ, where τ = 7.22 d corresponds to the merger. The cavity size decreases and the BHs accrete efficiently during the entire inspiral. A spectrum computed at the beginning of the inspiral is also shown in Fig. 2. In this figure, we additionally show three distinct components of the spectrum, arising from outer disc (r > 3a0), the shocked streams inside the cavity and shocked hot gas near the cavity wall (a0 < r < 3a0) and minidiscs (r < a0). At frequencies below ∼1keV, the emission from the viscously heated outer disc dominates, while at higher frequencies, the emission mostly comes from the shock-heated gas in the minidiscs, streams, and near the cavity wall. At frequencies between ∼(1–20) keV, the total spectrum shows a plateau with a small depression near ∼3 keV, resembling the ‘notch’ discussed in Roedig, Krolik & Miller (2014). The spectrum computed in Farris et al. (2015a) has a larger magnitude at the highest frequencies, and does not show a plateau. This is because in our simulation we resolve the event horizon, and gas near the black hole is efficiently cleared out, as it falls rapidly into the BH. In Farris et al. (2015a), a much slower sink prescription was implemented, causing hot, shock-heated gas to accumulate around each black hole, and brightening the high-frequency emission. Figure 2. View largeDownload slide Thermal spectrum computed from a simulation snapshot at the beginning of the inspiral (t = 0). The full spectrum is shown by the blue curve. The distinct components arising from gas in the minidiscs (r < a), streams and cavity wall (a < r < 3a) and from the outer regions of the circumbinary disc (r > 3a) are represented by the orange dashed line, green dotted line, and red dash–dotted line, respectively. Figure 2. View largeDownload slide Thermal spectrum computed from a simulation snapshot at the beginning of the inspiral (t = 0). The full spectrum is shown by the blue curve. The distinct components arising from gas in the minidiscs (r < a), streams and cavity wall (a < r < 3a) and from the outer regions of the circumbinary disc (r > 3a) are represented by the orange dashed line, green dotted line, and red dash–dotted line, respectively. In Fig. 3, we show the total spectrum at three different snapshots, t = 0 (blue; same as in Fig. 2), 0.8τ (orange), and 0.99τ (green). The low-frequency part of the spectrum, which is dominated by the outer regions of the circumbinary disc, shows little change throughout the simulation. By comparison, the higher frequencies (≳ 1 keV), which arises from the minidiscs or streams in the cavity, declines noticeably by 0.8τ; the decline in the hard X-ray band becomes significant during the end of inspiral (0.99τ). Note that in this paper we chose Mbin = 2 × 106 M⊙, so that the GW signal is in the LISA band, and 0.99τ corresponds to the last ≈2 h of the inspiral. Figure 3. View largeDownload slide Total composite thermal spectra computed from simulation snapshots at t = 0 (blue curve), t = 0.8τ (orange dashed), and t = 0.99τ (green dot–dashed). Figure 3. View largeDownload slide Total composite thermal spectra computed from simulation snapshots at t = 0 (blue curve), t = 0.8τ (orange dashed), and t = 0.99τ (green dot–dashed). The scaling relationships in our simulation can be calculated by equating the viscous heating qvis( ∝ αΩP) to the radiative cooling qcool( ∝ Σ−1T4). Note that we fix a0 = 60GMbh/c2, so that the orbital velocity is a constant, independent of Mbin, and the angular velocity $$\Omega \propto M_{\rm bh}^{-1}$$. Since we have further fixed $$\mathcal {M}_{\rm a}=10$$, the mid-plane temperature T and ratio P/Σ are also fixed, i.e. P ∝ Σ. Finally, with α fixed as well, we have qvis ∝ Σ/Mbh and qcool ∝ Σ−1. This then yields the scaling $$\Sigma _{0}\propto M_{\rm bh}^{1/2}$$ and the effective temperature Teff ∝ Σ−1/4 scales as $$M_{\rm bh}^{-1/8}$$. Note that this implies that the photon frequency scales with total binary mass as $$\nu \propto T_{\rm eff}\propto M_{\rm bin}^{-1/8}$$, as well. For a steady disc around a solitary BH, it would further scale with Mach number as $$h\nu \propto \mathcal {M}_{\rm a}^{-5/4}$$ ($$\Sigma \propto \mathcal {M}^{-3}$$ and $$T\propto \mathcal {M}^{-2}$$). This suggests that for a thinner disc, the frequencies could be lower than shown for our fiducial $$\mathcal {M}_{a}=10$$ in Figs 2 and 3. However, for gas near the cavity wall and minidiscs, the variation of the disc thickness caused by its interaction with the binary BHs requires further study (e.g. Ragusa et al. 2016). In Fig. 4, we show h/r of our simulated disc. In the shock-heated region, h/r is increased to ∼0.2, while in the outer region h/r stays at ∼0.1. Figure 4. View largeDownload slide Local disc thickness h/r throughout the disc at t = 0. Figure 4. View largeDownload slide Local disc thickness h/r throughout the disc at t = 0. In Fig. 5, we illustrate the effects of the relativistic Doppler shift and gravitational redshift on the observed spectra. As mentioned in Section 2, we assume the disc is observed edge-on, so that the Doppler effect is maximized. Note the spectra in Fig. 5 do not include the inclination angle. The observed spectrum will be dimmer by a factor of cos (θ), where θ is the angle of inclination to the line of sight. We include four viewing angles in the plane of the disc, corresponding to the $$\pm \hat{x}$$ and $$\pm \hat{y}$$ directions. As the figure shows, at the low-frequency end (hν ≲ 1 keV), the spectrum is not modified by either the Doppler or gravitational shifts. At higher frequencies (1 keV ≲ hν ≲ 20 keV), the Doppler shift causes an overall dimming (i.e. red and blue curves), while at the highest frequencies (hν ≳ 20 keV), the spectrum is brightened significantly. Interestingly, we find that the dimming and the brightening retain their sign in all four different viewing directions; we offer a simple interpretation of this result in Section 4 below. For reference, the figure also shows a case with only the GR redshift (green curve), which is independent of viewing angle, and always reduces the observed flux. Figure 5. View largeDownload slide Spectra illustrating the impact of the special relativistic Doppler shift and the gravitational red/blue shift. For reference, the black curve shows the spectrum without considering either of these effects. The green curve shows the spectrum including only the gravitational redshift, while the other four spectra, shown in blue and red, include both relativistic effects (for the four different viewing directions $$\pm \hat{x}$$ and $$\pm \hat{y}$$, as labelled in the inset). All spectra are calculated from the snapshot at t = 0, when the orbital velocity of the BHs is 0.091c. Figure 5. View largeDownload slide Spectra illustrating the impact of the special relativistic Doppler shift and the gravitational red/blue shift. For reference, the black curve shows the spectrum without considering either of these effects. The green curve shows the spectrum including only the gravitational redshift, while the other four spectra, shown in blue and red, include both relativistic effects (for the four different viewing directions $$\pm \hat{x}$$ and $$\pm \hat{y}$$, as labelled in the inset). All spectra are calculated from the snapshot at t = 0, when the orbital velocity of the BHs is 0.091c. In order to study the variability of the spectrum at different frequencies, we calculate the light curves in the soft (at 2 keV) and hard (at 10 keV) X-ray bands. In Fig. 6, we show the light curve at the beginning of the simulation and right before merger, with and without the Doppler effect. In the 2 keV band, we find clear periodicity in the light curve throughout the inspiral. In the 10 keV band, the inspiral light curve at the beginning is much noisier. We can still see some periodicity by eye, but it is harder to distinguish than in the 2keV case. However, in the late stages of the inspiral, the 10 keV band also shows clear periodicity at the orbital period. The light curves also show an overall see-saw modulation, with a period approximately eight times longer than torb. Furthermore, the light curves show a clear a behaviour analogous to the GW chirp signal, in that the periodicity tracks the binary's increasing orbital frequency. Figure 6. View largeDownload slide Gravitational wave strain (top panels) and X-ray light curves (middle and bottom) calculated at the beginning (left) and end (right) of the simulated inspiral in the LISA band. The top panel shows the GW strain observed edge-on from a binary at redshift z = 1. The middle panel corresponds to the 2 keV band and the bottom panel to 10 keV band. The blue curves show the light curves neglecting the relativistic Doppler shift or gravitational redshift. The orange curves show the light curves including these effects. While the Doppler effect imposes a small additional periodic modulation due to the binary's orbital motion, its dominant effect is a nearly uniform dimming of the light curve (this results from the asymmetry between the large blue and redshifts due to the internal motion of gas in the minidiscs; see the text for discussion). The x-axis shows the look-back time in the source's rest frame, and the number of binary orbits (in parentheses). Figure 6. View largeDownload slide Gravitational wave strain (top panels) and X-ray light curves (middle and bottom) calculated at the beginning (left) and end (right) of the simulated inspiral in the LISA band. The top panel shows the GW strain observed edge-on from a binary at redshift z = 1. The middle panel corresponds to the 2 keV band and the bottom panel to 10 keV band. The blue curves show the light curves neglecting the relativistic Doppler shift or gravitational redshift. The orange curves show the light curves including these effects. While the Doppler effect imposes a small additional periodic modulation due to the binary's orbital motion, its dominant effect is a nearly uniform dimming of the light curve (this results from the asymmetry between the large blue and redshifts due to the internal motion of gas in the minidiscs; see the text for discussion). The x-axis shows the look-back time in the source's rest frame, and the number of binary orbits (in parentheses). In order to identify the origin of the periodicity in the light curves, in Fig. 7 we decompose them into emission arising from r < a (dominated by gas in/near the minidiscs) and from r > a (dominated by the circumbinary disc). This figure reveals that the emission from the circumbinary disc is clearly periodic from the beginning of the simulation all the way to the merger in both bands. This is because each BH periodically interacts with the cavity wall and creates a stream which hits the cavity wall, shock-heats the gas, and leads to enhanced X-ray emission. By comparison, the emission from the minidiscs is noisy in the beginning, and develops clear periodicity only closer to the merger. One explanation for this phenomenon is that in the earlier stage the minidiscs act as buffers of the accretion streams. Unlike steady-state α discs, the internal shocks in the minidiscs create significant noise; their overall shape often differs significantly from a circular disc. In the late stage of the inspiral (in the last 60 h before the merger), as the truncation radius decreases, this buffer effect declines, and the narrow accretion streams eventually connect the cavity wall and event horizon almost directly (see the right-hand panel of Fig. 1). Therefore, the noise produced by the minidiscs gradually disappears. Figure 7. View largeDownload slide Decomposition of the light curves. The orange and blue curves correspond to emission arising from r < a (dominated by gas in/near the minidiscs) and from r > a (dominated by the circumbinary disc). The green curves correspond to the total emission. The upper two panels show the light curves in the 2 keV band, and the bottom two panels in the 10 keV band. For clarity of display, the orange curves, which show sub-dominant emission from the minidiscs in the soft band in the upper two panels, use the y-axis labels on the right side. Figure 7. View largeDownload slide Decomposition of the light curves. The orange and blue curves correspond to emission arising from r < a (dominated by gas in/near the minidiscs) and from r > a (dominated by the circumbinary disc). The green curves correspond to the total emission. The upper two panels show the light curves in the 2 keV band, and the bottom two panels in the 10 keV band. For clarity of display, the orange curves, which show sub-dominant emission from the minidiscs in the soft band in the upper two panels, use the y-axis labels on the right side. In the 2keV band, emission from the circumbinary disc dominates, and the light curve therefore also shows very clear periodicity, during the entire inspiral process. The opposite situation occurs in the 10keV band, in which the total emission is dominated by the minidiscs. The light curve is noisy in the beginning, and shows much clearer periodicity in the end stage. From the bottom-most panel of Fig. 7, we see that the emission from the minidiscs and circumbinary discs has the same period, but with different phases. Therefore, in the later stages, when the strength of the emission from these two regions is comparable, the total light curve sometimes shows two maximas in a single period. Haiman (2017) and D'Orazio & Haiman (2017) proposed a toy model for the minidiscs as two ‘light bulbs’, both of which are steady in their own rest frame. Therefore, unless observed perfectly face-on, the Doppler modulation will lead to periodicity at the orbital frequency. Note that the light curve would still be modulated, even in the case of an equal-mass binary, because the red and blueshifts of the thermal emission from the two BHs, with equal but opposite line-of-sight velocities, will not exactly cancel. However, in our simulation, we find that the internal structure of the minidiscs is far from stationary, with shocks forming and streams of gas entering and leaving the fiducial Hill sphere. The emission from the minidisc regions is therefore quite noisy. As a result, the Doppler modulation may be difficult to extract from the observed light curve. We expect that this may still be feasible, given a GW template from LISA, with which the X-ray light curve can be cross-correlated, but this requires further analysis. Fortunately, we find that the shock-heated cavity wall provides a clear periodicity throughout the entire inspiral process. Moreover, since this periodicity does not originate from the Doppler effect (see Fig. 6), it should exist independently of the inclination angle. In order to assess how well the periodicity in the light  curves track the binary's orbital period, we rescale the light-curve time series with the evolving binary orbital period. In Fig. 8, we plot the Lomb–Scargle periodograms with and without this time rescaling. The rescaled periodograms show peaks at ≈2 × Ωorb, both at the beginning and the end of the inspiral, implying that this periodicity is able to follow the accelerating binary orbital frequency (i.e. the ‘chirp’) throughout the entire inspiral process. The first peak in the periodogram without rescaling arises from a dense lump in the lopsided cavity wall, as reported in a series of previous works (e.g. MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Krolik 2010; Roedig et al. 2011; Noble et al. 2012; Shi et al. 2012; D'Orazio et al. 2013; Farris et al. 2014; Bowen et al. 2017a). This feature is smeared out and disappears after time rescaling, because it does not follow the binary's frequency chirp. More specifically, the frequency of this low-frequency peak increases from ≈0.25Ωorb(0) to ≈0.4Ωorb(0), but this increase is slower than the evolving binary orbital frequency Ωorb(t). Just prior to the merger, the period of this low-frequency modulation is around 15 orbits. Similarly, in Fig. 1, we see that the overall size of the cavity is shrinking, but cannot keep up with the binary separation, which decreases more rapidly. Nevertheless, this low-frequency modulation persists throughout the inspiral process, as we can clearly see in Fig. 6. In principle, the evolution of this low-frequency modulation could be measured in the X-ray light curves, and would probe the gas morphology near the central cavity. Figure 8. View largeDownload slide Lomb–Scargle periodograms of the light curves. The four panels correspond to those in Fig. 6, namely to the 2 and 10 keV bands (top and bottom, respectively) and to the beginning and end of the simulated inspiral (left and right, respectively). In each panel, the red curve corresponds to the periodogram of the original light curves shown in Fig. 6. The blue curves correspond to scaled light curves, in which time is measured in units of the instantaneous binary orbital period. Figure 8. View largeDownload slide Lomb–Scargle periodograms of the light curves. The four panels correspond to those in Fig. 6, namely to the 2 and 10 keV bands (top and bottom, respectively) and to the beginning and end of the simulated inspiral (left and right, respectively). In each panel, the red curve corresponds to the periodogram of the original light curves shown in Fig. 6. The blue curves correspond to scaled light curves, in which time is measured in units of the instantaneous binary orbital period. 4 DISCUSSION One of the motivations of this work was to compute the modulation of the light curves due to relativistic Doppler effect. Naively, such a modulation inevitably arises from the same orbital motion that produces GWs, and therefore the GW and X-ray chirp signals should track one another with a known phase (Haiman 2017). This would be of particular interest for LISA, allowing a robust measurement of the speed difference of photons and gravitons, and more generally, a secure identification of the EM counterpart of a LISA binary. Our results suggest that it may be difficult to extract this Doppler modulation from the X-ray light curve. We find efficient fuelling of both BHs through the inspiral process, with gas funnelled inside each BH's Hill radius, forming tidally truncated ‘minidiscs’. These minidiscs persist nearly all the way to the merger, until the tidal truncation radius shrinks to a size comparable to that of the ISCOs (roughly 60 h prior to merger). However, the effective temperature, density, and velocity fields of the gas inside these minidiscs are not stationary, and the net Doppler shift is dominated by these large ‘internal’ gas motions, rather than the orbital velocity of the BHs. Interestingly, as Figs 5 and 6 show, the most prominent effect of the Doppler shift in the 2–10 keV range is an overall dimming of the light curve, rather than a sinusoidal modulation. This is because the hot gas dominating this emission is located very close to each BH (near or inside its Hill radius), and the velocity of this gas is larger than the orbital speed of the BHs. Regardless of the viewing direction, there is always gas with both large blue and redshifts along the line of sight, and the total composite spectrum includes both of these shifts. The net result can be either an overall dimming or brightening, depending on the local slope and curvature of the spectrum. In other words, for the steeply curved blackbody spectral shape, the Doppler effect is not symmetric under switching the sign of the line-of-sight velocity. Fig. 5 shows that in the 2–10 keV range, the redshift and the corresponding dimming from receding gas dominate over the blueshift and brightening from approaching gas. We emphasize that this result is sensitive to the spectral shape. One caveat is that, as mentioned above, due to numerical constraints our simulated disc is thicker than a typical active galactic nucleus (AGN) disc of a solitary BH. If we were to scale the simulation to correspond to a higher Mach number, this would shift the normalization of the photon energy. For example, the soft versus hard X-ray bands could become FUV and EUV bands for a disc with a Mach number of $$\mathcal {M}_{a} \sim 1000$$, instead of the $$\mathcal {M}_{a}=10$$ we adopted (recall that $$h\nu \propto \mathcal {M}_{\rm a}^{-5/4}$$ for single black hole accretion disc). On the other hand, the gas in the inner regions is being heated primarily by shocks, rather than viscosity, and we find that h/r in the inner regions of the flow remains ∼ 0.2 throughout the simulation, justifying the scaling to lower Mach numbers (see Fig. 4). However, it requires further study to understand how this depends on the choice of $$\mathcal {M}_{a}$$. A higher $$\mathcal {M}_{a}$$ could lead to a larger cavity and reduce the gas fuelling from the cavity wall (Ragusa et al. 2016). This may reduce the luminosity of the minidiscs. We also note that in equation (2), we assumed that the fluid is gas-pressure dominated; this assumption is likely violated in the shock-heated streams and minidiscs. Another caveat in our conclusions is that we focused on thermal emission. We have verified that the effective absorptive opacity, defined as τeff = [τabs(τabs + τscat)]1/2, where τabs is the optical depth to free–free absorption and τscat is the optical depth to electron scattering, remains τeff ∼ 104 ≫ 1 throughout our simulation. This justifies considering a thermalized spectrum. However, empirically, the X-ray emission from an AGN consists of a combination of thermal disc emission and a hot corona, and in this paper, we have not included the latter component. If the emission from a corona turns out to be more stable (i.e. if the corona lacks large internal bulk motions), the Doppler modulation in the X-ray band may be dominated by the BH's orbital motion, and much more prominent. Moreover, in this paper, we have focused on an equal-mass binary. We expect that the Doppler modulation would be more prominent in the case of low-mass ratio binaries, because the circum-primary and -secondary minidiscs would be more stable, and there would be less cancellation between their emission. The overall shape of the light curves we found show a clear ‘X-ray chirp’, caused by the interaction between the binary and the gas discs. The periodicity of the chirp tracks the shrinking binary orbital period torb, in tandem with the GW signal. Since the chirp is produced by hydrodynamical modulation, its absolute phase is not tied directly to the GW emission. Nevertheless, the speed difference between photons and gravitons should still be measurable, in principle, given a precise measurement of the phase-evolution of both the GW and the X-ray chirp signals. In the future, it would be interesting to compute the accuracy to which the phase difference can be measured in practice, if the absolute phase is unknown ab-initio, given realistic observational S/N from LISA, and from an X-ray instrument such as Athena,2 Lynx,3 ULTRASAT,4 or EXTP (Zhang et al. 2016). In this paper, we neglect spin of the BHs, but the spin would change the GW signal significantly, and would further complicate a direct comparison of the X-ray light curve with the LISA signal. The ISCO for a spinning black hole is smaller than a non-spinning one, which may enhance the luminosity at high frequencies. We have made several simplifying assumptions throughout this paper, which we intend to relax in our future work. We have focused here on an equal-mass binary, but we intend to study the dependence of the spectra and light curves on the binary mass ratio. Similarly, we intend to determine the dependence of these features on the disc thickness. We also used an α-prescription for the viscosity. We intend to perform magnetohydrodynamic (MHD) simulations in the future (e.g. Krolik 2010). In this work, we have also assumed that the binary BHs are on circular orbits, with zero eccentricity. On the other hand, eccentricity is likely to develop as a result of strong gas torques during the earlier stages of the inspiral (e.g. Cuadra et al. 2009; Roedig et al. 2011; Sesana, Gualandris & Dotti 2011; Roedig et al. 2012; Miranda, Muñoz & Lai 2017). In the future, we intend to study live binaries, and follow the development of eccentricity, and to assess the impact of eccentricity on the accretion rates and the predicted spectra and light curves. In the present simulation, we followed a pseudo-Newtonian approach; a full relativistic simulation of the gas in the vicinity of the event horizon could modify our results. Finally, we started our current simulations somewhat past the nominal ‘decoupling’ stage. In the cartoon picture of the binary inspiral (e.g. Milosavljevic & Phinney 2005), once the binary is sufficiently compact that the GW inspiral time-scale is shorter than the viscous time-scale in the nearby disc, the BHs outpace the disc and leave the disc gas behind. The decoupling separation adec can be calculated with tgw = tvis, where $$t_{\rm gw}=\frac{5c^{5}a^{4}}{64G^{3}m_{1}m_{2}(m_{1}+m_{2})}$$ is the GW inspiral time-scale (Peters 1964) and $$t_{\rm vis}=\frac{2(2a)^2}{3\nu }$$ is the viscous time in an unperturbed disc at the radius r = 2a. With the parameters (α = 0.1; $$\mathcal {M}_{a}=10$$) adopted in our simulation, we find the decoupling separation adec = 128GMbh/c2, which is somewhat larger than our initial separation of 60GMbh/c2. In this picture, our procedure of holding the binary's orbit initially fixed, and allowing it to reach a steady state would not be justified. However, we first note that even in a simple 1D toy model, the decoupling radius can only be taken as an order-of-magnitude estimate. For example, torques modify the density profile outside the binary's orbit, causing a pile-up, increasing the pressure, steepening the radial gradients, shortening the viscous time, and decreasing adec (e.g. Milosavljevic & Phinney 2005; Kocsis, Haiman & Loeb 2012). More importantly, from the results we find here, as well as from our earlier work (Farris et al. 2015b), we see that the cartoon picture is inaccurate, and the cavity can, in fact, shrink and follow the inspiraling binary, even well past the nominal decoupling. This is because the gas dynamics near the black holes is dominated by gravity and shocks, and not by viscosity, and the gas configuration is neither steady nor axisymmetric. More specifically, in Farris et al. (2015b), a simulation was started at the separation corresponding to tgw = 5tvis, i.e. safely before the nominal decoupling. Nevertheless, we found that the gas flows inwards, following the binary, and the accretion rate is only halved, even well past the decoupling, 0.1tvis before the merger. Similarly, in our simulation here, the spectrum is only slightly modified even at 0.8τ (see Fig. 3). We therefore expect our way of setting up the initial condition to be reasonably accurate. However, one caveat to note here is that the nominal decoupling separation is larger for a disc with smaller aspect ratios h/r (i.e. large Mach numbers). This means that the reduction of accretion rate and luminosity in more realistic, thinner discs may be more significant than the discs in our simulations. We intend to study this in future work, by beginning simulation near or prior to the decoupling stage and with larger Mach numbers. 5 CONCLUSIONS We have performed 2D simulations of accretion on to an inspiraling and merging binary BH system, and studied the corresponding EM signatures. We have computed the evolving multicolour blackbody emission from snapshots of our simulations, including the effects of the relativistic Doppler shift and gravitational redshift. The most important conclusions of this work can be summarized as follows: The light curves show clear periodicity. Accretion on to the individual BHs has long been known to be periodic (Artymowicz & Lubow 1994), and it is not surprising that any luminosity produced by gas near the individual BHs is periodic. However, we have found that the emission arising from gas farther out, near the cavity wall, which dominates the softer X-ray bands, shows a similar (and even clearer) periodicity. This is because streams of gas, flung outwards by the BHs, periodically hit and shock the cavity wall. The frequency is twice the binary's instantaneous orbital frequency. Distinct behaviour in different bands. In the beginning of the simulation, the emission from the minidiscs is noisy, but develops clear periodicity during the later stages. In the soft X-ray band, the emission from the circumbinary disc dominates, and the light curve shows a clear periodicity from the beginning of the simulation (a = 60GM/c2). By comparison, in the hard X-ray band, the emission in the beginning is dominated by the minidiscs, and periodicity in the early-inspiral light curve is obscured by noise, with clearer periodicity developing only in the late stages. Doppler modulation is sub-dominant. We find that the Doppler modulation of the light curve does not follow a simple toy model of two moving light bulbs in orbit. This is because the effective temperature, density, and velocity field of the gas in the minidiscs are not stationary, and the line-of-sight velocity of this gas is dominated by the large internal gas motions inside the minidiscs, rather than the clean orbital motion of the BHs. We find that the most conspicuous effect of the Doppler shift is an overall dimming (1keV≲ hν ≲ 20keV) or brightening (hν ≳ 20keV) of the light curves. We caution that these results are sensitive to the spectral shape and emission mechanism. In particular the Doppler effect could more closely track the binary's orbital motion for X-ray emission from more stable coronae around the individual BHs which do not have large internal motions. The EM chirp follows the GW inspiral. We have found that the minidiscs persist until the very end of the inspiral process. The discs are truncated at the gradually shrinking size of the Hill radius. The thermal emission from the minidiscs declines gradually, but remains overall significant all the way to the merger. The overall periodicities in both the soft and hard X-ray bands track the shrinking orbital period of the binary. No ‘second decoupling’. Fontecilla, Chen & Cuadra (2017) suggested that the inner (circum-primary) disc will become geometrically thick (h/r ∼ 1) in the late stages of the inspiral, as a result of tidal heating. They find, using 1D analytical models and numerical simulation, that the viscous time in the inner disc decreases to tvis ∼ tGW, and it effectively decouples from the binary. By comparison, in our 2D simulation, the size and surface density of the two minidiscs are continuously adjusted, and eventually, in the late stages, the minidiscs disappear, and narrow accretion streams instead fall inside the event horizon directly. The inspiral process does thicken the minidiscs, but the aspect ratios remain below h/r ∼ 0.2. The simulations presented here are based on many simplifications. Nevertheless, our results suggest that an EM chirp signal may accompany the GW emission, with the EM period tracking that of the GW signal right up to the merger of an SMBH binary. A detection of this EM chirp could provide a potential ‘smoking gun’ evidence to uniquely identify the EM counterpart of an SMBH binary detected by LISA. Acknowledgements We thank Daniel D'Orazio, Paul Duffell, Julian Krolik, and Geoffrey Ryan for useful discussions. Financial support was provided by NASA through ATP grant NNX15AB19G and ADAP grant NNX17AL82G and by NSF grants 1715661 and 1715356. ZH also gratefully acknowledges support from a Simons Fellowship in Theoretical Physics and hospitality by NYU during his sabbatical leave when this work began. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. Footnotes 1 We note there is a typographical error in equation (2) of Farris et al. (2015a); the correct exponent of the Mach number is $$\mathcal {M}_a^{8}$$, as in equation (2) above. 2 See http://www.cosmos.esa.int/web/athena. 3 See http://www.astro.msfc.nasa.gov/lynx. 4 See http://www.weizmann.ac.il/ultrasat. REFERENCES Abbott B. P. et al.  , 2017, ApJ , 848, L13 https://doi.org/10.3847/2041-8213/aa920c CrossRef Search ADS   Amaro-Seoane P. et al.  , 2017, preprint (arXiv:1702.00786) Artymowicz P., Lubow S. H., 1994, ApJ , 421, 651 https://doi.org/10.1086/173679 CrossRef Search ADS   Bode T., Haas R., Bogdanović T., Laguna P., Shoemaker D., 2010, ApJ , 715, 1117 https://doi.org/10.1088/0004-637X/715/2/1117 CrossRef Search ADS   Bowen D. B., Mewes V., Campanelli M., Noble S. C., Krolik J. H., Zilhao M., 2017a, ApJ , 853, L17 https://doi.org/10.3847/2041-8213/aaa756 CrossRef Search ADS   Bowen D. B., Campanelli M., Krolik J. H., Mewes V., Noble S. C., 2017b, ApJ , 838, 42 https://doi.org/10.3847/1538-4357/aa63f3 CrossRef Search ADS   Corrales L. R., Haiman Z., MacFadyen A., 2010, MNRAS , 404, 947 https://doi.org/10.1111/j.1365-2966.2010.16324.x CrossRef Search ADS   Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, MNRAS , 393, 1423 https://doi.org/10.1111/j.1365-2966.2008.14147.x CrossRef Search ADS   D'Orazio D. J., Di Stefano R., 2017, MNRAS , 474, 2975 https://doi.org/10.1093/mnras/stx2936 CrossRef Search ADS   D'Orazio D. J., Haiman Z., 2017, MNRAS , 470, 1198 https://doi.org/10.1093/mnras/stx1269 CrossRef Search ADS   D'Orazio D. J., Haiman Z., MacFadyen A., 2013, MNRAS , 436, 2997 https://doi.org/10.1093/mnras/stt1787 CrossRef Search ADS   Dotti M., Sesana A., Decarli R., 2012, Adv. Astron. , 2012 Duffell P. C., 2016, ApJS , 226, 2 https://doi.org/10.3847/0067-0049/226/1/2 CrossRef Search ADS   Duffell P. C., MacFadyen A. I., 2012, ApJ , 755, 7 https://doi.org/10.1088/0004-637X/755/1/7 CrossRef Search ADS   Farris B. D., Gold R., Paschalidis V., Etienne Z. B., Shapiro S. L., 2012, Phys. Rev. Lett. , 109, 221102 https://doi.org/10.1103/PhysRevLett.109.221102 CrossRef Search ADS PubMed  Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ , 783, 134 https://doi.org/10.1088/0004-637X/783/2/134 CrossRef Search ADS   Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015a, MNRAS , 446, L36 https://doi.org/10.1093/mnrasl/slu160 CrossRef Search ADS   Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015b, MNRAS , 447, L80 https://doi.org/10.1093/mnrasl/slu184 CrossRef Search ADS   Fontecilla C., Chen X., Cuadra J., 2017, MNRAS , 468, L50 Giacomazzo B., Baker J. G., Miller M. C., Reynolds C. S., van Meter J. R., 2012, ApJ , 752, L15 https://doi.org/10.1088/2041-8205/752/1/L15 CrossRef Search ADS   Haiman Z., 2017, Phys. Rev. D , 96, 023004 https://doi.org/10.1103/PhysRevD.96.023004 CrossRef Search ADS   Haiman Z., Kocsis B., Menou K., 2009, ApJ , 700, 1952 https://doi.org/10.1088/0004-637X/700/2/1952 CrossRef Search ADS   Kocsis B., Haiman Z., Menou K., 2008, ApJ , 684, 870 https://doi.org/10.1086/590230 CrossRef Search ADS   Kocsis B., Haiman Z., Loeb A., 2012, MNRAS , 427, 2680 https://doi.org/10.1111/j.1365-2966.2012.22118.x CrossRef Search ADS   Kormendy J., Ho L. C., 2013, ARA&A , 51, 511 https://doi.org/10.1146/annurev-astro-082708-101811 CrossRef Search ADS   Krolik J. H., 2010, ApJ , 709, 774 https://doi.org/10.1088/0004-637X/709/2/774 CrossRef Search ADS   Lang R. N., Hughes S. A., 2008, ApJ , 677, 1184 https://doi.org/10.1086/528953 CrossRef Search ADS   MacFadyen A. I., Milosavljević M., 2008, ApJ , 672, 83 https://doi.org/10.1086/523869 CrossRef Search ADS   Milosavljevic M., Phinney E. S., 2005, ApJ , 622, L93 https://doi.org/10.1086/429618 CrossRef Search ADS   Miranda R., Muñoz D. J., Lai D., 2017, MNRAS , 466, 1170 https://doi.org/10.1093/mnras/stw3189 CrossRef Search ADS   Noble S. C., Mundim B. C., Nakano H., Krolik J. H., Campanelli M., Zlochower Y., Yunes N., 2012, ApJ , 755, 51 https://doi.org/10.1088/0004-637X/755/1/51 CrossRef Search ADS   Paczyńsky B., Wiita P. J., 1980, A&A , 88, 23 Peters P. C., 1964, Phys. Rev. , 136, B1224 https://doi.org/10.1103/PhysRev.136.B1224 CrossRef Search ADS   Phinney E. S., 2009, Astro2010: The Astronomy and Astrophysics Decadal Survey, Science White Papers , National Research Council (2010) Washington, 235 Ragusa E., Lodato G., Price D. J., 2016, MNRAS , 460, 1243 https://doi.org/10.1093/mnras/stw1081 CrossRef Search ADS   Roedig C., Dotti M., Sesana A., Cuadra J., Colpi M., 2011, MNRAS , 415, 3033 https://doi.org/10.1111/j.1365-2966.2011.18927.x CrossRef Search ADS   Roedig C., Sesana A., Dotti M., Cuadra J., Amaro-Seoane P., Haardt F., 2012, A&A , 545, A127 https://doi.org/10.1051/0004-6361/201219986 CrossRef Search ADS   Roedig C., Krolik J. H., Miller M. C., 2014, ApJ , 785, 115 https://doi.org/10.1088/0004-637X/785/2/115 CrossRef Search ADS   Schnittman J. D., Dal Canton T., Camp J., Tsang D., Kelly B. J., 2017, ApJ , 853, 123 https://doi.org/10.3847/1538-4357/aaa08b CrossRef Search ADS   Sesana A., Gualandris A., Dotti M., 2011, MNRAS , 415, L35 https://doi.org/10.1111/j.1745-3933.2011.01073.x CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Shi J.-M., Krolik J. H., Lubow S. H., Hawley J. F., 2012, ApJ , 749, 118 https://doi.org/10.1088/0004-637X/749/2/118 CrossRef Search ADS   Tang Y., MacFadyen A., Haiman Z., 2017, MNRAS , 469, 4258 https://doi.org/10.1093/mnras/stx1130 CrossRef Search ADS   Will C. M., 2006, Living Rev. Relativ. , 9, 3 https://doi.org/10.12942/lrr-2006-3 CrossRef Search ADS PubMed  Zhang S. N. et al.  , 2016, in Jan-Willem A., den H., Tadayuki T., Marshall B., eds, Proc. SPIE Conf. Ser. Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray . SPIE, Bellingham, p. 99051Q © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations