Uncovering Hidden Massive Black Hole Companions with Tidal Disruption Events

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

Tidal disruption events (TDEs) occur when a star passes close enough to a black hole that the tidal force from the supermassive black hole (SMBH) overwhelms the self-gravity of the star (e.g.Rees 1988;Evans & Kochanek 1989).The star is ripped apart and the stellar debris is accreted by the SMBH, producing a flare similar in brightness to supernovae that encodes properties of the star and SMBH in its light curve and spectra (e.g.Lodato et al. 2009;Ramirez-Ruiz & Rosswog 2009;Guillochon & Ramirez-Ruiz 2013;Yang et al. 2017;Law-Smith et al. 2019;Mockler et al. 2019Mockler et al. , 2022)).Theoretical tidal disruption event rates have most often been calculated using the dynamics of two-body relaxation (e.g., Magorrian & Tremaine 1999;Wang & Merritt 2004;MacLeod et al. 2012;Stone & Metzger 2016).In this picture, a star will experience random kicks to its angular momentum through interactions with other stars, and if it is very unlucky, one of these kicks will send it radially toward the SMBH, where it will be tidally disrupted.However, in a SMBH binary, 3-body interactions can produce TDEs much faster than two-body relaxation processes (e.g.Ivanov et al. 2005;Chen et al. 2011;Wegg & Nate Bode 2011).
In general, SMBH binary systems are very difficult to uncover observationally.If one or both of the black holes is accreting as an AGN, it is sometimes possible to measure the shift in velocity in the emission lines due to the binary orbit (Gaskell 1983;Boroson & Lauer 2009;Decarli et al. 2010;Tsalmantza et al. 2011;Eracleous et al. 2012;Runnoe et al. 2015;Wang et al. 2017).There are also predictions for how the overall shape of the broad emission lines might differ from a solitary AGN (Decarli et al. 2010;Nguyen & Bogdanović 2016).Unfortunately, measuring the velocity shift in emission lines is only possible over a limited range of SMBHB separations.The SMBHs must be close enough such that the velocity shifts in the broad line regions (BLR) are measureable, but not so close that the BLRs are truncated (see overview in De Rosa et al. 2019).
If we are lucky enough to catch both SMBHs actively accreting, this is likely the surest sign of a SMBH binary (SMBHB) system, however, resolving these systems requires the SMBHs to still be at large separations (before the system has formed a bound binary) and are also mostly limited to binaries with mass ratios near q = 1 (e.g., Rodriguez et al. 2006;Comerford et al. 2018;Dey et al. 2019;Goulding et al. 2019;Bhattacharya et al. 2020;Foord et al. 2020;Kim et al. 2020;Reines et al. 2020;Monageng et al. 2021;Severgnini et al. 2021;Shen et al. 2021).Additionally, even though the merger process can funnel gas to the center of galaxies and trigger AGN activity, most SMBHBs in the local universe do not currently host AGN (Casey-Clyde et al. 2022).
When stars are disrupted by the smaller SMBH in the binary, the black hole mass measured from the TDE light curve or spectra will correspond to the smaller SMBH (e.g.Mockler et al. 2019;Ramsden et al. 2022;Wen et al. 2020;Ryu et al. 2020).Given that the potential of the galactic nucleus will be dominated by the larger SMBH, typical scaling relations are likely to produce estimates close to the mass of the larger black hole in the binary (Izquierdo-Villalba et al. 2022).Therefore, if TDE light curve measurements suggest a much smaller black hole mass than scaling relations such as M − σ or the bulge-mass relation, this is an indication that the galaxy might host a SMBH binary1 .Similarly, if a TDE is observed in a galaxy that appears to be too massive to host an observable flare (i.e. a system with an expected black hole mass ≳ 10 8 M ⊙ ), this also suggests the existence of a second, smaller SMBH in the galactic nucleus (e.g.Coughlin & Armitage 2018).
In this paper, we study how hierarchical 3-body interactions between the SMBH binary and the stars orbiting the smaller black hole can excite the eccentricities of the stars and increase the rates of tidal disruption events around the smaller black hole.In particular, we focus on the effect of the secular orbital perturbations known as the eccentric Kozai-Lidov (EKL) mechanism (Kozai 1962;Lidov 1962;Naoz 2016).It was shown that SMBH binaries and the EKL mechanism can lead to increased rates of TDEs.However, most previous work on this subject has either focused on the stars surrounding the larger SMBH (e.g.Ivanov et al. 2005;Chen et al. 2011;Wegg & Nate Bode 2011;Fragione et al. 2021), on IMBH systems (Fragione & Leigh 2018), or on extreme mass ratio inspiral (EMRI) transients instead of TDEs (Naoz et al. 2022;Mazzolari et al. 2022).Past work has also shown that other dynamical mechanisms suppress EKL oscillations for stars around the larger black hole, and the majority of TDEs that are produced around the larger SMBH come from chaotic orbital crossings (e.g., Chen et al. 2011).This is generally not the case for the stars around the smaller SMBH, largely because the EKL perturbations due to the more massive SMBH produce stronger excitations in the stars' orbital eccentricities.Here we expand on work by Li et al. (2015) and explore the effect of the EKL mechanism on the rates of TDEs around the smaller SMBH in SMBH binaries.
We propose that every SMBH binary system will undergo a burst (short, enhanced period) of TDEs from the smaller black hole once the binary separation is of order the sphere of influence (so long as one of the SMBHs is below the Hills mass for the star, ≈ 10 8 M ⊙ ).At these separations, the binary will not be resolved spatially, and the TDE will appear to be at the center of the galaxy.TDEs from the smaller black hole provide a unique avenue to discover SMBHBs at much closer separations than the kiloparsec scales at which dual AGN are observed.This is especially significant when the more massive SMBH is unable to disrupt sun-like stars outside its event horizon (as also discussed in Fragione & Leigh 2018) .
Cartoon of simulation setup (not to scale).We study stars orbiting the smaller of the two SMBHs, and focus on the region within the hierarchical radius (defined in Equation 2), as orbits in this region are less likely to be chaotically scattered by the second black hole.

Physical processes
We consider a SMBH binary m 1 < m 2 , separated by a semi-major axis a bin , and with an eccentricity e bin .Throughout the paper, the subscripts '1', '2', and '*', denote the inner black hole, the outer black hole, and the star, respectively.The subscript 'bin' refers to the black hole binary.Surrounding each SMBH is a nuclear star cluster, and for the purposes of this paper we will focus on the stars surrounding m 1 .Because m 1 < m 2 , we can think of the nuclear star cluster of m 1 as embedded in the larger cluster around m 2 .We expect the smaller SMBH m 1 to retain stars within its Roche lobe (e.g.O' Leary & Loeb 2012;Rantala et al. 2018;Greene et al. 2021).Each star in the cluster around m 1 can be modeled as part of a hierarchical three-body system made up of an inner binary pair -m * and m 1 -and an outer binary pair -m 1 and m 2 (see illustration in Figure 1).We assume the stars orbiting m 1 follow approximately Keplerian orbits, but over time their orbits can be secularly perturbed by m 2 through the eccentric Kozai-Lidov mechanism (e.g. Naoz 2016).The dashed horizontal line denotes our critical radius for partial disruption (β = 0.5), and the dotted vertical lines denote t quad given the initial orbital conditions.As you can see, t quad is of order the period of one eccentricity oscillation, as expected.This simulation is part of the set of runs with m1 = 10 6 M⊙, m2 = 10 7 M⊙, α = 2.
We restrict the initial parameter space to stable configurations, adopting a hierarchical limit that ensures the orbit of the star stays well within the orbit of the binary.
where ϵ is the pre-factor in front of the octupole expansion term of the Hamiltonian (e.g., Naoz et al. 2013) 2 This defines a 'hierarchical radius' limit: The EKL mechanism describes the coherent perturbations from m 2 on the stars orbiting m 1 .These perturbations exchange angular momentum between the outer SMBH binary (m 1 and m 2 ) and the inner binaries made up of stars orbiting m 1 .This produces excitations in the eccentricity and inclination of the stars orbiting m 1 over a wide range of the parameter space (e.g., Naoz et al. 2011;Katz et al. 2011;Lithwick & Naoz 2011;Naoz et al. 2013;Li et al. 2014;Will 2017).Importantly, the eccentricity excitations can lead to the production of a large number of tidal disruption events (e.g., Li et al. 2015;Fragione & Leigh 2018).To model these systems and determine the rate of tidal disruption events, we employ the secular approximation, where the inner and outer orbits are treated as phase-averaged wires that torque each other.We solve the secular, hierarchical three-body equations following Naoz et al. (2013), including up to the octupole order terms (O(a * /a bin ) 3 ).The timescale associated with the quadrupole level of approximation is (e.g., Antognini 2015): where G is the gravitational constant.This timescale is associated with the period of oscillations of the eccentricity and inclination (see Figure 2).Li et al. (2015) considered the relevant dynamical timescales at the center of galaxies and compared them to the EKL timescale.They found that the dynamical processes competing with the EKL mechanism for the largest effect on the stellar orbital parameters are general relativistic precession (due to the potential of the disrupting black hole), and Newtonian precession (due to the potential of the stellar cusp).The timescale for general relativistic precession is, where c is the speed of light.The Newtonian (NT) precession timescale is given by: where ψ is the true anomaly of the inner orbit, r(ψ) = a * (1 − e 2 * )/(1 + e * cosψ) from Kepler's equation.The function M * (r) is the integrated stellar mass within r(ψ) -we will discuss the density profiles used to calculate this in Section 2.2.If precession due to GR or the stellar potential is strong enough, it can prevent EKL perturbations from building up coherently and exciting the eccentricity and inclination of the star (e.g., Li et al. 2015;Chen et al. 2011).
Motivated by a similar analysis in Li et al. (2015), in Figure 3, we look at how changing different parameters of the system changes the number of stars likely to experience large eccentricity oscillations from the EKL mechanism.For each point in Figure 3, we calculate the number of stars within our hierarchical limit (Equation 2) whose shortest dynamical timescale is t quad , and plot density contours displaying this information.The dynamical timescales we compare against are described in the figure caption (and included in Appendix A).Like Li et al. (2015), we find that t GR1 and t NT are the two shortest timescales apart from t quad , implying that general relativistic precession due to the disrupting black hole and Newtonian precession due to the stellar cusp are the two dynamical processes most likely to interfere with EKL oscillations.This parameter space exploration also further motivates our focus on the smaller black hole in the binary.The number of stars susceptible to EKL oscillations drops off sharply as the disrupting black hole grows to a comparable size to the outer black hole (previous work focusing on the larger black hole found that the majority of disruptions were caused by chaotic orbital crossings, not the EKL mechanism, e.g., Wegg & Nate Bode 2011;Chen et al. 2011).
Orbital precession due to GR and the stellar cusp can be important at small radii, however the strength of EKL perturbations increase with a star's radius from m 1 as stars move closer to the perturbing black hole (m 2 ).Therefore it is more difficult for the EKL mechanism to excite the eccentricities of stars at smaller semimajor axes, but at larger radii, where there are many more stars, the EKL mechanism often dominates and can produce many TDEs.For this work we chose to isolate the effects of the EKL mechanism and used the parameter space exploration described by Figure 3 to choose systems to model where the EKL timescale t KL is the shortest dynamical timescale for the majority of stars surrounding the smaller black hole.
We include GR precession in our modeling, and while it is possible that it suppresses the number of disruptions at small radii close to m 1 (see Figure 6), we find that the EKL mechanism acts much faster than GR precession for the majority of the stars around m 1 in the systems we are modeling, as expected (see for similar results Naoz & Silk 2014;Li et al. 2015;Naoz et al. 2022).Unlike Li et al. (2015), we do not include the NT precession from the stellar cusp in our models, and instead chose our setups such that the timescale for NT precession is longer than the timescale for EKL excitations for the majority of the stars.Thus, our results are consistent with Li et al. (2015) (see Section 3).Finally, we note that we do not include any star-star collisions or compact object collisions.We expect that these processes will generally become important on longer timescales (e.g., Rose et al. 2022Rose et al. , 2023)).
Another possible dynamical mechanism that can affect EKL perturbations is two-body relaxation.Naoz et al. (2022) showed that even though the canonical timescale for two-body relaxation to significantly change stellar orbits is many orders of magnitude larger than the Kozai-Lidov timescale in the relevant parameter space (consistent with our calculations for Figure 3), individual (incoherent) kicks from two-body scattering The white stars represent the parameter combinations chosen for numerical simulation runs.This was calculated by determining the range of semi-major axes for each parameter combination where the EKL timescale is the shortest dynamical timescale.All of our systems have e bin = 0.5, a bin = 0.5r h, disrupting .While we compared against a variety of dynamical timescales (including those for two-body relaxation, vector resonant relaxation, general relativistic precession, and binary hardening through gravitational wave emission, see Appendix A), the two competing timescales were tGR, 1 and tNT.These are the timescales for general relativistic (GR) precession due to the disrupting black hole and Newtonian precession due to the stellar potential, respectively.
Table 1.Simulation parameters and results.See Section 2.3 for a description of how upper and lower limits were calculated for the rightmost three columns.can affect the build up of Kozai-Lidov perturbations.In this paper we focus on the effects of the EKL mechanism alone, however, in a companion paper we explore the effect of two-body relaxation on the EKL mechanism (see Melchor et al. 2023).

Monte Carlo dynamical simulations
Using Figure 3 as our guide, we chose a variety of parameter combinations for our simulations that we expected to produce observable tidal disruption events.The parameters we varied were: the mass of the disrupting and perturbing black holes (m 1 and m 2 ), their mass ratio (q = m 2 /m 1 ), and the power law constant of the density profile for the stellar cusp (α).The mass of m 1 was set to a range of reasonable values for tidal disruption events: 10 5 , 10 6 , and 10 7 solar masses.Black holes above ≈ 10 8 M ⊙ will swallow most stars whole, and disrupt them within their event horizons.Black holes below 10 5 M ⊙ can and do disrupt stars (Law- Smith et al. 2017), however these disruptions may be more difficult to observe (e.g.due to slower circularization and lower Eddington limits, Guillochon & Ramirez-Ruiz 2015), and almost all confirmed TDEs have occurred around black holes above 10 5 M ⊙ (e.g., Ramsden et al. 2022;Hammerstein et al. 2022, with the exception of Angus et al. 2022).Therefore, we focus on this mass range for the purposes of this paper (see Fragione & Leigh 2018 for simulations of similar systems focusing instead on IMBH-SMBH binaries).We use mass ratios q = m 2 /m 1 of 10 and 100, and this determines the mass of m 2 .
The stellar cusp is modeled using the density distribution described in O' Leary et al. (2009): and our sphere of influence is defined where the mass enclosed is equal to 2 × m 1 , The density profile is written as a power law of semimajor axis, with the integrated stellar mass normalization fixed by the M − σ relation.The constants M 0 = 1.3 × 10 8 M ⊙ , σ 0 = 200 km s −1 , and k = 4 are the best-fit values for the M − σ relation in Tremaine et al. (2002) We vary the power law exponent 'α' between 1-2 to explore how the steepness of the stellar cusp affects TDE rates.While our own galactic center appears to have a flatter cusp (α ∼ 1, e.g., Schödel et al. 2018), TDE host galaxies tend to be very centrally concentrated and to have stellar light profiles indicative of steep stellar cusps down to the size scales that are observable (e.g., Law-Smith et al. 2017;French et al. 2020;Dodd et al. 2021).  .The final distributions of the stars' semi-major axes (a1) and maximum eccentricities (emax, plotted as 1 − emax) for 3 runs with varying inner black hole mass (m1 = 10 5 , 10 6 , 10 7 M⊙).The critical radius for disruption used here is Rp ≤ 2RT (allowing for both partial and full disruptions), and denoted by dashed lines in the plot.Disruptions are plotted with a star-shaped marker, and plunges (disruptions that occur with RT < RS) are shaded in gray.The maximum semi-major axis included for stars in a given run is defined by the hierarchical condition (Equation 1).
We set the eccentricity of the SMBH binary to an intermediate eccentricity of e bin = 0.5 and its semi-major axis to half of the sphere of influence of the disrupting black hole (a bin = 0.5r h, 1 ).This puts the SMBHs at separations of ∼ 1 pc (a bin = 0.2 − 1.9 pc for m 1 = 10 5 − 10 7 M ⊙ ) -a regime where the EKL mechanism is highly effective, however the SMBH binary's hardening timescale is still long enough that we do not have to worry about the binary's orbit changing significantly over the course of our simulations.At these radii, the SMBH binary's hardening timescale is mainly dependent on the effectiveness of loss-cone diffusion, which is likely ≳ 100 million years (e.g.Kelley et al. 2017), much longer than the EKL timescales for the majority of the stars in our simulation setups.
For each set of parameter combinations, we ran a 1000 numerical simulations of the 3-body interactions with different initial conditions for the star.The priors for the star's initial conditions were determined as follows: The semi-major axis of the star was drawn randomly from the density distribution, with a maximum radius determined by the hierarchical condition of the system (see Equation 1).The star's eccentricity was drawn from an isothermal distribution, and the mutual inclination was drawn from an isotropic distribution (uniform in cos(i)).The argument of periapsis and the longitude of the ascending nodes were drawn from a uniform distribution between 0 − 2π.See Table 1 for an overview of the simulation parameters.
The simulations were run until either the star was disrupted or the system was evolved for 10 9 years (by which point the TDE rate from the EKL mechanism had long since dropped off, see Figure 5).Our criteria for a tidal disruption to occur was that the star passed within twice the tidal radius (R p = 2R T ) of the black hole This meant that our minimum impact parameter for disruption was β = R T /R p = 0.5.Setting our TDE radius condition to 2R T instead of R T ensured that we captured partial disruptions as well as full disruptions 3 .Each set of a 1000 simulations provides a random Monte Carlo sampling of the stars around the smaller black hole in each binary system, given the priors described above.We plot the final eccentricity distributions of a subset of our runs in Figure 4, where the runs that end in disruption are clearly visible as a buildup of points at high eccentricities (low values of 1 − e).
Our simulations give us the fraction of stars that result in TDEs for each parameter combination, as well as 3 Below β ≈ 0.5 very little mass will return to the black hole and it is unlikely the TDE will be observable (e.g., Guillochon 2013) the time of disruption for each random sample.We also find that a number of stars plunge directly into the black hole, with pericenter radii less than the Schwarszchild radii of the black holes (shaded in gray in Figure 4).We note that if the SMBH is spinning, the star may zoom around the ergosphere instead of plunging directly into the black hole (e.g., Glampedakis & Kennefick 2002;Healy et al. 2009;Schnittman 2015).This can continue for a few tens to thousands of times the pericenter passage (see figure 3 in Naoz et al. 2019).
Finally, we can scale our results for the fraction of stars disrupted by the expected total number of stars around m 1 within the hierarchical radius (using Equation 6) to compare to real galaxies.The simulation parameters and the fraction and number of stars that become TDEs are quoted in Table 1 for each of our runs.We can also calculate a TDE rate as a function of time by binning our simulation results as a function of the time of disruption.We plot our scaled, time-dependent rates for a subset of the runs in Figure 5 (using Gaussian processes to smooth the distributions).

Limits and stability
For each set of simulations, we calculate an upper and lower limit for the TDE fraction.The lower limit assumes that any star whose orbit evolves such that its apocenter moves outside the Hill radius of the smaller black hole is scattered away from the smaller black hole and is not disrupted.We define the Hill radius at the black hole binary's pericenter to ensure our lower limit is sufficiently conservative, The upper limit assumes that all stars (within our hierarchical radius limit, see Equation 2) on disrupting orbits become TDEs, even those whose apocenters move outside the Hill radius before disruption.While the Hill radius is a good approximation of the stability boundary on long timescales (e.g., Grishin et al. 2017;Bhaskar et al. 2021;Tory et al. 2022), stars at larger radii have shorter EKL timescales, and can easily be disrupted before their orbits become unstable.Stars whose apocenters are just outside the Hill radius are not instantaneously scattered -the impulse felt on the orbit from the second black hole must build up over time before the star's orbit is significantly affected (Zhang et al. 2023).
It is also possible for stars to remain "Lagrangestable" even when they are no longer "Liapunov-stable" (e.g., Hayashi et al. 2022), meaning that stars can perform chaotic figure eight orbits around both black holes but remain bound to the system.In this situation, a star can become unbound from the first SMBH only to return later and get disrupted.
Finally, three-body scattering experiments have shown that even when stars are far outside the Hill radius, significant fractions of stars can be disrupted through chaotic orbital crossings with the second SMBH (although this has mostly been studied for stars orbiting the larger black hole, e.g.Chen et al. 2011).Given that we are focused here on studying TDEs disrupted through the EKL mechanism, we only consider disruptions within the hierarchical radius, which is within a factor of two of the Hill's radius in our simulations with q = 10 (and within a factor of four for simulations with q = 100), even for eccentricities approaching 1.
Thus, we conclude that using the Hill radius as a lower limit produces estimates on the TDE rate that are too conservative, and the true rate is somewhere in between the limits used in this paper.We note that a similar approach was taken in Naoz et al. (2022) and Melchor et al. (2023).See the former study's appendix for an exploration of the particles' trajectories.

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

Rates
As you can see in Figure 5, we obtain the expected 'burst' of TDEs in all of our simulations.The bursts peak around ∼ 10 5 − 10 6 years but remain higher than the global observed rate (≲ 10 −4 tdes/galaxy/year) for ∼ 10 6 − 10 8 years.This is comparable to or longer than the burst timescales predicted in Chen et al. (2011) for the more massive black hole companion (e.g.comparing our timescales for 10 6 M ⊙ secondaries to Chen et al. (2011) for 10 7 M ⊙ primaries).The burst is limited by the eccentric Kozai-Lidov timescale in that it will not begin until the stars' eccentricities have grown large enough for disruptions to occur.However, the burst can extend much longer than the eccentric Kozai-Lidov timescale (which is generally of order 10 4 − 10 6 years for these systems), as stars often go through many EKL cycles before they are disrupted.
We also note that we do not include the evolution of the black hole binaries in these simulations.If we modeled them as they evolved, we would expect the EKL 'loss-cone' to move from larger radii to smaller radii, as stars at larger radii are disrupted or scattered and the perturbations grow stronger at small radii.Therefore, we would likely expect a more extended rate if we included this binary evolution.This is what Chen TDE rate (dn/dt) m 1 = 10 6 , = 1.5 q = 10.0 q = 100.0observed rate observed rate PSB 2-body cuspy 2-body core Figure 5. TDE rates calculated for the EKL mechanism compared to rates calculated for 2-body relaxation and to observed rates.The 2-body relaxation rates plotted for 'cuspy' galaxies (dashed lines) correspond to galaxies with an inner density profile (ρ ∝ r −α ) α > 0.5, while the rates plotted for 'core' galaxies (dotted lines) correspond to α < 0.3 (Stone & Metzger 2016, note that density profiles are estimated observationally from radii > r h , and therefore are not a direct comparison to the density profiles used in simulations in this paper for radii interior to r h .).The observed rate for all TDEs (shaded in gray) is from van Velzen et al. ( 2020), and the observed rate in post-starburst galaxies (hatched region) is from French et al. (2020).
et al. ( 2011) found in their n-body simulations which included binary hardening due to stellar scattering (and Mazzolari et al. 2022, for similar simulations focused on EMRIs).Initially, while the disruption timescale was comparable to the hardening timescale, the peak of the burst was extended, but then when the hardening timescale grew longer than the tidal disruption timescale and the loss cone could no longer be replenished, the rate dropped off (although the authors note that loss-cone diffusion processes and gravitational wave emission were not included, and would likely result in further hardening at later times and a longer-lived plateau in the rate).

TDE fraction
The fraction of stars that become tidal disruption events depends on the number of stars within the hierarchical limit and the strength of the EKL mechanism i.e., the maximum value of the eccentricity excitations.In the test particle limit, the eccentricity can be excited to extreme values (e.g., Li et al. 2014;Naoz & Silk 2014), naturally driving the stars to low angular momentum orbits.We find that for the majority of our simulations, the upper limit on the percentage of stars that are disrupted is ≈ 30 − 45% (within the maximum radius defined by our hierarchical condition, Equation 2).This limit does vary slightly with the parameters chosen.
As depicted in Figure 6, the fraction of stars disrupted decreases as m 1 increases, particularly at small radii.This is mostly because GR precession is faster than the EKL precession close to the black hole.On the other hand, the fraction of stars disrupted at small radii increases when m 2 increases because the strength of the EKL oscillations increases with the mass of the perturber.The number of disruptions at large radii (close to the hierarchical radius) does not vary significantly because we are effectively running out of stars at those radii.Near the hierarchical radius, most of the stars in the inclination range susceptible to EKL perturbations are getting excited and disrupted, so increasing the strength of the perturbations is not noticeably increasing the number of disruptions.As there are much fewer stars at small radii (M * ,contained ∝ r 3−α ), increasing the mass ratio to 100 only increases the upper limit on the total fraction of disrupted stars in our simulations by ≈ 2% (see Table 1).
As expected, the density profile power law (α) affects the maximum fraction of TDEs, because changing α changes the relative fraction of stars near the hierarchical radius, where most of the disruptions originate.The number of stars disrupted is proportional to 1/t KL , which is, in turn, proportional to a 3/2 * .This can be used to estimate the upper limit on the total fraction of stars disrupted.The number of stars disrupted as a function of radius is dependent on both the eccentric Kozai-Lidov timescale and the density profile, and therefore goes as Then, the total fraction of stars disrupted goes as the integral of n disrupt (r) divided by the total number of stars and so the upper limit on the total fraction of stars disrupted (integrated over the radius range) is proportional to: TDE fraction (upper limit) While changing the mass of both black holes and the steepness of the stellar cusp affects the fraction of stars disrupted, it does not change the maximum TDE fraction by more than ∼ 0.1 (see Table 1).However, changing these parameters does dramatically change the number of stars within the hierarchical radius (our maximum radius, see Equation 1), and therefore the number of stars that can become TDEs.Because of this, we find that the strongest determinant of the maximum number of TDEs in the system is simply the number of stars within the hierarchical radius (N * (r max ) ∝ r 9/2−α max )4 .This is shown clearly in Figure 7.
On the other hand, the lower limit on the number of TDEs is strongly dependent on the mass ratio q.The lower limit on the number of TDEs is proportional to the number of stars that are not susceptible to scattering by the perturber -stars that remain within the Hill radius (r Hill ∝ q −1/3 ).The number of stars within the Hill radius follows the proportionality N * (r Hill ) ∝ r 9/2−α Hill .Combined with the fact that r Hill ∝ q −1/3 , the lower limit on the TDE fraction goes as5 , TDE fraction (lower limit) ∝ q It drops from ≈ 10% for q = 10 to ≈ 1% for q = 100 (m 1 = 10 6 , α = 1.5).We note that we find similar qualitative results to that of Li et al. (2015) (see following section), however, here we explore a wider part of the parameter space.
As we discussed in the previous section, our lower limit on the stability radius is likely too conservative.Our lower limit is approximately equal to the strongest limits found in Grishin et al. (2017) and Tory et al. (2022) (a * ≲ 0.5r H ).However, some of the stars outside this radius are disrupted on very short timescales, likely before their orbits would become unstable.Recently, Zhang et al. (2023) determined an analytical stability criterion for 3-body systems that they found to be consistent with n-body simulations.Adopting the Zhang et al. (2023) stability timescale definition means that stars from larger radii can undergo EKL eccentricity excitations before becoming unstable.This would increase the lower limit on the number of TDEs by a factor of ≈ 2 across all of our simulations.Therefore, it is likely that many of these stars will be disrupted before their orbits become unstable.It is also likely that many of the stars on unstable orbits will also become TDEs, even though analytical approaches are unable to model their evolution (as discussed in Section 2.3).

Comparison with previous work
We compared our results with Li et al. (2015), as they used a similar setup to our own.In their simulation, m 1 = 10 7 M ⊙ , m 2 = 10 8 M ⊙ , α = 1.75, e 2 = 0.5 and a 2 = 0.5pc, and they find 5.7% of stars within the Hill's radius are disrupted (using the more lenient criteria that a star must only pass within 3R T to be disrupted).This is comparable to the results for our simulations numbered 5-7, which were run with similar parameters and found that ≈ 5% of stars within the Hill's radius are disrupted.To compare directly, we also ran one additional test simulation with the same parameters (and disruption criteria) as Li et al. (2015), and found that 4.2% of stars within the Hill's radius were disrupted.We note that Li et al. (2015) include the effect of orbital precession due to the stellar cusp in their simulations but that our simulations appear consistent nonetheless, justifying our exclusion of this effect for this set of simulations (see Section 2 for further explanation).
We also compared our results to the 3-body scattering simulations of the eccentric Kozai-Lidov mechanism around the larger SMBH in Chen et al. (2011) and Wegg & Nate Bode (2011).Both of these works showed that the EKL oscillations induced by a smaller SMBH on stars surrounding a larger companion are suppressed for large parts of the parameter space.In Figure 5, we show that we get a comparable number of TDEs from the smaller black hole in our q = 10 simulations as predicted for the larger black hole by scaling relations calculated from the results in Chen et al. (2011) and Wegg & Nate Bode (2011).

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

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

Caveats
To begin with, only a fraction of galaxies will host close SMBH binaries at any given time.Given the uncertainties on SMBHB hardening timescales (e.g.Kelley et al. 2017), it is not straightforward to estimate the current population of SMBH binaries from the rate of galaxy mergers, however, an overview is still useful for understanding possible scenarios.
The current day merger fraction (as defined by the fraction of galaxies with close companions) is ≲ 3% for galaxies ≲ 10 11 M ⊙ (this is the major merger fraction, the minor merger fraction is likely similar, Bundy et al. 2009;López-Sanjuan et al. 2012).However, it is possible that the current-day SMBHB population reflects a higher merger fraction from earlier times, depending on the hardening timescale of the binary (Bluck et al. 2009, 2012 find a major merger fraction of up to ≈ 25 − 40% at z = 3, which is ≈ 10 10 year lookback time).If the SMBHB population is reflective of these higher merger rates, upwards of 25% of galaxies might host two SMBHs at a given time (Rodriguez-Gomez et al. 2015), but their relevance to this work would still depend on their hardening timescales.
In addition to requiring that the host galaxies and their SMBHs be at a specific evolutionary stage, the actual burst of TDEs from the EKL mechanism is shortlived.If the stars are not replenished rapidly through star formation or dynamical mechanisms, the rate will drop off relatively quickly and the elevated rate will only be observable for tens of millions of years for anything m 1 = 10 5 ; = 1.50 m 1 = 10 6 ; = 1.00 m 1 = 10 6 ; = 1.50 m 1 = 10 6 ; = 2.00 m 1 = 10 7 ; = 1.00 m 1 = 10 7 ; = 1.50 m 1 = 10 7 ; = 2.00 q = 10 q = 100 . LHS: The number of TDEs as a function of m1, with different density profiles represented as different marker shapes (q=10 in all cases).The points in orange are the lower limits, excluding all stars that move outside r Hill , while the points in blue are the upper limits and include these stars.The blue and orange shaded regions are plotted by determining the slopes between the fraction of TDEs for the α = 1 & 2 runs at the black hole masses used in the simulations and then extrapolating out to higher and lower masses.The black hatched region is the number of TDEs disrupted by m2 in the same system from chaotic scatterings (the upper limit is from Liu et al. 2013 and the lower limit is scaled down to match Wegg & Nate Bode 2011).The gray shaded region is the number of TDEs predicted from 2-body relaxation over 10 8 years for a single SMBH system with mass m1 (from Stone & Metzger 2016).RHS: The number of TDEs as a function of the number of stars within the hierarchical radius.Each m1 and α combination is plotted as a different color, and the marker shape is determined by the mass ratio q.
The upper and lower limits are determined in the same way as for the LHS.
but the most massive systems (see Figure 5).Therefore, it is possible that even in galaxies with close SMBH binaries, we will miss the 'burst' of TDEs predicted here.Finally, we note that although we focus on the tidal disruptions around the smaller of the black holes in the SMBH binary, the global rate will be limited by the number of larger mass SMBHs available to form binaries because the black hole mass function decreases with increasing mass (e.g.Shankar et al. 2016;Gallo & Sesana 2019).

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

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

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

SUMMARY OF KEY FINDINGS
Our key results are summarized as follows: • The EKL mechanism produces a burst of TDEs around the smaller black hole in a SMBHB at a rate that is orders of magnitude higher than the rate from two-body relaxation for timescales of millions to tens of millions of years (see Figure 5).Unlike EKL oscillations around the larger SMBH, the EKL oscillations around the smaller SMBH are not significantly suppressed by precession due to general relativity or the stellar cusp.This is mostly because of the increased strength of the perturber.
• The number of disruptions is primarily dependent on the number of stars within the hierarchical (or Hill) radius.This is because the fraction of stars disrupted does not change dramatically with most parameters (so long as we choose regions of parameter space such that the EKL mechanism is the dominant dynamical mechanism, see Figure 7, Figure 3, Table 1).The timescale is most strongly dependent on the black hole mass and mass ratio (see Figure 5).
• The EKL mechanism generates a larger fraction of disruptions around the smaller black hole (∼ 3% − 45% of stars within the hierarchical radius are disrupted) than chaotic orbital interactions are able to produce around the larger black hole (e.g.Chen et al. 2011;Wegg & Nate Bode 2011).This can lead to a comparable number of disruptions around both black holes (see Figure 7).
• Disruptions around the smaller black hole can provide evidence of the existence of the SMBH binary as their light curve timescales will likely be shorter than what is expected from the larger black hole, which dominates the potential (Mockler et al. 2019).
We have shown that tidal disruptions from the secondary black hole in a SMBHB contribute meaningfully to the rate from the binary system.We emphasize that this provides a novel way to find SMBHB candidates at separations that are too small to resolve distinct sources. We

2013) APPENDIX
The equations above represent timescales for the following physical processes: t quad is the quadrupole eccentric kozailidov timescale (Naoz et al. 2013;Antognini 2015); t GR, 1 & t GR, 2 are the first-order post-Newtonian GR precession timescales for the inner and outer black holes; and t GR, int is the timescale of the first order post-Newtonian GR interaction between the inner and outer orbits.As defined in Kocsis & Tremaine (2011), t NT is the Newtonian precession timescale caused by the potential of the stellar cusp, t RR, v is the vector resonant relaxation timescale, and t LT is the Lense-Thirring precession timescale.Finally, t GW is the binary SMBH orbital decay timescale due to gravitational wave emission.In equation A6, Ω is the orbital frequency of the star and f vrr is estimated to be 1.2 by Kocsis & Tremaine (2015).In equation A5, ψ is the true anomaly of the inner orbit, r(ψ) = a 1 (1 − e 2 1 )/(1 + e 1 cosψ) from Kepler's equation.The function M * (r) is the integrated stellar mass within r(ψ).In equation A7, lnΛ is the Coulomb logarithm, and σ * is the velocity dispersion of the bulge.In equation A8, the constant (G 2 m 2 0 s)/c is the spin angular momentum of the inner SMBH (Peters 1964;Kocsis & Tremaine 2011;Naoz et al. 2013).

Figure 2 .
Figure 2. Example run showing EKL cycles as a function of time.The dashed horizontal line denotes our critical radius for partial disruption (β = 0.5), and the dotted vertical lines denote t quad given the initial orbital conditions.As you can see, t quad is of order the period of one eccentricity oscillation, as expected.This simulation is part of the set of runs with m1 = 10 6 M⊙, m2 = 10 7 M⊙, α = 2.

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

Figure 4
Figure4.The final distributions of the stars' semi-major axes (a1) and maximum eccentricities (emax, plotted as 1 − emax) for 3 runs with varying inner black hole mass (m1 = 10 5 , 10 6 , 10 7 M⊙).The critical radius for disruption used here is Rp ≤ 2RT (allowing for both partial and full disruptions), and denoted by dashed lines in the plot.Disruptions are plotted with a star-shaped marker, and plunges (disruptions that occur with RT < RS) are shaded in gray.The maximum semi-major axis included for stars in a given run is defined by the hierarchical condition (Equation1).
t ∝ M thank the participants and organizers of the fall 2023 Como meeting on the dynamics of TDEs and EMRIs for fruitful discussions and helpful clarification.B.M. is grateful for support from the U.C. Chancellor's Postdoctoral fellowship and the Mani L. Bhaumik Institute for Theoretical Physics as well as from Swift (80NSSC21K1409).D.M. is grateful for support from the NSF graduate fellowship.S.N. acknowledges the partial support from NASA ATP 80NSSC20K0505 and from NSF-AST 2206428 grant as well as thanks Howard and Astrid Preston for their generous support.E.R.R. thanks the Heising-Simons Foundation, NSF (AST-1615881 and AST-2206243), Swift (80NSSC21K1409, 80NSSC19K1391) and Chandra (22-0142) for support.The calculations for this work were carried out in part on the UCSC lux supercomputer (supported by NSF MRI grant AST-1828315).