Massive gaseous discs around SMBH binaries: Binary decay and tidal disruptions

We investigate the evolution of black hole binaries embedded within geometrically thin gas discs. Our results imply that such discs can produce black hole mergers for relatively low-mass binaries, and that a significant population of eccentric binaries might exist at separations of a few 0.01 pc. These binaries may be detectable due to the time-variable accretion on to the black holes. If the disc fragments, then the newly-born stars will continue driving the binary to its coalescence, although at a slower rate. Interestingly, our preliminary analysis shows that these stars will be disrupted at a rate of ∼10−4–2 · 10−5 events per


INTRODUCTION
Galaxies form via mergers of lower mass progenitors, most of which host supermassive black holes (SMBHs).Following a galaxy merger, the SMBHs sink toward the centre of the merger product and form a binary.What happens next is less clear: angular momentum loss to either stars or gas is required to bring the black holes into the gravitational radiation inspiral regime, in which the coalescence occurs rapidly.
The dissipative nature of gas allows inflow down to pc scales in galactic nuclei after galaxy mergers (Springel et al. 2005;Mayer et al. 2007).Prior work has established that nuclear gas discs are likely to be important for SMBH binary mergers ( Here we revisit the interaction between a SMBH binary and a relatively low mass circumbinary gas disc.Our main innovation is to use simulations that directly resolve the transport of angular momentum due to the disc self-gravity.We study the rate at which the gas drives the binary toward merger, the effect that the gas disc has on the eccentricity of the binary, the time-variable signal associated with accretion on to the black holes, and the effects that disc fragmentation could have in the evolution of the system.For further details refer to our full-length papers (Cuadra et al. 2009;Roedig et al. 2011; Amaro-Seoane, Brem & Cuadra, in prep.).

SIMULATIONS OF NON-FRAGMENTING DISCS
We present models that follow a binary with mass ratio q = M 1 /M 2 = 3. Results are expressed in terms of the initial binary separation a 0 , mass M = M 1 + M 2 , and orbital frequency Ω 0 = (GM/a 3 0 ) 1/2 .Around the binary we set an aligned, co-rotating circular disc of gas.The disc has a total mass of M d = 0.2M, extends initially from 2a 0 to 5a 0 , and has a thickness of H/R = 2M d /M. a e-mail: jcuadra@astro.puc.clWe evolve the system with a modified version of the SPH code G-2 (Springel 2005).The code solves for the hydrodynamics and Newtonian gravitational forces of the system.We take into account radiative cooling defining a cooling time proportional to the orbital time of the gas, t cool = β/Ω.For most simulations, we take β = 10.This value of β ensures that the disc transports angular momentum instead of fragmenting (Gammie 2001;Rice et al. 2005;Nayakshin et al. 2007;Alexander et al. 2008).In our latest study (Amaro-Seoane et al., in prep.), however, we follow the evolution of a disc that does fragment ( §3).

Orbital evolution
The evolution of the disc for our fiducial model (β = 10; Cuadra et al. 2009) is shown in Fig. 1.After a transient phase, the gas reaches a quasi-steady state characterised by a spiral pattern that transfers the angular momentum of the binary outward, making its orbit shrink.A thorough study of the nature and strength of the different torques involved in the process was presented by Roedig et al. (2012).
Figure 2 shows the evolution of the binary orbital elements from our simulations.The orbit of the binary shrinks initially at a rate da/dt ∼ −10 −4 a 0 Ω 0 , in good agreement with analytical estimates (Ivanov et al. 1999).Later the rate of decay slows down due to the disc expansion, which is a consequence of the finite disc extension in our numerical models.
We also measured the binary eccentricity, finding that it initially grows at a rate de/dt ∼ 10 −4 Ω 0 when starting from low e. Figure 3 shows the eccentricity evolution from a series of simulations that have different initial eccentricities.It is clear that this quantity tends to a value in the range e ∼ 0.6-0.8.We interpret this behaviour as the result of the secondary-disc gravitational interaction at apocentre.During that orbital phase the secondary drives an instantaneous overdensity in the inner part of disc.If the eccentricity is low, then the angular velocity of the disc is lower than that of the secondary and the overdensity decelerates the secondary, increasing its eccentricity.If This is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.the eccentricity is high, the opposite happens and the eccentricity goes down.Equilibrium is achieved when the angular velocities match, at a binary eccentricity of e ∼ 0.7, which is what we observe in the simulations.Even though this eccentricity would be mostly lost at later stages dominated by gravitational wave emission, we still expect a small residual eccentricity to be detectable by a LISAkind mission (Roedig et al. 2011).

Accretion and its variability
Gas that gets within 0.1a 0 of either black hole is taken away from the simulation to prevent the very short timesteps that would be necessary to follow it.This numerical trick is also useful to study the accretion onto each black hole (Cuadra et al. 2006), in particular its variability and behaviour as a function of the orbital phase.Figure 4 shows the accretion rate onto each black hole as a function of time.Overall, more gas is accreted by the secondary black hole, which is closer to the gaseous disc and has a larger specific angular momentum.Both curves show a large peak at the beginning, which is a result of the initial conditions -gas that was situated inside the tidal truncation radius.The accretion rates increase again at t ≈ 400Ω −1 0 , after the disc loses its ring-like shape and a lot of radial mixing occurs.
The accretion rates are strongly variable, especially on the time-scale of the orbital frequency.Moreover, as the binary gains eccentricity, the periodicity of its signal becomes more pronounced and additional preferred timescales show up: higher harmonics of the binary frequency, the orbital frequency of the densest part of the disc, and the difference (beat) frequency between both orbits (Roedig et al. 2011).Once a M ∼ 3 × 10 6 M binary reaches a separation of 0.01 pc, the orbital frequency gets to 05003-p.2  an interesting observable range ∼0.1 yr −1 .Observation of variability on these frequencies will thus help identify the presence of binaries in such close orbits, which otherwise are extremely difficult to detect.

DISC FRAGMENTATION
In our latest study we have concentrated on the evolution of discs that cool fast enough for its gas to fragment (Amaro-Seoane et al. in prep.).Numerically, this means that we choose a cooling time β ≤ 5.As a result, the disc now forms clumps, which grow in a runaway fashion.Modelling this with a pure SPH code is not feasible, as the growing densities require ever shorter time-steps.To circumvent this problem, we follow Nayakshin et al (2007) and introduce sink particles to model the proto-stars that we expect would form in these large density regions.In the models, these particles grow by mergers with other stars and by accreting the surrounding gas.We run the SPH simulations for several hundred binary dynamical times.During this time interval the gaseous disc fragments and the vast majority (≥90%) of the gas is turned into stars, as expected.The system then reaches a quasisteady state in which stars very slowly accrete the tenous left-over gas (see Nayakshin et al. 2007).
The fragmentation rate is set by the cooling time of the disc, thus discs with lower values of β will evolve faster.We can see this in Fig. 5, which shows the mass in stars as a function of time for all the simulations.Moreover, for otherwise equal runs, shorter cooling times result in larger amounts of stars, as expected.As the total stellar mass is approximately constant, the typical stellar masses will be lower for shorter cooling times.

Long-term N-body evolution and tidal disruption events
To continue our study of the evolution of the SMBHs and circumbinary disc system, we take the masses, positions, and velocities of all sink particles and use them as input 05003-p.3EPJ Web of Conferences in direct-summation N-body simulations.These runs are performed using a modified version of the N6 code (Aarseth 1999).
As the system evolves, some of the stars interact with the binary, typically taking away some of its orbital energy.Figure 6 shows the evolution of the binary orbit for different simulations.The shrinking is about two orders of magnitude slower than in the case of the gas-driven evolution, proceeding at a rate da/dt ∼ −10 −6 a 0 Ω 0 .
During the direct-summation N-body runs any star entering the tidal radius of one of the SMBHs is considered to be tidally disrupted and its mass is added to the mass of the hole.Our preliminary analysis shows interesting rates of ∼10 −4 -2 • 10 −5 events per year per system.

Figure 1 .
Figure 1.Logarithmic maps of the disc column density (in of Ma −2 0 ) at different times.The panel at t = 0 shows the smooth initial conditions.Until t ≈ 350Ω −1 0 , material piles up at R ≈ 3a 0 , forming a dense ring.The ring breaks due to its self-gravity, developing a spiral pattern.The disc stays in that state until at least t ≈ 1200Ω −1 0 , when the simulation is stopped.Figure from Cuadra (2011).

Figure 2 .
Figure 2. Evolution of the binary orbit elements a (left panel) and e (right) for different simulations.Solid lines show runs with the binary in an initially circular orbit, while dashed lines show binaries with e = 0.1 initially.In both cases the eccentricity grows, and the orbital decay proceeds at the same rate.Figures from Cuadra et al (2009).

Figure 3 .
Figure 3. Eccentricity evolution of binaries starting with different eccentricities.This quantity tends to a value in the range e ∼0.6-0.8. Figure from Roedig et al (2011).

Figure 4 .
Figure 4. Accretion rate onto each black hole as a function of time.Notice the variability and the fact that the accretion rate is higher for the secondary.Figures from Cuadra et al (2009).

Figure 5 .
Figure 5. Accumulated stellar mass formed in the disc for our fragmentation simulations, scaled to a binary of 3.5 × 10 6 M .The shorter the cooling time, the fastest the disc turns its mass into stars.Figure from Amaro-Seoane et al. in prep.

Figure 6 .
Figure 6.Evolution of the semi-major axis of the binary driven by its interaction with the stars formed in the circumbinary disc.