Migration of massive black hole binaries in self--gravitating accretion discs: Retrograde versus prograde

We study the interplay between mass transfer, accretion and gravitational torques onto a black hole binary migrating in a self-gravitating, retrograde circumbinary disc. A direct comparison with an identical prograde disc shows that: (i) because of the absence of resonances, the cavity size is a factor a(1+e) smaller for retrograde discs; (ii) nonetheless the shrinkage of a circular binary semi--major axis, a, is identical in both cases; (iii) a circular binary in a retrograde disc remains circular while eccentric binaries grow more eccentric. For non-circular binaries, we measure the orbital decay rates and the eccentricity growth rates to be exponential as long as the binary orbits in the plane of its disc. Additionally, for these co-planar systems, we find that interaction (~ non--zero torque) stems only from the cavity edge plus a(1+e) in the disc, i.e. for dynamical purposes, the disc can be treated as a annulus of small radial extent. We find that simple 'dust' models in which the binary- disc interaction is purely gravitational can account for all main numerical results, both for prograde and retrograde discs. Furthermore, we discuss the possibility of an instability occurring for highly eccentric binaries causing it to leave the disc plane, secularly tilt and converge to a prograde system. Our results suggest that there are two stable configurations for binaries in self-gravitating discs: the special circular retrograde case and an eccentric (e~ 0.6) prograde configuration as a stable attractor.


INTRODUCTION
Supermassive black hole binaries (BHBs) 1 have been predicted to be ubiquitous intermediate products of comparable mass galaxy mergers (e.g. Begelman et al. 1980). We consider the possibility of a BHB that has migrated to reach parsec separation and is now evolving in a retrograde selfgravitating disc. Without making any a-priori statements about the astrophysical likelihood of this scenario, we compare the migration in a retrograde disc to the prograde case. Following the seminal work of Goldreich & Tremaine (1980) in the context of planetary dynamics, several studies of sub-parsec BHBs in gaseous discs (Artymowicz & Lubow 1994;Syer & Clarke 1995;Artymowicz & Lubow 1996;Ivanov et al. 1999;Gould & Rix 2000;Armitage & Natarajan 2002, 2005Hayasaki et al. 2007Hayasaki et al. , 2008MacFadyen & Milosavljević 2008;Hayasaki 2009;Cuadra et al. 2009;Haiman et al. 2009;Lodato et al. 2009;Lin & Papaloizou 2011;Nixon et al. 2011;Roedig et al. 2011;Shi et al. 2011) employing various techniques, have established the picture sketched below. Initially the binary migrates inwards following the 'fluid elements' of the discs (type I) migration. When the local mass of the disc (see, e.g., Kocsis et al. 2012, for a thoughtful discussion) is of the order of the mass of the lighter hole, torques exerted by the binary are strong enough to evacuate the central region of the disc carving a central cavity. Accretion continues onto both holes through the cavity. In the prograde case (studied by the aforementioned papers), the cavity itself has a radius which is about twice the binary separation, and is held up by outer Lindblad resonances. The gas at the inner edge of the circumbinary disc leaks through the cavity streaming towards either holes. The streams bend in front of each hole and are partially gravitationally captured forming minidiscs. However, a large amount of gas is flung back out impacting the disc edge, efficiently mediating the disc-binary energy and angular momentum exchange (MacFadyen & Milosavljević 2008;Shi et al. 2011). The semi-major axis, a, shrinks and the eccentricity, e, grows and saturates around e∼ 0.6 (Roedig et al. 2011).
Recently Nixon et al. (2011) investigated the dynamical interaction of BHB embedded in light, non-selfgravitating circumbinary retrograde discs. The absence of outer Lindblad resonances (OLR) allows the circumbinary disc to survive unperturbed outside the orbit of the binary, with the inner edge impacting directly onto the secondary BH. So, the retrograde cavity is smaller and there is neither room to accommodate large minidiscs nor will there be stream bending nor will there be sling-shots. The gas at the inner edge is either accreted directly by the secondary hole, or loses its energy and angular momentum in the impact being funneled onto the primary. We study here the evolution of BHB embedded in massive, self gravitating retrograde discs. The first goal of this paper is to compare the evolution of the orbital elements of an initially circular BHBs migrating inside pro-versus retrograde discs. The pro-retro run comparison allows us to better understand the physics of the BHB-disc interaction. We demonstrate that, in both configurations, the evolution of the binary elements is well reproduced if the BHB-disc interplay is described in terms of purely gravitational interactions, in which the gas particles are effectively treated as dust. In this picture, the hydrodynamical properties of the disc is only relevant in determining the pace at which gas is supplied for direct interaction with the BHB. At very high eccentricity we find hints of an instability that leads to the binary tilting out of the rotation plane of the disc. We pursue this possibility with a dedicated higher resolution simulation and point out the most likely cause of this behavior. Combining these finding with those of Roedig et al. (2011), we discuss that BHBs in selfgravitating retrograde discs have three possible end-stages: (1) a stable circular motion coupled with a orbit decay similar to the aligned state; depending on the stochastic occurrence of the instability either (2) a fast inspiral if the exponential eccentricity growth and semi-major axis decay bring the BHB into the gravitational wave regime or (3) realignment of the BHB with its disc and evolution towards the eccentric, prograde standard scenario.
The paper is organized as follows. In § 2, we describe the numerical method, initial conditions and parameters use in our simulations. In § 3, we provide a compilation of our numerical findings, that we try to explain coherently with a simple model presented in §4. The aforementioned tilting instability is studied in § 5. Concluding, we discuss and summarize the overall findings in § 6.

METHOD
The model and numerical setup of this work is closely related to that of Roedig et al. (2011Roedig et al. ( , 2012(the latter hereafter pa-perI) and Cuadra et al. (2009), hence we only outline the key aspects in the following, and refer the reader to these papers for further details.
We simulate a self-gravitating gaseous disc of mass MDisc around two BHs of combined mass MBHB = M1+M2, mass ratio q = M2/M1 ∼ 1/3, eccentricity e and semi-major axis a, using the SPH-code Gadget-2 (Springel 2005) in a Table 1. Run labels and initial parameters in code units. The prograde run is taken from paperI Roedig et al. (2012). On a side note: N is the resolution, * * denotes that additional crosschecks were preformed to ensure that all conservation laws were numerically satisfied.  Cuadra et al. 2006). The disc, which is counter-rotating with the BHs, radially extends to about 7a, contains a circumbinary cavity of radial size ∼ > a and is numerically resolved by at least 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). The gas in the disc is allowed to cool on a time scale proportional to the local dynamical time of the disc t dyn = 2π/Ω0, where Ω0 = (GM0/a 3 0 ) 1/2 is the initial orbital frequency of the binary 2 . To prevent it from fragmenting, we force the gas to cool slowly, setting β = t cool /t dyn = 10 (Gammie 2001;Rice et al. 2005). This choice of β is discussed in previous work (e.g. paperI). All figures are shown in code units G = M0 = 1.

Discussion of initial conditions
We start the simulations using two different approaches for obtaining initial conditions (ICs), the relevant quantities of which are listed in Table 1. In the standard case, labeled retro throughout, a relaxed prograde disc is taken and the orbital angular momentum vector of the BHB, LBHB, is flipped and the eccentricity is changed. To be specific, we took the ICs of model adia05 in paperI, in order to facilitate a direct prograde vs retrograde comparison. This prograde run is labeled pro e = 0. Additionally, to measure how fast the memory of any out-of-equilibrium IC is lost and in order to avoid any spurious systematics, the following test was performed: the cavity was filled with gas to form a solid accretion disc lacking any low density region around the BHB. A circular BHB was placed around the center of mass (CoM) and integrated to measure how many Figure 1. Retrograde runs, average surface density Σ(r) -Top: "retro runs. Bottom: constant (1 + e)a runs "count ′′ runs, Additionally, the circular prograde case is overlayed, in which the surface density peaks at a radius that is 1 a larger than in the retrograde case. Note that M Disc in the count runs is smaller than in the retro runs. See text for more details; quantities are in code units.
orbits the system would need to relax. We observed an initial transient marked by high accretion and BHB shrinkage that lasted 45 orbits. Thereafter the system converged to a perfectly relaxed behavior, the evolution of which in terms of time derivative of the semi-major axisȧ, the eccentricityė and mass,ṀBHB, could not be distinguished from the retro e = 0 setup. We thus only show data that exclude the first 45 orbits, to avoid any potential initial transient. Unless otherwise stated, for each setup we use the orbits 45 − 90 for the analysis.
The second set of ICs, labeled count, uses a disc that was left to relax onto a retrograde, circular BHB for 750 orbits, which is why the disc is lighter as it got depleted over time. After 750 orbits, we modify the binary semi-major axis and eccentricity, keeping the quantity a(1+e) approximately constant. This is physically motivated by the previous observation that the gap size in a prograde disc is, at zero order, determined by the apocenter of the BHB orbit 3 . We then let the system evolve for 90 orbits and, as in the retro case, we keep orbits 45 − 90 for the analysis.
In Fig. 1, we plot the surface densities of all runs (averaging over orbits 45 − 90). In the top panel we show the retro profiles as a function of radius in units of a. As e increases, the edge of the cavity is pushed outwards, and the profile of the disc inside the cavity becomes shallower. The lower panel shows the profiles of the four count setups in comparison with the prograde, circular setup, as a function 3 In reality, in the prograde case, we observed the gap size to be slightly larger than 2a(1+e) due to non-accreted gas being flung out in a sling-shot, impacting at the gap edge, shocking and pushing the edge further outwards. In the retrograde setup, we do not expect sling-shots to play any role.  of radius now in units of a(1+e). For the count runs Σ(r) are very similar and the radial displacement in Σmax between prograde and retrograde is one semi-major axis (i.e., the disc extends inwards because of lack of resonances). In Fig. 2, the logarithmic surface density inside 5a is shown for the two runs retro e = 0 and retro e = 08, with color code and scale identical to paperI, Fig. 9. We can appreciate here the much smaller and well defined cavity in the circular case compared to the eccentric one.
The question of the necessary resolution is addressed both by changing the total number of SPH particles and by changing the mass of each SPH particle (thus giving a different disc thickness, which must be sufficiently sampled) As can be seen by inspecting table 1, the retro runs use 2, the count runs use 3 and the retro * * runs use 4 Million particles (where we state the resolution after discard-ing relaxation). We measured differential torque-densities, dTG,z/dr (see §3.2), per unit mass in retro e = 0 versus count00, and found small, but no significant differences between the averaged torques for the circular runs. This assures us that the resolution is sufficient in the retro runs, that the results are not sensitive to the precise choice of disc mass and that the relaxation of the discs to the new BHB orbit is fast (the count00 run effectively discards the first 795 orbits of relaxation at which time it has a resolution of 3 Million). To be conservative, we use 4 Million particles and monitor closely all conservation laws for the more complex dynamical geometry discussed in the end of the paper.

Binary evolution: general features
We start with a qualitative description of the several outputs of the simulations. In Fig. 3, we show the measured evolution of the individual accretion ratesṀ1,Ṁ2, semi-major axis a (for clarity normalized to a(t=45)) and the eccentricity e for all the retro runs. The circular prograde case is shown in the top left panel for comparison. For the highest e, we add a panel showing the inclination i 4 in degrees; 180°implies perfect antialignment (an additional line shows the inclination i2.5 of only the material inside of r < 2.5a(1 + e)). In all cases ȧ < 0 (here and in the following * indicates time average) and its absolute value increases with eccentricity; i.e., more eccentric counterrotating binaries appear to shrink faster. ė is positive for e > 0, but it is almost exactly zero for the quasi-circular run (retro e = 0). This is in contrast to the circular prograde case, which features ė > 0 (top left panel). The accretion rate onto each hole is comparable in the retro e = 0 and it is smaller than the prograde case. With increasing e, the accretion onto the secondary hole becomes larger and highly periodic (modulated on the orbital period of the binary). The retro e = 0.8 case shows a different general behavior. The binary leaves the disc orbital plane, making accretion on the primary hole easier (the secondary hole does not 'screen' the primary anymore) and the eccentricity growth has a turnover at e=0.9. We will return on this latter run extensively in §5.

Gravitational and accretion torques
As elucidated in paperI the evolution of the binary elements is driven by the sum of the gravitational plus the accretion torques:LBHB = TG + Tacc 5 . The total gravitational torque exerted by the sum of the gaseous particles on each individual BH can be computed as where r is the position vector with respect to the system CoM. The accretion torque Tacc is measured from the conservation of linear momentum in the accretion process of each particle, i.e, Tacc = r k × mjvj (here r k is the position Throughout the paper we use boldface for vectors. vector of the accreting BH, and mj and vj are the mass and velocity of the accreted gas particle). In the following, all torques are averaged both in time and azimuth. In all retro runs, we use a reference frame in which the angular momentum of the disc, L disc , is initially aligned to the z axis and points in the positive z-direction, whereas LBHB points in the negative z-direction. Assuming that the binary and the disc are coplanar, the system is effective two-dimensional from the dynamical standpoint. We can therefore safely take L ≈ Lz and T ≈ Tz, and unless otherwise specified, we will always refer to the z components of these quantities 6 . As discussed in the previous section, the natural scale of the system is the apocenter distance of the BHB. We thus integrate the differential torques both in-and outwards of a(1 + e) as: In Fig. 4, we show the differential, gravitational torque dTG,z/dr on the primary BH (red), on the secondary BH (green), the sum of the two (blue), and the total integrated torque (black). In opposition to the prograde case, the figure clarifies that the interaction between the disc and the binary takes place mostly at the location of the secondary, supporting a direct 'BH-gas impact' model, as described by Nixon et al. (2011). The differential torque density, dTG,z/dr , of the runs pro e = 0 and retro e = 0 is shown in Fig. 5, where the leftmost panels show the total dTG,z/dr on the binary, the middle panel the torque density onto the primary BH and the right, onto the secondary BH.
Here, we can clearly appreciate the absence of a gap between a and 2a in the retrograde run. The lower left panel highlights the unbalanced positive torque (blue) at the location of the secondary hole, which is responsible for the binary evolution (Note, that the orbital angular momentum of the binary is negative in the retro runs). Differentiating LBHB with respect to the relevant binary elements, angular momentum conservation yields: where µ = M1M2/M is the binary reduced mass. Note that all quantities in Eq. (3) can be directly measured from the simulation outputs. We can therefore use it to get an empirical measurement on how the angular momentum exchange is redistributed among the relevant binary elements. We have cast Eq. (3) in terms of the shrinkage rate and show in Fig. 7 that the two major contributions are Figure 3. Evolution of the accretion rates, semi-major axis, eccentricity and inclination -All retro runs are shown, as indicated above each panel. We also show the prograde run in the top left panel. In each panel we plot the individual accretion ratesṀ 1 , M 2 (multiplied by 10 5 ), the semi-major axis a and the eccentricity e. In the retro e = 0.8 panel we also include the relative BHB-disc inclination i in degrees. Note the different scales on the y axis; quantities are in code units.

Comparison of circular runs: prograde vs retrograde phenomenology
As we chose the setup of the retro runs to be the same as those of the adia05 run of paperI, we can compare the evolution of a and e directly. In Fig. 6, we thus plot the circular runs retro e = 0 vs pro e = 0. The evolution of a is not distinguishable over the 45 orbits that we show; on the other hand, the eccentricity behavior is different, such thaṫ e|retro = 0 whereasė|pro > 0. Upon inspecting Fig. 3, we find that also the accretion rates for the prograde setup are larger and that while the individual rates are equal for the retrograde, the secondary BH accretes and varies more in the prograde. Considering now the very left panels of Fig. 4, focusing only on the gravitational inwards torque, Tin,z, we find them to be very similar, but opposite in sign with respect to LBHB,z. However, the integrated Tin,z is slightly smaller in the prograde and there are three clear differences: (i) the radial location of largest (differential) torque amplitude onto the individual holes are displaced by a factor a ; (ii) the largest peak of the (differential) torque onto the binary (blue dashed line) is around r ∼ 1.7a(1 + e), (coinciding with the first OLR) where the torques onto the individual holes are in phase, is absent in the retrograde case; (iii) the substantial outwards torques found in the prograde case are completely absent in the retrograde case; Note that, integrated across the whole r range, the total torques on the binary always have opposite sign with respect to the orientation of LBHB, causing a net transfer of angular momentum from the binary to the disc both in the prograde and in the retrograde case. Turning to the redistribution of the angular momentum transfer among the binary elements as detailed by equation (3), we show the individual terms from the perspective ofȧ in Fig. 7 and first compare the two leftmost (i.e. circular) runs 7 . We observe that (i) the contribution of 2e 1−e 2ė (filled pentagons) is negligible for the retrograde asė = 0 (ii) 2 Tacc,z/Lz (open squares) contributes negatively for the retrograde, whereas positively (and three times stronger) for the prograde. The sign simply reflects the sign of Lz.
(iii) the mass transfer factors 2μ/µ andṀ /M are larger by a factor ∼ 2 in the prograde case.
(iv) summing up all components gives a very similarȧ Taken together, these points show that for the retrograde case, the gravitational torque is indeed directly related to the binary semi-major axis evolution, given that the distribution of angular momentum among the other elements of the binary (mass, mass ratio and eccentricity) is negligible. . We show the differential torque on the primary (red), on the secondary (green), the sum of the two (dashed-blue), and the total torques integrated inwards and outwards (black). Note that the angular momentum of the BHB is negative for retrograde runs and positive for prograde runs. See text for details.

Eccentric runs
We turn now to the description of the phenomenology of eccentric runs. In the following, we provide some analytical fits to the few measured data points provided by our set of 5 simulations ({retro e = 0 e = 0.1 e = 0.3 e = 0.5 e = 0.8}). We note that sometimes either the circular or the tilting run (e = 0.8) had to be excluded, and we caution that the fitting functions have only few degree of freedom less than the number of sample points, which is why we do not pursue any extrapolation or other use of the numerical fit-coefficients.
In Fig. 4, we plot the evolution of the torque structure from circular to eccentric, and in Fig. 7 the contributions to the shrinkage rate,ȧ are given, also using the total integrated torque components. Clear trends with increasing e are: (i) with the exception of the run retro e = 0.8, the integrated inward torque grows with e as Tin,z(e) ≈ Tin,z(e = 0) + 10 · Tin,z(e = 0) e 2 ; (4) (ii) the integrated outward torque grows linearly with e , for e = 0, but is always negligible with respect to the inward torque; (iii) a positive bump appears in the outward torque around r ∼ 1.5a(1 + e), for e = 0; (iv) the amplitude of the differential torques, dTG,z/dr , grows with e; (v) the maximum radius outside of which dTG,z/dr ∼ 0 is approximately rmax ∼ < 2.5a(1 + e) (cf. dashed blue lines).
(vi) retro e = 0.8 shows a qualitatively different torque structure, namely there is a large differential torque outside of r = a(1 + e), there is a high and positive net torque resulting from the fact, that now pairs of two subsequent individual, differential torques (either red or green) alternate in sign, rather than single peaks alternating after every extreme (all other panels). We conjecture, that this reflects the BHB leaving the plane of the disc; an effect that will be discussed further in § 5.
Measuring the derivatives of the orbital elements as a function of e, we additionally find (cf . Figs 8 and 9): • the dependence ofṀ2/Ṁ1(e) is a concave parabola with a best fit maximum around e ∼ 0.43;

Accretion rates
In Fig. 3, we can spot clear trends in bothṀBHB anḋ M2/Ṁ1 with e. We explore this further in Fig. 9, where we show the accretion rates, both their sum and their ratio, as a function of eccentricity. The combined accretion rate increases with increasing e, in contrast to what was observed for prograde discs in paperI. The explanation is straightforward: for prograde runs the cavity is forced to expand with higher e due to resonances and slingshots whereas this is not the case for retrograde runs where instead the binary penetrates further and into the disc edge and thus accretes more. Note also that the accretion rate growth is proportional to the enhancement of the inwards torque in Fig. 4 (with the exception of the retro e = 0.8 case), establishing a clear link between the two quantities. The ratio of secondary to primary accretion rate reaches a maximum somewhere in the range 0.4 < e < 0.6 and declines for higher e. Conversely, in the prograde case (Roedig et al. 2011), the ratio decreases to aroundṀ2/Ṁ1 ∼ 1 only for very high e ∼ 0.8 and there is no turnover.

ANALYTICAL MODELING
We now show that the phenomenology described in the previous section can be interpreted in the light of a simple 'dust model' in which particles experience a pure gravitational interaction with the BHB .

Retrograde runs
The retrograde case has been already described by Nixon et al. (hereafter N11, 2011). Following the same argument, we assume that a pure kinetic interaction between the secondary BH and the gas occurs at each apoastron passage. Since in our case q ≈ 1/3, we relax the assumption q ≪ 1, but we keep the Ansatz that the inner boundary of the disc interacts with the secondary hole only. A close inspection of the disc-BHB dynamics shows that the gas that 'impacts' on the secondary hole is partly accreted and partly slowed down to the point that is easily captured by the primary. We therefore assume that when a mass ∆M interacts with the secondary, a fraction 0 < α < 1 is accreted by MBH 2 , and the remaining 1 − α is captured by MBH 1 . We can therefore impose conservation of linear momentum in the interaction to get (see N11 for the mathematical derivation) F(e, q, α), where and the 1+q factors (absent in N11) stem from the fact that we relaxed the condition q ≪ 1. If we now average over the  binary orbit, and we take the continuum limit, the above expressions result in two coupled differential equations foṙ a andė, as a function of the accretion rateṀ and the parameter α. A complete solution to the dynamics requires an additional knowledge ofṀ and α as a function of a and e. However Fig. 3 indicates that on the short period of the simulations those are basically constant. We also note that e enters in the RHS of the equations only through the combinations 1 − e and 1 + e, providing, over small timescales, only small corrections that we neglect. We can therefore set  (6) and (8), perform the integral via separation of variables and linearize the solution to get: The results are shown in table 2 for t = 2π × 45 (this corresponds to t = 90 in Fig. 3). A visual comparison with Fig.  3 indicates that, although there are differences, the overall quantitative agreement is quite good. Discarding the run retro e = 0.8 (the BHB leaves the plane and our analysis does not apply), we notice however that the model slightly underestimate both the binary shrinking and the eccentricity growth for moderate e (runs retro e = 0.3 , e = 0.5).
If we now approximate the torque exerted on the gas as T G = (dL/dt) G ≈ −Ṁ ∆v G a (v G refers to the gas velocity), and consider that the gas is partly accreted by the secondary (i.e, ∆v G ≈ −2vM 2 ) and partly slowed down to the point is accreted by M1 (i.e., ∆v G ≈ −vM 2 ), we obtain T G ≈Ṁ (in simulation units), as observed in Fig. 4. This also trivially explains the behavior of T (e) total : for eccentric binaries, the secondary penetrates deeper in the inner rim of the disc, ripping off much more gas, and the torque is simply proportional toṀ . Note that, compared to N11 the systems we consider are substantially different: N11 use 3D SPH with an isothermal equation of state and no selfgravity. A key difference, is that they find minidiscs around the individual BH, which we do not. This must be attributed to the fact, that our gas is allowed to cool proportional to the local dynamical time scale, whereas their simulations are globally isothermal. Nonetheless the BHB-disc interaction is described by pure gravitational model, which can be applied to a large variety of situations, independently on the particular structure of the disc.

Prograde run
We put forward a 'dust' model also for the prograde case. As a matter of fact, Lindblad resonances are purely gravitational effects reflecting the fact that no stable hierarchical triplets exist in Newtonian dynamics unless the ratio of the semi-major axis of the inner and the outer binary is larger than a certain value. In our case we can think of the BHB as the inner binary and a gas particle in the disc as the outer body orbiting the center of mass of the BHB. The transition between stable and unstable configurations corresponds to the location of the inner edge of the cavity. We can therefore suppose that outside the inner edge, the disc-BHB interaction is weak (as certified by the lower left panel in Fig. 4). At the inner edge, particles are ripped off and stream toward the binary suffering a slingshot. In a three body scattering, the average BHB-intruder energy exchange in a single encounter is given by (Quinlan 1996;Sesana et al. 2006) where E = −GM1M2/(2a) is the binary total energy and m is the mass of the intruder. C is a numerical coefficient of order 1.5, which is determined via extensive scattering experiments. Recasting equation (11) in terms of a and differentiating, yieldsȧ However, hereṁ is not the accretion rate onto the BHs, but the 'streaming rate' of gas into the cavity, flung back to the disc. We give this number in table 2 of paperI 8 , and in the adia05 case we getṁ ≈ 5 × 10 −5 . Again, the linear proxy to the solution of equation (12) is simply at = a0[1−2Cṁt/M ]. If we substitute all the relevant quantities and fix t = 2π×45, we get a45 ≈ 0.965a. This is about 30% larger than the observed shrink, which is acceptable for such a simple model. The origin of the discrepancy is likely due to dissipation not taken into account here. Part of the streaming gas is actually captured by the hole: in this case, for a particle of mass m the binary acquires a binding energy −GM m/(2Rcavity) rather than losing an energy given by equation 11 (i.e., GmµC/a). In our case, µ ≈ 0.2. so that the two terms are comparable. Given that ≈ 25% of the streaming material is accreted, this accounts for the discrepancy. This model also provides an intuitive explanation of the negative torque building up in the region a < r < 2a in the prograde case. The energy exchange given by equation (11) corresponds to a velocity increment of the gas of the order ∆v G ≈ C 2Gµ/(M a). Providing a torque on the gas T G = (dL/dt) G ≈ −Ṁ ∆v G a (we assume an impulsive interaction at r = a and v G ·r = 0). This gives, in simulation units a torque acting on the binary of the order TBHB = (dL/dt)BHB = −(dL/dt) G ≈ 2 × 10 −5 , as observed in the lower left panel of Fig. 4. In our simulations, therefore, it makes no significant difference to the shrinking rate whether the BHB is aligned or antialigned with its disc. However this appears to be a coincidence residing in the ratio of the ripped-off gas in the co-rotating case, versus the accreted gas in the retrograde case. Our purely dynamical model does not provide an explanation for this, which stems from the hydrodynamical properties of the gas flow. However, it can be used to quickly determine the evolution of a binary once the in-streaming and accretion rates are known.

Accretion rates
We turn now to the observed trends in the accretion rates highlighted in Fig. 9. Those can be understood by comparing the effect of the mass ratio q and e on the mean distance and relative velocity of the two BHs with respect to the gas at the disc edge. We expect the primary to start accreting much more efficiently as soon as its influence radius at apocenter is equal to the apocenter of the secondary, since in this case the secondary can no longer completely shield the primary from the infalling material. At apocenter, the distances of the BHs to the CoM are ri = a(1+e)/(1+Mi/Mj ). If indeed accretion is Bondi-like, then the influence radius is given by: ra ≡ GM/(v 2 ∞ +c 2 s,∞ ), where v∞ = vi +v edge and the speeds at apocenter are (13) If we impose a turnover in the relative accretion when r1(e)+ ra 1 (e) = r2(e), and we solve for e, we get: which for the average sound speed of our disc, cs ≈ 0.2, and a ≈ 1 gives We should mention that the dependence on sound speed, cs, is not very strong as long as the discs are supersonic. Our discs have Mach numbers between 5 and 10. Additionally, the fact that in the retrograde case for e ∼ > 0.8 the ratiȯ M2/Ṁ1 becomes smaller than unity can be attributed to the fact that the binary leaves the disc plane and starts tilting, so the gas can come closer to the primary without being captured by the secondary.

ALIGNMENT BETWEEN THE BHB AND THE DISC PLANE: THE TILTING INSTABILITY
When looking at the bottom right panel Fig. 3, right, we find that the inclination of the BHB tilts mildly away from exact antialignment (i = 180) in the retro e = 0.8 run. At the same time e saturates its growth and turns over to become again more circular. In order to study this effect in more detail, we setup the additional run "retro * * e = 0.9" with two times the resolution of the default runs. Additionally, we monitored the numerical drift of the total angular momentum vector L to be small and, averaged over an orbit, consistent with a numerical zero. Also, all conservation laws are monitored very closely, and we stop trusting the evolution of "retro * * e = 0.9" after orbit 600. As indicators of the onset of tilting, we identify the following criteria observed in the low-resolution runs: Figure 10. Retrograde high eccentricity run, tilting of the BHB -The "retro * * e = 0.9" run; from the top to the bottom, panels show the individual accretion ratesṀ 1 ,Ṁ 2 (×1000), the total inclination i in degrees, the inclination i 2.5 towards only the inner rings, semi-major axis a and the eccentricity e. See text for more details; quantities are in code units.
• the ratio of the individual accretion rates would drop belowṀ2/Ṁ1 < 1; • the outward torque would become significantly positive; • the individual differential z−torques onto the BHs would stop alternating sign between each peak and rather be pairwise in phase; • the angular momentum (modulus) of the relevant parts of the disc be larger than that of the BHB, i.e. LDisc > LBHB.
We define the relevant angular momentum of the disc to be Because outwards of this radius, the amplitude of the individual differential torque is essentially zero 9 . While it is true that in order to transport outwards the angular momentum, a radially larger disc is required, the dynamical influence onto the BHB, i.e. the strength and orientation of the torque is given only by the mass and distance of the matter that orbits at radii where the torque is non-zero. For e > 0.9 the angular momentum of the BHB drops below 0.5LDisc, which is the condition under which a counteraligned stable configuration is not possible anymore King et al. (2005). At this point, it should be mentioned, that the two CoM, the one of the BHB and the global CoM exhibit increasing amount of relative movement with e.

Physical reasons for the tilting
The retrograde setup differs in several physical properties from the prograde, such that the BHB grazes the disc without the "stream bending" within the cavity. The appearance of minidiscs is suppressed in the retrograde case, because firstly the cavity is too small to accommodate the discs and secondly because the gas inflows onto the BHs is much more radial than in the prograde case. So the gas is accreted mostly in the wake of the BHs in the retrograde case. Every direct BHB-disc interaction is thus much stronger than in the prograde case, but takes place in more closely confined space. As e increases (and with it LDisc/LBHB) this induces an increasingly large relative movement between the CoM of the BHB with respect to the global CoM. It should be mentioned that the retrograde disc remains circular even if the BHB has high e ∼ 0.8 (cf. Fig. 11, top panel). Thus, it can be understood that the BHB, which has increasingly small angular momentum, has less and less rotational support feeling the impacts into its disc more strongly with increasing e. The details for a BHB to become secularly unstable to this tilting instability will strongly be disc-model dependent, which is why we only list under which conditions we expect a system to go unstable: (i) a retrograde disc of i = 180°is an unstable equilibrium configuration if LDisc > 2LBHB as is the case here: with growing e the BHB is increasingly interacting with the retrograde disc, it is not off-hand to expect perturbations that are not perfectly in-plane; (ii) the corotating case is stable attractor, there should be a threshold value of perturbation that leads the BHB to leave the retro-configuration in favor of the prograde one; (iii) if the accretion disc is not perfectly cylindrically symmetric, but rather clumpy and of non-zero vertical extent, perturbations out-of-plane and out-of-synchronization 10 are well possible.
It is thus not surprising that simulations where discs are light and self-gravity is not considered (e.g, Nixon 2012;Nixon et al. 2013) would not find such an instability.

The stages of tilting
We now further explore the tilting itself, assuming that it is a physical effect. The time evolution of this run covers almost 3000 orbits (when taking into account the very significant shrinkage of a(t)). After an initial increase of e there follows circularisation which again is followed by an increase in e. At the same time the inclination i is decreasing then saturating and oscillating until it finally decreases. a experiences several phases of strong shrinkage 11 , whereas the accretion rate onto the primary exceeds that onto the secondary until very late times. This behavior of exchanging inclination and eccentricity is somewhat reminiscent of the 10 Out-of-synchronization means that the accretion events do not happen coincidentally at the two BHs 11 Note, that while energy of the binary shrinks drastically, this is not due to numerical integration error but to the cooling function imposed on the gas. We have verified that at any given moment in the simulation, the numerical energy errors are (about a factor 10) smaller than the prescribed energy loss of the system.  . However, it is not straightforward to apply these few body calculations for the BHB-disc scenario studied here: a gaseous disc is not a rigid body and can dissipate angular momentum; nor is the disc a well located point mass with clearly defined eccentricity and inclination in itself.
The stages that the BHB undergoes on its way towards alignment can be subdivided into the following which are supported by Fig. 10 and 11: (i) Fast circularisation: the BHB tilts out of plane, e quickly shrinks to e ∼ 0.2, while at the same time a shrinks almost by a factor of two.
(ii) Disc fallback I : while the BHB has shrunk very fast and circularized and tilted out of plane, the disc had disintegrated into many annuli of various inclinations due to conservation of angular momentum. Now, however, the inner part of the gas forms a flat inner disc again. Meanwhile, the gravity torques from the precessing gaseous annuli cause the BHB eccentricity to grow. It should be noted, that the CoM of the BHB and the (inner) disc have drifted apart significantly by now. The accretion rates are low and the semi-major axis remains almost stationary.
(iii) Inner asymmetric ring collapse: once the gas has reformed into one (inner) planar disc, the new discconfiguration has lost the formerly large BHB upkeeping the cavity size. So, at this point, it re-adjusts to the new BHB "size" and a hot, dense ring forms (cf. third panel Fig. 11). This ring has very low Toomre parameter, becomes extremely clumpy and the CoM of BHB and disc-ring drift even further apart. Note, that the cooling in the ring depends on the local orbital timescale. The fallback of gas is accompanied by large accretion rates and thus induced semi-major axis shrinkage. The BHB eccentricity remains very high.
(iv) Circularisation II : when finally the ring was able to cool enough for its inner part to reach the BHB, the BHB starts circularizing and tilting into the disc. The resulting meta-stable configuration has an inclination of about i ∼ 90°, is mildly eccentric, a stable semi-major axis and comparatively low accretion rates.
(v) Disc fallback II : the initial disc has left behind some outer annuli that now start interacting with the outer part of the (inner) disc. The inner disc and the outer ring are now aligning which is accompanied by a series of precessions of the outer ring. Meanwhile the BHB increases its e and stays in the same orientation while the disc tilts towards the BHB.
(vi) Outer ring collapse: Similar to stage (3), now the outer ring becomes unstable under its own gravity as it has lost too much rotational energy by aligning with the inner 12 The Kozai effect (Kozai 1962) is a well known, well studied phenomenon of a secular three-body interaction when a binary interacts with a perturber that is on an inclined orbit with respect to the binary plane. The perturber can then oscillate in its inclination, pericentre, longitude of ascending node and eccentricity (Takeda et al. 2008) in so-called Kozai-cycles. The perturber feels secularly torqued by the binary and starts precessing. Depending on the time-scale of this precession, the binary eccentricity can be excited. For four-body systems, i.e. a binary and two perturbers, the same effect can be found (e.g. Murray & Dermott 1999) where now the two perturbers can additionally exchange inclination and eccentricity. Figure 11. The tilting instability: logarithmic column density, face-on xy, edge-on xz, edge-on yz -Time evolution (top to bottom) of the high accuracy, tilting "retro * * e = 0.9" run. Panels show the successive stages of the binary first tilting out of plane and then subsequently aligning with the disc. Note the spread into various annuli that the disc undergoes in the middle panels.
disc. (Most of) the outer ring collapses onto the inner disc. The accretion rates are high, but a only mildly shrinks.
(vii) CoM alignment: Once the inner disc and the BHB are the only two major pieces reformed and the inclination is below 60°, the BHB and its new, much more compact disc start aligning their CoMs. This is accompanied by a circularisation of the BHB, as during the interaction the BHB pushes the disc towards one side by impacting into it.
(viii) Standard e-growth: when finally the BHB and the disc have re-united their CoMs (apart from a tenuous outer ring), the excitation of e similar to what was described in Roedig et al. (2011) starts occurring, while the inclination furthermore decreases.
At the final state of the simulation, the BHB has an inclination of i ∼ 30°, an eccentricity of about e ∼ 0.3 and the secondary accretion rate is exceeding the primary. During this evolution, the drift of the total angular momentum vector was about 0.8°, and the β−cooling has induced a loss of L of about 1.2%. However, we do not trust the evolution past orbit 600, because there appears noise in both the total energy and the total angular momentum and the ratio between r sink and a are severely weakening the inwards torque, thus vitiating the secular effects drastically. However, from a physical point of view, the aligned, prograde state should be energetically preferred and once the disc has aligned with the BHB, it will be forced to expand its cavity analogously to the prograde case. The numerical problem of performing several thousands of orbits is challenging and we argue here that studying the entire final alignment up to e ∼ 0.6 (cf. Roedig et al. (2011)) would be a waste of CPU time, since the physical arguments are quite strong if indeed the instability occurs in the first place. Also, it is our suspicion that different choices for β, cooling prescriptions, and global disc thermodynamics will alter the 8 stages identified here. The cooling, i.e. the loss of angular momentum seems to be the most crucial numerical parameter in setting the timescales for the "ring collapse" and the dissipation of the several annuli.

SUMMARY
We presented numerical work on the evolution of the orbital elements of a comparable mass black hole binary surrounded by a retrograde selfgravitating disc. First, it was found that for quasi-circular BHB-disc configurations, the difference between the retrograde and prograde scenario is chiefly in the eccentricity and mass ratio evolution. The eccentricity remains close to zero in the retrograde case whereas it grows in the prograde case. The semi-major axis shrinkage was found to be the same for the circular BHB irrespective of disc (anti)alignment. For increasingly eccentric set-ups, the growth of e was found to be exponential in time and the coupled evolution of the semi-major axisȧ was observed to be proportional to e(t), thus making the decay very efficient. Additionally, it was found that the individual accretion rates onto the two BHs grow, but that their ratio secondary over primaryṀ2/Ṁ1 assumes a maximum at eccentricities around e ∼ 0.45.
We demonstrated that all the observed features can be satisfactorily explained using a simple 'dust' model in which the BHB-disc interaction is purely gravitational. In the retrograde case, the absence of resonances together with the high BH-disc relative velocity (i.e., a small cross section) imply that the gas interacts with the binary only via direct impact onto the secondary. This is confirmed by the fact, that the turnover ofṀ2/Ṁ1(e) could be explained by Bondi-Hoyle like accretion assuming ballistic speeds of discedge and binary at apocenter. Generalizing the model of N11, we find that this simple 'impact model' is able to catch of the relevant features of the binary evolution at different eccentricities. In the prograde case, the evolution is determined by the gas inflows flung back to the disc. Treating the gas as a dust particle, we found that a simple model based on three body scattering theory fully accounts for the BHBdisc interplay. In those models, however, the mass accretion rate and instreaming rate stem from the hydrodynamical properties of the gas flow and are input parameters.
The retrograde setup shows non-negligible relative motion between the BHB center of mass and the global center of mass (CoM), even in circular runs. In the prograde case, this had been a very small effect only. There are indications that the BHB eccentricity will not grow to e ∼ 1 in the retrograde case, but that the impacts into the disc inner edge destabilize the BHB until it starts moving out of the binary-disc plane. The most important physical conditions for this instability to occur be that the angular momentum of the disc exceed twice that of the BHB and that the disc be sufficiently nonaxisymmetric for out-of-plane perturbations to occur. The details of the subsequent alignment are likely to be sensitive to the cooling prescription of the disc, but it is to be expected that the disc disintegrates into several annuli of various inclinations which precess within each other. It was observed that precession of the disintegrated disc annuli causes the BHB eccentricity to grow while the BHB CoM and the disc CoM do not overlap. There were stages wherein the multiple annuli aligned again to try form a coherent, flat disc. Due to loss of rotational support from the shrunk BHB, this was preceded by a short period of fragmentation-like behavior where the disc collapsed shortly into a dense ring, to spread out into spiral arms later. The BHB eccentricity exhibited few oscillations and it is expected that the final state will be a perfectly aligned disc-BHB system of eccentricity around e ∼ 0.6.
We note, that it is well conceivable that, depending on the physical separation and mass, during any point of the highly eccentric phases the BHB might enter the relativistic regime and rapidly lose energy via gravitational radiation. At this point the circularization effect of gravitational radiation would compete with the eccentricity growth enforced by the retrograde disc, possibly causing a distinctive eccentricity evolution of the system that would be possibly detectable and discernible by low-frequency gravitational wave interferometers (Amaro-Seoane et al. 2013;Consortium 2013), offering an interesting observational test of the BHB-counterrotating disc coupling.

ACKNOWLEDGMENTS
CR is indebted to Kareem Sorathia and Sarah Buchanan without whom this paper would not have been written. CR is grateful to Julian Krolik and Massimo Dotti for helpful sug-gestions and was happy to have informative conversations with Cole Miller, Sanchayeeta Borthakur and Jens Chluba. All computations were performed on the datura cluster of the AEI and CR thanks Nico Budewitz for his HPC support. This work was funded in part by the International Max-Planck Research School, NSF grants AST-0908336, AST-1028111 and the SFB Transregio 7.