© ESO, 2012

1. Introduction

In the currently favoured hierarchical framework of structure formation (White & Rees 1978), galaxies evolve through a complex sequence of merger and accretion events, and the existence of massive black holes (BHs) at their centres is nowadays a well established observational fact (see Gültekin et al. 2009, and references therein). By combining together these two pieces of information, the formation of a large number of massive black hole binaries (BHBs), following galaxy mergers throughout cosmic history, is a natural consequence of the structure formation process (Begelman et al. 1980). Although this is corroborated by several observed quasar pairs at a  ~ 100 kpc projected separation (Hennawi et al. 2006; Myers et al. 2007, 2008; Foreman et al. 2009; Shen et al. 2011), and by few  ≲ kpc dual accreting BHs embedded in the same galaxy (e.g. Komossa et al. 2003; Fabbiano et al. 2011), identification of gravitationally bound BHBs in galaxy centres remains elusive (for an up-to-date review on the candidates see Dotti et al. 2012, and references therein).

Even though observationally there is little evidence for their existence, much theoretical work has focused lately on Keplerian massive BHBs residing in galactic nuclei. One reason is that a deep understanding of the interplay between BHBs and their dense (stellar and gaseous) environment is required to predict robust signatures that may allow their identification. Additionally, it is still unclear how Nature bridges the gap between the two theoretically well understood stages of BHB evolution: (i) the dynamical friction driven stage, when the two BHs spiral in toward the centre of the merger remnant down to pc separations and (ii) the final inspiral driven by gravitational waves (GWs), which become efficient when the two BHs are at a separation  ≲ 10^-2 pc. Both dense stellar and gaseous environments have been shown to be effective in extracting the binary energy and angular momentum (see, e.g., Escala et al. 2005; Dotti et al. 2007; Cuadra et al. 2009; Khan et al. 2011; Preto et al. 2011), likely driving the system to final coalescence (an extensive discussion on the fate of sub-parsec BHBs can be found in Dotti et al. 2012). Scenarios involving cold gas are particular appealing not only because they might produce distinctive observational signatures, but also because cold gas dominates the baryonic content in most galaxies at redshifts higher than one, providing a natural reservoir of energy and angular momentum to drive the BHB towards coalescence.

Numerical simulations of wet galaxy mergers (e.g., Mihos & Hernquist 1996; Mayer et al. 2007) have led to the following picture for the post merger evolution. Cold gas is funnelled to the central  ≈ 100 pc by gravitational instabilities, where it forms a puffed, rotationally supported circumnuclear disc. Disc-like structures are actually observed in ULIRGs which are thought to be gas rich post-merger star-forming galaxies (e.g., Sanders & Mirabel 1996; Downes & Solomon 1998; Davies et al. 2004a,b; Greve et al. 2009). In the models, the two nuclear BHs efficiently spiral to sub-pc scales owing to dynamical friction against the massive circumnuclear disc (Dotti et al. 2007), eventually opening a cavity (or hollow) in the gas distribution (Goldreich & Tremaine 1980). The subsequent evolution of the system is determined by the efficiency of energy and angular momentum transfer between the BHB and its outer circumbinary disc.

The investigation of coupled disc-binary systems has a long-standing tradition in the context of planetary dynamics (Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Ward 1997; Bryden et al. 1999; Lubow et al. 1999; Nelson et al. 2000), where the focus usually lies on the extreme mass ratio situation, i.e., a star surrounded by a circumstellar disc, with a planetary companion embedded in it. The comparable mass case has also been extensively investigated in the context of binary star formation (Artymowicz & Lubow 1994; Bate & Bonnell 1997; Günther & Kley 2002; Günther et al. 2004), where a boost of activity has been triggered by imaging of nearby young binary stars embedded in hollow circumbinary discs (Dutrey et al. 1994). More recently, the techniques adopted in these fields have been applied to the BHB case. In the last decade, several investigations were devoted to the study of comparable mass BHBs evolving in circumbinary discs, exploiting a variety of analytical and numerical techniques (Ivanov et al. 1999; Armitage & Natarajan 2002, 2005; Hayasaki et al. 2007, 2008; MacFadyen & Milosavljević 2008; Hayasaki 2009; Cuadra et al. 2009; Lodato et al. 2009; Nixon et al. 2011; Roedig et al. 2011; Shi et al. 2012). However only few of them focused on the details of the dynamical disc-binary interplay. MacFadyen & Milosavljević (2008, hereafter MM08) made use of two-dimensional grid-based hydrodynamical simulation to study a BHB embedded in a thin α-disk (Shakura & Sunyaev 1973). They showed that the gas flowing through the cavity increases the energy and angular momentum transfer between the disc and the binary. Shi et al. (2012, hereafter S11) confirmed that result using full 3D magnetohydrodynamics (MHD) simulations. However, in both studies the BHB was on a fixed circular orbit, the central region of the disc (i.e., within twice the binary separation) was excised from the computational domain, and the disc self-gravity was neglected.

We employ full 3D smoothed particle hydrodynamical (SPH) simulations to study the interaction between BHBs and their surrounding discs. Our goal is to give a detailed description of the coupled disc-binary dynamics, paying particular attention to the competing effects of disc-binary gravitational torques and gas accretion onto the BHs in the evolution of the binary angular momentum budget. We simultaneously evolve the disc and the two BHs in a self-consistent way. This enables us to separately investigate the effect of gravitational torques coming from different disc regions and to directly link different physical mechanisms (gravitational torques and accretion) to the evolution of individual quantities describing the binary (mass, mass ratio, semi-major axis and eccentricity). This is different than previous studies (MM08, S11), where the binary was modelled as a fixed forcing quadrupolar potential and the central region was excised. To make sure that our results are not spuriously dependent on particular choices of sensitive parameters, we performed different runs varying physical and numerical prescriptions that may play a relevant role in the disc-binary energy and angular momentum exchanges. Specifically, we checked the possible dependence of our results on the “numerical size” of the two BHs (i.e., their sink radii, see next section), and on the adopted equation of state (EoS) for the gas within the central cavity.

The paper is structured as follows: we first describe the numerical setup and initial conditions in Sect. 2, then we show the importance of both accretion and gravitational torques on the binary evolution in Sect. 3. We study in detail the origin and strength of the mutual disc-binary gravitational torques in Sect. 4 highlighting the influence of the thermodynamics prescription. In Sect. 5 we briefly discuss the temporal behaviour and the spatial geometry of the accretion, speculating on the evolution of the binary mass ratio and BH individual spins. Finally, we compare our results to previous work in Sect. 6 and discuss the implications of our findings and give our conclusions in Sect. 7.

2. Simulations

Fig. 1 Relaxed meridional density maps of the disc for the adia (top) and iso (bottom) runs; the black dots indicate the BHs, the axis are in units of the binary semi-major axis a. Note that this is not an azimuthal average but a vertical slice of the disc.

Fig. 2 Time and azimuth-averaged values of the disc scale-height H / r and the disc eccentricity edisc as a function of radius for the different simulations.

The model and numerical setup of this work is closely related to that of Roedig et al. (2011) and Cuadra et al. (2009), hence we only outline the key aspects in the following, and refer the reader to these two papers for further details.

We simulate a self-gravitating gaseous disc of mass M_d = 0.2  M around two BHs of combined mass M = M₁ + M₂, mass ratio q = M₂ / M₁ = 1 / 3, eccentricity e and semi-major axis a, using the SPH-code Gadget-2 (Springel 2005) in a modified version that includes sink-particles which model accretion on to the BHs (Springel et al. 2005; Cuadra et al. 2006). Moreover, the orbit of the BHB is followed very accurately by using a fixed small time-step and summing up directly the gravitational force from every other particle in the simulation (Cuadra et al. 2009). The disc, which is co-rotating with the BHs, radially extends to about 7a, contains a circumbinary cavity of radial size  ~ 2a and is numerically resolved by 2 million¹ particles. The numerical size of each BH is denoted as r_sink, the radius below which a particle is accreted, removed from the simulation, and its momentum is added to the BH (Bate et al. 1995). We do not scale r_sink by BH-mass, but use one fixed value for both BHs. For all runs, the size of r_sink is smaller than the Roche lobe of both BHs, respectively. A particle is accreted if its separation from either BH is smaller than r_sink, no other conditions are imposed. The gas in the disc is allowed to cool on a time scale proportional to the local dynamical time of the disc , where is the initial orbital frequency of the binary². To prevent it from fragmenting, we force the gas to cool slowly, setting β = t_cool / t_dyn = 10 (Gammie 2001; Rice et al. 2005). The choice of β is motivated by the requirement of maintaining the disc in a configuration of marginal stability (Toomre parameter Q in the range [1, 2], see Sect. 4.3.2 and Fig. 10). Assuming β = 10, we thus effectively create an environment where self gravity acts as a viscosity that can transport angular momentum outwards.

We use two different treatments for the thermodynamics inside the cavity, denoted as adia and iso, respectively. In the adia runs, the gas within the cavity is allowed to both cool via the β prescription and to heat up adiabatically, as it is the remainder of the disc (Cuadra et al. 2009). In the iso runs, we define a threshold radius r_cavity = 1.75a, below which the gas is treated isothermally, meaning that the internal energy per unit mass is set to be u ≈ 0.14(GM / R) (Roedig et al. 2011). These two numerical experiments were chosen to be highly idealised, allowing us to study how torques behave in two opposite scenarios: one in which hydrodynamics is important in the cavity – the gas can cool efficiently and settles down in the binary plane forming mini-discs (iso); and one in which hydrodynamical forces have very little influence inside the cavity and the dynamics is dominated by the gravity of the BHs (adia).

As initial conditions we use a relaxed snapshot from Cuadra et al. (2009) taken at their (see Roedig et al. 2011, Sect. 2.1)³. The nomenclature of the runs and the parameters used are listed in Table 1. Each run is evolved for  ≈ 90 binary orbits, and we store the output in single precision  ~ 6 times per orbit. We also closely track accretion by storing the position and velocity of the two BHs and of each accreted particle at the time the latter crosses the sink radius. For both prescriptions (adia and iso) we show the relaxed density configurations sliced along the meridional plane in Fig. 1 (see also Figs. 8, 9 for a face-on view)⁴ and the measured, time-averaged values for both the (semi) scale height of the disc H / r, and the eccentricity of the disc e_disc in Fig. 2. H is taken to be the height above and below the disc-orbital plane in which 70% of the mass inside the annulus (r,r + dr) is found. Generally, we can consider the physics to be resolved if the smoothing length h⁵ is smaller than the characteristic lengthscale: radially the criterion h / r ≪ 1 is achieved well: h / r ∈ (0.005,0.2); whereas vertically the criterion h < (c_s / Ω) is well satisfied throughout the simulation domain, except very close to the primary hole: h / (c_s / Ω) ∈ (0.1,0.5) _{| r ≥ 0.5}. Finally, our resolution ensures that h ≲ r_sink close to the BHs, a condition which is necessary to properly resolve accretion onto a sink particle.

Unless otherwise stated, we will present all the relevant quantities in the natural units of the simulation by setting G = M₀ = a₀ = 1. It follows that also the initial circular velocity of the binary is V_c,0 = 1 and its initial period P₀ = 2π. We will then discuss the astrophysical implications of our findings by scaling them to fiducial astrophysical BHB systems.

3. Binary evolution: causes

For each of the four runs, we measure the two primary orbital elements: eccentricity e and semi-major axis a. Their evolution is shown by the red lines in Fig. 3. In all four runs, we find ȧ(t) < 0 and ė(t) > 0. While the eccentricity evolution is largely independent on the sink radius value and the adopted EoS within the cavity, the orbital shrinking is much faster in the adia runs, in which the binary shrinks by 4 − 5% over 90 orbits. Conversely, in the iso runs, the two BHs get only about 1% closer. Linear extrapolation of such instantaneous shrinking rates at face value would imply binary coalescence on a time-scale of  ~ 2000P₀ and  ~ 10⁴P₀ for the adia and the iso runs, respectively. Whereas, assuming a constant eccentricity growth rate, the limiting value predicted by Roedig et al. (2011) would be reached in  ≲ 1000P₀. As a note of caution, for a very similar numerical setup, Cuadra et al. (2009) showed that the secular evolution of a and e departs from being linear for longer timescales. This has to be expected given the limited energy and angular momentum reservoir provided by an isolated circumbinary disc, which cannot transfer further out the angular momentum extracted by the binary. However, in a more realistic situation (e.g., a binary interacting with a steady inflow of gas), such slow- down may well not happen. How the long-term secular evolution is influenced by external replenishing of the disc is beyond the scope of this paper.

Fig. 3 Semimajor axis (top panel in each plot) and eccentricity (bottom panel in each plot) evolution for all the four runs. In each panel, the red line represents the full evolution directly taken from the simulation whereas the black line is the equivalent evolution when considering the energy and angular momentum exchanges due to accreted particles only (i.e., basically due the evolution of the mass and mass ratio of the binary).

3.1. Gravitational torque and accretion contribution to the angular momentum budget

Being interested in the physical mechanisms driving the binary evolution, we firstly identify two distinct processes:

• 1.

  the gravitational torques exerted by the gaseous particles ontoeach individual BH, T_G (gravity torque);

• 2.

  the accretion of instreaming particles crossing either BH sink radius, (dL / dt)_acc (accretion torque)⁶.

Conservation of the total angular momentum in the simulations implies (1)where L is the BHB orbital angular momentum vector⁷. All vectors are computed with respect to a Cartesian reference frame centred in the BHB centre of mass (CoM)⁸; the binary initially lies in the x–y plane with angular momentum oriented along the positive z axis. In our SPH simulations, T_G can be computed at each snapshot by direct summation of individual particle torques onto each BH yielding (2)where j runs over all N gas particles, k identifies the two BHs, r are position vectors, and m and M denote particle and BH masses respectively. On the other hand, (dL / dt)_acc can be computed by assuming instantaneous linear momentum conservation of each accreted particle yielding ΔL = r_k,j × m_jv_j. Here r_k,j is the position vector of the accreting BH at the moment of swallowing the particle j, and v_j is the particle velocity vector. Over the interval between two subsequent snapshots, we evaluate the binary angular momentum change, ΔL, by splitting Eq. (1) into two parts: (3)where Δt is the interval between the two snapshots, and the sum runs over all the accretion events j occurring during this time lapse at their respective positions. Note, that we assume an accretion event to be instantaneous, i.e. the accreting BH k does not move during the accretion.

Fig. 4 Consistency check of the evolution of the binary angular momentum components for all the simulations. In each plot, the Lx, Ly and Lz component evolution is shown from the top to the bottom. In each panel, the dashed line is the ΔL exchange balancing the gravity torques at each time step, the dotted line is the ΔL exchange computed from the accreted particles, and the solid black line is the sum of those two. For comparison, the red line is the angular momentum evolution taken directly from the raw simulation data. All ΔL are in units of  [M0a0Vc,0] .

We use Eq. (3) to assess the relative importance of gravitational torques and accretion in the evolution of the system. In Fig. 4, we plot the evolution of the x,y,z components of the angular momentum L as directly evaluated from the positions and velocities of the two BHs stored in the simulation outputs (red), together with the relative changes due to accretion (black dotted) and gravity torques (black dashed) according to Eq. (3). The solid black line is the sum of the latter two which should overlap with the red line if our decomposition is sufficiently accurate and these are the only two physical processes determining the binary evolution. For comparison, in simulation units, the initial binary angular momentum is  ≈ 0.18. The L_x and L_y components are almost constant, showing fluctuations at the 10^-4 level. The mismatch between the red and the solid black lines here is because simulation outputs are in single precision, i.e. accurate to the fourth digit, which is the fluctuation level of the computed quantity. On the contrary, the L_z component changes up to a 10^-2 absolute level (i.e., about 5% of the initial value), and in this case we see very good agreement between the two lines. Note that in three of the runs, the overall angular momentum grows over time, even though a shrinks and e increases. This counter intuitive result is due to the dominance of accretion in the angular momentum budget, and will be further discussed in Sect. 3.2 below (see also S11).

As a general rule, the accretion contribution to the binary L evolution is at least comparable to the effect of the gravitational torques. It is therefore also interesting to see the accretion contribution to the evolution of the binary orbital elements. Having stored all the accretion events, this can be done separately by simply evolving the binary from a₀ and e₀ only by adding one accreted particle after the other, imposing linear momentum conservation. This “would-be” evolution of e(t) and a(t), accounting only for the accretion effect, is shown in black in all panels of Fig. 3. It is clear that e(t) is mainly driven by gravitational torques, with accretion playing a negligible role. This is in line with our previous findings (Roedig et al. 2011), where the evolution of the eccentricity was attributed mainly to the gravitational interaction between the binary and the overdensities excited by its quadrupolar potential in the inner rim of the circumbinary disc. Conversely, the a(t) plots highlight the importance of accretion for the semimajor axis evolution. In the iso05 case, the binary shrinking due to accretion is actually larger than the total one, meaning that the disc torques would force the binary to expand; an effect similar to that described by Lin & Papaloizou (2012) in the context of planetary migration in self-gravitating discs. Note that in the adia cases ȧ is dominated by the effect of the disc. This explains the match between such simulations and the Ivanov et al. (1999) analytical model based on angular momentum transport through the disc (Cuadra et al. 2009).

3.2. Dissection of the binary evolution into its components

So far, we have separated the relative contribution of gravitational torques and accretion to the binary angular momentum budget. We now investigate how the angular momentum change is distributed among the relevant binary quantities. As clearly shown in Fig. 4, L_x and L_y show only small fluctuations (and therefore remain small compared to L_z, since the binary initially orbits in the x–y plane). From here on, we will therefore concentrate on the dominant L_z component. The binary angular momentum is (4)where μ = M₁M₂ / M is the binary reduced mass. Equation (4) can be differentiated to separate the contribution of each single relevant quantity to the angular momentum change: (5)By directly measuring a, e, M, μ from the simulation, we can evaluate each single term and decompose the L_z evolution accordingly. This is shown in Fig. 5 for all the four runs. Note that the sum of the four contributions (solid black line in each panel) closely matches the measured L_z evolution directly measured from the simulation (red), validating our first-order expansion.

It is worth noting that the increase of L_z is perfectly compatible with semi-major axis shrinkage and eccentricity growth. In fact, Eq. (5) can be inverted to obtain (6)where we used the fact that . Even if the total torque is positive, and the eccentricity grows, the binary can shrink because of the negative contribution of the M and the μ terms. In particular, as we will see below, the higher accretion onto the secondary hole results in a large increase of μ (long-dashed line in Fig. 5) which is the main driver of the binary angular momentum growth. Equations (5) and (6) highlight the importance of accretion in determining the binary temporal evolution.

Fig. 5 Decomposition of the Lz evolution for the four runs. The ΔLz contributions stem from a (dotted), e (dashed), M (dot dashed) and μ (long dashed) variations separately. The sum of the four contributions is given by the solid black line, whereas the solid red line is, as in Fig. 4, the angular momentum evolution taken directly from the raw simulation data. All ΔL are in units of  [M0a0Vc,0] .

4. Gravitational torques

As pointed out in the previous section, understanding secular dynamics of BHBs in gaseous environments is closely tied to studying the interplay between long distance gravitational forces, hydrodynamics and accretion. In this section, we focus on the torques coming from the self-gravitating gaseous disc, investigating the influence of regions at different distances from the BHB.

4.1. Time averaged torque profiles

Gravitational torques exerted by the disc onto the BHB can be directly calculated at each snapshot from Eq. (2). It is of great interest to understand where such torques originate, and given the cylindrical nature of the problem, it is natural to investigate their azimuthally-averaged radial distribution. We take r to be the projected radial coordinate (in the x–y plane) from the binary CoM, and we define dT / dr to be the differential torque, integrated over the azimuthal coordinate and the disc height. The total average torque exerted by gas in the projected distance interval  [a,b]  from the binary CoM is (7)where ⟨·⟩ denotes temporal averaging over the entire simulation.

In Fig. 6, we show the time averaged differential torque acting on the binary, ⟨dT / dr⟩ (blue), together with its integral according to Eq. (7) (black). We also decompose the former in the two components acting on the primary (green) and the secondary (red) BH. All torques are null at the binary corotation radius we therefore show the total torque by integrating from this point inwards (⟨T _[1,0] ⟩, binary region) and outwards (⟨T _[1,∞] ⟩, disc region).

Fig. 6 Differential torques dT / dr in units of averaged over the entire simulations. In each panel we show the differential torque on the primary (green), on the secondary (red), the sum of the two (blue), and the total integrated torque (black) according to Eq. (7). This latter is integrated starting from a inwards and outwards. Notably, the inward torque is positive, whereas the outward torque is negative.

Fig. 7 Each plot depicts the torques exerted by the disc on the BHB in the time domain (left panels) and in the frequency domain (right panels). In each plot, from the top to the bottom, pairs of panels refer to the five domains discussed in the text: total torque (i), torque exerted by the material located at r < a (ii), a < r < 1.8a (iii), 1.8a < r < 2.5a (iv) and r > 2.5a (v). In each left panel, the red line is the raw oscillating torque, and the blue one is the torque smoothed over three periods to show the average behaviour. In the associated right panel, we plot the power spectrum as a function of frequency in units of the initial binary orbital frequency. All torques strongly oscillate, showing several characteristic frequencies, however, note the different scales of the panels.

In all simulations, the local average torque shows an oscillatory behaviour with a sharp maximum at the location of the secondary BH, r ~ 0.75, and a deep minimum in the cavity region at 1 < r < 2. In the body of the disc (r > 2) positive and negative peaks alternate, but they are shallow and almost cancel out, giving a negligible contribution to the total torque (as witnessed by the fact that the integrated torque is basically constant for r > 2). Note that the torque on the secondary BH is always larger than the torque on the primary, due to its proximity to the outer disc resulting in a stronger interaction. The general behaviour is qualitatively the same in the adia and iso runs. In this latter case, however, we found a much sharper negative peak at r ≈ 1.7 followed by a smaller secondary negative bump at r ≈ 1.2. This appears to be an artifact of the sudden change in EoS of the gas inside the cavity, at r = 1.75. The net result is a smaller negative ⟨T _[1,∞] ⟩ which has important consequences on the binary evolution. We see in fact that in the iso runs ⟨T _[1,0] ⟩ and ⟨T _[1,∞] ⟩ almost cancel out, meaning that, overall, gravitational torques do not change the binary angular momentum. Conversely, in the adia runs ⟨T _[1,0] ⟩ is much smaller (in absolute value) than ⟨T _[1,∞]  ⟩, implying an efficient angular momentum transfer from the binary to the gas (MM08; Cuadra et al. 2009).

4.2. Spectral analysis

We now consider the torque evolution in time. We separately discuss torques coming from different disc regions by showing both the time series and their associated power spectra. We cut the spatial domain into the following radial annuli:

• (i)

  0     a < r < 10   a: the entire domain;

• (ii)

  0     a < r < 1      a: the “binary region”;

• (iii)

  1     a < r < 1.8a: the “cavity region”;

• (iv)

  1.8a < r < 2.5a: the “rim region”;

• (v)

  2.5a < r < 10   a: the “disc region”.

The associated time series and power spectra are shown respectively in the left and right panels of Fig. 7.

In all the simulations, the overall torque (i), shows a clear periodic oscillation, much larger in amplitude than its average value. The power spectrum unveils several distinctive peaks, with relative amplitudes that can vary significantly for different simulations. In particular, peaks in the iso simulations are much sharper and better defined. This is because in the adia simulations the BHB shrinks significantly, resulting in a broadening of the characteristic frequencies. Moreover, as we shall see in the next section, the disc sub-structures are much better defined in the iso runs, giving rise to neater features.

In the binary region (ii) the torque is mostly coherent and positive. Because the torque strength in (ii) is regulated by the mass inflow, it shows periodicities that are related to the accretion flow: a disc component around (corresponding to the disc peak density at r ≈ 2.5a), the forcing frequency of the binary at , and the beat between these two at (see Sect. 5 and Roedig et al. 2011, Sect. 5.1). In the cavity region (iii) the torque is negative on average but strongly oscillating. Several periodicities are detectable, the most striking being a peak at which appears again to be directly related to the binary period. In the rim region (iv) the torque is highly oscillating, and the strongest feature is a sharp peak at , whereas in the disc region (v) the only significant spectral component is at . As a general trend, moving from the inner region to the disc body, torques become incoherent (i.e. they average to zero) and strongly oscillating (compare the power spectra scales in the different panels of each plot).

4.3. Interpretation: torque origin and location in the disc

In this Section, we provide a global interpretation of the features observed both in the radial distribution of the time averaged torques and in their temporal evolution. The arguments discussed below are supported by Figs. 8–10 and by Table 2.

4.3.1. Origin of the positive and negative peaks

It seems natural to compare the torque radial profiles obtained by our simulations with linear perturbation theory, in which torque minima and maxima are connected to outer Lindblad resonances (OLRs). It is in fact tempting to associate the torque minimum at r ≈ 1.6a with the 2:1 OLR. We should however be careful in pushing this interpretation too far, since, as already pointed out by MM08 and S11, the assumptions of linear theory are not satisfied in this context. Most importantly, looking at the upper panels of Figs. 8 and 9, we notice that the region r < 2a is almost devoid of gas, and the streams are almost radial. Mass fluxes reported in Table 2 clearly show that, at any radius, there are always fluxes of ingoing and outgoing mass, resulting in a steady net inward flux consistent with the accretion rate onto the two BHs. A strict OLR interpretation of the torque minimum at r ≈ 1.6a would instead require particles in circular orbit at that radius, experiencing a secular effect due to the phase-coherent periodic forcing of the binary; this is not what happens within the cavity region. The strongest 2:1 OLR is certainly responsible for the evacuation of the gas close to the binary and for the formation and maintenance of a cavity, however cannot be directly responsible for the coherent torque seen in the cavity region. This is also supported by the fact that MM08 and S11 find a minimum at the same location (r ≈ 1.6 − 1.7a) for an equal mass binary, where the 2:1 resonance is absent (because of the symmetry of the forcing potential), and the location of the strongest OLR (3:2) would be at r ≈ 1.3a. The strong negative torque in the cavity region has a purely kinetic origin: material ripped off the disc edge forms well defined streams following the two BHs, which are clearly distinguishable in both the surface density plots shown in Figs. 8 and 9. The streams are responsible for the yellow tails following the two BHs at  ~1.5a in the torque density panels, which lead to a net negative torque. Conversely, at r ≳ 2a, we have a well defined, almost circular disc, and the torque density peaks at r ≈ 2a and r ≈ 2.5a can be identified with the loci of the strong 3:1 and 4:1 OLRs (Artymowicz & Lubow 1994).

Fig. 8 Top panel: typical instantaneous logarithmic surface mass density of the iso05 simulation. Bottom row: combined surface torque density exerted by the disc onto the binary (left panel), onto the primary BH (central panel) and onto the secondary BH (right panel). All plots are in code units: G = M0 = a0 = 1.

Fig. 9 Same as Fig. 8 but for the adia05 run. Here the features highlighted above are even more evident (especially the symmetric overdensities at the inner edge of the disc).

Our simulations also allow us to investigate the torques within the binary corotation radius at r < a, a region often excised in grid-based simulations (see MM08 and S11). Here we find strong positive torques on both BHs. This is because the infalling material approaches the BHs at super-Keplerian velocities, and bends in a horseshoe fashion, exerting a net positive torque in front of them. In fact, the maximum positive torque basically coincides with the location of the two BHs (sharp peak at r ≈ 0.75a for the secondary and a broader peak around r ≈ 0.3a for the primary, see Fig. 6). The very same effect, in the context of planetary migration, was discussed by Lin & Papaloizou (2012)⁹. Note that the positive torque is related to this “stream bending”, and not directly to the small discs of gas orbiting each BH; its magnitude is in fact similar in the iso and in the adia simulations, even though in the former, the mass in the minidiscs is significantly larger as shown by the time and azimuthally averaged surface density profiles in the upper panel of Fig. 10.

Fig. 10 Disc properties: average surface density (top panel), and Toomre parameter Q (bottom panel) as a function of r. Line style denotes iso10 (long dashed), iso05 (dot dashed), adia10 (solid) and adia05 (dashed).

4.3.2. Disc structure and characteristic torque frequencies

The surface density profiles in Figs. 8 and 9 highlight some clear differences between the iso and the adia runs: (i) the appearance and shape of the minidiscs; (ii) the amount of gas feeding the streams; (iii) the definition of the inner disc edge. In the iso runs, we find a large amount of gas forming mini-discs around both BHs and an almost empty cavity except for two tenuous streams connecting the outer edge to the mini-discs. The edge of the disc is quite circular and well defined. Conversely, due to the gas inside the cavity being hotter in the adia runs, the streaming activity is more violent with thick spirals that do not settle into well defined mini-discs, but rather form a diluted, ∞-shaped cloud around the binary. The edge of the disc is strongly disturbed by the ripped out streams and is usually not circular. The different streaming activity is also confirmed by numbers in Table 2, where we collect the average inwards and outwards mass fluxes at r = 1.5a and r = a. Mass fluxes are generally larger in the adia case. Note that this does not necessarily result in larger accretion, as outgoing fluxes are also larger in this case. Most of the material in the streams, does not end up in an accretion flow, but is accelerated back to the disc (a sort of “slingshot”) impacting on the inner edge, an effect already seen and extensively described by MM08 and S11. The amount of impacting gas is 50% larger in the adia case, possibly contributing to the inner edge destabilization and to the larger perturbation in the disc structure. The appearance of spirals in the disc is related to the behaviour of the Toomre Q parameter, shown in the bottom panel of Fig. 10, which quantifies the degree of self-gravity of a disc. Not surprisingly the disc develops a spiral pattern in the region 3a   ≲ r   ≲ 5a, where, in fact, we find that the Toomre parameter reaches its minimum value Q ≈ 1.5 (Cuadra et al. 2009). The shape of the spiral arms varies much in time, and their definition is different from simulation to simulation. However, at any time we find spiral structures with 2, 3 or 4 arms, which remain confined between 3a < r < 5a.

The disc structures highlighted above are reflected in the torques depicted in the lower panels of Figs. 8 and 9. The total torque shows a typical quadrupolar structure, which leads to rapid oscillations averaging to zero in the disc body. Conversely, as already discussed, in the inner cavity we can appreciate the yellow tails following the two BHs at  ~ 1.5a, leading to a net negative torque.

Fig. 11 Simple test model for the periodicity generated in the outside disc, i.e. at r > a. In the top panel we show the temporal evolution of the torque, whereas in the bottom panel we show the Fourier transform. In each panel, the blue curve is obtained by placing three point masses at distance 1.6a, 2.1a, 3.5a from a circular binary, and the one red is the torque found in the iso10 simulation.

The main peaks observed in the torque power spectra must therefore arise from the structures we just described. In particular we identified the overdensities at r ≈ 2a in the inner rim of the disc, and the spiral structures at r ≳ 3a in the main body of the disc. The main torque frequencies must come from the interaction between the forcing binary quadrupolar potential and such overdensities propagating in the disc. To test this hypothesis, we mimic the situation by placing test particles in circular orbits at r = 2a and r = 3.5a around a circular binary. We compute the torques and show them as blue lines in Fig. 11. The natural frequency between a quadrupolar potential oscillating with frequency f₁ and an overdensity orbiting at frequency f₂ is the beat frequency 2f₁ − f₂, and this is what we see in the Fourier spectrum. The beats between the binary and the particles at r = 2a and at r = 3.5a give rise to sharp peaks seen at and , respectively. The fact that such peaks are much sharper in the iso simulations stems from two facts: (i) the binary does not significantly shrink in the iso runs, limiting the broadening of the spectral feature and (ii) the disc features (disc edge and spiral arms) are much more defined in this case, and the torques are well localized. This can be also seen by comparing the much neater torque structure in the bottom panel of Fig. 8 to the blurry one of Fig. 9. The peak of the blue line in Fig. 11 is obtained by placing a third particle at r = 1.6a, a separation corresponding to the maximum negative torque inside the cavity. However, Fig. 7 shows that the peak at is stronger at r > 1.8a than in the a < r < 1.8a region, meaning that it must be mostly directly related to the periodicity at which the binary rips gas off the disc edge (i.e. the binary period), and not to a specific radius in the cavity. Finally, as already noticed, torques in the binary region are related to the mass fuelling the two BHs, and therefore show periodicities related to the binary (), to the inner rim of the disc () and to the beat between the two (). Note that the latter two are much more evident in the adia runs, where the disc rim overdensities feeding the accretion streams are more pronounced.

5. Accretion

In Table 2 we also show the average mass accretion rates. We recall here that a particle is considered accreted when it crosses the sink radius of one of the BHs, In this sense, r_sink defines a small excision region around each BH. However, almost all the accreted particles are bound and lie well within the Roche lobe of the accreting BH. The mass accretion rate on both BHs is almost independent on the sink radius in the iso simulations, where instreaming gas settles in well defined accretion discs progressively loosing angular momentum. In the adia simulations, this is also true for M₁, but the accretion rate onto M₂ scales linearly with the sink radius, suggesting a direct relation to the BH cross section (defined by the sink radius itself). As a consequence, we are likely overestimating accretion on this hole for the adia runs. We should state here, that we are not interested in whether the particle is accreted immediately or is dragged along in a minidisc: for the purpose of linear momentum transfer from the disc to the BH, the two processes are equivalent. The physics of the minidiscs, beyond resolution issues, is likely to be radiation dominated, and is certainly not well described in our simulations.

It is particularly interesting to compare the accretion rates to the mass flow rates. Firstly we notice that in all simulations the accretion rate and the mass flow rates at r = a and r = 1.5a are nearly the same, meaning that the system is in a steady state configuration. Interestingly, the mass accretion rate is  ~ 25% lower than the inflow rate at r = a. This means that not all the material crossing r = a is accreted, an assumption commonly adopted in grid simulations where the r < a region in excised. Part of the gas suffers a slingshot impacting back to the disc, affecting further the disc-binary angular momentum transfer.

Fig. 12 Accretion onto the two BHs. Top box: accretion rate as a function of time. The red line is the total accretion rate while the blue and black lines are the accretion rates onto the secondary and the primary respectively. Lower box: Fourier transform of the total accretion rate (lower panel) and of the torques (upper panel). In the upper panel, the green line is the total torque (highest peak normalized to 1) and the red line is the torque exerted by the material at r < a (multiplied by ten).

In the top panel of Fig. 12 we show the temporal evolution of the accretion rate on the primary hole, Ṁ₁, on the secondary hole, Ṁ₂, and the sum of the two, Ṁ. As already found in Cuadra et al. (2009); Roedig et al. (2011), Ṁ₂ is a factor  ~ 2 larger than Ṁ₁, due to the closer interaction between the secondary BH and the disc. The relative mass growth δM / M is therefore about six time faster for the secondary BH, implying the binary mass ratio will tend toward q ≈ 1 if this condition persists over the entire evolution of the system. Accretion rates are strongly modulated in time, and notably, Ṁ₂ shows a striking periodicity related to the binary orbital period, whereas Ṁ₁ is dominated by longer timescale fluctuation, related to the overdensities developing at the inner rim of the circumbinary disc.

The Fourier transform of Ṁ is given in the lower panel of Fig. 12. We observe the usual three distinctive peaks observed in the torques coming from the binary region, as shown by the central panel of the figure. It is clear that the accretion and inner-torque periodicities are intimately connected, both reflecting the temporal evolution of the mass inflow in the binary region (ii).

Finally we can check the orientation of the angular momentum vector of the accreted material in the reference frame of the accreting BH. This number quantifies the level of “coherence” of the accretion flow, and has important consequences for the individual BH spin evolution and the magnitude of gravitational recoil at the coalescence (Bogdanović et al. 2007; Dotti et al. 2010; Kesden et al. 2010; Volonteri et al. 2010; Lousto & Zlochower 2011; Lousto et al. 2012). At each snapshot we therefore compute the average L_z and L of the accreted particles in the accreting BH reference frame, and the angle θ = arccos(L_z / L) defining the degree of misalignment with respect to the z axis. Our results are similar to those already found by Dotti et al. (2009), although in the cited study the spin evolution was studied only before the gas scouring. More specifically, an isothermal EoS leads to extremely coherent accretion flows, with fluctuations in the angular momentum direction of the mini-discs of few degrees. In the adiabatic case, the orientation of the minidisc orbiting the secondary changes by at most  ≈ 20 degrees, and has an average excursion of 7 degrees, 50% bigger than the oscillations in the primary disc. Assuming that the BH spins have efficiently aligned with the circumbinary disc before the opening of the cavity (Dotti et al. 2010) this implies that the efficient spin-up of the binary will continue following cavity opening.

6. Comparison with previous works

The binary evolution and torque structure have been previously investigated by MM08 and S11, as mentioned in the introduction. In particular, we can compare results regarding the average gravitational torques and the mass accretion rates. While both of them can be expressed in several ways, we find it practical to normalize these quantities to the disc mass, M_d, initial binary period, P₀, and initial binary angular momentum magnitude, L₀. In these units we can write: \arraycolsep1.75pthere is the disc mass computed using the peak surface density¹⁰ Σ_p. In both MM08 and S11, r_p ≈ 3a, close to what we find in our discs (see Fig. 10). Before proceeding with this comparison we should mention several issues that have to be considered. Firstly L₀ = (M₀ / 4)(GM₀a₀)^1 / 2 for the circular equal mass binary simulated by MM08 and S11, whereas in our simulation with q = 1 / 3 we have L₀ = (3M₀ / 16)(GM₀a₀)^1 / 2 (assuming a circular binary). Secondly, in MM08 and S11, the concept of “accretion rate” is defined by measuring the mass flux through the inner boundary, placed at r = a for MM08 and r = 0.8a for S11. As already discussed (see Table 2)  ≈ 25% of the mass crossing r = a is accelerated back to the disc edge extracting angular momentum from the binary (an effect already seen and described extensively by MM08 and S11 for the material inflowing in the cavity but not reaching the excised region). We can therefore consider MM08 and S11 to overestimate the accretion rate by a factor of 25%. Lastly, MM08 and S11 can compute gravitational torques only for r ≳ a, outside the binary corotation radius¹¹. As shown by both studies (and independently recovered by us), gravitational torques exerted by the disc elements on the binary are negative in this region.

For the sake of the comparison we use the iso05 and the adia05 runs. In the torque computation, we find Γ₁ =  −3.5 × 10^-3 for the former and Γ₁ =  −5.3 × 10^-3 for the latter. This can be compared to  − 1.26 × 10^-3 of MM08 and  − 1.8 × 10^-2 of S11. Our gravity torques are a factor 3–4 larger than MM08 and a factor 4–5 smaller than S11. Concerning the accretion, we get Γ₂ = 1.1 × 10^-3 for the iso05 run and Γ₂ = 6 × 10^-4 for the adia05 run, to be compared to Γ₂ = 1.2 × 10^-4 of MM08 and Γ₂ = 6   × 10^-3 for S11. Our accretion rates are a factor 5–10 larger than MM08 and a factor 5–10 smaller than S11. The similar scaling between ⟨T _[1,∞] ⟩ and Ṁ is not surprising, since, as discussed in Sect. 4, the former is directly linked to the streams fuelling the BHs.

7. Discussion and conclusions

The evolution of BHBs in circumbinary discs remains largely an open issue. Here we carried out a detailed study of the torques acting on the BHB, including the effect of the outer disc, the mass flowing in the cavity and, for the first time, of the gas entering the BHB region and eventually accreted onto the two BHs. We paid particular attention to investigating the individual contribution of different disc regions and we aimed at separating the effect of accretion from the effect of gravitational torques. The emerging picture is, in general, more complex than the often adopted simple assumption that resonant torques arising at the disc inner edge transfer angular momentum from the BHB to the disc, causing the orbital decay (e.g. Papaloizou & Pringle 1977; Goldreich & Tremaine 1980; Ivanov et al. 1999; Armitage & Natarajan 2002; Haiman et al. 2009).

Firstly, as already pointed out by a number of authors (e.g. MM08; Cuadra et al. 2009; Roedig et al. 2011; S11), the presence of a BHB clearing a cavity in the disc does not prevent gas inflows and eventual accretion onto the two BHs. Such inflowing gas has a non negligible role in defining the BHB-disc mutual evolution; it forms streams that follow the BHs, eventually exerting a net negative torque in the cavity region. The inflows catch-up with the BHs at a super-Keplerian speed, bend in front of them, thus exerting a net positive torque inside the binary corotation radius (i.e. at r < a). We therefore conclude that the origin of the dominating gravitational torques on the binary is purely kinematic, and can be understood without directly invoking resonances of the forcing binary potential with the disc.

The strong positive gravitational torque is related to the gas inflows only, and not to the formation of minidiscs around the two BHs nor to the eventual accretion. This becomes clear by inspecting Figs. 6, 10 and Table 2; although there is a factor up to ten difference in the mass bound to the BHs in minidiscs, and a factor of three difference in the accretion rate among the simulations, the positive torque is almost the same. This is because the gravitational torque is only related to the inflowing mass bending in front of the BHs, which is of the same order in all simulations. Note that different thermodynamics lead to different negative/positive torque balance. Although this might be an artifact of the sudden change of EoS in the iso simulations, such a balance has to be better understood to assess the fate of the binary.

Torques exerted by gas in the main body of the disc (excluding the disc edge) are instead negligible for the long term evolution of the system, since they average to zero (Fig. 6). Local torque maxima can be associated to the 3:1 and 4:1 OLRs (Artymowicz & Lubow 1994). Given the nature of the forcing potential, a quadrupolar torque pattern naturally develops (Figs. 8 and 9) causing a time oscillating torque (Fig. 7) with characteristic frequencies related to the beat between twice the binary frequency and the characteristic orbital frequencies of inhomogeneities developing in the disc in form of inner edge overdensities and spiral arms (Fig. 11).

Accretion plays an important role in the binary angular momentum budget. The binary can actually gain angular momentum even if it shrinks (see also S11) and its eccentricity increases. This is because accretion deposits a significant amount of angular momentum in the system by increasing its mass and especially its mass ratio (see Eq. (5)). The peculiar case of iso05 is quite interesting: here the total gravitational torque exerted by the gas is almost null (Fig. 6). However, the BHB-disc mutual torques are responsible for the eccentricity growth, meaning that, if accretion were neglected, the BHB would be forced to expand, as shown in the upper-right panel of Fig. 3. In this particular case, accretion is responsible for the shrinkage. Although whether this happens in Nature is questionable, it highlights the complexity of the BHB-disc interplay.

We should, nonetheless, be careful in drawing any conclusion about accretion from our simulations. If we scale our system to the fiducial binary of Roedig et al. (2011) and set M₁ = 2.6 × 10⁶   M_⊙ and a₀ = 0.039   pc, then the accretion rate we find corresponds to about 20–40 Ṁ_Edd for the secondary BH, and 4–7 Ṁ_Edd for the primary BH. In such situation, it is likely that most of the gas which is numerically accreted in our simulations will be expelled through outflows and winds. In this case, if the gas binds to the BHs before being expelled, it transfers its linear momentum to the holes, even if not accreted (Nixon et al. 2011). However, its mass does not add-up to the binary, making the question of angular momentum transfer more delicate and worthy of deeper and more refined investigations.

We numerically investigated different EoSs and sink radii r_sink. The numerical size of r_sink has a minor impact on the general behaviour of the system, affecting the mass bound to the BH in minidiscs and the accretion rate in the adia case. Conversely, the adopted EoS has a major impact on the results. In the iso runs, tenuous streams leak from an almost circular disc edge and fuel two well defined mini-discs, whereas streaming activity is more violent in the adia case, with thick spirals leaking from a deeply distorted disc edge, forming a diluted, ∞-shaped cloud around the binary. This difference in streaming activity affects both the overall structure of the outer disc and the long-term binary evolution (as mentioned above). Although the extrapolation of the results of all our four simulations would imply a fast BHB coalescence (in less than 10⁷ yr, for the fiducial binary parameters considered above), the almost perfect cancellation of the net gravitational torque in the iso simulations suggest that more realistic models for gas thermodynamics will be necessary to make clear-cut statements on this issue.

Acknowledgments

We thank the anonymous referee for suggesting necessary clarifying comments. We are indebted to Monica Colpi for lively and inspiring discussions. C.R. is grateful to Nico Budewitz for HPC support. The computations were performed on the datura cluster of the AEI. J.C. acknowledges support from FONDAP (15010003), FONDECYT (11100240), Basal (PFB0609), and VRI-PUC (Inicio 16/2010). M.D. and J.C. appreciate the warm hospitality at the AEI. This work has, in part, been supported by the Transregio 7 “Gravitational Wave Astronomy” financed by the Deutsche Forschungsgemeinschaft DFG (German Research Foundation).

References

All Tables

All Figures

	Fig. 1 Relaxed meridional density maps of the disc for the adia (top) and iso (bottom) runs; the black dots indicate the BHs, the axis are in units of the binary semi-major axis a. Note that this is not an azimuthal average but a vertical slice of the disc.
In the text

	Fig. 2 Time and azimuth-averaged values of the disc scale-height H / r and the disc eccentricity edisc as a function of radius for the different simulations.
In the text

Fig. 3 Semimajor axis (top panel in each plot) and eccentricity (bottom panel in each plot) evolution for all the four runs. In each panel, the red line represents the full evolution directly taken from the simulation whereas the black line is the equivalent evolution when considering the energy and angular momentum exchanges due to accreted particles only (i.e., basically due the evolution of the mass and mass ratio of the binary).

In the text

Fig. 4 Consistency check of the evolution of the binary angular momentum components for all the simulations. In each plot, the Lx, Ly and Lz component evolution is shown from the top to the bottom. In each panel, the dashed line is the ΔL exchange balancing the gravity torques at each time step, the dotted line is the ΔL exchange computed from the accreted particles, and the solid black line is the sum of those two. For comparison, the red line is the angular momentum evolution taken directly from the raw simulation data. All ΔL are in units of  [M0a0Vc,0] .

In the text

	Fig. 5 Decomposition of the Lz evolution for the four runs. The ΔLz contributions stem from a (dotted), e (dashed), M (dot dashed) and μ (long dashed) variations separately. The sum of the four contributions is given by the solid black line, whereas the solid red line is, as in Fig. 4, the angular momentum evolution taken directly from the raw simulation data. All ΔL are in units of  [M0a0Vc,0] .
In the text

Fig. 6 Differential torques dT / dr in units of averaged over the entire simulations. In each panel we show the differential torque on the primary (green), on the secondary (red), the sum of the two (blue), and the total integrated torque (black) according to Eq. (7). This latter is integrated starting from a inwards and outwards. Notably, the inward torque is positive, whereas the outward torque is negative.

In the text

Fig. 7 Each plot depicts the torques exerted by the disc on the BHB in the time domain (left panels) and in the frequency domain (right panels). In each plot, from the top to the bottom, pairs of panels refer to the five domains discussed in the text: total torque (i), torque exerted by the material located at r < a (ii), a < r < 1.8a (iii), 1.8a < r < 2.5a (iv) and r > 2.5a (v). In each left panel, the red line is the raw oscillating torque, and the blue one is the torque smoothed over three periods to show the average behaviour. In the associated right panel, we plot the power spectrum as a function of frequency in units of the initial binary orbital frequency. All torques strongly oscillate, showing several characteristic frequencies, however, note the different scales of the panels.

In the text

	Fig. 8 Top panel: typical instantaneous logarithmic surface mass density of the iso05 simulation. Bottom row: combined surface torque density exerted by the disc onto the binary (left panel), onto the primary BH (central panel) and onto the secondary BH (right panel). All plots are in code units: G = M0 = a0 = 1.
In the text

	Fig. 9 Same as Fig. 8 but for the adia05 run. Here the features highlighted above are even more evident (especially the symmetric overdensities at the inner edge of the disc).
In the text

	Fig. 10 Disc properties: average surface density (top panel), and Toomre parameter Q (bottom panel) as a function of r. Line style denotes iso10 (long dashed), iso05 (dot dashed), adia10 (solid) and adia05 (dashed).
In the text

	Fig. 11 Simple test model for the periodicity generated in the outside disc, i.e. at r > a. In the top panel we show the temporal evolution of the torque, whereas in the bottom panel we show the Fourier transform. In each panel, the blue curve is obtained by placing three point masses at distance 1.6a, 2.1a, 3.5a from a circular binary, and the one red is the torque found in the iso10 simulation.
In the text

Fig. 12 Accretion onto the two BHs. Top box: accretion rate as a function of time. The red line is the total accretion rate while the blue and black lines are the accretion rates onto the secondary and the primary respectively. Lower box: Fourier transform of the total accretion rate (lower panel) and of the torques (upper panel). In the upper panel, the green line is the total torque (highest peak normalized to 1) and the red line is the torque exerted by the material at r < a (multiplied by ten).

In the text