∎

11institutetext: T. Bogdanović 22institutetext: School of Physics and
Center for Relativistic Astrophysics
Georgia Institute of Technology
837 State St. NW, Atlanta GA 30332, USA
22email<EMAIL_ADDRESS>33institutetext: M.C. Miller 44institutetext:
Department of Astronomy
University of Maryland
4296 Stadium Dr., College Park, MD 20742, USA
44email<EMAIL_ADDRESS>55institutetext: L. Blecha 66institutetext:
Department of Physics
University of Florida
2001 Museum Rd., Gainesville, FL 32611, USA
66email<EMAIL_ADDRESS>
# Electromagnetic counterparts to massive black hole mergers

Tamara Bogdanović  M. Coleman Miller  Laura Blecha

(Received: date / Accepted: date)

###### Abstract

The next two decades are expected to open the door to the first coincident
detections of electromagnetic (EM) and gravitational wave (GW) signatures
associated with massive black hole (MBH) binaries heading for coalescence.
These detections will launch a new era of multimessenger astrophysics by
expanding this growing field to the low-frequency GW regime and will provide
an unprecedented understanding of the evolution of MBHs and galaxies. They
will also constitute fundamentally new probes of cosmology and would enable
unique tests of gravity. The aim of this Living Review is to provide an
introduction to this research topic by presenting a summary of key findings,
physical processes and ideas pertaining to EM counterparts to MBH mergers as
they are known at the time of this writing. We review current observational
evidence for close MBH binaries, discuss relevant physical processes and
timescales, and summarize the possible EM counterparts to GWs in the
precursor, coalescence, and afterglow stages of a MBH merger. We also describe
open questions and discuss future prospects in this dynamic and quick-paced
research area.

###### Keywords:

accretion, accretion disks black hole physics gravitational waves galaxies:
nuclei radiation mechanisms: general quasars: supermassive black holes

††journal: Living Reviews in Relativity

###### Contents

  1. 1 Introduction
  2. 2 Observational evidence for close massive black hole binaries
    1. 2.1 Emission of gravitational waves
    2. 2.2 Direct imaging of double nuclei
    3. 2.3 Photometric measurements of quasi-periodic variability
    4. 2.4 Spectroscopic measurements of offset broad emission lines
  3. 3 Physical processes and relevant timescales
    1. 3.1 Rates and typical redshifts
    2. 3.2 Effects before MBHs are bound to each other
      1. 3.2.1 Stellar and gas dynamical friction
      2. 3.2.2 Implications for formation and properties of gravitationally bound MBHs
    3. 3.3 Effects after MBH are gravitationally bound to each other
      1. 3.3.1 Three-body scattering of stars and interactions with additional MBHs
      2. 3.3.2 Interactions with the circumbinary disk
    4. 3.4 The gravitational-wave phase
  4. 4 Emission signatures before, during, and after coalescence
    1. 4.1 Electromagnetic radiation before the merger
      1. 4.1.1 Emission from the circumbinary disk
      2. 4.1.2 Emission from the mini-disks and accretion streams
      3. 4.1.3 Other modulation possibilities before the merger
    2. 4.2 Signatures during or immediately after the merger
      1. 4.2.1 Reaction of matter to dynamic gravity
      2. 4.2.2 Effect of dynamical spacetime on magnetic fields and formation of jets
      3. 4.2.3 Retreat of the disk due to mass-energy loss
    3. 4.3 Afterglows and other post-merger signatures
      1. 4.3.1 Afterglows from circumbinary gap refilling
      2. 4.3.2 Afterglows from GW recoil kicks
      3. 4.3.3 Post-merger signatures of AGN jets
      4. 4.3.4 Long-lived accretion signatures of GW recoil
      5. 4.3.5 Stellar signatures of GW recoil
    4. 4.4 Questions uniquely addressed with detection of EM and GW signals from MBHBs
  5. 5 Conclusions

##  1 Introduction

Coincident detections of electromagnetic (EM) and gravitational wave (GW)
signals from coalescences of massive black hole binaries111Black hole binaries
with mass in the range of $10^{5}$–$10^{10}\,M_{\odot}$. (MBHBs) have the
potential to provide unparalleled understanding of the evolution of massive
black holes (MBHs) in the context of large-scale structure (Kormendy and Ho,
2013; Heckman and Best, 2014). Furthermore, since direct detection of GWs can
yield accurate luminosity distances and EM counterparts can provide redshifts
of coalescing binary systems, their combination yields fundamentally new
multimessenger probes of the Universe’s expansion (Schutz, 1986; Holz and
Hughes, 2005; Abbott et al., 2017). Coincident EM and GW detections can
measure or place limits on differences in the arrival times of photons and
gravitons from the same cosmological source, which would in turn inform us
about the mass of the graviton and possible violations of the equivalence
principle and Lorentz invariance in the gravitational sector (Kocsis et al.,
2008; Hazboun et al., 2013; Yagi and Stein, 2016; Abbott et al., 2016b).

The outcome of these scientific endeavors directly depends on our ability to
identify observed MBHB systems characterized by both messengers. At the time
of this writing such detections are anticipated but are not a foregone
conclusion. Their feasibility will be determined by the properties of MBHBs
and environments in which they reside, as well as by the technical
capabilities of EM and GW observatories. Coincident EM and GW detections can
provide… $\bullet$ Understanding of the evolution of MBHs in the context of
large-scale structure.
$\bullet$ Fundamentally new probes of the universe’s expansion. $\bullet$
Unique tests of the theory of gravity.

A basic requirement for a successful detection using any messenger is that at
least some fraction of galaxy mergers lead to the formation of close MBH pairs
that coalesce within the age of the Universe. Observations of dual and
multiple AGNs with kiloparsec separations confirm that galactic mergers are
natural sites for formation of wide MBH pairs, and in some cases multiplets.
Theoretical models and simulations of galactic mergers suggest that subsequent
evolution of widely separated MBHs to smaller scales should be common.
However, the timescales for this process remain uncertain, owing to the
inherent challenges in modeling MBHB evolution from galactic to milliparsec
scales. This, combined with a lack of concrete observational evidence for
close, gravitationally bound MBH pairs with sub-parsec separations, currently
precludes a definitive conclusion that galactic mergers result in MBHB
coalescences.

The link between merging galaxies and merging MBHs is nontrivial to establish
in observations, because the two events are separated by hundreds of millions
to billions of years and because the MBH coalescence timescale is very short
relative to all other evolutionary stages of the binary. Consequently, MBH
coalescences are difficult to identify in EM observations alone, without prior
knowledge about their occurrence. This is where a specialized messenger, such
as GWs, will play an unparalleled role in pinpointing the instance of a MBH
merger. In this case, the temporal coincidence of a prompt EM counterpart with
a gravitationally timed MBH merger may offer the best chance of identifying a
unique host galaxy. Combined together, the EM and GWs can provide crucial new
information about the link between hierarchical structure formation and MBH
growth.

Identification of a galaxy that is a host to a MBH coalescence also relies on
uniqueness of the associated EM counterpart and the ability of astronomers to
recognize it in observations. This entails a priori knowledge about the
observational appearance of the EM counterpart as well as the wavelength band
optimal for its detection. This is because there may be more than one, or
indeed many, plausible host galaxies that are EM bright in the area on the sky
associated with the GW signal.

To be useful beyond detection, the EM and GW signatures must ultimately encode
in a predictable way the information about MBHBs and environments in which
they reside. If a MBHB heading for coalescence were to be detected within the
next few years, the limiting factor in interpretation of its coincident GW and
EM signatures would likely be the uncertainties related to its EM signatures.
Indeed, the LIGO-Virgo GW community has demonstrated that mergers of Open
questions $\bullet$ Do MBHBs form and coalesce and on what timescales?
$\bullet$ Are there unique EM signatures associated with MBHB coalescences?
$\bullet$ Do they encode the properties of MBHBs and their environments?
stellar origin black hole binaries can be successfully detected and
interpreted even in the absence of any EM counterparts. The EM counterparts
associated with MBH coalescences nevertheless remain an active and quick-paced
research area with realistic potential for major breakthroughs before the
first multimessenger detection is made.

The aim of this Living Review is to provide an introduction to this research
topic by presenting a summary of key findings, physical processes and ideas
pertaining to EM counterparts to MBH mergers as they are known at the time of
this writing. We review current observational evidence for close MBHBs in
Sect. 2 and discuss the relevant physical processes and timescales in Sect. 3.
In Sect. 4 we summarize the possible precursor, coalescence, and afterglow EM
counterparts to GWs associated with MBHBs. We describe open questions, discuss
future prospects in this research area, and conclude in Sect. 5.

##  2 Observational evidence for close massive black hole binaries

Galactic mergers are an important driver of galaxy evolution and a natural
channel for formation of MBH pairs and, possibly, multiplets (Begelman et al.,
1980). Most of the scientific attention was initially focused on the galaxies
in the act of merging, and their central MBHs were considered passive
participants, taken for a ride by their hosts. The realization that MBHs play
an important role in the evolution of their host galaxies (Kormendy and
Richstone, 1995; Magorrian et al., 1998; Ferrarese and Merritt, 2000; Gebhardt
et al., 2000; Tremaine et al., 2002) spurred a number of observational and
theoretical studies, undertaken over the last few decades, with an aim to
determine what happens to MBHs once their host galaxies merge.

From an observational point of view, this question is pursued through EM
searches for dual and multiple MBHs with a variety of separations, ranging
from $\sim$ tens of kpc to sub-parsec scales (see De Rosa et al., 2019, for a
review). The multi-wavelength searches for MBH systems with large separations,
corresponding to early stages of galactic mergers, have so far successfully
identified a few dozen dual and offset active galactic nuclei (AGNs; Komossa
et al., 2003; Koss et al., 2011, 2016; Liu et al., 2013; Comerford et al.,
2015; Barrows et al., 2016, and others). MBHs with even smaller (parsec and
sub-parsec) separations are representative of the later stages of galactic
mergers in which the two MBHs are sufficiently close to form a gravitationally
bound pair.

A MBH pair is gravitationally bound if the amount of gas and stars enclosed
within its orbit is smaller than its own mass. For a wide range of MBH masses
and host galaxy properties, this happens when the two holes reach separations
of $\sim 1-10$ pc (see Sect. 3.3). Throughout this paper, we refer to such
systems as _massive black hole binaries_ and use the term “dual MBHs” for
unbound pairs of MBHs in the same general galactic vicinity. A key
characteristic of gravitationally bound MBHBs is that they are Observational
evidence for MBHBs (from more to less direct) $\bullet$ Emission of
gravitational waves.
$\bullet$ Direct imaging of double nuclei with Very Long Baseline
Interferometry. $\bullet$ Quasi-periodicity in light curves of AGNs and
quasars. $\bullet$ Doppler-shifted emission lines in spectra of AGNs and
quasars.
observationally elusive and expected to be intrinsically rare. While the
frequency of binaries is uncertain and dependent on their unknown rate of
evolution on small scales (see Sect. 3), theorists estimate that a fraction of
$<10^{-2}$ AGNs at redshift $z<0.6$ may host MBHBs (Kelley et al., 2019b),
although only a small fraction can be identified as such (see for example the
estimates in Volonteri et al., 2009).

These results have two important implications: (a) any EM search for MBHBs
must involve a large sample of AGNs, making the archival data from large
surveys of AGNs an attractive starting point and (b) observational techniques
used in the search must be able to distinguish signatures of binaries from
those of AGNs powered by single MBHs. Observational techniques used to search
for such systems in large archival datasets have so far largely relied on
direct imaging, photometry, and spectroscopic measurements. They have recently
been complemented by observations with pulsar timing arrays (PTAs), which
could detect GWs from massive MBHBs. In the remainder of this section we
summarize the outcomes of and future prospects for these different
observational approaches. We also direct the reader to Dotti et al. (2012) and
Schnittman (2013) for reviews of a broad range of MBHB signatures proposed in
the literature.

###  2.1 Emission of gravitational waves

Figure 1: Top: Sky-averaged sensitivity (shown as the upper limit on the GW
strain amplitude) as a function of the GW frequency from the NANOGrav 5-year
data set (green), 9-year data set (orange), and 11-year data set (blue). The
data is presently most sensitive at the frequency of 8 nHz and least sensitive
around $3\times 10^{-8}\,{\rm Hz}=1{\rm year}^{-1}$, corresponding to the
orbital frequency of Earth around the Sun. Any MBHBs with GW strain above the
11-year sensitivity curve would at this point be detected by NANOGrav. Figure
from Aggarwal et al. (2019). Bottom: Angular correlation of common-spectrum
pulsar timing noise versus angular separation of pulsars on the sky, from the
NANOGrav 12.5 year dataset. The violin plots (gray/black) show the
marginalized posteriors for the interpulsar spatial correlations. For
comparison, the orange dashed line shows the flat correlation signature
between pulsars expected as a signature of error in the timescale (e.g.,
drifts in clock standards), and the blue dashed curve shows the Hellings &
Downs correlation predicted to result from a stochastic background of GWs from
MBHBs. Figure from Arzoumanian et al. (2020).

Detection of gravitational radiation emitted by a MBHB heading for coalescence
constitutes the most direct evidence for the existence of such systems. At the
time of this writing, GWs from MBHBs have not yet been definitively detected.
However, expectations for such detections of gravitational radiation have been
raised by the success of the LIGO-Virgo Collaboration (e.g., Abbott et al.,
2016a), by the selection of the Laser Interferometer Space Antenna (LISA) for
a large-class mission in the European Space Agency science program, and by the
success of the LISA Pathfinder mission (Armano et al., 2016) and the TianQin-1
experimental satellite (Luo et al., 2020). Recently, three PTA experiments,
the North American Nanohertz Observatory for Gravitational Waves (NANOGrav),
the Parkes Pulsar Timing Array (PPTA) and the European Pulsar Timing Array
(EPTA) reported strong evidence for a signature of the common-spectrum,
stochastic process (Arzoumanian et al., 2020; Goncharov et al., 2021; Chen et
al., 2021). Even though the signal does not show sufficient evidence of
quadrupolar correlations needed to claim detection of a stochastic background
of GWs from MBHBs (see Fig. 1 and the text below), it provides a useful
testing ground for a variety of theoretical models, and adds to the sense that
an era of low-frequency GW astronomy is imminent.

PTAs seek to detect GWs by searching for correlations in the timing
observations of a network of millisecond pulsars (Hellings and Downs, 1983;
Foster and Backer, 1990). Specifically, the deviation of pulse arrival times
from different pulsars due to an intervening nHz GW is predicted to exhibit a
characteristic correlation with the angle between the pulsars on the sky; this
is the so-called Hellings and Downs (1983) curve. The three mentioned PTA
experiments currently in operation, NANOGrav (McLaughlin, 2013), EPTA
(Desvignes et al., 2016), and PPTA (Hobbs, 2013), together form the
International Pulsar Timing Array (IPTA; Verbiest et al., 2016). PTAs are
sensitive to GWs with frequencies between a few and few hundred nHz, and their
sensitivity increases for longer observational time baselines. This is
illustrated in Fig. 1, which shows the sensitivity as function of the GW
frequency of NANOGrav based on the 5, 9 and 11 years of timing observations of
about forty pulsars (Aggarwal et al., 2019). PTA sensitivity is also being
steadily improved as more pulsars are identified and added to the network.

The expected signals for PTAs include (i) the aforementioned stochastic GW
background, which is produced by an ensemble of unresolved MBHBs emitting GWs
at different frequencies and distances throughout the Universe, and (ii)
continuous GWs from individual MBHBs, which would exhibit negligible frequency
evolution on human timescales. At this point it is unclear whether the
stochastic GW background will be detected first (Rosado et al., 2015;
Arzoumanian et al., 2020) or whether it will be detected contemporaneously
with the individual MBHBs (Mingarelli et al., 2017; Kelley et al., 2018). In
either case, a wide range of theoretical predictions agree that GWs from MBHBs
are expected to be observed by the PTAs within the next decade. If so, their
discovery and analysis will be well timed to inform expectations for GW
detections with the space-based observatory LISA, which is planned to launch
in the mid-2030s. Because PTAs are sensitive to the most massive MBHB systems
out to $z\sim 2$ (Sesana et al., 2008; Chen et al., 2017), they stand a chance
of being associated with some of the most EM luminous AGNs in the local
Universe. This bodes well for the prospects for multimessenger studies of
MBHBs in the near future.

PTAs are sensitive to individual MBHBs with high total mass, $\sim
10^{8}$–$10^{10}\,M_{\odot}$, and orbital periods on the order of months to
years. For example, the NANOGrav collaboration performed a search in their
11-year dataset for GWs from individual MBHBs on circular orbits. While no
statistically significant evidence for GWs has been found in the data, a
physically interesting limit has been placed on the distances to, or
equivalently masses of, individual MBHBs (Aggarwal et al., 2019). The 11-year
data can rule out the existence of MBHBs within 120 Mpc with chirp mass222The
chirp mass is defined as
$\mathcal{M}\equiv(M_{1}\,M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}=q^{3/5}M/(1+q)^{6/5}$,
where $M_{1}$ and $M_{2}$ are the masses of the primary and secondary MBH,
$M=M_{1}+M_{2}$ is the binary mass, and $q=M_{2}/M_{1}\leq 1$ is the mass
ratio. Note that $\mathcal{M}\approx q^{3/5}M$ for $q\ll 1$ and
$\mathcal{M}\approx M/2$ for $q=1$. $\mathcal{M}=10^{9}\,M_{\odot}$ emitting
GWs with frequency $8$ nHz, which falls in the most sensitive part of the
spectrum (see Fig. 1). The existence of binaries with even larger chirp masses
of $\mathcal{M}=10^{10}\,M_{\odot}$, emitting at the same frequency, can be
ruled out to even higher, cosmological distances of 5.5 Gpc.

Figure 2: Contours of constant signal-to-noise ratio for the LISA observatory
plotted as a function of redshift (or equivalently, cosmic time) and MBHB
mass. The contours are calculated for systems of non-spinning MBHs with mass
ratio $q=0.5$. Overlaid is the illustration of an evolutionary track, which
ends with a formation of a $\sim 10^{7}\,M_{\odot}$ MBH merger remnant at
$z\approx 3$. Black dots and arrows represent the MBHs and their spins,
respectively. MBHs are embedded in galactic halos (illustrated by white-yellow
circles) and grow through episodes of accretion (black lines) and mergers.
Black stars mark the redshifts of the most distant quasar, long gamma-ray
burst and galaxy detected so far. Figure from Colpi et al. (2019).

In contrast to PTAs, the planned LISA and TianQin space-based missions will
detect higher GW frequencies and MBHBs with lower total masses (Amaro-Seoane
et al., 2017; Mei et al., 2021). LISA will be sensitive to GWs with
frequencies in the range from about $100\,\mu$Hz to 100 mHz. The frequency of
GW radiation emitted by MBHBs with total masses of
$10^{4}$–$10^{7}\,M_{\odot}$ falls squarely within this bandwidth in the last
stages of their evolution (Klein et al., 2016). Figure 2 illustrates that for
many binary configurations in this mass range, high signal-to-noise merger
detections will be possible out to redshift of $z\sim 20$, which will in turn
allow measurements of the masses and spins of coalescing MBHs to a few $\%$
accuracy. LISA will therefore be able to place the evolution of MBHs in the
context of important cosmic epochs that coincide with the formation of MBH
seeds (marked as Cosmic Dawn in Fig. 2), the most distant quasars
(Reionization), and the epoch of growth through accretion and mergers (Cosmic
Noon). TianQin will operate over a similar (albeit not identical) parameter
space and, like LISA, is expected to detect GWs from actual MBHB merger events
(Feng et al., 2019; Wang et al., 2019). Both in terms of the mass range and
redshift distribution, the space-based GW observatories will therefore unearth
a population of MBHBs that is different from and complementary to those
discovered by PTAs.

Figure 3: Illustration of the evolutionary tracks (multicolor) through the
LISA frequency band for three equal mass MBHBs at redshift $z=3$ with total
intrinsic masses $10^{5}$, $10^{6}$ and $10^{7}\,M_{\odot}$. The remaining
time before coalescence is marked on the tracks. Also shown are the signals
from the extreme mass ratio inspiral (EMRI), stellar origin black hole
binaries of the type discovered by the LIGO-Virgo Collaboration, and from
resolved galactic binaries (including the already known verification
binaries). The solid green curve shows LISA’s sensitivity and black dashed
line is the sensitivity with an additional confusion signal at $\sim 1$ mHz
contributed by the unresolved galactic binaries. The sensitivity curves
correspond to a configuration with six links, 2.5 million km arm length and a
mission lifetime of 4 years. Figure from the LISA L3 mission proposal (Amaro-
Seoane et al., 2017).

From all MBHBs detected with the space-based GW observatories, the population
that stands the best chance of also being detected in EM observations includes
systems at intermediate and low redshifts. Figure 3 shows the tracks for
several representative MBHBs with equal mass ratios at redshift $z=3$ as they
evolve through the LISA frequency band. MBHBs with mass
$10^{5}-10^{7}\,M_{\odot}$ merging in galaxies at $z\leq 4$ are sufficiently
loud GW sources to be detected by LISA weeks before coalescence. For systems
in this mass and redshift range, the sky localization determined from the GW
observations alone is expected to reach accuracy of about $10-100\,{\rm
deg}^{2}$ in the weeks to hours before coalescence and $\approx 0.1\,{\rm
deg}^{2}$ at merger (Haiman, 2017; Mangiagli et al., 2020). As long as these
systems are also sufficiently EM luminous and can maintain emission regions
until very close to the coalescence, this localization accuracy provides an
opportunity for subsequent identification of an EM counterpart to GW
detection. For example, if the culprit source (an AGN on the sky) exhibits
characteristic EM variability that correlates with the GW chirp, this could
lead to a convincing EM identification of the binary headed for coalescence.
Under such conditions, Dal Canton et al. (2019) find that given a fiducial
LISA detection rate of 10 mergers per year at $z<3.5$, a few detections of
modulated X-ray counterparts are possible over the nominal, 5 year duration of
the LISA mission, with an X-ray telescope with a relatively large (1 deg2)
field of view. Whether MBHBs have unique EM signatures in this or any other
stages of their evolution that can be distinguished from regular AGNs is
presently a subject of active investigation.

###  2.2 Direct imaging of double nuclei

A clear EM manifestation of a parsec-scale MBHB is an image of a binary AGN
that forms a gravitationally bound system (as opposed to an accidental
projection on the sky). A practical obstacle in the detection of such objects
arises from their small angular separation: for example, a parsec-scale binary
at a redshift of $z=0.1$ subtends an angle of only $\sim 0.5$ mas on the sky
(neglecting the projection effects). Such scales are below the angular
resolution of most astronomical observatories operating at present time,
except the radio and millimeter observatories using the Very Long Baseline
Interferometry (VLBI) technique (e.g., D’Orazio and Loeb, 2018).

Figure 4: The MBHB candidate in the radio galaxy 0402+379. Two compact radio
cores in the middle of the figure are separated by 7.3 pc on the sky and are
thought to harbor MBHs. Both cores are characterized by flat radio spectra and
have been identified as possible AGNs based on this signature. Radio lobes
from a jet associated with one of the cores are also apparent. The color
illustrates the spectral index distribution between 8 and 22 GHz from the VLBA
observations. The inset panels show flux density as a function of the
logarithm of frequency for each core. Figure from Rodriguez et al. (2006).

This approach has been used to identify the most convincing MBHB candidate
thus far, in the radio galaxy 0402+379 (Rodriguez et al., 2006, 2009; Morganti
et al., 2009). This system hosts a pair of two compact radio cores at a
projected separation of 7.3 pc on the sky (Fig. 4). Both cores are
characterized by flat radio spectra and have been identified as possible AGNs
based on this signature (Maness et al., 2004; Rodriguez et al., 2006). Bansal
et al. (2017) subsequently reported that long term VLBI observations reveal
relative motion of the two cores, consistent with orbital motion, lending
further support to the MBHB hypothesis. While the MBHB candidate in the galaxy
0402+379 was discovered serendipitously, it demonstrated the power of radio
interferometry in imaging of small-separation MBHs.

In a subsequent investigation, Burke-Spolaor (2011) searched for binaries in
the archival VLBI data. The search targeted spatially resolved, double radio-
emitting nuclei with a wide range of orbital separations ($\sim 3\,{\rm
pc}-5\,$kpc) among 3114 radio-luminous AGNs in the redshift range $0<z\leq
4.715$. Another investigation, in addition to MBHBs, searched for the
recoiling MBHs spatially offset from the centers of their host galaxies
(Condon et al., 2011). The latter study is based on Very Long Baseline Array
(VLBA) observations at a frequency of 8GHz of a sample of 834 nearby radio-
luminous AGNs with typical distances of $\sim 200$ Mpc. Radio sources brighter
than 100 mJy were selected from the 1.4GHz NVSS catalogue data, under the
hypothesis that the most massive galaxies have undergone a major merger in
their history and may therefore contain a binary or recoiling MBH. More
recently, Tremblay et al. (2016) performed multi-frequency imaging of a flux-
limited sample of $\sim$1100 AGN included in the VLBA Imaging and Polarimetry
Survey, all of which were pre-selected to be bright and flat-spectrum. None of
these searches unearthed new instances of double radio nuclei, leaving several
possible interpretations: (a) there is a true paucity of MBHBs, (b) MBHBs may
be present but have low radio brightness, or (c) only one component of the
binary is radio-bright and it may or may not show a detectable spatial offset
relative to the center of the host galaxy. This highlights the difficulties in
using radio imaging as the primary technique for selecting gravitationally
bound MBHB candidates or their progenitors, given their unknown radio
properties.

Figure 5: Projected nuclear separation as a function of redshift for the
observed population of candidate AGN pairs and binaries, compared with the
spatial resolution of the proposed ngVLA observatory. Orbital periods
associated with each nuclear separation have been calculated neglecting the
orbital projection effects. Each symbol represents a candidate discovered in
observations in the radio (circles), optical / near-IR (dots), or X-ray
($\times$) part of the spectrum. Horizontal lines mark approximate
evolutionary stages of MBH pairs and binaries. The black curve marks the
spatial resolution limit at 10 GHz of the VLA observatory in the
A-configuration. The red and purple curves indicate the spatial resolution
limit in the 10 GHz and 120 GHz bands, respectively, for the ngVLA 150 km
baseline array (dashed) and 1000 km extended-baseline array with VLBI
capabilities (solid). Longer baselines and higher frequencies have the
potential to resolve MBHBs in the gravitational wave regime in the nearby
universe (indicated by the curves that cross the dark yellow region). Figure
from Burke-Spolaor et al. (2018).

Notwithstanding the challenges, VLBI imaging still represents a powerful
technique. One promising avenue for VLBI may be to test candidate MBHBs
selected by other techniques, especially if longer baseline arrays operating
at higher frequencies become available in the future. Figure 5 summarizes the
capabilities of one such array, proposed as an extension of the Very Large
Array (VLA), and named the Next Generation VLA (ngVLA; Murphy, 2018; Burke-
Spolaor et al., 2018). The ngVLA could resolve gravitationally bound binaries
at $z\lesssim 0.1$ when operating at 10 GHz, with the main array baseline of
1000 km and subarrays extending to $\sim 9000$ km used in VLBI mode. The same
array configuration operating at 120 GHz could in principle resolve MBHBs at
even higher redshifts, as well as those that are inspiraling due to the
emission of GWs at $z\lesssim 0.01$. Therefore, VLBI searches with longer
baselines and higher frequencies have the potential to enable both multi-
wavelength and multimessenger searches for MBHBs at a variety of separations.
If funded, the construction of the ngVLA could commence as early as 2024 with
full operations starting in 2034, enabling it to operate contemporaneously
with another next generation radio observatory currently under development in
the Southern hemisphere, the Square Kilometer
Array333https://www.skatelescope.org/technical/info-sheets, as well as with
the low-frequency GW detectors described in the previous section.

###  2.3 Photometric measurements of quasi-periodic variability

MBHBs at $\sim$ mpc-scale separations have orbital periods of $\sim$ years,
making any associated EM variability accessible on human timescales. Many MBHB
candidates have been identified via signatures of periodic or nearly periodic
variability in the lightcurves of quasars, which are interpreted as a
manifestation of binary orbital motion. The possible physical origins of
periodic EM variability from MBHBs are discussed in Sect. 4.1.1. A well-known
example of a MBHB candidate in this category is the blazar OJ 287, which
exhibits outburst activity in its optical light curve with an observed period
close to 12 years (see Fig. 6 from Valtonen et al., 2008). It is worth noting
however that OJ 287 is unique among photometrically selected binary candidates
because the first recorded data points in its light curve extend as far back
as the 19th century. Similarly, OJ 287 has received an unprecedented level of
observational coverage in modern times, yielding a high frequency of sampling
(from 1970s onwards), in addition to long observational baseline. While the
presence of quasi-periodicity has been claimed in other objects, it is usually
less pronounced and recorded over much shorter time span than in the case of
OJ 287, thus preventing an ironclad case for existence of MBHBs from being
made (for e.g., Fan et al., 1998; Rieger and Mannheim, 2000; De Paolis et al.,
2002; Liu et al., 2014a). In 2019, a bright flare was observed from OJ 287
that occurred within 4 hours of its predicted arrival time, which provides
further support for the MBHB hypothesis for this source (Laine et al., 2020).

Figure 6: Historical light curve of OJ 287 in V-magnitude recorded over more
than 100 years. Quasi-periodic outbursts with an observed period of $\sim 12$
yr in the light curve of this object have been interpreted as a signature of
the MBHB orbital motion. Figure from Valtonen et al. (2008).

Recently, there have been a number of systematic searches for MBHBs in
photometric surveys, such as the Catalina Real-Time Transient Survey, the
Palomar Transient Factory and others (e.g., Graham et al., 2015a; Charisi et
al., 2016; Liu et al., 2016). Although some candidates have been identified,
there is an ongoing discussion about the level of evidence provided by the
data. For example, as pointed out by Vaughan et al. (2016) and later
demonstrated by Liu et al. (2018), the risk of false binary identifications in
photometric surveys is especially high when the evidence is based on only a
few apparent orbital cycles. If the current candidate fraction is treated as
an upper limit, then the incidence of MBHBs among quasars is $<5\times
10^{-4}$, which is consistent with theoretical predictions for low redshift
($<10^{-2}$; Volonteri et al., 2009; Kelley et al., 2019b).

Because of the finite temporal extent of the surveys, which must record at
least several orbital cycles of a candidate binary, most photometrically
identified MBHB candidates have relatively short orbital periods, of the order
a few years or less. While systematic searches like these signal more
efficient approach to finding MBHB candidates in surveys, a confirmation of
their identity requires additional evidence. This is because stochastically
variable lightcurves of “normal” AGNs and quasars, powered by single MBHs, can
be mistaken for periodic sources.

Photometric searches for MBHBs will expand greatly in the next decade, as the
capabilities of the optical time-domain surveys will be dramatically enhanced
with the Zwicky Transient Facility (ZTF) and the advent of the Legacy Survey
of Space and Time (LSST) with the Vera Rubin Observatory in 2023. Over time,
photometric searches for MBHBs will also become increasingly sensitive to a
wider range of orbital periods, as the existing surveys continue to collect
data on longer timescales. For example, Kelley et al. (2019b) find that if the
host galaxy does not outshine the AGN, up to $\sim 100$ such MBHB candidates
could be identifiable after 5 years of LSST observations with the Rubin
Observatory. They also suggest existence of a “sweet spot” for binary mass
ratios of $q\sim 0.1$, where differential accretion onto the secondary MBH
could enhance variability. This could potentially improve the prospects for
observing some MBHB systems with lower total mass. Xin and Haiman (2021)
furthermore predict that from all quasi-periodic AGNs detected by LSST (of
which there may be 20 to 100 million), there will be $10-150$ powered by MBHBs
with short orbital periods ($\lesssim 1$ day) that could in $\sim 5-15$ years
be detected by LISA. Such sources would thus serve as “LISA verification
binaries”, analogous to short-period Galactic compact-object binaries.

Already, some multimessenger constraints have been placed on the MBHB
population by combining these photometrically variable candidates with upper
limits from PTAs on their GW signals. PTAs are capable of probing the GW
background at nanoHertz frequencies from sources such as MBHBs with orbital
periods of a few years (Lentati et al., 2015; Shannon et al., 2015;
Arzoumanian et al., 2016, 2020). Specifically, the upper limit placed by the
sensitivity of PTAs largely rules out the amplitude of GW background resulting
from the so far identified $\sim 150$ photometric binary candidates, implying
that some fraction of them are unlikely to be MBHBs (Sesana et al., 2018).
Subsequently, a similar approach was used to place limits on the presence of
MBHBs in periodic blazars (Holgado et al., 2018) and in ultraluminous infrared
galaxies (Inayoshi et al., 2018). More recently, Liu and Vigeland (2021)
illustrated how detection and parameter estimation of individual MBHBs in PTA
data could be improved by up to an order of magnitude with the addition of EM
priors. These studies provide examples of the great potential of
multimessenger techniques, even prior to a low-frequency GW detection.

###  2.4 Spectroscopic measurements of offset broad emission lines

The principal assumption made by all spectroscopic searches is that some
fraction of MBHBs at sub-parsec orbital separations reside in emission regions
comparable in size to the broad line regions (BLRs) of regular AGNs444The term
“regular AGNs” refers to the AGNs that do not host MBHBs.. If so, the
dynamical perturbation of a BLR by the binary gravitational potential can in
principle be reflected in the low-ionization broad emission-line profiles
(e.g., Gaskell, 1983, 1996; Bogdanović et al., 2008; Montuori et al., 2011,
2012). The broad emission lines of particular interest are H$\alpha~{}\lambda
6563$, H$\beta\;\lambda 4861$ and Mg ii$\;\lambda 2798$, because they are
prominent in AGN spectra and are commonly used as tracers of dense, low-
ionization gas in BLRs at low (H$\alpha$ at $z<0.4$) and high redshift (Mg ii
at $z<2.5$).

Spectroscopic searches rely on the detection of the Doppler shift of broad
emission lines in the spectrum of a MBHB host caused by the binary orbital
motion. This approach is reminiscent of a well established technique for
detection of single- and double-line spectroscopic binary stars. In both
classes of spectroscopic binaries, the lines are expected to oscillate about
their local rest frame wavelength on the orbital time scale of a system. In
the context of the MBHB model, the spectral emission lines are assumed to be
associated with gas accretion disks that are gravitationally bound to the
individual black holes (Gaskell, 1983, 1996; Bogdanović et al., 2009a). Given
the high velocity of the bound gas, the emission-line profiles from the MBH
disks are expected to be Doppler-broadened, similar to the emission lines
originating in the BLRs of regular AGNs. Moreover, several theoretical studies
have shown that in unequal-mass binaries, accretion occurs preferentially onto
the lower mass object (Artymowicz and Lubow, 1994; Günther and Kley, 2002;
Hayasaki et al., 2007; Roedig et al., 2011; Farris et al., 2014), rendering
the secondary AGN potentially more luminous than the primary. If so, this
indicates that some fraction of MBHBs may appear as single-line spectroscopic
binaries. This scenario is illustrated in Fig. 7, where the Doppler-shifted
broad component of the H$\beta$ emission line traces the motion of the smaller
MBH. (See however Nguyen et al., 2020a, who suggest that the motion traced may
actually be that of the primary).

Figure 7: Illustration of the spectroscopic radial velocity technique used to
search for MBHBs. In the context of the MBHB model adopted by spectroscopic
surveys, the broad profile (shown in the top panels as the H$\beta$ line) is
attributed to the emission from the accretion disk of the smaller MBH. The
offset of the broad emission-line profile with respect to its rest frame
wavelength (marked by the vertical dotted line) is attributed to the binary
orbital motion. The center of mass of the MBHB is assumed to be at rest with
respect to the host galaxy, and the rest frame of the host galaxy is defined
by the narrow H$\beta$ and [O iii] lines. Figure from Guo et al. (2019).

This realization led to the discovery of a number of MBHB candidates based on
the criterion that the culprit sources exhibit broad optical lines offset with
respect to the rest frame of the host galaxy (Bogdanović et al., 2009b; Dotti
et al., 2009a; Boroson and Lauer, 2009; Tang and Grindlay, 2009; Decarli et
al., 2010a; Barrows et al., 2011; Tsalmantza et al., 2011; Eracleous et al.,
2012; Tsai et al., 2013)555In an alternative approach, anomalous line ratios
have been used to flag MBHB candidates with perturbed BLRs (Montuori et al.,
2011, 2012).. Because this effect is also expected to arise in the case of a
recoiling MBH receding from its host galaxy, the same approach has been used
to flag candidates of that type (Komossa et al., 2008; Shields et al., 2009b;
Civano et al., 2010; Robinson et al., 2010; Lusso et al., 2014). The key
advantage of the method is its simplicity, as the spectra that exhibit
emission lines shifted relative to the galaxy rest frame are fairly
straightforward to select from large archival data sets, such as the Sloan
Digital Sky Survey (SDSS). Its main complication, however, is that the
Doppler-shift signature is not unique to these two physical scenarios; for
example, broad line offsets can also be produced by AGN outflows.
Complementary observations are therefore required in order to determine the
nature of the observed candidates (e.g., Popović, 2012; Barth et al., 2015).

To address this ambiguity, a new generation of time-domain spectroscopic
searches has been designed to monitor the offset of the broad emission-line
profiles over multiple epochs and identify sources in which modulations in the
offset are consistent with binary orbital motion (Bon et al., 2012, 2016;
Eracleous et al., 2012; Decarli et al., 2013; Ju et al., 2013; Liu et al.,
2014b; Shen et al., 2013; Runnoe et al., 2015, 2017; Li et al., 2016; Wang et
al., 2017; Guo et al., 2019). For example, Eracleous et al. (2012) searched
for $z<0.7$ SDSS AGNs with broad H$\beta$ lines offset from the rest frame of
the host galaxy by $\gtrsim 1000\,{\rm km\,s^{-1}}$. Based on this criterion,
they selected 88 quasars for observational follow-up from an initial group of
about 16,000. The follow-up observations of this campaign span a temporal
baseline from a few weeks to 12 years in the observer’s frame. After multiple
epochs of follow-up, statistically significant changes in the velocity offset
have been measured for 29/88 candidates (Runnoe et al., 2015, 2017), all of
which remain viable as MBHB candidates. This implies that $\lesssim
10^{-3}-10^{-2}$ of quasars host MBHBs at $z<0.7$, a limit that is so far
consistent with theoretical predictions (Volonteri et al., 2009; Kelley et
al., 2019b).

Figure 8: H$\beta$ emission-line profiles of MBHB candidates selected in a
spectroscopic survey by Eracleous et al. (2012). Only the broad emission-line
profiles are shown; the narrow line components were removed for clarity. Peaks
of the asymmetric profiles are offset from the rest frame of the host galaxy
(marked by the dashed vertical line) by $\sim{\rm few}\times 10^{3}\;{\rm
km\,s^{-1}}$. If measured velocity offsets are associated with the orbital
motion of the secondary MBH, they imply orbital periods in the range of
$\sim{\rm few}\times 10-{\rm few}\times 100$ yr for MBH mass ratios larger
than 0.1. The dotted lines mark velocity offsets of $\pm 5000\;{\rm
km\,s^{-1}}$.

Figure 8 shows several representative broad $H\beta$ profiles selected in this
search. The profiles are asymmetric and have peaks offset by $\sim{\rm
few}\times 10^{3}\;{\rm km\,s^{-1}}$ from the rest frame of the galaxy,
inferred from the wavelength of the narrow emission lines (the narrow
emission-line components were subsequently removed from these profiles for
clarity). If measured velocity offsets of this magnitude are associated with
the orbital motion of the secondary MBH, as described above, then it is
possible to show that searches of this type (with a yearly cadence of
observations) are in principle sensitive to a subset of MBHBs with orbital
separations $\lesssim{\rm few}\times 10^{4}$ Schwarzschild radii and orbital
periods in the range of a few tens to a few hundreds of years (Pflueger et
al., 2018). In comparison, current spectroscopic observational campaigns
typically span a baseline about 10–15 years, shorter than the average period
of MBHBs targeted by this technique.

These relatively short monitoring baselines still allow one to reject MBHB
candidates with velocity curves inconsistent with the MBHB model, even if they
have been monitored for less than a full orbital cycle (Runnoe et al., 2017;
Guo et al., 2019). At the same time, they preclude measurement of multiple
orbital cycles, which have traditionally been used as a criterion for binarity
in stellar systems. Therefore, this technique alone cannot uniquely confirm
the identity of a MBHB; it must be supplemented by additional, independent
observational evidence that can help to further elucidate their nature. Some
of the more promising complementary approaches include direct VLBI imaging at
millimeter and radio wavelengths of nearby MBHB candidates ($z\lesssim 0.1$)
with orbital separations $a\gtrsim 0.01\,$pc (e.g., D’Orazio and Loeb, 2018;
Burke-Spolaor et al., 2018; Breiding et al., 2021).

As in the case of the photometric searches, a statistical test of the nature
of spectroscopic MBHB candidates can be performed using GW constraints from
PTAs. While hypothetical MBHBs detected spectroscopically do not emit
significant amounts of GWs because they are at relatively large separations,
they imply some number of MBHBs inspiraling toward coalescence, whose GW
signal is reaching Earth at this very moment. For example, a comparison of
this GW signal with the PTA sensitivity limit indicates that all the MBHB
candidates in the spectroscopic sample of Eracleous et al. (2012) could be
binaries without violating the observational constraints on the GW background
(Nguyen et al., 2020b).

The risk of binary misidentifications in spectroscopic surveys can also be
mitigated by continued monitoring and longer observational baselines. One
example of a survey that can provide such opportunity in the next decade is
SDSS-V, which is planned for the period
2020–2025666https://www.sdss.org/future. As a part of SDSS-V, the Black Hole
Mapper multi-object spectroscopic survey will provide multi-epoch optical
spectra for more than 400,000 X-ray sources to probe the evolution of MBHs
over cosmic time. This and similar surveys will therefore provide a
considerable sample of active galaxies that are potential hosts to MBHBs.

##  3 Physical processes and relevant timescales

###  3.1 Rates and typical redshifts

In the standard model of structure formation in the Universe, bound halos of
dark matter (along with their associated baryons) assemble in a hierarchical
fashion, with small structures forming first and combining to make larger
structures. If halos contain central MBHs, then halo mergers can lead to
eventual mergers of the black holes. Such a coalescence will typically not
happen if one halo is much less massive than the other, because in that case
dynamical friction is inefficient and the lower-mass halo will either have an
inspiral time longer than the Hubble time, or it will dissolve before it
settles in the center. Indeed, Taffoni et al. (2003) find that below mass
ratios of 1:10, merging of dark matter halos is ineffective. As we discuss in
the next section, a similar result is found in simulations of isolated (i.e.,
non-cosmological) galaxy mergers and cosmological simulations that follow
pairing of MBHs.

Analyses of halo mergers (Press and Schechter, 1974; Lacey and Cole, 1993;
Sheth et al., 2001; Menou et al., 2001; Volonteri et al., 2002, 2003; Wyithe
and Loeb, 2003; Sesana et al., 2004; Yoo and Miralda-Escudé, 2004; Islam et
al., 2004a, b; Shapiro, 2005; Sesana et al., 2005; Micic et al., 2007, and
others) indicate that most of the halo mergers that might lead to MBH mergers
occur at redshifts of several. One therefore finds that, if they follow
promptly after their host halos merge, most MBH mergers would happen at
redshifts $z>3$ and perhaps at Physical mechanisms leading to orbital
evolution of MBHBs $\bullet$ Dynamical friction (during pairing of MBHBs).
$\bullet$ Stellar “loss-cone” scattering. $\bullet$ Interactions with
additional MBHs. $\bullet$ Torques from a circumbinary gas disk. $\bullet$
Emission of gravitational waves.
even higher redshifts. Such mergers would be difficult to detect
electromagnetically because of their distance and the weakening of the signal
due to redshift, and it will be even more difficult to identify their host
galaxies. There is however evidence that delays due to orbital evolution of
MBHs may shift the peak of the MBH coalescence rate to $z\sim 1-2$, a much
more promising outcome for the EM followup (Volonteri et al., 2020; Li et al.,
2022). The space-based gravitational wave detectors (like LISA and TianQin)
will nonetheless be able to detect mergers out to $z\approx 20$ and will
therefore provide unique insight into the era of structure assembly.

As discussed in Sect. 2.1, EM counterparts to MBH mergers would be much easier
to detect at redshifts $z<1$. The rate of such mergers, which may be
accompanied by EM counterparts, is difficult to determine directly, but if we
assume that major galaxy mergers (with mass ratios $>1/4$) inevitably lead to
MBH mergers, then because galaxies comparable to or more massive than the
Milky Way have typically had at least one major merger since $z\sim 1$ (Bell
et al., 2006; Lotz et al., 2008; López-Sanjuan et al., 2009; Bundy et al.,
2009; Hopkins et al., 2010), and there are tens of billions of such galaxies
in the universe, one might expect a few such mergers to happen per year.

In the last few decades it has become clear that the mass of a central MBH is
correlated with properties of its host galaxy (e.g., Kormendy and Richstone,
1995; Magorrian et al., 1998; Ferrarese and Merritt, 2000; Gebhardt et al.,
2000; Merritt and Ferrarese, 2001; Tremaine et al., 2002; Ferrarese, 2002;
Ferrarese et al., 2006; Gültekin et al., 2009; McConnell and Ma, 2013).
Perhaps the tightest correlation is between the black hole mass, $M$, and the
stellar velocity dispersion of the bulge, $\sigma$. A recent assessment of
this correlation (McConnell and Ma, 2013) gives

$\log(M/M_{\odot})=\alpha+\beta\log(\sigma/200~{}{\rm km~{}s}^{-1})\,,$ (1)

with $\alpha=8.32\pm 0.05$, $\beta=5.64\pm 0.32$, and there is an intrinsic
scatter in the logarithm

of the MBH mass of $\epsilon_{0}=0.38$. This correlation implies that moderate
to massive galaxies typically contain black holes with masses $M\sim
10^{6}$–$10^{8}\,M_{\odot}$. Mergers between such black holes release energy
in GWs that is comparable to the total energy emitted by the stars in the
galaxy over a Hubble time. As we discuss later in this section, some of that
energy is emitted asymmetrically and hence, the remnant can receive a kick
sufficient to eject it from the merged host galaxy. Characteristic timescales
$\bullet$ Galaxies merge in a few Gyr.
$\bullet$ Bound MBHB forms in $\sim 10^{7-9}\,$yr due to dynamical friction.
$\bullet$ MBHB shrinks in $\sim 10^{6-9}\,$yr due to interactions with stars
and possibly circumbinary disk. $\bullet$ On the smallest scales MBHB
coalesces rapidly due to the emission of GWs.

There are several important stages between a galaxy merger and a black hole
merger, defined in terms of the physical mechanisms that determine the orbital
evolution of the MBH pair. In brief:

  1. 1.

The interacting galaxies settle together in a few orbital times. The bulges of
the galaxies, which contain the MBHs, and are several hundred times more
massive than the holes, are initially dense enough to act dynamically as
single objects. They therefore sink in the overall star-gas system of the
remnant galaxy. This can take a few Gigayears.

  2. 2.

When the cores get within several hundred parsecs of each other, they interact
dynamically and the MBHs shed their stars. The holes now spiral separately in
the potential of the remnant galaxy, and their orbit evolves due to dynamical
friction against the stars and gas. This stage can take between $\sim 10^{7}$
to $10^{9}\,$yr.

  3. 3.

When the MBHs get close enough so that their total mass significantly exceeds
the mass of the gas and stars within their orbits (typically, $\sim{\rm
few}\,$pc), their orbital evolution can slow, or even stall in principle, if
much of the interacting stars and gas are ejected. As a consequence, the MBHB
evolution in this stage is a sensitive function of the properties of its
environment and spans a wide range of timescales ($\sim 10^{6}-10^{9}\,$yr).

  4. 4.

If the holes can bridge this gap and reach separations of $10^{-2}-10^{-3}$
pc, then gravitational radiation by itself suffices to bring the MBHs to
coalescence in less than a Hubble time. This phase can be very short, compared
to the the other evolutionary stages mentioned above.

The last three evolutionary stages are illustrated in Fig. 9. We now discuss
them in more detail, in the context of stellar and gas dynamical processes.

Figure 9: Diagram of timescales associated with different stages of MBHB
evolution: dynamical friction ($t_{\rm df}$), loss-cone scattering ($t_{\rm
h}$), gas processes ($t_{\rm gas}$), and GW emission ($t_{\rm GR}$). The
timescales correspond to a relatively massive binary, $M\approx
10^{8}\,M_{\odot}$, residing in a massive elliptical galaxy with stellar
velocity dispersion $\sigma=300\,{\rm km\,s^{-1}}$. Unless the loss cone can
be replenished on a timescale much less than the stellar relaxation time
(dotted line), the MBHB orbital evolution takes longer than a Hubble time
(solid). Gas dynamical processes can in principle lead to faster MBHB
evolution (short dashed). For binaries that reach smaller separations,
gravitational radiation takes over as the dominant mechanism at $\sim
10^{-2}\,$pc (long dashed). Figure from Begelman et al. (1980).

###  3.2 Effects before MBHs are bound to each other

####  3.2.1 Stellar and gas dynamical friction

Dynamical friction is thought to be the principal mechanism responsible for
orbital evolution of MBHs in the aftermath of galactic mergers and an
important channel for formation of gravitationally bound MBH binaries. When
two galaxies merge, their MBHs are initially tightly bound to their individual
stellar bulges, so that the bulge-MBH system dynamically acts as one object.
The parent bulges continue to sink towards each other until they reach $\sim
1\,$kpc. Beyond this point the bulges merge and gradually dissipate until the
MBHs are on their own. As long as the mass of the stars and gas between the
MBHs is significantly larger than the total mass of the holes, the MBHs are
dragged towards each other by dynamical friction exerted by the background
stars and gas in the newly-formed remnant galaxy. Thus, for a wide range of
galactic merger scenarios, this process is expected to play an important role
for MBHs with separations in the range from $\sim{\rm few}\times 100\,$pc to
$\sim{\rm few}\,$pc, when they become gravitationally bound. Dynamical
friction is also expected to be important for MBHs that are ejected to
comparable distances in dynamical interactions with other MBHBs, or in the
aftermath of the GW recoil, the effects that we discuss in Sect. 3.3.1 and
Sect. 3.4.

In this section we summarize the key aspects of dynamical friction in stellar
and gaseous media, as applicable to MBHs in those scenarios.777We do not
consider dynamical friction due to dark matter, which may play an important
role for dynamics of intermediate mass black holes ($M\sim
10^{3-5}\,M_{\odot}$) in mergers of dwarf and dark-matter dominated galaxies
(see for example, Di Cintio et al., 2017; Tamfal et al., 2018). For a general
introduction to this topic, we direct the reader to the seminal work of
Chandrasekhar (1943), who developed the analytic theory for gravitational drag
in collisionless systems, and to Ostriker (1999), who performed an equivalent
treatment for a gaseous medium.

Dynamical friction arises when a massive perturber moving through some
background (stellar or gaseous) medium creates in it a density wake. The wake
trails the perturber and exerts on it gravitational force in the direction
opposite to the motion, thus acting as a brake and earning this interaction a
name: “dynamical friction.” Consider, for example, an encounter of a MBH with
mass $M$ moving at a velocity ${\bf v}$ relative to an isotropic background of
stars with much lower mass, and a total mass density $\rho$. If the stars have
a three-dimensional Maxwellian velocity dispersion $\sigma$, the drag
experienced by the MBH is (e.g., Binney and Tremaine, 1987)

${d{\bf v}\over{dt}}=-{4\pi\ln\Lambda\,G^{2}M\rho\over{v^{3}}}\left[{\rm
erf}(X)-{2X\over{\sqrt{\pi}}}e^{-X^{2}}\right]{\bf v}$ (2)

where $X\equiv v/(\sqrt{2}\sigma)$, $\ln\Lambda=\ln(b_{\rm max}/b_{\rm
min})\sim 10-20$ is the Coulomb logarithm, and $b_{\rm min}$ and $b_{\rm max}$
are the smallest and largest effective impact parameter, respectively. The
characteristic timescale for the MBH to decelerate due to multiple encounters
with stars is then

$t_{\rm df}=v/(dv/dt)\sim v^{3}/(4\pi G^{2}M\rho)\;.$ (3)

An important implication of this expression is that MBHs of larger mass,
immersed in denser environments, can sink faster toward the center of the
remnant galaxy.

The exact timescale for decay of the MBH’s orbit in the remnant stellar bulge
will depend on the properties of the stellar bulge, such as the distribution
of stars and the stellar velocity dispersion $\sigma$, which provides a
characteristic speed of stars relative to the MBH. For example, following
Antonini and Merritt (2012, see their Eq. 19), who calculate the effect of
stellar dynamical friction in nuclear star clusters around MBHs, we estimate
the inspiral time for a $10^{6}\,M_{\odot}$ MBH initially located at
$r=100\,$pc from the cluster center

$t_{\rm df}\approx 45\,{\rm Myr}\,\left(\frac{r}{100\,{\rm
pc}}\right)^{2}\left(\frac{\sigma}{100\,{\rm
km\,s^{-1}}}\right)\left(\frac{M}{10^{6}\,M_{\odot}}\right)^{-1}\left(\frac{\ln\Lambda}{15}\right)^{-1}\,.$
(4)

Note that in the fiducial example considered in Eq. 4 inspiral from 1 kpc due
to stellar dynamical friction takes 4.5 Gyr and is shorter for more massive
black holes. In the estimate above, the stellar background is described by the
mass density profile $\rho(r)=\rho_{0}(r/1\,{\rm pc})^{-1.8}$, where
$\rho_{0}=1.5\times 10^{5}\,M_{\odot}\,{\rm pc^{-3}}$ is the normalization.
This profile is an extrapolation of the observed distribution of stars in the
central 10 pc of the Galactic Center (e.g., Oh et al., 2009). Stellar density
profiles that decrease less rapidly with radius would result in even shorter
timescales. See for example an estimate in Eq. (5) in Souza Lima et al.
(2017), who instead assume a stellar distribution consistent with $\rho\propto
r^{-1}$ (i.e., the Hernquist bulge; Hernquist, 1990).

Figure 10: Magnitude of the dynamical friction force, as a function of the
Mach number, $\mathcal{M}=v/c_{s}$. The force is exerted either by the gas
(solid lines) or by the collisionless background (dashed) of the same density,
and with a sound speed of the gas equal to the velocity dispersion of the
collisionless system. From top to bottom, the lines correspond to the
dynamical friction force contributed by regions of increasing size around the
perturber. Figure from Ostriker (1999).

If the remnant galaxy contains gas, the gas can also play an important role
for pairing of MBHs. For example, Ostriker (1999) evaluated the gravitational
drag exerted by a uniform-density gas background on a massive perturber on a
linear trajectory888See also Kim and Kim (2007) for calculation with a massive
perturber on a circular orbit.. Figure 10, adapted from this work, shows a
comparison of the magnitude of the dynamical friction force, exerted by the
gas with the same density as a purely collisionless background, and with a
sound speed $c_{s}$, equal to the velocity dispersion $\sigma$ of the
collisionless system. Ostriker (1999) found the following results as a
function of the Mach number, ${\cal M}\equiv v/c_{s}$, where $v$ is the speed
of the massive object relative to the local standard of rest of the gas:

  1. 1.

For ${\cal M}\gg 1$, the gas drag is the same as the drag against a
collisionless background, for the same average density $\rho$ and a velocity
dispersion $\sigma$ of the collisionless background equal to $c_{s}$.

  2. 2.

For ${\cal M}<1$, the drag force is larger in a collisionless medium than in a
gaseous medium. This is because the gas pressure forces result in a more
symmetric distribution around the massive object.

  3. 3.

For ${\cal M}\sim 1$, the dynamical friction force due to the gas can exceed
the collisionless drag by a factor of four.

####  3.2.2 Implications for formation and properties of gravitationally bound
MBHs

Because of its efficiency relative to the collisionless case, gaseous drag has
been extensively investigated in simulations that follow pairing of MBHs in
gas-rich mergers of galaxies (see Mayer, 2013, 2017, for reviews). When
pairing is successful, gaseous dynamical friction is capable of transporting
the MBHs from galactic radii of few hundred pc to $\sim 1$ pc on timescales as
short as $\sim 10^{7}$ yr (for details see, e.g., Escala et al., 2004, 2005;
Dotti et al., 2006; Mayer et al., 2007; Van Wassenhove et al., 2014), and in
many cases, tens of times faster than stellar dynamical friction (e.g.,
Berczik et al., 2006; Preto et al., 2011; Khan et al., 2011; Holley-Bockelmann
and Khan, 2015; Vasiliev et al., 2015; Pfister et al., 2017). The propensity
of gas to effectively couple with the orbiting MBHs can be understood by
recalling that unlike collections of stars, gas can radiate energy. Thus,
again unlike stars, the same parcels of gas can interact multiple times with a
black hole pair. As a result, in the extreme case that only gas is present,
the binary can in principle be induced to spiral inwards without severely
depleting the gas reservoir.

Simulations of isolated galaxy mergers, however, show that MBH pairs with mass
ratios $q<0.1$ are less likely to form gravitationally bound binaries within a
Hubble time at any redshift, primarily due to the inefficiency of dynamical
friction on the lower-mass MBH (Callegari et al., 2009, 2011). They find that
in moderately gas rich galaxies MBH pairs with initially unequal masses tend
to evolve toward equal masses, through preferential accretion onto a smaller
MBH. This trend is also consistent with that found in cosmological simulations
Illustris (Vogelsberger et al., 2014; Kelley et al., 2017), ROMULUS25 (Tremmel
et al., 2018; Ricarte et al., 2021), and a suite of cosmological zoom-in
simulations (Bellovary et al., 2019, 2021); moreover, evidence for wandering
MBHs in nearby dwarf galaxies has recently been found (Reines et al., 2020).
These works show that low-mass galaxies, which lack a dense stellar core, are
more likely to become tidally disrupted early in the process of merging. As a
consequence, they deposit their MBHs at large radii, resulting in a population
of ‘wandering’ MBHs. If so, this signals a preference for $q>0.1$ mergers of
MBHs. In the circumstance that the opposite is suggested by observations of GW
mergers, this would point to an existence of a mechanism other than dynamical
friction that can promote pairing of MBHs of disparate masses.

Another important consideration for the pairing timescale of MBHs is the
overall rotation of the stars and gas in the merged galaxy. For example, if a
dominant component of overall speed of stars is in bulk motion, such as
rotation, the velocity dispersion of the stars is then less than if their
motion were random. Equation (2) indicates that the impact of dynamical
friction (and more generally, dynamical relaxation) is stronger in systems
with lower velocity dispersion $\sigma$. Studies of rotating systems confirm
their fast evolution (Einsel and Spurzem, 1999; Kim et al., 2002, 2004;
Berentzen et al., 2009; Amaro-Seoane et al., 2010a). Consequently, if the
collision of galaxies that leads to the formation of dual MBHs also imparts
significant rotation, as is likely in merger remnants, than the orbital
evolution time of MBHs is much shorter than in non-rotating systems. Holley-
Bockelmann and Khan (2015) for example find that for all other things being
the same, rotating stellar bulges drive MBHs more efficiently through the
dynamical friction and three-body scattering phase, resulting in timescales
that are between 3 (for corotating systems) and 30 times shorter (for
counterrotating systems) than the non-rotating bulges.

The latter timescale, obtained for systems in which a majority of stars are
rotating in the opposite sense from the MBH pair, stems from the ability of
counterrotating stellar orbits to effectively sap the orbital angular momentum
of the MBHs. Therefore, the same effect should produce high eccentricity MBH
pairs and bound binaries. This is important because if the eccentricity can
truly increase arbitrarily close to unity, then the MBHB can get to a point
where gravitational radiation is important and it will merge. This expectation
is confirmed by studies that follow the evolution of binary eccentricity with
respect to the host galaxy and find low eccentricities ($e\lesssim 0.1$) when
the MBH binary and the galaxy are corotating, and high eccentricities
($e\lesssim 1$) and nearly radial orbits in the case of counterrotation
(Amaro-Seoane et al., 2010a; Sesana et al., 2011; Gualandris et al., 2012;
Wang et al., 2014; Holley-Bockelmann and Khan, 2015; Mirza et al., 2017;
Rasskazov and Merritt, 2017). Furthermore, these studies find that MBH pairs
and binaries, whose orbital angular momentum is initially misaligned with
respect to that of the stellar cusp, tend to realign their orbital planes with
the angular momentum of the cusp. Similar outcomes have been observed for MBHs
evolving in gaseous circumnuclear disks. Specifically, the corotation of MBHs
with the disk results in efficient circularization of the black holes’ orbits
(Dotti et al., 2007; Li et al., 2020a). In the counterrotating case, however,
gaseous dynamical friction drives the increase in MBH eccentricities until
they experience “orbital angular momentum flip” and start to corotate with the
disk before a MBHB forms (Dotti et al., 2009b).

Simulations of isolated galaxy mergers have also revealed that pairing of MBHs
may be inefficient in clumpy galaxies, such as the massive star-forming
galaxies detected at redshifts $z\sim 1-3$. This is because gravitational
scattering of MBHs by giant molecular clouds and stellar clusters can cause
them to exhibit a random walk, rather than a smooth inspiral. Consequently,
the orbital evolution of MBHs in such environments can be substantially slower
(by one or two orders of magnitude relative to the timescale quoted above), or
can stall permanently in some cases, at distances $\sim 10-10^{3}\,$pc from
the center of the remnant galaxy (Fiacconi et al., 2013; Roškar et al., 2015;
Tamburello et al., 2017). The long MBH pairing timescales obtained in these
works leave room for the remnant galaxy to experience multiple mergers in this
period of time. If so, the existence of triple, or even quadruple, MBHs in a
single merger galaxy may be common, highlighting the need to understand their
frequency of occurrence, and importance of dynamical interactions of multi-MBH
systems. For example, using the Galaxy Zoo project catalog, Darg et al. (2011)
find that a fraction of galaxy mergers involving more than two galaxies at a
time is at least 2 orders of magnitude less than that of mergers of galaxy
pairs. Liu et al. (2011) provide the incidence of multiple MBH in galaxy
mergers, by investigating systems in which constituent MBHs accrete and shine
as AGNs. They report the fraction of kpc scale AGNs pairs as $\gtrsim 3\times
10^{-3}$, and kpc scale AGN triples as $\gtrsim 5\times 10^{-5}$, both among
optically selected AGNs at low redshift.

Interestingly, gravitational interaction of the MBH with the surrounding gas
is quite local. For example, for a black hole with mass $\sim
10^{6}\,M_{\odot}$ most of the gravitational drag force is contributed by the
gas that resides within only a few parsecs from the MBH (Chapon et al., 2013).
This proximity implies that the gas dynamical friction wake can be strongly
affected, and possibly obliterated, by irradiation and feedback from an
accreting MBH. Indeed, studies of dynamical evolution of MBHs find that
dynamical friction can be significantly reduced due to the “wake evacuation”
caused by purely thermal feedback from MBHs in simulations that follow
gravitationally recoiled MBHs (see Fig. 11; Sijacki et al., 2011) and black
hole pairing (Souza Lima et al., 2017; del Valle and Volonteri, 2018). These
results bring into question the assumed efficacy of gaseous dynamical
friction. In this context, Park and Bogdanović (2017), Toyouchi et al. (2020)
and Gruzinov et al. (2020) investigated the efficiency of gaseous dynamical
friction in the presence of radiative feedback. The first two groups used
local radiation hydrodynamic simulations and found that ionizing radiation,
which emerges from the innermost parts of the MBH’s accretion flow, strongly
affects the dynamical friction wake and renders gas dynamical friction
inefficient for a range of MBH masses and gas densities. They find that MBHs
in this regime tend to experience positive net acceleration, meaning that they
speed up, contrary to the expectations for gaseous dynamical friction in the
absence of radiative feedback. This effect is more severe at the lower mass
end of the MBH spectrum which, compounded with inefficiency of the gas drag
for lower mass objects in general, suggests that $\lesssim 10^{7}\,M_{\odot}$
MBHs have fewer means to reach the centers of merged galaxies and form bound
pairs with their counterpart MBHs. Indeed, Li et al. (2020b, 2022) show a
significant suppression of pairing probability of such MBHs in merger galaxies
with gas, in the presence of radiative feedback. This is of importance because
systems of MBHs in this mass range are direct progenitors of merging binaries
targeted by the space-based GW observatory LISA.

Figure 11: Mass-weighted gas density maps from hydrodynamic simulations of
recoiling MBHs in massive gas-rich galaxies. Left (right) panel shows the
results of the simulation without (with) feedback exerted by the accreting
MBH. The scenario with feedback results in the evacuation of gas from the wake
following the MBH. Both panels show the location of the MBHs approximately 40
Myr after the recoil and the dotted lines trace their orbit from the start of
each simulation. Figure adapted from Sijacki et al. (2011).

Another important effect for the gaseous media is that when the gas accreting
onto a MBH has much greater angular momentum than the hole itself, the MBH
tends to align its spin axis with the angular momentum of the gas. This effect
was originally proposed by Bardeen and Petterson (1975), and several papers
have investigated important aspects of this process (Papaloizou and Pringle,
1983; Natarajan and Pringle, 1998; Natarajan and Armitage, 1999; Sorathia et
al., 2013a, b; Tremaine and Davis, 2014; Morales Teixeira et al., 2014;
Zhuravlev et al., 2014; Nealon et al., 2015; Polko and McKinney, 2017; Hawley
and Krolik, 2018; Liska et al., 2019). Given that the radius of the accretion
disk around an individual MBH can be thousands of times, or more, greater than
the radius of the black hole, it follows that the angular momentum of the disk
can dominate if it has more than a few percent of the mass of the MBH. An
important implication of this statement is that it is, in principle, easy for
the accretion flow to reorient the spin of the MBH. Note that this is
different than the amount of gas required to significantly change the spin
_magnitude_ of the MBH, as conservation of angular momentum in this case
requires accretion of the amount of gas comparable the original black hole
mass.

Figure 12: Density distribution of the initial ($\theta_{i}$) and final
($\theta_{f}$) angles of MBH spins with the vector of the binary orbital
angular momentum. Left (right) panels refer to MBHs embedded in cold (hot)
circumnuclear gas disks, described with a polytropic equation of state with
$\gamma=7/5$ ($\gamma=5/3$). Upper (lower) panels refer to the spin of the
secondary (primary) MBH, both with mass $M=4\times 10^{6}\,M_{\odot}$. Dark,
medium and light-grey surfaces refer to high density regions encompassing
68.3, 95.5 and 99.7 per cent of the realizations. $\theta_{i}$ has been
sampled isotropically, and the magnitudes of the dimensionless spin parameters
(not shown) have been sampled from a constant probability distribution,
between 0 and 1. Figure from Dotti et al. (2010).

A major reason that this matters astrophysically is that, as we discuss in
§3.4, mergers between holes with coaligned spins lead to moderate GW kicks,
whereas mergers between comparable-mass black holes with randomly oriented
spins can produce kicks of thousands of kilometers per second, sufficient to
eject the remnant from any galaxy. This led Bogdanović et al. (2007) to
propose that gas alignment was effective in reorienting MBH prior to their
becoming a gravitationally bound binary. Calculations by Miller and Krolik
(2013) further find that gas torques are effective in aligning the MBH spins
and their orbit with the large scale gas flow not only at large orbital
separations but also well into the final parsec, where the two form a
gravitationally bound binary. Simulations by Dotti et al. (2009b, 2010) and
Perego et al. (2009) support the spin-alignment hypothesis, with the
refinement that because the large gas disk in the inner few hundred parsecs of
a galaxy will be somewhat turbulent, accretion is not always in the same
direction and hence, rather than being exactly aligned, the spins of MBHs will
be mostly aligned. This effect is illustrated in Fig. 12, which shows the
efficiency of alignment in circumnuclear disks that correspond to more
turbulent (hot) or less turbulent (cold) gas flows. This level of alignment
reduces the probability of a large kick but does not eliminate it, and
depending on specific assumptions one every few thousand or few tens of
thousands of galaxy mergers might result in a kick exceeding 2,000 km s-1
(Lousto et al., 2012).

It is worth noting that there are also works proposing that the packets of gas
that arrive at MBHs are insufficient to reorient their spin axes. In this
picture, the black holes accrete small quantities of gas stochastically (as
opposed to from a disk with a well defined angular momentum) and as a
consequence, the majority of MBHs are expected to spin slowly (King et al.,
2005; King and Pringle, 2006, 2007; King et al., 2008; Wang et al., 2009;
Nayakshin et al., 2012; Lodato and Gerosa, 2013). The stochastic accretion
hypothesis so far appears in tension with the spins of MBHs inferred from the
X-ray observations of AGNs with relativistically broadened Fe K$\alpha$ lines
(Iwasawa et al., 1996; Fabian et al., 2002; Reynolds and Nowak, 2003;
Brenneman and Reynolds, 2006, 2009; Schmoll et al., 2009; Brenneman et al.,
2011; Reynolds et al., 2012). If the properties of these observed MBHs are
representative of the broader population in the Universe, these measurements
imply that the dimensionless spin magnitudes of most MBHs exceed $s=0.6$. If
not by the EM measurements, this hypothesis will ultimately be tested by the
LISA GW observatory, which is expected to measure the distribution of MBH
spins with exquisite precision.

###  3.3 Effects after MBH are gravitationally bound to each other

If MBH inspiral due to the stellar and gas dynamical interactions described in
previous section is efficient, it eventually results in formation of a
gravitationally bound binary. Qualitatively, a gravitationally bound MBHB
forms when the amount of gas and stars enclosed within its orbit becomes
comparable to its mass, $M$. This orbital separation is comparable to the
radius of gravitational influence of a single MBH, where the circular velocity
around the black hole equals the stellar velocity dispersion, $r_{\rm
inf}\approx GM/\sigma^{2}$. Combining this with the $M-\sigma$ relationship in
Eq. (1) (McConnell and Ma, 2013) gives

$r_{\rm inf}\approx 3.2\,{\rm pc}\,\,M_{7}^{0.645}=6.6\times
10^{6}r_{g}\,\,M_{7}^{-0.355}\;.$ (5)

where $M_{7}=M/10^{7}M_{\odot}$ and $r_{g}=GM/c^{2}\approx 0.1~{}{\rm
au}~{}(M/10^{7}~{}M_{\odot})$ is the gravitational radius. Therefore, bound
binaries form at separations of $\sim 1-10$ pc for for a wide range of MBHB
masses and host galaxy properties.

####  3.3.1 Three-body scattering of stars and interactions with additional
MBHs

The simplest picture one can consider for a MBHB in this stage is that of a
gravitationally bound binary that evolves due to interactions with individual
stars in a completely spherically symmetric potential. For stars that come
close enough to the binary, three-body interactions lead to ejection of stars
with greater energy than they had initially, a process that leads to gradual
hardening of the MBHB. In the limit where the MBHB masses are much larger than
the masses of the interacting stars, the work of Quinlan (1996) is especially
relevant. He showed that, given an unlimited supply of stars, the binary would
harden at a rate given by the dimensionless parameter

$H={\sigma\over{G\rho}}{d\over{dt}}\left(1\over a\right)\,,$ (6)

where he found $H\approx 10-20$ in the cases most relevant to us, in which the
stellar speeds at infinity are much less than the binary orbital speed. In
this case, the parameter $H$ provides a dimensionless characterization of the
rate of change of the binding energy of the binary. The associated hardening
timescale can then be estimated as

$t_{\rm h}=\left|\frac{a}{\dot{a}}\right|=\frac{\sigma}{G\rho aH}=1.5\times
10^{6}\,{\rm yr}\,\left(\frac{\sigma}{100\,{\rm
km\,s^{-1}}}\right)\left(\frac{10^{4}\,M_{\odot}\,{\rm
pc^{-3}}}{\rho}\right)\left(\frac{0.1\,{\rm
pc}}{a}\right)\left(\frac{15}{H}\right)\,,$ (7)

and it corresponds to the time that a MBHB of mass $\sim 10^{7}\,M_{\odot}$
spends at orbital separation of $a=0.1\,$pc in a nuclear stellar cluster with
velocity dispersion of $\sigma=100\,{\rm km\,s^{-1}}$ and stellar mass density
of $\rho=10^{4}\,M_{\odot}\,{\rm pc^{-3}}$. Note also that the MBHB spends
increasingly longer times at smaller orbital separations and thus, its
evolution slows down progressively.

The implication of Eq. (7) is that three-body interactions with stars can in
principle lead to very efficient orbital evolution of MBHBs, with a wide range
of masses, given an unlimited supply of stars. This assumption is reasonable
in the stages of evolution when the MBHB orbital speed is comparable to the
velocity dispersion of the stars, because stars are not completely ejected and
can return for another interaction. When the MBHB orbital speed is
significantly greater than $\sigma$ however, interacting stars are ejected. In
that scenario, the supply of stars that can interact with the MBHB is depleted
and must, ultimately, be replenished if the binary is to continue its
evolution. The timescale on which the stars are replenished therefore
determines the rate of the binary orbital evolution.

This timescale is ordinarily determined by gravitational two-body interactions
between stars in the nuclear star cluster around the MBHB, which can change
the energy of a star by the order of itself. This is referred to as the
_energy relaxation time_ (Binney and Tremaine, 1987; Spitzer, 1987) and is
directly related to the dynamical friction time shown in Eq. (4), as it
describes the same type of gravitational interaction between the two point
masses. For example, using Eq. (4) one can estimate a relaxation time of about
41 Gyr for a $M=1\,M_{\odot}$ star residing in the nuclear star cluster with
$\sigma=100\,{\rm km\,s^{-1}}$, from an initial radius $r\approx 3\,$pc,
corresponding to the influence radius of a $10^{7}\,M_{\odot}$ binary. Such a
MBHB would thus be unable to shrink its orbit much below the influence radius
by scattering of stars and would in principle have to rely on other processes,
like emission of GWs, to continue its evolution. At this point, however,
gravitational radiation driven inspiral would also take longer than a Hubble
time for such a binary, leading to a hang-up in the evolution of the MBHB.
This phenomenon, first described by Begelman et al. (1980), has been dubbed
the “final parsec problem” (Milosavljević and Merritt, 2003).

Thus, an important question for MBHBs evolving due to three-body interactions
with stars is: are there physical scenarios in which the stars are supplied at
a sufficiently high rate, so as to guarantee the delivery of the MBHBs to the
GW regime and coalescence in less than a Hubble time? The answer depends on
the orbital properties of the stars. If the interacting star is on an orbit
with eccentricity $e$, then its specific angular momentum is
$\ell=\sqrt{GMa(1-e^{2})}$, where $M$ is the mass of the central object
(either a single MBH or a MBHB). Two-body relaxation is a diffusive process,
so the time required for the angular momentum to change by the order of itself
is $(1-e^{2})$ times the energy relaxation time. Thus, orbits with high
eccentricities can change their angular momenta rapidly. The rapidity with
which stars can change their angular momenta leads to the distinction between
full and empty loss cones (for discussion of stellar dynamics in the centers
of galaxies, see, e.g., Peebles, 1972; Frank and Rees, 1976; Bahcall and Wolf,
1976, 1977; Lightman and Shapiro, 1977; Wang and Merritt, 2004; Hopman and
Alexander, 2005).

In general, the “loss cone” is the solid angle inside of which a particular
star can have an important interaction with a target object. For example, for
tidal disruption around a single MBH, the loss cone at a given distance from
the black hole is the set of angles such that if a star’s orbit is inside that
solid angle, it will be disrupted. In the case of interactions with a MBHB,
the loss cone is the solid angle such that the star will pass within $\sim 2$
binary semimajor axes of the binary, where the star can have a significant
interaction with it. If the angular momentum can change by of order itself or
more in a single orbit, the orbit is in the _full_ loss cone regime, so-called
because any star that is ejected or destroyed from such an orbit will be
replaced immediately by another star. In contrast, if the typical change in
the angular momentum in one orbit is less than the magnitude of the angular
momentum, then the orbit is in the _empty_ loss cone regime, and loss of stars
must be replaced diffusively, on the timescale discussed in the previous
paragraph. In such scenarios, there is a decrease in the phase space density
inside empty loss cones, leading to significant slow down or complete hang-up
in the orbital evolution of the MBHB. These two regimes are illustrated in
Fig. 9, which shows different evolutionary tracks for MBHBs that find
themselves in full and empty loss cones.

Given the importance of the final parsec problem for plausibility of
coalescences of MBHBs, there have been many papers devoted to its solution.
One approach has been to consider the evolution of MBHBs in non-spherical
potentials. The motivation for this has been twofold. Firstly, galaxies that
are remnants of major mergers are not expected to be spherically symmetric.
This is because the relaxation timescales in their nuclei (discussed above)
are expected to exceed a Hubble time in many relevant physical scenarios.
Secondly, if the underlying stellar potential is triaxial then, because it
lacks the azimuthal or spherical symmetry which would guarantees conservation
of angular momentum for individual orbits, individual stars can replenish the
loss cone much more rapidly than would be possible via two-body relaxation.
Indeed, for galaxies that can be described as triaxial ellipsoids, there are a
substantial number of stars on “centrophilic” orbits that can plunge
arbitrarily close to their center (Holley-Bockelmann et al., 2002; Merritt and
Poon, 2004). Studies of these orbits have demonstrated that even moderate
triaxiality suffices to allow a MBHB to merge within a Hubble time (e.g.,
Holley-Bockelmann and Sigurdsson, 2006; Berczik et al., 2006; Preto et al.,
2011; Khan et al., 2011; Holley-Bockelmann and Khan, 2015; Vasiliev et al.,
2015; Pfister et al., 2017).

For example, recent studies show that in all cases, a MBHB forms and hardens
on timescales of at most 100 Myr, coalescing on another few-hundred-megayear
timescale, depending on the binary orbital eccentricity and the characteristic
density of stars (Khan et al., 2018). Higher central densities provide a
larger supply of stars that can efficiently extract energy from the MBHB
orbit, and shrink it to the phase where GW emission becomes dominant, leading
to the coalescence of the MBHs. Overall, the sinking of the MBHB in such
studies takes no more than $\sim 0.5\,$Gyr after the merger of the two
galaxies is complete, but this can in principle be much faster for plunging
binary orbits.

In addition to non-spherical potentials, a number of studies have considered
the evolution of MBHBs in rotating galactic nuclei (see the discussion in
Sect. 3.2.2 and Amaro-Seoane et al., 2010a; Sesana et al., 2011; Gualandris et
al., 2012; Wang et al., 2014; Holley-Bockelmann and Khan, 2015; Mirza et al.,
2017; Rasskazov and Merritt, 2017). An important result of these studies is
that in general, rotation has been found to increase the hardening rate of
MBHB even more effectively than galaxy geometry alone. In this case too, the
corotation (counterrotation) of a MBHB with stars in the galactic nucleus
results in relatively low (high) binary orbital eccentricity. For MBHBs that
end up on relatively eccentric orbits, a residual eccentricity may be
preserved all the way into the GW regime, potentially allowing to distinguish
different evolutionary scenarios from GW observations.

It is worth noting that the same types of interactions that lead to efficient
binary hardening also result in some number of stellar tidal disruptions. As a
result, the evolution of the tidal disruption rates happens in tandem with the
binary hardening rates, since both are driven by the depletion of loss-cone
orbits (Chen et al., 2011; Lezhnin and Vasiliev, 2019). Specifically, the
stellar tidal disruptions of unequal mass MBHBs can reach a rate $\sim
10^{-1}\,{\rm yr^{-1}}$ (three orders of magnitude higher than the
corresponding value for a single MBH) soon after the gravitationally bound
binary is formed (Chen et al., 2009). Depletion of loss-cone orbits leads to a
precipitous drop in the disruption rate after $\sim 10^{6}$ yr, which
asymptotes to the value typical of single MBHs ($\sim 10^{-5}$ to $\sim
10^{-4}\,{\rm yr^{-1}}$) for binaries in the GW regime. While the MBHB
hardening rate is only weakly influenced by tidal disruptions (Chen et al.,
2011), they can nevertheless give rise to bright EM flares (Coughlin et al.,
2017; Coughlin and Armitage, 2018), that may allow identification of a galaxy
that is the host to a gravitationally bound binary, as long as this scenario
can be distinguished from a disruption by a single MBH and variability of
regular AGNs.

While stellar dynamical processes seem capable of delivering many MBHBs into
the GW regime, it is prudent to consider what happens to those that experience
ineffective stellar scattering and stall as a result at sub-parsec orbital
separations. Given enough time, the remnant galaxy could undergo subsequent
mergers with other galaxies and acquire additional (i.e., third, fourth, etc.)
MBHs in the process (Mikkola and Valtonen, 1990; Heinämäki, 2001; Blaes et
al., 2002). As noted before, studies of galaxy pairs lead to the conclusion
that, on average, all massive galaxies have experienced a merger event in the
last 10 billion years (Bell et al., 2006; Lotz et al., 2008). Assuming an
average binary lifetime of 1 billion years implies that about 10 percent of
MBHBs may form a triplet. Indeed, Kelley et al. (2017) and Bonetti et al.
(2018) find that about 20–30 percent of all MBHBs could involve the presence
of a third MBH.

Some of the consequences of the intrusion of another MBH into the MBHB system
are: ($i$) the enhanced rate of loss cone refilling due to perturbation of
stellar orbits by the intruding MBH (Perets and Alexander, 2008), ($ii$)
chaotic, non-hierarchical three-body interactions, which tend to shrink the
binary semimajor axis and to increase its orbital eccentricity (Valtonen and
Mikkola, 1991), and ($iii$) if they form a hierarchical triple, the merger
time of the inner binary can be dramatically reduced due to eccentricity
oscillations induced by the Kozai-Lidov mechanism (Blaes et al., 2002). The
latter two effectively increase the binary eccentricity to values close to
unity, resulting in the GW losses and eventual MBHB coalescence on a timescale
of $\sim{\rm few}\times 10^{8}\,$yr. All of them are likely to accelerate the
MBHB coalescence and increase the ejection rate of (typically less massive)
MBHs (Hoffman and Loeb, 2007; Amaro-Seoane et al., 2010b; Ryu et al., 2018;
Bonetti et al., 2019). This bodes well for GW searches for MBHBs, as it
provides multiple means for them to reach this evolutionary stage.

####  3.3.2 Interactions with the circumbinary disk

The evolution of sub-parsec MBHBs in predominantly gaseous environments has
also been a topic of a number of theoretical studies in the past 15 years
(Armitage and Natarajan, 2005; MacFadyen and Milosavljević, 2008; Bogdanović
et al., 2008, 2011; Cuadra et al., 2009; Haiman et al., 2009; Hayasaki, 2009;
Roedig et al., 2012; Shi et al., 2012; Shi and Krolik, 2015; Noble et al.,
2012; Kocsis et al., 2012a, b; D’Orazio et al., 2013; Farris et al., 2014;
Rafikov, 2016; Schnittman et al., 2016; Tang et al., 2017; Miranda et al.,
2017; Bowen et al., 2017, 2018, 2019; Tang et al., 2018; d’Ascoli et al.,
2018; Moody et al., 2019, etc.). Most have focused on a particular physical
scenario, commonly referred to as the _binary in circumbinary disk_ , in which
a MBHB is embedded within a geometrically thin and optically thick gas disk.
This setup is of interest because a dense and relatively cold circumbinary
disk can, in principle, act as an effective conveyor belt for angular
momentum. Such disks are also thought to be radiatively efficient and similar
to disks that power AGNs – thus, they can, as an added bonus, provide a
luminous EM tracer for inspiraling MBHBs.

As noted in Sect. 3.2.2, simulations of isolated galaxy mergers indicate that
MBH pairs with mass ratios considerably lower that $q\approx 0.1$ are less
likely to form gravitationally bound binaries within a Hubble time, primarily
due to the inefficiency of dynamical friction on the lower-mass MBHs. If, at
the point of formation of the gravitationally bound binary, the two MBHs have
a mass ratio larger than a few percent, theoretical studies indicate that the
binary torques can truncate a sufficiently cold circumbinary disk and create
an inner low-density cavity with a radius corresponding to approximately twice
the binary semimajor axis (see Lin and Papaloizou, 1979; D’Orazio et al.,
2016, and references above). MBHs in this phase can accrete by capturing gas
from the inner rim of the circumbinary disk and can in this way maintain mini-
disks bound to individual holes. The left panel of Fig. 13 illustrates the
geometry of the circumbinary region for MBHBs in this evolutionary stage,
found in simulations that have investigated this scenario.

Figure 13: Illustration of the geometry of circumbinary region after the
binary with comparable mass MBHs has created a low density region in the disk.
Left: Simulations indicate that MBHs in this evolutionary phase can accrete by
capturing gas from the inner rim of the circumbinary disk and can maintain
mini-disks bound to individual holes. Adapted from Bogdanović (2015). Right:
Snapshot of gas surface density from a 2D hydrodynamic simulation of
circumbinary disk around a MBHB with mass ratio 1:4. Orbital motion is in the
counterclockwise direction. Superposed green arrows trace the fluid velocity.
Adapted from Farris et al. (2014).

Simulations also indicate that despite strong binary torques, accretion into
the central cavity continues unhindered relative to the single MBH case (e.g.,
D’Orazio et al., 2013; Farris et al., 2014; Shi and Krolik, 2015). If so, this
has two important implications: (1) that necessary conditions exist for
efficient transport of angular momentum through the circumbinary disk and
potentially, orbital evolution of the MBHBs (see the discussion below,
however) and (2) by analogy with single MBHs, binary MBHs powered by accretion
can also shine as AGNs. Another property that many simulations of MBHBs in
circumbinary disks identify robustly is that in unequal-mass binaries,
accretion occurs preferentially onto the smaller of the two objects, because
it orbits closer to the inner edge of the circumbinary disk (e.g., Artymowicz
and Lubow, 1994; Günther and Kley, 2002; Hayasaki et al., 2007; Roedig et al.,
2011; Nixon et al., 2011; Farris et al., 2014). This indicates that the lower-
mass secondary MBH may grow at a rate higher than the primary and that, given
enough time, the binary may evolve toward an equal-mass configuration.

Thus, the simplest scenario that one can consider in this case is that of a
MBHB orbit that decays under the influence of accretion torques. In this case
the inner rim of the circumbinary disk follows the MBHB inward until the
timescale for orbital decay by gravitational radiation becomes shorter than
the viscous timescale999The timescale on which angular momentum is transported
outward through the disk. of the disk (Milosavljević and Phinney, 2005;
Armitage and Natarajan, 2005). At that point, the rapid loss of orbital energy
and angular momentum through gravitational radiation can cause the binary to
detach from the circumbinary disk and to accelerate towards coalescence. If
the orbital evolution of the binary is driven by the circumbinary disk, the
rate of shrinking of the MBHB orbit can be described in terms of the viscous
inflow rate at the disk inner edge as

$\displaystyle\dot{a}_{\rm visc}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\,\alpha\left(\frac{h}{r}\right)^{2}v_{\rm Kep}$
$\displaystyle=$ $\displaystyle-1.1\times
10^{-8}c\left(\frac{\alpha}{0.1}\right)^{4/5}\left(\frac{\dot{M}}{0.1\dot{M}_{\rm
E}}\right)^{2/5}\left(\frac{10^{4}r_{g}}{a}\right)^{2/5}\left(\frac{10^{7}M_{\odot}}{M}\right)^{1/5},$

where $\alpha$ is the dimensionless viscosity parameter and $h/r$ is the
geometric aspect ratio of the disk in the Shakura and Sunyaev (1973) accretion
disk model. $v_{\rm Kep}=(GM/2a)^{1/2}$ is the circular orbital speed of the
circumbinary disk with the inner edge at $r=2a$, around a central object, here
assumed to be a binary with semimajor axis $a$ and mass $M$. Here
$\dot{M}_{\rm E}=L_{\rm E}/\eta c^{2}$ is the Eddington mass accretion rate,
$\eta$ is the radiative efficiency, $L_{\rm E}=4\pi GMm_{p}c/\sigma_{\rm T}$
is the Eddington luminosity, $\sigma_{\rm T}$ is the Thomson cross section,
and other constants have their usual meaning. We obtain the second line of Eq.
(3.3.2) by using the expression for $h/r$ for a gas pressure dominated disk
(see Eq. 2.16 in Shakura and Sunyaev, 1973).

The associated timescale for the inner edge of the disk to move inward is then

$t_{\rm visc}=-\frac{a}{\dot{a}_{\rm visc}}\approx 1.7\times 10^{6}\,{\rm
yr}\,\left(\frac{0.1}{\alpha}\right)^{4/5}\left(\frac{0.1\dot{M}_{\rm
E}}{\dot{M}}\right)^{2/5}\left(\frac{a}{10^{4}r_{g}}\right)^{7/5}\left(\frac{M}{10^{7}M_{\odot}}\right)^{6/5},$
(9)

where we assumed a binary with orbital separation $a=10^{4}\,r_{g}$ surrounded
by a disk with a steady state accretion rate of 10% of the Eddington limit. It
is worth noting that the Shakura-Sunyaev solution for an accretion disk around
a single MBH implies that for each $M$ and $\dot{M}$ there is a critical
radius beyond which the disk is gravitationally unstable to fragmentation and
star formation. This truncation radius in general arises somewhere between
$10^{2}\,r_{g}$ and $10^{6}\,r_{g}$, depending on the properties of the disk
and the mass of the central object (see Eqs. 15 and 16 in Haiman et al., 2009,
for example). This is commonly considered as the outer edge of the accretion
disk, beyond which the transport of angular momentum transitions from
accretion torques to some other mechanism, either stellar or gaseous. Given
that the structure and stability of circumbinary accretion disks is an area of
active research, the existence and location of such an outer truncation radius
are still open questions.

Equation (9) implies that if accretion torques are responsible for evolution
of the MBHB orbit, its semimajor axis should monotonically decrease with time
until the binary reaches the GW regime. Since the MBHB-circumbinary-disk
system has more angular momentum than the nominal Shakura-Sunyaev MBH-disk
system, one may also expect the timescale for orbital evolution to be longer
than estimated above by some factor that depends on the properties of the MBHB
and circumbinary disk. The expectation of systematic shrinking of binary orbit
has been confirmed by many theoretical works leading to 2017 (see the
references above), based on hydrodynamic and magnetohydrodynamic simulations
of MBHBs in circumbinary disks. Works after 2017, however, show that
significant uncertainties in the MBHB migration still remain.

The basic challenge with determining whether and how fast the MBHB orbit
shrinks arises from the fact that angular momentum loss by gravitational
torques is nearly canceled by angular momentum gain in accretion (Shi et al.,
2012), making quantitative determination of the net sensitive to details of
the calculation (Dittmann and Ryan, 2021). For example, using 2D hydrodynamic
simulations Tang et al. (2017) found that the torques that are responsible for
MBHB orbital evolution emerge from the gas very close to the MBHs, and are
contributed by the small asymmetries in the shape of the gas distribution of
the mini-disks. The implication of this result is that the timescale for
orbital evolution is sensitive to both numerical effects (numerical resolution
and dimensionality of the model) and numerical implementation of the physical
processes (such as accretion, thermodynamic properties of the disk, etc.).
Surprisingly, the simulations by Tang et al. (2017) showed that when MBHs are
modeled as rapid sinks101010These simulations do not resolve the spatial
scales close to the event horizon of the MBHs and instead model them as sinks
of larger size. and accrete gas on a timescale comparable to or shorter than
the local viscous time, the binary experiences positive torques, leading to
expansion of its orbit. A few other groups, which performed close studies of
the angular momentum transfer to the MBHB and the resulting binary orbital
evolution, also report positive orbital torques: (a) Miranda et al. (2017),
based on 2D simulations of circular binaries with different mass ratios and
eccentricities, (b) Muñoz et al. (2019) for equal-mass binaries with different
eccentricities, modeled in 2D, and (c) Moody et al. (2019) based on their 2D
and 3D simulations of circular binaries with equal-mass ratios, coplanar and
misaligned with the circumbinary disk.

More recently, Tiede et al. (2020) modeled equal-mass MBHBs on circular orbits
and found that binary orbital evolution switches from outspiralling for warm,
thicker disks (with aspect ratios $h/r\sim 0.1$), to inspiralling for cooler,
thinner disks at a critical value of $h/r\sim 0.04$. They explain that the
reason why a group of earlier studies find positive orbital torques is because
they model disks with $h/r=0.1$. Dittmann and Ryan (2022) arrive to a similar
conclusion and find an additional dependence of the orbital torque on the
assumed value of kinematic viscosity. These results underscore that the
evolution of MBHBs is sensitive to the viscous and thermodynamic properties of
the disk (reflected in its temperature and $h/r$). They also raise additional
questions regarding the degree to which the calculation of the MBHB orbital
torques may be affected by other approximations (like the description of the
accretion flow equation of state) or by effects not captured by most present-
day simulations (such as loss of angular momentum due to magnetically or
radiatively driven outflows). At the time of this writing, some insight is
provided by a study by Williamson et al. (2022), the first one to investigate
orbital evolution of binary AGNs at parsec-scale separations affected by
radiation pressure from the AGNs themselves. The study shows that radiation
pressure cyclically destroys MBH mini-disks, leading to reduced gravitational
torques and accretion, and significantly increased coalescence time scale.

The uncertainties mentioned above should also affect the evolution of the MBHB
orbital eccentricity, which is a function of both the orbital angular momentum
and energy. The works mentioned above (Tang et al., 2017; Miranda et al.,
2017; Muñoz et al., 2019; Moody et al., 2019; Tiede et al., 2020; Dittmann and
Ryan, 2021, 2022) assume a MBHB with fixed orbital elements, and thus make no
predictions about the evolution of eccentricity. Earlier works that reported
the inspiral of MBHBs found that in prograde binaries, the exchange of angular
momentum with the circumbinary disk drives a steady increase in binary
eccentricity (Armitage and Natarajan, 2005; Cuadra et al., 2009), which
saturates in the range of $0.45-0.6$ (Roedig et al., 2011; Zrake et al.,
2021). This finding seemed intuitive, as the comparable mass MBHBs were found
to produce circumbinary disks with non-axisymmetric, lopsided cavities. These
were in turn expected to induce the increase of binary orbital eccentricity
through gravitational interactions. For retrograde MBHBs, the cancelation in
angular momentum with the circumbinary disk was found to lead to growth in the
orbital eccentricity, which is slow for initially circular binaries, and
fairly rapid and leading to $e\approx 1$ for initially eccentric binaries
(Nixon et al., 2011; Roedig and Sesana, 2014; Bankert et al., 2015; Schnittman
and Krolik, 2015; Amaro-Seoane et al., 2016). The question of MBHB orbital
eccentricity is important for the final outcome of the binary evolution
(stalling or coalescence), as very eccentric MBHBs that reach small enough
pericentric distances may evolve directly into the GW regime and coalesce
rapidly (see the next section). Similarly, the orbital eccentricity of the
MBHB may also be imprinted in its EM and GW signatures at coalescence, as we
discuss in Sect. 4.

As noted in §3.2.2, the MBHs are expected to achieve an effective alignment of
their spin axes with the rotation axis of the circumnuclear disk (and thus,
with each other) by the time they form a gravitationally bound binary
(Bogdanović et al., 2007). This alignment is expected to be further enforced
after they form a binary, if their evolution continues in a circumbinary disk
with a well defined orbital angular momentum (Miller and Krolik, 2013).
Alternatively, there is no a priori reason to expect the alignment of the MBH
spins or high spins in scenarios characterized by the absence of coherent
accretion or absence of gas altogether (King et al., 2005; King and Pringle,
2006, 2007; King et al., 2008; Wang et al., 2009; Nayakshin et al., 2012;
Lodato and Gerosa, 2013).

###  3.4 The gravitational-wave phase

When the binary separation is small enough, the evolution of the MBHB is
driven primarily by gravitational radiation, rather than by interactions with
stars or gas. The subsequent evolution and coalescence of the black holes is
commonly divided into three approximately defined phases, illustrated in Fig.
14. The most protracted phase is _inspiral_ , in which the black holes are
sufficiently separated that their inward radial migration is slow compared to
their orbital motion. The second phase is _merger_ , in which the black holes
plunge together and merge. The third phase is _ringdown_ , in which the
initially asymmetric combined horizon radiates away modes and settles into a
Kerr configuration. The motivation for this classification originates from
different techniques used to describe these evolutionary stages. Specifically,
inspiral and ringdown phases can be treated semi-analytically and are modeled
using the post-Newtonian technique and perturbation method, respectively. The
merger phase, however, involves full nonlinear gravity and must usually be
treated numerically (in some cases the Effective One Body approach can also be
used to compute the full GW signal; Buonanno and Damour, 1999, 2000; Damour et
al., 2000).

Figure 14: The loss of energy and angular momentum via the emission of
gravitational radiation drives a binary to coalescence, which proceeds in
three different phases, illustrated here in terms of the GW strain amplitude
as a function of time. Figure from Baumgarte and Shapiro (2011).

Emitted gravitational radiation carries away some fraction of the MBHB energy
and angular momentum, and it is illustrative to examine how this affects the
properties of the holes in the three evolutionary stages mentioned above.
Consider a binary with component masses $m_{1}$ and $m_{2}$, total mass
$M=m_{1}+m_{2}$, and reduced mass ratio $\mu=m_{1}m_{2}/M$. Since the inspiral
phase includes the evolution of the MBHB from infinity to the innermost stable
circular orbit111111The smallest circular orbit in which a test particle can
stably orbit a massive object in general relativity. (ISCO), the energy
released is $\mu$ times the specific binding energy at the ISCO, so that
$E_{\rm inspiral}\sim\mu$. The total energy emitted in the merger and ringdown
phases may be estimated using the order of magnitude rule that the strain
amplitude from a binary with semimajor axis $a$, at a luminosity distance
$d_{L}$ from us, is $h\sim(G\mu/ac^{2})(GM/d_{L}\,c^{2})$. For the merger and
ringdown phases, the semimajor axis is comparable to the size of the event
horizon of the binary121212Here we adopt $GM/c^{2}\equiv M$ as a natural unit
of length, expressed in geometrized units when $G=c=1$, and refer to it as the
_gravitational radius_., $a\sim M$, so $h\sim\mu/d_{L}$. Since the GW
luminosity is $L\sim d_{L}^{2}\,h^{2}f^{2}$, where $f$ is the GW frequency,
then $L\sim\mu^{2}f^{2}$. If the merger and ringdown phase last a time $\tau$,
then the total energy released is $E\sim\mu^{2}f^{2}\tau$. But the
characteristic GW frequency at merger is $f\sim 1/M$ and the characteristic
time is $\tau\sim M$, resulting in $E_{\rm merger}\sim\mu^{2}/M$, or a factor
$\sim\mu/M$ times the energy released in the inspiral. As a consequence, for a
comparable-mass binary the energy emitted in the merger and ringdown could be
a few percent of the rest mass energy of the binary, and of the same order as
the energy emitted in the inspiral. For systems with extreme mass ratios
however, the inspiral is likely to be the only detectable phase.

In most of the inspiral phase, down to the separation $\sim 10-20\,M$, the
pure gravitational radiation driven evolution is well described by quadrupolar
radiation (Peters, 1964). In this approximation, the orbital semimajor axis,
$a$, and eccentricity, $e$, of the binary evolve as

$\biggl{\langle}{da\over{dt}}\biggr{\rangle}=-{64\over 5}{G^{3}\mu
M^{2}\over{c^{5}a^{3}(1-e^{2})^{7/2}}}\left(1+{73\over{24}}e^{2}+{37\over{96}}e^{4}\right)$
(10)

and

$\biggl{\langle}{de\over{dt}}\biggr{\rangle}=-{304\over{15}}e{G^{3}\mu
M^{2}\over{c^{5}a^{4}(1-e^{2})^{5/2}}}\left(1+{121\over{304}}e^{2}\right)$
(11)

where the angle brackets indicate an average over an orbit. One can show that
these rates imply that the quantity

$ae^{-12/19}(1-e^{2})\left(1+{121\over{304}}e^{2}\right)^{-870/2299}$ (12)

is constant throughout the inspiral, meaning that for moderate to small
eccentricities we have $e\propto a^{19/12}$. These equations imply that the
characteristic time to inspiral for high-eccentricity orbits is

$t_{\rm gw}\approx 10^{9}~{}{\rm yr}\left(\frac{10^{18}~{}M_{\odot}^{3}}{\mu
M^{2}}\right)\left(\frac{a}{0.001~{}{\rm pc}}\right)^{4}(1-e^{2})^{7/2}\;,$
(13)

and that the eccentricity also decreases on that timescale. The dependence of
the characteristic inspiral time on $e$, as well as circularization due to
emission of gravitational radiation can be understood by recalling that the GW
luminosity depends very strongly on the orbital frequency (and thus speed) of
two masses. For high eccentricities the orbital speed is significantly larger
at pericenter than it is at most places in the orbit. If one simplifies by
assuming that GWs are _only_ emitted at pericenter (as in the impulse
approximation), the loss of energy and angular momentum at pericenter means
that the apocenter distance decreases with time, and consequently, so does the
eccentricity.

The strong dependence of the inspiral time, $t_{\rm gw}\propto a^{4}$, means
that at large distances the binary inspirals slowly. This timescale should be
compared to the characteristic timescale for the inner edge of the
circumbinary disk to move inward, $t_{\rm visc}\propto a^{7/5}$, given by Eq.
(9). This implies that for MBHBs surrounded by circumbinary disks, the radial
diffusion of the circumbinary disk initially keeps the inner edge of the disk
at a radius $r\approx 2a$ as the binary shrinks. When $a$ becomes small enough
however, the GW inspiral time can become shorter than the disk radial
diffusion time, so that the MBHB may decouple from the disk (Pringle, 1991;
Artymowicz and Lubow, 1994, 1996; Liu et al., 2003; Milosavljević and Phinney,
2005). Comparison of the two timescales indicates that the radius of
decoupling can be $\sim 100\,M$ for comparable-mass binaries and accretion at
rates comparable to Eddington. It is reduced if one member of the MBHB is much
less massive than the other, because gravitational radiation is weaker in that
case. The implied decoupling radius is however larger at lower accretion
rates, as long as the disks are in the radiatively efficient regime (with
luminosity above $\sim 0.01\,L_{\rm E}$; see, e.g., Ichimaru, 1977; Rees et
al., 1982; Narayan and Yi, 1994; Blandford and Begelman, 1999). This is
because the disks then have smaller vertical thickness and hence smaller
inward radial speeds (see Eq. (3.3.2)).

The question of decoupling of the MBHB from the circumbinary gas disk before
coalescence is important, because coalescences that occur in gas poor
environments may not be accompanied by detectable EM counterparts. For that
reason, this question has been examined in hydrodynamic simulations of
inspiraling MBHBs in circumbinary disks. For example, Noble et al. (2012)
pointed out that magnetic fields in the circumbinary disk are stretched as the
binary shrinks, and this can lead to much larger stresses at its inner edge
than previously realized. They find that the disk can follow the binary down
to $\sim 20M$, which has encouraging implications for the detectability of EM
counterparts. Similarly, Tang et al. (2018) find that the MBHs capture gas
from the circumbinary disk and continue to accrete efficiently via their own
mini-disks, well after their inspiral outpaces the viscous evolution of the
disk.

In addition to effects on the semimajor axis and eccentricity, Schnittman
(2004) showed that the directions of the spins of the two black holes can
change systematically during GW inspiral. In particular, if the more massive
black hole is at least partially aligned with the orbital axis and the mass
ratio is close to unity, then the spins of both holes tend to align with the
orbital axis. This has been proposed by Kesden et al. (2010) and Berti et al.
(2012) as an important mechanism to further align already partially aligned
spins, which (as we discuss below) in turn reduces the kicks produced by the
MBH merger. Table 1 of Kesden et al. (2010) shows that alignment is only
effective if the mass ratio is very close to unity and that even for mass
ratios of 1/3, alignment is ineffective. Thus, this alignment mechanism might
play a role in a limited set of mergers. Interestingly, however, Gerosa et al.
(2015a) show that even initially aligned spin configurations can be brought
out of alignment as a consequence of an instability that exists for a wide
range of spin magnitudes and mass ratios, and can occur in the strong-field
regime near the merger. They show that when the spin of the higher-mass black
hole is aligned with the orbital angular momentum and the spin of the lower-
mass black hole is antialigned, binaries in these configurations are unstable
to precession to large spin misalignment. If prevalent in real MBHBs, this
instability provides a way to circumvent gas-driven astrophysical spin
alignment that may be acquired at large binary separations (see discussion in
Sect. 3.2.2 and 3.3.2), allowing significant spin precession prior to merger.

The orientation of spins, as well as the asymmetry in masses of the two black
holes, play an important role for the resulting magnitude of the “GW kick”,
received by the daughter black hole that forms in the coalescence. The merger
phase is the origin of the majority of this kick, relative to the original
center of mass of the binary. The small pitch angle of the orbit during MBHB
inspiral means that net emission of momentum in GWs nearly self-cancels over
an orbit. Analytical calculations of the kick during inspiral (Peres, 1962;
Bekenstein, 1973; Fitchett, 1983; Fitchett and Detweiler, 1984; Redmount and
Rees, 1989; Kidder, 1995; Favata et al., 2004; Blanchet et al., 2005; Damour
and Gopakumar, 2006; Kidder et al., 2007; Schnittman and Buonanno, 2007;
Schnittman et al., 2008) suggest that prior to the final orbit the net speed
relative to the original center of mass is at most tens of km s-1. However,
the momentum radiated in the last half-orbit before merger is _not_ self-
cancelled, coinciding with moments before coalescence when the emission of
gravitational radiation is strongest.

The breakthroughs that ushered in stable numerical simulations of black hole
coalescence through ringdown (Pretorius, 2005; Campanelli et al., 2006; Baker
et al., 2006) have allowed GW kicks to be computed for many configurations,
from nonspinning black holes to fully generic mergers of comparable-mass black
holes without special spin orientations. They showed that if the spins are
aligned with each other (as opposed to antialigned or at some angle) and the
orbital angular momentum, then the maximum kick is moderate: less than 200 km
s-1 (González et al., 2007). However, starting in 2007 other configurations
were found that can lead to much greater kicks. Campanelli et al. (2007a, b)
showed that for equal-mass maximally spinning black holes with oppositely-
directed spins in the orbital plane (configuration also known as the “super
kick”), the kick can reach almost 4,000 km s-1. A different, “hangup kick”
configuration, in which the spins are partway between the orbital plane and
orbital axis, was introduced in Lousto and Zlochower (2011) and Lousto et al.
(2012). The maximum velocity predicted for this configuration was shown to be
close to 5,000 km s-1.

As numerical simulations of black hole mergers are extremely computationally
expensive, several groups have developed approximate analytic formulae,
calibrated on simulations, that predict kick velocity as a function of mass
ratio and spin magnitudes and directions (Baker et al., 2007, 2008; Lousto and
Zlochower, 2008, 2009; van Meter et al., 2010a; Lousto et al., 2010, 2012). We
show the one by Lousto et al. (2012) that includes the super kick
configurations

$\displaystyle{\bf v_{\rm kick}}$ $\displaystyle=$ $\displaystyle
v_{m}\,\hat{e}_{1}+v_{\perp}(\cos{\xi}\,\hat{e}_{1}+\sin{\xi}\,\hat{e}_{2})+v_{\parallel}\,\hat{L}\,,$
(14) $\displaystyle v_{m}$ $\displaystyle=$ $\displaystyle
A_{m}\,\frac{\eta^{2}(1-q)}{1+q}(1+B_{m}\eta)\,,$ (15) $\displaystyle
v_{\perp}$ $\displaystyle=$ $\displaystyle
H\frac{\eta^{2}}{1+q}(s_{1}^{\parallel}-qs_{2}^{\parallel})\,,$ (16)
$\displaystyle v_{\parallel}$ $\displaystyle=$
$\displaystyle\frac{16\eta^{2}}{1+q}(V_{11}+V_{A}S_{z}+V_{B}S_{z}^{2}+V_{C}S_{z}^{3})|{\bf
s_{1}}^{\perp}-q\,{\bf s_{2}}^{\perp}|\cos{(\phi_{\Delta}-\phi_{1})}\,.$ (17)

The velocity component $v_{m}$ refers to the kick velocity contribution from
unequal masses and $v_{\perp}$, $v_{\parallel}$, are contributions from the
spins. The indices $\perp$ and $\parallel$ refer to the components
perpendicular and parallel to the orbital angular momentum, respectively.
Here, $\hat{e}_{1}$, $\hat{e}_{2}$ are the orthogonal unit vectors in the
orbital plane, $\xi$ is the angle that determines the projection of the spin
contribution to the kick velocity in the orbital plane ($v_{\perp}$) on
$\hat{e}_{1}$ and $\hat{e}_{2}$, and $\hat{L}$ is the unit vector in the
direction of the orbital angular momentum. $s_{1}=(cS_{1}/Gm_{1}^{2})$ is the
dimensionless spin angular momentum of the more massive, primary black hole,
and $s_{2}$ is that of the secondary, with similar definition. The parameter
$q=m_{2}/m_{1}\leq 1$ is the mass ratio, $\eta=q/(1+q)^{2}$ is the symmetric
mass ratio, and $S_{z}=2(s_{1}^{\parallel}+q^{2}s_{2}^{\parallel})/(1+q)^{2}$.
$\phi_{1}$ is the phase angle associated with the primary black hole and
$\phi_{\Delta}$ is the angle between the in-plane vector, defined in Lousto et
al. (2012) as a linear combination of the two spins, and a fiducial direction
at merger. The best fit coefficients are $A_{m}=1.2\times 10^{4}$ km s-1,
$B_{m}=-0.93$, $H=6.9\times 10^{3}$ km s-1, $V_{11}=3677.76$ km s-1,
$V_{A}=2481.21$ km s-1, $V_{B}=1792.45$ km s-1, $V_{C}=1506.52$ km s-1, and
$\xi=145^{\circ}$.

It is worth noting that improved, albeit more complex versions of the analytic
formulae are now available in the literature (e.g., Healy and Lousto, 2018).
In addition to the analytic formulae, a family of non-analytic, so called
“surrogate models”, has been developed more recently. Surrogate models are
trained against numerical simulations and were shown to be almost
indistinguishable from simulations in the portion of parameter space where
$s_{1}$ , $s_{2}\leq 0.8$ and $q\leq 1/4$ (Varma et al., 2019a, b).

Figure 15: Central escape speeds in units of ${\rm km\,s^{-1}}$ from different
types of stellar systems: core ellipticals (E$\cap$), power-law ellipticals
(E$\backslash$), dwarf ellipticals (dE), dwarf spheroidals (dSph), and
globular clusters (GC). The solid line is the mean escape speed from the DM
halos associated with the luminous matter. The dashed line is the escape speed
from the combined luminous+mean DM potentials for E galaxies. Figure from
Merritt et al. (2004).

As discussed by Merritt et al. (2004) and illustrated in Fig. 15, GW kicks
exceeding $\sim 2,000$ km s-1 suffice to eject the remnant MBH from any
galaxy. Kicks exceeding several hundred km s-1 are enough to eject remnant
MBHs from galaxies comparable to or smaller than the Milky Way. Thus, if such
kicks are produced in some galactic mergers, there may exist some number of
“empty nest” galaxies without a central MBH or with the MBH offset from the
galaxy center. We discuss possible observational signatures associated with
this phenomenon further in Sect. 4.3.

##  4 Emission signatures before, during, and after coalescence

In this section we discuss the emission signatures that arise before, during
and after coalescence. Because our goal is to review the EM counterparts to
merging MBHs that are sources of GWs, we focus on stages when the binary
evolution is driven by emission of gravitational radiation (inspiral, merger
and ringdown) rather than interactions with stars or gas. As discussed in 3,
this may happen when the MBHs are within few hundred to a thousand
gravitational radii of each other, or $\sim 10^{-3}(M/10^{8}\,M_{\odot})$ pc,
and hence their orbital period is a few years or less. At such separations,
the EM signatures are most likely to be associated with gas accretion flows
surrounding the MBHBs. As noted in previous section, stellar tidal disruptions
are not likely to play an important role for MBHBs in the GW-dominated phase
of their inspiral, since their enhanced rate would have long since asymptoted
to the value typical of single MBHs. We therefore focus attention on the MBHBs
in circumbinary disks and assess pathways to EM signatures by considering
several questions:

  * •

What is the semimajor axis at which MBHB evolution driven by stellar and gas
interactions transitions to binary evolution driven by gravitational
radiation?  • What is the orbital eccentricity of the MBHB once it transitions
to GW-driven evolution? Because emission of GWs leads to circularization, this
transition determines the eccentricity of the MBHB at smaller separations and
thus affects possibilities for time-variability of its EM signatures. • What
fraction of the matter in the circumbinary disk eventually accretes EM
counterparts of merging MBHBs can originate from $\bullet$ The circumbinary
disk.
$\bullet$ Mini-disks and accretion streams. $\bullet$ Magnetized jets.
$\bullet$ If all else fails: emission from the host galaxy.
onto either of the MBHs? Even if this fraction is small, the high efficiency
of energy production in the mini-disks around each BH could mean that they
dominate the overall luminosity.

  * •

Down to what spatial scale can the circumbinary disk follow the GW-driven
shrinking of the MBHB? This could have a large impact on the simultaneity of
EM and GW emission during the merger, as well as the EM luminosity of any
merger-related event.

  * •

What is the role of periodic or quasi-periodic modulation of EM emission, as
opposed to secular variability? Nearly periodic modulation, if observed over
enough cycles, will be more easily distinguished from background noise than a
gradual increase in intensity. Another important question is whether the
fractional variation in intensity for periodic MBHB sources will be large
enough to be detected, especially when combined with emission from the host
galaxy.

We will keep these questions in mind as we lay out analytic estimates and
discuss possibilities for production of characteristic EM signatures from
merging MBHBs presented in the literature.

###  4.1 Electromagnetic radiation before the merger

####  4.1.1 Emission from the circumbinary disk

The possible locations of EM emission before the MBHB merger are the
circumbinary disk, the mini-disks around the individual black holes, and the
streams feeding the mini-disks (see Fig. 13). Let us focus first on the
circumbinary disk, where we will present some simple analytic estimates to
guide the intuition and compare them with the results from simulations, paying
attention to situations where they differ. The analytic estimates are obtained
under the assumption that the circumbinary disk luminosity is produced by
gravitational energy and angular momentum transport rather than by
perturbations induced by the binary torques. We will also assume that the mass
accretion rate through the disk is constant, and that it does not change as
the MBHB inspirals and the inner rim of the circumbinary disk shrinks. As the
binary inspirals, the disk then becomes more luminous and its spectrum peaks
at higher photon energies.

This picture is similar to the scenario envisioned by many researchers
(Pringle, 1991; Artymowicz and Lubow, 1994, 1996; Liu et al., 2003;
Milosavljević and Phinney, 2005). They noted that the GW-driven inspiral time
decreases more rapidly with semimajor axis ($t_{\rm gw}\propto a^{4}$, Eq.
(13)) than does the viscous inflow time ($t_{\rm visc}\propto a^{7/5}$, Eq.
(9)). Therefore, when the binary semimajor axis decreases sufficiently, the
binary may decouple from the disk. As noted in Sect. 3.4, comparison of the
two timescales indicates that the decoupling may happen when $a\sim 100\,M$
for comparable-mass binaries and accretion rates comparable to Eddington. If
so, this implies that the gradual refilling of the hole after MBHB merger
could provide a distinct signature of the merger (as suggested by
Milosavljević and Phinney, 2005). It also implies that the circumbinary disk
will not be very luminous. The reason is that if the inner edge of the
circumbinary disk is at $r_{\rm in}\approx 2a\sim 200\,M$ when the binary
decouples, then the energy is released with only $M/r_{\rm in}\sim 0.5$%
efficiency. In comparison, the radiative efficiency of any matter accreting
all the way to ISCO of either MBH is about 10%, so the luminosity of the
circumbinary disk is at least a factor of twenty less. Based on these simple
arguments, the emission luminosity from the MBH mini-disks is expected to
dominate that from the circumbinary disk in a steady-state accretion flow some
time before merger.

These considerations also imply that when the MBHBs merge, the circumbinary
disk may be left behind at some appreciable distance. If so, there may be
little matter around the MBHB at the moment of coalescence and hence, few
opportunities for direct interaction of that matter with the dynamically
changing spacetime near the MBHs that could produce a luminous EM counterpart.
This early concern was mitigated by simulations (Noble et al., 2012; Farris et
al., 2015a; Bowen et al., 2017; Tang et al., 2018), which indicate that the
inner rim of the circumbinary disk is able to follow the binary to radii of
$\sim 10M$. This is also true for the spinning MBHBs, as the spins do not
strongly affect the bulk properties of the circumbinary disk (Lopez Armengol
et al., 2021). For the circumbinary disk with an inner rim radius of $\sim
20M$, the energy release could then be $\sim 5$% and substantial enough to
produce an EM luminous signature. Similarly, the inspiral time of the gas is
then much shorter than previously assumed, thus allowing the MBHs to accrete
all the way to coalescence. This is illustrated in Fig. 16, where the two MBHs
are fed by the narrow gas streams till coalescence.

Figure 16: Surface density of the gas around an equal mass $10^{6}\,M_{\odot}$
MBHB evolving from the initial separation of $60\,M$ (left panel) to merger
(right), on a timescale of about a week. Despite the presence of an asymmetric
low density cavity, the two MBHs continue to accrete gas carried by the narrow
gas streams all the way through coalescence. The blue and green (orange and
red) colors indicate low (high) surface densities. Figure adapted from Tang et
al. (2018).

The prime difficulty with gravitationally-powered circumbinary disk emission
as a signature of binary MBHB coalescence is that the gradual increase in
luminosity and change in spectrum is smooth and thus, might not be easily
distinguishable from other gradual changes that exist naturally in disks
around single AGN. As a reminder, we expect at most a few MBHB mergers per
year out to a redshift $z\sim 1$, compared to the millions of AGNs that exist
to that redshift. Thus without a contemporaneous GW signal and decent sky
localization to guide a search, one needs to identify a signature that exists