Tidal Disruption Events from the Combined Effects of Two-Body Relaxation and the Eccentric Kozai-Lidov Mechanism

Tidal disruption events (TDEs) take place when a star ventures too close to a supermassive black hole (SMBH) and becomes ruptured. One of the leading proposed physical mechanisms often invoked in the literature involves weak two-body interactions experienced by the population of stars within the host SMBH's sphere of influence, commonly referred to as two-body relaxation. This process can alter the angular momentum of stars at large distances and place them into nearly radial orbits, thus driving them to disruption. On the other hand, gravitational perturbations from an SMBH companion via the eccentric Kozai-Lidov (EKL) mechanism have also been proposed as a promising stellar disruption channel. Here we demonstrate that the combination of EKL and two-body relaxation in SMBH binaries is imperative for building a comprehensive picture of the rates of TDEs. Here we explore how the density profile of the surrounding stellar distribution and the binary orbital parameters of an SMBH companion influence the rate of TDEs. We show that this combined channel naturally produces disruptions at a rate that is consistent with observations and also naturally forms repeated TDEs, where a bound star is partially disrupted on multiple orbits. Recent observations show stars being disrupted in short-period orbits, which is challenging to explain when these mechanisms are considered independently. However, the diffusive effect of two-body relaxation, combined with the secular nature of the eccentricity excitations from EKL, is found to drive stars on short eccentric orbits at a much higher rate.

1. INTRODUCTION Tidal Disruption Events (TDEs) occur when a star passes near a supermassive black hole (SMBH) and gets torn apart by tidal forces (e.g., Hills 1975;Rees 1988;Guillochon & Ramirez-Ruiz 2013).As the star begins to be torn apart, a fraction of it may form an accretion disk which may result in an electromagnetic signature (Rees 1988;Evans & Kochanek 1989;Ulmer 1999;Guillochon et al. 2014).Thus, TDEs are promising signatures for understanding stellar population around SMBHs as well as accretion processes onto SMBHs (e.g., Dai et al. 2021;Gezari 2021;Mockler et al. 2022).The rate of these events seems also to have many repercussions for SMBH and host galaxy demographics (e.g.van Velzen & Farrar 2014;Kochanek 2016;Gezari 2021;Law-Smith et al. 2017;French et al. 2020;Dodd et al. 2021).For instance it has been suggested that the TDE rate can be used to discriminate between SMBH formation scenarios (Stone & Metzger 2016) and help probe the spin distribution of massive SMBHs (Kesden 2012;Law-Smith et al. 2020).Despite these exciting prospects, many challenges exist in effectively estimating TDE rates.
Previous studies calculating rates of TDEs have often focused on the two-body relaxation process.In this process, weak gravitational interactions, kicks from neighboring stars over long periods of time are able to place stars in near radial orbits around an SMBH (Frank & Rees 1976;Rees 1988;Rauch & Tremaine 1996).This channel receives a lot of attention, both analytically and numerically (Magorrian & Tremaine 1999;Wang & Merritt 2004;Brockamp et al. 2011;MacLeod et al. 2012;Stone & Metzger 2016).Notably, these studies estimate the rate at which stars can be positioned in orbit with pericenter distances comparable to or smaller than their corresponding tidal radius (Frank & Rees 1976).The estimated rates, both numerical and analytical, are in agreement that about 10 −5 − 10 −4 stars per year will undergo a TDE in a typical galaxy.
Gravitational perturbations from an SMBH companion can significantly modify the orbits of surrounding stars (e.g., Chen et al. 2008Chen et al. , 2009;;Chen et al. 2011;Chen & Liu 2013).Most notably, the EKL mechanism (Kozai 1962;Lidov 1962;Naoz 2016) has been shown to excite the eccentricities of stars to high values (Li et al. 2014;Li et al. 2014bLi et al. , 2015)).The EKL-channel is expected to result in a burst-like TDE rate, (Mockler et al. 2023), of tens to hundreds of times higher than two-body relaxation alone for ≈ 10 7 years.
Recent observations of repeating tidal disruption events (rTDE) have raised the question of why none of the aforementioned channels, two-body relaxation or EKL, predicted a sizable number of rTDEs.In an rTDE, a star is partially disrupted and may experience multiple disruptive events (MacLeod et al. 2013(MacLeod et al. , 2014;;Campana et al. 2015;Payne et al. 2021Payne et al. , 2022)).A popular example is ASASSN-14ko, which is a periodically flaring transient, every ≈ 115 days (Payne et al. 2021;Liu et al. 2023).How do these stars migrate to such close distances around the SMBH without becoming fully dis-rupted?The two-body relaxation process drives stars on large separations from the SMBH onto nearly radial orbits (e.g., Fragione & Sari 2018;Sari & Fragione 2019); therefore, these stars have large semi-major axes and low angular momentum.The EKL mechanism changes the angular momentum of an orbit and can drive the eccentricity to extreme values (e.g., Naoz et al. 2013a) while keeping the star separation constant.At face value, both of these channels seem to face challenges in explaining rTDEs.Specifically, for repeated TDEs, the star separation needs to be small, 10 −3 to 10 −2 pc for SMBH with masses in the range of 10 6.5−7.86M in order to have an orbital period ranging from 115 days to 30 years (e.g., Payne et al. 2021;Liu et al. 2023;Malyali et al. 2023;Wevers et al. 2023).Another possible formation channel for rTDEs is the widely discussed binary disruption and capture via the Hills mechanism.Interestingly, eccentric stellar disks can yield promising breeding grounds for a high rate of binary disruptions.Having said this, the disruption of highly bound binaries is required in order to explain a bound star with an orbital period ≈ 115 days.
Following Naoz et al. (2022), we propose a novel mechanism combining gravitational perturbations from a faraway SMBH companion with weak two-body interactions from the overall population of stars around the primary SMBH.The combination of these two mechanisms, EKL, and two-body relaxation, is necessary for the formation of rTDEs.While EKL, on its own, conserves the energy of any individual orbit, two-body relaxation describes the diffusive changes of the star's energy and angular momentum.Therefore, diffusive effects can drive a star's orbit to a part of the parameter space that is highly sensitive to EKL, thus triggering eccentricity and inclination excitations, which can more easily result in a rTDE.Here we show that this combined mechanism not only creates a more continuous TDE rate, in contrast to the burst-like EKL channel found in Mockler et al. (2023), but also naturally forms a large relative number of rTDEs.The paper is organized as follows: We describe the basic physical concepts in Section 2.Then, we present our numerical simulations, both for a representative system and then for a large number of runs in Section 3. In Section 4, we describe the TDE rate and formation of rTDEs using this novel mechanism.Finally, we end with a discussion of our findings in Section 5.

BASIC CONCEPTS AND CHARACTERISTIC TIMESCALES
We consider an SMBH with mass m 1 and an SMBH companion m 2 with a semi-major axis a bin and eccen- tricity e bin .In this system, we select m 1 < m 2 , where the primary SMBH was fixed to mass m 1 = 10 7 M , while the secondary SMBH mass was varied between m 2 = 10 8 M and 10 9 M .Surrounding m 1 is a population of stars with mass m ≈ 0.8 M , each, and separation r = a (1 − e 2 )/(1 + e cos f ), where a , e and f are the star's semi-major axis, eccentricity, and true anomaly, respectively.The density profile ρ (r) of the stellar components is calibrated by the M-σ relation (e.g.Tremaine et al. 2002): , where M 0 = 10 8 M and σ 0 = 200 km sec −1 , are scaling factors.We study two nominal density profiles, α = 1 and α = 2, corresponding to a core and cusp stellar distribution, respectively.The maximum separation of the stars is defined by the hierarchical edge set by (e.g., Lithwick & Naoz 2011): as depicted in Figure 1.We solve the secular, hierarchical three-body equations up to the octupole level of approximation following Naoz et al. (2013a).The EKL mechanism can excite the stars' eccentricity and inclination as a time function.These eccentricity excitations can lead to TDEs if the stars' pericenter a (1 − e ), crosses the tidal threshold where R is the radius of the star.Here we assume a population of stars with 0.8 M , which correspond to R = 0.7 R .As mentioned, previous studies showed that the effect of EKL is efficient in creating TDEs for stars distributed around the less massive SMBH (e.g., Li et al. 2015;Mockler et al. 2023).The relevant timescale for these events is (e.g., Antognini 2015) estimated by: This timescale is shown in Figure 2. It has been suggested that the EKL mechanism is less efficient in driving high eccentricity of the stars around Here, we normalize the timescales by the two-body relaxation timescale (Eq.5) and depict the EKL and the 1pN timescales (Eq. 3 and Eq. 4, respectively) as a function of the mass ratio q.To highlight the dependency on the mass ratio, we explicitly denote m disruptor as the one that forms the TDE and m perturber as the far away SMBH companion.To generate this figure, we adopted m perturber = 10 8 M and vary m disruptor between 10 6 M and 10 10 M with e bin = 0.7 and e = 0.9.The binary is set at a bin = 1 pc and the star at a = 0.07 pc.The dashed line denotes q = 1.To the right of the dashed line we have the regime where m disruptor > m perturber , and the 1pN precession suppresses the EKL eccentricity excitations.For m disruptor < m perturber , on the other hand, EKL excitations are dominant.In the presence of two-body relaxation, high eccentricities can be excited to larger values (Figure 2).Throughout the paper, we thus focus on the left regime, which can excite eccentricities more efficiently.
the more massive SMBH because general relativistic (GR) precession can suppress the eccentricity excitation (e.g.Ford et al. 2000;Naoz et al. 2013b;Li et al. 2015;Mockler et al. 2023).The 1st post-Newtonian precession takes place on the following timescale: where c is the speed of light.The relevant timescale is shown in Figure 2. Large eccentricity excitations will take place when the EKL timescale is shorter than the GR precession timescale.We illustrate this in Figure 3, which depicts the EKL and GR precession timescales normalized to the relaxation timescale as a function of the mass ratio.We define the mass ratio as q = m disruptor /m perturber .As shown in the Figure, for m disruptor > m perturber , GR precession dominates the EKL eccentricity excitations, thus suppressing TDE formation.On the other hand, when m disruptor < m perturber , the EKL eccentricity excitation timescale is faster than GR precession.This trend motivates the choice of having m 1 < m 2 in our analysis.We note that the SMBH binary orbital timescale is comparable to the EKL timescale in some parts of the parameter space.At face value, this may suggest that the double average approach adopted here leads to misleading results.However, as was shown by Antonini et al. (2014); Antognini et al. (2014), andLuo et al. (2016), the main difference between our double average approach and an N-body approach (for these setups) is that N-body simulations can result in even higher eccentricity excitations.Thus, we estimate that the eccentricity values achieved below represent a lower limit.
As mentioned before, two-body relaxation processes were proposed as a promising channel to produce TDEs.The typical timescale to change the orbit's angular momentum (h) and the orbital energy by an order of themselves is estimated by (e.g., Binney & Tremaine 2008): where m scat is the mass of the average star that acts as a scatterer, σ is the velocity dispersion of stars around the SMBH where α is the slope of the stellar density profile.Lastly, the coulomb logarithm is: We adopt m scat = m = 0.8 M .We show this timescale (Eq.( 5)), in Figure 2, for different density profiles α, between 1 and 2. The weak two-body relaxation processes are often neglected in the literature because the relevant timescale is larger compared to the EKL timescale (see Figure 2).However, as shown recently, comparing the timescales can be misleading, and instead, the change applied to the angular momentum should be compared (Naoz et al. 2022).As depicted in Figure 2, bottom panel, when comparing the h/∆h, in a large part of the parameter space, the two-body relaxation for a cusp-like profile, α = 2, yields changes to the angular momentum h/∆h comparable, or even shorter than the EKL's.
Following Naoz et al. (2022), we model the change to the stellar orbit's velocity v = Gm 1 (2/r − 1/a ) due to one encounter as a random walk with isotropically oriented kicks to the stellar velocity once per orbit around the SMBH.The kick is assumed to be instantaneous at some random phase of the orbit about the An example of a repeated TDE.The time evolution of a supermassive black hole binary is shown.The primary black hole is of mass m1 = 10 7 M with a secondary of m2 = 10 9 M .There is a stellar mass of m = 0.8 M that is gravitationally bound to the primary at a separation of 4.7 × 10 3 au.The top panel compares the changes to the semi-major axis and pericenter radius of the inner binary when the system undergoes only EKL (+GR) (yellow) and EKL (+GR) with two-body relaxation (blue).The tidal radius is shown in red, where in crossing this line, we consider the system to have become a TDE.The semimajor axis is shown to drop considerably in the EKL and two-body relaxation simulation, indicating that this system could potentially be a repeated TDE.
SMBH.Each component of this three-dimensional kick is drawn from a Gaussian distribution with a zero average and a standard deviation of ∆v/ √ 3, where (see also Bradnick et al. 2017).See Naoz et al. (2022) for the complete set of equations.

Dynamical Evolution of a Representative Example
In Figure 4, we consider our fiducial model of m 1 = 10 7 M , and m 2 = 10 9 M , set on a bin = 4 × 10 5 au, and e bin = 0.5.The star is initialized with a = 4711 au, e = 0.53, Ω = 168 • and ω = 38 • , where Ω and ω are the longitude of ascending nodes and the argument of periapsis, respectively.For this example, we depict two cases, one which only includes EKL + GR, yellow, and EKL + GR + two-body relaxation processes, blue.We note that while we limit the EKL + GR presentation to 15 Myr, to provide a comprehensive comparison, the system never had its eccentricity excited to cause a TDE.However, in the EKL + GR + two-body relaxation run, the star's pericenter crossed the tidal radius signifying a TDE.
The final separation of m 1 and m decreases to 877 au, crossing the tidal radius.At around 10 7 years the system begins eccentricity oscillations which change the inclination and semi-major axis of the primary's orbit with the stellar mass.By 15 Myr, the system has reached an eccentricity of 0.998, resulting in a nearly 90 • rotation in the inclination.The stellar mass can potentially become a repeated tidal disruption event by m 1 and continue to orbit the SMBH, where this rTDE would be seen every P = 9 yr.

Monte-Carlo simulations
We ran a total of 12,000 runs for 12 realizations varying the SMBH binary eccentricity, mass ratio, and the power law of the stellar density profile.The system was initialized with m 1 = 10 7 M and m = 0.8 M .
We chose two representative masses of the companion of m 2 = 10 8 M and 10 9 M .We explore three outer orbit eccentricities of e bin = 0.3, 0.5, 0.7.The semi-major axis of the companion was set to be half the sphere of influence of the primary for all runs.
For the stellar population, we chose two representative density profiles, α = 1 and 2, following Equation 1. Furthermore, we adopted a thermal eccentricity distribution for the stars while the argument of periapsis and longitude of ascending nodes were chosen from a uniform distribution between 0 − 2π.The mutual inclination was chosen from an isotropic distribution, i.e., uniform in cos i. See Table 1 for the list of parameters chosen for the 12 realizations.
Figure 5 shows a representative example of the runs, where the grey points mark the initial conditions.In total, we see four distinct outcomes (the color code in Figure 5, as well as 12 and 13 from the Appendix, is mentioned below): 1.During the evolution the star's pericenter crossed R T , see Eq. ( 2).We mark these systems as stars that become tidally disrupted.These systems are marked in green points below the a (1 − e ) = R T line.The fraction of systems that we mark TDEs is shown as the upper limit in the 6th column (f TDE ) of Table 1.The lower limit of the TDE fraction includes only TDEs that are within the Hill radius of the primary SMBH.
2. Systems that have crossed the R T threshold and have periods shorter or equal to 30 yrs we mark as repeated TDEs (f rTDE ).These systems are highlighted as green stars.
3. System evolved up to 1 Gyr, without crossing the R T , we consider as survived.These systems are marked in dark blue.
4. Finally, in some cases, the system violated the hierarchical condition, Eq. ( 1).These systems are marked as pink and reside to the right of the vertical hierarchical line.The fate of these systems is somewhat unclear, but as discussed in Naoz et al. (2022), we expect these stars to have their eccentricity excited to large values (e.g., Bhaskar et al. 2021), and thus, perhaps all of them will end as TDEs.We consider this case as an upper limit for the TDE fraction (see Table 1, 7th column (f TDE + f non−hierarchical )).
The evolution of systems with a stellar density distribution of α = 1 (leftmost panels of Figure 5) showcase stars ending in the three main outcomes: stars surviving, stars ending beyond the hierarchical edge, and stars becoming TDEs.Systems with α = 2 (rightmost panels of Figure 5) experience most systems having either the stellar component leaving the hierarchical edge or stellar components becoming tidally disrupted.
Not surprisingly, fewer systems survive when α = 2, since this density profile results in a shorter timescale, i.e., larger change to the angular momentum ∆h/h, as shown in Figure 2. See Naoz et al. (2022) for a detailed discussion.For the α = 1 case, more systems survive since this density profile yields a longer timescale, which results in a smaller change to the angular momentum, as depicted in Figure 2.However, compared to the EKLonly scenario (Mockler et al. (2023) in prep.results), we still find the inclusion of two-body relaxation enhances the rate of systems that result in TDEs.

TDE rate
A key factor in estimating the TDE rate is the number of stars in the sphere inside the hierarchical edge N (r ≤ r max ).This number is sensitive to the density profile as well as to the eccentricity of the SMBH binary.To estimate the number of stars within this range, N (≤ r max ), we use the M − σ relation: , and r max = a e bin (1 − e 2 bin ) , where we take = 0.1.Given this number of stars and the fraction of systems that become TDEs for all simulation runs (see Table 1), we can estimate the number of stars that become TDEs, N TDE .Figure 6, depicts N TDE as a function of N (≤ r max ), for the different density profile values (see arrows), eccentricity, and mass ratios.As shown in Figure 6, the cusp profile, α = 2, which has a larger constriction of stars closer to the SMBH, compared to α = 1, results in a larger number of stars.
We calculate the rate of TDEs for the different values of the power law of the stellar distribution (α), the outer orbit eccentricity (e bin ), and the mass ratio (q), as outlined in Table 1.The calculated TDE rates from our simulations are then compared to the average observed TDE rate as taken from (van Velzen et al. 2020) and the observed post-starburst (PSB) TDE rate from (French et al. 2020), shaded regions in Figure 9.The calculated rates span from a minimum which is chosen by including all the hierarchical systems that end up as TDEs (indicated as green dots and green stars below the solid critical line in Figure 5).The maximum rate additionally includes all the non-hierarchical systems (indicated as the light pink dots to the right of the dashed lines in Figure 5).These non-hierarchical systems are likely to be excited to high eccentricities (e.g., Bhaskar et al.Comparing the number of stars to the number of tidally disrupted events.Here, the eccentricity of the outer binary is varied on values e bin = 0.3, 0.5, 0.7.Mass ratios of q = 10 and q = 100 are also compared and are indicated by circles and "×'s," respectively.Pairs of e bin and q are indicated by different colors, as seen in the legend.The stellar density profile is changed, and a dashed line shows α = 1 values to the left and α = 2 to the right.A linear trend is seen in the number of TDEs as e bin decreases and α increases.

2021
) and may eventually become TDEs (similarly to the maximum rate of EMRIs Naoz et al. 2022).
We point out that some of the stars reside beyond the Hill radius of the SMBH binaries.However, as was shown by (e.g., Zhang et al. 2023), this configuration does not produce an instantaneous destabilization of the star's orbit.Rather the system can undergo eccentricity excitations before changing the energy of the stellar orbit (i.e., the semi-major axis).We find that in most of the systems, the destabilization timescale (following Zhang et al. 2023) is larger than the time to become a TDE.We, thus, adopt the lower limit of the rate to be the TDEs within the hierarchical limit (set by ≤ 0.1, see Eq. ( 1)).
Consider first the comparison between the calculated rates as a function of α.As shown above, α is one of the key parameters that affect the combined effect of EKL with 2-body relaxation (see Figure 2).Thus, in Figure 9, we compare the results between α = 1 and α = 2, left and right column, respectively, for the simulations with e bin = 0.3, 0.5, and 0.7 (see labels), and for q = 10 and q = 100, top and bottom panels, respectively.As depicted, the α = 2 case yields a slightly higher rate, with somewhat less dependency on the eccentricity of the SMBH binary eccentricity e bin .The α = 2 case represents a cusp density profile for which the two-body Comparing the number of repeating tidally disrupted events to the total number of tidal disruption events.The eccentricity of the outer binary is varied on values e bin = 0.3, 0.5, 0.7.Mass ratios of q = 10 and q = 100 are represented as circles and "×'s," respectively.The colors represent the pairs of e bin and α.Here, unlike in Fig. 6, only α = 2 is shown as only this cuspy profile produced rTDEs as seen in Fig. 5, 12, 13.
relaxation angular momentum changes are significant over a larger part of the parameter space.Thus, this behavior assists the EKL eccentricity excitations nearly independent of the binary eccentricity.
However, several studies suggested that the stellar distribution in our own Milky Way's galactic center is corelike rather than cusp, with α closer to unity (e.g., Lu & Naoz 2019;Schödel et al. 2017Schödel et al. , 2018Schödel et al. , 2020;;Gallego-Cano et al. 2020).Thus, if our galactic center is representative of galaxy nuclei's stellar distribution, the left column of Figure 9 may be a more representative scenario.In this case, the weak kicks due to two-body relaxation are important but do not wash out the EKL sensitivity to the outer orbit eccentricity.Thus, as depicted in this Figure, lower SMBH eccentricity yields a higher TDE rate (simply more stars).We see that the e bin = 0.5 TDE rate is consistent with the post-star burst TDE observed rate.Thus this may highlight a preferred SMBH formation channel (e.g., Dotti et al. 2012).
We remind the reader that the SMBH binary was set at half the distance to the sphere of influence.The latter dictates the number of stars, within that sphere, via the M − σ relation.Thus, setting the number of stars that can potentially undergo TDEs.Therefore the dependency on the eccentricity is degenerate with the dependency of the SMBH binary separation, which is beyond the scope of this paper.

Formation of repeated TDEs
rTDE observations present many intriguing puzzles.For example, how do these stars get only partially disrupted?This question was addressed in earlier studies (Guillochon & Ramirez-Ruiz 2013;Coughlin & Nixon 2015) and is still under investigation.Here we focus on the question of how these stars get to such short separations from the SMBH without having their eccentricity excited to high values at larger distances from the SMBH.
The combined effect of EKL with two-body relaxation slightly changes the star's SMA due to kicks.Suppose the star migrates to a regime where EKL is more efficient (for example, a higher or an inclination closer to 90 • ).In that case, the star will undergo large eccentricity excitations due to EKL.This behavior is depicted in Figure 4, where the star undergoes faster EKL oscillations, and its eccentricity is excited to higher values as its SMA slightly increases.Furthermore, since the kicks are proportional to the star's orbital velocity (Eq.( 8)), the kick can be larger for a star initially closer (larger velocity).The example system shown in Figure 4 highlights this behavior, and we outline this signature below.
We arbitrarily choose a period of 30 yrs to mark the rTDEs candidates.As for the non-repeating TDEs, the number of expected rTDEs is proportional to the number of stars in the sphere of influence (see Figure 7).As depicted in Figure 5, the initial cuspier density distribution (i.e., α = 2) results in a larger fraction of rTDEs.Thus, these results suggest that detections of rTDEs may be used as an indicator for the stars underlying density profile.
Figure 8 shows the expected rTDE period distribution adopting the M-σ normalization explained above.As shown, this mechanism results in a period distribution consistent with known observations of potential repeated TDEs with periods ranging hundreds of days up to 30 years, (e.g.Payne et al. 2022;Wevers et al. 2023;Liu et al. 2023;Malyali et al. 2023).Note that our simulations do not suggest a preferred initial semimajor axis that rTDEs originate from in the star cluster based on this channel.It was suggested that the Hill's mechanism may contribute to the formation of rTDEs as recently proposed by Cufari et al. (2022).However, one can relate the period of the captured star from the Hill's mechanism as a function of the separation of the inner binary (e.g., Hills 1975;Yu & Tremaine 2003).A binary undergoes many weak gravitational interactions with neighboring stars during its lifetime.These encounters tend to widen and unbind the binary (e.g., Binney & Tremaine 2008).Thus, in order for a binary to remain bound before the Hill's mechanism takes place, its semi-major axis should be smaller than a critical value and larger than at least two times the radius of the binary members (Rose et al. 2020).Taking the minimum separation between a binary (i.e., assuming an extremely hard binary that does not evolve), we can find the corresponding minimum period for a repeated TDE.We show this period in Figure 8, suggesting the long period J1331 could have originated from the Hills mechanism, while the other potential rTDEs seem to be consistent with the channel described here.

Comparison with TDE without two-body relaxation
Combined with EKL, two-body relaxation efficiently drives the stars onto the SMBH.As discussed above, the efficiency depends on the number of stars within the hierarchical limit (thus the SMBH binary separation of their eccentricity) and the density profile.The latter is largely unknown within the inner parsec of observed galaxies.Currently, some of the best estimations can reach as close as 30 pc to a few kiloparsecs (e.g., French et al. 2020), finding a cusp-like distribution.On the other hand, the stellar distribution of our own galactic center seems to follow a core-like distribution (e.g., Genzel et al. 2003;Gallego-Cano et al. 2018) In Figure 10, we compare the combined effect with an EKL-only approach adopted from Mockler et al. (2023).The Figure shows the time-dependent TDE rate for two examples of m 1 = 10 7 M , m 2 = 10 8 M , m = 0.8 M , e bin = 0.5, a bin ≈ 2 pc with the violet shade indicating the system with α = 1 and the cyan shade α = 2.We note that in the case of EKL-only, there is no system that is pushed beyond the hierarchical limit (systems beyond the vertical line in Figure 5).Thus, the maximum value calculated differs between the two runs.However, as can be seen, the possible additional nonhierarchical TDEs do not yield a significant difference, especially for the shallow, α = 1 profile.It is worth noting that the combined effect results in a longer, extended time-dependent rate, perhaps allowing star formation to occur and replenish the TDEs.
As expected, the EKL-only channel produces a burstlike event, depicted as a shorter rate in Figure 10.This behavior was noted in Naoz & Fabrycky (2014) and Naoz et al. (2022).In particular, we suggest that systems with a core-like density distribution, thus a longer two-body relaxation timescale, may undergo a burst-like rate.In other words, these systems may have already formed their TDEs and, without newly formed stars, are less likely to produce a TDE at a high rate.The rates suggested here (and in Mockler et al. 2023) are much higher than previously considered when considering single SMBHs (e.g., Stone & Metzger 2016).Given the right conditions, in terms of density profile and the SMBH binary separation and eccentricity, Figures 9 and 6 suggest that a galaxy may have multiple TDEs.Similarly, for the right conditions, a galaxy may have repeated TDEs (as suggested by Figures 5 and 7).Beyond TDEs, it was recently suggested that the similar physical processes described here could yield much higher extreme mass ratio inspiral (EMRIs) rates (orders of magnitude than estimated before Naoz et al. 2022).Thus, by combing the two findings, we speculate that a combined TDE and EMRI event may have a non-negligible rate.

DISCUSSION
We demonstrate that the combined physical processes of two-body relaxation and EKL in SMBH binaries significantly affect the dynamics of TDEs.Two-body relaxation has been proposed as one of the most promising physical processes to form TDEs efficiently.In this pro-cess, weak two-body kicks from the population of stars surrounding the SMBH can change the star's orbit over time, plunging it into the SMBH from large distances.Perturbations from an SMBH companion via the EKL mechanism can also excite the star to high eccentricities, providing another channel for forming TDEs.
Here we demonstrated that the two body relaxation, combined with the EKL-induced eccentricity, plays a crucial role in forming TDEs, and rTDEs.In this case, stars experience high eccentricities due to the two-body relaxation process, which drives the stars into a more EKL-sensitive regime as demonstrated in Figure 2.This combination not only leads to more TDEs than EKL on its own, but it also naturally produces repeating TDEs.Additionally, two-body relaxation leads to more stars scattering beyond the hierarchical radius (Figure 5, where we expect many to be disrupted (e.g., Grishin et al. 2018;Bhaskar et al. 2021).
In addition to combining the aforementioned two effects, we have also explored the effect of the star's density profile by comparing two extreme cases.One of a The rate for all systems is calculated and compared to both the average observed rate and the observed PSB rate.The left columns show systems with the density distribution power-law α = 1 while the right column panels show that of α = 2.The top panels have a mass ratio of q = 10, and the bottom panels q = 100.As seen, α = 2 systems produce TDEs at an enhanced rate with little dependence on eccentricity.
shallow profile, corresponding to α = 1, and the second of a cusp profile of α = 2 (see Figure 5).Observations of the galactic center suggest that the stellar distribution in the inner 0.1 pc may be shallow, i.e., closer to α = 1 (e.g., Genzel et al. 2003;Schödel et al. 2017Schödel et al. , 2018Schödel et al. , 2020)).While it may be that other galactic nuclei have similar conditions to our Galactic center, TDEs have been preferentially found in post-starburst galaxies, which may have different properties (e.g., Yang et al. 2008;French et al. 2016;Law-Smith et al. 2017;Dodd et al. 2021); thus, we also considered a cusp profile (i.e., α = 2).A shallower density profile has two main consequences.First, it creates a longer two-body relaxation timescale, thus, the angular momentum changes via the 2-body scattering are smaller compared to the EKL (see bottom panel of Figure 2).Second, it effectively changes the number of stars inside the hierarchical limit, which has a direct consequence on the number of stars resulting in TDEs and rTDEs (as highlighted in Figures 6 and  7).
We also investigated the effect of the SMBH binary's mass ratio and eccentricity.We show representative results of our Monte Carlo simulations in Figure 5 (see Figures 12 and 13 in the Appendix for the configurations involving e bin = 0.3 and e bin = 0.7).As seen in Figure 9, the rate is nearly independent of the binary eccentricity.However, e bin = 0.3 tended toward a higher and more extended rate compared to the higher eccen- The shaded rates are for systems that underwent both twobody relaxation and EKL (+GR), while the hatched greylike rates are for those with only EKL (+GR) processes as adapted from Mockler et al. (2023).The adapted rates showcase high peaks on shorter timescales, while rates in this work have lower peaks and are extended further in time.The minimum for all rates is calculated only for our conservative lower limit for TDEs, which are constrained within the Hill radius.The maximum for the EKL-only rates includes the systems that are beyond the Hill and the Roche limits.This limit corresponds to the solid line in the Figure .The maximum value for the rates for the EKL (+GR) + 2 body Relaxation includes the systems beyond the hierarchical limit as well.
tricities.This is a direct result of the fact that lower eccentricity SMBH binary configuration allows for more stars to reside in the hierarchical limit (see for reference Figure 6).Thus, as shown in 9, the rate is largely insensitive to the mass ratio and the orbital eccentricity of the SMBH binary1 .
We note that we have neglected possible stellar collisions or star-compact object collisions, which may have significant consequences on the masses of stars in galactic nuclei (e.g., Rose et al. 2021Rose et al. , 2023)).Additionally, a substantial fraction of stars is expected to be in binaries (see Raghavan & Steprāns (2010), even in the galactic center Naoz et al. (2018)).Collisions, two-body relaxation, and weak perturbations that may unbind the binary stars may also affect the orbital configuration of these binaries (e.g., Rose et al. 2020).The EKL mechanism can also merge binary stars together (e.g., Stephan et al. 2016;Stephan et al. 2019).In general, we expect that binaries that undergo the combined effect of twobody relaxation and EKL from an SMBH companion will yield a possible Hills (e.g., Hills 1975) process or, in some cases, even double TDE (e.g., Mandel & Levin 2015).The details of such a process are beyond the scope of this paper.
The main contributor to the TDE rate is the underlying stellar distribution and, of course, the existence of the SMBH binary.Therefore comparing to observations, this may help constrain the stellar cusp-or corelike distribution and the frequency of SMBH binaries.Note that increasing the cusp of the density profile yields a shorter two-body relaxation timescale compared to a core-like profile.While this was noted before for a single SMBH system, Stone et al. (2018), here we report an opposite dependency as a function of the SMBH (disruptor) mass and thus the number of stars.In other words, since the two-body relaxation timescale is increasing as a function of the SMBH mass, a single SMBH system predicts a decreasing rate as the SMBH mass.In the case of the combined two-body relaxation with EKL (for a given mass ratio), we find the opposite trend, i.e., an increasing rate with the disruptor mass (see Figure 6).A similar result was obtained for the EMRIs rate as a function of the SMBH mass; see figure 5 in Naoz et al. (2022).
As shown by Mockler et al. (2023), the EKL mechanism leads to a unique signature of TDEs on the smaller SMBH companion.Combining the two-body relaxation processes and the EKL mechanism leads to higher eccentricity excitations, yielding an extended time-dependent rate (see Figures 9 and 10).Additionally, we demonstrated that combining EKL and two-body relaxation can be used as a signature of the underlying cusp of the less massive SMBH in an SMBH binary configuration.Further, we suggested that this combined dynamical effect can naturally give rise to rTDEs.Lastly, these results can be used to constrain the SMBH binary fraction.In other words, we suggest that since these predictions only rely on one main assumption, i.e., the existence of binary SMBHs, they can be used to constrain SMBH binary fractions.D.M. acknowledges the partial support from NSF graduate fellowship DGE-2034835, the Eugene V. Cota-Robles Fellowship, and the NASA ATP 80NSSC20K0505.B.M. is grateful for support from the U.C. Chancellor's Postdoctoral fellowship and from Swift (80NSSC21K1409).S.N. acknowledges the partial support from NASA ATP 80NSSC20K0505 and from NSF-AST 2206428 grant as well as thanks Howard and Astrid Preston for their generous support.S.R. 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.E.R.R. thanks the Heising-Simons Foundation, NSF (AST-1615881 and AST-2206243), Swift (80NSSC21K1409, 80NSSC19K1391) and Chandra (22-0142) for support.The EKL mechanism is a resonant system at which the argument of periapsis (ω) and the longitude of ascending nodes (Ω) can undergo libration or rotation (e.g., Li et al. 2014;Hansen & Naoz 2020).The diffusive nature of the two-body relaxation can drive the system from libration to rotation (as shown in Naoz et al. 2022).In Figure 11 we illustrate the evolution of ω and Ω, for both EKL (+GR) and the combined EKL(+GR) and two-body relaxation effects.While the evolution of ω and Ω appear solid, they oscillate rapidly for the EKL (+GR) and two-body relaxation case.We remind the reader that this system for the combined EKL(+GR) and two-body relaxation case is a repeated TDE system.

B. ADDITIONAL MONTE CARLO RESULTS
In Table 1, we described the full Monte Carlo simulations.As an example, in Figure 5, we showed the results from simulations with binary eccentricity e b in = 0.5.Here we show additional systems undergoing two-body relaxation and EKL in Figure 12 and Figure 13.Like in 5, these additional figures show the final evolution of several systems with different pairings of mass ratio, stellar density distribution, and binary eccentricity.In Figure 12, the binary eccentricity is set to e bin = 0.3, and in Fig. 13, it is set to e bin = 0.7.Similar to 5, repeated TDEs are only produced in the cuspier stellar density environments.

Figure 1 .
Figure1.The hierarchical triple system.This hierarchical triple includes an SMBH binary of masses m1 and m2 on an outer orbit, and stellar cluster with components m on an inner orbit around m1.

Figure 2 .
Figure 2. Comparing characteristic timescales.Top panel: In orange we show the EKL timescale for an m1 = 10 7 M and m2 = 10 9 M SMBH binary and a population of stars of m = 0.8 M .In blue we show the relaxation timescale for the inner binary of m1 and m , allowing α to range between 1 and 2. The dotted grey line is the star's orbital period around m1, while the dashed line is the period of the SMBH binary.Bottom panel: Here we compare the relative change in angular momentum h.In orange h/∆h ≈ tEKL/P for the EKL mechanism and in blue h/∆h ≈ t rlx /P for two-body relaxation.Further description is provided at the end of Section 2. SeeNaoz et al.  (2022)  for a similar analysis.

Figure 3 .
Figure 3. Principal timescales in a SMBH binary.Here, we normalize the timescales by the two-body relaxation timescale (Eq.5) and depict the EKL and the 1pN timescales (Eq. 3 and Eq. 4, respectively) as a function of the mass ratio q.To highlight the dependency on the mass ratio, we explicitly denote m disruptor as the one that forms the TDE and m perturber as the far away SMBH companion.To generate this figure, we adopted m perturber = 10 8 M and vary m disruptor between 10 6 M and 10 10 M with e bin = 0.7 and e = 0.9.The binary is set at a bin = 1 pc and the star at a = 0.07 pc.The dashed line denotes q = 1.To the right of the dashed line we have the regime where m disruptor > m perturber , and the 1pN precession suppresses the EKL eccentricity excitations.For m disruptor < m perturber , on the other hand, EKL excitations are dominant.In the presence of two-body relaxation, high eccentricities can be excited to larger values (Figure2).Throughout the paper, we thus focus on the left regime, which can excite eccentricities more efficiently.
Figure 4.An example of a repeated TDE.The time evolution of a supermassive black hole binary is shown.The primary black hole is of mass m1 = 10 7 M with a secondary of m2 = 10 9 M .There is a stellar mass of m = 0.8 M that is gravitationally bound to the primary at a separation of 4.7 × 10 3 au.The top panel compares the changes to the semi-major axis and pericenter radius of the inner binary when the system undergoes only EKL (+GR) (yellow) and EKL (+GR) with two-body relaxation (blue).The tidal radius is shown in red, where in crossing this line, we consider the system to have become a TDE.The semimajor axis is shown to drop considerably in the EKL and two-body relaxation simulation, indicating that this system could potentially be a repeated TDE.

Figure 5 .
Figure5.Results of the Monte Carlo runs for e bin = 0.5.We show a collection of systems that have undergone two-body relaxation and EKL.In grey, we indicate the initial conditions for the inner binary.Each panel has distinct initial conditions: the top panels have a BH mass ratio of 10, while the bottom panels have a BH mass ratio of 100; left panels have α = 1 while right panels have α = 2. Systems that end near their starting point are marked in dark blue, while those to the right of the dashed light-blue line have crossed the threshold for maintaining their hierarchical triple system configuration and thus are marked as non-hierarchical in pink.Systems below the red line, resulting in the crossing of the tidal radius, are marked as TDEs in green.Green stars on the RT line are considered to be repeated TDEs on orbital periods P ≤ 30 yrs.As shown, more systems become tidally disrupted or non-hierarchical for an α of 2.
Figure 6.Comparing the number of stars to the number of tidally disrupted events.Here, the eccentricity of the outer binary is varied on values e bin = 0.3, 0.5, 0.7.Mass ratios of q = 10 and q = 100 are also compared and are indicated by circles and "×'s," respectively.Pairs of e bin and q are indicated by different colors, as seen in the legend.The stellar density profile is changed, and a dashed line shows α = 1 values to the left and α = 2 to the right.A linear trend is seen in the number of TDEs as e bin decreases and α increases.
Figure 7.Comparing the number of repeating tidally disrupted events to the total number of tidal disruption events.The eccentricity of the outer binary is varied on values e bin = 0.3, 0.5, 0.7.Mass ratios of q = 10 and q = 100 are represented as circles and "×'s," respectively.The colors represent the pairs of e bin and α.Here, unlike in Fig.6, only α = 2 is shown as only this cuspy profile produced rTDEs as seen in Fig.5, 12, 13.

Figure 8 .
Figure 8. Period distribution of expected rTDEs.The expected number of rTDEs, normalized by the M − σ relation, is shown for different pairs of binary eccentricity and mass ratio, top panel q = 10 and bottom panel q = 100.Arrows show the current observations of potential rTDEs.The shaded grey region spans the parameter space for which the Hill's mechanism produces rTDEs (see text for details).

Figure 9 .
Figure9.Calculated TDE rates.The rate for all systems is calculated and compared to both the average observed rate and the observed PSB rate.The left columns show systems with the density distribution power-law α = 1 while the right column panels show that of α = 2.The top panels have a mass ratio of q = 10, and the bottom panels q = 100.As seen, α = 2 systems produce TDEs at an enhanced rate with little dependence on eccentricity.

Figure 10 .
Figure10.TDE rates with and without two-body relaxation.Rates shown here are for systems with two distinct stellar distributions α = 1 in violet and 2 in cyan.The shaded rates are for systems that underwent both twobody relaxation and EKL (+GR), while the hatched greylike rates are for those with only EKL (+GR) processes as adapted fromMockler et al. (2023).The adapted rates showcase high peaks on shorter timescales, while rates in this work have lower peaks and are extended further in time.The minimum for all rates is calculated only for our conservative lower limit for TDEs, which are constrained within the Hill radius.The maximum for the EKL-only rates includes the systems that are beyond the Hill and the Roche limits.This limit corresponds to the solid line in the Figure.The maximum value for the rates for the EKL (+GR) + 2 body Relaxation includes the systems beyond the hierarchical limit as well.

Figure 11 .
Figure11.Repeated TDE evolution including Ω and ω.Additional orbital elements of the same system in Figure4illustrate the evolution of the argument of periapsis ω, third panel, and the longitude of ascending nodes Ω, bottom panel, both of which are taken from a uniform distribution between 0 and 2π.We compare the evolution of these orbital parameters as they undergo EKL (+GR) (yellow) to the combined effects of EKL (+GR) and two-body relaxation (blue).

Figure 12 .
Figure 12. Results of the Monte Carlo runs for e bin = 0.3.Similarly to Figure 5, we show a collection of systems that have undergone two-body relaxation and EKL, with the only change being the binary eccentricity reduced to e bin = 0.3.As shown, only systems with α = 2 produce rTDEs on orbital periods of 30 years or less.

Figure 13 .
Figure 13.Results of the Monte Carlo runs for e bin = 0.7.Similarly to Figure 5, we show a collection of systems that have undergone two-body relaxation and EKL with the only change being the binary eccentricity increased to e bin = 0.7.As shown, only systems with α = 2 produce rTDEs on orbital periods of 30 years or less.