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 populations around SMBHs as well as accretion processes onto SMBHs (e.g., Dai et al. 2021;Gezari 2021;Mockler et al. 2022). The rates of these events have recently been studied as possible probes in the search for ultra-light bosons (e.g., Du et al. 2022). These rates seem also to have many repercussions for SMBH and host-galaxy demographics (e.g., van Velzen & Farrar 2014;Kochanek 2016;Law-Smith et al. 2017;French et al. 2020;Dodd et al. 2021;Gezari 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 has received 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.

Another channel often considered in the literature involves the presence of an SMBH binary. The hierarchical nature of galaxy formation, combined with SMBHs residing in the centers of almost every galaxy, implies that SMBH binaries are prevalent in our Universe (e.g., Begelman et al. 1980;Di Matteo et al. 2005;Hopkins et al. 2006;Robertson et al. 2006;Callegari et al. 2009;Li et al. 2020). Furthermore, observations of dual active galactic nuclei, typically a few kiloparsecs or more apart, seem to suggest that these configurations will eventually lead to a tight SMBH binary (e.g., Komossa et al. 2003;Bianchi et al. 2008;Comerford et al. 2009Comerford et al. , 2018;;Green et al. 2010;Liu et al. 2010;Smith et al. 2010;Foord et al. 2020;Li et al. 2020;Stemo et al. 2021). Closer to home, theoretical arguments combined with observational campaigns insinuate that our Galactic center may also host a massive black hole companion (e.g., Hansen & Milosavljević 2003;Maillard et al. 2004; Gürkan & Rasio 2005;Gualandris & Merritt 2009;Chen & Liu 2013;Fragione et al. 2020;Generozov & Madigan 2020;GRAVITY Collaboration et al. 2020;Naoz et al. 2020;Zheng et al. 2020).

Gravitational perturbations from an SMBH companion can significantly modify the orbits of surrounding stars (e.g., Chen et al. 2008Chen et al. , 2009Chen et al. , 2011;;Chen & Liu 2013). Most notably, the Kozai-Lidov (EKL) mechanism (Kozai 1962;Lidov 1962;Naoz 2016) has been shown to excite the eccentricities of stars to high values (Li et al. 2014a(Li et al. , 2014b(Li 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 yr.

Recent observations of repeating tidal disruption events (rTDEs) have raised the question of why none of the aforementioned channels, two-body relaxation or EKL, have predicted a sizable number of rTDEs. In an rTDE, a star is partially disrupted and may experience multiple disruptive events (e.g., MacLeod et al. 2013MacLeod 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 disrupted? 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 semimajor 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 -10 -2 pc for an SMBH with mass in the range of 10 6.5-7.86 M e in order to have an orbital period ranging from 115 days to 30 yr (e.g., Payne et al. 2021;Evans et al. 2023;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 (e.g., Antonini et al. 2011;Cufari et al. 2022). 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 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 an 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 semimajor axis a bin and eccentricity e bin .

In this system, we select m 1 < m 2 , where the primary SMBH was fixed to m 1 = 10 7 M e , while the secondary SMBH mass was chosen as either m 2 = 10 8 or 10 9 M e for mass ratios of q = m 2 /m 1 of 10 or 100. Surrounding m 1 is a population of stars, each with mass m å ≈ 0.8 M e . The initial mass function in our Galactic center is suggested to be top heavy (peaking at higher masses), as well as possibly exhibiting a depletion of smaller-mass stars (e.g., Lu et al. 2013). Thus, our mass choice reflects our Galactic centerʼs nuclear star cluster. The stellar mass is treated as a test particle; therefore, changes in mass choice will exhibit no major dynamical effects. The stars are set on separation

), where a å , e å , and f å are the stars' semimajor 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 e and σ 0 = 200 km s -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)   = -a a e e 1 , 2 bin bin bin 2 ( )

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 e , which corresponds to R å = 0.7 R e . As mentioned, previous studies have shown 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. Throughout this section, we introduce relevant timescales to this system (as shown in Figure 2), and whose comparison is further discussed toward the end of the section. It has been suggested that the EKL mechanism is less efficient in driving high eccentricity of the stars around 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 first 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 Figure 3, 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. As shown in Bhaskar et al. (2022) and Grishin et al. (2018), overall the secular approximation underpredicts the eccentricity. 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 its circular angular momentum is estimated by (e.g., Binney & Tremaine 2008 6)) and depict the EKL and the 1 pN timescales (Equations ( 4) and (5), 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 faraway SMBH companion. To generate this figure, we adopted m perturber = 10 8 M e and vary m disruptor between 10 6 M e and 10 10 M e 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 1 pN 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.

where 〈m scat 〉 is the mass of the average star that acts as a scatterer, and σ å 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 e . We show this timescale (Equation ( 6)), in Figure 2, for different density profiles α, between 1 and 2. This timescale is explicitly for stellar circular orbits; noncircular stellar orbits are expected to be significantly shorter (e.g., Sari & Fragione 2019).

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. 2022a). 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 cusplike profile, α = 2, yields changes to the angular momentum h/Δh comparable, or even shorter than, the EKL's.

Following Naoz et al. (2022a), we model the change to the stellar orbit's velocity

) 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 SMBH. Each component of this threedimensional kick is drawn from a Gaussian distribution with a zero average and a standard deviation of Dv 3 , where

(see also Bradnick et al. 2017). See Naoz et al. (2022a) for the complete set of equations.

Numerical Simulations

Dynamical Evolution of a Representative Example

In Figure 4, we consider our fiducial model of m 1 = 10 7 M e and m 2 = 10 9 M e , 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 (in yellow), and EKL + GR + two-body relaxation processes (in 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 yr the system begins eccentricity oscillations, which change the inclination and semimajor 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 TDE 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 e and m å = 0.8 M e . The primary SMBH choice of m 1 = 10 7 M e is motivated by the inferred SMBH masses of the potential rTDE observations, which average to ∼10 7 M e (e.g., Payne et al. 2021;Evans et al. 2023;Liu et al. 2023;Malyali et al. 2023;Wevers et al. 2023). The main consequence of choosing a different mass will be to change the number of stars available in the TDE process, as was highlighted in Naoz & Fabrycky (2014) There is a stellar mass of m å = 0.8 M e that is gravitationally bound to the primary at a separation of 4.7 × 10 3 au. The top panel compares the changes to the semimajor 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.

We chose two representative masses for the companion of m 2 = 10 8 M e and 10 9 M e . We explore three outer-orbit eccentricities of e bin = 0.3, 0.5, and 0.7. The semimajor 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), which yields the semimajor axis distribution of the stars; see Figure 5. Moreover, we require that ò 0.1; see Equation (2), as well as the Mardling & Aarseth (2001) stability criterion. Motivated by our own Galactic center (e.g., Yelda et al. 2014;Gillessen et al. 2017), we set the minimum semimajor axis of the stars to 500 au. 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 and 2π. The mutual inclination was chosen from an isotropic distribution, i.e., uniform in i cos ; see Table 1 for the list of parameters chosen for the 12 realizations.

Figure 5 shows a representative example of the runs, where the gray points mark the initial conditions. In total, we see four distinct outcomes (the color code in Figure 5, as well as Figures 12 and 13 from Appendix B, is described below):

1. During the evolution the star's pericenter crossed 2R T ; see Equation (3). We mark these systems as stars that become tidally disrupted. These systems are marked in green points below the a å (1 -e å ) = 2R T line. The fraction of systems in which 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 disruption threshold of R peri = 2R T (β = 0.5) and have periods shorter or equal to 30 yr we mark as repeated TDEs ( f rTDE ).foot_1 These systems are highlighted as green stars. 3. Systems that have not crossed 2R T within 1 Gyr are considered as "survived" and marked in dark blue.foot_2 4. Finally, in some cases, the system violated the hierarchical condition, Equation (2). These systems are marked as magenta 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, seventh column ( f TDE + f nonhierarchical )).

Simulation Results

As discussed in the previous section, Figure 5 shows the trajectory of stars at the end of our simulations (see also Figures 12 and 13, Appendix B). We have summarized the variable parameters in our simulations and their impact on the fraction of TDEs and rTDEs in Table 1. Here we examine the impact of the stellar density distribution, the SMBH mass ratio, and the SMBH binary eccentricity on the evolution of our simulations.

Two-body relaxation is sensitive to the underlying density profile. This is further highlighted in Figure 5, which shows simulation results for α = 1 and α = 2, left and right panels, respectively, as described above. A shallow density profile (i.e., α = 1) yields a longer two-body relaxation timescale, thus minimizing the angular momentum changes per interaction (see Figure 2). Thus, in this case, we find less migration inwards (compared to the steeper, α = 2 case). The immediate consequence of the combined EKL and two-body relaxation is that repeated TDEs are more likely to take place for a cusplike (α = 2) distribution, as discussed in Section 4.2. Compared to the EKL-only scenario (Mockler et al. 2023), we still find that including two-body relaxation is required for the formation of rTDEs.

As shown in Figures 5, 12, and 13, the dynamics is nearly independent of the mass ratio and the eccentricity of the SMBH binary. In other words, the same qualitative TDE and rTDE results are duplicated across these parameters. The main contributor to the rate and the number of TDEs and rTDEs is the "available" number of stars within the hierarchical limit. The M-σ relation (Equation ( 1)) correlates to the number of stars within the sphere of influence with the mass of the SMBH. On the other hand, the hierarchical limit (Equation ( 2)) determines the number of stars based on the SMBH eccentricity, which ultimately influences the TDE rate. This is discussed below, in Section 4.1.

Predictions and Observational Signatures

Tidal Disruption Event 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:

1 0 0 2 3 ( ) and   =r a e e 1 , 10 max bin bin 2 ( ) ( )

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.

Table 1 Simulation Parameters Run # q e bin α f TDE f TDE + f rTDE  N r max ( ) f nonhierarchical 0 10 0.3 1 0.18-0.33 0.85 L 6e5 1 10 0.3 2 0.22-0.31 0.98 0.0156 4e6 2 100 0.3 1 0.035-0.37 0.82 L 6e5 3 100 0.3 2 0.15-0.47 0.99 0.0101 4e6 4 10 0.5 1 0.24-0.41 0.84 L 1e5 5 10 0.5 2 0.21-0.30 0.99 0.0260 2e6 6 100 0.5 1 0.056-0.51 0.87 L 1e5 7 100 0.5 2 0.18-0.48 1.0 0.0122 2e6 8 10 0.7 1 0.20-0.48 0.81 L 3e4 9 10 0.7 2 0.14-0.29 1.0 0.0238 9e5 10 100 0.7 1 0.034-0.62 0.87 L 3e4 11 100 0.7 2 0.087-0.50 1.0 0.0196 9e5 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), shown by the shaded regions in Figure 7. As seen, for α = 1, the spread of rates for the different eccentricities are consistent with both the average and PSB observed rates. In particular, the α = 2 rates are at least an order of magnitude larger than the observed rates. Note that the plateau in the figure is a result of adopting only a single star formation episode for the purposes of this exercise. Hopefully, in the near future, surveys such as the Vera Rubin Observatory and the LISA space mission will close the gap between our rate estimates and observations. 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 nonhierarchical systems (indicated as the light magenta dots to the right of the dashed lines in Figure 5). These nonhierarchical systems are likely to be excited to high eccentricities (e.g., Bhaskar et al. 2021) and may eventually become TDEs (similarly to the maximum rate of extreme mass ratio inspirals, EMRIs; Naoz et al. 2022).

We point out that while we initialize the particles within the hierarchical limit, the two-body relaxation may drive them beyond this limit (e.g., magenta dots to the right of the vertical line in Figure 5), thus potentially getting too close to the companion m 2 . However, as was shown by 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 semimajor 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 Equation ( 2)).

Consider first the comparison between the calculated rates as a function of α. As shown above, α is one of the key The rate for all systems is calculated and compared to both the average observed rate and the observed PSB rate. The left column panels 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.

parameters that affect the combined effect of EKL with twobody relaxation (see Figure 2). Thus, in Figure 7, we compare the results between α = 1 and α = 2, left and right columns, 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 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 have suggested that the stellar distribution in our own Milky Way's Galactic center is corelike rather than cusp-like, with α closer to unity (e.g., Schödel et al. 2018Schödel et al. , 2020;;Lu & Naoz 2019;Gallego-Cano et al. 2020). Thus, if our Galactic center is representative of galaxy nuclei's stellar distribution, the left column of Figure 7 may be a more representative scenario. In this case, the weak kicks due to twobody 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 PSB 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 Tidal Disruption Events

rTDE observations present many intriguing puzzles. For example, how do these stars get only partially disrupted? This question has been 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 semimajor axis 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 semimajor axis slightly increases. Furthermore, since the kicks are proportional to the star's orbital velocity (Equation ( 9)), 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 yr to mark the rTDEs candidates. As for the nonrepeating TDEs, the number of expected rTDEs is proportional to the number of stars in the sphere of influence (see Figure 8). 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 9 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 from hundreds of days up to 30 yr (e.g., Wevers et al. 2019Wevers et al. , 2023;;Payne et al. 2022;Evans 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 has been suggested that the Hills 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 Hills 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, for a binary to remain bound before the Hills mechanism takes place, its semimajor axis should be smaller than a critical value and larger than at least 2 times the radius of the binary members (Rose et al. 2020). Taking the minimum separation between a binary, i.e., assuming a hard binary in which its binding energy is larger than the kinetic energy of neighboring stars (Binney & Tremaine 2008) and more specifically 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 9, 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. Note that the minimum period is, in fact, a contact binary; thus, in practice, the shaded area in Figure 9 begins further to the left.

Comparison with Tidal Disruption Event without Twobody 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 EKLonly approach adopted from Mockler et al. (2023). The figure shows the time-dependent TDE rate for two examples of m 1 = 10 7 M e , m 2 = 10 8 M e , m å = 0.8 M e , e bin = 0.5, and a bin ≈ 2 pc, with the violet shaded region indicating the system with α = 1 and the cyan shaded region α = 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 timedependent rate, perhaps allowing star formation to occur and replenish the TDEs.

As expected, the EKL-only channel produces a burst-like event, depicted as a shorter rate in Figure 10. This behavior was noted in Naoz & Fabrycky (2014) and Naoz et al. (2022a). 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. Figure 10. 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 two-body relaxation and EKL (+GR), while the hatched gray-like 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, 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) + twobody relaxation includes the systems beyond the hierarchical limit as well.

Combined Tidal Disruption Event and Extreme Mass Ratio Inspiral Events

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 regarding density profile and the SMBH binary separation and eccentricity, Figures 6 and 7 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 8). Beyond TDEs, it was recently suggested that the similar physical processes described here could yield much higher EMRIs rates (orders of magnitude than estimated before; Mazzolari et al. 2022;Naoz et al. 2022a). Thus, by combining the two findings, we speculate that a combined TDE and EMRI event may have a nonnegligible 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 process, 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 stars' density profile by comparing two extreme cases, one of a 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. 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 PSB 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, and thus the angular momentum changes via the two-body scattering are smaller than via 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 8).

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 Appendix B for the configurations involving e bin = 0.3 and e bin = 0.7). As seen in Figure 7, 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 eccentricities. This is a direct result of the fact that the lowereccentricity SMBH binary configuration allows for more stars to reside in the hierarchical limit (see for reference Figure 6). Thus, as shown in Figure 7, the rate is largely insensitive to the mass ratio and the orbital eccentricity of the SMBH binary. 7We 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. 2022Rose et al. , 2023)). Additionally, a substantial fraction of stars are expected to be in binaries (see Raghavan & Steprāns 2012), 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. 2016Stephan et al. , 2019)). In general, we expect that binaries which undergo the combined effect of two-body 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 core-like 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 EMRI 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 7 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.

Appendix 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 bin = 0.5. Here we show additional systems undergoing two-body relaxation and EKL in Figures 12 and 13. Like in Figure 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 Figure 13, it is set to e bin = 0.7. Similar to Figure 5, repeated TDEs are only produced in the cuspier stellar density environments.