On the orbital evolution of supermassive black hole binaries with circumbinary accretion discs

TL;DR

This study uses hydrodynamical simulations to analyze how circumbinary accretion discs influence the orbital evolution of supermassive black hole binaries, highlighting the importance of gas torques in merger timescales.

Contribution

It provides detailed insights into the gravitational torques from circumbinary discs on SMBH binaries and their dependence on accretion modeling parameters.

Findings

Gas torques are mainly gravitational and from gas near the black holes.

Slower sink rates lead to more negative torques and faster mergers.

Binary inspiral timescale is about 3 million years, affecting GW signals.

Abstract

Gaseous circumbinary accretion discs provide a promising mechanism to facilitate the mergers of supermassive black holes (SMBHs) in galactic nuclei. We measure the torques exerted on accreting SMBH binaries, using 2D, isothermal, moving-mesh, viscous hydrodynamical simulations of circumbinary accretion discs. Our computational domain includes the entire inner region of the circumbinary disk with the individual black holes (BHs) included as point masses on the grid and a sink prescription to model accretion onto each BH. The BHs each acquire their own well-resolved accretion discs ("minidiscs"). We explore a range of mass removal rates for the sink prescription removing gas from the central regions of the minidiscs. We find that the torque exerted on the binary is primarily gravitational, and dominated by the gas orbiting close behind and ahead of the individual BHs. The torques from the distorted circumbinary disc farther out and from the direct accretion of angular momentum are subdominant. The torques are sensitive to the sink prescription: slower sinks result in more gas accumulating near the BHs and more negative torques, driving the binary to merger more rapidly. For faster sinks, the torques are less negative and eventually turn positive (for unphysically fast sinks). When the minidiscs are modeled as standard alpha discs, our results are insensitive to the choice of sink radius. Scaling the simulations to a binary orbital period tbin = 1yr and background disc accretion rate Mdot = 0.3MEdd in Eddington units, the binary inspirals on a timescale of 3X10^6 years, irrespective of the SMBH masses. For binaries with total mass <10^7Msun, this is shorter than the inspiral time due to gravitational wave (GW) emission alone, implying that gas discs will have a significant impact on the SMBH binary population and can affect the GW signal for Pulsar Timing Arrays.

Figures

Tables

Table 1. Table 1: Summary of our simulation runs. The columns, from left to right, show the (i) sink time-scale, (ii) the corresponding effective viscosity parameter, (iii) sink radius, and (iv) total dimensionless specific torque.

${\tau }_{\mathrm{sink}}$	effective $\alpha$	$r_{\mathrm{sink}}$	$T_{\mathrm{tot}}/\dot{M}$
0.0125	238.51	0.1	0.433
0.5	5.96	0.1	0.373
1.25	2.38	0.1	0.269
2.5	1.19	0.1	-0.0271
5.0	0.59	0.1	-0.625
29.84	0.1	0.1	-5.77
9.621	0.1	0.05	-5.52

Equations

$$\Sigma (r)={\Sigma }_0\exp [-(r/r_0{)}^{-10}]r^{-0.5}.$$

$$c_s=(h/r)\left(\frac{M_1}{r_1}+\frac{M_2}{r_2}\right).$$

$$t_{\mathrm{vis}}=\frac{2}{3}{\left[\alpha {\left(\frac{h}{r}\right)}^2\Omega (r)\right]}^{-1}\sim 300{\left(\frac{r}{2a}\right)}^{3/2}{t}_{\mathrm{bin}},$$

$${\left(\frac{d\Sigma }{dt}\right)}_{\mathrm{sink}}\equiv -\frac{\Sigma }{{\tau }_{\mathrm{sink}}}.$$

$$t_{\mathrm{vis}}=\frac{2}{3}\frac{r^2}{\nu }=\frac{1}{3\pi }\frac{1}{\alpha (h/r{)}^2}{\left(\frac{r}{a}\right)}^{3/2}{u}^{-1/2}{t}_{\mathrm{bin}}.$$

$$T_{\mathrm{tot}}=T_{\mathrm{gr}}+T_{\mathrm{acc}}.$$

$$T_{\mathrm{gr}}=\sum _{\mathrm{cells}}{r}_{\mathrm{bh}}\times f_{\mathrm{gr}}.$$

$$T_{\mathrm{acc}}=\sum _{\mathrm{sink}}{r}_{\mathrm{bh}}\times \dot{m_i}(v_i-v_{\mathrm{bh}}),$$

$$M_{\mathrm{bin}}\times \frac{d\ell }{dt}=T_{\mathrm{tot}},$$

$$T_{\mathrm{tot}}=(-0.209{\tau }_{\mathrm{sink}}+0.437)\dot{M}.$$

$$\Sigma (r)=A/\sqrt{r}-B/r,$$

$$B/A=\dot{J}/\dot{M}.$$

$$\frac{⟨\dot{J}⟩}{⟨\dot{M}⟩}=1.02\left(\frac{B}{A}\right)+0.21,$$

Full text

††pubyear: 2017††pagerange: On the orbital evolution of supermassive black hole binaries with circumbinary accretion discs–References

1 Introduction

Binaries embedded in gas discs occur frequently in nature. Examples include planet+star systems with protoplanetary discs (e.g. Ward, 1997), young stellar binaries (e.g. Orosz et al., 2012), and supermassive black hole binaries (SMBHBs) expected in many galactic nuclei (e.g. Dotti et al., 2012; Mayer, 2013). The nature of the binary-disc interaction is of compelling interest for understanding the long-term evolution of these systems. In the case of SMBHBs, in particular, a long-standing question is whether such binaries can reach orbital separations as small as $\sim 1$pc, so that gravitational radiation would cause an inevitable merger (e.g. Begelman et al., 1980). Gas discs could help facilitate orbital decay down to these separations (e.g. Escala et al., 2005). Understanding the interaction between SMBHBs and their gas discs is also crucial for predicting the electromagnetic (EM) signatures of such systems (e.g. Tanaka et al., 2012).

Significant work has been done to understand the linear binary-disc interaction in the limit of small mass ratio, $q\equiv M_{2}/M_{1}\ll 1$. This so-called Type I planetary migration is facilitated by linear spiral density waves launched at resonant disc locations (e.g. Goldreich & Tremaine, 1980). The rate and even the direction (inward or outward) of Type I migration is sensitive to thermodynamics (e.g. Paardekooper & Mellema, 2006). The regime 10^{-4}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}q\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.05 has also been widely explored. The secondary opens a narrow annular gap in the disc, resulting in so-called Type II migration (e.g. Ward, 1997). The interaction is non-linear, generally causing inward migration on a time-scale comparable to the viscous time-scale of the disc near the binary (e.g. Lin & Papaloizou, 1986). Deviations from this time-scale can, however, can be significant, especially when the mass of the smaller BH exceeds that of the nearby disc (e.g. Syer & Clarke, 1995; Ivanov et al., 1999; Duffell et al., 2014; Dürmann & Kley, 2017).

By comparison, relatively little analogous work has been done for SMBH binaries in the range 0.1\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}q\leq 1. Studies of the demography and assembly of SMBHs in hierarchical cosmologies, based on the merger histories of dark matter haloes, find broad distributions between 10^{-2}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}q<1, generally peaking in the range $q\sim 0.1-1$ (e.g. Volonteri et al., 2003; Sesana et al., 2005; Lippai et al., 2009; Sesana et al., 2012). MacFadyen & Milosavljević (2008), D’Orazio et al. (2013) and Miranda et al. (2017) have made use of 2D grid-based hydrodynamical simulations to study SMBHBs embedded in thin $\alpha$ discs (Shakura & Sunyaev, 1973), and included a discussion of the disc torques. However, these simulations excised the innermost region surrounding the binary, and imposed an inner boundary condition, neglecting potentially important gas dynamics and torque contributions from inside the excised region. Farris et al. (2014) included the innermost region in their simulated domain, but their study did not focus on the binary-disc interaction, and did not present measurements of the gas torques from that region. Closest to the present paper are the previous studies by Cuadra et al. (2009) and Roedig et al. (2012), both of which followed the interaction between SMBHBs with $q=1/3$ and a self-gravitating circumbinary disc, using 3D SPH simulation. We will present a comparison to these studies below.

In the present work, we use a 2D grid-based code to evolve a thin (aspect ratio $h/r=0.1$) $\alpha$ disc for a few viscous times, until the system has relaxed to a quasi-steady state. We use DISCO, a high-resolution, finite-volume moving-mesh code (Duffell & MacFadyen, 2011, 2012, 2013), which allows us to not impose an inner boundary condition, and therefore to include the BHs in the computational domain. Our main goal is to accurately measure the disc-binary torques, including contributions from gravitational torques, as well as effective torques from the accretion of mass and angular momentum by the BHs. Accretion is implemented with a sink prescription; we perform numerical experiments to study how the mass removal rate in the sinks affects the global structure of the circumbinary disc and the torques exerted on the binary.

This paper is organized as follows. In § 2, we summarise our numerical methods. In § 3, we present our main results on the disc torques. We also determine the dependence of the disc morphology and the disc torque on the sink time scale $\tau_{\rm sink}$. In § 4, we discuss the implications of our results, list some caveats, and outline topics for future work. Finally, in § 5, we summarise our conclusions. Throughout this paper, we adopt units in which Newton’s constant, the binary mass, and the binary separation are $G=M_{\rm bin}=2\times M_{\rm bh}=a=1$, and the binary’s orbital period is $t_{\rm bin}=2\pi$.

2 Numerical Setup

2.1 Summary of Simulation Setup

The model and numerical setup of this work closely resembles that in our previous work (Farris et al., 2014, 2015). We only briefly outline the key aspects and some modifications here, and refer the reader to the above papers for further details. All simulations are performed using the DISCO code in 2D; a recent version of this code, including 3D magnetohydrodynamics functionality, is publicly available (Duffell, 2016a, b).

Our initial disc configuration is similar to that in Farris et al. (2015). We modify the standard steady-state Shakura–Sunyaev disc profile by exponentially reducing the density and pressure at small radii by a factor of $\exp[-(r/r_{0})^{-10}]$ to create a hollow cavity within $r_{0}\equiv 2.5a$,

[TABLE]

We do not include self-gravity; this assumption is justified for compact binaries in the inner regions of the disc. As a result, the normalization of the surface density $\Sigma_{0}$ does not affect the simulation and is arbitrary. For convenience, we set $\Sigma_{0}=1$ in code units.

We simulate the hydrodynamical evolution of a thin circumbinary disc by solving the vertically averaged, 2D Navier–Stokes equations. The pressure is set to enforce a locally isothermal equation of state of the form

[TABLE]

Throughout this work, we choose $h/r=constant=0.1$. We employ a standard $\alpha$-law viscosity (Shakura & Sunyaev, 1973) of the form $v=\alpha c_{s}^{2}\Omega_{k}^{-1}$, with $\alpha=constant=0.1$ (see Appendix A in D’Orazio et al. 2013 and Farris et al. 2014 for detailed descriptions of the viscosity implementation in polar coördinates).

In each of our simulations, we assume that the disc mass is small compared to the binary mass, $M_{\rm disc}\ll M_{\rm bin}$ so that we may neglect the self gravity of the gas, as well as changes in the binary’s orbital elements due to torques from the gas. In this work, we further focus on the equal-mass binary case ($q=1$).

In order to ensure that a quasi-steady state has been established over the spatial region of interest in our simulations, we evolve the system for 2.5 viscous timescales at the radius of the cavity wall at $r\approx 2a$,

[TABLE]

where $\Omega(r)$ is the Keplerian angular frequency in the disc at the radial distance $r$ from the binary’s center of mass and $t_{\rm bin}$ is the binary’s orbital period.

2.2 Accretion Prescription

To mimic accretion onto the individual BHs, we follow the sink prescription used in Farris et al. (2014) and Farris et al. (2015), but with some refinements. At each time-step, the gas surface density in cells within a distance $r_{\rm sink}$ (measured from each individual BH) is reduced uniformly by a fraction $\frac{dt}{\tau_{\rm sink}}$, where $dt$ is the time-step, and the sink time-scale $\tau_{\rm sink}$ is a parameter to control the mass removal rate, defined as

[TABLE]

In this work, we use $r_{\rm sink}=0.1a$ as the fiducial value. In our previous work (Farris et al., 2014), $\tau_{\rm sink}$ was set at the local viscous time based on the $\alpha$ viscosity law, applied to the individual minidiscs,

[TABLE]

Here $r$ is the distance to each individual accreting BH (analogous to Eq. 3 for the circumbinary disc), $u=1/2$ is the ratio of the BH mass to the total binary mass, and $t_{\rm bin}$ is the binary orbital period (which is $2\pi$ in code units).

In this work we choose $\tau_{\rm sink}$ to be a constant inside the sink radius, and perform a suite of simulations with a range of values of $\tau_{\rm sink}$. This allows us to examine the relationship between the binary-disc torques and the mass removal rate in the sink. When we model the minidisc as a thin $\alpha$ disc, we adopt $\tau_{\rm sink}=t_{\rm vis}(r_{\rm sink})$ as the sink time-scale.

2.3 Measuring Torques

In our simulations, the mass and velocity of the black hole binary are both fixed. In other words, we calculate the torque exerted on the BHs by the circumbinary disc, but do not include the impact of these torques on the binary’s orbital elements. This is justified as long as the torques are small enough so that the binay’s orbit does not evolve during the simulation.

Similar to Roedig et al. (2012) and Dürmann & Kley (2017), the total torque on the binary is computed by tracking both the gravitational and accretion torques,

[TABLE]

The gravitational torque is given by summing the gravitational force $f_{\rm gr}$ exerted by the fluid elements in each computational cell on the individual BHs,

[TABLE]

For the accretion torque, we measure the linear momentum of the gas removed in each simulation cell inside the sink, relative to the accreting BHs,

[TABLE]

where $v_{i}$ and $v_{\rm bh}$ are the velocities of the gas and of an individual BH (as measured in the inertial frame). With the above definitions, $T_{\rm tot}$ denotes the time derivative of the specific angular momentum $\ell$ of the binary,

[TABLE]

where $M_{\rm bin}=M_{1}+M_{2}$ is the total binary mass ($M_{\rm bin}\equiv 1$ in code units, as stated above).

3 Results

Snapshots of the surface density from simulations with fast ($\tau_{\rm sink}=0.0125$, or 0.002$t_{\rm bin}$) and slow ($\tau_{\rm sink}=5.0$, or 0.64$t_{\rm bin}$) mass removal rate are shown in Figure 1. As discussed below, the former sink is unphysically fast and represents a numerical experiment; the latter value is chosen to correspond to $\alpha=0.1$. In both cases, we confirm the findings in many previous studies (beginning with Artymowicz & Lubow 1994, 1996 and with recent works including, e.g. MacFadyen & Milosavljević 2008, Cuadra et al. 2009, D’Orazio et al. 2013, Shi et al. 2012, Noble et al. 2012, Roedig et al. 2012, and Farris et al. 2014), namely that a low-density, lopsided cavity forms around the binary, while narrow accretion streams connect the cavity wall to the minidiscs around the individual BHs. Most of the material in the narrow streams does not accrete onto a BH, but rather is flung outward. The outward going streams create an overdense lump at the cavity wall. As the figure shows, the more rapid sink (top panel) depletes the gas near the BHs, as well as in the entire inner region of the circumbinary disc, compared to the slower sink (bottom panel).

In Figure 2 we show the time-averaged surface density profiles in a frame co-rotating with the binary. We find again that with a faster sink, the surface density of the minidiscs, of the accretion streams and of the circumbinary discs near the cavity wall, are all visibly reduced. The three dashed circles in the top panel have radii of 0.6, 1.05 and 1.95$a$. In the annulus between the latter two circles, the azimuthal phase of the accretion stream trails behind the binary. As a result, the gravitational force between the binary and these portions of the accretion streams exert a negative torque on the binary (as further illustrated in Figure 4 below). By comparison, in the annulus between the first two circles, the tails of the narrow streams are ahead of the individual BHs, and therefore exert a positive gravitational torque on the binary. Comparing the shape of the minidiscs in the top and bottom panels, we find that the minidiscs in the fast-sink case are less dense overall; they are furthermore less round and symmetric than in the slow-sink case. Finally, an important conclusion from this figure is that the fast sink preferentially depletes the gas density behind (versus ahead of) the individual BHs (i.e. in the annulus $0.6<r<1.05$).

As discussed in several previous works (MacFadyen & Milosavljević, 2008; Roedig et al., 2012; Shi et al., 2012; D’Orazio et al., 2013; Farris et al., 2014), most of the material in the narrow streams gains angular momentum from the faster moving BHs and is flung back towards the cavity wall. In Figure 3, we plot contours of specific angular momentum in runs with both fast and slow sinks (top and bottom panel, respectively). These confirm that the specific angular momentum of the material forming the narrow streams is larger than that of the cavity wall, and also that overall, the streams have a net outward-going velocity. This suggest that the narrow streams could extract angular momentum from the binary and deposit it into the circumbinary disc, effectively acting as media to carry angular momentum from the binary to the gas disc.

In Figure 4, we plot the 2D distribution of the gravitational torque surface density, and in Figure 5, we plot the radial distribution of the azimuthally-averaged gravitational torque density. The previous works by MacFadyen & Milosavljević (2008) and Farris et al. (2014) calculated the torques only in the region $r>1.0a$. Our results are qualitatively similar to these studies in this region. In particular, we find significant negative torques from the disc material in the annulus between $a<r<2a$, corresponding to the extended faint blue streams in Figure 4. The total torque from the outer regions $r>a$ is also negative, with $r>3a$ contributing negligibly to this total. However, we find a significant positive torque in the annulus $0.6a<r<a$. Figure 4 further quantifies the point mentioned above, namely that this positive torque can be attributed to the compact tails of the narrow streams. These can be seen as the prominent red clumps in the figure, located near and ahead of the individual BHs. In the fast sink case, the contribution from this region dominates, making the total torque positive. In the slow sink case, the negative torques from the more massive minidiscs roughly cancel the positive torques from the tails of the narrow streams, and the result is a net negative torque.

In Figure 6, we plot the radial angular momentum flux profiles for both fast and slow sinks. We show, separately the fluxes corresponding to advection (blue), and to the viscous (green) and gravitational torques (cyan). The accretion torques are not shown explicitly, but correspond to a source term near the BHs (they are negligible in the lower panel, but account for the visible constant offset between the red and cyan curves at $r\approx 0$ in the upper panel). With the fast sink, the net torque exerted on the binary is positive, implying that angular momentum is transported from the disc to the binary. With the slow sink, the situation is the reverse, and angular momentum is transported outward from the binary to the disc.

We find that inside the cavity, the viscous angular momentum flux is close to zero, and advection is the only effective way to transfer angular momentum. With a slow sink, the advected angular momentum flux turns positive in the cavity. This implies that the narrow streams in the cavity are able to carry angular momentum outward to the cavity wall. When compared to the analogous results from a single-BH reference disc (shown by the dashed curves), we find that both the viscous and advected fluxes shift up or down based on the direction of the net angular momentum flux. However, at large radii, the viscous flux converges to the single-BH solution faster than the advected flux.

Overall, we conclude that when the net torque on the binary is negative (bottom panel), the angular momentum of the binary ends up being deposited into the circumbinary disc, and transported outward first by gravitational torques and advection (in the inner regions; r\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}3a) and then by the elevated viscosity (farther out; r\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}3a). For the fast sink (top panel), the behavior is the opposite, with reduced viscosity enhancing the inward advection of angular momentum, which ultimately results in a net transfer of angular monentum from the disc onto the binary.

In Figure 7, we show the net specific torque on the binary, in units of $\dot{M}$, as a function of the sink time scale $\tau_{\rm sink}$. We find a clear linear relationship with longer sink timescales yielding more negative specific torque. The total torque exerted on the binary is well approximated by the equation:

[TABLE]

When $\tau_{\rm sink}$ is below a critical value of $\tau_{\rm sink,crit}=2.1$, the binary feels a net positive torque from the disc, and is driven apart by the gas.111We here ignore the negative torques from gravitational waves, which could overwhelm the positive gas torque for sufficiently compact binaries. In the bottom pannel of Figure 7, we show the total torque as a function of the corresponding effective $\alpha$ (using Eq. (5) to convert $\tau_{\rm sink}=t_{\rm visc}$ to $\alpha$). The effective viscosity parameter corresponding to $\tau_{\rm sink,crit}=2.1$ is $\alpha_{\rm crit}=1.42$. We conclude that for sink time-scales corresponding to physically realistic choices of viscosity, binaries are always driven towards merger, rather than pushed apart, by the presence of the circumbinary gas disc.

Kocsis et al. (2012a, b), Rafikov (2013) and recently Rafikov (2016) have used simplified 1D models to study how the surface density profile of a steady accretion disc responds when angular momentum is injected to (or extracted from) an $\alpha$ disc at some small radius by a binary companion (see also Syer & Clarke 1995 for related earlier work). These analytic 1D steady-state solutions resemble standard Shakura-Sunyaev discs, but the disc adjusts itself to transfer angular momentum outward (or inward) at the same rate as the injection (or extraction) rate imposed at the bottom (i.e. at a small radius, comparable to the binary separation). This modified angular momentum transfer rates requires a surface density ’pile-up’ (or depletion) in the inner disc regions.

The form of this modified disc profile can be approximated by combining equations (3) and (11) of Rafikov (2016) (see also his Figure 2). The modified surface density follows the form

[TABLE]

where $A$ and $B$ are constants, related to one another as

[TABLE]

Here $\dot{J}$ is the angular momentum flux absorbed by the binary (which can be either negative or positive), and $\dot{M}$ is the mass accretion rate which is always positive.

In Figure 8, we show the time- and azimuthally-averaged surface density profiles from our simulations, for both a fast (top panel) and a slow sink (bottom panel). We find that outside the cavity wall, $\Sigma(r)$ accurately follows the form given in Eq. (12). In Figure 8, we further test the 1D predictions, and plot the empirical relationship between $\langle\dot{J}\rangle/\langle\dot{M}\rangle$ and the $B/A$ value obtained from fitting the profiles found in the simulations. We find

[TABLE]

which matches the prediction (eq. 12) quite well. Nevertheless, we caution that in order to reach a steady-state throughout the disc, we would need to run the simulations for several viscous times measured at the outermost disc radius. Thus, Figure 8 should be interpreted only as evidence for a “temporary” steady state in the inner regions, and for the sink time scale affecting the global surface density structure as expected, rather than a prediction for the long-term steady-state surface density profile, which would be established on a much longer time-scale.

Ryan & MacFadyen (2017) have performed a general relativistic local simulation of an isolated minidisc, subject to the tidal field of a companion. Their simulation does not include gas viscosity, but shows gravitational torques arising from non-axisymmetric structures (spiral waves) inside the minidisc, induced by the tidal torques of the companion BH. They find that the minidisc has an effective “gravitational Shakura-Sunyaev $\alpha$” on the order of a few $\times 10^{-2}$. Since this effect does not significantly enhance the effective $\alpha$ parameter of the minidisc, we proceed to assume the physical sink time scale $\tau_{\rm sink}$ equals the viscous time at $r_{\rm sink}$.

To investigate whether the sink radius $r_{\rm sink}$ affects our results, we have performed simulations with $r_{\rm sink}=0.1a$ and $0.05a$. In each case, we have set the sink time scale to $\tau_{\rm sink}=t_{\rm vis}(r_{\rm sink})$. Figure 9 shows the specific torque, in units of $\dot{M}$, as a function of time in these two simulations. The figure shows that the torques converge rapidly to a value that is insensitive to the choice of $r_{\rm sink}$. In both runs, we find that the gravitational torque dominates over the accretion torque, with $|T_{\rm acc}|/|T_{\rm tot}|\approx 0.01$.

In summary, in Table 1 we list the simulation parameters and the gas disc torques exerted on the binary in each of our simulations.

4 Discussion

4.1 Implications for the SMBH binary population

In our simulations the binary orbital period is $2\pi$, and the accretion rate $\dot{M}$ is arbitrary, depending linearly on the surface density normalization $\Sigma_{0}$. By scaling the simulation to $t_{\rm bin}=1$yr, and $\dot{M}=0.3\dot{M}_{\rm Edd}$ (with $\dot{M}_{\rm Edd}\equiv L_{\rm Edd}/\epsilon c^{2}$, $\epsilon=0.1$ i.e. the accretion rate corresponding to the Eddington luminosity), we obtain the binary’s residence time $t_{\rm res}\equiv a/[da/dt]$ in physical units. We find that this time-scale due to the gas torque is $t_{\rm res}=3.26\times 10^{6}$yr. In our code units ($M_{\rm bin}=a=1$, $t_{\rm bin}=2\pi$) the conversion constant for the accretion rate is $t_{\rm bin}/(2\pi M_{\rm bin})$. Therefore if $\dot{M}/\dot{M}_{\rm Edd}$ is fixed, in our code units $\dot{M}$ is proportional to $t_{\rm bin}$ and independent of $M_{\rm bin}$. In code units the residence time $\widetilde{t}_{\rm res}\varpropto T_{tot}^{-1}\varpropto\Sigma_{0}^{-1}\varpropto\dot{M}^{-1}\varpropto t_{\rm bin}^{-1}$ and the physical residence time $t_{\rm res}=\widetilde{t}_{\rm res}\times t_{\rm bin}/2\pi$. Thus the inferred residence time in physical units does not depend on $M_{\rm bin}$ or $t_{\rm bin}$. The same conclusion would follow by noting that the residence time in physical units would scale as $t_{\rm res}\propto$ (angular momentum/torque) $\propto a/(\Sigma_{0}t_{\rm bin})$, with $a$, $\Sigma_{0}$, and $t_{\rm bin}$ in physical units. However, the condition $\dot{M}=const\times\dot{M}_{\rm Edd}$ implies $\Sigma_{0}\propto a/t_{\rm bin}$, making $t_{\rm res}$ a constant. This counter-intuitive conclusion is a result of our scale-free simulation set-up, and would no longer hold in the presence of realistic cooling/heating or other physics introducing a physical scale.

Further choosing a total binary mass $M_{\rm bin}=10^{6-9}~{}{\rm M_{\odot}}$, the inspiral time-scale due to gravitational wave emission is $t_{\rm GW}=5.2\times 10^{8-3}$yr (Peters, 1964). We conclude that the scaled gas-driven inspiral time, implied by our simulations, is comparable to (or shorter than) the GW-driven inspiral time for binaries with total mass M_{\rm bin}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}10^{7}~{}{\rm M_{\odot}}.

Pulsar timing arrays are sensitive to $\sim$nHz-frequency gravitational waves generated by SMBHBs with orbital periods between $\sim$0.1-10 yr. Because gas torques can promote orbital decay of the binary much more efficiently than pure GW emission, the expected number of binaries, and the GW background (GWB) at these separations should be reduced accordingly (Kocsis & Sesana, 2011; Kelley et al., 2017). Although we find that for the $10^{8-9}~{}{\rm M_{\odot}}$ binaries which dominate the GWB, the gas torques are subdominant, we expect the impact of the gas torques to increase for smaller mass ratios (because gas torques increase, whereas GW torques decrease, for smaller $q$; e.g. Haiman et al. 2009). Syer & Clarke (1995) (see also Kocsis et al. 2012a, b) have proposed analytical models analogous to Type-II migration, in which migration rates are tied to the disc viscous time, but are slowed down by a factor related to the ratio $M_{\rm 2}/M_{\rm disc}$, when this ratio is above unity.

For a steady-state $\alpha$ disc with $M_{\rm bin}=2\times 10^{6-9}M_{\odot}$, and accretion rate $\dot{M}=0.3\dot{M}_{\rm Edd}$, the viscous time $t_{\rm vis}$ at $r=0.5a$ is approximately $(8\times 10^{4}$-$10^{6}$) yr, and the predicted residence time $t_{\rm res}$ is approximately (50-3)$t_{\rm vis}\approx(4-3)\times 10^{6}$ years (see Haiman et al. 2009 for a detailed calculation). This turns out to be very close to the result $t_{\rm res}\sim 3\times 10^{6}$ years we infer here by scaling the torques measured in our simulations. We emphasize that this good agreement appears to be a coincidence, because in the 1D models, the torques are assumed to arise from the overdense cavity wall of the circumbinary disk, whereas in our simulations, the torques are dominated by the gas near the individual BHs.

We note further that a steady-state $\alpha$ disc with an accretion rate of $0.3\dot{M}_{\rm Eedd}$ has a scale-height of $H/r\sim 10^{-3}$, and is much thinner than the disc in our simulations $(H/r=0.1)$. While simulating discs numerically with $H/r\sim 10^{-3}$ is currently prohibitively challenging, our choice of $H/r=0.1$ is typical for the values employed in previous works (Artymowicz & Lubow, 1994; MacFadyen & Milosavljević, 2008; D’Orazio et al., 2013; Farris et al., 2014) and allows for a direct comparison with those works. We expect that the dimensionless torques measured in the simulations will depend on $H/r$, and we intend to investigate this dependence in future work.

We have found that the sink prescription and time-scale affect the magnitude (and even the sign) of the torques on the binary. Slower sinks result in more negative torques on the binary, while faster sinks result in diminished negative (or even positive) torques. During the late stages of the merger, the Schwarzschild radius $r_{S}$ is not much smaller than the binary seperation. A more physically realistic sink prescription during this stage would be to adopt the dynamical time as the sink time-scale $\tau_{\rm sink}$, and the innermost stable circular orbit ($r_{\rm ISCO}$) as the sink radius $r_{\rm sink}$. However a more careful treatment of relativistic effects, e.g. through the use of the Paczynski-Wiita potential (Paczyńsky & Wiita, 1980) is required in this case.

4.2 Comparison with previous works

The nature of the torques and the corresponding binary evolution, arising from the disc-binary interaction, has been previously investigated by Cuadra et al. (2009), Roedig et al. (2012) and Miranda et al. (2017) as mentioned in the Introduction. Cuadra et al. (2009) and Roedig et al. (2012) have both simulated self-gravitating gaseous discs of mass $M_{\rm disc}=0.2M_{\rm bin}$ and BH mass ratio $q=1/3$. The disc is Toomre unstable, and self-gravity acts as a source of viscosity that can transport angular momentum outwards. In their SPH simulations, gas particles in the sink are absorbed by the black hole immediately, which corresponds to a sink time scale $\tau_{\rm sink}=0$. Similar to our results, Roedig et al. (2012) find that the gravitational torque changes sign at $r=a$ (positive at $r<a$ and negative at $r>a$). They also find that the accretion torque alone would shrink the binary seperation significantly, while in our work the accretion torque drives the binary apart, but is generally very weak, and unimportant compared to the gravitational torque (a \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}1\% correction). Furthermore, in Roedig et al. (2012), the binary is driven together by the disc, and the residence time is approximately 8000 binary orbits (taking their run closest to ours, i.e. “iso10” in their Fig. 3).

These differences are partly explained by the fact that in Roedig et al. (2012), the typical accretion rate is much higher than the Eddington rate (20-40 $\dot{M}_{\rm Edd}$ for the secondary and 4-7 for the primary when $M_{\rm primary}=2\times 10^{6}{\rm M_{\odot}}$ and $a=0.039$pc). We have found that the total torque exerted on the binary sensitively depends on the sink timescale, and with a fast sink the binary is driven apart by the gas. If we re-scale our simulations to the same $\dot{M}$, $M_{\rm bin}$, and seperation $a$, and we adopt $\tau_{\rm sink}=t_{visc}$, the residence time we obtain is approximately 250 binary orbits, which is around 30 times smaller. This difference may come from the different binary mass ratio, equation of state, disc thickness, as well as from the effect of self-gravity. Moreover, we suspect that for a more physical (i.e. less rapid) sink prescription was employed, the migration rate $-\frac{da}{dt}$ in Cuadra et al. (2009) and Roedig et al. (2012) would be faster, and closer to our value.

Finally, we note that Miranda et al. (2017) perform a series of 2D viscous hydrodynamics simulations of circumbinary discs, but with the inner regions ($r<a$) excised from the simulation domain. They conclude that the binary feels a net positive torque. This positive torque arises from accretion, and is a direct result of the assumption that the specific angular momentum of the dic gas near $r\approx a$ ends up being accreted by the binary. As we have seen in our simulations which resolve the inner regions, this is not the case. Rather than being deposited into the binary, the large angular momentum of this gas is returned to the disk, via narrow streams, in a ’gravitational slingshot’ mechanism.

4.3 Caveats and future work

In future work, we intend to relax several assumptions, while expanding the parameter space of our simulations. In this paper, we have focused on a mass ratio $q=1$, viscosity $\alpha=0.1,$ and disc aspect ratio $H/r=0.1$. In the future, we intend to to study how these parameters affect the torques and the evolution of the binary’s orbital parameters. In this paper we have further assumed a locally isothermal equation of state, with the temperature set to the same value for the circumbinary disc and for both minidiscs. This assumption should be relaxed by incorporating a more sophisticated treatment of the thermodynamics, so that the minidiscs are allowed to heat or cool independently (e.g. Farris et al., 2015). Realistic circumbinary discs will also have a range of thicknesses, and we intend to investigate the dependence of gas torque on $H/r$ (e.g. Ragusa et al., 2016). In this paper, the binary moves on a fixed, prescribed circular orbit. In future work, this assumption should also be relaxed, in order to measure the growth of eccentricity due to the binary-torque interaction (e.g. Cuadra et al., 2009) and self-consistently couple the evolution of the binary and the disc (e.g. Dürmann & Kley, 2017). If significant eccentricity is imparted to the binary, it may have observational consequences for GW observations, especially for the GWB measured by PTAs. Finally, our simulations are in 2D. We expect that the 3D vertical structure will modify the structure of the minidisc and the accretion streams, and inevitably have a strong affect the gas torques.

5 Conclusions

We have performed hydrodynamic simulations of thin, circumbinary accretion discs using a moving-mesh, finite volume code with the BHs present on the simulated grid. Our discs are locally isothermal, co-rotating in the plane of the binary, and we employ an $\alpha$ viscosity prescription. The binaries follow a fixed, prescribed circular orbit. We carried out a detailed study of the gas torques exerted on the BHB. We performed a numerical experiment, to study the effect of the mass removal rate in the sink prescription used to model accretion onto the individual BHs, on both the gas torques and on the global circumbinary disc structure. The most important conclusions of this work can be summarized as follows:

(1) We find a strong reciprocal relationship between the mass-removal time scale and the torques on the binary. The total torque is positive when the mass-removal time scale $\tau_{\rm sink}$ is $\lesssim 2.5$. For slower sinks, gas accumulates near the BHs and generates large negatives torques, making the net total torque exerted on the binary negative. The relationship between $\tau_{\rm sink}$ and total gas torque $T_{\rm tot}$ can be well approximated by a linear regression. When we model the minidiscs as thin $\alpha$ discs, and adopt a viscous time as the sink time-scale, the gas torque exerted on the binary is always negative for $\alpha\leq 1$, and insensitive to $r_{\rm sink}$.

(2) The gas torques are predominantly gravitational, and determined by gas dynamics near individual BHs. The gas torques come from threee regions: the narrow streams, tails of the narrow streams near the minidiscs, and the minidiscs themselves. The tail of the narrow stream contributes positive torque, while the other regions contribute negative torque. Rapid sinks make the tail of the narrow streams dominate, and result in a net positive total torque. For slower sinks, the negative torque from the material in the minidiscs behind the BHs cancel the positive torques from the gas ahead of the BHs, and result in a net negative torque. When we adopt the viscous time as the sink time-scale, the accretion torques are at most a few percent of the total torque.

(3) When are simulations are scaled to physical units, we find that binaries are driven together on a time scale of $t_{\rm res}=a/|da/dt|\sim 3\times 10^{6}$ years. This time-scale would be independent of $M_{\rm bin}$ or $t_{\rm bin}$, if the dimensionless torques were independent of the disc aspect ratio $H/R$. However, this is unlikely, and the migration rates will in general depend on $H/R$ (and then also on $M_{\rm bin}$ or $t_{\rm bin}$). Our current results suggest that the scaled gas-driven inspiral time is comparable to (or shorter than) the GW-driven inspiral time for binaries with total mass M_{\rm bin}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}10^{7}~{}{\rm M_{\odot}}. However, the above dependencies, as well as the caveats listed in the previous section, will have to be explored in future work, in order to accurately assess the impact of gas torques on the binary population.

Acknowledgements

We thank Daniel D’Orazio, Paul Duffell, and Geoff Ryan for useful discussions. Financial support was provided from NASA through the ATP grant NNX15AB19G (ZH). ZH also gratefully acknowledges support from a Simons Fellowship for Theoretical Physics.