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

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.


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).
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 relax- 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•), 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 bin , and the resulting EKL timescale is the red line labeled t 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/∆h (the relative change in angular momentum 'h'), due to both two-body relaxation and EKL, where h/∆h relx ∼ t relx /P• for two-body processes, and h/∆h EKL ∼ t EKL /P• for EKL. See text ( §3.1) for more details. ation 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 2.1. Fiducial System We consider a system of SMBH binary with masses m 1 and m 2 , and orbital period P bin . Surrounding the primary m 1 is a sphere of compact objects at distance r • and masses m • , where for simplicity we assume the same masses 6 . 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 • ) density profile ρ(r • ) is calibrated by the M − σ relation : where M 0 = 10 8 M , and σ 0 = 200 km sec −1 , are scaling factors. Below we adopt a Bahcall & Wolf (1976)  Each BH (m • ) 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.

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 EKL ∼ 16 30π where P bin and e bin are the period and eccentricity of the SMBH binary respectively, and P • is the period of the stellar mass black hole around m 1 . We show this timescale in Figure 1.

General relativity and Gravitational Waves
The 1st post Newtonian effects induced by m 1 cause m • to precess on a characteristic timescale 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 • 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: where f (e • ) is a function of e • and for all values of e • is between 0.979 and 1.81, (Blaes et al. 2002). We show this timescale for our fiducial system in Figure 1 for e • = 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): where m scat is the mass of the average scatterer, σ is the velocity dispersion of BHs around the SMBH where α is the slope of the density profile. The coulomb logarithm is: For simplicity we adopt m scat ≈ m • . However, if m scat << m • , 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 • = Gm 1 (2/r • − 1/a • ) due to one encounter is:  We also consider a population of stellar mass black holes around m 1 , following a Bahcall & Wolf (1976) profile (i.e., α = 1.75). We normalize the density profile according to the m − σ relation, (see Equation (1)), which result in two-body relaxation timescale of t relx ∼ 3.5 × 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 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 × 10 11 yrs (by assuming scatter masses of 5 × 10 −3 M ). As depicted this system qualitatively follows the EKL (+GR) behaviour.
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 where f • is the true anomaly 7 . Thus, the vector r • in the invariable plane 8 can be considered constant during the encounter. See appendix A for full set of the two-body relaxation equations.

DYNAMICAL EVOLUTION
3.1. Example system and revisiting the time-scale argument 7 Note that we choose the Eccentric anomaly from a uniform distribution and finding the true anomaly from there. 8 Note 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).
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 EKL ∼ 1.4 × 10 4 yr) is shorter than the GR precession timescale (t 1pN ∼ 6 × 10 6 yr).
However, as can be seen in Figure 2, left column, a twobody 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., α = 7/4). The two-body relaxation timescale from Equation (5) is t relx ≈ 3.5 × 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 ≈ 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 relx in this case). However, the combined effect of twobody 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 (e.g., Lithwick & Naoz 2011a Thus, as a • increases, so does , 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 relx . The answers to both of these questions can be understood by examining Equation (8), which suggests that h/∆h| relx ∼ √ t relx , where h is the angular momentum and δ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/∆h| EKL ∼ t EKL (e.g., Naoz et al. 2013a). Thus, effectively, we should compare t relx /P • to t EKL /P • . We show this comparison in Figure 1, bottom panel, where we compare h/∆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/∆h| relx > h/∆h| 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., α = 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 twobody 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 relx ≈ 4.3 × 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/∆h| relx ≈ 700 × h/∆h| 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 (+1pN), as shown by the thick red lines. Note that the apparent drift in ω in this case is due to the 1pN 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 1pN precession does not take place. To guide the eye we have outline the t EKL = t 1pN line for a e • = 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, ω, 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.

Monte Carlo Proof-of-concept
As mentioned, two-body relaxation processes are often neglected when analyzing the EKL like systems. On the other 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 bin = 1 pc and eccentricity of e 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 + twobody 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 crit , thus a system that ended up below the line is marked as a potential GW source, i.e., EMRI.
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 bin = 1 pc and eccentricity of e bin = 0.7. We populate the area of m 1 with 1000 stellar mass BHs, adopting a Bahcall & Wolf (1976), i.e., α = 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π.
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: 3. The BH semi-major axis changed due to two-body relaxation such that > 0.1 (pink points, to the right of the dashed line). This is only possible when the twobody relaxation is turned on.
While it is clear that the systems whose pericenter crossed R sc , are GW sources (i.e., EMRI candidates), it may be less obvious to understand what is the outcome of those with > 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 well 9 . We find that between ≈ 50 − 100% of the BHs (corresponding to pericenter smaller than R crit , to > 0.1) become a GW source.
For peri-centers smaller than R 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., α, 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.).

EMRI RATE ESTIMATION
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., Thus, here we aim to highlight the efficiency of the proposed mechanism by utilizing the M − σ 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 9 Note 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. (2008Chen et al. ( , 2009Chen et al. ( , 2011. Similar arguments were done for systems for which > 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 + twobody relaxation processes may continue to occur for BHs for which > 0.1, until resulting in possible EMRIs (see Appendix B).  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 EKL , constant mass ratio, and that the maximum distance for stellar mass BHs corresponds to = 0.1, see Eq. (10). We estimate the number of stellar mass BHs using the M − σ, relation (see text for more details).
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 × 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×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 place 10 .
To estimate the number of black holes, n BH (≤ r • ), within a distance r max we use the M − σ relation: where M(≤ r max ) = rmax 0 ρ(r )4πr 2 dr and ρ is the density profile form Equation (1). Furthermore, m is the average mass 10 Note 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.
of the stars and f BH is the fraction of BHs from the overall stellar population, where we adopt f BH = 3.2 × 10 −3 (e.g., Aharon & Perets 2016). In our fiducial system r max = 0.07 pc, which corresponds to the = 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) constant 11 and considering the number of BHs up to r max , for = 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 max is then mass dependent and it takes the following form: 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: where f EMRI is the fraction of systems that may become an EMRI rather than a plunged orbit, f EKL is the fraction of systems that have their eccentricity excited to cross R sch and Γ EKL is the rate estimated in the simulation. We estimate the latter by calculating the average accretion time and estimating ±68% of it from our fiducial simulations (i.e., taking 1σ 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 EKL ∼ 0.5 − 1.
In Appendix B, we estimate the fraction of systems that are likely to appear within the LISA band ( f 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 2-bod y relx 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 (14). 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 rmax (sphere of influence), shown as dashed (solid) line. 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. (11), the rate is proportional to the mass of the SMBH primary in the following way: Thus, for the scaling relation chosen in this proof of concept, where α = 7/4, q = 0.01, e 2 = 0.7, and t EKL are constant, the rate is proportional to m 7/8 1 . We show this rate in Figure 5, as the shaded band for the following limits: f EKL × f EMRI = 1 − 0.2, where f EKL corresponds to having 50% (100%) from the total number of available BHs become EMRIs and f 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×10 7 yr. The EKL (+GR)-only scenario may represent a shallow density distribution (α ≈ 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.
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 max [i.e., n BH (≤ r max ), Eq. (11)]. 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: where δ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 85 3072 where for consistency we evaluate this critical value at r max (as well as t 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 max , or the sphere of influence), the rate depends on the SMBH mass weakly: Γ relx ≈ m −1/4 1 , (e.g., Hopman & Alexander 2005). In Figure 5, we show this rate for r max (the sphere of influence), dashed (solid) line.

DISCUSSION
EMRIs are the result of an SMBH that captures a stellarmass 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 twobody 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/∆h| relx needs to be within a couple of orders of magnitude (or close to) h/∆h| EKL . If h/∆h| relx >> h/∆h| EKL , the angular momentum change ∆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/∆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 ∆h/h RR ∼ t res /P • , 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 ∆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 EKLfavorable 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 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 1pN contribution to the inner and outer orbits 12 . 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 EM-RIs, within a few×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 EM-RIs rate is highly uncertain and is sensitive to the number of BHs as a function of time, we used the M − σ 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 bin = 1 pc(a bin = 2.2 pc). Furthermore, the number of BHs inside a sphere at which ≤ 0.1 varies accordingly, for m 1 = 10 7 M (m 1 = 10 8 M ), N BH ≈ 331 (N BH ∼ 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 The associated velocity vector at the plane of the ellipse is: v • = h/a • (− sin f • , e • + cos f • , 0)/(1 − e 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, ω, longitude of ascending nodes, Ω and inclination i. For example, and similarly for the velocity vector, we have: where the subscript "inv" and "ell" refer to the invariable and ellipse coordinate systems, respectively. Given a rotation angle θ, the rotation matrices R z and R x are R z (θ) = cos θ − sin θ 0 sin θ cos θ 0 0 0 1 and R x (θ) = 1 0 0 0 cos θ − sin θ 0 sin θ cos θ . (A4) A two-body encounter can change its velocity by: We model this change as a random walk, applying a single isotropic, instantaneous, kick to the BH velocity once per P • . Each directional component of this 3D kick is drawn from a Gaussian distribution with a zero average and a standard deviation of ∆v j / √ 3, where j is 1, 2, 3 for the three components of the velocity vector. The instantaneous assumption means that r • is kept constant during the kick (see Kalogera 2000).
Thus, post-kick, the new velocity vector (in the invariable plane) is: v •,p = ∆v + v • , where we dropped the subscript "inv" to avoid clutter and the subscript "p" means post-kick. The angular momentum post-kick is h p = r • ×v •,p . Thus, it is straightforward to find the orbital parameters. Specifically, the semi-major axis of the BH post-kick is: •,p Figure 6. An example of the characteristic strain for 150 (chosen randomly) of the runs that reached P• ≤ 10 yr, and a•(1 − e•) < 1 au. We consider the case which includes EKL (+GR) + two-body relaxation + GW.
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.