11affiliationtext: Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095, USA22affiliationtext: Mani L. Bhaumik
Institute for Theoretical Physics, Department of Physics and Astronomy, UCLA,
Los Angeles, CA 90095, USA33affiliationtext: Department of Astronomy and
Astrophysics, University of California,Santa Cruz, CA 95064,
USA44affiliationtext: NASA Goddard Space Flight Center, Greenbelt, MD 20771,
USA55affiliationtext: Maryland Joint Space-Science Institute, College Park, MD
20742, USA

#  The Combined Effects of Two-Body Relaxation Processes and the Eccentric
Kozai-Lidov Mechanism on the EMRI Rate

Smadar Naoz1,2, Sanaea C. Rose1,2, Erez Michaely1,2, Denyz Melchor1,2, Enrico
Ramirez-Ruiz3, Brenna Mockler3, Jeremy D. Schnittman4,5

###### Abstract

Gravitational wave (GW) emissions from extreme-mass-ratio inspirals (EMRIs)
are promising sources for low-frequency GW-detectors. They result from a
compact object, such as a stellar-mass black-hole (BH), captured by a
supermassive black hole (SMBH). Several physical processes have been proposed
to form EMRIs. In particular, weak two-body interactions over a long time
scale (i.e., relaxation processes) have been proposed as a likely mechanism to
drive the BH orbit to high eccentricity. Consequently, it is captured by the
SMBH and becomes an EMRI. Here we demonstrate that EMRIs are naturally formed
in SMBH binaries. Gravitational perturbations from an SMBH companion, known as
the eccentric Kozai-Lidov (EKL) mechanism, combined with relaxation processes,
yield a significantly more enhanced rate than any of these processes operating
alone. Since EKL is sensitive to the orbital configuration, two-body
relaxation can alter the orbital parameters, rendering the system in a more
EKL-favorable regime. As SMBH binaries are expected to be prevalent in the
Universe, this process predicts a substantially high EMRI rate.

##  1\. Introduction

Extreme-mass-ratio inspirals (EMRIs) arise from the capture of a stellar mass
compact object by a supermassive black hole (SMBH). The Gravitational Wave
(GW) emission from such a system is expected to be at the mHz band, thus a
promising signal for the Laser Interferometer Space Antenna (LISA), as well as
other mHz detectors, such as TianQin (e.g., Amaro-Seoane et al., 2017; Robson
et al., 2019; Baker et al., 2019; Mei et al., 2020). Thus, the rate estimation
of EMRIs is of high importance for these detections.

EMRI rate estimation studies often focused on the ”loss cone” mechanism, in
which stellar mass black holes (BHs) undergo weak two-body scatterings and
over time are able to reach high eccentricities (e.g., Hopman & Alexander,
2005; Aharon & Perets, 2016; Amaro-Seoane, 2018; Sari & Fragione, 2019).
Additionally, weak two-body interactions can also lead to mass segregation if
the BH is more massive than the surrounding population of stars (e.g., Hopman
& Alexander, 2006; Alexander & Hopman, 2009; Preto & Amaro-Seoane, 2010;
Amaro-Seoane & Preto, 2011; Chen & Han, 2018). Other physical processes have
also been suggested to contribute to the formation of EMRIs, for example, the
tidal separation of BH binaries by SMBHs was suggested to form a low
eccentricity LISA event (e.g., Miller et al., 2005; Raveh & Perets, 2021, the
latter also include the effects of mass segregation). Furthermore, accretion
disks around SMBHs in active galactic nuclei (AGN) have been suggested to
further increase the EMRI rate (e.g., Pan & Yang, 2021).

Of particular interest here is the formation of EMRIs in SMBH binaries. Thanks
to the hierarchical nature of galaxy formation, and since almost every galaxy
hosts a SMBH at its center, SMBH binaries are expected to be a common
phenomenon (e.g., Di Matteo et al., 2005; Hopkins et al., 2006; Robertson et
al., 2006; Callegari et al., 2009; Li et al., 2020). Observations of AGN
pairs, which are typically few kpc (and more) apart, suggest that these
configurations may lead to the formation of SMBH binaries with parsec to sub-
parsec separations (e.g., Komossa et al., 2003; Bianchi et al., 2008;
Comerford et al., 2009, 2018; Green et al., 2010; Liu et al., 2010; Smith et
al., 2010; Stemo et al., 2020). Moreover, observations of SMBH binaries on
wide orbits, as well as some sub-parsec candidates, seem to support this idea
(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).
Lastly, a combination of theoretical and observational studies suggested that
our own galactic center may also host a companion (albeit a small one) (e.g.,
Hansen & Milosavljević, 2003; Maillard et al., 2004; Gürkan & Rasio, 2005;
Gualandris & Merritt, 2009; Chen & Liu, 2013; Generozov & Madigan, 2020;
Fragione et al., 2020; Zheng et al., 2020; Naoz et al., 2020; GRAVITY
Collaboration et al., 2020).

A SMBH companion gravitationally perturbs the orbit of a stellar-mass BH via
the Eccentric Kozai-Lidov mechanism (EKL, e.g., Kozai, 1962; Lidov, 1962;
Naoz, 2016, see latter for review). These perturbations can result in extreme
eccentricities (e.g., Li et al., 2014b, a; Naoz & Silk, 2014), which can lead
to the formation of EMRIs (e.g., Bode & Wegg, 2014; Haster et al., 2016). A
similar process is often considered for the production of tidal disruption
events (e.g., Chen et al., 2008, 2009, 2011; Chen & Liu, 2013; Li et al.,
2015, Mockler et al. in prep.).

The EKL approach often neglects collective dynamical interaction because these
interactions operate on much longer timescales. In particular, relaxation by
gravitational encounters typically takes place on such long time scales,
compared to other physical processes (see below), and thus often is neglected
when considering EKL processes. Here we show that the combined effect of EKL
and relaxation processes enhances the EMRI formation efficiency more than any
of these processes operating alone. Furthermore, the two-body relaxation
processes overcome the general relativistic precession that suppresses EKL
resonances. We begin by describing the methodology of the system in Section 2.
We then consider an example system and present proof of concept Monte Carlo
results of a fiducial system in Section 3. We show that this mechanism can
potentially result in a much higher EMRI rate in Section 4, and we offer our
discussion in Section 5.

##  2\. Methodology and system set up

Figure 1.—  Top panel: An illustration of the system. Middle panel: An example
of the timescales in the system. We consider a SMBH of mass $m_{1}=10^{7}$ M⊙
with a population of $10$ M⊙ BHs. The period of the BHs around $m_{1}$ is
shown by the dashed line (labeled $P_{\bullet}$), and the associated 1pN
precession is shown as the blue line, labeled $t_{1pN}$, according to Equation
(3). The weak interactions between the BHs results in the two-body relaxation
timescale (see Eq. (5)), shown by the black line. We also consider a SMBH
companion with $m_{2}=10^{9}$ M⊙, at $1$ pc separation (note that
$m_{1}<m_{2}$ in this configuration). The period of the SMBH binary is shown
as the dashed line labeled $P_{\rm bin}$, and the resulting EKL timescale is
the red line labeled $t_{\rm EKL}$, see Equation (2. Bottom panel: Proposed
alternative to quantify the relative importance of the two-body relaxation
processes compared to EKL. We consider $h/\Delta h$ (the relative change in
angular momentum ‘$h$’), due to both two-body relaxation and EKL, where
$h/\Delta h_{\rm relx}\sim\sqrt{t_{\rm relx}/P_{\bullet}}$ for two-body
processes, and $h/\Delta h_{\rm EKL}\sim{t_{\rm EKL}/P_{\bullet}}$ for EKL.
See text (§3.1) for more details.

###  2.1. Fiducial System

We consider a system of SMBH binary with masses $m_{1}$ and $m_{2}$, and
orbital period $P_{\rm bin}$. Surrounding the primary $m_{1}$ is a sphere of
compact objects at distance $r_{\bullet}$ and masses $m_{\bullet}$, where for
simplicity we assume the same masses111Different population may result in
slightly different density profiles, see Aharon & Perets (2016).. Note that
here $m_{1}<m_{2}$. We emphasize that the physical processes described below
are scalable beyond the fiducial system adopted here. In particular, we expect
that two-body relaxation will play a critical role in the EKL process of a
population of stars and a wide range of compact object masses surrounding
SMBHs. The stellar mass BHs ($m_{\bullet}$) density profile
$\rho(r_{\bullet})$ is calibrated by the $M-\sigma$ relation (Tremaine et al.,
2002):

$\rho(r_{\bullet})=\frac{3-\alpha}{2\pi}\frac{m_{1}}{r_{\bullet}^{3}}\left(\frac{G\sqrt{m_{1}M_{0}}}{\sigma_{0}^{2}r_{\bullet}}\right)^{-3+\alpha}\
,$ (1)

where $M_{0}=10^{8}$ M⊙, and $\sigma_{0}=200$ km sec-1, are scaling factors.
Below we adopt a Bahcall & Wolf (1976) profile, i.e., $\alpha=1.75$. Note that
these values have been slightly modified recently (e.g., van den Bosch, 2016;
McConnell & Ma, 2013). However, it does not affect the underlying physical
processes described below, and only may slightly change the relaxation
timescale (see Section 2.4).

Each BH ($m_{\bullet}$) undergoes eccentricity and inclination excitations due
to the far away SMBH companion ($m_{2}$) according to the EKL mechanism.
Additionally, general relativity effects induce precession and can also
circularize and shrink the orbit through gravitational wave (GW) emission.
Finally, collective relaxation interactions with the sea of objects in the
sphere of influence tend to change the angular momentum and energy of the
orbit by an order of themselves over long timescales. Below we specify these
different physical processes and outline the methodology of including them in
our analysis.

###  2.2. Three body secular analysis

We solve the hierarchical three body secular equations up to the octupole-
level of approximation (see for complete set of equations Naoz et al., 2013a).
The timescale of the lowest order of approximation, namely the quadrupole, is
(e.g., Antognini, 2015) estimated by

$t_{\rm
EKL}\sim\frac{16}{30\pi}\frac{m_{1}+m_{\bullet}+m_{2}}{m_{2}}\frac{P_{\rm
bin}^{2}}{P_{\bullet}}(1-e_{\rm bin}^{2})^{3/2}\ ,$ (2)

where $P_{\rm bin}$ and $e_{\rm bin}$ are the period and eccentricity of the
SMBH binary respectively, and $P_{\bullet}$ is the period of the stellar mass
black hole around $m_{1}$. We show this timescale in Figure 1.

###  2.3. General relativity and Gravitational Waves

The 1st post Newtonian effects induced by $m_{1}$ cause $m_{\bullet}$ to
precess on a characteristic timescale

$t_{\rm 1pN}\sim\frac{P_{\bullet}c^{2}a_{\bullet}(1-e_{\bullet}^{2})}{6\pi
G(m_{1}+m_{\bullet})}$ (3)

where $c$ is the speed of light. When this timescale is shorter than the
quadrupole timescale from Equation (2) eccentricity excitations are typically
suppressed (e.g., Ford et al., 2000; Naoz et al., 2013b; Will & Maitra, 2017;
Lim & Rodriguez, 2020). However, when these two timescales are similar, the
precession may excite eccentricities and even re-trigger the EKL behaviour of
extreme eccentricity and inclination flips (Naoz et al., 2013b), by
destabilizing the quadrupole level resonance (see, Hansen & Naoz, 2020). The
timescales of $m_{\bullet}$ for a fiducial system are shown in Figure 1.

We include in our calculations both 1st post Newtonian effects from the
primary $m_{1}$ and the secondary $m_{2}$. As mentioned in Naoz & Silk (2014)
and Li et al. (2015), we choose to focus on the BHs around the less massive
SMBH to minimize the part of the parameter space in which 1st pN precession
suppresses the EKL’s eccentricity excitations. As we highlight below, in the
presence of two-body relaxation this suppression is minimized.

In addition to 1pN precession we also include the shrinking and
circularization of the stellar BH orbit due to gravitational wave emission
following Peters & Mathews (1963). The characteristic timescale to merge an
EMRI is:

$\displaystyle t_{\rm GW}$ $\displaystyle\sim$ $\displaystyle 5.8\times
10^{9}~{}{\rm yr}\left(\frac{m_{1}}{10^{6}~{}{\rm
M}_{\odot}}\right)^{-2}\left(\frac{m_{\bullet}}{10~{}{\rm
M}_{\odot}}\right)^{-1}\left(\frac{a_{\bullet}}{10^{-4}~{}{\rm
pc}}\right)^{4}$ (4) $\displaystyle\times$ $\displaystyle
f(e_{\bullet})(1-e_{\bullet}^{2})^{7/2}\ ,$

where $f(e_{\bullet})$ is a function of $e_{\bullet}$ and for all values of
$e_{\bullet}$ is between $0.979$ and $1.81$, (Blaes et al., 2002). We show
this timescale for our fiducial system in Figure 1 for $e_{\bullet}=0.9$.

###  2.4. Two-body relaxation

Scattering relaxation interactions of a target black hole with the sea of
objects are modeled by considering the two-body relaxation timescale (e.g.,
Binney & Tremaine, 2008):

$t_{\rm relx}=0.34\frac{\sigma^{3}}{G^{2}\rho\langle m_{\rm
scat}\rangle\ln\Lambda}\ ,$ (5)

where $\langle m_{\rm scat}\rangle$ is the mass of the average scatterer,
$\sigma$ is the velocity dispersion of BHs around the SMBH

$\sigma^{2}=\frac{Gm_{1}}{r_{\bullet}(1+\alpha)}\ ,$ (6)

where $\alpha$ is the slope of the density profile. The coulomb logarithm is:

$\Lambda=\frac{r_{\bullet}\sigma^{2}}{2G\langle m_{\rm scat}\rangle}\ .$ (7)

For simplicity we adopt $\langle m_{\rm scat}\rangle\approx m_{\bullet}$.
However, if $\langle m_{\rm scat}\rangle<<m_{\bullet}$, mass segregation may
migrate the BHs inwards.

The relaxation time from Equation (5), is the timescale for a change of energy
of the stellar mass BH around the SMBH $m_{1}$ by an order of its orbital
energy, or a change in angular momentum by an order of its circular angular
momentum. We show the relaxation timescale in Figure 1 (solid black line on
top), which for large part of the parameter space is much larger than the EKL
timescale. As mentioned, this motivated many studies to ignore the
contribution of two-body relaxation when considering EKL effects.

The typical change in the BH’s velocity
$v_{\bullet}=\sqrt{Gm_{1}(2/r_{\bullet}-1/a_{\bullet})}$ due to one encounter
is:

$\Delta v=v_{\bullet}\sqrt{\frac{P_{\bullet}}{t_{\rm relx}}}\ .$ (8)

We model this change as a random walk, applying a single isotropically
oriented kick to the BH velocity once per orbit around the SMBH. Each
directional component of this 3D kick is drawn from a Gaussian distribution
with a zero average and a standard deviation of $\Delta v/\sqrt{3}$ (see
Bradnick et al., 2017, for a similar approach for binaries around a single
SMBH). We assume that the kick is instantaneous at some random phase of the
BH’s orbit

$r_{\bullet}=\frac{a_{\bullet}(1-e_{\bullet}^{2})}{1+e_{\bullet}\cos
f_{\bullet}}\ ,$ (9)

where $f_{\bullet}$ is the true anomaly222Note that we choose the Eccentric
anomaly from a uniform distribution and finding the true anomaly from there..
Thus, the vector $\vec{r}_{\bullet}$ in the invariable plane333Note that the
system evolves due to EKL and thus we need to project the separation vector on
the invariable plane. For similar analysis see (e.g., Lu & Naoz, 2019). can be
considered constant during the encounter. See appendix A for full set of the
two-body relaxation equations.

Figure 2.—  Time Evolution of an example system in the presence of different
physical processes. We show, from top to bottom, a stellar mass black hole
separation around an SMBH, inclination with respect to the outer perturber,
argument of perihapsis, and longitude of ascending nodes. Left side: We
consider a stellar mass black hole ($m_{\bullet}=10$ M⊙) orbiting an SMBH
($m_{1}=10^{7}$ M⊙), at $a_{\bullet}=8000$ au, initially with
$e_{\bullet}=0.02$, $\omega_{\bullet}=45^{\circ}$,
$\Omega_{\bullet}=110^{\circ}$. We also consider a population of stellar mass
black holes around $m_{1}$, following a Bahcall & Wolf (1976) profile (i.e.,
$\alpha=1.75$). We normalize the density profile according to the $m-\sigma$
relation, (see Equation (1)), which result in two-body relaxation timescale of
$t_{\rm relx}\sim 3.5\times 10^{8}$ yrs. We show the resulting orbital
evolution of the stellar black hole in the thick red line. We also introduce a
binary SMBH with mass $m_{2}=10^{9}$ M⊙ set on $1$ pc separation, with
eccentricity of $e_{\rm bin}=0.7$. The evolution that includes both the two-
body relaxation and the EKL from the outer orbit (as well as GR precession on
the inner orbit) is shown in thin blue line. As depicted this system reached
extreme eccentricities induced by a combination of two-body relaxation and EKL
and pushed toward the SMBH, producing a GW source. We also consider the case
of which we ignore the contribution of two-body relaxation processes and
consider only the EKL (+GR) in light grey. This system never reached high
eccentricity to become an EMRI. At the right side we consider the same system,
only this time we arbitrary increased the relaxation timescale to $4.3\times
10^{11}$ yrs (by assuming scatter masses of $5\times 10^{-3}$ M⊙). As depicted
this system qualitatively follows the EKL (+GR) behaviour.

##  3\. Dynamical evolution

###  3.1. Example system and revisiting the time-scale argument

The EKL mechanism tends to excite high eccentricities and inclination.
However, only about $30\%$ of the parameter space in the aforementioned
configuration is available to reach the extreme eccentricities needed to drive
an object into the black hole, and cross its Schwarzschild radius (e.g., Li et
al., 2015; Naoz & Silk, 2014; Naoz et al., 2019). As an example, we consider
in Figure 2 a system whose EKL eccentricity excitations do not result in
values sufficient to cross the SMBH’s Schwarzschild radius (grey lines in both
columns). For this system, the EKL timescale ($t_{\rm EKL}\sim 1.4\times
10^{4}$ yr) is shorter than the GR precession timescale ($t_{\rm 1pN}\sim
6\times 10^{6}$ yr).

However, as can be seen in Figure 2, left column, a two-body relaxation
process combined with EKL results in aggravated EKL eccentricity and
inclination excitations. We note that we include GR precession for the inner
and outer orbit. The former suppresses the EKL eccentricity excitations when
two-body relaxation is not included (gray lines). However, in this example we
do not include GW emission. To avoid clutter, GW is included in the Monte
Carlo analysis below. In this example (Figure 2, left column) we consider a
black hole population with a Bahcall & Wolf (1976) distribution (i.e.,
$\alpha=7/4$). The two-body relaxation timescale from Equation (5) is $t_{\rm
relx}\approx 3.5\times 10^{8}$ yrs, well above the the EKL timescale (see also
Figure 1, for this case it’s about four orders of magnitude larger). By
definition, over the $\approx 1.5$ Myrs run, the relaxation timescale is
insufficient to change the angular momentum by an order of itself (because the
timescale is shorter than $t_{\rm relx}$ in this case). However, the combined
effect of two-body relaxation and EKL results in higher eccentricity and
inclination amplitude modulations.

In fact, the eccentricity excitations were large enough to drive this stellar
mass BH onto the SMBH, thus forming an EMRI. The higher eccentricity values
reached are correlated with the BH semi-major-axis slightly drifting to higher
values, due to two-body relaxation, thus getting closer to the secondary SMBH
($m_{2}$). This process yields a shorter EKL timescale (recall Eq. (2)
dependency on the inner orbital period). Furthermore, as the inner orbit gets
closer to the secondary SMBH, the octupole-level of approximation dominates
more. This behaviour is expressed by the pre-factor of the octupole-level
Hamiltonian $\epsilon$ (e.g., Lithwick & Naoz, 2011a):

$\epsilon=\frac{a_{\bullet}}{a_{\rm bin}}\frac{e_{\rm bin}}{1-e_{\rm
bin}^{2}}\ .$ (10)

Thus, as $a_{\bullet}$ increases, so does $\epsilon$, which excites the
eccentricity of the BH toward larger values (e.g., Li et al., 2014b, a).

The obvious questions from this result are why these diffusion processes
create such a large effect, and will it always happen regardless the value of
$t_{\rm relx}$. The answers to both of these questions can be understood by
examining Equation (8), which suggests that $h/\Delta h|_{\rm
relx}\sim\sqrt{t_{\rm relx}}$, where $h$ is the angular momentum and $\delta
h$ is the change of the angular momentum due to a small kick over the particle
orbit around $m_{1}$. However, the angular momentum changes due to the EKL are
$h/\Delta h|_{\rm EKL}\sim{t_{\rm EKL}}$ (e.g., Naoz et al., 2013a). Thus,
effectively, we should compare $\sqrt{t_{\rm relx}/P_{\bullet}}$ to $t_{\rm
EKL}/P_{\bullet}$. We show this comparison in Figure 1, bottom panel, where we
compare $h/\Delta h$, due to the different processes. Using this picture, it
is clearer that two-body relaxation is relevant to a large part of the
parameter space.

In the example depicted in the left column of Figure 2, even though $h/\Delta
h|_{\rm relx}>h/\Delta h|_{\rm EKL}$, it is only by a factor of $20$, which
yields this cumulative effect (examining the bottom panel in Figure 1, helps
clarify the comparison between the two effects). About an order of magnitude
difference can still lead to a significant cumulative effect. This behavior is
similar to the way that GR precession destabilizes the quadrupole resonance,
even when it’s timescale is much longer than the quadrupole level (e.g., Naoz
et al., 2017; Hansen & Naoz, 2020). We note of course that for this system,
two-body relaxation effects would have eventually change the energy and
angular momentum of the orbit by an order of themselves, regardless of EKL.
However, this does not guaranteed an orbit that will plunge onto $m_{1}$. In
our case we have adopted a Bahcall & Wolf (1976), i.e., $\alpha=7/4$, which
results in zero net flux, thus, the BHs are expected to undergo diffusion, but
not preferentially migrate.

We emphasize that the two-body relaxation effect on the orbital configuration
is indeed small compared to the long-term EKL eccentricity excitation. This is
highlighted in Figure 2 for the two-body relaxation-only case (red lines),
which does not excite the eccentricity to any meaningfully high values during
the simulation run-time. Instead, the BH simply undergoes diffusion in its
energy and angular momentum. However, since the EKL is sensitive to the
orbital configuration, the diffusion in energy and angular momentum due to
two-body relaxation can still contribute to large effects on the BH orbit. If
the small changes in the orbit’s energy and angular momentum can cause a
change of the angular momentum of about $10-15\%$, the effects on EKL are
substantial.

For comparison, we consider the same system in Figure 2 (right column), only
this time we artificially increased the relaxation timescale, for illustration
purposes. In this example $t_{\rm relx}\approx 4.3\times 10^{11}$ yrs, which
is also longer than the lifetime of the system, and the BH simply undergoes
diffusion. As clearly depicted in the Figure, the diffusion in this system is
insignificant and does not trigger larger EKL effects. Furthermore, in this
example, we find that $h/\Delta h|_{\rm relx}\approx 700\times h/\Delta
h|_{\rm EKL}$. Thus, the relaxation effects, according to this comparison,
results in a negligible change. In this panel, we again over-plot the two-body
relaxation-only effect (+$1$pN), as shown by the thick red lines. Note that
the apparent drift in $\omega$ in this case is due to the $1$pN precession, a
similar drift, is depicted in the left column, only modulated by the diffusion
processes.

As depicted in the bottom two panels in Figure 3, the nominal suppression of
eccentricity excitations due to $1$pN precession does not take place. To guide
the eye we have outline the $t_{\rm EKL}=t_{1pN}$ line for a
$e_{\bullet}=2/3$. Indeed, without two-body relaxation processes, eccentricity
excitations are suppressed in the presence of GR precession (e.g., Ford et
al., 2000; Naoz et al., 2013b). However, the small kicks result in a
diffusion, thus allowing the eccentricity excitation to take place over a wide
range of the parameter space.

Lastly, a striking feature of Figure 2 is that in the presence of two-body
relaxation the system moves in and out libration regime, not in-sync with EKL.
The resonant angle, $\omega$, is known to change from libration to circulation
in EKL (e.g., Li et al., 2014a; Hansen & Naoz, 2020). However, as depicted,
the diffusion process changes these processes, even when the two-body
relaxation effects are insignificant. These small kicks allow the (already
chaotic) system to transfer zones.

###  3.2. Monte Carlo Proof-of-concept

As mentioned, two-body relaxation processes are often neglected when analyzing
the EKL like systems. On the other hand, EKL is often neglected when
considering objects sinking onto a SMBH. Here, we qualitatively show the
importance of combining these two processes. We consider the system
highlighted in Figure 1, of $m_{1}=10^{7}$ M⊙ and $m_{2}=10^{9}$ M⊙ with a
binary separation of $a_{\rm bin}=1$ pc and eccentricity of $e_{\rm bin}=0.7$.
We populate the area of $m_{1}$ with $1000$ stellar mass BHs, adopting a
Bahcall & Wolf (1976), i.e., $\alpha=7/4$, profile. We also adopt a thermal
distribution for the stellar mass black holes, and a mutual inclination that
is taken from isotropic distribution (i.e., uniform in $\cos i$). The argument
of perihapsis and longtitue of ascending nodes are taken from a uniform
distribution between $0-2\pi$.

Figure 3.—  Monte-Carlo results. As a proof of concept, we consider a system
composed out of $m_{1}=10^{7}$ M⊙ and $m_{2}=10^{9}$ M⊙ with a binary
separation of $a_{\rm bin}=1$ pc and eccentricity of $e_{\rm bin}=0.7$. We
present three runs, of $1000$ particles each. We consider the following
processes: (top) EKL + GR, (middle) EKL + GR + two-body relaxation, (bottom)
EKL + GR + two-body relaxation + GW. The initial conditions are the same at
each run and are shown in grey in each panel. Red line marks the limit of
crossing $R_{\rm crit}$, thus a system that ended up below the line is marked
as a potential GW source, i.e., EMRI.

In Figure 3 we present the results of $3$ set of simulations of $10^{3}$
particles each, while adopting the following physical processes (top) EKL +
GR, (middle) EKL + GR + two-body relaxation, (bottom) EKL + GR + two-body
relaxation + GW. The light grey point in each panel represent the initial
conditions (which are identical in each panel). We have three stopping
conditions:

  1. 1.

The simulation reaches $10^{9}$ yrs (depicted as black small points). This
result only happen in the EKL + GR run (top panel), where about $69\%$ of the
systems have survived throughout the EKL + GR simulation (this is consistent
with the results for dark matter particles, by Naoz & Silk, 2014).

  2. 2.

The stellar mass BH pericenter crossed a critical distance, which we adopt as
$R_{\rm sch}=8Gm_{1}/c^{2}$, (following Naoz & Silk, 2014; Naoz et al., 2019),
which is inside the inside the Kerr black hole’s inner-most retrograde stable
orbit. These are represented by red points below the solid line. We label them
as “GW sources.” In the EKL + GR, about $31\%$ of all systems crossed the the
critical radius, while $50\%$ ($53\%$) of all systems in the EKL + GR + two-
body relaxation (+GW) run have ended up as GW sources.

  3. 3.

The BH semi-major axis changed due to two-body relaxation such that
$\epsilon>0.1$ (pink points, to the right of the dashed line). This is only
possible when the two-body relaxation is turned on.

While it is clear that the systems whose pericenter crossed $R_{\rm sc}$, are
GW sources (i.e., EMRI candidates), it may be less obvious to understand what
is the outcome of those with $\epsilon>0.1$. We emphasize that this condition
for hierarchy is based on the octupole pre-factor and therefore is somewhat
arbitrary (e.g. Lithwick & Naoz, 2011a). Furthermore, it was suggested in
Bhaskar et al. (2021) that violating this role often results in even higher
eccentricities. Thus, we refer to those systems as possible EMRIs candidate as
well444Note that systems that crossed the Roche limit (or the Hill Sphere) of
the secondary may also be considered as systems that descend toward the SMBH
(either the primary or secondary) following Chen et al. (2008, 2009, 2011).
Similar arguments were done for systems for which $\epsilon>0.1$ (e.g.,
Bhaskar et al., 2021). Furthermore, Zhang et al. in prep. showed that even in
the case of this Roche limit crossing of a tertiary the system may not change
its energy or angular momentum at the order of itself for long timescales. In
other words, the system may still considered “stable” and the eccentricity may
continue to increase via EKL. Thus, the combined effect of EKL + two-body
relaxation processes may continue to occur for BHs for which $\epsilon>0.1$,
until resulting in possible EMRIs (see Appendix B). . We find that between
$\approx 50-100\%$ of the BHs (corresponding to pericenter smaller than
$R_{\rm crit}$, to $\epsilon>0.1$) become a GW source.

For peri-centers smaller than $R_{\rm crit}$, Kerr geometry may cause the BHs
to spend a lot of time on the SMBH’s ergosphere (Schnittman, 2015) where GW
emission may shrink their separations. Furthermore, special relativity effects
should also be taken into account (Yunes et al., 2008; Berry & Gair, 2013).

As can be seen from Figure 3, the combination of EKL with two-body relaxation
allows the system to access a larger part of the parameter space, thus
triggering the EKL mechanism. In general, the number of objects that undergo
high eccentricity excitation depend on the density distribution (e.g., Li et
al., 2015, Mockler et al. in prep. ). Moreover, since the two-body relaxation
timescale is highly sensitive to the density profile (i.e., $\alpha$, see for
example, Figure 1 in Rose et al., 2020) we expect that the efficiency of the
combined system will depend on the underlying density distribution (see
Melchor et al. in prep.).

##  4\. EMRI Rate estimation

Figure 4.—  Scaling relation proof of concept. Top panel, shows the
probability density of the BHs that cross $R_{\rm crit}$ and thus becomes a GW
source, as a function of time, in each of the three simulations from Figure 3.
Note that the statistical difference between the two probability densities
that include the two-body relaxation is negligible. In particular, the two-
sample Kolmogorov-Smirnov test does not rejects the null hypothesis, that both
has the dame distribution, at the $20\%$ significance level. Bottom panel
shows the number of stellar mass BHs that became GW sources as a function of
time, for a range of primary masses between $10^{5}-10^{8}$ M⊙ (from bottom to
top). In generating this estimate we have assumed constant $t_{\rm EKL}$,
constant mass ratio, and that the maximum distance for stellar mass BHs
corresponds to $\epsilon=0.1$, see Eq. (10). We estimate the number of stellar
mass BHs using the $M-\sigma$, relation (see text for more details).

The rate estimation is very sensitive to the steady state number of BHs around
the SMBH. It varies over three orders of magnitude between the various
assumptions for EMRIs formation processes (e.g., Freitag, 2001; Hopman &
Alexander, 2005; Hopman, 2009; Amaro-Seoane & Preto, 2011; Aharon & Perets,
2016; Bar-Or & Alexander, 2016; Babak et al., 2017). Thus, here we aim to
highlight the efficiency of the proposed mechanism by utilizing the $M-\sigma$
relation for the number of BHs. We then compare to similar approaches in the
literature for the two-body relaxation process.

The EKL-only runs compared to the ones with two-body relaxation processes
yield a significantly different flux of GW source formation. This is shown in
the top panel of Figure 4, where a striking feature is the EKL (+GR)-only
result. This feature is consistent with a “burst”-like behavior that depletes
the stellar mass BHs, which could otherwise become GW sources (a similar
behavior was found for TDEs and dark matter participle depletion Li et al.,
2015; Naoz & Silk, 2014). Thus, for a relatively short time ($6\times 10^{5}$
yr, corresponding to the width of the distribution), the rate is high, but on
the timescale it takes to replenish the stellar mass BH population, the rate
is low. Replenishment of BHs can take place via mass segregation, which brings
BHs in from the sphere of influence (e.g., Hopman & Alexander, 2006). The
corresponding timescale at the order of the two-body relaxation timescale up
to a factor of the mass ratio between the BHs and the stars. Another source of
replenishment is star formation, which for our galactic center is estimated to
occur every few$\times 10^{6}$ yr (Lu et al., 2013). Unlike the EKL (+GR)-only
result, the inclusion of two-body relaxation expands the timescales at which
GW sources can form, thus allowing for the replenishment of stellar-mass BHs
to take place555Note that in these cases, during the long timescales the SMBH
binary’s separation is expected to shrink, yielding an enhancement to the EMRI
rate (e.g., Iwasa & Seto, 2016). The inclusion of this effect is beyond the
scope of this paper. .

To estimate the number of black holes, $n_{\rm BH}(\leq r_{\bullet})$, within
a distance $r_{\rm max}$ we use the $M-\sigma$ relation:

$\displaystyle n_{\rm BH}(\leq r_{\rm max})$ $\displaystyle=$ $\displaystyle
f_{\rm BH}\frac{M(\leq r_{\rm max})}{\langle m_{\star}\rangle}$
$\displaystyle=$ $\displaystyle 2f_{\rm BH}\frac{m_{1}}{\langle
m_{\star}\rangle}\left(\frac{G\sqrt{m_{1}M_{0}}}{\sigma_{0}^{2}r_{\bullet}}\right)^{-3+\alpha}\
,$

where $M(\leq r_{\rm max})=\int_{0}^{r_{\rm max}}\rho(r^{\prime})4\pi
r^{\prime 2}dr^{\prime}$ and $\rho$ is the density profile form Equation (1).
Furthermore, $\langle m_{\star}\rangle$ is the average mass of the stars and
$f_{\rm BH}$ is the fraction of BHs from the overall stellar population, where
we adopt $f_{\rm BH}=3.2\times 10^{-3}$ (e.g., Aharon & Perets, 2016). In our
fiducial system $r_{\rm max}=0.07$ pc, which corresponds to the
$\epsilon=0.1$, and the number of BHs within this radius is about $330$.

As highlighted in previous studies, it is straightforward to scale the system
to a wide range of primary masses, for a constant mass ratio, while holding
the quadrupole timescale (Eq. 2) constant666Note that we limit our analysis to
systems for which $t_{\rm EKL}<t_{\rm relx}$, to allow for a the behavior
outlined in Figure 3 to take place. and considering the number of BHs up to
$r_{\rm max}$, for $\epsilon=0.1$, (e.g., Naoz & Silk, 2014; Naoz et al.,
2019). Thus, in Figure 4, bottom panel, we show the number of stellar mass BHs
that are sunk onto the SMBH, for the run that includes all of the
aforementioned physical processes. In this scaling, proof of concept, $r_{\rm
max}$ is then mass dependent and it takes the following form:

$\displaystyle r_{\rm max}$ $\displaystyle=$
$\displaystyle\left(\frac{15}{16}\right)^{2/3}t_{\rm
EKL}^{2/3}q^{-2/3}\left(\frac{\epsilon}{e_{\rm bin}}\right)^{2}(1-e_{\rm
bin}^{2})(Gm_{1})^{1/3}$ (12) $\displaystyle\approx$ $\displaystyle
0.07~{}{\rm pc}\left(\frac{t_{\rm EKL}}{930~{}{\rm
yr}}\right)^{2/3}\left(\frac{q}{0.01}\right)^{-2/3}\left(\frac{\epsilon}{0.1}\right)^{2}$
$\displaystyle\times$ $\displaystyle\left(\frac{e_{\rm
bin}}{0.7}\right)^{2}\left(1-\bigg{[}\frac{e_{\rm
bin}}{0.7}\bigg{]}^{2}\right)\left(\frac{m_{1}}{10^{7}~{}{\rm
M}_{\odot}}\right)\ ,$

where $q=m_{1}/m_{2}$ is the mass ratio. We note that both in Figure 3 and
below we refer to these objects as GW sources, and EMRIs.

The EMRI rate is then estimated by:

$\Gamma\approx\Gamma_{\rm EKL}\times f_{\rm EKL}\times f_{\rm EMRI}\times
n_{\rm BH}(\leq r_{\rm max})\ ,$ (13)

where $f_{\rm EMRI}$ is the fraction of systems that may become an EMRI rather
than a plunged orbit, $f_{\rm EKL}$ is the fraction of systems that have their
eccentricity excited to cross $R_{\rm sch}$ and $\Gamma_{\rm EKL}$ is the rate
estimated in the simulation. We estimate the latter by calculating the average
accretion time and estimating $\pm 68\%$ of it from our fiducial simulations
(i.e., taking $1\sigma$ of the accretion time, estimated from Figure 4) and
normalized to the range of primary masses as described above (see Figure 4
bottom panel). As highlighted in Figure 3, a large fraction of systems sink
onto the SMBH when both EKL and two-body relaxation operate, i.e., $f_{\rm
EKL}\sim 0.5-1$.

In Appendix B, we estimate the fraction of systems that are likely to appear
within the LISA band ($f_{\rm EMRI}$). Roughly speaking one divides between
plunging orbits which may be characterized with a short GW burst and EMRIs
that have many to a few cycles before merging with the SMBH (e.g., Rubbo et
al., 2006; Yunes et al., 2008; Berry & Gair, 2013, for further discussion). In
the former case, special relativity correction may need to be included (e.g.,
Yunes et al., 2008). Additionally, we note that the pN treatment unitized here
may break down around a rotating SMBH, because the stellar mass BHs are
expected to spend a lot of time close to the SMBH’s ergosphere, before
continuing on their original trajectory (Schnittman, 2015). Thus, GW emission
may alter their orbit. Therefore, the distinction between plunging and cycling
orbit represents a larger problem in this field.

Based on the above distinction (see Appendix B for more details), we find that
about $40\%$ of the systems may be defined as EMRIs. In Appendix B we also
present possible SNR of an example system. Note that the fraction of systems
that may end up in the LISA band may depend on the distance of the source.

Using our scaling relation, and the number of BHs from Eq. (4), the rate is
proportional to the mass of the SMBH primary in the following way:

$\displaystyle\Gamma$ $\displaystyle\approx$ $\displaystyle\Gamma_{\rm
EKL}\times f_{\rm EKL}\times f_{\rm EMRI}\times 2f_{\rm BH}\times
m_{1}^{(3+\alpha)/6}$ $\displaystyle\times$
$\displaystyle\left(\frac{\left(e_{\rm bin}/\epsilon\right)^{2}}{1-e_{\rm
bin}^{2}}\frac{G^{2/3}M_{0}^{1/2}}{\sigma_{0}^{2}}\right)^{-3+\alpha}\left(\frac{15\times
t_{\rm EKL}}{16q}\right)^{2(3-\alpha)/3}\ .$

Thus, for the scaling relation chosen in this proof of concept, where
$\alpha=7/4$, $q=0.01$, $e_{2}=0.7$, and $t_{\rm EKL}$ are constant, the rate
is proportional to $m_{1}^{7/8}$. We show this rate in Figure 5, as the shaded
band for the following limits: $f_{\rm EKL}\times f_{\rm EMRI}=1-0.2$, where
$f_{\rm EKL}$ corresponds to having $50\%$ ($100\%$) from the total number of
available BHs become EMRIs and $f_{\rm EMRI}=1-0.4$ (see Appendix B).

We also depict the EKL (+GR) - only case, during burst (thin dashed line), and
the average over replenishment time, taken to be few$\times 10^{7}$ yr. The
EKL (+GR)-only scenario may represent a shallow density distribution
($\alpha\approx 1$) for the BHs, where the two-body relaxation effect is
longer and thus can be neglected. However, the density distribution of BH is
expected to be steep (e.g., Bahcall & Wolf, 1976), and therefore, as
highlighted here two-body relaxation processes cannot be ignored. We thus,
predict the shaded band as the rate from SMBH binary.

Figure 5.—  A comparison of EMRI formation rate for consistent number of BHs.
We consider the case which includes EKL (+GR) + two-body relaxation estimated
rate from Equation (4). We compare to the EKL (+GR) only runs, where we
consider during burst (light dashed line), or over sufficient replenishment
time (dark dashed line). The latter is loosely estimated by assuming star
formation episode and life time of stars to be about $50$ Myr. Finally, we
depict the EMRI rate for the number of BH limited up to $r_{\rm max}$ (sphere
of influence), shown as dashed (solid) line.

For comparison, we examine the EMRI rate due to only two-body relaxation. For
consistent comparison, we only consider the rate due to the “available” BHs up
to $r_{\rm max}$ [i.e., $n_{\rm BH}(\leq r_{\rm max})$, Eq. (4)]. The rate is
proportional to the number of BHs over the two-body relaxation timescale.
However, as highlighted by Hopman & Alexander (2005), the onset of GW
dissipation does not necessarily correspond with the emission of detectable GW
emission. Thus, following Hopman & Alexander (2005, equation 31.) we write the
two-body relaxation EMRI rate as:

$\Gamma_{\rm relx}(\leq r_{\rm max})\approx\frac{n_{\rm BH}}{t_{\rm
relx}(r_{\rm max})\ln{\delta J}}\left(\frac{a_{c}}{r_{\rm
max}}\right)^{9/2-2\alpha}\ ,$ (15)

where $\delta J$ is the ratio of the maximal circular angular momentum at
$a_{c}$, compared to the angular momentum at the loss cone. Finally, $a_{c}$
is the critical semimajor axis at which the angular momentum relaxation time
is equal to the GW emission decay time (Hopman & Alexander, 2005)

$a_{c}=r\left(\frac{85}{3072}\right)^{1/(3-\alpha)}\left[\frac{1}{r}\left(\frac{Gm_{\rm
BH}t_{\rm relx}}{\sqrt{GM_{\bullet}}}\right)^{2/3}\right]^{3/(6-2\alpha)}\ ,$
(16)

where for consistency we evaluate this critical value at $r_{\rm max}$ (as
well as $t_{\rm relx})$, but in the literature this and the rate from Eq. (15)
are evaluated at the sphere of influence. Regardless of the distance we choose
(i.e., either $r_{\rm max}$, or the sphere of influence), the rate depends on
the SMBH mass weakly: $\Gamma_{\rm relx}\approx m_{1}^{-1/4}$, (e.g., Hopman &
Alexander, 2005). In Figure 5, we show this rate for $r_{\rm max}$ (the sphere
of influence), dashed (solid) line.

##  5\. Discussion

EMRIs are the result of an SMBH that captures a stellar-mass compact object,
such as BH. Thus, these are some of the promising GW signals for low-frequency
GW detectors such as LISA. Different channels have been suggested to form
EMRIs. In particular, two-body relaxation has been proposed as one of the
likely physical processes to form EMRIs efficiently. In this process, weak
two-body kicks from the population of stars and compact object that surrounds
the SMBH can change the BH’s orbit over time, driving it into the SMBH. On the
other hand, perturbations from SMBH companions, via the EKL mechanism, can
excite the SMBH to high eccentricities, thereby forming EMRIs. Here we
demonstrated that EMRIs are naturally formed in SMBH binaries with higher
efficiency than either of these processes considered alone.

In the presence of an SMBH companion, the EKL mechanism can excite the BHs
eccentricity to high values. However, the EKL mechanism’s efficiency depends
to some extent on the initial conditions (e.g., Li et al., 2014a). Therefore,
the small kicks due to two-body relaxation do not need to accumulate to change
the angular momentum by order of itself. Instead, they can change the orbital
parameters of the stellar mass BH, such as eccentricity, semi-major axis, and
argument of periapsis, rendering it in a favorable EKL regime. We show an
example of such as system in Figure 2. Even if the two-body relaxation
timescale is orders of magnitude longer than the EKL timescale (see Figure 1),
the small-kicks are effective as long as they result in a change of angular
momentum comparable to that due to EKL. In particular, we suggest that
$h/\Delta h|_{\rm relx}$ needs to be within a couple of orders of magnitude
(or close to) $h/\Delta h|_{\rm EKL}$. If $h/\Delta h|_{\rm relx}>>h/\Delta
h|_{\rm EKL}$, the angular momentum change $\Delta h$ due to the two-body
relaxation can be neglected (see for example Figure 2). In Figure 1 we
highlight the proposed comparison between the two-body relaxation process and
EKL, using $h/\Delta h$ rather than timescales.

In general, other collective processes may also be considered. For example,
resonant relaxations (Rauch & Tremaine, 1996), which arise from orbit-averaged
mass distribution of the objects around the primary, can be added as well
(e.g., Eilon et al., 2009; Kocsis & Tremaine, 2011; Sridhar & Touma, 2016;
Touma et al., 2019). However, scalar and vector resonant relaxation processes
modify the angular momentum $\Delta h/h_{\rm RR}\sim t_{\rm res}/P_{\bullet}$,
thus using their timescales to estimate their contribution may not be as
misleading as the aforementioned timescale analysis of the two-body relaxation
(instead of using $\Delta h/h$). Vector resonant relaxation processes have
been added recently to the EKL context and were shown to drive low-inclination
configurations to a more EKL favorable regime (e.g., Hamers et al., 2018).
However, the latter study concluded that overall, the combined effect is not
very efficient in the context of BH-BH mergers. In contrast, as highlighted
here, two-body relaxation results in populating EKL-favorable regimes very
efficiently.

As a proof of concept, we choose a fiducial system composed of an SMBH binary
($m_{1}=10^{7}$ M⊙, and $m_{1}=10^{9}$ M⊙) on an eccentric orbit $e_{\rm
bin}=0.7$, at $1$ pc separation. We begin by considering the effect of the EKL
mechanism on stellar black holes around $m_{1}$ (Figure 3 top panel). Note
that all runs include the $1$pN contribution to the inner and outer
orbits777Note that we do not include crossing terms (e.g., Naoz et al., 2013b;
Lim & Rodriguez, 2020), because their overall effect should be minimal in this
configuration. . Stellar-mass BHs whose pericenter distance passed a critical
value are considered as EMRIs. As noted in previous studies, the efficacy of
this mechanism is about $30\%$ (e.g., Naoz & Silk, 2014). We then
systematically add two-body relaxation (middle panel in Figure 3) and
gravitational wave emission (bottom panel). As a result, the efficacy
increased to $50-100\%$, meaning nearly all of the stellar mass BHs ended up
descending into the SMBH, thereby possibly forming EMRIs, within a few$\times
10^{8}$ yr, after a single star formation burst, i.e., not including
replenishment. .

To highlight the efficiency of this scenario, we extrapolate the EMRI
formation rate to different SMBHs. Since EMRIs rate is highly uncertain and is
sensitive to the number of BHs as a function of time, we used the $M-\sigma$
relation. Moreover, we rescale our fiducial example by keeping the quadrupole-
level of the EKL approximation constant. This means a constant power law and a
constant mass ratio and the SMBH binary separation varies accordingly, for
example, for $m_{1}=10^{7}$ M⊙ ($m_{1}=10^{8}$ M⊙), $a_{\rm bin}=1$ pc$(a_{\rm
bin}=2.2$ pc). Furthermore, the number of BHs inside a sphere at which
$\epsilon\leq 0.1$ varies accordingly, for $m_{1}=10^{7}$ M⊙ ($m_{1}=10^{8}$
M⊙), $N_{\rm BH}\approx 331$ ($N_{\rm BH}\sim 1979$). We depict the rescaling
in Figure 4.

Even for this simple scaling, it is clear that having the entire population of
BHs, (or even just $50\%$) becoming EMRIs has large implications on the EMRI
rate. We compare the predicted EMRI rate from this scenario to the prediction
from two-body relaxation only in Figure 5. As depicted in this Figure, the
EMRI rate in SMBH binaries is orders of magnitude larger than in isolated
SMBHs. Additionally, the dependency on the SMBH mass is different, offering a
potential way to disentangle between the different scenarios. Furthermore,
because SMBH binaries are expected to be ubiquitous in the Universe, our
results suggest that the EMRI rate may be much higher than nominal
estimations. In particular, post star burst galaxies may be interesting
candidates for enhanced EMRIs formation as they possibly host a SMBH binary.
Moreover, this result suggests that the observed EMRI rate may be used to
constrain the prevalence of SMBH binaries in the Universe.

We thank the referee for useful comments. SN acknowledges the partial support
from NASA ATP 80NSSC20K0505 and thanks Howard and Astrid Preston for their
generous support. SR thanks the Nina Byers Fellowship, the Charles E Young
Fellowship, and the Michael A. Jura Memorial Graduate Award for support, as
well as partial support from NASA ATP 80NSSC20K0505. EM acknowledges the
support of s Howard and Astrid Preston, the Mani L. Bhaumik Institute for
Theoretical Physics, and as well as partial support from NASA ATP
80NSSC20K0505. DM acknowledges the partial support from NSF graduate
fellowship, the Eugene Cota-Robles Fellowship, and the NASA ATP 80NSSC20K0505.
B.M. is grateful for the AAUW American Fellowship, and the UCSC Presidents
Dissertation Fellowship. E.R.-R. and B.M. are grateful for support from the
Packard Foundation, Heising-Simons Foundation, NSF (AST-1615881, AST-1911206
and AST-1852393), Swift (80NSSC21K1409, 80NSSC19K1391) and Chandra
(GO9-20122X).

##  Appendix A A: The post kick orbital parameters

Consider a BH orbiting SMBH. In the plane of the ellipse we can define the
separation vector as ${\bf r}_{\bullet}=r_{\bullet}(\cos f_{\bullet},\sin
f_{\bullet},0)$, where $f_{\bullet}$ is the true anomaly and

$r_{\bullet}=\frac{a_{\bullet}(1-e_{\bullet}^{2})}{1+e_{\bullet}\cos
f_{\bullet}}\ .$ (A1)

The associated velocity vector at the plane of the ellipse is: ${\bf
v}_{\bullet}=h/a_{\bullet}(-\sin f_{\bullet},e_{\bullet}+\cos
f_{\bullet},0)/(1-e_{\bullet}^{2})$. These vectors are projected onto the
invariable plane, where in the case of test-particle EKL is simply the plane
of the outer orbit (e.g., Lithwick & Naoz, 2011b). Thus, we rotate the the
separation and velocity vectors at each time-step given their argument of
perihapsis, $\omega$, longitude of ascending nodes, $\Omega$ and inclination
$i$. For example, and similarly for the velocity vector, we have:

${\bf r}_{\bullet,\textnormal{inv}}=R_{z}(\Omega)R_{x}(i)R_{z}(\omega){\bf
r}_{\bullet,\textnormal{ell}}\ ,$ (A2)

where the subscript “inv” and “ell” refer to the invariable and ellipse
coordinate systems, respectively. Given a rotation angle $\theta$, the
rotation matrices $R_{z}$ and $R_{x}$ are

$R_{z}(\theta)=\left(\begin{array}[]{ccc}\cos\theta&-\sin\theta&0\\\
\sin\theta&\cos\theta&0\\\ 0&0&1\end{array}\right)$ (A3)

and

$R_{x}(\theta)=\left(\begin{array}[]{ccc}1&0&0\\\ 0&\cos\theta&-\sin\theta\\\
0&\sin\theta&\cos\theta\end{array}\right)\ .$ (A4)

A two-body encounter can change its velocity by:

$\Delta v=v_{\bullet}\sqrt{\frac{P_{\bullet}}{t_{relx}}}\ .$ (A5)

We model this change as a random walk, applying a single isotropic,
instantaneous, kick to the BH velocity once per $P_{\bullet}$. Each
directional component of this 3D kick is drawn from a Gaussian distribution
with a zero average and a standard deviation of $\Delta v_{j}/\sqrt{3}$, where
$j$ is $1,2,3$ for the three components of the velocity vector. The
instantaneous assumption means that ${r}_{\bullet}$ is kept constant during
the kick (see Kalogera, 2000).

Thus, post-kick, the new velocity vector (in the invariable plane) is: ${\bf
v}_{\bullet,p}={\bf\Delta v}+{\bf v}_{\bullet}$, where we dropped the
subscript “inv” to avoid clutter and the subscript “p” means post-kick. The
angular momentum post-kick is ${\bf h}_{p}={\bf r}_{\bullet}\times{\bf
v}_{\bullet,p}$. Thus, it is straightforward to find the orbital parameters.
Specifically, the semi-major axis of the BH post-kick is:

$a_{\bullet,p}=\left(\frac{2}{r_{\bullet}}-\frac{v_{\bullet,p}^{2}}{Gm_{1}}\right)^{-1}\
,$ (A6)

the post-kick eccentricity is:

$e_{\bullet,p}=\sqrt{1-\frac{h_{p}^{2}}{Gm_{1}a_{\bullet,p}}}\ .$ (A7)

Because the $z$ axis is defined by the outer orbit, the new inclination is
$\cos i_{p}=h_{p,z}/h_{p}$, where $h_{p,z}$ is the $z$ component of the post-
kick angular momentum.

The post-kick longitude of acsending nodes is:

$\Omega_{p}=\arctan_{2}\left(\pm\frac{h_{p,1}}{h_{p}\sin
i_{p}},\mp\frac{h_{p,2}}{h_{p}\sin i_{p}}\right)\ .$ (A8)

The post-kick true anomaly is:

$f_{p}=\arctan_{2}\left(\frac{a_{\bullet,p}(1-e_{\bullet,p}^{2})}{h_{p}e_{\bullet,p}}\dot{R},\frac{1}{e_{\bullet,p}}\bigg{[}\frac{a_{\bullet,p}(1-e_{\bullet,p}^{2})}{r_{\bullet}}-1\bigg{]}\right)\
,$ (A9)

where $\dot{R}=\pm\sqrt{v^{2}_{\bullet,p}-h_{p}^{2}/r_{\bullet}^{2}}$, where
the sign is defined by the sign of ${\bf r}_{\bullet}\cdot{\bf v}_{\bullet,p}$
(e.g. Murray & Dermott, 2000). The post-kick argument of pericenter is then:

$\displaystyle\omega_{p}$ $\displaystyle=$
$\displaystyle\arctan_{2}\left(\frac{r_{\bullet,3}}{r_{\bullet}\sin
i_{p}},\bigg{[}\frac{r_{\bullet,1}}{r_{\bullet}}+\frac{r_{\bullet,3}\sin\omega_{p}\cos
i_{p}}{r_{\bullet}\sin i_{p}}\bigg{]}\sec\Omega_{p}\right)-f_{p}$ (A10)

##  Appendix B B: Plunging orbits and an Example of Signal to Noise in the
LISA band

We first differentiate between plunging orbits and EMRIs, where the former is
described as a burst associated with their peri-center passage. Our adopted
stopping condition of $R_{\rm sch}=8Gm_{1}/c^{2}$ means that beyond this
threshold the BH trajectory will be modified by Kerr geometry and special
relativity (e.g., Schnittman, 2015; Schnittman et al., 2018; Yunes et al.,
2008; Berry & Gair, 2013). The specific trajectories are beyond the scope of
this study. Nonetheless, in the presence of GW emission, we can roughly
estimate the fraction of systems that are more likely to appear as EMRIs
rather than GW bursts. For that, we first confirmed that all of the systems in
the EKL (+GR) + two-body relaxation indeed reach the Schwarzschild radius by
integrating all the systems below the solid line in the bottom panel of Figure
3.

Second, examining the integration prior to the threshold we found that $\sim
40\%$ of the system reach a configuration for which $P_{\bullet}\leq 10$ yr,
and $a_{\bullet}(1-e_{\bullet})<1$ au. This specific configuration is chosen
such that the characteristic strain will appear in the LISA band, resulting in
mHz signals (see below). Assuming LISA lifetime to be about $10$ yr. We
emphasize that the $40\%$ estimation is rather conservative because, as
mentioned, even the plunged BHs trajectories may spend a long time zooming in
the SMBH’s ergosphere, where GW emission may alter their separation can result
in an EMRI-like signal. Note that even when the BH period is smaller than $10$
years (roughly equivalent to S0-2’s orbital period), two-body relaxation may
still result in small kicks, about $0.0003$ of the BH velocity, according to
Eq. (8). Thus, overall the orbit will not substantially change over the BH
period888Note that we are not taking into account star-BH collisions and tidal
interactions that may result in electromagnetic signatures or larger BHs
(e.g., Metzger et al., 2021; Rose et al., 2021; Kremer et al., 2022). .

To estimate the signal to noise we follow Robson et al. (2018, 2019). The
strain and thus the SNR depend on the the orbital period, the eccentricity,
and the luminosity distance. As a proof of concept we depict in Figure 6 the
characteristic strain for all of the runs that crossed $R_{\rm sch}$ in our
nominal system (i.e., all the point below the line in Figure 3). For this
example, we adopt a luminosity distance of $0.7$ Mpc, and LISA observation
time of $10$ years. We find that $52\%$ of the systems have a SNR $>5$. Out of
these systems $3\%$ have GW dissipation timescale which is shorter than $10$
years, which implies that a more careful analysis of the characteristic strain
should be conducted for them (e.g., Barack & Cutler, 2004). Eccentricity
oscillations due to the EKL signature on the characteristic strain (e.g.,
Hoang et al., 2019; Deme et al., 2020) are unlikely to be detected in this
configuration.

Figure 6.—  An example of the characteristic strain for $150$ (chosen
randomly) of the runs that reached $P_{\bullet}\leq 10$ yr, and
$a_{\bullet}(1-e_{\bullet})<1$ au.  We consider the case which includes EKL
(+GR) + two-body relaxation + GW.

## References

  * Aharon & Perets (2016) Aharon, D., & Perets, H. B. 2016, ApJ, 830, L1, 1609.01715
  * Alexander & Hopman (2009) Alexander, T., & Hopman, C. 2009, ApJ, 697, 1861, 0808.3150
  * Amaro-Seoane (2018) Amaro-Seoane, P. 2018, Living Reviews in Relativity, 21, 4, 1205.5240
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P. et al. 2017, arXiv e-prints, arXiv:1702.00786, 1702.00786
  * Amaro-Seoane & Preto (2011) Amaro-Seoane, P., & Preto, M. 2011, Classical and Quantum Gravity, 28, 094017, 1010.5781
  * Antognini (2015) Antognini, J. M. O. 2015, MNRAS, 452, 3610, 1504.05957
  * Babak et al. (2017) Babak, S. et al. 2017, Phys. Rev. D, 95, 103012, 1703.09722
  * Bahcall & Wolf (1976) Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214
  * Baker et al. (2019) Baker, J. et al. 2019, arXiv e-prints, arXiv:1907.06482, 1907.06482
  * Bar-Or & Alexander (2016) Bar-Or, B., & Alexander, T. 2016, ApJ, 820, 129, 1508.01390
  * Barack & Cutler (2004) Barack, L., & Cutler, C. 2004, Phys. Rev. D, 69, 082005, gr-qc/0310125
  * Batcheldor et al. (2010) Batcheldor, D., Robinson, A., Axon, D. J., Perlman, E. S., & Merritt, D. 2010, ApJ, 717, L6, 1005.2173
  * Berry & Gair (2013) Berry, C. P. L., & Gair, J. R. 2013, MNRAS, 433, 3572, 1306.0774
  * Bhaskar et al. (2021) Bhaskar, H., Li, G., Hadden, S., Payne, M. J., & Holman, M. J. 2021, AJ, 161, 48, 2008.04335
  * Bianchi et al. (2008) Bianchi, S., Chiaberge, M., Piconcelli, E., Guainazzi, M., & Matt, G. 2008, MNRAS, 386, 105, 0802.0825
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Blaes et al. (2002) Blaes, O., Lee, M. H., & Socrates, A. 2002, ApJ, 578, 775, arXiv:astro-ph/0203370
  * Bode & Wegg (2014) Bode, J. N., & Wegg, C. 2014, MNRAS, 438, 573
  * Bogdanović et al. (2009) Bogdanović, T., Eracleous, M., & Sigurdsson, S. 2009, ApJ, 697, 288, 0809.3262
  * Boroson & Lauer (2009) Boroson, T. A., & Lauer, T. R. 2009, Nature, 458, 53, 0901.3779
  * Bradnick et al. (2017) Bradnick, B., Mandel, I., & Levin, Y. 2017, MNRAS, 469, 2042, 1703.05796
  * Callegari et al. (2009) Callegari, S., Mayer, L., Kazantzidis, S., Colpi, M., Governato, F., Quinn, T., & Wadsley, J. 2009, ApJ-Lett, 696, L89, 0811.0615
  * Chen & Han (2018) Chen, X., & Han, W.-B. 2018, Communications Physics, 1, 53, 1801.05780
  * Chen & Liu (2013) Chen, X., & Liu, F. K. 2013, ApJ, 762, 95, 1211.4609
  * Chen et al. (2008) Chen, X., Liu, F. K., & Magorrian, J. 2008, ApJ, 676, 54, 0712.0246
  * Chen et al. (2009) Chen, X., Madau, P., Sesana, A., & Liu, F. K. 2009, ApJ, 697, L149, 0904.4481
  * Chen et al. (2011) Chen, X., Sesana, A., Madau, P., & Liu, F. K. 2011, ApJ, 729, 13, 1012.4466
  * Comerford et al. (2009) Comerford, J. M., Griffith, R. L., Gerke, B. F., Cooper, M. C., Newman, J. A., Davis, M., & Stern, D. 2009, ApJ-Lett, 702, L82, 0906.3517
  * Comerford et al. (2018) Comerford, J. M., Nevin, R., Stemo, A., Müller-Sánchez, F., Barrows, R. S., Cooper, M. C., & Newman, J. A. 2018, ApJ, 867, 66, 1810.11543
  * Deane et al. (2014) Deane, R. P. et al. 2014, Nature, 511, 57, 1406.6365
  * Deme et al. (2020) Deme, B., Hoang, B.-M., Naoz, S., & Kocsis, B. 2020, ApJ, 901, 125, 2005.03677
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604, astro-ph/0502199
  * Dotti et al. (2009) Dotti, M., Montuori, C., Decarli, R., Volonteri, M., Colpi, M., & Haardt, F. 2009, MNRAS, 398, L73, 0809.3446
  * Eilon et al. (2009) Eilon, E., Kupi, G., & Alexander, T. 2009, ApJ, 698, 641, 0807.1430
  * Ford et al. (2000) Ford, E. B., Joshi, K. J., Rasio, F. A., & Zbarsky, B. 2000, ApJ, 528, 336, arXiv:astro-ph/9905347
  * Fragione et al. (2020) Fragione, G., Loeb, A., Kremer, K., & Rasio, F. A. 2020, ApJ, 897, 46, 2002.02975
  * Freitag (2001) Freitag, M. 2001, Classical and Quantum Gravity, 18, 4033, astro-ph/0107193
  * Generozov & Madigan (2020) Generozov, A., & Madigan, A.-M. 2020, ApJ, 896, 137, 2002.10547
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration et al. 2020, A&A, 636, L5, 2004.07187
  * Green et al. (2010) Green, P. J., Myers, A. D., Barkhouse, W. A., Mulchaey, J. S., Bennert, V. N., Cox, T. J., & Aldcroft, T. L. 2010, ApJ, 710, 1578, 1001.1738
  * Gualandris & Merritt (2009) Gualandris, A., & Merritt, D. 2009, ApJ, 705, 361, 0905.4514
  * Gürkan & Rasio (2005) Gürkan, M. A., & Rasio, F. A. 2005, ApJ, 628, 236, astro-ph/0412452
  * Hamers et al. (2018) Hamers, A. S., Bar-Or, B., Petrovich, C., & Antonini, F. 2018, ApJ, 865, 2, 1805.10313
  * Hansen & Milosavljević (2003) Hansen, B. M. S., & Milosavljević, M. 2003, ApJ-Lett, 593, L77, arXiv:astro-ph/0306074
  * Hansen & Naoz (2020) Hansen, B. M. S., & Naoz, S. 2020, MNRAS, 499, 1682, 2011.07103
  * Haster et al. (2016) Haster, C.-J., Antonini, F., Kalogera, V., & Mandel, I. 2016, ApJ, 832, 192, 1606.07097
  * Hoang et al. (2019) Hoang, B.-M., Naoz, S., Kocsis, B., Farr, W. M., & McIver, J. 2019, ApJ, 875, L31, 1903.00134
  * Hopkins et al. (2006) Hopkins, P. F., Hernquist, L., Cox, T. J., Di Matteo, T., Robertson, B., & Springel, V. 2006, ApJS, 163, 1, astro-ph/0506398
  * Hopman (2009) Hopman, C. 2009, Classical and Quantum Gravity, 26, 094028, 0901.1667
  * Hopman & Alexander (2005) Hopman, C., & Alexander, T. 2005, ApJ, 629, 362, astro-ph/0503672
  * Hopman & Alexander (2006) ——. 2006, ApJ, 645, L133, astro-ph/0603324
  * Iwasa & Seto (2016) Iwasa, M., & Seto, N. 2016, Phys. Rev. D, 93, 124024, 1508.05762
  * Kalogera (2000) Kalogera, V. 2000, ApJ, 541, 319, astro-ph/9911417
  * Kocsis & Tremaine (2011) Kocsis, B., & Tremaine, S. 2011, MNRAS, 412, 187, 1006.0001
  * Komossa et al. (2003) Komossa, S., Burwitz, V., Hasinger, G., Predehl, P., Kaastra, J. S., & Ikebe, Y. 2003, ApJ-Lett, 582, L15, astro-ph/0212099
  * Komossa et al. (2008) Komossa, S., Zhou, H., & Lu, H. 2008, ApJ-Lett, 678, L81, 0804.4585
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Kremer et al. (2022) Kremer, K., Lombardi, James C., J., Lu, W., Piro, A. L., & Rasio, F. A. 2022, arXiv e-prints, arXiv:2201.12368, 2201.12368
  * Li et al. (2014a) Li, G., Naoz, S., Holman, M., & Loeb, A. 2014a, ApJ, 791, 86, 1405.0494
  * Li et al. (2014b) Li, G., Naoz, S., Kocsis, B., & Loeb, A. 2014b, ApJ, 785, 116, 1310.6044
  * Li et al. (2015) ——. 2015, MNRAS, 451, 1341, 1502.03825
  * Li et al. (2020) Li, K., Bogdanović, T., & Ballantyne, D. R. 2020, ApJ, 896, 113, 2006.08520
  * Lidov (1962) Lidov, M. L. 1962, planss, 9, 719
  * Lim & Rodriguez (2020) Lim, H., & Rodriguez, C. L. 2020, Phys. Rev. D, 102, 064033, 2001.03654
  * Lithwick & Naoz (2011a) Lithwick, Y., & Naoz, S. 2011a, ApJ, 742, 94, 1106.3329
  * Lithwick & Naoz (2011b) ——. 2011b, ApJ, 742, 94, 1106.3329
  * Liu et al. (2010) Liu, X., Greene, J. E., Shen, Y., & Strauss, M. A. 2010, ApJ-Lett, 715, L30, 1003.3467
  * Lu & Naoz (2019) Lu, C. X., & Naoz, S. 2019, MNRAS, 484, 1506, 1805.06897
  * Lu et al. (2013) Lu, J. R., Do, T., Ghez, A. M., Morris, M. R., Yelda, S., & Matthews, K. 2013, ApJ, 764, 155, 1301.0540
  * Maillard et al. (2004) Maillard, J. P., Paumard, T., Stolovy, S. R., & Rigaut, F. 2004, A&A, 423, 155, arXiv:astro-ph/0404450
  * McConnell & Ma (2013) McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184, 1211.2816
  * Mei et al. (2020) Mei, J. et al. 2020, Progress of Theoretical and Experimental Physics, 2021, https://academic.oup.com/ptep/article-pdf/2021/5/05A107/37953035/ptaa114.pdf, 05A107
  * Metzger et al. (2021) Metzger, B. D., Stone, N. C., & Gilbaum, S. 2021, arXiv e-prints, arXiv:2107.13015, 2107.13015
  * Miller et al. (2005) Miller, M. C., Freitag, M., Hamilton, D. P., & Lauburg, V. M. 2005, ApJ, 631, L117, astro-ph/0507133
  * Murray & Dermott (2000) Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics, ed. Murray, C. D. & Dermott, S. F.
  * Naoz (2016) Naoz, S. 2016, ARA&A, 54, 441, 1601.07175
  * Naoz et al. (2013a) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013a, MNRAS, 431, 2155, 1107.2414
  * Naoz et al. (2013b) Naoz, S., Kocsis, B., Loeb, A., & Yunes, N. 2013b, ApJ, 773, 187, 1206.4316
  * Naoz et al. (2017) Naoz, S., Li, G., Zanardi, M., de Elía, G. C., & Di Sisto, R. P. 2017, AJ, 154, 18, 1701.03795
  * Naoz & Silk (2014) Naoz, S., & Silk, J. 2014, ApJ, 795, 102, 1409.5432
  * Naoz et al. (2019) Naoz, S., Silk, J., & Schnittman, J. D. 2019, ApJ, 885, L35, 1905.03790
  * Naoz et al. (2020) Naoz, S., Will, C. M., Ramirez-Ruiz, E., Hees, A., Ghez, A. M., & Do, T. 2020, ApJ, 888, L8, 1912.04910
  * Pan & Yang (2021) Pan, Z., & Yang, H. 2021, Phys. Rev. D, 103, 103018, 2101.09146
  * Pesce et al. (2018) Pesce, D. W., Braatz, J. A., Condon, J. J., & Greene, J. E. 2018, ApJ, 863, 149, 1807.04598
  * Peters & Mathews (1963) Peters, P. C., & Mathews, J. 1963, Physical Review, 131, 435
  * Preto & Amaro-Seoane (2010) Preto, M., & Amaro-Seoane, P. 2010, ApJ, 708, L42, 0910.3206
  * Rauch & Tremaine (1996) Rauch, K. P., & Tremaine, S. 1996, NA, 1, 149, astro-ph/9603018
  * Raveh & Perets (2021) Raveh, Y., & Perets, H. B. 2021, MNRAS, 501, 5012, 2011.13952
  * Robertson et al. (2006) Robertson, B., Bullock, J. S., Cox, T. J., Di Matteo, T., Hernquist, L., Springel, V., & Yoshida, N. 2006, ApJ, 645, 986, astro-ph/0503369
  * Robson et al. (2018) Robson, T., Cornish, N., & Liu, C. 2018, arXiv e-prints, 1803.01944
  * Robson et al. (2019) Robson, T., Cornish, N. J., & Liu, C. 2019, Classical and Quantum Gravity, 36, 105011, 1803.01944
  * Rodriguez et al. (2006) Rodriguez, C., Taylor, G. B., Zavala, R. T., Peck, A. B., Pollack, L. K., & Romani, R. W. 2006, ApJ, 646, 49, astro-ph/0604042
  * Rose et al. (2020) Rose, S. C., Naoz, S., Gautam, A. K., Ghez, A. M., Do, T., Chu, D., & Becklin, E. 2020, ApJ, 904, 113, 2008.06512
  * Rose et al. (2021) Rose, S. C., Naoz, S., Sari, R., & Linial, I. 2021, arXiv e-prints, arXiv:2201.00022, 2201.00022
  * Rubbo et al. (2006) Rubbo, L. J., Holley-Bockelmann, K., & Finn, L. S. 2006, ApJ, 649, L25
  * Runnoe et al. (2017) Runnoe, J. C. et al. 2017, MNRAS, 468, 1683, 1702.05465
  * Sari & Fragione (2019) Sari, R., & Fragione, G. 2019, ApJ, 885, 24, 1907.03312
  * Schnittman (2015) Schnittman, J. D. 2015, ApJ, 806, 264, 1506.06728
  * Schnittman et al. (2018) Schnittman, J. D., Dal Canton, T., Camp, J., Tsang, D., & Kelly, B. J. 2018, ApJ, 853, 123, 1704.07886
  * Sillanpaa et al. (1988) Sillanpaa, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, ApJ, 325, 628
  * Smith et al. (2010) Smith, K. L., Shields, G. A., Bonning, E. W., McMullen, C. C., Rosario, D. J., & Salviander, S. 2010, ApJ, 716, 866, 0908.1998
  * Sridhar & Touma (2016) Sridhar, S., & Touma, J. R. 2016, MNRAS, 458, 4143, 1509.02401
  * Stemo et al. (2020) Stemo, A., Comerford, J. M., Barrows, R. S., Stern, D., Assef, R. J., Griffith, R. L., & Schechter, A. 2020, arXiv e-prints, arXiv:2011.10051, 2011.10051
  * Touma et al. (2019) Touma, J., Tremaine, S., & Kazandjian, M. 2019, Phys. Rev. Lett., 123, 021103, 1907.01555
  * Tremaine et al. (2002) Tremaine, S. et al. 2002, ApJ, 574, 740, astro-ph/0203468
  * van den Bosch (2016) van den Bosch, R. C. E. 2016, ApJ, 831, 134, 1606.01246
  * Will & Maitra (2017) Will, C. M., & Maitra, M. 2017, Phys. Rev. D, 95, 064003, 1611.06931
  * Yunes et al. (2008) Yunes, N., Sopuerta, C. F., Rubbo, L. J., & Holley-Bockelmann, K. 2008, ApJ, 675, 604, 0704.2612
  * Zheng et al. (2020) Zheng, X., Lin, D. N. C., & Mao, S. 2020, arXiv e-prints, arXiv:2011.04653, 2011.04653