Evolution of supermassive black hole binaries and tidal disruption rates in non-spherical galactic nuclei Free

ABSTRACT

1 INTRODUCTION

It is commonly accepted that supermassive black holes (SMBHs) reside in the centres of most galaxies (Kormendy & Ho 2013; Graham 2016). In the conventional cosmological paradigm, larger galaxies are formed through hierarchical mergers of smaller galaxies. Thus the merger remnants may contain multiple SMBHs, which will eventually sink towards the galaxy centre due to dynamical friction, and form a binary (or a higher multiplicity) system. The subsequent evolution of a binary SMBH is driven by a combination of several effects, which all lead to a gradual shrinking of the binary orbit (Begelman, Blandford & Rees 1980, see also more recent reviews by Merritt & Milosavljević 2005; Colpi 2014). Three-body interactions between the pair of SMBHs and stars of the galactic nucleus typically result in an ejection of the lighter of the three bodies (the star), extracting the energy from the binary orbital motion (Quinlan 1996). Torques from the circumbinary gaseous discs also lead to the shrinking of the binary (e.g. Mayer 2013; Rafikov 2016; Tang, MacFayden & Haiman 2017), although this process appears to be effective only when the separation becomes small enough (Lodato et al. 2009). Finally, when the separation between the two SMBHs reaches ≲10⁻³ pc, the emission of gravitational waves (GWs) becomes efficient enough and leads to the ultimate coalescence of the SMBHs.

Observational evidence for the existence of binary SMBHs is scarce, and to date only a few candidate objects have been found through direct imaging (Rodriguez et al. 2006; Kharb, Lal & Merritt 2017). Among various indirect observational signatures (see Komossa 2006; Dotti, Sesana & Decarli 2012 for a review) are tidal disruption events (TDEs), when a star passing too close to one of the SMBH is torn apart by the tidal forces. The presence of a second SMBH modifies the light curve (e.g. Ricarte et al. 2016; Coughlin et al. 2017; Vigneron, Lodato & Guidarelli 2018), and rates of TDEs in galactic nuclei hosting single and binary SMBHs may also differ substantially. In particular, during the early stage of binary evolution (∼10⁵ to 10⁶ yr after its formation), TDE rates could be dramatically elevated, reaching 10⁻² to 1 events per year (Ivanov, Polnarev & Saha 2005; Chen et al. 2009, 2011; Liu & Chen 2013; Li et al. 2017), while for the remaining part of its evolution the rates could be much lower (Chen, Liu & Magorrian 2008). On the other hand, high TDE rates could also result from other factors, such as triaxial geometry of the galactic nucleus (Merritt & Poon 2004; Vasiliev 2014; Lezhnin & Vasiliev 2016), or steep density cusps (Stone & van Velzen 2016; Stone et al. 2018). Finally, a double TDE within a single galaxy (not yet observed), separated by an interval of a few hundred days, would be a possible signature of two stars in a stellar binary being captured by two components of an SMBH binary (Coughlin et al. 2018; Wu & Yuan 2018). If TDE by binary SMBHs can be reliably distinguished from those by single SMBHs, they would serve as a valuable probe of the cosmic population of low-mass binary SMBHs in upcoming large-scale surveys (Thorp, Chadwick & Sesana 2018).

The lifetime of an SMBH binary in a gas-poor galactic nucleus is determined by the rate of its interaction with stars that come to its vicinity: some of them would be tidally disrupted, but most will be eventually ejected with high velocity, carrying away the energy from the binary. Early studies identified a possible bottleneck related to the eventual depletion of the reservoir of interacting stars (so-called final-parsec problem, Milosavljević & Merritt 2001). However, now it is commonly believed to be an artefact of the assumption of spherical geometry. N-body simulations of merger remnants produce less symmetric systems and do not exhibit the stalling of the binary evolution (Khan, Just & Merritt 2011; Preto et al. 2011); simulations of isolated non-spherical systems confirm this conclusion (Berczik et al. 2006; Vasiliev, Antonini & Merritt 2014). The disappearance of the final-parsec problem has been linked to the existence of a sufficiently large population of so-called centrophilic orbits (those which may approach the binary closely enough, although not at every pericentre passage), sustaining the binary hardening and leading to its coalescence on a time-scale ≲1 Gyr (Vasiliev, Antonini & Merritt 2015; Gualandris et al. 2017). Fully cosmological simulations cannot explicitly follow the dynamical evolution of SMBH binaries, but can still be used to estimate their population under reasonable assumptions about unresolved physical processes (e.g. Kelley, Blecha & Hernquist 2017; Ryu et al. 2018); the estimated lifetimes lie in the range from 10⁸ to few ×10⁹ yr.

The stellar-dynamical processes in galactic nuclei with single or binary SMBH can be studied with various approaches, but for several reasons discussed later in this paper, we choose the Monte Carlo simulation method, implemented in the code raga (Vasiliev 2015). It has been used previously to study the TDE rates around single SMBH (Vasiliev 2014; Lezhnin & Vasiliev 2016) and the evolution of binary SMBH (Vasiliev et al. 2015; Vasiliev 2016), but not simultaneously. In this study, we extend it to include both processes, and conduct a large suite of simulations, exploring the influence of several properties of the system: the total mass of the binary SMBH, its mass ratio, initial eccentricity of its orbit, the geometry of the stellar distribution (spherical, axisymmetric, or triaxial), its mean rotation, and the inclusion or neglect of stellar captures. We confirm the earlier conclusions of Vasiliev et al. (2015) and Gualandris et al. (2017) that the final-parsec problem can be avoided only in triaxial systems, and that the hardening rate of the binary due to interactions with stars is almost independent of its eccentricity, but is higher for binaries with a larger mass ratio. We also find that the occasional captures of stars do not have any influence on the hardening rate. The difference between the capture rates in galactic nuclei hosting a binary SMBH or a single SMBH, but identical other properties, depends on the evolutionary stage and the geometry of the stellar distribution. At early times, the capture rates by binary SMBHs are much higher than by single SMBHs, but subsequently they drop to lower values. In spherical or axisymmetric systems they remain very low throughout the evolution, in agreement with earlier studies (Chen et al. 2008), but in triaxial systems they may be comparable to or only moderately lower than the single SMBH capture rates. We examine the properties of stellar orbits before and after interaction with the binary. We demonstrate that the loss cone remains well populated in triaxial systems, but significantly depleted otherwise. Hypervelocity stars ejected from the galaxy by the slingshot interaction have a roughly exponential distribution in velocity, reaching few thousand km s⁻¹. Finally, we study the impact of the binary on the structure of the galactic nucleus, in particular, the erosion of the stellar cusp. We argue that while it definitely takes place, it is very hard to measure observationally the amount of missing mass or to relate it to the properties of the binary.

The paper is organized as follows. We start by summarizing the classical loss-cone theory and its extension to non-spherical galactic nuclei in Section 2, referring the reader to chapter 6 in Merritt (2013) for a more detailed treatment. Then in Section 3 we review various approaches used for studying the evolution of binary SMBHs and TDE rates, and justify our choice of the Monte Carlo method. Section 4 presents our initial conditions for single and binary SMBH simulations. The results of Monte Carlo simulations and the various aspects of the evolution of galactic nuclei are discussed in Section 5, and finally Section 6 summarizes our findings.

2 LOSS-CONE THEORY FOR SINGLE AND BINARY SMBH

The largest of the two values is called the loss-cone radius r_LC. A more rigorous treatment (e.g. Servin & Kesden 2017, fig. 1) shows that this simple approximation only slightly underestimates the exact relativistic solution (the largest error |${\sim } 20{{\ \rm per\ cent}}$| is at r_t = r_c, which for η = 1 corresponds to |$M_\bullet \simeq 1.5\times 10^7\, \mathrm{ M}_\odot$|⁠), and that stars can still be disrupted (and not captured) by SMBH up to ∼4.5 times more massive than this critical value. It is also common to formulate the loss-cone boundary in terms of the critical angular momentum |$L_\mathrm{LC}\equiv \sqrt{2\, G\, M_\bullet \, r_\mathrm{LC}}$|⁠.

The concept of loss cone may be extended to the case of a binary SMBH (e.g. Yu 2002; Milosavljević & Merritt 2003), describing the distance of closest approach to the binary centre of mass at which the incoming star undergoes a complex three-body scattering process and is ultimately ejected (if it was not captured by either of the two SMBHs in this process). Naturally, the critical distance is of order a – the semimajor axis of the binary, which is much larger than the loss-cone radius of an isolated black hole.

Regardless of the physical nature of the loss cone, the stars on the orbits with pericentre distances smaller than the critical distance are essentially lost within one dynamical time (either destroyed or ejected, in the latter case they may still return back at some later time, but this may well be viewed as an incomplete scattering event). The important question is whether the time-scale of repopulation of these loss-cone orbits is longer or shorter than the dynamical time. In the former case (empty loss cone), the rate of interactions between the stars and the SMBH(s) is supply-limited and depends on the mechanism of repopulation, while in the latter case (full loss cone), this repopulation is efficient enough that the rate no longer depends on it, and is simply proportional to the geometrical size of the loss cone.

In the (idealized) case of a spherically symmetric stellar system and the SMBH sitting in its centre, only the orbits with low enough angular momentum can enter the loss cone. Since in this case the angular momentum of a star can only be modified by two-body (collisional) relaxation, which is very inefficient in galactic nuclei with large enough number of stars (≫10⁶), the loss cone is typically quite empty. On the other hand, in non-spherical systems the angular momentum changes due to the torques from the global stellar potential (and not because of individual encounters), i.e. through collisionless processes. Even in this case, the loss cone needs not be full, but is typically at least partially refilled, so that the rate of interactions is a significant fraction of the ‘limiting’ full-loss-cone rate. A clear signature of this regime is that a star can enter the loss cone with any value of the angular momentum smaller than L_LC. By contrast, in the opposite (empty loss cone) regime, the typical change of angular momentum per one orbital period is small, hence a star is likely to enter the loss cone with the angular momentum only slightly smaller than the critical value L_LC. By studying the angular-momentum distribution of stars entering the loss cone, one may draw conclusions about the dominant mechanism of the loss-cone repopulation (collisional or collisionless). Collisional processes are believed to be unimportant in all but the densest galactic nuclei.

3 METHODS

The dynamical evolution of a stellar system surrounding a single or a binary SMBH may be studied using various approaches. N-body simulations offer the most straightforward way of incorporating all relevant physical processes (gravitational encounters of stars with the components of the binary and between stars themselves, stellar evolution, gas dynamics, etc.), but they are very computationally expensive. Broadly speaking, existing simulation methods fall into two categories (see Dehnen & Read 2011 for a review).

‘Collisional’ codes accurately compute the gravitational force, typically using the direct-summation approach (with the cost proportional to N²), and employ high-accuracy time integration schemes, often combined with regularization techniques to handle close encounters between particles, or tightly bound subsystems. Because of this, they cannot presently be used for systems with the number of particles significantly in excess of 10⁶, even when running on special-purpose hardware (Makino 1996) or modern GPU-accelerated parallel supercomputers (e.g. Wang et al. 2015; Li et al. 2017).¹ Real galaxies, of course, contain a much larger number of stars, so any simulation with a lower number of particles must be appropriately scaled, which is tricky when the various physical processes have different dependence on N. In particular, the two-body relaxation rate is |$\propto N^{-1}\, \log N$|⁠, being significantly higher in scaled-down N-body systems than in real galaxies. This distorts the interplay between relaxation-driven and collisionless repopulation of the loss cone and makes it difficult to draw firm conclusions about the long-term evolution of real galaxies (Vasiliev et al. 2015).

On the other hand, ‘collisionless’ codes sacrifice the accuracy of force computation and time integration to allow a significantly larger number of particles, reducing (but not eliminating) the impact of two-body relaxation. The gravitational force is computed approximately, using grid, tree-code, or fast-multipole methods, which scale as |$\mathcal {O}(N\, \log N)$| or even |$\mathcal {O}(N)$|⁠; the force is also softened at small distances, which additionally suppresses the two-body relaxation (but only moderately). We stress that these methods are called ‘collisionless’ only because they cannot properly handle two-body encounters (collisions), not because they somehow completely mitigate their influence: at a fixed number of particles, the two-body relaxation rate in these approaches is only a factor of few lower than in collisional codes (Hernquist & Barnes 1990). It is only the use of a substantially larger N that allows them to be nearly relaxation-free. However, in their standard form these methods are unsuitable for studying close encounters between stars and SMBH(s) by design.

It is not impossible to imagine a hybrid method that would combine a large number of particles and hence negligible two-body relaxation with an accurate treatment of close encounters between the particles and SMBHs. Indeed, the fast-multipole method can achieve the accuracy of force computation comparable to the direct-summation approach, while still retaining its |$\mathcal {O}(N)$| scaling (Dehnen 2014). Gualandris et al. (2017) used this method with up to N ∼ 10⁸ particles to study the rate of loss-cone repopulation in non-spherical galaxies, but their simulations did not actually contain SMBHs, because the leap-frog time integration scheme (standard for collisionless codes) and the need for gravitational softening prevented the possibility of accurately tracking the close encounters with SMBHs. Karl et al. (2015) and Rantala et al. (2017) introduced hybrid codes that combine the tree-code approach with regularization techniques suitable for accurate computation of SMBH orbits and stellar encounters with SMBHs. Although the force errors in the conventional tree-code schemes are too large for reliable analysis of loss-cone orbits, as demonstrated in the appendix of Gualandris et al. (2017), these approaches seem promising. A combination of the fast-multipole method and regularization techniques is conceivable, and would be the best method for studying the evolution of such stellar systems.

A very different approach to studying the dynamical interactions between stars and a binary SMBH is offered by scattering experiments, which determine the statistical properties of a large ensemble of independent three-body interactions (e.g. Quinlan 1996; Sesana, Haardt & Madau 2006, 2008; Chen et al. 2008). This information, augmented by suitable assumptions about the distribution of incoming stars, could be used to approximately follow the evolution of the binary and estimate the TDE rates (Sesana, Haardt & Madau 2007; Sesana 2010; Chen et al. 2011; Darbha et al. 2018). The main limitation of this approach is the lack of self-consistency regarding the distribution of stars: it either remains fixed or is evolved using simplified prescriptions, in contrast to a self-consistent evolution in N-body simulations (and, of course, in real galaxies).

In this work we use a hybrid Monte Carlo approach (Vasiliev 2015; Vasiliev et al. 2015) that combines the benefits of scattering experiments and N-body simulations, while largely avoiding their respective limitations. The evolving stellar system is still represented by N particles that do not correspond to individual stars, but rather sample their distribution function probabilistically. These particles move in the smooth gravitational potential, self-consistently computed from the ensemble of particles, and can interact with the binary SMBH.

Unlike conventional N-body simulations, the gravitational potential is expressed in terms of a spherical-harmonic expansion, centred on the binary centre of mass, with the radial dependence of each multipole term in this expansion described by a smooth (spline) curve with a relatively small (∼20−30) number of grid points in log r. This resembles the self-consistent field method of Hernquist & Ostriker (1992), with the difference that we use spline interpolation instead of a weighted combination of basis functions in radial coordinate. The ‘MEX’ method of Meiron et al. (2014) also computes the radial dependence of each multipole term directly from particle coordinates, without decomposing it into the basis-set expansion, but only at each particle’s position, without constructing a global smooth approximation for the potential. In our implementation, we first create a smooth approximation to the density profile from particle positions, which is also expressed in terms of a spherical-harmonic expansion with radially interpolated coefficients, and then solve the Poisson equation by numerically evaluating the 1D integrals in radius for each multipole term (for details see Vasiliev 2018, Section A.3.1).

The potential of the entire system is updated not at every time-step, as in ordinary N-body simulations, but much less frequently – at the end of each ‘episode’, which should be much shorter than the characteristic time-scale of the global evolution (e.g. the time needed to significantly shrink the binary), but still could be longer than the dynamical time, at least for the most bound particles in the central part of the system. Accordingly, the particles move in a smooth static stellar potential for the entire episode, but they do experience the time-dependent forces from the two SMBHs, which themselves follow a Keplerian orbit with fixed parameters during each episode. The time-dependent (non-monopole) part of the binary potential is only important when a particle is close enough to the centre. We record the changes in energy and angular momentum of all particles entering and exiting the sphere of radius r_enc = 5a (where a is the semimajor axis of the binary). At the end of each episode, these changes are summed up for all particles that had close encounters with the binary, and the binary orbit (a and eccentricity e) is adjusted by subtracting the equivalent amount of energy and angular momentum, and optionally taking into account the GW losses. Hence our approach may be viewed as a superposition of individual three-body scattering experiments, but with the parameters of the incoming stars drawn from their actual distribution function, which evolves in the course of the simulation.

Because each particle’s orbit is computed independently from the others, the computational cost scales as |$\mathcal {O}(N)$|⁠, and the method is trivially parallelized (we use the openmp approach to distribute the load to all cores of a single workstation). Despite this favourable scaling, the number of particles in our runs is typically ≲10⁶, similar to direct N-body simulations. However, unlike the latter, our method has a dramatically lower (∼100×, see fig. 1 in Vasiliev 2015) level of two-body relaxation for the same number of particles, due to several factors. First, the smooth global potential, with a low number of spherical-harmonic terms and a rather coarse radial grid, acts as a low-pass filter (similar to the use of a relatively large but spatially variable softening length, avoiding a loss of force resolution at small radii). Secondly, in updating the potential, we do not just use the instantaneous positions of particles at the end of each episode, but rather take N_samp ≫ 1 points sampled from each particle’s orbit during the entire episode, further suppressing the discreteness noise. Finally, the rather long update intervals for the potential act as a temporal smoother.

Since the loss-cone problems involve both collisional and collisionless processes, we also need to simulate the effect of two-body relaxation. This is achieved by adding perturbations to particle velocities at every time-step, with the magnitude determined by diffusion coefficients computed from the distribution function of all stars, which is also updated after each episode. Importantly, the amplitude of these perturbations can be assigned manually to mimic the level of relaxation expected in a stellar system with N_⋆ stars, which is unrelated to the actual number of particles N in the simulation. Various tests performed in previous studies demonstrate a good agreement between the Monte Carlo and N-body simulations of the same system, when the level of relaxation in the former approach is set to match the latter.

Quinlan & Hernquist (1997) and Hemsendorf, Sigurdsson & Spurzem (2002) developed hybrid codes combining the representation of the stellar potential in terms of a multipole expansion with a direct N-body simulation of a small number of particles in the very centre, including the two SMBHs. Our approach resembles these methods (especially in the first part), but with important differences (additional temporal smoothing and oversampling to reduce the level of noise, and the explicit account for two-body relaxation).

The Monte Carlo code raga, introduced in Vasiliev (2015), was used to study the TDE rates by single SMBHs in Vasiliev (2014) and Lezhnin & Vasiliev (2016), and applied to the binary SMBH evolution in Vasiliev et al. (2015). The new version of the code, used in this paper, is built on the same principles, but redesigned almost from scratch; it is more robust and computationally efficient, and can simultaneously deal with both the binary SMBH evolution and the captures of stars by each SMBH. The code is publicly available as part of the agama library for galactic dynamics (Vasiliev 2019).

The main advantages of the Monte Carlo approach in the context of dynamics of galactic nuclei are

• correct balance between collisional and collisionless effects (unlike N-body simulations with a scaled-down number of particles);

• self-consistent treatment of the evolution of the stellar distribution and the loss-cone repopulation (unlike scattering experiments).

• possibility of using physically correct size of the loss cone [unlike Li et al. (2017), who had to rely on extrapolation].

• Moderate computational cost – simulations with 10⁶ particles run for a significant fraction of the Hubble time only take about a day on a conventional 16-core workstation.

The limitations of the method are as follows:

• It can only be applied to systems with a well-defined centre and an already formed hard SMBH binary, so we cannot reliably simulate the early stages of evolution, when the TDE rates briefly reach a peak (e.g. Chen et al. 2009; Li et al. 2017).

• We only follow the evolution of binary semimajor axis and eccentricity, and neglect the orientation; however, the latter is not expected to change dramatically, unless the stellar cluster is strongly rotating and the binary orbital plane is misaligned with it (Rasskazov & Merritt 2017).

• We also assume that the binary resides in the galaxy centre, neglecting the Brownian motion. Bortolas et al. (2016) found that allowing the SMBH to wander has little effect on the evolution of systems where the loss cone is replenished mainly by collisionless processes. On the other hand, Holley-Bockelmann & Khan (2015) and Mirza et al. (2017) discovered a regime where the SMBH centre of mass exhibits coherent motion around the galactic nucleus, which occurs when the stars in the nucleus are on orbits corotating with the binary; however, this does not seem to significantly affect the hardening rate.

• Our approach for simulating two-body relaxation does not account for its possible enhancement due to resonant relaxation (Rauch &Tremaine 1996); however, we do not expect it to play any role in highly chaotic three-body interactions.

• Stellar binaries are absent in our simulations; in reality they may be an important source of hypervelocity stars, as predicted by Hills (1988) for single SMBH–binary star systems and by Wang et al. (2018) for binary SMBH–binary star systems.

• We do not explore the effect of mass segregation (all particles have the same dynamical mass); however, any pre-existing segregation is expected to be erased during the merger (Gualandris & Merritt 2012), and is unlikely to be re-generated within the lifetime of the binary (although this is difficult to avoid in N-body simulations, see e.g. Khan, Berczik & Just 2018).

4 VARIANTS OF MODELS

We explore a large suite of models, varying the parameters of the binary SMBH and the stellar system. We only consider ‘isolated’ galaxies, initially constructed to be in equilibrium, not the merger products: Vasiliev et al. (2015) and Gualandris et al. (2017) demonstrated that the long-term evolution of the binary depends qualitatively on the geometry of the stellar potential (spherical, axisymmetric, or triaxial), but only moderately on its particular properties such as axial ratios (Bortolas et al. 2018b), and is similar between isolated and merger simulations. Therefore, we examine three series of models, with initial profiles following the Dehnen (1993) model having the inner slope γ = 1 and constant axial ratios – 1|$\colon$|1|$\colon$|1 for spherical, 1|$\colon$|1|$\colon$|0.8 for axisymmetric and 1|$\colon$|0.9|$\colon$|0.8 for triaxial models. These models are designed to be in a self-consistent equilibrium with a single SMBH in the centre, which has a mass 0.01 times the total stellar mass. The initial N-body snapshots are prepared with the Schwarzschild orbit-superposition code smile (Vasiliev 2013) and contain 0.5 × 10⁶ particles. Next we take 25 per cent of particles with lowest values of angular momentum, and replace each one with five identical particles of a correspondingly smaller mass, bringing the total number of particles to N = 10⁶. This static mass refinement scheme (similar to the one in Lezhnin & Vasiliev 2016) improves the statistics of TDE rates, since the more numerous smaller particles better sample the underlying distribution function (due to independent velocity perturbations, they quickly spread out in the phase space). Thus our mass resolution is equivalent to a 2.5 × 10⁶ particle snapshot in the low-angular-momentum part of the phase space, which is responsible for the interaction with the SMBH binary and for the TDE rates. The presence of different mass groups does not lead to mass segregation, since we use the same the drift coefficient in the treatment of two-body relaxation for all particles (i.e. particle have different tracer masses but identical dynamical masses).

The loss-cone radius for each SMBH is set by equation (1); particles reaching a smaller distance to the SMBH are captured and removed from the simulation, with their mass added to the respective SMBH mass (although the total added mass is always significantly smaller than the initial SMBH mass). The relaxation level is set by the number of stars in the physical system being modelled, using the value of Coulomb logarithm ln Λ = ln (0.3M_•/M_⊙). We stress that the number of particles in the simulation is always the same (10⁶), but we quote all physical values (e.g. the galaxy mass) according to the chosen dimensional scaling factors. The simulations are run for several thousand N-body time units, corresponding to a physical time of up to 3 Gyr.

To compare the TDE rates between galactic nuclei hosting single and binary SMBHs, we also conducted simulations with single SMBHs. The complication in such a comparison is that a binary SMBH disrupts the pre-existing stellar cusp at the early stage of its evolution, resulting in a shallower density profile. To account for this, we started the single SMBH runs from the snapshots of corresponding simulations with the binary SMBH after the initial stage of relatively fast evolution (∼10⁷ yr). After this early stage, the density profile changes very little in the simulations with binary SMBHs, so the comparison is more adequate.

Given the large number of parameters (total mass of the binary M_•, mass ratio q, initial eccentricity e, three variants of geometry, inclusion, or absence of captures), we do not discuss all possible combinations of them, but instead focus on the main trends arising from changing one parameter at a time. In general, we believe the triaxial geometry to be the most relevant physically, given that a realistic galaxy merger would not produce a precisely axisymmetric (much less spherical) system, but we include the other two cases for comparison.

5 RESULTS

5.1 Hardening rate of the binary

where the integration is carried over the range of stellar energies where the stars are expected to efficiently repopulate the loss cone. One may take the value of the stellar potential at origin Φ_⋆(r = 0) as the lower limit of integration, which is equivalent to considering only the stars that are unbound to the SMBH (but of course still bound to the entire galaxy). These two definitions typically agree to within a factor ≲1.5.

Of course, the stellar density profile itself evolves with time (see Section 5.6). The early stage of binary formation and hardening leads to a rapid erosion of the pre-existing stellar cusp (e.g. Milosavljević & Merritt 2001), and afterwards the density continues to decrease, but very gradually. Most of the stars that interact with the binary at the later stage are scattered into the loss cone from relatively large radii, hence their disappearance has little effect on the density profile. Accordingly, we compute the value of the ‘full-loss-cone’ hardening rate right after the formation of the hard binary, and neglect its subsequent evolution.

Fig. 1 shows the evolution of 1/a for models with M_• = 10⁸ M_⊙ (similar results were obtained for other M_• values) and different q. The hardening rate is nearly independent of the eccentricity (the figure shows the runs with e = 0) and is insensitive to whether the stars are captured by individual SMBHs or not. The likely explanation is that the geometric size of the loss-cone of a single SMBH is much smaller than the cross-section for strong three-body interactions with the binary, therefore the probability of being captured during a single three-body interaction is rather minor, and most stars are ejected and contribute to the binary hardening.

Figure 1.

Evolution of inverse semimajor axis of the binary 1/a with total mass M_• = 10⁸ M_⊙ and mass ratio q = 1 (blue), q = 1/3 (red), q = 1/9 (cyan), q = 1/27 (magenta), and q = 1/81 (orange). Solid lines show the triaxial systems, which happily continue hardening, and dashed – axisymmetric systems, which virtually stall after a short period of initial hardening, due to the depletion of the loss cone (spherical models, not shown on this plot, stall even sooner). Dot–dashed lines show equivalent simulations with GW emission turned on, which reach coalescence in ∼2−2.5 Gyr in the case of zero eccentricity; for a higher initial eccentricity this would occur sooner. Dashed black line illustrates the hardening rate expected in the initial model if the loss cone were full (equation 4), which is clearly much higher than the actual rate measured in simulations. This reference full-loss-cone rate is somewhat reduced at later stages due to a partial destruction of the cusp, but the actual hardening rate still remains considerably lower because of the depletion of centrophilic orbits, which is less prominent for q ≪ 1.

Open in new tabDownload slide

For spherical and axisymmetric systems, the evolution of 1/a slows down dramatically after the initial phase, and a barely reaches |$0.1\, a_\mathrm{h}$|⁠, which is still far larger than a_GW. Therefore, unless e₀ is close to unity or other non-stellar-dynamical processes are considered, the binary SMBHs would not merge in these cases, in agreement with Vasiliev et al. (2015); Gualandris et al. (2017). Using direct N-body simulations, Khan & Holley-Bockelmann (2013) argued that the hardening rate remains sufficiently high even in axisymmetric systems to reach coalescence within the Hubble time, which seems to contradict our findings. However, their different conclusion is likely to result from the enhanced two-body relaxation in the simulations with N ∼ 10⁶ compared to real galaxies with N_⋆ ∼ 10⁸−10¹⁰, as discussed earlier and corroborated by the analysis of N-dependence of the hardening rate in the two papers quoted above.

5.2 Eccentricity evolution

The eccentricity of the binary changes rather moderately in the case of a (nearly-)isotropic distribution of field stars (Fig. 2). For models with a rather high initial value (e ≳ 0.6), it tends to increase further, reaching e ≳ 0.98 for an initial value e = 0.9, which, of course, substantially shortens the lifetime of the binary before the GW-induced coalescence. On the other hand, for lower initial values the eccentricity remains roughly constant. The inclusion or neglect of captures has little influence on the eccentricity evolution. However, the geometry of the stellar potential seems to play some role, in particular, the eccentricity grows faster in spherical than in the triaxial cases, though it is clear from Fig. 1 that the degrees of binary hardness in different geometries are different.

Figure 2.

Eccentricity evolution over time for simulations with initial eccentricities equal to e = 0.01 (red), 0.3 (green), 0.6 (blue), 0.9 (black) in a triaxial potential with M_• = 10⁹ M_⊙ and mass ratios q = 1 (upper panel) and q = 1/27 (lower panel). We also demonstrate how the eccentricity evolves in the spherical potential (dot–dashed line). Dashed lines correspond to simulations without captures. In all cases with e > 0.3 the eccentricity growth and saturation is evident. However, the overall evolution of e(t) is quite stochastic, especially for q ≪ 1.

Open in new tabDownload slide

Observationally, there is growing evidence that nuclear star clusters have a substantial amount of rotation (e.g. Seth et al. 2008; Schödel, Merritt & Eckart 2009; Hartmann et al. 2011). This is important for the eccentricity evolution, since the scattered stars coming from prograde (retrograde) orbits tend to increase (decrease) the angular momentum of the binary. Iwasawa et al. (2011) and Sesana, Gualandris & and Dotti (2011) studied the evolution of very unequal mass ratio binaries (q ∼ 10⁻²) and found that the eccentricity strongly grows or declines when the stellar cluster is counterrotating or corotating with the binary orbit. Holley-Bockelmann & Khan (2015) and Mirza et al. (2017) confirmed these trends with higher-resolution N-body simulations for a wider range of parameters, and additionally found a regime where the binary centre of mass settles into a nearly circular orbit with a radius ∼r_infl in a corotating stellar cluster. The latter study also considered the evolution of binary orbital plane orientation in their simulations, while Rasskazov & Merritt (2017) addressed this question using Fokker–Planck formalism.

To explore the impact of the net rotation of the stellar cluster, we conducted an auxiliary set of simulations with the initial distribution of stars modified in the following way. To bias the stellar distribution towards corotation (counterrotation), we change the sign of the velocity for particles with negative (positive) |$z$|-component of angular momentum, with a probability ranging from 0 to 100 per cent. This preserves the orbital structure of the stellar system (required for it to remain in a self-consistent equilibrium) and a nearly-isotropic distribution of orbit eccentricities (an assumption used in computing the two-body relaxation rate in the Monte Carlo approach).

Fig. 3 shows the evolution of binary eccentricity in a range of triaxial M_• = 10⁹, q = 1 or 1/9 systems with different degree of rotation, starting from the initial value e = 0.3. In agreement with previous studies, the eccentricity rapidly decreases or increases in corotating or counterrotating clusters, correspondingly. However, this stage lasts only around 10⁷ yr, after which the eccentricity reaches a nearly stationary value for q = 1, or keeps increasing only moderately for q = 1/9. The likely reason is that the loss cone is repopulated mainly from centrophilic orbits on longer time-scales, and they do not conserve the sign of angular momentum, hence the orbits of incoming stars are eventually isotropized.

Figure 3.

Eccentricity evolution over time for simulations with the initial eccentricity equal to e = 0.3, mass ratios q = 1 (upper panel) and q = 1/9 (lower panel), and various amount of rotation of the stellar cluster. Solid line represents an isotropic cluster in L_|$z$|, dashed – counterrotating clusters with specific fractions of the co- and counterrotating orbits, dot–dashed – corotating clusters. After the short initial phase of rapid eccentricity growth (in counterrotating clusters) or decay (in corotating ones), the subsequent evolution is much more gradual in both q = 1 and q = 1/9 cases, although the long-term growth of eccentricity is more prominent for q = 1/9.

Open in new tabDownload slide

In the absence of GW losses, the hardening rates are not significantly affected by eccentricity. The capture rate is moderately (less than a factor of two) suppressed for corotating systems. The inclusion of captures slightly enhances the increase of eccentricity in the majority of our simulations, as predicted by Sesana et al. (2008) and Chen et al. (2011). In the case of very unequal-mass mergers, the eccentricity at the time of the formation of a hard binary is expected to be higher, having grown at the previous stage of dynamical friction (Dosopoulou & Antonini 2017).

5.3 Capture rates

We now discuss the dependence of capture rates on the parameters of the binary SMBH, and compare them with the rates expected in galactic nuclei hosting a single SMBH. As explained in Section 2, black holes more massive than a few|${}\times 10^7\, \mathrm{ M}_\odot$| swallow most stars entirely (except very extended giants), without producing an observable TDE flare. We refer to all cases when a star passes at r < r_LC (equation 1) as captures, regardless of the outcome.

Fig. 4 shows the time-dependent capture rates obtained in our simulations of spherical, axisymmetric, and triaxial clusters, for three different SMBH masses (10⁷, 10⁸, and 10⁹ M_⊙) and two values of mass ratio (q = 1 and q = 1/27). For comparison, we also plot the average rates in equivalent simulations with a single SMBH, started from the initial configuration shortly after the formation of a hard binary (when the original density cusp has already been destroyed, but the population of loss-cone orbits has not yet been depleted). These ‘reference rates’ are typically around few |${}\times 10^{-5}~\mathrm{ M}_\odot \, \mbox{yr}^{-1}$| in spherical systems (e.g. Wang & Merritt 2004; Stone & Metzger 2016), but higher in non-spherical cases (e.g. fig. 4 in Vasiliev 2014 or fig. 5 in Lezhnin & Vasiliev 2016).

Figure 4.

Time-dependent capture rates in simulations of binary SMBH in various geometries: triaxial (left-hand panel), axisymmetric (centre), spherical (right). Colour corresponds to the SMBH mass, from top to bottom: 10⁹ (blue), 10⁸ (green), and 10⁷ M_⊙ (red); solid lines show the equal-mass case (q = 1), and dashed – for the mass ratio q = 1/27. Average capture rates in otherwise equivalent galaxies with a single SMBH are shown by symbols in the right side of each panel (they are approximately constant over time). Since the density cusp is less depleted in unequal-mass cases, the rates are higher for q ≪ 1 than for q = 1. In simulations with binary SMBH, the rates are initially higher, but decline over time (faster in more symmetric potentials – in a spherical geometry they are essentially zero after a fraction of Gyr, note the different scale of the horizontal axis in the right-hand panel).

Open in new tabDownload slide

As seen from the right-hand panel (spherical case), the capture rates in simulations with binaries are initially much higher than the steady-state rates in galaxies with a single SMBH, but they rapidly decline with time. We caution that this early transient period is rather sensitive to the particular way of setting up the initial conditions, so one should not compare these values quantitatively, but the difference of more than an order of magnitude is a robust feature. This enhancement has been found previously in several studies (Chen et al. 2011; Wegg & Bode 2011; Li et al. 2015), but lasts only a short time (≲10⁷ yr).

The initial transient period is similar in all three geometries, but the late-time evolution is rather different. In the spherical case, both the capture rates and the binary hardening rates drop essentially to zero after this initial period, because the loss cone is nearly empty. In the axisymmetric case the situation is similar but less dramatic; nevertheless, the capture rates drop by more than an order of magnitude compared to the single-SMBH case. In the triaxial case, however, the long-term capture rates are only moderately smaller or even comparable to those found in single-SMBH systems, scaling nearly linearly with the SMBH mass. This proportionality is an indication of the major role of collisionless loss-cone refilling mechanism, whose rate is independent of the number of stars (recall that the geometric size of the loss cone is proportional to M_• for direct-capture events).

The capture rates are insensitive to the binary eccentricity, varying by less than 30 per cent. In the equal-mass binary (q = 1), each SMBH captures roughly half of the stars, but the proportion of stars captured by the primary (more massive) component of the binary tends to 100 per cent very quickly as the mass ratio q decreases (even in the case q = 1/3, this fraction is |${\gtrsim }90{{\ \rm per\ cent}}$|⁠). On the other hand, since a massive enough primary SMBH would swallow stars entirely, the secondary remains the only source of observable TDE (Fragione & Leigh 2018).

Recently Darbha et al. (2018) studied the rates of TDE by binary SMBHs, using scattering experiments (drawing from a uniform distribution of incoming stars in angular momentum) under the assumption of a full loss cone. They found that the rates are increased compared to the case of a single SMBH by a factor of few, weakly depending on the mass ratio (as long as it is higher than few × 10⁻²). The likely explanation is that the time-dependent potential of the binary perturbs the incoming orbits and forces a larger fraction of them to enter the loss cone of either of the two SMBH, especially in the case of chaotic resonant scattering events with multiple pericentre passages. We ran a series of isolated scattering experiments, using a similar setup to Darbha et al. (2018), and confirmed their conclusions regarding the enhancement of capture rates in binary systems.

At a first glance, this seems to contradict our findings that the capture rates by binaries are smaller than or comparable to those by single SMBHs in simulations of triaxial systems, where the loss cone is at least partially full (the most relevant case for comparison). However, the analysis of time-dependent capture rates in our simulations indicates a different behaviour in single- and binary-hosting nuclei. At the early stage of evolution, the rates are considerably higher for systems with binary SMBHs, but subsequently they decrease; their time evolution is similar to that of the hardening rates, and the slowdown is due to the loss-cone depletion in both cases. By contrast, in systems with a single SMBH they stay at a roughly constant level, or even increase with time, as the phase-space gap at low angular momenta, carved out by a pre-existing binary, gradually fills up (see Merritt & Wang 2005; Lezhnin & Vasiliev 2015, 2016). In the case of a continuously evolving binary, this gap is not refilled but rather deepens with time, because most centrophilic orbits that still exist in the system are scattered away by the binary, rather than find their way into the loss cone of a single SMBH. The loss cone of the binary is gradually emptied at any given energy, progressing from inside out to less bound orbits (see the discussion in section 4.3 and fig. 6 in Vasiliev et al. 2015). The fact that we find comparable capture rates in triaxial single- and binary-hosting nuclei appears to be a coincidence arising from the cancellation of opposite trends – enhancement of captures by binaries compared to single SMBH for a fixed supply rate of incoming stars, and the decrease of this rate with time in galactic nuclei where the binary continues to evolve.

Darbha et al. (2018) also found that the fraction of stars captured by the secondary SMBH scales roughly linearly with the mass ratio q, whereas we observe a much steeper drop in the secondary disruptions as q decreases. The difference likely arises from a shallower dependence of capture radius on the mass of a single SMBH (equation 1) for less massive SMBH (⁠|${\lesssim } 10^6\, \mathrm{ M}_\odot$|⁠) studied in their work.

5.4 Loss-cone orbits

We analyse the properties of stars that interact with the binary in the following way. For each particle entering or leaving the sphere of radius 5a from the binary centre of mass, we record the initial and final values of energy and angular momentum, and also the direction of the velocity on exit. Particles that leave the interaction zone with a positive total energy are ejected from the galaxy as hypervelocity stars (see next section); otherwise a particle may return later and experience a ‘secondary slingshot’ (Milosavljević & Merritt 2003). Since a three-body interaction may last a long time and consist of multiple close approaches with the binary, we merge the data recorded for successive interactions of the same particle if the time between them is ≤1 orbital period at the given energy. The majority of particles experienced more than one interaction, and a substantial fraction of interactions with the initial value of angular momentum L² ≲ 2GM_•a resulted in the ejection of the particle with a positive total energy (Fig. 5).

Figure 5.

Properties of stars in the loss cone of the binary SMBH. Shown is the distribution of stars in squared angular momentum: if a star were orbiting a single SMBH instead of a binary, this quantity would be a proxy for the pericentre distance. Different lines show the outcome of the interaction: blue are all stars that have been scattered by the binary, but remained in the galaxy; green are hypervelocity stars that were ejected from the galaxy with positive total energy; and red are stars that were captured by one of the SMBH. Solid lines are for the triaxial case and dashed – for the axisymmetric case; top panel is for equal-mass binary, and bottom – for the mass ratio q = 1/27.

Open in new tabDownload slide

In triaxial systems, the distribution of initial values of squared angular momentum (dN/dL²) is approximately flat, which is a characteristic sign of a full-loss-cone regime. Note that this does not mean that the interaction rate corresponds to the full-loss-cone rate (4): rather, the stars close to the loss-cone boundary are equally likely to cross it with any value of angular momentum, not just the value only slightly smaller than the boundary (as would happen in the empty-loss-cone regime). However, in more symmetric systems the situation is different: the dashed lines illustrate that in the axisymmetric case, the probability distribution increases as one moves outwards from the loss-cone boundary (to the right in that plot). The distribution of captured stars (red lines) is concentrated at low L² in unequal-mass cases (q ≪ 1), because the more massive component of the binary resides closer to its centre of mass.

The orbits of stars in a time-dependent potential cannot be characterized rigorously (but see Li, Holley-Bockelmann & Khan F. 2018 for an attempt); however, it is likely that most of the orbits that bring stars into the vicinity of the binary are centrophilic, i.e. can attain arbitrarily small values of angular momentum (e.g. boxes, pyramids, some minor resonant families, or chaotic orbits fulfil these criteria). At least in the case of a single SMBH in non-spherical system, the majority of captured stars arrive from these orbits (fig. 5 in Vasiliev 2014), and we expect this to hold also for the loss cone of a binary (see section 4.2 in Vasiliev et al. 2015).

5.5 Hypervelocity stars

The stars involved in the three-body interaction events involving a hard binary are accelerated to high velocities (Yu & Tremaine 2003; Sesana et al. 2006; Darbha et al. 2019), occasionally exceeding the escape velocity from the galaxy centre (∼10³ km s⁻¹) by large factors. These stars would quickly leave the host galaxy and travel in the intergalactic space with large enough speeds that their proper motion could be detectable even at large distances (Guillochon & Loeb 2015).

Fig. 6 shows the velocity distribution function of stars after 1 Gyr of evolution, illustrating the emergence of an unbound population of stars reaching velocities up to 10⁴ km s⁻¹. Their total number is proportional to the mass of the binary, with a weak dependence on the mass ratio, as explained in the next section. The shape of the velocity profile is roughly the same in all cases and is close to exponential, with a characteristic scale 1000−1500 km s⁻¹ comparable to the orbital speed of the binary during most of its evolution. In the case of eccentric binary, a small fraction of stars happens to be ejected while the binary is at the pericentre of its orbit, acquiring much higher speeds (dubbed ‘semirelativistic stars’ by Guillochon & Loeb 2015).

Figure 6.

Distribution of stars in velocity after 1 Gyr of evolution, for the case |$M_\bullet =10^8\, \mathrm{ M}_\odot$|⁠, e = 0 and various values of mass ratio: q = 1 (blue), q = 1/3 (red), q = 1/9 (cyan), q = 1/27 (magenta), q = 1/81 (orange). In the top panel, all stars are shown (both bound and unbound), and the additional black line corresponds to the initial model. In the bottom panel, only the unbound stars are shown, which roughly follow an exponential distribution; additionally, the case of equal-mass binary on an initially eccentric orbit (e = 0.9) is shown in dashes, illustrating a longer tail at high ejection velocities.

Open in new tabDownload slide

Fig. 7 illustrates the angular distribution of hypervelocity stars with |$v$| > 1300 km s⁻¹, in an equal-area projection. The binary orbit lies in the equatorial plane, and the semimajor axis is directed at longitude 0°. Most of the stars are ejected in directions close to the orbital plane of the binary, in agreement with Zier & Biermann (2001) and Sesana et al. (2006). The azimuthal distribution is naturally uniform in the case of small eccentricity, but becomes more biased towards the direction of the binary semimajor axis in the high-eccentricity case (peaks at longitudes 0° and 180°). This appears to disagree with the statement made in Sesana et al. (2006); however, they do not plot the angular distribution itself, only quote its moments (mean azimuthal angle and its dispersion), and these numbers are in good agreement with our results. Recently Rasskazov et al. (2018) independently found an enhanced proportion of ejected stars near the orbital plane and along the direction of the binary semimajor axis (their fig. 4).

Figure 7.

Angular distribution of hypervelocity stars for mass ratios q = 1 (top row), q = 1/3 (middle row), q = 1/9 (bottom row) and initial eccentricities e = 0.01 (left column) or e = 0.9 (right column). The majority of stars are directed close to the orbital plane of the binary, and in the case of an eccentric binary, along its semimajor axis.

Open in new tabDownload slide

To verify our calculations, we also performed a suite of scattering experiments with isolated binary SMBH interacting with an isotropically distributed population of incoming stars; the angular distribution of ejected stars was similar to the one obtained in the fully self-consistent simulations of binary evolution. In our approach, the orientation of the binary semimajor axis is assumed to be fixed, but in reality it precesses in the orbital plane (even when the plane itself is largely conserved), hence the azimuthal distribution is averaged over time to a nearly uniform one. Finally, we note that the angular distribution is nearly independent of the radius at which it is measured, indicating that the triaxial galactic potential is unable to significantly deflect these high-velocity stars.

5.6 Evolution of the stellar distribution

Stars that have been scattered by the binary are ejected from the galaxy, leaving a ‘mass deficit’ in its central part. The depleted density profile in the centre has long been considered a signature of the binary being present in the galaxy currently or at some time in the past (e.g. Milosavljević et al. 2002; Volonteri, Madau & Haardt 2003; Graham 2004; Merritt 2006). While we confirm that some depletion takes place, it does not appear to have a pronounced and unambiguously detectable effect on the density profile.

The central mass deficit is usually defined as the difference between the actual surface density profile Σ(R) and the reconstructed initial profile before the binary evolution had taken place. The trouble is, of course, that this initial profile is not known, and the inference about missing mass strongly depends on the assumptions about the original un-depleted profile. For instance, Milosavljević et al. (2002) deprojected Σ(R) of a sample of galaxies to obtain the intrinsic density profiles ρ(r), then for each galaxy located the radius r₂ where the negative logarithmic slope γ ≡ −dρ/dr attains the value 2, and extrapolated the initial profile inwards as |$\rho _\mathrm{init} = \rho (r_2)\, (r/r_2)^{-2}$|⁠. This assumption of an initial steep cusp is adequate for low-mass galaxies, but if the actual density was shallower, this procedure will overestimate the amount of missing mass; indeed their inferred mass deficits change by a factor of few when the inner slope is reduced to 3/2. Moreover, there is no fundamental reason why the initial density profile would follow a power law with any particular slope.

Graham (2004) assumes that the un-depleted galaxy follows the Sérsic density profile, and the deviations from it may be parametrized by a core-Sérsic model (equation 5 in Graham et al. 2003). Thus the surface density profile is fitted by the core-Sérsic model, and its structural parameters (Sérsic index n and effective radius R_e) are assumed to characterize the un-depleted profile. He argues that the Sérsic model provides a smoothly varying logarithmic slope, and is hence able to describe Σ(R) of some galaxies without the need to introduce a core (fig. 3b in Graham 2004). However, the (non-cored) Sérsic model implies a rigid link between the density slopes in the central parts and in the outskirts of the galaxy: the larger is the parameter n, the steeper is the central profile and the shallower is the outer profile. In the core-Sérsic model, this relation is removed, and the index n is mainly determined by the outer profile, while the size of the core and the steepness of the transition are determined by the density profile in central parts. So far there is nothing wrong about such fitting procedure; however, when the mass deficit is inferred as the integrated difference between the cored and non-cored profiles with the same n and R_e, it depends strongly on the parameter n, i.e. on the properties of the galaxy in the outer parts, which have no physical relation to the region around the binary. Fig. 1 in Hopkins & Hernquist (2010) illustrates this point: if a galaxy with an initial Sérsic profile experiences a merger (without any secondary SMBH), the accreted material in the outer parts makes the density profile shallower, and hence increases the Sérsic index. Even though there is no change in the actual density profile in the inner part, the core-Sérsic model would imply a significant mass deficit w.r.t. a steeper cusp resulting from the extrapolation of the outer profile inwards.

Fig. 8 shows the evolution of surface density profile in a simulation of a triaxial galaxy with M_• = 10⁸ and q = 1/9. After the formation of a hard binary (semimajor axis a ≃ a_h, equation 2), the density is somewhat reduced within R ≲ r_infl), and after 1 Gyr the semimajor axis shrinks by another factor of 30, and the density is further depleted within a region a few times larger. However, there is no dramatic difference between density profiles at any stage of evolution. The reason is that the stars interacting with the binary are supplied from relatively large radii by centrophilic orbits, and hence their depletion has a very weak effect (⁠|${\lesssim } 10{{\ \rm per\ cent}}$|⁠) on the density profile at these radii. For comparison, we show the reconstructed ‘un-depleted’ profiles using either the γ = 2 inward extrapolation or the core-Sérsic fit: both approaches grossly misrepresent the actual initial profile. The inferred mass deficit is actually not too far off when using the γ = 2 extrapolation, because an overestimate of the initial density at R ≲ r_infl is partly compensated by a slight underestimate at R ≳ r_infl, while the core-Sérsic fit implies a 10 × larger mass deficit than the true value. Unsurprisingly, fitting a core-Sérsic model to the initial profile also implies a presence of a core with roughly the same amount of missing mass. Moreover, varying the strength of the transition between the core and the outer Sérsic profile leads to a few-fold change in the estimated mass deficit, even though the variation of the density between the fitted profiles is only of order 1 per cent.

Figure 8.

Surface density profiles of triaxial models with M_• = 10⁸, q = 1/9, at different stages of evolution: initial (Dehnen γ = 1 profile, blue), shortly after the formation of a hard binary (semimajor axis a ≃ 1 pc, green), and after 10⁹ yr of evolution (a ≃ 0.03 pc, red). For the latter profile, we additionally show two models which aim at reconstructing the ‘original’ (non-depleted) density under various assumptions. Dashed magenta line shows the density profile continued as ρ ∝ r⁻² inwards from R ≃ 100 pc (the radius at which the logarithmic slope of the 3D density profile equals −2, as in Milosavljević et al. 2002). Dot–dashed cyan line shows the Sérsic profile with shape parameter n ≃ 4.5, which is extrapolated from the outer part of the actual density profile. More precisely, a core–Sérsic model is fitted to the actual profile (it matches the red curve reasonably well), and then the core is undone, as in Dullo & Graham (2012, 2014) and Rusli et al. (2013). Both these extrapolated reconstructions bear little resemblance to the actual initial profile, and they yield overestimated mass deficits (⁠|${\sim } 4\, M_\bullet$| for γ = 2 and |${\sim } 18\, M_\bullet$| for the Sérsic models), while the actual missing mass is |${\sim } 0.4\, M_\bullet$| at t = 10⁷ yr and |${\sim } 2\, M_\bullet$| at t = 10⁹ yr.

Open in new tabDownload slide

The bottom line is that estimating the mass evacuated from the central part of the galaxy based on its present-day density profile is nearly impossible. Not only the original density profile is unknowable, but also the mass deficit is highly sensitive to the behaviour of the fitted profile near the transition radius. Even though this could not be done for any individual galaxy because of large scatter (e.g. table 4 in Rusli et al. 2013), nevertheless the general trends could be explored by analysing a large ensemble of galaxies. Hopkins & Hernquist (2010) assess the typical amount of missing mass by comparing the averaged profiles of cored and non-cored galaxies, and infer it to be comparable to M_• (note however that M_• in their list of galaxies is not measured directly but estimated from existing relations with other galaxy structural parameters).

On the theory side, Merritt (2006) explored the evolution of unequal-mass binaries (1/40 ≤ q ≤ 1/2) in spherical galaxies with initial density profiles following the Dehnen model with γ = 1/2, 1, 3/2, using direct-summation N-body simulations with N ≤ 2 × 10⁵. In a spherical system, the binary evolution is expected to stall at a semimajor axis a_stall that is only moderately smaller than a_h. He computes the mass deficit as the difference between the (known) initial density profile and the one at the moment when a = a_stall. He finds that the mass deficit roughly equals half the total mass of the binary, and weakly depends on the mass ratio: |$M_\mathrm{def} \simeq 0.7\, M_\bullet \, q^{0.2}$|⁠. If the binary continues to shrink due to the interaction with stars (and not because of other mechanisms such as gas drag or GWs), the evacuated mass is expected to increase, roughly proportional to |$J\, M_\bullet \, \ln [a_\mathrm{init}/a(t)]$|⁠, where J ≃ 0.5 is the dimensionless mass ejection rate coefficient (Quinlan 1996; Sesana et al. 2006).

Fig. 9 shows the evolution of evacuated mass M_evac in our simulations of triaxial galaxies with various values of M_• and q, confirming a nearly-linear growth of M_evac with ln 1/a. In this plot, M_evac is computed as the total mass of particles with positive energy (ejected from the entire galaxy). This is likely an underestimate of the mass deficit, because not all stars scattered by the binary acquire velocities greater than the escape speed (although most of them do at later stages, when a ≪ a_h). Another approach is to compute the difference between the initial density profile and the one at the given time (e.g. Khan et al. 2012, section 4), which produces up to ∼2 × higher estimate of the ejected mass, in agreement with a similar analysis of Milosavljević & Merritt (2001, their equation 41). We note in passing that in triaxial systems, stars ejected with less-than-escape speed are unlikely to return back to the centre and experience a secondary slingshot, being deflected by the large-scale torques in the potential.

Figure 9.

Evolution of evacuated mass as a function of inverse binary semimajor axis. Solid lines correspond to the total mass of the binary |$M_\bullet = 10^7\, \mathrm{ M}_\odot$|⁠, circles – |$10^8\, \mathrm{ M}_\odot$|⁠, crosses – |$10^9\, \mathrm{ M}_\odot$|⁠; different colours correspond to mass ratio q = 1 (red), 1/3 (green), 1/9 (blue), and 1/27 (black). Total evacuated mass at the end of simulation (3 Gyr of evolution) in triaxial case scales with q as |$M_{\rm evac} \approx M_\bullet \, q^{1/4}$|⁠. Dashed line represents the simulation with |$M_\bullet =10^7\, \mathrm{ M}_\odot$| and q = 1 in an axisymmetric stellar potential, in which case the evolution of semimajor axis does not progress too far; the case of a spherical potential is not shown, as the evolution is negligible.

Open in new tabDownload slide

Since the stars ejected after interacting with the binary arrive from nearly-radial orbits, their disappearance creates a tangential velocity anisotropy in the central parts of the galaxy. Fig. 10 shows the radial profile of the anisotropy coefficient, |$\beta \equiv 1 - \sigma _t^2/(2\sigma _r^2)$|⁠, where σ_r and σ_t are the radial and tangential velocity dispersions. The initial models were constructed to be nearly-isotropic, but the depletion of radial orbits leaves a pronounced trace in the anisotropy profile in the same range of radii where the density profile is modified (up to a few influence radii). The magnitude of tangential anisotropy is in qualitative agreement with the simulations of Quinlan & Hernquist (1997, their fig. 7), Milosavljević & Merritt (2001, their fig. 16), and Rantala et al. (2018, their fig. 7). We stress that this velocity anisotropy coefficient is not directly observable, but has to be inferred from dynamical models (at least when its value is not assumed or kept constant in the modelling procedure, as often done in the Jeans approach). In fact, it is hard to find any direct signature of core depletion in the observable values – the most apparent one is the lack of a central spike in the line-of-sight velocity dispersion profile in tangentially anisotropic systems. Using the Schwarzschild method, Thomas et al. (2014) did find an indication for tangential anisotropy in central parts of several galaxies with cored density profiles, which may be regarded as a more substantial evidence for the scouring effect of the binary than the core alone.

Figure 10.

Radial profiles of velocity anisotropy coefficient β for triaxial models with |$M_\bullet =10^8\, \mathrm{ M}_\odot$| and q = 1 (solid) or q = 1/9 (dashed lines). The initial model (blue) is close to isotropic, but as the binary ejects stars preferentially from radial orbits, the velocity distribution becomes more tangentially anisotropic in the central region in the course of evolution.

Open in new tabDownload slide

Finally, given that the continued shrinking of the binary depends crucially on the existence of centrophilic orbits in a triaxial potential, it is worthwhile to check whether it remains triaxial over the entire evolution. Fig. 11 plots the radial profile of the axial ratio (y/x, |$z$|/x) at the initial moment and after 1 Gyr of evolution (for a |$M_\bullet =10^8\, \mathrm{ M}_\odot$|⁠, q = 1 model, but the results are similar in other cases). It demonstrates that the system evolves towards a more axisymmetric shape (y/x ∼ 0.95 in the central parts, higher than the initial value 0.9), but remains sufficiently triaxial over its lifetime.

Figure 11.

Radial profiles of the axial ratios for triaxial models with |$M_\bullet =10^8\, M_\odot$| and q = 1 (they are almost insensitive to q). The initial model (blue line) had a constant axial ratio x|$\colon$|y|$\colon$||$z$| = 1|$\colon$| 0.9|$\colon$| 0.8, and after 1 Gyr of evolution (red line) still remains substantially triaxial despite the presence of the binary. For comparison, we show the shape of a model with a single black hole after 1 Gyr of evolution (dashed cyan line), which also did become somewhat more axisymmetric.

Open in new tabDownload slide

6 DISCUSSION AND CONCLUSIONS

In this paper, we study the long-term co-evolution of binary SMBHs and galactic nuclei, using a large suite of stellar-dynamical simulations conducted with the Monte Carlo code raga. We explore the influence of several parameters of the binary (total mass, initial eccentricity, mass ratio, inclusion or neglect of stellar captures, and GW losses) and the stellar cluster (shape and rotation). We start the simulations from equilibrium initial models for the stellar distribution, at the moment of formation of a hard binary, and conduct them for a duration of time sufficient for the binary to reach the GW-dominated regime (up to a few Gyr). Our main conclusions are the following:

• The stellar-dynamical hardening rate of the binary is almost independent of the initial eccentricity, and is nearly the same whether we allow the stars to be captured or not. Of course, if one takes the GW losses into account, a highly eccentric binary would reach coalescence sooner.

• The hardening rate depends weakly on the binary mass ratio q, but much more strongly on the geometry of the stellar cluster (spherical, axisymmetric, or triaxial). Only in the latter case the binary is able to reach the GW-dominated regime in a time much shorter than the Hubble time, regardless of its initial parameters, in agreement with Vasiliev et al. (2015) and Gualandris et al. (2017). The hardening rate is significantly lower (by a factor 3−10) than the so-called ‘full loss-cone rate’ due to the depletion of centrophilic orbits, however, their reservoir is large enough to sustain the long-term evolution of the binary. The q-dependence of the hardening rate seen in Fig. 1 arises primarily from the fact that a smaller secondary SMBH leads to a less dramatic destruction of the stellar cusp at the early stages of evolution. However, the coalescence time (equation 6) is nearly independent of q.

• The eccentricity evolution depends on its initial value and on the rotation and geometry of the stellar cluster. We confirm previous results regarding the strong tendency of eccentricity to increase or decrease in counterrotating or corotating stellar clusters, correspondingly (e.g. Sesana et al. 2011; Holley-Bockelmann & Khan 2015). However, in triaxial clusters this growth or decay saturates more quickly, because the principal supply source of stars into the loss cone comes from centrophilic orbits, which may not have net rotation. We also find that the eccentricity evolution in non-rotating clusters is rather stochastic, which has been observed in many previous studies (e.g. Wang et al. 2014; Khan et al. 2018), but on average it tends to increase if started from a sufficiently high initial value (≳0.3−0.5) or in systems with more unequal mass ratios.

• The capture rates of stars by binary SMBH are higher than in otherwise equivalent systems with a single SMBH at the early stage of its evolution, in agreement with Chen et al. (2011), Li et al. (2015), and Darbha et al. (2018). However, subsequently they drop precipitously, and in the spherical and axisymmetric cases remain well below the rates expected in single-SMBH systems, again confirming the earlier analysis of Chen et al. (2008). The decrease of capture rates parallels that of the binary hardening rates, since both are driven by the depletion of loss-cone orbits. In triaxial systems this depletion is less prominent, and the capture rates remain only slightly lower than or comparable to those in single-SMBH systems throughout the entire evolution. The capture rates do not appreciably depend on eccentricity.

• Most of the stars that interacted with the binary at later stages of evolution are ejected from the galaxy, forming a population of hypervelocity stars with an approximately exponential distribution in speed (with a characteristic scale ≳10³ km s⁻¹). In the case of eccentric binaries, the tail of this distribution extends to higher velocities, exceeding 10⁴ km s⁻¹. Their angular distribution is concentrated near the orbital plane of the binary.

• The total mass of ejected stars is comparable to the mass of the binary, with a weak dependence on the mass ratio (∝ q^1/4), in agreement with Merritt (2006) or Khan et al. (2012). At the early stage of evolution, the initially cuspy density profile is destroyed (cf. Milosavljević & Merritt 2001); however, at later stages the stars entering the loss cone arrive from larger distances, and the density profile is depleted only very slightly there. We argue that the impact of the binary on the density profile (known as core depletion or mass deficit) is very hard to assess from observations, as the difference between the (unknown) initial profile and the one modified by the binary is smaller than the uncertainties in estimating the initial profile. Most importantly, any attempt to infer the mass deficit from fitting a particular parametric model to the observed density profile cannot be accurate even to within an order of magnitude.

Our focus in this study was on the long-term evolution of stellar systems that occurs quasi-statically (the characteristic time-scales for changes in the binary orbit or the stellar distribution are much longer than the dynamical time). The secular Monte Carlo approach is well-suited for such problems, although it has a number of limitations, which we believe to be insignificant for most of our conclusions:

• We consider only stellar-dynamical processes, neglecting any gas physics, thus our analysis is applicable only to binary SMBH formed in dry mergers.

• We assume that the evolution proceeds through a sequence of nearly steady-state configurations. Thus we are not able to adequately simulate the very early stage when the binary just forms, or consider highly dynamic situations such as an infall of a globular cluster (Arca Sedda et al. 2019; Bortolas, Mapelli & Spera 2018a), a third SMBH (Iwasawa, Funato & Makino 2006; Hoffman & Loeb 2007; Ryu et al. 2018), or massive perturbers (Perets & Alexander 2008). The momentarily elevated capture rates at the early stage of binary evolution (Chen et al. 2011; Li et al. 2017) are also observed in our simulations, but cannot be quantified robustly. The duration of this early stage is short (a few Myr) compared to the overall lifetime of a binary (∼1 Gyr).

• We also assume that the binary orbit is aligned in a specific plane (x−y) and does not change its inclination throughout the simulation, which is reasonable to expect especially in non-spherical stellar potentials (Cui & Yu 2014). A detailed Fokker–Planck study of diffusive changes in the binary orbital parameters (eccentricity and inclination) has been conducted by Rasskazov & Merritt (2017), who showed that the inclination and eccentricity evolve independently, and the expected change in the inclination angle over the lifetime of the binary is moderate.

• In our approach, the centre of mass of the binary is pinned down at the origin. It was suggested by Chatterjee, Hernquist & Loeb (2003) that the Brownian motion of the binary may enhance the hardening rate. However, analytical estimates by Milosavljević & Merritt (2003) and N-body simulations by Bortolas et al. (2016) demonstrate that the stochastic wandering of the binary SMBH centre of mass has only a minor effect (⁠|${\lesssim } 10{{\ \rm per\ cent}}$|⁠) on the hardening rates. On the other hand, Holley-Bockelmann & Khan (2015) and Mirza et al. (2017) found that in stellar clusters with prograde rotation, the binary centre-of-mass settles on to a finite-size orbit around the galaxy centre; however, this global motion did not significantly affect the hardening rate, although the duration of their simulations is much shorter than considered in this paper.

• We do not consider the role of binary stars, which may influence the dynamical evolution of the single/binary SMBHs in various ways. Wang et al. (2018) explored various effects in four-body interactions between binary stars and binary SMBH, showing that these may lead to elevated rates of tidal disruptions or ejection of hypervelocity stars, possibly providing another observational signature of a binary SMBH.

• The Monte Carlo approach to two-body relaxation uses the diffusion coefficients computed in a spherical isotropic background, while the systems under study are in general non-spherical and anisotropic (especially in cases of strong rotation). Nevertheless, the effect of bulk rotation and the global shape of the potential are adequately taken into account at the level of stellar orbits (in the collisionless sense). Since the evolution of galactic nuclei in the mass regime considered in this paper (⁠|$M_\bullet \ge 10^7\, \mathrm{ M}_\odot$|⁠) is mainly caused by collisionless processes, this slight inconsistency in the treatment of collisional relaxation should not have any impact on the results.

The observational evidence of binary SMBHs is scarce, although this will hopefully change when the planned space-based GW observatory lisa (Amaro-Seoane et al. 2017) becomes operational. In this paper, we reaffirmed the expectation that the lifetime of SMBH binaries is much shorter than the Hubble time in realistic conditions, and explored their effects on the stellar distribution in the galactic nuclei (rates of stellar captures or tidal disruptions, ejection of hypervelocity stars and erosion of density cusps). Unfortunately, none of these effects constitute an observable smoking-gun evidence for the existence of a binary SMBH in any given galaxy.

ACKNOWLEDGEMENTS

This work has made use of the MIPT-60 cluster, hosted by Moscow Institute of Physics and Technology. KL thanks S.V. Ermakov for helpful discussions.

Footnotes

REFERENCES