# Uncovering Hidden Massive Black Hole Companions with Tidal Disruption Events

Brenna Mockler The Observatories of the Carnegie Institution for Science,
Pasadena, CA 91101, USA  Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095, USA  Denyz Melchor Smadar Naoz Department
of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and
Astronomy, UCLA, Los Angeles, CA 90095, USA  Enrico Ramirez-Ruiz Department of
Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064,
USA

###### Abstract

Dynamical perturbations from supermassive black hole (SMBH) binaries can
increase the rates of tidal disruption events (TDEs). However, most previous
work focuses on TDEs from the heavier black hole in the SMBH binary (SMBHB)
system. In this work, we focus on the lighter black holes in SMBHB systems and
show that they can experience a similarly dramatic increase in their TDE rate
due to perturbations from a more massive companion. While the increase in TDEs
around the more massive black hole is mostly due to chaotic orbital
perturbations, we find that, around the smaller black hole, the eccentric
Kozai-Lidov (EKL) mechanism is dominant and capable of producing a comparably
large number of TDEs. In this scenario, the mass derived from the light curve
and spectra of TDEs disrupted by the lighter SMBH companion is expected to be
significantly smaller than the SMBH mass estimated from galaxy scaling
relations, which are dominated by the more massive companion. This apparent
inconsistency can help find SMBHB candidates that are not currently accreting
as active galactic nuclei (AGN) and that are at separations too small to be
resolved as two distinct sources. In the most extreme cases, these TDEs
provide us with the exciting opportunity to study SMBHBs in galaxies where the
primary SMBH is too massive to disrupt sun-like stars.

stars: black holes — stars: tidal disruption events — galaxies: nuclei –
galaxies: active — galaxies: supermassive black holes

††software: astropy (Astropy Collaboration et al., 2013)

##  1 Introduction

Galaxy mergers are a natural consequence of hierarchical structure formation
(Barnes, 2002; Springel & Hernquist, 2005; Robertson et al., 2006; Hopkins et
al., 2008). Since nearly every galaxy hosts a supermassive black hole (SMBH),
binary SMBHs (SMBHBs) should be common (e.g., Sillanpaa et al., 1988;
Rodriguez et al., 2006; Komossa et al., 2008; Bogdanović et al., 2009; Boroson
& Lauer, 2009; Dotti et al., 2009; Batcheldor et al., 2010; Deane et al.,
2014; Runnoe et al., 2017; Pesce et al., 2018). Observations of merging
galaxies and of dual active galactic nuclei (AGN) give credence to this idea
(e.g., Komossa et al., 2004; Bianchi et al., 2008; Comerford et al., 2009,
2018; Greene et al., 2010; Liu et al., 2010; Smith et al., 2010; Stemo et al.,
2021; Di Matteo et al., 2005; Li et al., 2020), and it has even been suggested
that our Galaxy might host a small companion to our own SMBH (e.g., Hansen &
Milosavljević, 2003; Maillard et al., 2004; Gürkan & Rasio, 2005; Gualandris &
Merritt, 2009; Chen & Liu, 2013; Fragione et al., 2020; GRAVITY Collaboration
et al., 2020; Generozov & Madigan, 2020; Naoz et al., 2020; Zheng et al.,
2020)

After two galaxies merge, the central supermassive black holes (SMBHs) will
eventually sink to the bottom of the merged galaxy’s potential, where they
will form a binary system in the nucleus of the galaxy. While the two SMBHs
coexist, they will induce orbital perturbations in each others’ nuclear star
clusters. These perturbations can send stars on high eccentricity, nearly
radial orbits towards one of the black holes and lead to increased rates of
tidal disruptions events (e.g. Ivanov et al., 2005; Chen et al., 2011; Wegg &
Nate Bode, 2011; Li et al., 2015, 2019; Fragione & Leigh, 2018).

Tidal disruption events (TDEs) occur when a star passes close enough to a
black hole that the tidal force from the supermassive black hole (SMBH)
overwhelms the self-gravity of the star (e.g. Rees, 1988; Evans & Kochanek,
1989). The star is ripped apart and the stellar debris is accreted by the
SMBH, producing a flare similar in brightness to supernovae that encodes
properties of the star and SMBH in its light curve and spectra (e.g. Lodato et
al., 2009; Ramirez-Ruiz & Rosswog, 2009; Guillochon & Ramirez-Ruiz, 2013; Yang
et al., 2017; Law-Smith et al., 2019; Mockler et al., 2019, 2022). Theoretical
tidal disruption event rates have most often been calculated using the
dynamics of two-body relaxation (e.g., Magorrian & Tremaine, 1999; Wang &
Merritt, 2004; MacLeod et al., 2012; Stone & Metzger, 2016). In this picture,
a star will experience random kicks to its angular momentum through
interactions with other stars, and if it is very unlucky, one of these kicks
will send it radially toward the SMBH, where it will be tidally disrupted.
However, in a SMBH binary, 3-body interactions can produce TDEs much faster
than two-body relaxation processes (e.g. Ivanov et al., 2005; Chen et al.,
2011; Wegg & Nate Bode, 2011).

In general, SMBH binary systems are very difficult to uncover observationally.
If one or both of the black holes is accreting as an AGN, it is sometimes
possible to measure the shift in velocity in the emission lines due to the
binary orbit (Gaskell, 1983; Boroson & Lauer, 2009; Decarli et al., 2010;
Tsalmantza et al., 2011; Eracleous et al., 2012; Runnoe et al., 2015; Wang et
al., 2017). There are also predictions for how the overall shape of the broad
emission lines might differ from a solitary AGN (Decarli et al., 2010; Nguyen
& Bogdanović, 2016). Unfortunately, measuring the velocity shift in emission
lines is only possible over a limited range of SMBHB separations. The SMBHs
must be close enough such that the velocity shifts in the broad line regions
(BLR) are measureable, but not so close that the BLRs are truncated (see
overview in De Rosa et al., 2019).

If we are lucky enough to catch both SMBHs actively accreting, this is likely
the surest sign of a SMBH binary (SMBHB) system, however, resolving these
systems requires the SMBHs to still be at large separations (before the system
has formed a bound binary) and are also mostly limited to binaries with mass
ratios near $q=1$ (e.g., Rodriguez et al., 2006; Comerford et al., 2018; Dey
et al., 2019; Goulding et al., 2019; Bhattacharya et al., 2020; Foord et al.,
2020; Kim et al., 2020; Reines et al., 2020; Monageng et al., 2021; Severgnini
et al., 2021; Shen et al., 2021). Additionally, even though the merger process
can funnel gas to the center of galaxies and trigger AGN activity, most SMBHBs
in the local universe do not currently host AGN (Casey-Clyde et al., 2022).

When stars are disrupted by the smaller SMBH in the binary, the black hole
mass measured from the TDE light curve or spectra will correspond to the
smaller SMBH (e.g. Mockler et al., 2019; Ramsden et al., 2022; Wen et al.,
2020; Ryu et al., 2020). Given that the potential of the galactic nucleus will
be dominated by the larger SMBH, typical scaling relations are likely to
produce estimates close to the mass of the larger black hole in the binary
(Izquierdo-Villalba et al., 2022). Therefore, if TDE light curve measurements
suggest a much smaller black hole mass than scaling relations such as
$M-\sigma$ or the bulge-mass relation, this is an indication that the galaxy
might host a SMBH binary111Note that if the SMBHB is on a very close orbit
such that the orbital timescale is comparable to the fallback timescale of
tidal debris, other observables such as truncations in the accretion rate can
also become important (e.g. Mardling & Aarseth, 2001; Ricarte et al., 2016).
Similarly, if a TDE is observed in a galaxy that appears to be too massive to
host an observable flare (i.e. a system with an expected black hole mass
$\gtrsim 10^{8}{\rm M}_{\odot}$), this also suggests the existence of a
second, smaller SMBH in the galactic nucleus (e.g. Coughlin & Armitage, 2018).

In this paper, we study how hierarchical 3-body interactions between the SMBH
binary and the stars orbiting the smaller black hole can excite the
eccentricities of the stars and increase the rates of tidal disruption events
around the smaller black hole. In particular, we focus on the effect of the
secular orbital perturbations known as the eccentric Kozai-Lidov (EKL)
mechanism (Kozai, 1962; Lidov, 1962; Naoz, 2016). It was shown that SMBH
binaries and the EKL mechanism can lead to increased rates of TDEs. However,
most previous work on this subject has either focused on the stars surrounding
the larger SMBH (e.g. Ivanov et al., 2005; Chen et al., 2011; Wegg & Nate
Bode, 2011; Fragione et al., 2021), on IMBH systems (Fragione & Leigh, 2018),
or on extreme mass ratio inspiral (EMRI) transients instead of TDEs (Naoz et
al., 2022; Mazzolari et al., 2022). Past work has also shown that other
dynamical mechanisms suppress EKL oscillations for stars around the larger
black hole, and the majority of TDEs that are produced around the larger SMBH
come from chaotic orbital crossings (e.g., Chen et al., 2011). This is
generally not the case for the stars around the smaller SMBH, largely because
the EKL perturbations due to the more massive SMBH produce stronger
excitations in the stars’ orbital eccentricities. Here we expand on work by Li
et al. (2015) and explore the effect of the EKL mechanism on the rates of TDEs
around the smaller SMBH in SMBH binaries.

We propose that every SMBH binary system will undergo a burst (short, enhanced
period) of TDEs from the smaller black hole once the binary separation is of
order the sphere of influence (so long as one of the SMBHs is below the Hills
mass for the star, $\approx 10^{8}{\rm M}_{\odot}$). At these separations, the
binary will not be resolved spatially, and the TDE will appear to be at the
center of the galaxy. TDEs from the smaller black hole provide a unique avenue
to discover SMBHBs at much closer separations than the kiloparsec scales at
which dual AGN are observed. This is especially significant when the more
massive SMBH is unable to disrupt sun-like stars outside its event horizon (as
also discussed in Fragione & Leigh, 2018) .

##  2 Methodology

Figure 1:  Cartoon of simulation setup (not to scale). We study stars orbiting
the smaller of the two SMBHs, and focus on the region within the hierarchical
radius (defined in Equation 2), as orbits in this region are less likely to be
chaotically scattered by the second black hole.

###  2.1 Physical processes

We consider a SMBH binary $m_{1}<m_{2}$, separated by a semi-major axis
$a_{\rm bin}$, and with an eccentricity $e_{\rm bin}$. Throughout the paper,
the subscripts ‘1’, ‘2’, and ‘*’, denote the inner black hole, the outer black
hole, and the star, respectively. The subscript ‘bin’ refers to the black hole
binary. Surrounding each SMBH is a nuclear star cluster, and for the purposes
of this paper we will focus on the stars surrounding $m_{1}$. Because
$m_{1}<m_{2}$, we can think of the nuclear star cluster of $m_{1}$ as embedded
in the larger cluster around $m_{2}$. We expect the smaller SMBH $m_{1}$ to
retain stars within its Roche lobe (e.g. O’Leary & Loeb, 2012; Rantala et al.,
2018; Greene et al., 2021). Each star in the cluster around $m_{1}$ can be
modeled as part of a hierarchical three-body system made up of an inner binary
pair – $m_{*}$ and $m_{1}$ – and an outer binary pair – $m_{1}$ and $m_{2}$
(see illustration in Figure 1). We assume the stars orbiting $m_{1}$ follow
approximately Keplerian orbits, but over time their orbits can be secularly
perturbed by $m_{2}$ through the eccentric Kozai-Lidov mechanism (e.g. Naoz,
2016).

Figure 2: Example run showing EKL cycles as a function of time. The dashed
horizontal line denotes our critical radius for partial disruption
($\beta=0.5$), and the dotted vertical lines denote $t_{\rm quad}$ given the
initial orbital conditions. As you can see, $t_{\rm quad}$ is of order the
period of one eccentricity oscillation, as expected. This simulation is part
of the set of runs with $m_{1}=10^{6}{\rm M_{\odot}},\;m_{2}=10^{7}{\rm
M_{\odot}},\;\alpha=2$.

We restrict the initial parameter space to stable configurations, adopting a
hierarchical limit that ensures the orbit of the star stays well within the
orbit of the binary.

$\epsilon=\frac{a_{*}}{a_{\rm bin}}\frac{e_{\rm bin}}{1-e_{\rm bin}^{2}}<0.1\
,$ (1)

where $\epsilon$ is the pre-factor in front of the octupole expansion term of
the Hamiltonian (e.g., Naoz et al., 2013)222This limiting condition was found
by Naoz et al. (2013) to give results that are similar to the more complicated
Mardling & Aarseth (2001) stability criterion. This defines a ‘hierarchical
radius’ limit:

$a_{\rm hier}=0.1\times a_{\rm bin}\frac{1-e_{\rm bin}^{2}}{e_{\rm bin}}\ .$
(2)

The EKL mechanism describes the coherent perturbations from $m_{2}$ on the
stars orbiting $m_{1}$. These perturbations exchange angular momentum between
the outer SMBH binary ($m_{1}$ and $m_{2}$) and the inner binaries made up of
stars orbiting $m_{1}$. This produces excitations in the eccentricity and
inclination of the stars orbiting $m_{1}$ over a wide range of the parameter
space (e.g., Naoz et al., 2011; Katz et al., 2011; Lithwick & Naoz, 2011; Naoz
et al., 2013; Li et al., 2014; Will, 2017). Importantly, the eccentricity
excitations can lead to the production of a large number of tidal disruption
events (e.g., Li et al., 2015; Fragione & Leigh, 2018). To model these systems
and determine the rate of tidal disruption events, we employ the secular
approximation, where the inner and outer orbits are treated as phase-averaged
wires that torque each other. We solve the secular, hierarchical three-body
equations following Naoz et al. (2013), including up to the octupole order
terms ($\mathcal{O}(a_{*}/a_{\rm bin})^{3}$). The timescale associated with
the quadrupole level of approximation is (e.g., Antognini, 2015):

$t_{\rm quad}=\frac{16}{15}\frac{a_{\rm bin}^{3}(1-e_{\rm
bin}^{2})^{3/2}\sqrt{(m_{1}+m_{*})}}{\sqrt{G}a_{*}^{3/2}m_{2}}\ ,$ (3)

where $G$ is the gravitational constant. This timescale is associated with the
period of oscillations of the eccentricity and inclination (see Figure 2).

Figure 3: Analytical exploration of the parameter space showing the number of
stars susceptible to the EKL mechanism. The white stars represent the
parameter combinations chosen for numerical simulation runs. This was
calculated by determining the range of semi-major axes for each parameter
combination where the EKL timescale is the shortest dynamical timescale. All
of our systems have $e_{\rm bin}=0.5$, $a_{\rm bin}=0.5r_{h,\;\rm
disrupting}$. While we compared against a variety of dynamical timescales
(including those for two-body relaxation, vector resonant relaxation, general
relativistic precession, and binary hardening through gravitational wave
emission, see Appendix A), the two competing timescales were $t_{\rm GR,\;1}$
and $t_{\rm NT}$. These are the timescales for general relativistic (GR)
precession due to the disrupting black hole and Newtonian precession due to
the stellar potential, respectively.

Li et al. (2015) considered the relevant dynamical timescales at the center of
galaxies and compared them to the EKL timescale. They found that the dynamical
processes competing with the EKL mechanism for the largest effect on the
stellar orbital parameters are general relativistic precession (due to the
potential of the disrupting black hole), and Newtonian precession (due to the
potential of the stellar cusp). The timescale for general relativistic
precession is,

$t_{\rm GR1}=\frac{2\pi
a_{*}^{5/2}c^{2}(1-e_{*}^{2})}{3G^{3/2}(m_{1}+m_{*})^{3/2}}\ ,$ (4)

where $c$ is the speed of light. The Newtonian (NT) precession timescale is
given by:

$t_{\rm NT}=2\pi\left(\frac{\sqrt{Gm_{1}/a_{*}^{3}}}{\pi
m_{1}e_{*}}\int_{0}^{\pi}{\rm d}\psi\>M_{\ast}(r)\>{\rm cos}\psi\right)^{-1}\
,\\\ $ (5)

where $\psi$ is the true anomaly of the inner orbit,
$r(\psi)=a_{*}(1-e_{*}^{2})/(1+e_{*}{\rm cos}\psi)$ from Kepler’s equation.
The function $M_{\ast}(r)$ is the integrated stellar mass within $r(\psi)$ –
we will discuss the density profiles used to calculate this in Section 2.2. If
precession due to GR or the stellar potential is strong enough, it can prevent
EKL perturbations from building up coherently and exciting the eccentricity
and inclination of the star (e.g., Li et al., 2015; Chen et al., 2011).

Table 1: Simulation parameters and results. See Section 2.3 for a description of how upper and lower limits were calculated for the rightmost three columns. run # | $m_{1}$ | $m_{3}$ | $\alpha$ |  $a_{1}$ max | $N_{*}$ | TDE fraction | $N_{TDE}$ | $N_{plunge}$ |  |
---|---|---|---|---|---|---|---|---|---|---
| $\rm(log_{10}{\rm M}_{\odot})$ | $\rm(log_{10}{\rm M}_{\odot})$ |  | (au) |  |  |  |  |  |
0 | 5.0 | 6.0 | 1.50 | 6002 | 5134 | 0.11 - 0.44 | 576 - 2277 | 0 - 5 |  |
1 | 6.0 | 7.0 | 1.00 | 18980 | 14081 | 0.07 - 0.45 | 986 - 6379 | 0 - 70 |  |
2 | 6.0 | 7.0 | 1.50 | 18980 | 51307 | 0.09 - 0.41 | 4618 - 20984 | 51 - 257 |  |
3 | 6.0 | 7.0 | 2.00 | 18980 | 182705 | 0.09 - 0.33 | 15713 - 59562 | 0 - 731 |  |
4 | 6.0 | 8.0 | 1.50 | 18980 | 51366 | 0.01 - 0.43 | 462 - 22242 | 0 - 308 |  |
5 | 7.0 | 8.0 | 1.00 | 60021 | 140784 | 0.05 - 0.38 | 7180 - 54061 | 141 - 3801 |  |
6 | 7.0 | 8.0 | 1.50 | 60021 | 512509 | 0.05 - 0.32 | 23575 - 165028 | 1025 - 10250 |  |
7 | 7.0 | 8.0 | 2.00 | 60021 | 1860848 | 0.05 - 0.27 | 100486 - 504290 | 0 - 20469 |  |

Motivated by a similar analysis in Li et al. (2015), in Figure 3, we look at
how changing different parameters of the system changes the number of stars
likely to experience large eccentricity oscillations from the EKL mechanism.
For each point in Figure 3, we calculate the number of stars within our
hierarchical limit (Equation 2) whose shortest dynamical timescale is $t_{\rm
quad}$, and plot density contours displaying this information. The dynamical
timescales we compare against are described in the figure caption (and
included in Appendix A). Like Li et al. (2015), we find that $t_{\rm GR1}$ and
$t_{\rm NT}$ are the two shortest timescales apart from $t_{\rm quad}$,
implying that general relativistic precession due to the disrupting black hole
and Newtonian precession due to the stellar cusp are the two dynamical
processes most likely to interfere with EKL oscillations. This parameter space
exploration also further motivates our focus on the smaller black hole in the
binary. The number of stars susceptible to EKL oscillations drops off sharply
as the disrupting black hole grows to a comparable size to the outer black
hole (previous work focusing on the larger black hole found that the majority
of disruptions were caused by chaotic orbital crossings, not the EKL
mechanism, e.g., Wegg & Nate Bode, 2011; Chen et al., 2011).

Orbital precession due to GR and the stellar cusp can be important at small
radii, however the strength of EKL perturbations increase with a star’s radius
from $m_{1}$ as stars move closer to the perturbing black hole ($m_{2}$).
Therefore it is more difficult for the EKL mechanism to excite the
eccentricities of stars at smaller semi-major axes, but at larger radii, where
there are many more stars, the EKL mechanism often dominates and can produce
many TDEs. For this work we chose to isolate the effects of the EKL mechanism
and used the parameter space exploration described by Figure 3 to choose
systems to model where the EKL timescale $t_{\rm KL}$ is the shortest
dynamical timescale for the majority of stars surrounding the smaller black
hole.

We include GR precession in our modeling, and while it is possible that it
suppresses the number of disruptions at small radii close to $m_{1}$ (see
Figure 6), we find that the EKL mechanism acts much faster than GR precession
for the majority of the stars around $m_{1}$ in the systems we are modeling,
as expected (see for similar results Naoz & Silk, 2014; Li et al., 2015; Naoz
et al., 2022). Unlike Li et al. (2015), we do not include the NT precession
from the stellar cusp in our models, and instead chose our setups such that
the timescale for NT precession is longer than the timescale for EKL
excitations for the majority of the stars. Thus, our results are consistent
with Li et al. (2015) (see Section 3). Finally, we note that we do not include
any star-star collisions or compact object collisions. We expect that these
processes will generally become important on longer timescales (e.g., Rose et
al., 2022, 2023).

Another possible dynamical mechanism that can affect EKL perturbations is two-
body relaxation. Naoz et al. (2022) showed that even though the canonical
timescale for two-body relaxation to significantly change stellar orbits is
many orders of magnitude larger than the Kozai-Lidov timescale in the relevant
parameter space (consistent with our calculations for Figure 3), individual
(incoherent) kicks from two-body scattering can affect the build up of Kozai-
Lidov perturbations. In this paper we focus on the effects of the EKL
mechanism alone, however, in a companion paper we explore the effect of two-
body relaxation on the EKL mechanism (see Melchor et al., 2023).

###  2.2 Monte Carlo dynamical simulations

Using Figure 3 as our guide, we chose a variety of parameter combinations for
our simulations that we expected to produce observable tidal disruption
events. The parameters we varied were: the mass of the disrupting and
perturbing black holes ($m_{1}$ and $m_{2}$), their mass ratio
($q=m_{2}/m_{1}$), and the power law constant of the density profile for the
stellar cusp ($\alpha$). The mass of $m_{1}$ was set to a range of reasonable
values for tidal disruption events: $10^{5},10^{6},$ and $10^{7}$ solar
masses. Black holes above $\approx 10^{8}{\rm M}_{\odot}$ will swallow most
stars whole, and disrupt them within their event horizons. Black holes below
$10^{5}{\rm M}_{\odot}$ can and do disrupt stars (Law-Smith et al., 2017),
however these disruptions may be more difficult to observe (e.g. due to slower
circularization and lower Eddington limits, Guillochon & Ramirez-Ruiz, 2015),
and almost all confirmed TDEs have occurred around black holes above
$10^{5}{\rm M}_{\odot}$ (e.g., Ramsden et al. 2022; Hammerstein et al. 2022,
with the exception of Angus et al. 2022). Therefore, we focus on this mass
range for the purposes of this paper (see Fragione & Leigh 2018 for
simulations of similar systems focusing instead on IMBH-SMBH binaries). We use
mass ratios $q=m_{2}/m_{1}$ of 10 and 100, and this determines the mass of
$m_{2}$.

The stellar cusp is modeled using the density distribution described in
O’Leary et al. (2009):

$\rho_{\ast}(r)=\frac{3-\alpha}{2\pi}\frac{m_{1}}{r^{3}}\left(\frac{GM_{0}(m_{1}/M_{0})^{1-2/k}}{\sigma_{0}^{2}r}\right)^{-3+\alpha}\
,$ (6)

and our sphere of influence is defined where the mass enclosed is equal to
$2\times m_{1}$,

$r_{h}=\frac{GM_{0}(m_{1}/M_{0})^{1-2/k}}{\sigma_{0}^{2}}\ .$ (7)

The density profile is written as a power law of semi-major axis, with the
integrated stellar mass normalization fixed by the $M-\sigma$ relation. The
constants $M_{0}=1.3\times 10^{8}{\rm M}_{\odot}$, $\sigma_{0}=200~{}\rm
km~{}s^{-1}$, and $k=4$ are the best-fit values for the $M-\sigma$ relation in
Tremaine et al. (2002) ($M_{h}/{\rm M}_{\odot}=M_{0}(\sigma/\sigma_{0})^{k}$).
We vary the power law exponent ‘$\alpha$’ between 1-2 to explore how the
steepness of the stellar cusp affects TDE rates. While our own galactic center
appears to have a flatter cusp ($\alpha\sim 1$, e.g., Schödel et al., 2018),
TDE host galaxies tend to be very centrally concentrated and to have stellar
light profiles indicative of steep stellar cusps down to the size scales that
are observable (e.g., Law-Smith et al., 2017; French et al., 2020; Dodd et
al., 2021).

Figure 4: The final distributions of the stars’ semi-major axes ($a_{1}$) and
maximum eccentricities ($e_{\rm max}$, plotted as $1-e_{\rm max}$) for 3 runs
with varying inner black hole mass ($m_{1}=10^{5},10^{6},10^{7}{\rm
M}_{\odot}$). The critical radius for disruption used here is $R_{p}\leq
2R_{T}$ (allowing for both partial and full disruptions), and denoted by
dashed lines in the plot. Disruptions are plotted with a star-shaped marker,
and plunges (disruptions that occur with $R_{T}<R_{S}$) are shaded in gray.
The maximum semi-major axis included for stars in a given run is defined by
the hierarchical condition (Equation 1).

We set the eccentricity of the SMBH binary to an intermediate eccentricity of
$e_{\rm bin}=0.5$ and its semi-major axis to half of the sphere of influence
of the disrupting black hole ($a_{\rm bin}=0.5r_{h,\;1}$). This puts the SMBHs
at separations of $\sim 1$ pc ($a_{\rm bin}=0.2-1.9$ pc for
$m_{1}=10^{5}-10^{7}{\rm M}_{\odot}$) – a regime where the EKL mechanism is
highly effective, however the SMBH binary’s hardening timescale is still long
enough that we do not have to worry about the binary’s orbit changing
significantly over the course of our simulations. At these radii, the SMBH
binary’s hardening timescale is mainly dependent on the effectiveness of loss-
cone diffusion, which is likely $\gtrsim 100$ million years (e.g. Kelley et
al., 2017), much longer than the EKL timescales for the majority of the stars
in our simulation setups.

For each set of parameter combinations, we ran a 1000 numerical simulations of
the 3-body interactions with different initial conditions for the star. The
priors for the star’s initial conditions were determined as follows: The semi-
major axis of the star was drawn randomly from the density distribution, with
a maximum radius determined by the hierarchical condition of the system (see
Equation 1). The star’s eccentricity was drawn from an isothermal
distribution, and the mutual inclination was drawn from an isotropic
distribution (uniform in cos($i$)). The argument of periapsis and the
longitude of the ascending nodes were drawn from a uniform distribution
between $0-2\pi$. See Table 1 for an overview of the simulation parameters.

The simulations were run until either the star was disrupted or the system was
evolved for $10^{9}$ years (by which point the TDE rate from the EKL mechanism
had long since dropped off, see Figure 5). Our criteria for a tidal disruption
to occur was that the star passed within twice the tidal radius
($R_{p}=2R_{T}$) of the black hole

$R_{T}=R_{*}\Big{(}\frac{M_{h}}{M_{*}}\Big{)}^{1/3}\ .$ (8)

This meant that our minimum impact parameter for disruption was
$\beta=R_{T}/R_{p}=0.5$. Setting our TDE radius condition to $2R_{T}$ instead
of $R_{T}$ ensured that we captured partial disruptions as well as full
disruptions333Below $\beta\approx 0.5$ very little mass will return to the
black hole and it is unlikely the TDE will be observable (e.g., Guillochon,
2013).

Each set of a 1000 simulations provides a random Monte Carlo sampling of the
stars around the smaller black hole in each binary system, given the priors
described above. We plot the final eccentricity distributions of a subset of
our runs in Figure 4, where the runs that end in disruption are clearly
visible as a buildup of points at high eccentricities (low values of $1-e$).

Our simulations give us the fraction of stars that result in TDEs for each
parameter combination, as well as the time of disruption for each random
sample. We also find that a number of stars plunge directly into the black
hole, with pericenter radii less than the Schwarszchild radii of the black
holes (shaded in gray in Figure 4). We note that if the SMBH is spinning, the
star may zoom around the ergosphere instead of plunging directly into the
black hole (e.g., Glampedakis & Kennefick, 2002; Healy et al., 2009;
Schnittman, 2015). This can continue for a few tens to thousands of times the
pericenter passage (see figure 3 in Naoz et al., 2019).

Finally, we can scale our results for the fraction of stars disrupted by the
expected total number of stars around $m_{1}$ within the hierarchical radius
(using Equation 6) to compare to real galaxies. The simulation parameters and
the fraction and number of stars that become TDEs are quoted in Table 1 for
each of our runs. We can also calculate a TDE rate as a function of time by
binning our simulation results as a function of the time of disruption. We
plot our scaled, time-dependent rates for a subset of the runs in Figure 5
(using Gaussian processes to smooth the distributions).

###  2.3 Limits and stability

For each set of simulations, we calculate an upper and lower limit for the TDE
fraction. The lower limit assumes that any star whose orbit evolves such that
its apocenter moves outside the Hill radius of the smaller black hole is
scattered away from the smaller black hole and is not disrupted. We define the
Hill radius at the black hole binary’s pericenter to ensure our lower limit is
sufficiently conservative,

$r_{\rm Hill}=a_{\rm bin}(1-e_{\rm
bin})\left(\frac{m_{1}}{3m_{2}}\right)^{1/3}\ ,$ (9)

The upper limit assumes that all stars (within our hierarchical radius limit,
see Equation 2) on disrupting orbits become TDEs, even those whose apocenters
move outside the Hill radius before disruption. While the Hill radius is a
good approximation of the stability boundary on long timescales (e.g., Grishin
et al., 2017; Bhaskar et al., 2021; Tory et al., 2022), stars at larger radii
have shorter EKL timescales, and can easily be disrupted before their orbits
become unstable. Stars whose apocenters are just outside the Hill radius are
not instantaneously scattered – the impulse felt on the orbit from the second
black hole must build up over time before the star’s orbit is significantly
affected (Zhang et al., 2023).

It is also possible for stars to remain “Lagrange-stable” even when they are
no longer “Liapunov-stable” (e.g., Hayashi et al., 2022), meaning that stars
can perform chaotic figure eight orbits around both black holes but remain
bound to the system. In this situation, a star can become unbound from the
first SMBH only to return later and get disrupted.

Finally, three-body scattering experiments have shown that even when stars are
far outside the Hill radius, significant fractions of stars can be disrupted
through chaotic orbital crossings with the second SMBH (although this has
mostly been studied for stars orbiting the larger black hole, e.g. Chen et
al., 2011). Given that we are focused here on studying TDEs disrupted through
the EKL mechanism, we only consider disruptions within the hierarchical
radius, which is within a factor of two of the Hill’s radius in our
simulations with $q=10$ (and within a factor of four for simulations with
$q=100$), even for eccentricities approaching 1.

Thus, we conclude that using the Hill radius as a lower limit produces
estimates on the TDE rate that are too conservative, and the true rate is
somewhere in between the limits used in this paper. We note that a similar
approach was taken in Naoz et al. (2022) and Melchor et al. (2023). See the
former study’s appendix for an exploration of the particles’ trajectories.

##  3 Results

In this section we present the results of our dynamical simulations, focusing
on the rates and fraction of stars that become tidal disruption events.

###  3.1 Rates

Figure 5: TDE rates calculated for the EKL mechanism compared to rates
calculated for 2-body relaxation and to observed rates. The 2-body relaxation
rates plotted for ‘cuspy’ galaxies (dashed lines) correspond to galaxies with
an inner density profile ($\rho\propto r^{-\alpha}$) $\alpha>0.5$, while the
rates plotted for ‘core’ galaxies (dotted lines) correspond to $\alpha<0.3$
(Stone & Metzger, 2016, note that density profiles are estimated
observationally from radii $>r_{h}$, and therefore are not a direct comparison
to the density profiles used in simulations in this paper for radii interior
to $r_{h}$.). The observed rate for all TDEs (shaded in gray) is from van
Velzen et al. (2020), and the observed rate in post-starburst galaxies
(hatched region) is from French et al. (2020).

As you can see in Figure 5, we obtain the expected ‘burst’ of TDEs in all of
our simulations. The bursts peak around $\sim 10^{5}-10^{6}$ years but remain
higher than the global observed rate ($\lesssim 10^{-4}$ tdes/galaxy/year) for
$\sim 10^{6}-10^{8}$ years. This is comparable to or longer than the burst
timescales predicted in Chen et al. (2011) for the more massive black hole
companion (e.g. comparing our timescales for $10^{6}M_{\odot}$ secondaries to
Chen et al. (2011) for $10^{7}M_{\odot}$ primaries). The burst is limited by
the eccentric Kozai-Lidov timescale in that it will not begin until the stars’
eccentricities have grown large enough for disruptions to occur. However, the
burst can extend much longer than the eccentric Kozai-Lidov timescale (which
is generally of order $10^{4}-10^{6}$ years for these systems), as stars often
go through many EKL cycles before they are disrupted.

We also note that we do not include the evolution of the black hole binaries
in these simulations. If we modeled them as they evolved, we would expect the
EKL ‘loss-cone’ to move from larger radii to smaller radii, as stars at larger
radii are disrupted or scattered and the perturbations grow stronger at small
radii. Therefore, we would likely expect a more extended rate if we included
this binary evolution. This is what Chen et al. (2011) found in their n-body
simulations which included binary hardening due to stellar scattering (and
Mazzolari et al., 2022, for similar simulations focused on EMRIs). Initially,
while the disruption timescale was comparable to the hardening timescale, the
peak of the burst was extended, but then when the hardening timescale grew
longer than the tidal disruption timescale and the loss cone could no longer
be replenished, the rate dropped off (although the authors note that loss-cone
diffusion processes and gravitational wave emission were not included, and
would likely result in further hardening at later times and a longer-lived
plateau in the rate).

###  3.2 TDE fraction

The fraction of stars that become tidal disruption events depends on the
number of stars within the hierarchical limit and the strength of the EKL
mechanism i.e., the maximum value of the eccentricity excitations. In the test
particle limit, the eccentricity can be excited to extreme values (e.g., Li et
al., 2014; Naoz & Silk, 2014), naturally driving the stars to low angular
momentum orbits. We find that for the majority of our simulations, the upper
limit on the percentage of stars that are disrupted is $\approx 30-45\%$
(within the maximum radius defined by our hierarchical condition, Equation 2).
This limit does vary slightly with the parameters chosen.

As depicted in Figure 6, the fraction of stars disrupted decreases as $m_{1}$
increases, particularly at small radii. This is mostly because GR precession
is faster than the EKL precession close to the black hole. On the other hand,
the fraction of stars disrupted at small radii increases when $m_{2}$
increases because the strength of the EKL oscillations increases with the mass
of the perturber. The number of disruptions at large radii (close to the
hierarchical radius) does not vary significantly because we are effectively
running out of stars at those radii. Near the hierarchical radius, most of the
stars in the inclination range susceptible to EKL perturbations are getting
excited and disrupted, so increasing the strength of the perturbations is not
noticeably increasing the number of disruptions. As there are much fewer stars
at small radii ($M_{*,\rm contained}\propto r^{3-\alpha}$), increasing the
mass ratio to 100 only increases the upper limit on the total fraction of
disrupted stars in our simulations by $\approx 2\%$ (see Table 1).

Figure 6:  The fraction of stars that are disrupted as TDEs as a function of
semi-major axis (scaled by the hierarchical radius). Each panel varies one
simulation parameter. Top: $\alpha=1.5$, $q=10$. Middle: $m_{1}=10^{6}{\rm
M}_{\odot}$, $q=10$. Bottom: $m_{1}=10^{6}{\rm M}_{\odot}$, $\alpha=1.5$.

As expected, the density profile power law ($\alpha$) affects the maximum
fraction of TDEs, because changing $\alpha$ changes the relative fraction of
stars near the hierarchical radius, where most of the disruptions originate.
The number of stars disrupted is proportional to $1/t_{\rm KL}$, which is, in
turn, proportional to $a_{*}^{3/2}$. This can be used to estimate the upper
limit on the total fraction of stars disrupted. The number of stars disrupted
as a function of radius is dependent on both the eccentric Kozai-Lidov
timescale and the density profile, and therefore goes as

$n_{\rm disrupt}(r)\propto 1/t_{\rm KL}\times n_{*}(r)\propto
r^{3/2}r^{2-\alpha}.$ (10)

Then, the total fraction of stars disrupted goes as the integral of $n_{\rm
disrupt}(r)$ divided by the total number of stars

$\frac{\int{n_{\rm
disrupt}(r)dr}}{\int{n_{*}(r)dr}}\propto\frac{\int{r^{3/2}r^{2-\alpha}dr}}{\int{r^{2-\alpha}dr}}=\frac{r_{\rm
max}^{3/2}(3-\alpha)}{(9/2-\alpha)},$ (11)

and so the upper limit on the total fraction of stars disrupted (integrated
over the radius range) is proportional to:

$\rm{TDE\;fraction\;(upper\;limit)}\propto\frac{3-\alpha}{9/2-\alpha}.$ (12)

While changing the mass of both black holes and the steepness of the stellar
cusp affects the fraction of stars disrupted, it does not change the maximum
TDE fraction by more than $\sim 0.1$ (see Table 1). However, changing these
parameters does dramatically change the number of stars within the
hierarchical radius (our maximum radius, see Equation 1), and therefore the
number of stars that can become TDEs. Because of this, we find that the
strongest determinant of the maximum number of TDEs in the system is simply
the number of stars within the hierarchical radius ($N_{*}(r_{\rm max})\propto
r_{\rm max}^{9/2-\alpha}$)444We note that if we did not choose simulation
parameters such that $t_{KL}<t_{NT}$ and $t_{KL}<t_{GR1}$ for the majority of
stars, the steepness of the stellar cusp and the mass of the black holes might
have a more noticeable affect on the fraction of stars that become TDEs.. This
is shown clearly in Figure 7.

On the other hand, the lower limit on the number of TDEs is strongly dependent
on the mass ratio $q$. The lower limit on the number of TDEs is proportional
to the number of stars that are not susceptible to scattering by the perturber
– stars that remain within the Hill radius ($r_{\rm Hill}\propto q^{-1/3}$).
The number of stars within the Hill radius follows the proportionality
$N_{*}(r_{\rm Hill})\propto r_{\rm Hill}^{9/2-\alpha}$.

Combined with the fact that $r_{\rm Hill}\propto q^{-1/3}$, the lower limit on
the TDE fraction goes as555We note that the TDE fractions are scaled by the
number of stars within $r_{\rm hier}$, and $r_{\rm hier}$ is not dependent on
the mass ratio ‘q’. If we instead calculated the lower limit on the TDE
fraction by scaling by the number of stars within $r_{\rm Hill}$, the fraction
would go as $\propto r_{\rm Hill}^{3/2}\propto q^{-1/2}$ instead.,

$\rm{TDE\;fraction\;(lower\;limit)}\propto q^{\frac{\alpha}{3}-\frac{3}{2}}\
.$ (13)

It drops from $\approx 10\%$ for $q=10$ to $\approx 1\%$ for $q=100$
($m_{1}=10^{6}$, $\alpha=1.5$). We note that we find similar qualitative
results to that of Li et al. (2015) (see following section), however, here we
explore a wider part of the parameter space.

As we discussed in the previous section, our lower limit on the stability
radius is likely too conservative. Our lower limit is approximately equal to
the strongest limits found in Grishin et al. (2017) and Tory et al. (2022)
($a_{*}\lesssim 0.5r_{H}$). However, some of the stars outside this radius are
disrupted on very short timescales, likely before their orbits would become
unstable. Recently, Zhang et al. (2023) determined an analytical stability
criterion for 3-body systems that they found to be consistent with n-body
simulations. Adopting the Zhang et al. (2023) stability timescale definition
means that stars from larger radii can undergo EKL eccentricity excitations
before becoming unstable. This would increase the lower limit on the number of
TDEs by a factor of $\approx 2$ across all of our simulations. Therefore, it
is likely that many of these stars will be disrupted before their orbits
become unstable. It is also likely that many of the stars on unstable orbits
will also become TDEs, even though analytical approaches are unable to model
their evolution (as discussed in Section 2.3).

###  3.3 Comparison with previous work

Figure 7: LHS: The number of TDEs as a function of $m_{1}$, with different
density profiles represented as different marker shapes (q=10 in all cases).
The points in orange are the lower limits, excluding all stars that move
outside $r_{\rm Hill}$, while the points in blue are the upper limits and
include these stars. The blue and orange shaded regions are plotted by
determining the slopes between the fraction of TDEs for the $\alpha=1\;\&\;2$
runs at the black hole masses used in the simulations and then extrapolating
out to higher and lower masses. The black hatched region is the number of TDEs
disrupted by $m_{2}$ in the same system from chaotic scatterings (the upper
limit is from Liu et al. 2013 and the lower limit is scaled down to match Wegg
& Nate Bode 2011). The gray shaded region is the number of TDEs predicted from
2-body relaxation over $10^{8}$ years for a single SMBH system with mass
$m_{1}$ (from Stone & Metzger, 2016). RHS: The number of TDEs as a function of
the number of stars within the hierarchical radius. Each $m_{1}$ and $\alpha$
combination is plotted as a different color, and the marker shape is
determined by the mass ratio $q$. The upper and lower limits are determined in
the same way as for the LHS.

We compared our results with Li et al. (2015), as they used a similar setup to
our own. In their simulation, $m_{1}=10^{7}{\rm M}_{\odot}$, $m_{2}=10^{8}{\rm
M}_{\odot}$, $\alpha=1.75$, $e_{2}=0.5$ and $a_{2}=0.5$pc, and they find 5.7%
of stars within the Hill’s radius are disrupted (using the more lenient
criteria that a star must only pass within $3R_{T}$ to be disrupted). This is
comparable to the results for our simulations numbered 5-7, which were run
with similar parameters and found that $\approx 5\%$ of stars within the
Hill’s radius are disrupted. To compare directly, we also ran one additional
test simulation with the same parameters (and disruption criteria) as Li et
al. (2015), and found that 4.2% of stars within the Hill’s radius were
disrupted. We note that Li et al. (2015) include the effect of orbital
precession due to the stellar cusp in their simulations but that our
simulations appear consistent nonetheless, justifying our exclusion of this
effect for this set of simulations (see Section 2 for further explanation).

We also compared our results to the 3-body scattering simulations of the
eccentric Kozai-Lidov mechanism around the larger SMBH in Chen et al. (2011)
and Wegg & Nate Bode (2011). Both of these works showed that the EKL
oscillations induced by a smaller SMBH on stars surrounding a larger companion
are suppressed for large parts of the parameter space. In Figure 5, we show
that we get a comparable number of TDEs from the smaller black hole in our
$q=10$ simulations as predicted for the larger black hole by scaling relations
calculated from the results in Chen et al. (2011) and Wegg & Nate Bode (2011).

##  4 Discussion

We find that the EKL mechanism in SMBH binaries produces a burst of tidal
disruption events at rates much higher than those expected from two-body
relaxation (see Figure 5). This phenomenon has been previously explored for
tidal disruption events, but the majority of past work has focused on the
disruptions of stars orbiting the more massive black hole in the binary system
(e.g. Ivanov et al., 2005; Chen et al., 2009, 2011; Wegg & Nate Bode, 2011).
This meant that the majority of disruptions did not actually come from
eccentric Kozai-Lidov oscillations, but rather from chaotic orbital crossings
between the stars and the perturbing black hole (Chen et al., 2011; Wegg &
Nate Bode, 2011). As we are interested in determining when eccentric Kozai-
Lidov oscillations in particular are important, we followed Li et al. (2015)’s
lead and explored a parameter space centered on the smaller black hole.

###  4.1 Connections to simulated and observed TDE rates

In Figure 5 we compare the rate calculated from our simulations with the rates
estimated for two-body relaxation based on properties of observed galaxies
(dashed lines, from Stone & Metzger, 2016). We find that the rate from the EKL
mechanism is higher than the rate from two-body relaxation for timescales
between millions to tens of millions of years depending on the parameters of
the system. This means that during a burst the EKL rate is also much higher
than the overall observed rate (as two-body relaxation calculations generally
over-predict the observed rate). However, the scenario described in this paper
is only applicable to a subset of galaxies, and only over a relatively short
time period.

####  4.1.1 Caveats

To begin with, only a fraction of galaxies will host close SMBH binaries at
any given time. Given the uncertainties on SMBHB hardening timescales (e.g.
Kelley et al., 2017), it is not straightforward to estimate the current
population of SMBH binaries from the rate of galaxy mergers, however, an
overview is still useful for understanding possible scenarios.

The current day merger fraction (as defined by the fraction of galaxies with
close companions) is $\lesssim 3\%$ for galaxies $\lesssim 10^{11}M_{\odot}$
(this is the major merger fraction, the minor merger fraction is likely
similar, Bundy et al., 2009; López-Sanjuan et al., 2012). However, it is
possible that the current-day SMBHB population reflects a higher merger
fraction from earlier times, depending on the hardening timescale of the
binary (Bluck et al. 2009, 2012 find a major merger fraction of up to $\approx
25-40\%$ at $z=3$, which is $\approx 10^{10}$ year lookback time). If the
SMBHB population is reflective of these higher merger rates, upwards of 25% of
galaxies might host two SMBHs at a given time (Rodriguez-Gomez et al., 2015),
but their relevance to this work would still depend on their hardening
timescales.

In addition to requiring that the host galaxies and their SMBHs be at a
specific evolutionary stage, the actual burst of TDEs from the EKL mechanism
is short-lived. If the stars are not replenished rapidly through star
formation or dynamical mechanisms, the rate will drop off relatively quickly
and the elevated rate will only be observable for tens of millions of years
for anything but the most massive systems (see Figure 5). Therefore, it is
possible that even in galaxies with close SMBH binaries, we will miss the
‘burst’ of TDEs predicted here.

Finally, we note that although we focus on the tidal disruptions around the
smaller of the black holes in the SMBH binary, the global rate will be limited
by the number of larger mass SMBHs available to form binaries because the
black hole mass function decreases with increasing mass (e.g. Shankar et al.,
2016; Gallo & Sesana, 2019).

####  4.1.2 Comparison with two-body relaxation

Given these restrictions, it is reasonable that we do not observe dramatically
higher TDE rates due to the EKL mechanism. However, we still expect this
mechanism to contribute meaningfully to TDE rates, particularly for higher
mass systems with lower rates from two-body relaxation. In Figure 7, the
number of TDEs from two-body relaxation from a single galaxy is comparable to
the number from an EKL burst over 100 million years for a disrupting black
hole mass of $10^{6}$ (perturbing $M_{h}=10^{7}M_{\odot}$). At higher black
hole masses, the EKL mechanism dominates even over these relatively long
timescales. While this comparison is still during the bursting phase, Stone et
al. (2018) estimated the average TDE rates per galaxy for SMBHBs (assuming
that the larger black hole dominated the rates) and compared it to two-body
relaxation rates. They found that the average rate from SMBHBs approaches the
rate from two-body relaxation for black hole masses $\gtrsim
10^{7.5}M_{\odot}$. When we compare the number of disruptions from our
simulations with the number of disruptions expected around the larger black
hole (from chaotic orbital interactions), we find that the numbers are similar
for $q=10$ (see hatched region in Figure 7, Chen et al., 2011; Wegg & Nate
Bode, 2011). Assuming the smaller black hole contributes a similar number of
disruptions would mean that the average rate from SMBHBs would approach the
rate from two-body relaxation for primary SMBH masses around
$10^{7}M_{\odot}$. In addition to increasing the total number of disruptions
from each SMBHB pair, it would also allow for the inclusion of higher mass
galaxy mergers where the larger black hole is past the Hill mass turnover
point ($\approx 10^{8}M\odot$).

####  4.1.3 Comparison with post-starburst rates

We note that the rate we calculate during the EKL burst is much more
comparable to the observed rates calculated for post-starburst galaxies, which
are often considered candidates for recent mergers (Mihos & Hernquist, 1994;
Bekki et al., 2005; Hopkins et al., 2009, as mergers can trigger starbursts).
The reason why the TDE rate in post-starburst galaxies is $\sim 20-200\times$
higher than the galaxy averaged rate remains an open question in the field
(e.g., Arcavi et al., 2014; Law-Smith et al., 2017; French et al., 2020; Dodd
et al., 2021). Previous work has suggested that SMBH binaries might be able to
increase the rates of TDEs in these galaxies if the starburst is due to a
recent merger (e.g. Arcavi et al., 2014). However, as Stone et al. (2018)
point out, if the starburst is triggered by a merger, and the increased rates
are from the resulting SMBH binary, this requires the smaller black hole to
sink to the center of the galaxy on timescales short enough that the post-
starburst features are still observable. If the starburst is triggered when
the merger starts, this may be difficult, as the timescales for most SMBH
binaries to harden to distances of $\approx 1$pc are thought to be of order
$10^{9}$ years (Kelley et al., 2017). Additionally, the observed TDEs from
post-starburst galaxies appear to have SMBH masses measured from light curves
that are consistent with the SMBH masses from galaxy scaling relations
(Mockler et al., 2019; Ramsden et al., 2022). Therefore, even if the galaxies
host SMBH binaries, it is unlikely that the smaller black hole is responsible
for the disruptions and therefore that the EKL mechanism is producing these
TDEs.

###  4.2 Finding hidden SMBH binaries

Despite their rarity, bursts of TDEs produced by the EKL mechanism could help
us find hidden SMBH binaries. Tidal features from minor mergers dissipate
quickly (Conselice et al., 2000), and the gravitational potential of a
galactic nucleus hosting a SMBH binary with a mass ratio $\gtrsim 10$ will be
dominated by the mass of the larger black hole, making it difficult to find
SMBH binary candidates. However, as shown here, disruptions produced by the
EKL mechanism will produce a significant number of TDEs around the smaller
black hole in a SMBH binary, comparable in number to TDEs around the larger
black hole for mass ratios of $q=10$ (see Figure 7). Tidal disruption flares
contain information about the mass of the disrupting black hole, therefore
TDEs from the smaller black hole can expose its presence. For example, the
light curve timescale of a TDE is dependent on the mass of the black hole
(e.g. $t\propto M_{h}^{1/2}$, Mockler et al., 2019), therefore if a very short
flare occurred around what seemed to be a much larger black hole, that would
be an indication of a potential SMBH binary. If the mass of the larger black
hole is above $\gtrsim 10^{8}{\rm M}_{\odot}$, it would not be able to produce
an observable flare from the disruptions of most stars, and therefore seeing
any TDE at all from very massive galaxies would be a very strong indication of
a hidden companion SMBH (e.g. Coughlin & Armitage, 2018; Fragione & Leigh,
2018).

##  5 Summary of Key Findings

Our key results are summarized as follows:

  * •

The EKL mechanism produces a burst of TDEs around the smaller black hole in a
SMBHB at a rate that is orders of magnitude higher than the rate from two-body
relaxation for timescales of millions to tens of millions of years (see Figure
5). Unlike EKL oscillations around the larger SMBH, the EKL oscillations
around the smaller SMBH are not significantly suppressed by precession due to
general relativity or the stellar cusp. This is mostly because of the
increased strength of the perturber.

  * •

The number of disruptions is primarily dependent on the number of stars within
the hierarchical (or Hill) radius. This is because the fraction of stars
disrupted does not change dramatically with most parameters (so long as we
choose regions of parameter space such that the EKL mechanism is the dominant
dynamical mechanism, see Figure 7, Figure 3, Table 1). The timescale is most
strongly dependent on the black hole mass and mass ratio (see Figure 5).

  * •

The EKL mechanism generates a larger fraction of disruptions around the
smaller black hole ($\sim 3\%-45\%$ of stars within the hierarchical radius
are disrupted) than chaotic orbital interactions are able to produce around
the larger black hole (e.g. Chen et al., 2011; Wegg & Nate Bode, 2011). This
can lead to a comparable number of disruptions around both black holes (see
Figure 7).

  * •

Disruptions around the smaller black hole can provide evidence of the
existence of the SMBH binary as their light curve timescales will likely be
shorter than what is expected from the larger black hole, which dominates the
potential (Mockler et al., 2019).

We have shown that tidal disruptions from the secondary black hole in a SMBHB
contribute meaningfully to the rate from the binary system. We emphasize that
this provides a novel way to find SMBHB candidates at separations that are too
small to resolve distinct sources.

We thank the participants and organizers of the fall 2023 Como meeting on the
dynamics of TDEs and EMRIs for fruitful discussions and helpful clarification.
B.M. is grateful for support from the U.C. Chancellor’s Postdoctoral
fellowship and the Mani L. Bhaumik Institute for Theoretical Physics as well
as from Swift (80NSSC21K1409). D.M. is grateful for support from the NSF
graduate fellowship. S.N. acknowledges the partial support from NASA ATP
80NSSC20K0505 and from NSF-AST 2206428 grant as well as thanks Howard and
Astrid Preston for their generous support. E.R.R. thanks the Heising-Simons
Foundation, NSF (AST-1615881 and AST-2206243), Swift (80NSSC21K1409,
80NSSC19K1391) and Chandra (22-0142) for support. The calculations for this
work were carried out in part on the UCSC lux supercomputer (supported by NSF
MRI grant AST-1828315).

##  Appendix A Dynamical timescales

We present an overview of the relevant dynamical timescales for our SMBH
binary systems. These were used in the making of Figure 3.

$\displaystyle t_{\rm quad}=\frac{16}{15}\frac{a_{\rm bin}^{3}(1-e_{\rm
bin}^{2})^{3/2}\sqrt{(m_{1}+m_{*})}}{\sqrt{G}a_{*}^{3/2}m_{2}}$ (A1)
$\displaystyle t_{\rm GR,\>1}=\frac{2\pi
a_{*}^{5/2}c^{2}(1-e_{*}^{2})}{3G^{3/2}(m_{1}+m_{*})^{3/2}}$ (A2)
$\displaystyle t_{\rm GR,\>2}=\frac{2\pi a_{\rm bin}^{5/2}c^{2}(1-e_{\rm
bin}^{2})}{3G^{3/2}(m_{*}+m_{1}+m_{2})^{3/2}}$ (A3) $\displaystyle t_{\rm
GR,\>int}=\frac{16}{9}\frac{a_{\rm bin}^{3}c^{2}(1-e_{\rm
bin}^{2})^{3/2}}{\sqrt{a_{*}}e_{*}\sqrt{1-e_{*}^{2}}G^{3/2}m_{1}^{1/2}m_{2}}$
(A4) $\displaystyle t_{\rm NT}=2\pi\left(\frac{\sqrt{Gm_{1}/a_{*}^{3}}}{\pi
m_{1}e_{*}}\int_{0}^{\pi}{\rm d}\psi\>M_{\ast}(r)\>{\rm cos}\psi\right)^{-1}$
(A5) $\displaystyle t_{\rm RR,\>v}=\frac{2\pi f_{\rm
vrr}}{\Omega}m_{1}\frac{1}{\sqrt{M_{\ast}(r)m_{*}}}$ (A6) $\displaystyle
t_{\rm rel}=0.34\frac{\sigma_{\ast}^{3}}{G^{2}\rho m_{*}{\rm ln}\Lambda}$ (A7)
$\displaystyle t_{\rm
LT}=\frac{a_{*}^{3}c^{3}(1-e_{*}^{2})^{3/2}}{2G^{2}m_{1}^{2}s}$ (A8)
$\displaystyle t_{\rm GW}=\frac{a_{\rm
bin}^{4}}{4}\frac{5}{64}\frac{c^{5}}{G^{3}m_{1}m_{2}(m_{1}+m_{2})}$ (A9)

The equations above represent timescales for the following physical processes:
$t_{\rm quad}$ is the quadrupole eccentric kozai-lidov timescale (Naoz et al.,
2013; Antognini, 2015); $t_{\rm GR,\>1}$ & $t_{\rm GR,\>2}$ are the first-
order post-Newtonian GR precession timescales for the inner and outer black
holes; and $t_{\rm GR,\>int}$ is the timescale of the first order post-
Newtonian GR interaction between the inner and outer orbits. As defined in
Kocsis & Tremaine (2011), $t_{\rm NT}$ is the Newtonian precession timescale
caused by the potential of the stellar cusp, $t_{\rm RR,\>v}$ is the vector
resonant relaxation timescale, and $t_{\rm LT}$ is the Lense-Thirring
precession timescale. Finally, $t_{\rm GW}$ is the binary SMBH orbital decay
timescale due to gravitational wave emission. In equation A6, $\Omega$ is the
orbital frequency of the star and $f_{vrr}$ is estimated to be 1.2 by Kocsis &
Tremaine (2015). In equation A5, $\psi$ is the true anomaly of the inner
orbit, $r(\psi)=a_{1}(1-e_{1}^{2})/(1+e_{1}cos\psi)$ from Kepler’s equation.
The function $M_{\ast}(r)$ is the integrated stellar mass within $r(\psi)$. In
equation A7, $\rm ln\Lambda$ is the Coulomb logarithm, and $\sigma_{\ast}$ is
the velocity dispersion of the bulge. In equation A8, the constant
$(G^{2}m_{0}^{2}s)/c$ is the spin angular momentum of the inner SMBH (Peters,
1964; Kocsis & Tremaine, 2011; Naoz et al., 2013).

## References

  * Angus et al. (2022) Angus, C. R., Baldassare, V. F., Mockler, B., et al. 2022, Nature Astronomy, 6, 1452, doi: 10.1038/s41550-022-01811-y
  * Antognini (2015) Antognini, J. M. O. 2015, Monthly Notices of the Royal Astronomical Society, 452, 3610, doi: 10.1093/mnras/stv1552
  * Arcavi et al. (2014) Arcavi, I., Gal-Yam, A., Sullivan, M., et al. 2014, ApJ, 793, 38, doi: 10.1088/0004-637X/793/1/38
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Barnes (2002) Barnes, J. E. 2002, Monthly Notices of the Royal Astronomical Society, 333, 481, doi: 10.1046/j.1365-8711.2002.05335.x
  * Batcheldor et al. (2010) Batcheldor, D., Robinson, A., Axon, D. J., Perlman, E. S., & Merritt, D. 2010, ApJ, 717, L6, doi: 10.1088/2041-8205/717/1/L6
  * Bekki et al. (2005) Bekki, K., Couch, W. J., Shioya, Y., & Vazdekis, A. 2005, MNRAS, 359, 949, doi: 10.1111/j.1365-2966.2005.08932.x
  * Bhaskar et al. (2021) Bhaskar, H., Li, G., Hadden, S., Payne, M. J., & Holman, M. J. 2021, The Astronomical Journal, 161, 48, doi: 10.3847/1538-3881/abcbfc
  * Bhattacharya et al. (2020) Bhattacharya, A., Nehal C., P., Das, M., et al. 2020, Automated Detection of Double Nuclei Galaxies using GOTHIC and the Discovery of a Large Sample of Dual AGN, Tech. rep., doi: 10.48550/arXiv.2011.12177
  * Bianchi et al. (2008) Bianchi, S., Chiaberge, M., Piconcelli, E., Guainazzi, M., & Matt, G. 2008, Monthly Notices of the Royal Astronomical Society, 386, 105, doi: 10.1111/j.1365-2966.2008.13078.x
  * Bluck et al. (2009) Bluck, A. F. L., Conselice, C. J., Bouwens, R. J., et al. 2009, Monthly Notices of the Royal Astronomical Society, 394, L51, doi: 10.1111/j.1745-3933.2008.00608.x
  * Bluck et al. (2012) Bluck, A. F. L., Conselice, C. J., Buitrago, F., et al. 2012, The Astrophysical Journal, 747, 34, doi: 10.1088/0004-637X/747/1/34
  * Bogdanović et al. (2009) Bogdanović, T., Eracleous, M., & Sigurdsson, S. 2009, The Astrophysical Journal, 697, 288, doi: 10.1088/0004-637X/697/1/288
  * Boroson & Lauer (2009) Boroson, T. A., & Lauer, T. R. 2009, Nature, 458, 53, doi: 10.1038/nature07779
  * Bundy et al. (2009) Bundy, K., Fukugita, M., Ellis, R. S., et al. 2009, The Astrophysical Journal, 697, 1369, doi: 10.1088/0004-637X/697/2/1369
  * Casey-Clyde et al. (2022) Casey-Clyde, J. A., Mingarelli, C. M. F., Greene, J. E., et al. 2022, The Astrophysical Journal, 924, 93, doi: 10.3847/1538-4357/ac32de
  * Chen & Liu (2013) Chen, X., & Liu, F. K. 2013, The Astrophysical Journal, 762, 95, doi: 10.1088/0004-637X/762/2/95
  * Chen et al. (2009) Chen, X., Madau, P., Sesana, A., & Liu, F. K. 2009, ApJ, 697, L149.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2009ApJ...697L.149C&link_type=ABSTRACT
  * Chen et al. (2011) Chen, X., Sesana, A., Madau, P., & Liu, F. K. 2011, ApJ, 729, 13.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2011ApJ...729...13C&link_type=ABSTRACT
  * Comerford et al. (2009) Comerford, J. M., Griffith, R. L., Gerke, B. F., et al. 2009, The Astrophysical Journal, 702, L82, doi: 10.1088/0004-637X/702/1/L82
  * Comerford et al. (2018) Comerford, J. M., Nevin, R., Stemo, A., et al. 2018, The Astrophysical Journal, 867, 66, doi: 10.3847/1538-4357/aae2b4
  * Conselice et al. (2000) Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, The Astrophysical Journal, 529, 886, doi: 10.1086/308300
  * Coughlin & Armitage (2018) Coughlin, E. R., & Armitage, P. J. 2018, Monthly Notices of the Royal Astronomical Society, 474, 3857, doi: 10.1093/mnras/stx3039
  * De Rosa et al. (2019) De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, New Astronomy Reviews, 86, 101525, doi: 10.1016/j.newar.2020.101525
  * Deane et al. (2014) Deane, R. P., Paragi, Z., Jarvis, M. J., et al. 2014, Nature, 511, 57, doi: 10.1038/nature13454
  * Decarli et al. (2010) Decarli, R., Dotti, M., Montuori, C., Liimets, T., & Ederoclite, A. 2010, The Astrophysical Journal, 720, L93, doi: 10.1088/2041-8205/720/1/L93
  * Dey et al. (2019) Dey, L., Gopakumar, A., Valtonen, M., et al. 2019, Universe, 5, 108, doi: 10.3390/universe5050108
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604, doi: 10.1038/nature03335
  * Dodd et al. (2021) Dodd, S. A., Law-Smith, J. A. P., Auchettl, K., Ramirez-Ruiz, E., & Foley, R. J. 2021, ApJ, 907, L21, doi: 10.3847/2041-8213/abd852
  * Dotti et al. (2009) Dotti, M., Montuori, C., Decarli, R., et al. 2009, Monthly Notices of the Royal Astronomical Society, 398, L73, doi: 10.1111/j.1745-3933.2009.00714.x
  * Eracleous et al. (2012) Eracleous, M., Boroson, T. A., Halpern, J. P., & Liu, J. 2012, The Astrophysical Journal Supplement Series, 201, 23, doi: 10.1088/0067-0049/201/2/23
  * Evans & Kochanek (1989) Evans, C. R., & Kochanek, C. S. 1989, ApJ, 346, L13.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=1989ApJ...346L..13E&link_type=ABSTRACT
  * Foord et al. (2020) Foord, A., Gültekin, K., Nevin, R., et al. 2020, The Astrophysical Journal, 892, 29, doi: 10.3847/1538-4357/ab72fa
  * Fragione & Leigh (2018) Fragione, G., & Leigh, N. 2018, Monthly Notices of the Royal Astronomical Society, 479, 3181, doi: 10.1093/mnras/sty1600
  * Fragione et al. (2020) Fragione, G., Loeb, A., Kremer, K., & Rasio, F. A. 2020, The Astrophysical Journal, 897, 46, doi: 10.3847/1538-4357/ab94b2
  * Fragione et al. (2021) Fragione, G., Perna, R., & Loeb, A. 2021, Monthly Notices of the Royal Astronomical Society, 500, 4307, doi: 10.1093/mnras/staa3493
  * French et al. (2020) French, K. D., Arcavi, I., Zabludoff, A. I., et al. 2020, The Astrophysical Journal, 891, 93, doi: 10.3847/1538-4357/ab7450
  * French et al. (2020) French, K. D., Arcavi, I., Zabludoff, A. I., et al. 2020, ApJ, 891, 93, doi: 10.3847/1538-4357/ab7450
  * Gallo & Sesana (2019) Gallo, E., & Sesana, A. 2019, The Astrophysical Journal, 883, L18, doi: 10.3847/2041-8213/ab40c6
  * Gaskell (1983) Gaskell, C. M. 1983, 24, 473.  https://ui.adsabs.harvard.edu/abs/1983LIACo..24..473G
  * Generozov & Madigan (2020) Generozov, A., & Madigan, A.-M. 2020, The Astrophysical Journal, 896, 137, doi: 10.3847/1538-4357/ab94bc
  * Glampedakis & Kennefick (2002) Glampedakis, K., & Kennefick, D. 2002, Physical Review D, 66, 044002, doi: 10.1103/PhysRevD.66.044002
  * Goulding et al. (2019) Goulding, A. D., Pardo, K., Greene, J. E., et al. 2019, The Astrophysical Journal, 879, L21, doi: 10.3847/2041-8213/ab2a14
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2020, Astronomy and Astrophysics, 636, L5, doi: 10.1051/0004-6361/202037813
  * Greene et al. (2010) Greene, J. E., Peng, C. Y., Kim, M., et al. 2010, The Astrophysical Journal, 721, 26.  http://stacks.iop.org/0004-637X/721/i=1/a=26
  * Greene et al. (2021) Greene, J. E., Lancaster, L., Ting, Y.-S., et al. 2021, The Astrophysical Journal, 917, 17, doi: 10.3847/1538-4357/ac0896
  * Grishin et al. (2017) Grishin, E., Perets, H. B., Zenati, Y., & Michaely, E. 2017, Monthly Notices of the Royal Astronomical Society, 466, 276, doi: 10.1093/mnras/stw3096
  * Gualandris & Merritt (2009) Gualandris, A., & Merritt, D. 2009, The Astrophysical Journal, 705, 361, doi: 10.1088/0004-637X/705/1/361
  * Guillochon (2013) Guillochon, J. 2013, PhD thesis, University of California, Santa Cruz
  * Guillochon & Ramirez-Ruiz (2013) Guillochon, J., & Ramirez-Ruiz, E. 2013, The Astrophysical Journal, 767, 25, doi: 10.1088/0004-637X/767/1/25
  * Guillochon & Ramirez-Ruiz (2015) Guillochon, J., & Ramirez-Ruiz, E. 2015, ApJ, 809, 166, doi: 10.1088/0004-637X/809/2/166
  * Gürkan & Rasio (2005) Gürkan, M. A., & Rasio, F. A. 2005, The Astrophysical Journal, 628, 236, doi: 10.1086/430694
  * Hammerstein et al. (2022) Hammerstein, E., van Velzen, S., Gezari, S., et al. 2022, arXiv:2203.01461 [astro-ph].  http://arxiv.org/abs/2203.01461
  * Hansen & Milosavljević (2003) Hansen, B. M. S., & Milosavljević, M. 2003, The Astrophysical Journal, 593, L77, doi: 10.1086/378182
  * Hayashi et al. (2022) Hayashi, T., Trani, A. A., & Suto, Y. 2022, Lagrange vs. Liapunov stability of hierarchical triple systems: dependence on the mutual inclination between inner and outer orbits, Tech. rep.  https://ui.adsabs.harvard.edu/abs/2022arXiv220908487H
  * Healy et al. (2009) Healy, J., Levin, J., & Shoemaker, D. 2009, Physical Review Letters, 103, 131101, doi: 10.1103/PhysRevLett.103.131101
  * Hopkins et al. (2008) Hopkins, P. F., Cox, T. J., Kereš, D., & Hernquist, L. 2008, The Astrophysical Journal Supplement Series, 175, 390, doi: 10.1086/524363
  * Hopkins et al. (2009) Hopkins, P. F., Cox, T. J., Younger, J. D., & Hernquist, L. 2009, ApJ, 691, 1168, doi: 10.1088/0004-637X/691/2/1168
  * Ivanov et al. (2005) Ivanov, P. B., Polnarev, A. G., & Saha, P. 2005, MNRAS, 358, 1361.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2005MNRAS.358.1361I&link_type=ABSTRACT
  * Izquierdo-Villalba et al. (2022) Izquierdo-Villalba, D., Sesana, A., & Colpi, M. 2022, Unveiling the hosts of parsec-scale massive black hole binaries: morphology and electromagnetic signatures, Tech. rep.  https://ui.adsabs.harvard.edu/abs/2022arXiv220704064I
  * Katz et al. (2011) Katz, B., Dong, S., & Malhotra, R. 2011, Physical Review Letters, 107, 181101, doi: 10.1103/PhysRevLett.107.181101
  * Kelley et al. (2017) Kelley, L. Z., Blecha, L., & Hernquist, L. 2017, Monthly Notices of the Royal Astronomical Society, 464, 3131, doi: 10.1093/mnras/stw2452
  * Kim et al. (2020) Kim, D. C., Yoon, I., Evans, A. S., et al. 2020, The Astrophysical Journal, 904, 23, doi: 10.3847/1538-4357/abb9a0
  * Kocsis & Tremaine (2011) Kocsis, B., & Tremaine, S. 2011, Monthly Notices of the Royal Astronomical Society, 412, 187, doi: 10.1111/j.1365-2966.2010.17897.x
  * Kocsis & Tremaine (2015) —. 2015, Monthly Notices of the Royal Astronomical Society, 448, 3265, doi: 10.1093/mnras/stv057
  * Komossa et al. (2004) Komossa, S., Halpern, J., Schartel, N., et al. 2004, ApJ, 603, L17.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2004ApJ...603L..17K&link_type=ABSTRACT
  * Komossa et al. (2008) Komossa, S., Zhou, H., & Lu, H. 2008, The Astrophysical Journal, 678, L81, doi: 10.1086/588656
  * Kozai (1962) Kozai, Y. 1962, The Astronomical Journal, 67, 591, doi: 10.1086/108790
  * Law-Smith et al. (2019) Law-Smith, J., Guillochon, J., & Ramirez-Ruiz, E. 2019, The Astrophysical Journal, 882, L25, doi: 10.3847/2041-8213/ab379a
  * Law-Smith et al. (2017) Law-Smith, J., MacLeod, M., Guillochon, J., Macias, P., & Ramirez-Ruiz, E. 2017, ApJ, 841, 132, doi: 10.3847/1538-4357/aa6ffb
  * Law-Smith et al. (2017) Law-Smith, J., Ramirez-Ruiz, E., Ellison, S. L., & Foley, R. J. 2017, The Astrophysical Journal, 850, 22, doi: 10.3847/1538-4357/aa94c7
  * Li et al. (2014) Li, G., Naoz, S., Kocsis, B., & Loeb, A. 2014, The Astrophysical Journal, 785, 116, doi: 10.1088/0004-637X/785/2/116
  * Li et al. (2015) Li, G., Naoz, S., Kocsis, B., & Loeb, A. 2015, MNRAS, 451, 1341, doi: 10.1093/mnras/stv1031
  * Li et al. (2020) Li, K., Bogdanović, T., & Ballantyne, D. R. 2020, The Astrophysical Journal, 896, 113, doi: 10.3847/1538-4357/ab93c6
  * Li et al. (2019) Li, S., Berczik, P., Chen, X., et al. 2019, ApJ, 883, 132, doi: 10.3847/1538-4357/ab3e4a
  * Lidov (1962) Lidov, M. L. 1962, Planetary and Space Science, 9, 719, doi: 10.1016/0032-0633(62)90129-0
  * Lithwick & Naoz (2011) Lithwick, Y., & Naoz, S. 2011, The Astrophysical Journal, 742, 94, doi: 10.1088/0004-637X/742/2/94
  * Liu et al. (2013) Liu, S.-F., Guillochon, J., Lin, D. N. C., & Ramirez-Ruiz, E. 2013, ApJ, 762, 37\.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2013ApJ...762...37L&link_type=ABSTRACT
  * Liu et al. (2010) Liu, X., Greene, J. E., Shen, Y., & Strauss, M. A. 2010, The Astrophysical Journal, 715, L30, doi: 10.1088/2041-8205/715/1/L30
  * Lodato et al. (2009) Lodato, G., King, A. R., & Pringle, J. E. 2009, MNRAS, 392, 332.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2009MNRAS.392..332L&link_type=ABSTRACT
  * López-Sanjuan et al. (2012) López-Sanjuan, C., Le Fèvre, O., Ilbert, O., et al. 2012, Astronomy and Astrophysics, 548, A7, doi: 10.1051/0004-6361/201219085
  * MacLeod et al. (2012) MacLeod, M., Guillochon, J., & Ramirez-Ruiz, E. 2012, ApJ, 757, 134, doi: 10.1088/0004-637X/757/2/134
  * Magorrian & Tremaine (1999) Magorrian, J., & Tremaine, S. 1999, Monthly Notices of the Royal Astronomical Society, 309, 447, doi: 10.1046/j.1365-8711.1999.02853.x
  * Maillard et al. (2004) Maillard, J. P., Paumard, T., Stolovy, S. R., & Rigaut, F. 2004, Astronomy and Astrophysics, 423, 155, doi: 10.1051/0004-6361:20034147
  * Mardling & Aarseth (2001) Mardling, R. A., & Aarseth, S. J. 2001, MNRAS, 321, 398, doi: 10.1046/j.1365-8711.2001.03974.x
  * Mazzolari et al. (2022) Mazzolari, G., Bonetti, M., Sesana, A., et al. 2022, Extreme Mass Ratio Inspirals triggered by Massive Black Hole Binaries: from Relativistic Dynamics to Cosmological Rates, Tech. rep.  https://ui.adsabs.harvard.edu/abs/2022arXiv220405343M
  * Melchor et al. (2023) Melchor, D., Mockler, B., Naoz, S., Rose, S. C., & Ramirez-Ruiz, E. 2023, submitted to ApJ
  * Mihos & Hernquist (1994) Mihos, J. C., & Hernquist, L. 1994, ApJ, 425, L13, doi: 10.1086/187299
  * Mockler et al. (2019) Mockler, B., Guillochon, J., & Ramirez-Ruiz, E. 2019, The Astrophysical Journal, 872, 151, doi: 10.3847/1538-4357/ab010f
  * Mockler et al. (2022) Mockler, B., Twum, A. A., Auchettl, K., et al. 2022, The Astrophysical Journal, 924, 70, doi: 10.3847/1538-4357/ac35d5
  * Monageng et al. (2021) Monageng, I. M., Motta, S. E., Fender, R., et al. 2021, Monthly Notices of the Royal Astronomical Society, 501, 5776, doi: 10.1093/mnras/stab043
  * Naoz (2016) Naoz, S. 2016, Annual Review of Astronomy and Astrophysics, 54, 441, doi: 10.1146/annurev-astro-081915-023315
  * Naoz et al. (2011) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187, doi: 10.1038/nature10076
  * Naoz et al. (2013) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, MNRAS, 431, 2155, doi: 10.1093/mnras/stt302
  * Naoz et al. (2013) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, Monthly Notices of the Royal Astronomical Society, 431, 2155, doi: 10.1093/mnras/stt302
  * Naoz et al. (2013) Naoz, S., Kocsis, B., Loeb, A., & Yunes, N. 2013, ApJ, 773, 187, doi: 10.1088/0004-637X/773/2/187
  * Naoz et al. (2013) Naoz, S., Kocsis, B., Loeb, A., & Yunes, N. 2013, The Astrophysical Journal, 773, 187, doi: 10.1088/0004-637X/773/2/187
  * Naoz et al. (2022) Naoz, S., Rose, S. C., Michaely, E., et al. 2022, ApJL, 927, L18, doi: 10.3847/2041-8213/ac574b
  * Naoz & Silk (2014) Naoz, S., & Silk, J. 2014, The Astrophysical Journal, 795, 102, doi: 10.1088/0004-637X/795/2/102
  * Naoz et al. (2019) Naoz, S., Silk, J., & Schnittman, J. D. 2019, The Astrophysical Journal, 885, L35, doi: 10.3847/2041-8213/ab4fed
  * Naoz et al. (2020) Naoz, S., Will, C. M., Ramirez-Ruiz, E., et al. 2020, The Astrophysical Journal, 888, L8, doi: 10.3847/2041-8213/ab5e3b
  * Nguyen & Bogdanović (2016) Nguyen, K., & Bogdanović, T. 2016, The Astrophysical Journal, 828, 68, doi: 10.3847/0004-637X/828/2/68
  * O’Leary et al. (2009) O’Leary, R. M., Kocsis, B., & Loeb, A. 2009, MNRAS, 395, 2127, doi: 10.1111/j.1365-2966.2009.14653.x
  * O’Leary & Loeb (2012) O’Leary, R. M., & Loeb, A. 2012, Monthly Notices of the Royal Astronomical Society, 421, 2737, doi: 10.1111/j.1365-2966.2011.20078.x
  * Pesce et al. (2018) Pesce, D. W., Braatz, J. A., Condon, J. J., & Greene, J. E. 2018, The Astrophysical Journal, 863, 149, doi: 10.3847/1538-4357/aad3c2
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=1964PhRv..136.1224P&link_type=ABSTRACT
  * Ramirez-Ruiz & Rosswog (2009) Ramirez-Ruiz, E., & Rosswog, S. 2009, The Astrophysical Journal, 697, L77, doi: 10.1088/0004-637X/697/2/L77
  * Ramsden et al. (2022) Ramsden, P., Lanning, D., Nicholl, M., & McGee, S. L. 2022, Monthly Notices of the Royal Astronomical Society, 515, 1146, doi: 10.1093/mnras/stac1810
  * Rantala et al. (2018) Rantala, A., Johansson, P. H., Naab, T., Thomas, J., & Frigo, M. 2018, The Astrophysical Journal, 864, 113, doi: 10.3847/1538-4357/aada47
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523.  http://adsabs.harvard.edu/abs/1988Natur.333..523R
  * Reines et al. (2020) Reines, A. E., Condon, J. J., Darling, J., & Greene, J. E. 2020, The Astrophysical Journal, 888, 36, doi: 10.3847/1538-4357/ab4999
  * Ricarte et al. (2016) Ricarte, A., Natarajan, P., Dai, L., & Coppi, P. 2016, MNRAS, 458, 1712, doi: 10.1093/mnras/stw355
  * Robertson et al. (2006) Robertson, B., Bullock, J. S., Cox, T. J., et al. 2006, The Astrophysical Journal, 645, 986, doi: 10.1086/504412
  * Rodriguez et al. (2006) Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, The Astrophysical Journal, 646, 49, doi: 10.1086/504825
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, Monthly Notices of the Royal Astronomical Society, 449, 49, doi: 10.1093/mnras/stv264
  * Rose et al. (2022) Rose, S. C., Naoz, S., Sari, R., & Linial, I. 2022, The Astrophysical Journal, 929, L22, doi: 10.3847/2041-8213/ac6426
  * Rose et al. (2023) —. 2023, Stellar Collisions in the Galactic Center: Massive Stars, Collision Remnants, and Missing Red Giants, Tech. rep., doi: 10.48550/arXiv.2304.10569
  * Runnoe et al. (2015) Runnoe, J. C., Eracleous, M., Mathes, G., et al. 2015, The Astrophysical Journal Supplement Series, 221, 7, doi: 10.1088/0067-0049/221/1/7
  * Runnoe et al. (2017) Runnoe, J. C., Eracleous, M., Pennell, A., et al. 2017, Monthly Notices of the Royal Astronomical Society, 468, 1683, doi: 10.1093/mnras/stx452
  * Ryu et al. (2020) Ryu, T., Krolik, J., & Piran, T. 2020, The Astrophysical Journal, 904, 73, doi: 10.3847/1538-4357/abbf4d
  * Schnittman (2015) Schnittman, J. D. 2015, The Astrophysical Journal, 806, 264, doi: 10.1088/0004-637X/806/2/264
  * Schödel et al. (2018) Schödel, R., Gallego-Cano, E., Dong, H., et al. 2018, Astronomy and Astrophysics, 609, A27, doi: 10.1051/0004-6361/201730452
  * Severgnini et al. (2021) Severgnini, P., Braito, V., Cicone, C., et al. 2021, Astronomy and Astrophysics, 646, A153, doi: 10.1051/0004-6361/202039576
  * Shankar et al. (2016) Shankar, F., Bernardi, M., Sheth, R. K., et al. 2016, Monthly Notices of the Royal Astronomical Society, 460, 3119, doi: 10.1093/mnras/stw678
  * Shen et al. (2021) Shen, Y., Chen, Y.-C., Hwang, H.-C., et al. 2021, Nature Astronomy, 5, 569, doi: 10.1038/s41550-021-01323-1
  * Sillanpaa et al. (1988) Sillanpaa, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, The Astrophysical Journal, 325, 628, doi: 10.1086/166033
  * Smith et al. (2010) Smith, K. L., Shields, G. A., Bonning, E. W., et al. 2010, The Astrophysical Journal, 716, 866, doi: 10.1088/0004-637X/716/1/866
  * Springel & Hernquist (2005) Springel, V., & Hernquist, L. 2005, The Astrophysical Journal, 622, L9, doi: 10.1086/429486
  * Stemo et al. (2021) Stemo, A., Comerford, J. M., Barrows, R. S., et al. 2021, The Astrophysical Journal, 923, 36, doi: 10.3847/1538-4357/ac0bbf
  * Stone et al. (2018) Stone, N. C., Generozov, A., Vasiliev, E., & Metzger, B. D. 2018, Monthly Notices of the Royal Astronomical Society, 480, 5060, doi: 10.1093/mnras/sty2045
  * Stone et al. (2018) Stone, N. C., Generozov, A., Vasiliev, E., & Metzger, B. D. 2018, MNRAS, 480, 5060, doi: 10.1093/mnras/sty2045
  * Stone & Metzger (2016) Stone, N. C., & Metzger, B. D. 2016, 25
  * Stone & Metzger (2016) Stone, N. C., & Metzger, B. D. 2016, MNRAS, 455, 859, doi: 10.1093/mnras/stv2281
  * Tory et al. (2022) Tory, M., Grishin, E., & Mandel, I. 2022, Empirical Stability Boundary for Hierarchical Triples, Tech. rep.  https://ui.adsabs.harvard.edu/abs/2022arXiv220814005T
  * Tremaine et al. (2002) Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, The Astrophysical Journal, 574, 740, doi: 10.1086/341002
  * Tsalmantza et al. (2011) Tsalmantza, P., Decarli, R., Dotti, M., & Hogg, D. W. 2011, The Astrophysical Journal, 738, 20, doi: 10.1088/0004-637X/738/1/20
  * van Velzen et al. (2020) van Velzen, S., Holoien, T. W. S., Onori, F., Hung, T., & Arcavi, I. 2020, Space Sci. Rev., 216, 124, doi: 10.1007/s11214-020-00753-z
  * Wang & Merritt (2004) Wang, J., & Merritt, D. 2004, The Astrophysical Journal, 600, 149, doi: 10.1086/379767
  * Wang et al. (2017) Wang, L., Greene, J. E., Ju, W., et al. 2017, The Astrophysical Journal, 834, 129, doi: 10.3847/1538-4357/834/2/129
  * Wegg & Nate Bode (2011) Wegg, C., & Nate Bode, J. 2011, ApJ, 738, L8.  http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2011ApJ...738L...8W&link_type=ABSTRACT
  * Wen et al. (2020) Wen, S., Jonker, P. G., Stone, N. C., Zabludoff, A. I., & Psaltis, D. 2020, The Astrophysical Journal, 897, 80, doi: 10.3847/1538-4357/ab9817
  * Will (2017) Will, C. M. 2017, Physical Review D, 96, 023017, doi: 10.1103/PhysRevD.96.023017
  * Yang et al. (2017) Yang, C., Wang, T., Ferland, G. J., et al. 2017, The Astrophysical Journal, 846, 150, doi: 10.3847/1538-4357/aa8598
  * Zhang et al. (2023) Zhang, E., Naoz, S., & Will, C. M. 2023, A Stability Timescale for Non-Hierarchical Three-Body Systems, Tech. rep., doi: 10.48550/arXiv.2301.08271
  * Zheng et al. (2020) Zheng, X., Lin, D. N. C., & Mao, S. 2020, The Astrophysical Journal, 905, 169, doi: 10.3847/1538-4357/abc8e5