In currently favoured hierarchical cosmologies, the formation of massive black hole binaries (MBHBs) following galaxy mergers is unavoidable. Still, due the complex physics governing the (hydro)dynamics of the post-merger dense environment of stars and gas in galactic nuclei, the final fate of those MBHBs is still unclear. In gas-rich environments, it is plausible that turbulence and gravitational instabilities feed gas to the nucleus in the form of a series of cold incoherent clumps, thus providing a way to exchange energy and angular momentum between the MBHB and its surroundings. Within this context, we present a suite of smoothed-particle-hydrodynamical models to study the evolution of a sequence of near-radial turbulent gas clouds as they infall towards equal-mass, circular MBHBs. We focus on the dynamical response of the binary orbit to different levels of anisotropy of the incoherent accretion events. Compared to a model extrapolated from a set of individual cloud-MBHB interactions, we find that accretion increases considerably and the binary evolution is faster. This occurs because the continuous infall of clouds drags inwards circumbinary gas left behind by previous accretion events, thus promoting a more effective exchange of angular momentum between the MBHB and the gas. These results suggest that sub-parsec MBHBs efficiently evolve towards coalescence during the interaction with a sequence of individual gas pockets.


Clumpy accretion onto MBHBs]Accretion of clumpy cold gas onto massive black hole binaries: a possible fast route to binary coalescence Goicovic et al.]Felipe G. Goicovic1, Cristián Maureira-Fredes, Alberto Sesana,\newauthorPau Amaro-Seoane, and Jorge Cuadra
Heidelberg Institute for Theoretical Studies (HITS), Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mülenberg 1, 14476, Potsdam-Golm, Germany
School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom
Institute of Space Sciences (ICE, CSIC) & Institut d’Estudis Espacials de Catalunya (IEEC) at Campus UAB, Carrer de Can Magrans s/n 08193 Barcelona, Spain
Institute of Applied Mathematics, Academy of Mathematics and Systems Science, CAS, Beijing 100190, China
Kavli Institute for Astronomy and Astrophysics, Beijing 100871, China
Zentrum für Astronomie und Astrophysik, TU Berlin, Hardenbergstraße 36, 10623 Berlin, Germany
Max Planck Institute für extraterrestriche Physik (MPE), D-85748 Garching, Germany


[[ 2017

ccretion, accretion discs – black hole physics – hydrodynamics – galaxies: evolution – galaxies: nuclei

1 Introduction

Massive black holes (MBHs) with masses between (Kelly & Shen 2013) have been observed to inhabit the nuclei of massive galaxies. And even though these objects are tiny with respect to their hosts, they are believed to play a fundamental, yet not fully understood, role in galaxy evolution. In fact, the observed tight correlations between the MBH masses and the properties of their galaxies (e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000) strongly suggest a coevolution scenario through feedback processes (Kormendy & Ho 2013, and references therein).

Furthermore, according to the currently favoured CDM cosmology, galaxies assemble hierarchically through mergers following their parent dark matter haloes, building up from smaller structures (White & Rees 1978). A natural outcome of this process is the formation of massive black hole binaries (hereafter MBHBs). Following a merger of two galaxies, we expect both MBHs harboured in their nuclei to sink toward the centre of the merger remnant and eventually form a binary (Begelman et al. 1980). Therefore, understanding MBHB formation and dynamical evolution is a necessary piece for reconstructing the puzzle of the hierarchical growth of structures in the Universe. This is particularly interesting if we consider that MBHBs are powerful sources of gravitational waves (GWs), detectable by the future space-based laser interferometer space antenna (LISA, Consortium et al. 2013; Amaro-Seoane et al. 2017). For GWs to carry away enough energy, however, the two MBHs have to reach a separation of pc. This is extremely small compared to the pc separation at which dynamical friction, initially driving the evolution of the two MBHs in the merger remnant, ceases to be efficient (Binney & Tremaine 2008).

The exact evolution and the final fate of the MBHs strongly depends on the content and distribution of both gas and stars in the galactic nucleus, and different processes can change dramatically the merging timescale, even resulting in values higher than the Hubble time. Although there is some indirect evidence that these binaries coalesce into a single object (e.g., the low scatter in the relation, Jahnke & Macciò 2011), the processes involved are not completely understood.

For instance, due to its dissipative nature, gas can be very efficient in absorbing and transporting outwards the angular momentum of the pair, likely leading to a rapid evolution and eventual coalescence. From observations, as well as numerical simulations, it has been established that in gas-rich galaxy mergers there is a large inflow of gaseous material to the central kiloparsec of the galactic remnant, often resulting in a massive circumnuclear disc (Sanders & Mirabel 1996; Mayer et al. 2007). Driven by dynamical friction and global torques from this disc, the pair of MBHs decays very efficiently down to separations of the order of 1–0.1 pc, where it forms a gravitationally bound binary (Escala et al. 2004, 2005; Mayer et al. 2007; Fiacconi et al. 2013; Roškar et al. 2015; del Valle et al. 2015). At these sub-parsec scales, most theoretical and numerical studies have focused on the evolution of binaries surrounded by a gaseous circumbinary disc, often either co-rotating (see e.g. Ivanov et al. 1999; Armitage & Natarajan 2005; Cuadra et al. 2009; Haiman et al. 2009; Lodato et al. 2009; Nixon et al. 2011; Roedig et al. 2011, 2012; Kocsis et al. 2012; Amaro-Seoane et al. 2013; D’Orazio et al. 2013; Tang et al. 2017) or conuter-rotating (see e.g. Nixon & Lubow 2015; Roedig & Sesana 2014) with respect to the binary’s orbital motion. These discs are generally assumed to be well-defined, smooth, and relaxed, with no attempt to link their presence to the gaseous environment around the binary, nor to the fuelling mechanisms that bring gas to the nucleus. Furthermore, all these idealised scenarios are subject to the disc consumption problem, namely, if the disc dissolves through some process (e.g. star formation, AGN/supernovae feedback, Lupi et al. 2015), the evolution of the binary orbit stops. In fact, the evolution of MBHBs in gas-rich environments is intimately related to the unsolved problem of gas supply to the centre of galactic nuclei.

The fuelling of gaseous material onto galactic nuclei is currently a critical, yet uncertain piece in the galaxy formation puzzle. The wide range of physical scales involved makes the fuelling a very complex process – the gas needs to lose its angular momentum efficiently to be transported from galactic scales down to the nuclear region (Hopkins & Quataert 2011). A plausible mechanism that could break this angular momentum barrier is provided by chaotic feeding or ballistic accretion of cold streams when the gas is cold enough to form dense clouds or filaments (Hobbs et al. 2011).

Cold accretion onto black holes has long been predicted by both theory and simulations. The possible relevance of discrete and randomly oriented accretion events was first highlighted by King & Pringle (2006) (see also King & Pringle 2007; Nayakshin & King 2007; Nayakshin et al. 2012), which they referred to as \enquotechaotic accretion scenario. They argue that a series of gas pockets falling from uncorrelated directions could explain the rapid growth of MBHs across cosmic time, and can reproduce some of the observed properties (e.g. the relation). One of the most extensive efforts to study the evolution of multiphase gaseous haloes and its impact on MBH accretion is the work of Gaspari and collaborators (Gaspari et al. 2012; Gaspari et al. 2013, 2015; Gaspari & Sa̧dowski 2017; Gaspari et al. 2017a, b). By means of numerical simulations, they have shown that realistic turbulence, cooling, and heating affecting the hot halo, can dramatically change the accretion flow onto black holes, departing from the idealised picture of the Bondi prescription (Bondi 1952). Under certain physical conditions, cold clouds and filaments condense out of the hot phase due to thermal instabilities. Chaotic collisions then promote the funnelling of the cold phase towards the MBH, leading to episodic spikes in the accretion rate.

A likely example of this kind of discrete accretion events is the putative molecular cloud that resulted in the unusual distribution of stars orbiting our Galaxy’s MBH. Several numerical studies have shown that portions of a near-radial gas infall can be captured by a MBH to form one or more eccentric discs that eventually fragment to form stars (Bonnell & Rice 2008; Hobbs & Nayakshin 2009; Alig et al. 2011; Gaburov et al. 2012; Mapelli et al. 2012; Lucas et al. 2013; Mapelli & Trani 2016), roughly reproducing the observed stellar distribution (Paumard et al. 2006; Lu et al. 2009). An interesting feature of these observed stellar discs is that they appear to be misaligned with respect to the Galactic plane, indicating that the infalling gas had an angular momentum direction unrelated to that of the large scale structure of the Galaxy.

From an observational perspective, the development of multi-wavelength observations has started to unveil the multiphase structure of massive galaxies. These galaxies have been observed to harbour neutral and molecular gas down to 10 K (e.g. Combes et al. 2007). The Atacama Large Millimetre Array (ALMA) has opened the gate to high-resolution detection of molecular gas in early-type galaxies. For instance, there have been observations of giant molecular associations in the NGC5044 group within 4 kpc, inferred to have chaotic dynamics (David et al. 2014). In massive galaxy clusters, ALMA has detected molecular hydrogen in the core of the Abell 1835 galaxy cluster (McNamara et al. 2014), which may be supported in a rotating, turbulent disc oriented nearly face-on. In the Abell 1664 cluster, the molecular gas shows asymmetric velocity structure, with two gas clumps flowing into the nucleus (Russell et al. 2014). More importantly, in a recent study, Tremblay et al. (2016) reported observations that reveal a cold, clumpy accretion flow towards a massive black hole fuel reservoir in the nucleus of the Abell 2597 galaxy cluster, a nearby giant elliptical galaxy surrounded by a dense halo of hot plasma. They infer that these cold molecular clouds are within the innermost hundred parsecs, moving inwards at about 300 kilometres per second, likely fuelling the black hole with gaseous material. The authors claim that this is the first time a distribution of molecular clouds has been unambiguously linked with MBH growth.

Despite these several clues suggesting that stochastic condensation of cold gas and its accretion onto the central MBH is essential for active galactic nuclei, it has been barely explored as a possible source of gas for sub-parsec MBHBs. Notable exceptions are the studies of Dunhill et al. (2014), Goicovic et al. (2016, hereafter Paper I) and Goicovic et al. (2017, hereafter Paper II), where the Authors have numerically modelled the interaction of a single gaseous cloud with MBHBs for a variety of orbital configurations. In particular, Paper II extrapolated the results for these single cloud models to estimate the long-term evolution of a binary interacting with a sequence of near-radial infalling clouds, finding a very efficient shrinking of its orbit. However, the approach presented in Paper II is subject to a series of limitations, mainly, considering only the prompt evolution during the first impact of the clouds without the secular effects of the remaining gas. These effects can become very important after a series of accretion events due to the formation of circumbinary structures that can continue extracting angular momentum and increasing accretion onto the binary. With the goal of overcoming these limitations, we present a suite of hydrodynamical simulations of MBHBs interacting with a sequence of clouds, infalling incoherently onto the binary with different levels of anisotropy. These models now combine the ‘discrete’ evolution due to the prompt interaction of each individual cloud, with the secular torques exerted by the left-over gas.

This paper is organised as follows. We first describe our numerical model in Section 2, highlighting the main differences with respect to the studies presented in Paper I and Paper II. In Section 3 we show the dynamical evolution of the binary in response to the different levels of anisotropy of the accretion events, and by comparing with the single cloud extrapolations, we identify the effects of the non-accreted material. In Section 4 we study the same effects on a different suite of runs featuring an higher cloud supply rate. We finally summarise and discuss the implications of our results in Section 5. Unless stated otherwise, throughout this paper we use units such that , where is the gravitational constant, and , are the initial binary separation and mass, respectively.

2 The numerical model

We base our model of multiple clouds on the simulations presented in Paper I and Paper II. The new suite of simulations is fully described in the companion paper Maureira-Fredes et al. (2018). We model the interaction between the gas clouds and the MBHBs using the smoothed-particle-hydrodynamics (SPH) code gadget-3, an updated version of the public code gadget-2 (Springel 2005).

As described in Paper I, we modified the standard version of the code to include a deterministic accretion recipe, where the sink particles representing the MBHs accrete all bound SPH particles within a fixed radius () as implemented by Cuadra et al. (2006); we set this radius to . Additionally, we treat the thermodynamics of the gas using a barotropic equation of state (i.e. pressure is a function of density only) with two regimes; for densities lower than some critical value (), the gas evolves isothermally, while for higher densities the gas is adiabatic. This prescription captures the expected temperature dependence with density of the gas without an explicit implementation of cooling, and has the effect of halting the collapse of the densest gas, avoiding excessively small time-steps that would stall the simulations. For the simulations presented in this paper we have chosen , which is two orders of magnitude smaller than the value presented in Paper I and Paper II. This change is necessary to avoid stalling the simulations due to the formation of a large amount of clumps in some configurations. Since we are not following in detail the evolution of these clumps, stopping their collapse in a lower density should not affect greatly our results. This fact only becomes important when scaling our simulations to physical units – for the fiducial model presented in Paper I, where the binary has and pc, the value of the critical density is now g cm, which is not necessarily an accurate representation of the collapse of molecular hydrogen (Masunaga & Inutsuka 2000).

To improve the orbit integration accuracy we have made changes with respect to the single cloud models presented in Paper II. First, we have modified the code to remove the sink particles from the gravity tree approximation, and their forces are instead summed directly. In simple words, each sink particle sees all other particles individually. Additionally, the orbit of the black holes is integrated using a fixed time-step, set equal to .

A more critical modification with respect to the previous simulations is the frequency at which the code makes the domain decomposition and full tree reconstruction. The single cloud models of Paper II had sufficient resolution to study the transient evolution of the binary during the first impact of the gas. Conversely, the typically much smaller secular torques exerted by the remaining gas were not resolved. The main source of this numerical noise is the aforementioned frequency, more specifically, when the tree reconstruction is not done at the same time with the black hole accretion, the sudden disappearance of particles during an update creates spurious jumps in the total angular momentum, which are generally larger than the small torques that we are trying to resolve with these models. In order to improve the conservation of angular momentum, we modified gadget to enforce the domain decomposition and full tree reconstruction every time the sink particles are active. This implementation is complementary to the standard criterion based on the amount of active particles (Springel 2005). While this greatly increases the computational cost (typically by a factor of 10), the numerical fluctuations in the total angular momentum are vastly reduced, allowing a robust study of its dynamical evolution during all phases of the interaction (see below).

Considering the aforementioned modifications, and to make these simulations feasible, we have decreased the resolution elements of the individual clouds with respect to the single cloud models. We model the gas clouds using , , and particles (in contrast with the used before) to evaluate the convergence of our results with resolution. To identify these different resolutions throughout this paper we use 50k, 200k, 500k and 1m, respectively.

2.1 External potential

In the single cloud simulations presented in Paper I and Paper II, the system was treated in isolation, i.e. with no external forces acting upon it. In reality, such a binary would be placed deep in the potential well of a galaxy, where stars and dark matter exert their gravitational influence. To mimic the presence of the galactic spheroid we include an external potential with a Hernquist profile (Hernquist 1990) which can be written as


and an enclosed mass given by


where and are scaling constants that represent the total mass of the spheroid and the size of core, respectively. Based on a typical relation (Magorrian et al. 1998) and a radius-to-stellar-mass relation (Dabringhausen et al. 2008), we choose and , thus implying


As a consequence, this external force represents only a small perturbation to the dynamics of the interaction that is still dominated by the binary’s gravitational potential. Moreover, because this external potential is spherically symmetric, orbits experience precession only, where energy and angular momentum are conserved. In principle, this precession could enhance stream collisions, prompting more gas onto the binary, although the timescales are too long with respect to the frequency of cloud infall, which is the dominant effect in our models. Based on this, we do not expect the oscillations induced by this potential to have a secular effect in the overall evolution of the system during the stages modelled. In practice, this potential maintains the binary close to the centre of the reference frame by applying a restoring force when it drifts away due to the interaction with the incoming material. It is important to mention that the gas also feels this external force.

For the single cloud simulations this potential was not needed, as drifts produced by the interaction with the cloud or by numerical inaccuracies did not affect the prompt phase studied. For the simulations presented here, however, this restoring force guarantees that the orbit of each cloud remains close to the originally calculated values, as the pericenter distance for each cloud was computed assuming the binary’s centre of mass static at the origin of the Cartesian coordinates.

2.2 Initial conditions

The binary consists of two sink particles, initially having equal masses and a circular orbit in the plane. On the other hand, all the clouds are initially spherical with uniform density, a turbulent velocity field, and a total mass 100 times smaller than the binary (i.e. ). And similar to the single cloud models, the turbulent velocity field is drawn from a random distribution with power spectrum . It is worth mentioning that the random seed is exactly the same for all clouds. With these choices all clouds are identical except for their initial orbit, which simplifies the comparison between the different distributions.

Figure 1: Angular momentum orientations of the initial orbit for 20 randomly selected clouds from different distributions, as indicated at the top of each projection. The binary’s angular momentum initially points to the north pole on these projections. Note that and cloud orbits are generated by mirroring the inclination angle from with respect to the binary’s orbital plane (red solid line) to swap between one hemisphere to the other. Labelled numbers on each circle indicate the order in which we introduce the clouds.
Figure 2: Distribution of pericentre distances as a function of the time at which each cloud is included to the simulation. This distribution is used for the 3 levels of anisotropy shown in Fig. 1. The labelled numbers indicate the order of the clouds. The right plot shows the cumulative distribution of pericentre distances.

The initial setup of these simulations is motivated by the Monte Carlo models presented in Paper II, and is described in full detail in Maureira-Fredes et al. (2018). We use the -distributions to model the different levels of anisotropy of the cloud distribution, where the number represents the fraction of counter-rotating events with respect to the binary (southern hemisphere of the projections shown in Fig. 1). On the other hand, we use a flat distribution of pericentre distances, which derives from an uniform disposition of clouds and efficient gravitational focusing. We generate a total of 20 clouds with orientations shown in Fig. 1 and periapsis shown in Fig. 2. Note that the same pericentre distances are used for the 3 different levels of anisotropy .

We set the maximum value of the pericentre distance to , to ensure a transient evolution with the first impact independent of the cloud’s relative inclination (see Paper II). Finally, the time difference between events exhibited on figure 2 was determined from a Gamma distribution with shape 2 and a scale parameter of . The average time between successive injections is . For a binary with and pc, this rate of clouds corresponds to a mass inflow of  yr entering the inner parsec.

In order to make a robust comparison of the dynamical impact of each distribution it is necessary to reduce some of the stochasticity produced by the low number of events. With that goal in mind, the different angular momentum distributions were not sampled independently from each other. Only the angles from the isotropic distribution () are randomly drawn. We then generate the co-rotating distribution () by reflecting symmetrically the negative inclination angles () with respect to the binary’s orbital plane onto the northern hemisphere. The same is done to get the counter-rotating distribution (), but mirroring the angles onto the southern hemisphere. The azimuthal angles are kept unchanged. This means that the only difference between distributions is the inclination angle of the cloud’s orbits.

Because we are modelling 3 levels of anisotropy for the incoming clouds and we are performing runs at four different resolutions, we are presenting a total of 12 simulations, corresponding to RunB in Maureira-Fredes et al. (2018). We concentrate our analysis on this set of runs because we were able to evolve all the three 50k simulations up to the 20th cloud, thus providing a consistent dataset for the analysis. RunA of Maureira-Fredes et al. (2018) yields similar results, which we summarise in §4.

All the simulations described in this paper are run in different ‘segments’. Initially one cloud is placed with the corresponding initial position and velocity vectors. This system is evolved until the time set for the following cloud to enter has been reached; at this point we stop the simulation. The new cloud is added to the final snapshot to generate the initial conditions for the following section.

3 Dynamical evolution of the binary

3.1 Angular momentum

Figure 3: Angular momentum evolution for the lower resolution runs (50k). The upper panels of each plot show the evolution of the angular momentum magnitude of the entire system (red solid lines) and the binary (dashed blue lines). The dotted light blue lines represent the contribution of the semimajor axis evolution to the binary’s angular momentum, while the dotted-dashed orange lines are the contribution of the total mass change, as shown in equation (5). The lower panels display the relative change between consecutive snapshots of the same lines. The vertical shaded areas indicate the arrival time of the clouds. The extreme “jumps” seen in the red lines are caused by new clouds entering the system, thus providing additional angular momentum. The top row of plots show the evolution during the entire simulations, while the bottom row is a zoom-in on the first orbits, which shows in detail the interaction with 2 clouds.

The first step is to establish whether the simulations have the accuracy to measure the effects of the interaction between the binary and the sequence of clouds. With this goal in mind, we measure the conservation of angular momentum in the different models, as shown in Fig. 3. This corresponds to the evolution in the lower resolution simulation, but because the binary evolution is almost independent from the number of particles (see below), these examples illustrate the evolution of the higher resolution models as well.

Consistent with the single cloud models, the largest effect in the binary angular momentum occurs during the interaction with the first incoming material, depicted as vertical shaded areas. The fluctuations following the first encounter with each cloud rapidly decrease. It is important to clarify that the rather extreme jumps observed in the total angular momentum are not numerically driven, but produced by the addition of new clouds into the system. From the lower panels of figure 3 it is clear that this does not translate into spurious jumps in the binary evolution, as there are no such noticeable changes in its angular momentum. Moreover, simulation’s segments delimited by the addition of clouds are extremely flat, i.e. the numerical fluctuations in the total momentum are way below the observed binary evolution and consequently they can hardly be seen at this scale. This implies that the new models, thanks to the frequent tree updates, have the numerical accuracy to resolve the binary evolution beyond the prompt phase, as opposed to the old single cloud simulations.

From the top panels of figure 3 it is possible to observe that the evolution of the binary is somewhat following the injection of momentum in the form of clouds. In other words, the binary increases its angular momentum in the model, where as in the model the opposite occurs and it decreases. For the isotropic case () the angular momentum increases and decreases depending on whether the new cloud is co-rotating or counter-rotating, respectively. Although we expect a monotonically increasing angular momentum in the simulation, the small drops caused by the inclusion of clouds 11, 12, 16 and 19 occur because the binary inclination slightly changes during the simulations. Said clouds, being so close to the initial binary’s orbital plane (see Fig. 1, upper panel), are now located in the southern hemisphere of the binary at the time of inclusion.

The behaviour of the binary’s angular momentum is in agreement with our previous results where the main effect of the infalling material is to exchange its angular momentum through capture and accretion (Paper II). It is important to bear in mind that an increase in the binary angular momentum does not necessarily imply an expansion of the semimajor axis; this only means that the contribution of the total accreted mass, which is always positive, is much larger to that of the likely negative semimajor axis. This can be understood similarly to what is shown in Paper II, where we write the binary angular momentum in terms of its parameters as follows:


where is the gravitational constant. By differentiating we find that at first order this expression implies


Similar to the single cloud models, during all of our simulations the binary remains roughly equal-mass and circular. We can therefore calculate approximately the total change in angular momentum based on the change in and only,


The contributions of both semimajor axis and mass terms are displayed in the upper panels of Fig. 3. From there it is clear that the addition of these two terms is sufficient to reproduce the binary’s angular momentum evolution, justifying our approximation of ignoring the mass ratio and eccentricity contributions. It is also evident that the contribution to the MBHB angular momentum change is always negative, indicating that the binary shrinks also when .

3.2 Orbital elements

Figure 4: Evolution of the binary semimajor axis (top panels), eccentricity (upper-middle panels), total accreted mass (lower-middle panels), and accretion rate (bottom panels) for all the simulations. The corresponding distribution is indicated at the top of each plot. The different lines in each panel represent the cloud resolutions, as indicated in the legend. The vertical shaded areas indicate the arrival times of the clouds. Note that the non-solid lines are arbitrarily shifted either upwards or downwards with respect to the solid line (50k resolution) in the semimajor axis, accreted mass and accretion rate panels for clarity. The virtually identical evolution during these stages indicates convergence with the number of particles per cloud.

On the basis of the adequate angular momentum conservation, it is possible to study the evolution of the binary orbital elements during the entirety of the simulations. To properly include the influence of the external potential on the binary, we compute its orbital elements ( and ) using the periapsis and apoapsis, instead of the relative energy and angular momentum of the two MBHs. Figure 4 displays the time evolution of the semimajor axis, eccentricity, total mass and accretion rate comparing the different resolutions modelled. Because the different resolutions tend to lie on top of each other, we arbitrarily shifted the non-solid lines with respect to the solid line (50k resolution) for the semimajor axis, accreted mass and accretion rate, which is indicated by black arrows. The virtually identical orbital evolution during the modelled stages indicates convergence with the number of particles per cloud. This is because the dynamical evolution of the system is dominated by capture and gravitational slingshot, which are properly resolved even with a low number of particles. The largest differences are present in the eccentricity evolution. This is due to the typically low values () of this quantity that makes it much more sensitive to noise. For this reason we restrict the analysis to the semimajor axis and mass evolution. In any case, all our binaries remain essentially circular over the modelled timescales (although they consistently show a mild tendency to increase eccentricity), and the small discrepancies in eccentricity do not affect the overall evolution of the system.

() () () () ()
F0.0_50k 145.4 -0.160 0.158 -0.079 0.056
F0.0_200k 95.7 -0.157 0.110 -0.075 0.057
F0.0_500k 79.7 -0.094 0.062 -0.079 0.057
F0.0_1m 68.2 -0.083 0.057 -0.078 0.056
F0.5_50k 124.8 -0.327 0.165 -0.218 0.064
F0.5_200k 86.7 -0.266 0.094 -0.222 0.065
F0.5_500k 69.7 -0.233 0.069 -0.221 0.064
F0.5_1m 67.1 -0.227 0.070 -0.221 0.067
F1.0_50k 123.1 -0.713 0.164 -0.405 0.070
F1.0_200k 81.5 -0.555 0.099 -0.428 0.073
F1.0_500k 62.3 -0.338 0.059 - -
F1.0_1m 42.7 -0.207 0.033 - -
Table 1: Total change of the binary semimajor axis () and mass () for the different models at the final time () of the simulations. The model name is composed of the distribution followed by the resolution of the individual clouds. is the final time reached by each simulation. The values with subscript ‘10th’ correspond to the total change after the interaction with the 10th cloud, which is measured at . This point was not reached by the higher resolution runs of distribution.

In Table 1 we compile the total change in both binary semimajor axis and total mass for all the simulations. The first 3 columns show the values measured at the end of each run (), As expected, increasing the number of particles per cloud makes the simulations more computationally expensive, forcing us to stop them earlier in the evolution of the system. Additionally, the models are more expensive than the others because counter-rotating clouds tend to form more clumps. In order to facilitate the comparison between different resolutions, the values with subscript ‘10th’ correspond to the total change after the interaction with the 10th cloud, measured at , point that was reached by most of our simulations. Notice that, for a given distribution, the discrepancies on the binary parameters at this point are very small, being usually less than 5% with respect to the actual changes measured along the simulation.

Based on the results showed in Fig. 4 and the convergence of the simulation with resolution, most of the analysis in this paper is based on the 50k runs because more clouds have interacted with the binary. We expect the conclusions we will reach will apply for the higher resolution simulations.

Figure 5: Time evolution of the binary semimajor axis (top panel) and mass (bottom panel) in the 50k runs for the different cloud distributions, as indicated in the legend. The solid lines depict the evolution measured directly from the simulations, while the dashed lines correspond to predictions based on the results from the single cloud model presented in Paper II. Similar to previous figures, the shaded vertical areas illustrate the start of the interaction of the binary with a new cloud.

To highlight the differences between distributions, in the upper panel of Fig. 5 we display the semimajor axis evolution for the three levels of anisotropy, while the lower panel shows the mass evolution. The solid lines on this figure represent the evolution measured from the simulations. As expected, the binary shrinks more efficiently when we increase the fraction of retrograde events. For instance, after the infall of 10 clouds the orbit shrinks by in the simulation, by in the model and by a mere for (see Table 1). This can be understood in terms of the exchange of angular momentum with the captured material, as studied in detail in Paper II. Because retrograde clouds bring negative momentum to the binary, any captured and accreted material would immediately subtract angular momentum, while aligned encounters tend to actually expand the MBHB orbit. In contrast with the semimajor axis, the total mass evolution remains very similar for the three levels of anisotropy, as the capture of material does not depend strongly on the orbit of the incoming gas. However, because the gas exchanges its angular momentum with the accreting black holes, the net effect in the binary orbit does depend on the different relative inclinations.

In the Monte Carlo models presented in Paper II, we have extrapolated the results from the single cloud simulations to evolve a binary embedded in a clumpy environment where individual gas pockets infall with different levels of anisotropy. One of the critical assumptions made in this approach was ignoring the effects of the left-over material after the first stages of the interaction. To investigate this effect we compute the binary evolution following the same procedure explained in Paper II, i.e. estimating the total change in mass and semimajor axis based on the cloud’s initial orientation and pericentre distance, using the values extrapolated from the single cloud models. In this procedure we are ignoring any secular effects from the non-accreted gas (as well as cloud-cloud interactions), by considering only the evolution measured during the prompt accretion phase. The results obtained are shown with dashed lines in Fig. 5.

During the interaction with the first few clouds the single cloud predictions agree well with the measured binary evolution, but after few events this simple model consistently underestimates both semimajor axis and mass evolution. This departure allows us to directly measure the effects of the surrounding material, previously unresolved with the Monte Carlo models in Paper II. In the following subsections, we analyse two effects: i) enhancement in accretion due to the interaction of the new incoming material with the gas already in orbit, and ii) gravitational torques from long-lived circumbinary gas structures.

Enhanced Accretion

Figure 6: Time evolution of the accreted mass contributed by each cloud in the different runs. Note that the mass is presented in unit of , which corresponds to the initial mass of each cloud, thus the -axis can be interpreted as the fraction of the cloud that has been accreted by the binary. Here we observe that some of the incoming clouds produce spikes also in the previous clouds.

One striking difference with respect to the single cloud extrapolations is the total accreted mass by the binary. From the lower panel of Fig. 5 we notice that the single cloud model underestimates the final mass by a factor of 3. This occurs because there is an accumulation of material around the binary, and the new incoming clouds are able to drag some of it, increasing accretion and consequently the effect on the binary semimajor axis. Note that the deviation of the measured evolution from the predicted one is episodic rather then continuous.

To corroborate this point, we present the contribution from each individual cloud to the mass accreted by the binary in Fig. 6. Here it is possible to observe that some of the incoming clouds produce an accretion spike also in some of the previous clouds. For example, one of the most noticeable cases is the arrival of cloud 8, which brings to the binary a significant fraction of clouds 1, 3 and 6. Recall from Fig. 2 that the 8th cloud has the smallest pericentre distance, very close to zero. This event produces the largest accretion in all distributions. This extra-accretion triggered by interactions with leftover gas from previously infalling clouds is obviously not captured by the simple model derived from the single cloud simulations.

Figure 7: Time evolution of the ratio between the total accreted mass measured from the simulations and the value predicted from the single cloud models. The horizontal dotted line indicates . This plot shows that the departure from the predicted accreted mass is stronger when there are circumbinary structures present.

Finally, in order to quantify the total increase of the accretion due to the accumulation of gas around the binary, we compute the evolution of the ratio between the total accreted mass measured from the simulation and the value predicted using the single cloud models, which is shown in Fig. 7. As previously observed, the prediction and the actual evolution agree reasonably well during the first few clouds, and this ratio oscillates around 1. After the build up of material around the binary, however, there are large bursts of accretion, reaching values up to times the predicted level. This picture becomes stronger if we consider that around , following the infall of the 10th cloud, there are circumbinary discs formed in all the 3 distributions (see next subsection). These discs are disrupted at least partially by the new incoming clouds, promoting accretion onto the binary.

In a companion paper (Maureira-Fredes et al. 2018), we are using the same suite of simulations to study in detail the formation and evolution of gaseous structures arising from these sequences of infalling clouds onto a binary. In particular, we focus our analysis in the dependency with the level of anisotropy of the incoming clouds (expressed in terms of the -distributions) and the time difference between events. There we show that is typically very difficult to build a substantial circumbinary disc in this scenario of incoherent accretion events, especially with shorter times in between clouds, because the incoming gas constantly perturbs the arising structure. Nevertheless, we do observe the formation of transient circumbinary structures that exert some gravitational influence on the binary’s orbit, as we discuss below.

Circumbinary disc formation and secular evolution

As thoroughly discussed in Paper II, another effect expected from the left-over gas is the formation of circumbinary structures. This represents one of the most important limitations of the models based on the single cloud simulations. For instance, in the extreme case where every cloud comes exactly coplanar and corotating with the binary (sometimes referred to as the ‘prolonged accretion’ scenario), the extrapolation from the single cloud models predicts that the semimajor axis would only increase. In reality, however, following several of such accretion events a massive circumbinary disc is expected to form. This structure can continue extracting and transporting the binary angular momentum, evolving the binary towards coalescence.

Figure 8: Column density maps of an intermediate output from models (left), (centre) and (right). For reference, the binary lies on the -, rotating counter-clockwise. The top and bottom panels show a face-on and edge-on view of the circumbinary structure formed after the interaction with 10 gaseous clouds. For the and models, the disc is prograde with respect to the binary orbit, while for the model the structure is retrograde.

Following the arrival of the 10th cloud, there is enough time before the next event for the gas to settle in a well-defined coplanar circumbinary structure in the 3 distributions. This can be clearly seen on the density maps shown in Fig. 8. For the and models, we observe the formation of a corotating circumbinary disc, while for the model the structure is a retrograde and much narrower, ring-like circumbinary disc.

Figure 9: Density maps averaged over 5 binary orbits for the (left), (centre) and (right) models. These maps were computed in a reference frame corotating with the binary. Note that the large amount of clumps seen in these panels are due to few clumps caught at different orbital phases when creating the maps by superimposing the snapshots. For the and distributions, there are gaseous streams leaking through a cavity sustained by the gravitational influence of the binary. In contrast, for the distribution the gas is able to orbit much closer, directly impacting the black holes.

To study any persistent structures of these developing discs we present the column density maps averaged over five orbits (100 snapshots) in Fig. 9. These maps were computed in a reference frame corotating with the binary, with the black holes aligned on the -axis. It is clear that the gaseous structures are better defined in the simulation, where the binary has opened a clear cavity and each MBH is surrounded by prominent mini-discs. In contrast, the less coherent infall of material in the simulation delays the emergence of these features. Nevertheless, in both cases there is a dense, roughly circular ring of gas located at . In this region the material piles up due to resonances with the gravitational potential of the binary (Artymowicz & Lubow 1994). Although the gravitational forces keep most of the gas confined at this radius, some material leaks thorough the cavity in the form of narrow streams that reach the black holes, a result that has been extensively reported in the literature (e.g. Artymowicz & Lubow 1996; MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Roedig et al. 2012; D’Orazio et al. 2013; Farris et al. 2015; Dunhill et al. 2015; Tang et al. 2017). Such inflowing gas can exert a net negative torque inside the cavity region, shrinking the binary more effectively than the resonant torques.

In contrast, we observe the formation of a counter-rotating circumbinary ring in the run. Because of the absence of gravitational resonances for counter-rotating material, the gas is able to orbit much closer to the binary, directly impacting onto the black holes. As a consequence, the morphology of the gas is quite different, as there are no mini-discs, nor gas in the form of streams. These features are extensively discussed in Maureira-Fredes et al. (2018).

Figure 10: Torque radial profiles averaged over the last five binary orbits. In each panel, the green dashed line represents the differential torque, while the red solid line is the total integrated torque. This latter is integrated starting from the corotation radius (vertical dotted line) inwards and outwards. The quantities are presented in code units, namely, for the differential torque and for the integrated profile.

To establish the role of the discs in the binary evolution we compute the gravitational torques exerted by the gas. The total gravitational torque is given by directly summing the individual particle torques onto each MBH as follows


where the index runs over all the gas particles in the simulation, and indicates the two black holes. Using this expression, we compute the torque radial profiles shown in Fig. 10, averaged over 5 orbits of the simulation. As the binary’s inclination slightly evolves in our simulations, we compute these torques in a reference frame in which the binary angular momentum is exactly aligned with the axis. With this choice, we can directly link the component of the torque with the binary orbital evolution (eq. 5). Because the torques are null at the corotation radius, the integration of the differential profile is done from that point outwards and inwards.

From the profiles shown in figure 10 it is clear that the net effect of the circumbinary material is to produce a negative torque, hence to extract some angular momentum from the binary. For the co-rotating discs, the largest peak is located between and , i.e. inside the cavity wall. This peak basically determines the total strength of the torque, which is responsible for the shrinking observed in the binary orbit. Beyond the cavity region the torques oscillate between positive and negative, cancelling each other out, and consequently making the contribution of the disc negligible. For the counter-rotating disc, conversely, the gravitational torques come from the material inside the corotation radius, very close to the MBHs, and the outward torques are essentially absent because of the lack of Lindblad resonances.

Figure 11: Averaged surface density torque exerted by the gas onto the binary for the (left), (centre) and (right) models. Similar to the previous figures, these maps were averaged over five orbits. Note that the torque values for the minidiscs are off the scale in both and models. For the prograde discs this shows that the negative torques come from the gaseous streams reaching the binary. Conversely, in the counter-rotating disc these torques are located much closer to the binary.

To identify the location of the negative torques, we show the averaged surface density torque of the discs in Fig. 11. For the co-rotating circumbinary discs, and similar to the gas distribution shown in figure 9, the location of the torques are clearer in the model, where the cavity has been almost completely depleted of gas. Nonetheless, in both and cases the gaseous streams inside the cavity provide a strong source of negative torque. This confirms that the negative peak observed in the profiles of figure 10 comes from the material leaking through the cavity wall. Rather than being entirely captured by the binary, a fraction of this gas is flung back to the disc, via gravitational slingshot, carrying away angular momentum from the binary. This physical effect has been identified and extensively studied in the literature, mainly by means of numerical simulations (e.g. Roedig et al. 2012; D’Orazio et al. 2013; Farris et al. 2015; Dunhill et al. 2015; Tang et al. 2017).

On the other hand, gravitational torques in the counter-rotating disc are located much closer to the binary, with no features resembling the gaseous streams as in the previous cases. This material extracts some of the binary’s angular momentum, but because it is not able to transport it outwards, this material is eventually accreted by the MBHs. As this gas is counter-rotating with respect to the binary, its capture and subsequent accretion subtracts angular momentum much more efficiently than in the co-rotating case (Nixon et al. 2011).

The profiles shown in figure 10 can be directly compared with results from previous studies of \enquotestandard circumbinary discs. In particular, the simulations of a massive self-gravitating circumbinary disc by Roedig & Sesana (2014) find very similar torque profiles (cf. their Fig. 4), the main distinction being the strength of the peaks – the torques found in our simulations are noticeably larger. For instance, for prograde discs, Roedig et al. show a total outwards torque of , while in our simulations we find it to be , hence a factor of larger. In the retrograde case, Roedig et al. find a total integrated torque of , while we find , a factor of larger. Even though these values might not seem much higher, note that we are comparing very massive circumbinary discs with 20% the binary mass to the much lighter discs formed in our simulations, that reach at most a few percent of the MBHB mass (Maureira-Fredes et al. 2018).

These large torques are due to the transient nature of our discs; during the episodes of disc disruption, the binary clear its orbit of gaseous material, enhancing the torques. If these circumbinary discs were to settle in a more steady-state, such as the model of Roedig & Sesana (2014), the gaseous streams will be determined mainly by the internal viscous torques, and we expect much less mass interacting directly with the binary. This would decrease the efficiency of the gravitational torques, especially considering their low mass with respect to the MBHs. To confirm this point, we run a separate set of simulations of these binaries surrounded by their discs without adding any more clouds – we referred to these runs as ‘forked’ simulations in Maureira-Fredes et al. (2018). Because we are continuing further the simulations after the 10th cloud, we call here these the extended 10th cloud models (10x).

Figure 12: Time evolution of the integrated torque profile of the circumbinary discs, for the extended 10th cloud simulations. Each line corresponds to a time average over 5 orbital periods, taken further in the binary evolution from the lightest to the darkest line. From left to right we show the , , and distributions. The trend shows decreasing gravitational torques as a function of time.
Figure 13: Total integrated torque evolution measured during the extended 10th cloud simulations. Points correspond to the average of the profiles shown in Fig. 12. The light blue line shows the evolution of the distribution, the orange line is the distribution, and the green lines is the distribution. This shows that the slow relaxation of the discs translates into decreasing gravitational influence.

Using these simulations we measure how the total integrated torque profile evolves in time for each of the discs, as we show in Fig. 12. Similar to Fig. 10, each line is the radial profile of the torque’s component averaged over 5 binary orbital periods. From this figure it is clear that the gravitational influence of each disc tends to decrease as time passes (from lighter to darker lines), showing that the measured evolution driven by torques is still transient and not a long term effect. For the co-rotating circumbinary discs (left and middle plots of Fig. 12) the shape of the profile remains rather unperturbed, which means that the discs evolve smoothly. This is not the case for the counter-rotating disc (right plot of Fig. 12), where we observe strong changes in the profile outside the corotation radius. This occurs because this circumbinary ring-like structure has little angular momentum, and hence it can be perturbed much more easily by gaseous clumps and streams reaching the binary from previous interactions. Nevertheless, the profile inside the binary’s orbit remains roughly unperturbed, just decreasing its strength with time. Consequently, the evolution of the torque profiles displayed in Fig. 12 clearly shows the transient nature of the large torques found in our simulations.

The relaxation of the discs becomes even more apparent when we compute the total torque (averaged over five orbits) as a function of time, as shown in Fig. 13. The trend is clear: the net gravitational effect of the discs decreases in time. This damping of the torque is occurring somewhat slowly in our simulations, but this is likely due to the fact that we are not including any source of viscosity on top of the artificial term to capture shocks. This artificial viscosity acts as an effective viscosity within the disc, but is typically very small (Lodato & Price 2010). We would expect these discs to relax much faster with a more realistic prescription of the viscosity.

Nonetheless, we can use the formalism presented in Roedig & Sesana (2014) to make an order of magnitude estimation of total effect of such discs in the binary orbit. If we take the simple assumption that the gravitational torques evolve only the semimajor axis, we can write


where for our equal-mass and circular binary. By taking an average torque of , this expression yields


which over the dynamical times that span our 50k simulations translates onto a total evolution of


Comparing this value with the final of each distribution (Table 1), we find that corresponds to 40% for the , 20% for the , and 10% for the distribution. Even though this seems to indicate that the circumbinary structures contribute with an important fraction of the evolution, we arrived to this number by considering the presence of the circumbinary disc during the entirety of the run, which is certainly not true. As a matter of fact, these discs arose after the arrival of the 10th cloud, and only the one from the distribution survives the interaction during the rest of the simulation, the other two are destroyed by the interaction following a few clouds (Maureira-Fredes et al. 2018). Consequently, this estimate demonstrates that the extraction of angular momentum via gravitational torques from a circumbinary disc is subdominant with respect to capture and accretion of material, particularly for the and distributions, where these structures are unlikely to survive.

4 Additional run: the effects of shorter arrival times

Figure 14: Compilation of the initial conditions for RunA. This is similar to the distributions shown in Figs. 1 and 2, where on the three top plots we display the cloud’s angular momentum orientations for the different levels of anisotropy, while the bottom plot shows the distribution of pericentre distances and the as a function of the time of inclusion.

Due to the stochasticity intrinsic to the low number of injected clouds, significant scatter in the MBHB evolution can be expected in different realisations of the interacting cloud population. We thus analysed a second suite of simulations, corresponding to RunA in Maureira-Fredes et al. (2018). Initial conditions for this set of runs are shown in Fig. 14. As described in Maureira-Fredes et al. (2018), due to computational limitations we were unable to achieve the 20th cloud for all the distributions, and we integrated RunA up to cloud 15, RunA up to cloud 22 and RunA up to cloud 12. To facilitate the comparison with the previous results, here we present only up to cloud 20.

Besides featuring a different random seed for the generation of initial pericenter and angular momentum distributions, this set of simulations is characterised by a much shorter average elapsed time in between clouds. This can be clearly seen by comparing the time of the final cloud () in both runs – while in RunB the the 20th cloud is introduced at around 120  (Fig. 2), for RunA this time is 60  (Fig. 14, bottom plot). This translates onto an average time difference a factor of 2 smaller for the latter.

Figure 15: Same as Fig. 5, but for RunA.

One of the most significant results presented in the previous section is the deviation of the binary evolution from what predicted by the single cloud model. A difference that clearly shows the effects of the non-accreted material, which was not taken into account in that model. We repeated the analysis for RunA, results are shown in Fig. 15. Most of the conclusions drawn for RunB also apply to this case. For instance, although the the single cloud model is able to roughly reproduce the evolution driven by interaction with the first few clouds, it consistently underestimates the cumulative binary shrinking as further clouds are added to the system. A striking difference between the actual evolution and the prediction of the single cloud model is observed in the case, in which the binary slightly expands instead of shrinking its orbit.

Figure 16: Same as Fig. 7, but for RunA.

Similar to RunB, we observe that the total accreted mass measured in the multicloud simulation is about a factor 3 larger than the value predicted by the single cloud models. Consequently, the dominant effect of the surrounding gas on RunA also appears to be the enhanced accretion caused by the interaction with the new incoming clouds. To quantify this increase in accretion, we directly compare figures predicted by the single cloud model with those measured from the simulations in Fig. 16. This comparison clearly shows that the single cloud model underestimates the accreted mass onto the binary after the first few clouds even more consistently than in RunB. This occurs because the typically shorter times between events does not allow for material to settle into a disc. This is one of the main results shown in Maureira-Fredes et al. (2018) – in order for the gas to settle into a more stable circumbinary configuration it needs to gain some angular momentum, which takes a few binary orbits. When clouds interact frequently with the binary, new incoming clouds strongly perturb the surrounding gas structures, prompting a lot of accretion onto the binary. In agreement with RunB, this is the main source of deviation with respect to the single cloud predictions shown in Fig. 15.

Figure 17: Same as Fig. 10, but for the circumbinary structures observed during RunA. Please note that the circumbinary dics responsible for the represented torques were identified on different times depending on the distributions: after the 12th cloud for , 16th for and 10th for .

Despite the higher infall rate of clouds in RunA, we do observe also the formation of somewhat coplanar circumbinary discs, albeit they get continuously destroyed by the interaction with new incoming gas (Maureira-Fredes et al. 2018). Opposite to RunB, those discs are less prominent and on average much lighter. Nonetheless, we can compute the gravitational influence of these circumbinary structures through the torque that the gas exerts on the binary. The differential and integrated radial torque profiles are shown in Fig. 17.

Compared to the profiles shown in Fig. 10 we observe similar features, confirming the robustness of the dynamics against different realisations of the interacting cloud population. In particular, in the and runs, material outside the corotation radius is transporting some of the binary’s angular momentum outwards via gravitational slingshot, while for the distribution all the negative torques come from the material inside the binary orbit.

Interestingly, in the simulation the positive torque inside the corotation radius is larger than the outwards torque, which implies that the net effect of the gravitational forces is to add momentum to the binary, expanding its orbit. This is likely why the final semimajor axis during this simulation is larger than the initial value, in contrast to the predicted evolution which is a slight shrinking of the orbit. The continuous disruption of arising structures due to frequent infall of clouds, supplies a larger amount of gas inside the binary’s orbit, enhancing these positive torques. If the accretion events were to stop, or decrease in frequency, the surviving material would settle into a circumbinary disc, and as the binary opens up a cavity, gravitational slingshot of gaseous streams would becomes the dominant effect, shrinking its semimajor axis.

Overall, the main conclusions found in the previous section apply to this run. The dominant effect of the non-accreted material on the binary evolution comes from the enhanced accretion produced by stream interactions. Due to the higher infalling rate we find that is even more difficult to build massive and stable circumbinary structures (Maureira-Fredes et al. 2018), which further decreases the influence of the angular momentum transport by gravitational torques. The net result is that the dynamics of the system is mostly driven by capture and accretion of mass, which has the interesting effect of making the evolution of the binary even more sensitive to the details of the cloud distributions. In general, for isotropic or counter-rotating distributions, the mutual partial cancellation of the MBHB and the accreted gas angular momenta accelerates the binary shrinking, driving the binary to a faster coalescence. In the co-rotating case, the resulting binary evolution depends on a subtle balance between accreted mass and angular momentum, with the net effect that the MBHB can also slightly expand. A longer simulation is required to confirm whether the expansion seen in RunA is only transient or if it holds in the long term.

5 Summary and discussion

On this paper we have introduced a new suite of simulations to study the evolution of a binary embedded in a turbulent environment, where we expect several gaseous clumps to infall and dynamically interact with the binary. The main goal was to overcome the limitations of the Monte Carlo evolution presented in Paper II by directly modelling the interaction of a MBHB with a sequence of incoherent accretion events. These models combine the ‘discrete’ evolution due to the prompt interaction of each individual cloud, with the secular torques exerted by the left-over gas. The latter could not be numerically resolved in the single cloud models, but could have a role on the dynamics of both gas and binary.

In order to directly compare with the results of the Monte Carlo evolution of Paper II, the cloud orbits where sampled using the same probability density functions. In particular, we have used the -distributions for the angular momentum orientation which are an indication of the anisotropy levels in the surrounding medium. These distributions have been linked to the morphological and kinematical properties of the host galaxy (Sesana et al. 2014).

Our main findings can be summarised as follows:

  1. The MBHB angular momentum evolution follows the injection of momentum in the form of clouds. That is, the binary monotonically increases its angular momentum (even though it still shrinks) in the model where each new cloud adds to the total momentum budget, while continuously decreasing for where the opposite occurs. This is consistent with the scenario where the orbital evolution is dominated by direct exchange of angular momentum through capture and accretion, as demonstrated by the single cloud simulations in Paper II.

  2. The binary exhibits an almost identical orbital evolution at every investigated resolution. This is because the dynamical evolution of the system is dominated by capture and gravitational slingshot, which are properly resolved even with a low number of particles.

  3. The largest differences are present in the eccentricity evolution due to typically low values (), making it much more sensitive to noise. Because of this, we focused the analysis on the semimajor axis and mass evolution, considering that the binary remains circular throughout the whole simulation.

  4. Also concordant with the single cloud simulations, the evolution of the semimajor axis is strongly dependent on the distribution of the infalling gas, being fastest when it features more counter-rotating clouds. Conversely, the binary mass evolution is roughly independent of the cloud distribution during most of the interaction.

  5. By directly comparing the observed binary evolution with a semi-analytical model constructed on the basis of the single cloud simulations, we found that accretion onto the MBHB increases considerably (by about a factor of three). This occurs because new incoming material is able to drag inward a fraction of the surrounding gas leftover by previous accretion events.

  6. Enhanced accretion has a strong impact on the evolution of the binary. Because of angular momentum cancellation with a larger mass of infalling gas, the binary evolution is accelerated in the isotropic and counter-rotating cases (by 20-to-50% and a factor 2, respectively, see Figs. 5 and 15). In the co-rotating case, the evolution depends on a fine balance of mass and angular momentum accretion, and can lead, at least temporarily, to the expansion of the binary semimajor axis.

  7. In the and models we observe the intermittent formation of well-defined, coplanar and co-rotating circumbinary discs. The discs present well-known features, in particular, a cavity located at which is maintained by the resonances of the orbiting gas with the binary’s gravitational potential.

  8. Similarly, we observed the formation of intermittent coplanar and counter-rotating circumbinary discs in the model. Due to the absence of resonances, in this case the gas is able to orbit much closer to the binary.

  9. By computing the gravitational torques exerted by the surrounding gas, we have shown that these discs are able to transport some of the binary’s angular momentum via gaseous streams that are slingshotted away once they get close to the MBHs. This process can continue shrinking the binary orbit in the absence of new impacting clouds. However, due to the low mass of these discs, their effect on the binary orbital decay is subdominant with respect to the capture and subsequent accretion of the incoherent material deposited by clouds in the vicinity of the MBHB.

An interesting case is the isotropic scenario where the total angular momentum brought by the gas is close to zero, at least after a significant number of clouds. In this scenario, the intersection of highly misaligned streams tend to effectively cancel angular momentum, forcing material to plunge onto the binary. Therefore, we do not expect a steady circumbinary disc to survive under these conditions. Following a sequence of 10 clouds, however, there is a clear circumbinary structure rising from the non-accreted material, which is already producing a noticeable effect on the binary orbit. This implies that the formation of circumbinary discs is a relevant process even for the model, even if these discs are transient and with negligible self-gravity. One of the assumptions of the Monte Carlo models presented in Paper II is that there is no circumbinary disc during the isotropic scenario due to the low total angular momentum of the gas. Our results clearly show the inadequacy of this assumption. This occurs because corotating material is able to gain more angular momentum via resonances and gravitational slingshot compared to counter-rotating gas, thus forming meta-stable disc-like configurations.

At this point, it is important to stress that our simulations are subject to a number of limitations, that we now briefly consider. The mass of the circumbinary discs intermittently forming throughout the simulations is small compared to the binary and hence their self-gravity is negligible (Maureira-Fredes et al. 2018). In this situation the outwards transport of angular momentum throughout the disc will occur through an unresolved viscosity (see e.g. Syer & Clarke 1995; Ivanov et al. 1999). For geometrically thin and optically thick discs, such viscosity is commonly parametrised using the -disc model (Shakura & Sunyaev 1973), where represents the strength of the internal viscous stresses. In our simulations, however, the only source of viscosity is the artificial term introduced to capture shocks, hence the transport of angular momentum throughout the disc is likely not adequately modelled. Additionally, the evolution of the binary-disc system might be affected by the simple thermodynamics adopted in these models. As discussed in Roedig et al. (2012), the morphology of the streaming gas has a major impact in the binary evolution and the structure of the surrounding disc, thus it is directly influenced by the gas thermodynamics. Moreover, Tang et al. (2017) using 2D hydrodynamical simulations of circumbinary accretion discs found that the torques are sensitive to the sink prescription: slower sinks result in more gas accumulating near the MBHs, driving the binary to merge more rapidly. Nevertheless, all these effects become relevant for the long-term evolution of a MBHB embedded in a circumbinary disc, where torques mediated by resonances play an important role in the evolution of the system. Conversely, during the infall of a sequence of clouds, the relevant mechanisms driving the binary’s orbital response are almost purely dynamical, namely, gas capture and gravitational slingshot, which are well resolved with these models. Therefore, a more physical treatment for the viscosity, thermodynamics and accretion in required only if one wants to resolve the long-term evolution of the MBHB surrounded by a circumbinary disc arising in the aftermath of the accretion events.

Another delicate point is the simplified treatment of accretion, whereby every bound particle entering the sink radius is added to the MBHs. When scaled to our fiducial physical model (, pc), this results in sparse episodes of highly super-Eddington accretion. The exact fate of the gas during this ‘accretion spikes’ is uncertain. The extremely high density of the accretion flow might cause photon diffusion timescale to be longer than the accretion timescale, thus advecting photons into the MBHs, in turn allowing super-Eddington accretion rates. Conversely, radiation pressure in the proximity of the MBHs might cause copious mass loss through powerful winds, significantly decreasing the amount of accretion. As noticed in Paper II, the exchange of angular momentum is mostly mediated by capture of gas in bound orbits around the individual MBHs rather than subsequent accretion. Even assuming that the majority of the gas is accelerated away in form of isotropic winds, the evolution of the binary is affected at a level of 20-30%. This conclusion, however, does not take into account the interaction of the outgoing wind with the infalling clouds. If winds are isotropic and strong enough to clear the infalling gas on large scales, this would quench the supply of material onto the MBHB, eventually stalling the binary. A proper investigation of this potential issue requires the implementation of a suitable feedback prescription, and we plan to explore this in future work.

A somewhat different aspect that might be impacting our results is the resolution elements of the individual clouds. For instance, in their simulations Dunhill et al. (2015) found that the dynamics of gaseous streams going through the cavity is resolution dependent. As discussed above, the streams play an important role transporting angular momentum once the circumbinary disc has formed, thus the late evolution measured using the lower resolution runs is likely affected by this issue. However, we did not find any resolution-dependent feature in the binary evolution when modelling the clouds with up to 1 million particles. We can therefore claim that our main results are robust.

The results exhibited in this paper suggest that the secular effects of the left-over material further increase the efficiency of the incoherent infalling clouds to bring the MBHs towards coalescence respect to a scenario that only considers the evolution during the prompt accretion phase, as the one presented in Paper II. Consequently, provided there is a reasonable rate of events, these models confirm that infalling clouds present a viable mechanism to efficiently merge these binaries within less then a Gyr.


We thank Volker Springel for his suggestions on how to improve the conservation of angular momentum of these models. We thank Johanna Coronado for comments on the manuscript. The simulations were performed partially between the sandy-bridge nodes at HITS, and datura and minerva clusters at the AEI. FGG acknowledges support from the CONICYT-PCHA Doctorado Nacional scholarship, DAAD in the context of the PUC-HD Graduate Exchange Fellowship, the European Research Council under ERC-StG grant EXAGAL-308037, and the Klaus Tschira Foundation. CMF acknowledges support from the Transregio 7 “Gravitational Wave Astronomy” financed by the Deutsche Forschungsgemeinschaft DFG (German Research Foundation). CMF acknowledges support from the DFG Project “Supermassive black holes, accretion discs, stellar dynamics and tidal disruptions”, awarded to PAS, and the International Max-Planck Research School. PAS acknowledges support from the Ramón y Cajal Programme of the Ministry of Economy, Industry and Competitiveness of Spain, as well as the COST Action GWverse CA16104. This work has been partially supported by the CAS President’s International Fellowship Initiative. JC acknowledges support from CONICYT-Chile through FONDECYT (1141175) and Basal (PFB0609) grants. This work was partially developed while JC was on sabbatical leave at MPE. JC and FGG acknowledge the kind hospitality of MPE, and funding from the Max Planck Society through a “Partner Group” grant.



  1. thanks: E-mail: felipe.goicovic@gmail.com


  1. Alig C., Burkert A., Johansson P. H., Schartmann M., 2011, \mnras, 412, 469
  2. Amaro-Seoane P., Brem P., Cuadra J., 2013, \apj, 764, 14
  3. Amaro-Seoane P., et al., 2017, preprint, (arXiv:1702.00786)
  4. Armitage P. J., Natarajan P., 2005, \apj, 634, 921
  5. Artymowicz P., Lubow S. H., 1994, \apj, 421, 651
  6. Artymowicz P., Lubow S. H., 1996, \apjl, 467, L77
  7. Begelman M. C., Blandford R. D., Rees M. J., 1980, \nat, 287, 307
  8. Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  9. Bondi H., 1952, \mnras, 112, 195
  10. Bonnell I. A., Rice W. K. M., 2008, \sci, 321, 1060
  11. Combes F., Young L. M., Bureau M., 2007, \mnras, 377, 1795
  12. Consortium T. e., et al., 2013, arXiv:1305.5720,
  13. Cuadra J., Nayakshin S., Springel V., Di Matteo T., 2006, \mnras, 366, 358
  14. Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, \mnras, 393, 1423
  15. D’Orazio D. J., Haiman Z., MacFadyen A., 2013, \mnras, 436, 2997
  16. Dabringhausen J., Hilker M., Kroupa P., 2008, \mnras, 386, 864
  17. David L. P., et al., 2014, \apj, 792, 94
  18. Dunhill A. C., Alexander R. D., Nixon C. J., King A. R., 2014, \mnras, 445, 2285
  19. Dunhill A. C., Cuadra J., Dougados C., 2015, \mnras, 448, 3545
  20. Escala A., Larson R. B., Coppi P. S., Mardones D., 2004, \apj, 607, 765
  21. Escala A., Larson R. B., Coppi P. S., Mardones D., 2005, \apj, 630, 152
  22. Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015, \mnras, 447, L80
  23. Ferrarese L., Merritt D., 2000, \apjl, 539, L9
  24. Fiacconi D., Mayer L., Roškar R., Colpi M., 2013, \apjl, 777, L14
  25. Gaburov E., Johansen A., Levin Y., 2012, \apj, 758, 103
  26. Gaspari M., Sa̧dowski A., 2017, \apj, 837, 149
  27. Gaspari M., Ruszkowski M., Sharma P., 2012, \apj, 746, 94
  28. Gaspari M., Ruszkowski M., Oh S. P., 2013, \mnras, 432, 3401
  29. Gaspari M., Brighenti F., Temi P., 2015, \aap, 579, A62
  30. Gaspari M., et al., 2017a, preprint, (arXiv:1709.06564)
  31. Gaspari M., Temi P., Brighenti F., 2017b, \mnras, 466, 677
  32. Gebhardt K., et al., 2000, \apjl, 539, L13
  33. Goicovic F. G., Cuadra J., Sesana A., Stasyszyn F., Amaro-Seoane P., Tanaka T. L., 2016, \mnras, 455, 1989
  34. Goicovic F. G., Sesana A., Cuadra J., Stasyszyn F., 2017, \mnras, 472, 514
  35. Haiman Z., Kocsis B., Menou K., 2009, \apj, 700, 1952
  36. Hernquist L., 1990, \apj, 356, 359
  37. Hobbs A., Nayakshin S., 2009, \mnras, 394, 191
  38. Hobbs A., Nayakshin S., Power C., King A., 2011, \mnras, 413, 2633
  39. Hopkins P. F., Quataert E., 2011, \mnras, 415, 1027
  40. Ivanov P. B., Papaloizou J. C. B., Polnarev A. G., 1999, \mnras, 307, 79
  41. Jahnke K., Macciò A. V., 2011, \apj, 734, 92
  42. Kelly B. C., Shen Y., 2013, \apj, 764, 45
  43. King A. R., Pringle J. E., 2006, \mnras, 373, L90
  44. King A. R., Pringle J. E., 2007, \mnras, 377, L25
  45. Kocsis B., Haiman Z., Loeb A., 2012, \mnras, 427, 2680
  46. Kormendy J., Ho L. C., 2013, \araa, 51, 511
  47. Lodato G., Price D. J., 2010, \mnras, 405, 1212
  48. Lodato G., Nayakshin S., King A. R., Pringle J. E., 2009, \mnras, 398, 1392
  49. Lu J. R., Ghez A. M., Hornstein S. D., Morris M. R., Becklin E. E., Matthews K., 2009, \apj, 690, 1463
  50. Lucas W. E., Bonnell I. A., Davies M. B., Rice W. K. M., 2013, \mnras, 433, 353
  51. Lupi A., Haardt F., Dotti M., Colpi M., 2015, \mnras, 453, 3437
  52. MacFadyen A. I., Milosavljević M., 2008, \apj, 672, 83
  53. Magorrian J., et al., 1998, \aj, 115, 2285
  54. Mapelli M., Trani A. A., 2016, \aap, 585, A161
  55. Mapelli M., Hayfield T., Mayer L., Wadsley J., 2012, \apj, 749, 168
  56. Masunaga H., Inutsuka S.-i., 2000, \apj, 531, 350
  57. Maureira-Fredes C., Goicovic F. G., Amaro-Seoane P., Sesana A., 2018, preprint, (arXiv:1801.06179)
  58. Mayer L., Kazantzidis S., Madau P., Colpi M., Quinn T., Wadsley J., 2007, Science, 316, 1874
  59. McNamara B. R., et al., 2014, \apj, 785, 44
  60. Nayakshin S., King A., 2007, arXiv:0705.1686,
  61. Nayakshin S., Power C., King A. R., 2012, \apj, 753, 15
  62. Nixon C., Lubow S. H., 2015, \mnras, 448, 3472
  63. Nixon C. J., Cossins P. J., King A. R., Pringle J. E., 2011, \mnras, 412, 1591
  64. Paumard T., et al., 2006, \apj, 643, 1011
  65. Roedig C., Sesana A., 2014, \mnras, 439, 3476
  66. Roedig C., Dotti M., Sesana A., Cuadra J., Colpi M., 2011, \mnras, 415, 3033
  67. Roedig C., Sesana A., Dotti M., Cuadra J., Amaro-Seoane P., Haardt F., 2012, \aap, 545, A127
  68. Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn T. R., Wadsley J., 2015, \mnras, 449, 494
  69. Russell H. R., et al., 2014, \apj, 784, 78
  70. Sanders D. B., Mirabel I. F., 1996, \araa, 34, 749
  71. Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, \apj, 794, 104
  72. Shakura N. I., Sunyaev R. A., 1973, \aap, 24, 337
  73. Springel V., 2005, \mnras, 364, 1105
  74. Syer D., Clarke C. J., 1995, \mnras, 277, 758
  75. Tang Y., MacFadyen A., Haiman Z., 2017, \mnras, 469, 4258
  76. Tremblay G. R., et al., 2016, \nat, 534, 218
  77. White S. D. M., Rees M. J., 1978, \mnras, 183, 341
  78. del Valle L., Escala A., Maureira-Fredes C., Molina J., Cuadra J., Amaro-Seoane P., 2015, \apj, 811, 59
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description