The study of the dynamical evolution of massive black hole pairs in mergers is crucial in the context of a hierarchical galaxy formation scenario. The timescales for the formation and the coalescence of black hole binaries are still poorly constrained, resulting in large uncertainties in the expected rate of massive black hole binaries detectable in the electromagnetic and gravitational wave spectra. Here, we review the current theoretical understanding of the black hole pairing in galaxy mergers, with a particular attention to recent developments and open issues. We conclude with a review of the expected observational signatures of massive binaries and of the candidates discussed in literature to date.

1. Introduction

Understanding the formation and evolution of massive black holes (MBHs) is one of the most exciting goals of contemporary astrophysics and cosmology. It is now well established that MBHs are ubiquitous in nearby spheroids (e.g., [1]), most of them lurking in a quiescent accretion state, while during the cosmic history million-to-billion solar mass MBHs powered quasars. These objects have become in the last years central building blocks for all the proposed scenarios of galaxy formation (e.g., [2]), playing a major role in shaping galaxies through feedback processes springel [3]. Massive black hole binaries (MBHBs), formed during the galaxy merging process [4], promise to be among the most luminous gravitational wave (GW) sources for future space-borne interferometers like the proposed new gravitational wave observatory (NGO) (https://lisa-light.aei.mpg.de/bin/view/) and ongoing pulsar timing array (PTA) campaigns [5]. Theoretical modeling of MBHB dynamics is essential in addressing a number of fundamental astrophysical questions (such as the merger-quasar connection or the MBH-host relations) and in identifying putative signatures that may serve as a guidance for present and future observational campaigns.

Early stages of MBH pairing have been observed, from the initial phases of galaxy mergers, where two distinct but gravitationally bound galaxies are observable at separations of: ~100 kpc [611], down to unbound pairs of MBHs at separations of 1 kpc embedded in a single galaxy remnant [12, 13].

During this initial stage, the MBH pairing is driven by dynamical friction acting on the host galaxies. The two MBHs (here after and for the primary and secondary MBH, resp.) bind in a binary if they reach a relative separation where is the velocity dispersion of the host galaxy, is the total mass of the binary, and and are in units of 100 km s−1 and , respectively. If, in galaxy mergers, the mass of the MBHB scales with following the MBH mass versus relation (see, e.g., [14, 15] and references therein), (1) implies .

The efficiency of the process depends on the mass ratio between the merging galaxies. While equal mass (major) mergers result in a fast formation of a MBHB, in very unequal mass (minor) mergers, the satellite, tidally disrupted along the course of the encounter, leaves its naked MBH wandering in the outskirts of the most massive galaxy (see, e.g., [16]). Recent numerical studies have addressed the efficiency of the formation of a MBHB as a function of the mass ratio between the two galaxies, the gas fraction in the merging galaxies, and the redshift of the merger [17]. In particular, the presence of a significant gas component in the satellite helps the MBHB formation: during the first pericenters, the interaction between the two galaxies promotes the formation of bars, that convey a large fraction of the available gas in the center of the merging galaxies (already noticed in lower resolution merger simulations, see, e.g., [18, 19]). This new nuclear gas overdensity deepens the potential well of the secondary nucleus and prevents its tidal disruption. Callegari et al. [17] found a critical galaxy mass ratio for the formation of a MBHB of ~1/10 for gas-rich galaxies and 1/4 for gas poor galaxies. In zero-order approximation, assuming the MBH mass versus bulge mass relation [20, 21], we expect similar mass ratios (≤1) in MBHBs. However, can be increased for gas-rich galaxies. The strong gas inflows in the satellite (more perturbed) galaxy can result in a faster growth of the smaller MBH and in higher values of the expected MBH mass ratio ( [22]). (We caution that these ranges in implicitly assume similar morphology in the two merging galaxies, although minor mergers can involve very different morphological types. A study of MBH pairing and the formation of MBHBs in mixed mergers (e.g., disks merging with ellipticals) is not available to date.)

As a consequence, efficient dynamical friction promotes the formation of MBH binaries with similar mass ratios. The expected number of observable binaries and the rate of MBH coalescences, however, depend on the dynamical evolution after the binary formation. In order to coalesce through GW emission in less than an Hubble time, the two MBHs have to reach a separation where is a function of the binary eccentricity [23]. Circular equal mass binaries can coalesce shrinking by a factor while this factor decreases to 1 for binary with . In order to understand the final fate of a MBHB and to constrain theoretically its observability, it is fundamental to study at the same time the orbital decay and the eccentricity evolution of the binary. The final fate of the MBHs, that is, if they will coalesce in a single object or not, strongly depends on the amount of matter (stars and gas) they can interact with after the binary formation. A definite answer is not present to date. In this paper, different scenarios will be discussed, depending on the nature of the galaxy mergers and on the properties of the galaxy nuclei.

This paper is organized as follow. In Section 2, we review the dynamical evolution of MBHB in gas poor environments, while the effect of gas is discussed in Section 3. In Section 4 we describe the MBHB candidates observed to date. Finally, our conclusions are drawn in Section 5.

2. Dynamical Evolution in Gas-Poor Environment

In systems where dynamical friction is efficient in dragging the two MBHs to the center of the merger remnant, the now bound MBHB is inevitably embedded in a gas- and star-rich environment. Such rich ambient provides a variety of physical mechanisms to efficiently extract the energy and angular momentum of the MBHB, promoting its final coalescence. In this section, we focus on dynamical processes involving interactions with stars. MBHBs in pure stellar environments were the first to be examined (the basics going back to [4]), for the obvious reason that stars can be considered as point particles, affected by gravitational forces only. The MBHB-star interactions are therefore adequately described by Newton’s laws only, without all the complications involved in gas dynamics. Nonetheless, a single star-binary interaction is, by definition, a three-body problem, and the dynamics of the system is inevitably chaotic. Therefore, no simple analytical solutions are viable, and numerical studies (both involving three body scatterings and full N-body simulations) have been massively exploited to tackle the problem. The fate of the MBHB is determined by its semimajor axis and eccentricity evolution (see Section 1); in the following, we discuss them separately.

2.1. Shrinking of the Binary Semimajor Axis

The basic physical process driving the MBHB evolution in presence of stars is the slingshot mechanism. A star intersecting the MBHB orbit undergoes a complex three-body interaction being eventually ejected at infinity, carrying away energy and angular momentum from the binary. Extensive three-body scattering experiments [2426] have shown that ejected stars carry away an energy per unit mass of the order of where is the reduced mass of the binary. Assuming a classical interaction rate given by , where is the binary cross section, is the number density of the ambient stars, and is their typical velocity “at infinity” with respect to the binary (i.e., the velocity dispersion of the stellar system). Quinlan [25] showed that the evolution of the binary semimajor axis is simply given as where is a dimensionless hardening rate. If , independently on the mass, mass ratio and eccentricity of the system. In principle, given a stellar system with density and velocity dispersion , (5) predicts efficient coalescence of the MBHB.

However, the above simple treatment ignores the concept of loss cone depletion. In an extended stellar system, only a tiny fraction of the stellar phase space allows orbits intersecting the MBHB, commonly referred as “binary loss cone.” As stars are ejected, the loss cone is depleted, and the binary evolution is governed by the rate at which new stars are fed into the loss cone [27]. In typical stellar systems, the mass in stars in the binary loss cone is of the order of few times [28], insufficient to reach in most of the cases [29]. This is the origin of the “last parsec problem” [30]. In a spherical stellar system, the loss cone refilling proceeds on a two-body relaxation timescale [31], which is usually much longer than the Hubble time. In the last decade, this fact has been confirmed in N-body simulations [27, 30, 32]. In such simulations, after loss cone depletion, further hardening was provided by two-body relaxation. This is a process that depends on the “granularity” of the systems, and the result is an N-dependent hardening rate, with the binary evolution slowing down as the number of particle in the simulation increases. Extrapolating these results to a realistic N representative of a galactic bulge, the binary evolution would have stalled.

In recent years, evidence has emerged that the “last parsec problem” might be an artificial product of the “spherical cow” approximation which is often exploited in astronomy. Basically, the spherical systems studied in the simulations represent a worse case (and unrealistic) scenario. MBHBs are infact produced in galaxy mergers, in which the resulting stellar bulge is rotating, triaxial, and likely to undergo bar-like instabilities. In a triaxial potential, an orbiting star does not conserve any of its angular momentum components [31]. As a result, there is a vast family of orbits (called centrophilic) that are allowed to get arbitrarily close to the binary [3335], keeping the loss cone full during the MBHB hardening process. Recent N-body simulations have confirmed this scenario. Berczik et al. [36] studied the evolution of a MBHB in a rotating bulge. In this case, the stellar system experiences a bar instability resulting in a triaxial potential. The binary hardening rate was found to be N-independent; a proof that the hardening was not proceeding because of spurious two-body relaxation. More recently, the advent of GPU computing made possible to simulate “ab initio” the evolution of two interacting stellar bulges hosting MBHs; a first step toward a realistic galactic merger scenario. Several simulations were performed by Khan et al. [37], Preto et al. [38], and Gualandris and Merritt [39]. In all cases, the stellar remnant was triaxial and rotating, and the hardening rate, given by triaxiality driven loss cone replenishment, was found to be independent on N, implying coalescence timescales of 108 yr. Remarkably, when normalized to the merging galaxy properties, the binary hardening rates found in these simulations follow (5) where [39]. This is a consequence of the fact that, whatever is the geometry of the system, the average “quantum” of energy taken away from an interacting star is always the same and the evolution of the system is determined by the star-binary interaction rate only.

2.2. Eccentricity Evolution

As pointed out in the introduction, eccentricity plays an important role in driving the binary coalescence. However, addressing the eccentricity evolution of the system is more complicated because the caused by each individual interaction depends on a combination of energy and angular momentum exchanges. The angular momentum distribution of the interacting stars is therefore crucial. The eccentricity evolution can be described as Here is a dimensionless parameter that, differently than , depends on the binary mass ratio and eccentricity itself [25, 26]. In general is a positive number in the range 0–0.3 (the peak value occurs at ), meaning that the binary eccentricity grows during the shrinking process. Sesana [40] constructed complete binary evolutionary tracks by coupling three-body scattering experiments of bound and unbound stars to an analytical description of the stellar distribution and on the loss cone refilling. As a general trend, quasi-circular-equal mass MBHBs experience just a mild eccentricity growth, while systems which are already eccentric at the moment of pairing, or with significantly lower than 1, can evolve up to .

The eccentricity evolution in stellar environments has been tackled by several authors by means of full N-body simulations. However, the limited number of particles () in such simulations results in very noisy behavior for the binary eccentricity, and it is difficult to draw conclusions about the general trends behind the numerical noise. Milosavljevic and Merritt [30] carried out numerical integration of equal MBHBs embedded in two merging isothermal cusps ( with ). Starting with circular orbits, they find a mild eccentricity increase to a value of 0.2 during the stellar-driven hardening phase. Merritt et al. [32] considered equal MBHBs embedded in Dehnen density profiles [41] with with different initial eccentricities. Again, they find that circular binaries tend to stay circular, while eccentric binaries tend to increase their eccentricities in reasonable agreement with the prediction of scattering experiments. Simulations starting before the formation of MBHBs carried by Hemsendorf et al. [42] and Aarseth [43] produce binaries with at the moment of pairing, with subsequently increasing up to 0.95. Amaro-Seoane et al. [44] focused on binaries of intermediate MBHs () in massive star clusters. Coupling full N-body simulations to three-body scattering experiments, they find binaries with significant eccentricity (~0.5-0.6) at the moment of pairing, growing up to >0.9 during the hardening phase; similar conclusions are reported by Amaro-Seoane et al. [45]. On the small side, simulations were performed by Baumgardt et al. [46] and Matsubayashi et al. [47], assuming a stellar density profile , motivated by the analytical equilibrium solution for a dense relaxed stellar cusp around a massive object [48]. When properly rescaled, the eccentricity increase found in both papers agrees remarkably well with predictions based on the hybrid model by Sesana et al. [49]. Iwasawa et al. [50] investigated in detail the angular momentum exchanges between the binary and the stars responsible for the eccentricity growth, in bound stellar cusps. In particular, they showed that stars counterrotating with the binary tend to extract a lot of angular momentum from the MBHB, causing the eccentricity growth, whereas corotating stars do not. This is a simple consequence of angular momentum conservation during the ejection process, as shown by Sesana et al. [51].

The evolution of the binary eccentricity can be extremely different for nonisotropic systems. For example, Dotti et al. [52] showed that at large scales before the formation of a binary dynamical friction exerted by rotationally supported stellar disks tend to circularize the orbit of a MBH pair. At smaller separations, Sesana et al. [51] demonstrated, that in a rotating stellar system, the eccentricity evolution of unequal MBHB is dramatically affected by the level of co/counterrotation of the stellar distribution with respect to the binary, with corotating distributions promoting circularization rather than eccentricity growth. Nonetheless, most of the simulations involving rotating bulges [36, 53] or merging systems [3739] find quite eccentric binaries at the moment of pairing (ranging from 0.4 to 0.8), and the subsequent evolution leads to a general eccentricity growth, in good agreement of what predicted for an isotropic stellar distribution. This may be because the binary evolution is mostly driven by loss cone refilling of unbound stars on almost radial orbits, with negligible initial angular momentum.

Overall, the emerging general picture favors efficient coalescence of MBHBs in dense stellar merger remnants. The triaxial and rotating nature of the stellar distribution promotes efficient loss cone refilling, while large eccentricities (especially in unequal mass systems) shorten the gap between the binary pairing and the efficient GW emission stage. In the near future, massive N-body simulations with several million particles will offer a unique opportunity to confirm this scenario.

3. Dynamical Evolution in Gas-Rich Environment

3.1. Formation of a MBH Binary in a Circumnuclear Gas Disks

As discussed in the Introduction, in comparable mass, gas-rich galaxy mergers, the gravitational interaction drives strong inflows of gas toward the galactic centers. The numerical multiscale investigation of an equal mass galaxy merger discussed in Mayer et al. [54] revealed that, in advanced stages of the galaxy merger, the two MBHs, orbiting in the central 100 pc of the merger remnant, are embedded in a dense, rotationally supported, gas disk (see Figure 1). This circumnuclear disk is self-gravitating and can be up to ~500 times more massive than the MBH pair [54]. The dynamical evolution of the two MBHs is driven by dynamical friction, and, since the circumnuclear disk is the densest structure in the remnant nucleus, it is the main cause of their orbital decay. Mayer et al. [54] followed the evolution of the two MBHs from the initial stages of the galaxy merger down to 5 pc where they form a binary, as the mass of gas enclosed within their separation is less than the mass of the binary.

The evolution of MBHs in circumnuclear disks has been studied in details in dedicated simulations, in which the former evolution of the MBHs at distances 100 pc is not explored. This allows to achieve a better resolution in the central region of the remnant and to study the latter MBH pairing. Similar independent investigations discussed in Escala et al. [55] and Dotti et al. [52, 56] agree with Mayer et al. [54] on a rapid (on a timescale of 106 yr) formation of a MBHB at parsec separations. The higher resolution in these studies allow for a further decay of the binary, that reaches sub-pc separations comparable to the spatial resolution in these numerical studies.

During this intermediate stage (100 pc 0.1 pc), the interaction with the circumnuclear disk strongly affects the eccentricity of the MBH pair. Dotti et al. [52] first noted that, for MBHs with an initial significant eccentricity, the decay phase is preceded by a circularization of the orbit. This is due to the dynamical friction exerted by a rotating background onto the BHs at their apocenters. Since the circumnuclear disk is rotationally supported, the BHs at their apocenter move more slowly than the local gas and are dragged in the direction of their motion. This positive torque is exerted only in proximity of the apocenter, where the angular momentum of the MBHs can be maximally increased. This effect is a general feature of dynamical friction exerted by a rotating background. Dotti et al. [52] demonstrated that the same effect is present for MBHBs orbiting in a stellar circumnuclear disk (as discussed in Section 2.2), and Callegari et al. [17] found the same effect acting in unequal mass galaxy mergers at larger scales, where the disk of the primary galaxy circularizes the orbit of the satellite.

3.2. Evolution of Close MBH Binaries in Circumbinary Gas Disks

After the formation of a close binary, the MBHs, acting as a source of angular momentum, exerts a tidal torque that inhibits the gas from drifting inside its orbit. This creates a hollow density region, called gap that surrounds the binary (e.g., [5761]).

As a consequence of disk clearance, corotation and inner Lindblad resonances are reduced in power, drastically changing the dynamical evolution of the binary. The same transfer of angular momentum that keeps the disk from infalling onto the MBHs drives the shrinking of the binary. Analytical and numerical studies agree in predicting that the interaction between the binary and the surrounding material reduces the semimajor axis of the binary and increases the eccentricity of an initially low eccentricity binary (e.g., [6170]).

The MBHs/circumbinary disk interaction can be studied in two-limit cases: (i) assuming that the MBHB is embedded in a virtually infinite disk, as in the case in which it is continuously refilled by a long-lived larger-scale structure, or (ii) assuming that the disk has a finite mass (and, as a consequence, a finite reservoir of energy and angular momentum).

In the first limit of a constantly fuelled disk, analytical models of the evolution of the binary are available. The orbital decay timescale is where is an estimate of the disk mass inside the orbit of the secondary MBH. As long as is greater than , the MBH behaves as a fluid element, decaying on the viscous timescale onto the primary. As the binary orbit shrinks, decreases. When the enclosed mass becomes comparable with the secondary its decay timescale increases to resulting in very large migration timescales at small separations. However, as first noticed in Ivanov et al. [65] in a continuously refilled disk, the stalling of the binary would cause a steady increase of the density of the inner edge of the circumbinary disk, until the mass close to the binary becomes again comparable to the secondary. At this stage, fast migration starts again. Following Ivanov et al. [65] and Cuadra et al. [70] estimated that, for a disk on the verge to undergo fragmentation (i.e., as dense as possible), the coalescence timescale of a binary with , starting from an initial separation of <0.05 pc, is <109 yr. This timescale decreases with the decrease of the binary mass. Promisingly, the initial separation assumed here is close to the limit achieved in the larger scale simulations discussed in Section 3.1. Note, however, that such a timescale is comparable with the age of the Universe at . If the migration in a dense circumbinary disk is the fastest process driving the MBHs coalescence, no coalescences of binary with are expected at .

In a similar way, the study of the evolution of is possible. For simple alpha disks [71], Artymowicz and Lubow [59] found that the outer edge of the gap (i.e., the inner edge of the circumbinary disk) depends on the binary eccentricity, with more eccentric binaries opening larger gaps. However, the eccentricity itself increases as a consequence of the interaction with the expanding disk [62, 64, 67, 68, 70] resulting in a steady expansion of the gap. These two effect together would result in a limiting eccentricity of ~1, leading to a fast coalescence due to efficient gravitational wave emissions (see (2)). A smaller limiting eccentricity may result from a less efficient coupling between the gas and the binary, as expected if the disk moves farther out. Since the evolution of in an expanding gap has been performed with finite mass disks, we postpone its discussion to the following.

If the circumbinary disk is limited in mass, the evolution of the orbital separation and eccentricity can be quite different. Since, in this scenario, the disk is not continuously refilled from the outside, the gas mass initially available in the disk is the key parameter.(i)If the disk is ≫, the evolution is similar to the infinite-mass disk, with the binary coalescing on a short time-scale.(ii)If the mass in the circumbinary disk is less or of the order of , the interaction with the binary forces the whole circumbinary disk to move outward in few orbital periods. This expansion of the circumnuclear disk has been observed in simulations in which the components of the binary have similar masses (e.g., [62, 70]). In this case, the gas reaches distances 4a, at which the interaction with the binary is not efficient anymore. A small amount of gas can then fall again closer to the binary because of orbital angular momentum exchange with the bulk of the disk, but most of the gas would never get close enough to significantly alter the evolution of the binary. Note that the expansion of the circumbinary disk is most effective for eccentric binaries. Binaries do not coalesce because of the interaction with too small (), nonrefuelled disks.

In this simple description, the mass in gas within the disk is either accreted onto the MBHs or conserved. However, the gas in such a dense environment could be consumpted by star formation, decreasing the effective mass of the circumbinary disk. As a consequence, the disk could be initially ≫, but decreasing in mass with time, and could possibly fail in bringing the binary to the final coalescence. This scenario has been recently discussed in Lodato et al. [72]. In this investigation, stars are allowed to form in the disk whenever it becomes gravitationally unstable. The rate of new star formation is obtained requiring that they would inject enough energy to keep the disk on the verge of fragmentation (i.e., providing an heating term exactly equal to the cooling losses in the disk). Even considering such a simple “thermal” feedback from the newly formed stars, the disk loses so much mass that the binary cannot reach the final coalescence, unless its initial separation is 0.01 pc (for the MBH masses considered in the paper, and ).

Other possible feedback terms, not included in this model, that can help in preventing such a strong gas consumption have been suggested by Lodato and collaborators, such as momentum feedback from stellar winds and supernovae explosions. Furthermore, the interaction between the MBHs and the forming gas clumps and stars, not considered in the investigation, could help the binary decay. The consumption of gas may not be a problem if large inflows of gas are present, as in the case of a continuously refuelled disk discussed above.

In the finite mass disk scenario, the existence of a limiting eccentricity has been studied in Roedig et al. [62] through a suite of high-resolution SPH simulations. They find a critical value . In these simulations, the initial ratio between the gap size and the semimajor axis of the binary is 2 and can increase during the runs up to more than 4, when the interaction efficiency drops. The analytical model presented in the paper agrees with the simulations, predicting the limiting eccentricity to be

The initial choice of is somewhat arbitrary. In reality, the feeding of a MBH binary forming in a gas-rich galaxy merger can be a very dynamic process, and the interaction with a single circumbinary disk could be too idealized a picture. Larger scale simulations show episodic gas inflows due to the dynamical evolution of the nucleus of the remnant (see, e.g., [73, 74]). In this scenario, the binary can still interact with a disk and excavate a gap, but the size of it would be time dependent (as in the simulations presented here) and would depend on the angular momentum distribution of the inflowing streams, resulting in a range of .

Note that the discussion above implicitly assumes that the MBHB and the circumnuclear disk corotate with each other. This is the natural outcome of a evolutionary sequence in a gas-rich galaxy mergers, in which the two MBHs orbit in a large scale circumnuclear disk, are forced to corotate with it [56] and open a gap in the very central region of the gas distribution. This picture, however, could not apply in gas-poorer mergers, or even in a gas-rich scenario, if the circumnuclear disk formed during the merger fails in bringing the two MBHs to the final coalescence before it is consumed by star formation and/or MBH feedback. In one of these cases, an occasional small inflow of gas could happen with a random angular momentum and could form a retrograde circumbinary disk, counterrotating with respect to the binary.

The evolution of a MBHB in a retrograde disk has been discussed in Nixon et al. [75]. In this case, the gravitational interaction between the binary and the gas brakes both the components, so that, unlike in the prograde scenario, here the torques responsible for the binary shrinking and its eccentricity evolution cause the edge of the disk to move inwards. The binary and the secondary MBH in particular, experience the presence of a closer distribution of gas, which would imply a faster evolution, moving at a higher relative velocity, that results in a less effective interaction. Nixon and collaborators show that (i) if the binary is initially not exactly circular (, where is the aspect ratio of the disk) and (ii) if it interacts with the disk mainly at the apocenter, then the secondary evolves onto an almost radial orbit after interacting with a gas mass comparable to its own. (Note that assuming the interaction to take place mainly at the apocenter is in agreement with the circumbinary disk to form after the MBHs bind in a close binary, as a consequence of a randomly oriented accretion event.)

The increase of the eccentricity to in the retrograde case is due to direct accretion of linear momentum from counterrotating material. Since the secondary has null radial velocity at the apocenter, before and after the interaction, and far from the apocenter, the secondary is assumed to move on an unperturbed Keplerian orbit, the apocenter is constant. Interacting with counterrotating gas, the secondary decreases its angular momentum, reducing its semimajor axis (of up to a factor of 2) and, most importantly, increasing its eccentricity (up to 1). At very high eccentricities the emission of gravitational waves can bring the binary to the final coalescence in less than an Hubble time (see (2)).

Note that this scenario suffers of the same disk consumption problem as the prograde one. If the disk is consumed by star formation or evacuated by MBH or supernovae feedbacks, the process stops. This makes this process particularly interesting for very unequal mass binaries (, less likely to form from galaxy mergers, as discussed in the Introduction). Fast inflows of gas, on timescales shorter than the consumption time, would help the coalescence in both the prograde and the retrograde scenario.

4. Binary Candidates in the Realm of Observations

Despite being a natural outcome of galaxy mergers, MBH pairs are still elusive. Less than 20 systems with separations of ~10 pc to ~10 kpc pairs are of this kind are known to date. MBHs orbit in the common postmerger stellar environment, in-spiralling because of dynamical friction. They appear as a single galaxy (eventually, with disturbed morphology) with two active nuclei. Examples are the prototypical case of NGC 6240 [12]; the radio galaxy 3C75 [76]; the spiral galaxy NGC 3341 [77]; the ULIRG Mrk463 [78]; the interacting galaxy COSMOS J100043.15+020637.2 [79]; the quasar pair J1254+0846 [80]. All these objects have been discovered because of the presence of two resolved X-ray sources wandering in the merged galaxy. In order to look for these systems, an alternative approach is to search for objects with two systems of narrow lines at slightly different redshifts [81, 82]. Large spectroscopic surveys, like the Sloan Digital Sky Survey (SDSS [83]), have been used to search for these systems. Follow-up observations were then used to discriminate between dual AGN and single AGN with complex gas dynamics in the narrow line region [81, 84].

At separations of 10 pc, the two MBHs start experiencing their own gravitational interaction, binding in a binary. These systems cannot be spatially resolved in optical and X-ray observations, and radio interferometry is required. This has been successfully done only in the case of 0402+379 [85, 86]. The two flat-spectrum radio sources, corresponding to the two components of the candidate MBHB, have a projected separation of ~7 pc (few milliarcsec at ). We note, however, that radio interferometry at very high spatial resolution is not an efficient technique to search for rare objects as MBHBs, because of the limited field of views and the requirement that both the MBHs are radio-luminous (see, f.i., [87]).

Another approach to look for MBHBs is to study periodic variations in the luminosity of some AGN. The only MBHB candidate selected on these bases up to now is the BL Lac object OJ287 (see [88], and references therein). It shows a 11 yr periodic flaring. In the MBHB scenario, the secondary MBH orbits on a tilted plane with respect to the accretion disk of the primary. The flares are associated to the passages of the secondary MBH through the nodes. However, this is not the only available explanation of the peculiar behavior of this source [89].

More promisingly, signatures of MBHBs have been searched for in optical and NIR spectroscopic databases. According to the MBHB hypothesis, if at least one of the MBHs is active, the broad lines (BLs) emitted by gas bound to it may be red- or blue-shifted with respect to their host galaxy systemic recessional velocity, as a consequence of the Keplerian motion of the binary [4]. Therefore, looking for quasars with significant velocity shifts (>few hundreds km s−1) can be a valid approach to systematically search for MBHB candidates over large fields. This technique does not suffer of any angular resolution limitations. Actually, the closer (and more massive) the binary is, the more shifted/deformed the BLs are. Five MBHB candidates have been individually found in this way: SDSS J0927+2943 [9092], J1536+0441 [93], J1050+3456 [94], J1000+2233 [95], and J0932+0318 [96]. Tsalmantza et al. [97] recently applied this technique in a systematic way over the whole spectroscopic sample of SDSS, using a method developed for searching unresolved gravitational lenses [98]. This analysis resulted in 4 new MBHB candidates (J1012+2613, J1154+0134, J1539+3333, and J1714+3327) and the confirmation of all the previously known objects (see Figure 3). In a similar study, Eracleous et al. [99] searched for objects with anomalous H profiles in the SDSS quasar catalogue [100]. Among them, they identified 88 sources showing velocity shifts between the broad H line peak and the rest frame of the narrow emission lines.

It should be noted that alternative interpretations for the spectral properties of the known candidates are available.(i)Modest line shifts (500 km s−1) are often observed in “normal’’ AGN [101].(ii)Similarly, small velocity shifts (<4000 km s−1) can be associated to the remnant of a binary coalescence, recoiling because of anisotropic gravitational wave [90]. (However, if the galaxy merger is gas rich, the maximum recoil velocity is expected to be <100 km s−1  [102105].) (iii)An unobscured MBHB with both MBHs active could resemble the spectrum of a double peaked emitter (see, e.g., [106]), where broad double-peaked lines are emitted because of the almost edge-on, disk-like structure of the broad line region of a single MBH.(iv)Chance superposition of two AGNs (or an AGN-galaxy superposition) within the angular resolution of the used spectrograph can also mimic velocity shifts of different line systems [107].

The simplest way to discriminate between these scenarios and the MBHB hypothesis would be to look for a periodic oscillation of the broad line shifts around the host galaxy redshift. However, the orbital period of the binary could be too long to be easily observed [4]. Noticeably, Eracleous et al. [99] observed a variation in the shifts at two different epochs in 14 out of 88 candidates with resulting accelerations between −120 and +120 km/s/yr. Longer temporal baselines are needed to prove the MBH binary interpretation for these objects.

In order to increase the number of known MBHB candidates, and to confirm their interpretation, it is therefore of fundamental importance to identify new signatures of MBHBs. The simultaneous observation of various MBHB signatures could represent the only way to firmly validate the MBHB scenario in the known candidates. A possibility is to look for peculiar flux ratios between high- and low-ionization broad lines. The broad line region of each MBH in a binary can be perturbed and disrupted by the gravitational potential of the companion. External shells of the broad line region (where most of the low-ionization line flux is emitted) are affected first, resulting in peculiar flux ratios [108]. This criterion is particularly interesting for quasars at high redshift (), where high- and low-ionization lines are observable in large surveys as the SDSS, while the most prominent narrow lines, needed to measure a shift related to the orbital motions, are not present in the SDSS spectra anymore.

At even closer separations between the two MBHs, when the size of the BLR is significantly larger than the semimajor axis of the binary ( pc), the optical and UV spectral features discussed above become more complex and not directly related with the period of the binary [109]. However, in this case, typical MBHB periods are 10 yr, opening the interesting opportunity of directly detect periodic variability of the system, related to the periodicity of the accretion flows [110112]. Such close separations are particularly interesting, since they will be proven by the ongoing and future pulsar timing arrays (PTAs,[5]). In this context, Sesana et al. [112] estimate that up to few hundred MBHBs contributing to the GW signal in the PTA band may be identified through their periodicity in future X-ray all-sky surveys. Among those, few exceptionally bright sources may be resolved both in the GW and in the electromagnetic window through the detection of peculiar double K- fluorescence lines, offering a unique multimessenger astronomy opportunity. An alternative possibility, suggested by Tanaka et al. [113], is that the presence of a circumbinary cavity results in a suppression of the UV soft-X emission. MBHBs close to coalescence may therefore be identified as AGNs with exceptionally faint UV X-ray continuum.

5. Conclusions

We reviewed the most recent findings about the dynamical evolution of MBHBs and their detectability. Regarding the binary dynamics, in the last few years we recognized the importance of the medium/large scale galactic structures (100 pc) in the dynamical evolution of the binaries. In gas free environments, the shape of the bulge potential is directly related to the possibility of a MBHB to reach the final coalescence. In a spherical system, the stars that can interact with the binary (i.e., that orbit within its loss-cone) are evacuated from the center before the binary can coalesce (e.g., [30]). A fast refilling of the loss cone, that can result in the merger of the two MBHs, is possible in triaxial systems, in which the angular momentum components of the orbiting stars are not conserved (e.g., [33]). Recently, thanks to the advent of GPU computing, large-scale galaxy mergers proved the occurrence of such a replenishment in more realistic scenarios [3739]. The evolution of the binary eccentricity depends on the dynamical properties of the core as well. If the MBHs are embedded in a nonrotating stellar system, the general trend is towards an increase of with time (e.g., [25]). This promotes the coalescence of the MBHs, since gravitational wave emission is more efficient for eccentric binaries. In rotating systems, on the other hand, the evolution of depends on the orientation of the binary: a binary corotating with the stellar cusp tends to decrease its eccentricity, while, in the counterrotating case, grows up to 0.95 [51]. This clear cut scenario could be modified by the interaction with stars on quasiradial, centrophilic orbits. A study of the orbital properties of these stars has not been presented to date.

Galaxy mergers can easily promote strong inflows of gas towards the center of the galaxy remnant (e.g., [54]). Hence, it is fundamental to understand how the presence of massive gas structures in the cores modify the dynamical evolution of the forming MBHBs. It has been proven that the interaction with circumnuclear (~100 pc) disks can result in a fast (107 yr) formation of a MBHB (e.g., [55]). After this fast transient, the binary is thought to open a gap in the central gas distribution (e.g., [57]), and any further evolution is mediated by the interaction between the MBHs and the inner edge of the circumbinary disk.

Simulations and analytical studies about the interaction of MBHBs and circumbinary disks have improved our knowledge of the physical processes in play and their effect onto the binary evolution. However, a complete understanding of the binary/disk interaction is still to come. For example, in many investigations, the circumbinary disk is assumed to be cylindrically symmetric, that is, the study is reduced to an effective one-dimensional problem. The assumption of cylindrical symmetry removes any possibility of studying the effects of structures in the disk and, most importantly, of gas streams periodically inflowing from the disk onto the two MBHs, routinely observed in simulations (e.g., [62, 70, 110112]). The torques exerted by these inflowing streams have not been studied in detail yet (with the notable exception of [114]) and could provide additional help (or resistance) in bringing the binary to the final coalescence.

The final fate of a binary embedded in such a circumbinary disk is still debated. If the disk is continuously refueled from any larger-scale gas distribution, a fast coalescence can easily be achieved. However, if the binary cannot interact with enough gas (e.g., because it turns in stars), the circumbinary disk gets evacuated and fails in bringing the MBHs to coalescence (e.g., [72]). As a consequence, as in the stellar scenario, the final fate of the binary depends on the properties of larger scale structures and its ability to efficiently refuel the proximity of the MBHB with fresh gas. In principle, in presence of an intense inflow toward the center, the binary could fail in opening a gap at all, and the interaction between MBHs and a closer/denser gas structure could result in a faster coalescence of the MBHB [55]. Furthermore, if the angular momentum of gas can be efficiently reshuffled, inflowing streams could form counterrotating circumbinary disks, that can also promote the binary coalescence [75]. Only recently, the formation of a gap has been observed in large scale simulations [55, 115], in which the evolution of an extended (~100 pc) massive disk is followed (as massive as the binary in Escala et al., up to 10 in Dotti et al., see Figure 2). The spatial and mass resolution of these simulations do not allow yet a detailed study of the sub-pc evolution of the binary, down to a possible coalescence. Simulating the evolution of a sub-pc binary starting from large scale initial conditions, that can constrain the properties of the nuclear inflows together with the evolution of the binary, is the fundamental improvement needed to build a coherent picture of MBHBs in gas-rich environment.

To summarize the recent findings present in literature, we can draw a comparison between orbital decay timescales obtained considering different scenarios, for equal mass binaries.(i)In Dense Stellar Environments: if the loss cone of the MBHB is constantly refilled (see Section 2.1), a binary of with an initial separation of  pc will coalesce in years, while a binary 100 times more massive will inspire for about years before reaching the final coalescence [38, 40].(ii)In Gas Rich Environments: if the interaction with a steady, long-lived, corotating, maximally massive circumbinary disk is responsible for the MBHB orbital decay (see Section 3.2), the orbital decay timescale for a equal mass binary is less than the age of the Universe if the binary starts at an initial separation  pc (~10 times smaller than the separation at which the MBHs bind in a binary ) and is 2 orders of magnitude shorter for  pc [70]. For a equal mass binary, this timescale is less than the age of the Universe if  pc [70], ~100 times less than .(iii)In Gas Rich Environments: if the MBHB interacts with a continuous sequence of counterrotating accretion disks with an accretion rate corresponding to the MBH Eddington limit, the orbital decay timescale is of the order of 108 yr, regardless the binary mass, for . In this case, the coalescence timescale increases linearly with the inverse of the accretion rate.

Note that the estimates in gas-rich environment should be considered as lower limits, since they assume continuous accretion and, as stressed above, a continuous refuelling of the disks from larger scales. A single, not refuelled disk with enough mass and angular momentum to bring the MBHs to the final coalescence form  pc would undergo fragmentation and star formation, as discussed in [72].

The presence of gas close to the binary is necessary to its detection. If at least one of the two MBHs is active, the orbital velocity of the binary can result in a frequency shift between the broad emission lines and the narrow emission lines (e.g., [4]). This shift is expected to change periodically, on the orbital period. Using this technique, few tens of MBHB candidates have been selected (e.g., [97, 99]). For all the candidates discussed to date, possible explanations other than a MBHB have been proposed (see Section 4). Moreover, the orbital period expected for such binaries is often too long to be observed, thus the periodic variation of the velocity shift cannot be used (yet) to determine the real nature of the candidates. It is therefore fundamental to couple this with other (independent) signatures of MBHBs, to confirm their nature.

In the near future, space-based interferometers like NGO will detect the GWs emitted in the late inspiral and final coalescence of MBHBs, opening the fascinating prospective of multimessenger astronomy. The identification of an electromagnetic counterpart to a gravitational wave detection of merging MBHs can teach us about relativity, accretion physics, galaxy formation and evolution, and, more in general, cosmology. A large variety of potential EM signatures have been proposed (see, e.g., [116] and references therein). In particular, thanks to the recent quick progresses in numerical relativity, simulations of the coalescence process in presence of matter is now becoming possible [117120] and in the coming years may provide valuable insights about the nature of the associated electromagnetic signals. However, most of the counterparts discussed to date depend on some assumptions on the distribution of gas and/or stars in the immediate vicinities of the MBHB. A theoretical comprehension of the dynamics of matter close to the binary, necessary to understand the fate of a MBHB, is also fundamental to constrain the observability of an electromagnetic counterpart.

Tracing the formation, evolution, and final fate of MBHBs is certainly one of the open challenges of contemporary astrophysics. A better understanding of the interaction between MBHs and gas and the prediction of new peculiar observational features of MBHBs are needed to unambiguously constrain their properties and demography, adding an important missing piece to the galaxy evolution puzzle.


The authors are grateful to Jorge Cuadra, Monica Colpi, and Paraskevi Tsalmantza for useful discussions and comments.