Introduction

Galaxy mergers are a natural consequence of hierarchical structure formation (Barnes 2002;Springel & Hernquist 2005;Robertson et al. 2006;Hopkins et al. 2008). Since nearly every galaxy hosts a supermassive black hole (SMBH), binary SMBHs (SMBHBs) should be common (e.g., Sillanpaa et al. 1988;Rodriguez et al. 2006;Komossa et al. 2008;Bogdanović et al. 2009;Boroson & Lauer 2009;Dotti et al. 2009;Batcheldor et al. 2010;Deane et al. 2014;Runnoe et al. 2017;Pesce et al. 2018). Observations of merging galaxies and of dual active galactic nuclei (AGN) give credence to this idea (e.g., Komossa et al. 2004;Di Matteo et al. 2005;Bianchi et al. 2008;Comerford et al. 2009Comerford et al. , 2018;;Greene et al. 2010;Liu et al. 2010;Smith et al. 2010;Li et al. 2020;Stemo et al. 2021), and it has even been suggested that our Galaxy might host a small companion to our own SMBH (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).

After two galaxies merge, the central SMBHs will eventually sink to the bottom of the merged galaxy's potential, where they will form a binary system in the nucleus of the galaxy. While the two SMBHs coexist, they will induce orbital perturbations in each other's nuclear star clusters. These perturbations can send stars on high-eccentricity, nearly radial orbits toward one of the black holes and lead to increased rates of tidal disruption events (TDEs) (e.g., Ivanov et al. 2005;Chen et al. 2011;Wegg & Nate Bode 2011;Li et al. 2015Li et al. , 2019;;Fragione & Leigh 2018).

TDEs occur when a star passes close enough to a black hole that the tidal force from the 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 2019;Mockler et al. 2022). Theoretical TDE 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 an SMBHB, three-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, SMBHB systems are very difficult to uncover observationally. If one or both of the black holes are 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 that of 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 that the velocity shifts in the broad-line regions (BLRs) are measurable, 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 an SMBHB system; however, resolving these systems requires the SMBHs to still be at large separations (before the system forms a bound binary) and also to be 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;Foord et al. 2020;Kim et al. 2020;Reines et al. 2020;Monageng et al. 2021;Severgnini et al. 2021;Shen et al. 2021;Bhattacharya et al. 2023). Additionally, even though the merger process can funnel gas to the center of the 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 2019;Ryu et al. 2020;Wen et al. 2020;Ramsden et al. 2022). 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. 2023). Therefore, if TDE light-curve measurements suggest a much smaller black hole mass than do scaling relations such as M-σ or the bulgemass relation, this is an indication that the galaxy might host an SMBHB. 5 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 e ), this suggests the existence of a second, smaller SMBH in the galactic nucleus (e.g., Coughlin & Armitage 2018).

In this paper, we study how hierarchical three-body interactions between an SMBHB and the stars orbiting the smaller black hole can excite the eccentricities of the stars and increase the rates of TDEs 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 has been shown that SMBHBs and the EKL mechanism can lead to increased rates of TDEs. However, most previous work on this subject has 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 (Mazzolari et al. 2022;Naoz et al. 2022). Recent work by Li et al. (2019) as well as previous work by Li et al. (2015) does study TDEs around the smaller black hole in an SMBHB; however the first study focuses on lower mass ratios and uses N-body simulations,foot_1 while the second study only models one SMBHB system. 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 stars around the smaller SMBH, largely because 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 SMBHBs.

We propose that every SMBHB system undergoes a burst (a short, enhanced period) of TDEs from the smaller black hole once the binary separation is of the order of the sphere of influence (so long as one of the SMBHs is below the Hills mass for the star, ≈10 8 M e ). 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).

Methodology

Physical Processes

We consider an SMBHB m 1 < m 2 , separated by a semimajor axis a bin , and having 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 being 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 EKL mechanism (e.g., Naoz 2016).

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 prefactor in front of the octupole expansion term of the Hamiltonian (e.g., Naoz et al. 2013). 7 This defines a "hierarchical radius" limit:

1 1 . 2 hier bin bin 2 bin

The EKL mechanism describes the coherent perturbations from m 2 on the stars orbiting m 1 . These perturbations exchange angular momentum between the outer SMBHB (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., Katz et al. 2011;Lithwick & Naoz 2011;Naoz et al. 2011Naoz et al. , 2013;;Li et al. 2014;Will 2017). Importantly, the eccentricity excitations can lead to the production of a large number of TDEs (e.g., Li et al. 2015;Fragione & Leigh 2018). To model these systems and determine the rate of TDEs, 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 (( ) 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 the 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 (GR) precession (due to the potential of the disrupting black hole) and Newtonian precession (due to the potential of the stellar cusp). The timescale for GR precession is

where c is the speed of light. The Newtonian precession timescale is given by

where ψ is the true anomaly of the inner orbit,

1 cos 2 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 general relativity 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., Chen et al. 2011;Li et al. 2015).

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 the Appendix). Like Li et al. (2015), we find that t GR,1 and t NT are the two shortest timescales apart from t quad , implying that GR 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 size comparable to that of the outer black hole (previous work focusing on the larger black hole has found that the majority of disruptions are caused by chaotic orbital crossings, not by the EKL mechanism, e.g., Chen et al. 2011;Wegg & Nate Bode 2011).

Orbital precession due to general relativity and the stellar cusp can be important at small radii; however the strength of EKL perturbations increases with a star's radius from m 1 as it moves 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 choose to isolate the effects of the EKL mechanism and use the parameter space exploration described by Figure 3 to choose systems to model where the EKL timescale t quad 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 Newtonian precession from the stellar cusp in our models, and instead choose our setups such that the timescale for Newtonian 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 can affect the buildup 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 choose a variety of parameter combinations for our simulations that we expect to produce observable TDEs. The parameters we vary are 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 is set to a range of reasonable values for TDEs: 10 5 , 10 6 , and 10 7 solar masses. Black holes above ≈10 8 M e will swallow most stars whole, and disrupt them within their event horizons. Black holes below 10 5 M e 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 (e.g., Ramsden et al. 2022;Hammerstein et al. 2023; with the exception of that in 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 these determine 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 the semimajor axis, with the integrated stellar mass normalization fixed by the M-σ relation. The constants M 0 = 1.3 × 10 8 M e , σ 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 and 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).

We set the eccentricity of the SMBHB to an intermediate eccentricity of e bin = 0.5. There is debate in the community over the expected eccentricity of SMBHBs at this stage of evolution, as the order-of-magnitude difference in size scales between galaxy mergers and binary formation makes for a difficult computational challenge. For example, recent work by Vasiliev et al. (2022) studying the initial inspiral of a satellite galaxy (up to when the satellite's mass drops to 10 -3 × m initial , and therefore still well outside of the radius range studied here) found that the eccentricity evolution can depend on the mass ratio and also on the density profile of the galaxies. They found that less equal mass ratios and steeper density profiles lead to The white stars represent the parameter combinations chosen for numerical simulation runs. This is calculated by determining the range of semimajor axes for each parameter combination where the EKL timescale is the shortest dynamical timescale. All of our systems have e bin = 0.5 and a bin = 0.5r h,disrupting . While we compare against a variety of dynamical timescales (including those for two-body relaxation, vector resonant relaxation, GR precession, and binary hardening through gravitational wave emission; see Appendix), the two competing timescales are t GR,1 and t NT . These are the timescales for GR precession due to the disrupting black hole and Newtonian precession due to the stellar potential, respectively. more circularization in this initial inspiral phase. Recent work by Gualandris et al. (2022) studied the following phases of evolution, when the SMBHs at the centers of the galaxies become bound and start to harden. They also found that steeper density profiles lead to more circularization of the orbits; however, unlike the Vasiliev et al. (2022) paper, they found that less equal mass ratio systems end up with higher eccentricities than systems with more equal mass ratios.

If we use the predictions from Gualandris et al. (2022) for the radius range relevant to our systems (which are within the dynamical friction radius and outside of the hard-binary separation as defined in that work), we will find that the eccentricities predicted for the mass ratios and density profiles studied in this work range from ∼0.2 to 0.8. As long as e bin is high enough for the EKL mechanism to be effective and low enough that there are a significant number of stars within the hierarchical radius of the secondary black hole, we do not expect our results to change dramatically (see Melchor et al. 2023 for a follow-up analysis that varies the eccentricity of the SMBHB). The EKL mechanism can produce chaotic behavior even when the octupole expansion prefactor ò < 10 -2 (Katz et al. 2011). For the same separations used in this paper, e bin must decrease to 0.07 for ò to drop to 10 -2 at a hier . On the other end, the number of stars within the hierarchical radius of the black hole goes as ( )

1 3 2 (integrating the density profile), and so taking α = 1.5 as an example, e bin would have to increase to 0.9 for there to be a factor of ∼10 fewer stars available for disruption.

We set the semimajor axis of the binary 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 e )-a regime where the EKL mechanism is highly effective; however the SMBHB'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 SMBHB's hardening timescale is mainly dependent on the effectiveness of loss-cone diffusion, which is likely 100 Myr (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 run 1000 numerical simulations of the three-body interactions with different initial conditions for the star. The priors for the star's initial conditions are determined as follows: The semimajor axis of the star is drawn randomly from the density distribution, with the maximum radius determined by the hierarchical condition of the system (see Equation (1)). The star's eccentricity is drawn from an isothermal distribution, and the mutual inclination is drawn from an isotropic distribution (uniform in cos(i)). The EKL mechanism is mostly insensitive to the test particle eccentricity, and is capable of exciting stars to very high eccentricities even if they start on very low eccentricity orbits (Katz et al. 2011;Lithwick & Naoz 2011;Li et al. 2014). The argument of periapsis and the longitude of the ascending nodes are drawn from a uniform distribution between 0 and 2π. See Table 1 for an overview of the simulation parameters.

The simulations are run until either the star is disrupted or the system is evolved for 10 9 yr (by which point the TDE rate from the EKL mechanism has long since dropped off; see Figure 5). Our criterion for a tidal disruption to occur is that the star passes within twice the tidal radius (R p = 2R T ) of the black hole,

This means that our minimum impact parameter for disruption is β = R T /R p = 0.5. Setting our TDE radius condition to 2R T instead of R T ensures that we capture partial disruptions as well as full disruptions. 8 Each set of 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 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 Schwarzschild radii of the black holes (shaded 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 Note. See Section 2.3 for a description of how upper and lower limits are calculated for the rightmost three columns.

8 Below β ≈ 0.5 very little mass will return to the black hole and it is unlikely the TDE will be observable (e.g., Guillochon & Ramirez-Ruiz 2013).

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 the 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:

. 9

Hill bin bin 1 2 1 3

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).

Figure 4. The final distributions of the stars' semimajor axes (a 1 ) and maximum eccentricities (e max , plotted ase 1 max ) for three runs with varying inner black hole mass (m 1 = 10 5 , 10 6 , 10 7 M e ). The critical radius for disruption used here is R p 2R T (allowing for both partial and full disruptions), denoted by dashed lines in the plot. Disruptions are plotted with a star-shaped marker, and plunges (disruptions that occur with R T < R S ) are shaded gray. The maximum semimajor axis included for stars in a given run is defined by the hierarchical condition (Equation (1)).

Figure 5. TDE rates calculated for the EKL mechanism compared to rates calculated for two-body relaxation and to observed rates. The two-body relaxation rates (from Stone & Metzger 2016) plotted for "cuspy" galaxies (dashed lines) correspond to galaxies with an inner 2D surface brightness profile exponent of γ > 0.5 (I(R) ∝ R -γ ), while the rates plotted for "core" galaxies (dotted lines) correspond to γ < 0.3 (Lauer et al. 2007; the cuspy galaxies approximately correspond to density profile exponents α between 1 and 2; however note that these profiles are estimated observationally from radii >r h , and therefore are not a direct comparison to the profiles used in simulations in this paper for radii interior to r h ). The observed rate for all TDEs (shaded gray) is from van Velzen et al. (2020), and the observed rate in poststarburst galaxies (hatched region) is from French et al. (2020).

It is also possible for stars to remain "Lagrange-stable" even when they are no longer "Liapunov-stable" (e.g., Hayashi et al. 2022), meaning that stars can perform chaotic figure-8 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 them 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 2 of the Hill radius in our simulations with q = 10 (and within a factor of 4 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 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 TDEs.

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 yr but remain higher than the global observed rate (10 -4 TDEs/galaxy/yr) for ∼10 6 -10 8 yr. 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., when comparing our timescales for 10 6 M e secondaries to those of Chen et al. 2011 for 10 7 M e primaries). The burst is limited by the EKL 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 EKL timescale (which is generally of the order of 10 4 -10 6 yr 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 et al. (2011) found in their N-body simulations, which included binary hardening due to stellar scattering (and what Mazzolari et al. 2022 found 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 noted that loss-cone diffusion processes and gravitational wave emission were not included, and would likely have resulted in further hardening at later times and a longer-lived plateau in the rate).

TDE Fraction

The fraction of stars that become TDEs depends on the number of stars within the hierarchical limit and on 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 does not noticeably increase 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 quad , 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 within a radius r is dependent on both the Kozai-Lidov timescale and the density profile, and therefore goes as

The total fraction of stars disrupted within a given radius goes as 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

This matches the behavior we see in Table 1 when varying the density profile exponent α and holding all other parameters constant.

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. 9 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 dependent on 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 ). However, increasing the mass ratio q (by increasing m 2 ) while holding all other parameters constant also increases the strength of the EKL perturbations. It is, therefore, more difficult to analytically estimate the effect on the fraction of disruptions. If we ignore the increased perturbation strength that comes from increasing m 2 (which would only increase the lower limit), we can approximate the lower limit as purely proportional to the number of stars within the Hill radius at radii favorable for disruption:

( ) µ= aa

--N r r r r disrupt Hill Hill 3 2 Hill 3 Hill 9 2

. Combined with the fact that r Hill ∝ q -1/3 , the lower limit on the TDE fraction goes as 10

This predicts 10× fewer disruptions for q = 100 than for q = 10 for α = 1.5, and indeed we find that the lower limit on our TDE fraction drops from 9% for q = 10 to 1% for q = 100. We note that we find similar qualitative results to those 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 Hill ). 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 three-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 compare our results with those of Li et al. (2015), as they used a similar setup to our own. In their simulation, m 1 = 10 7 M e , m 2 = 10 8 M e , α = 1.75, e 2 = 0.5, and a 2 = 0.5 pc, and they found 5.7% of stars within the Hill radius are disrupted (using the more lenient criterion 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 are run with similar parameters and which find that ≈5% of stars within the Hill radius are disrupted. To compare directly, we run one additional test simulation with the same parameters (and disruption criteria) as those of Li et al. (2015), and find that 4.2% of stars within the Hill radius are disrupted. We note that Figure 6. The fraction of stars that are disrupted as TDEs as a function of the semimajor axis (scaled by the hierarchical radius). Each panel varies one simulation parameter. Top: α = 1.5, q = 10. Middle: m 1 = 10 6 M e , q = 10. Bottom: m 1 = 10 6 M e , α = 1.5. 9 We note that if we do not choose simulation parameters such that t quad < t NT and t quad < t GR,1 for the majority of stars, the steepness of the stellar cusp and the mass of the black holes might have a more noticeable effect on the fraction of stars that become TDEs. 10 We note that the TDE fractions are scaled by the number of stars within r hier , and r hier is not explicitly dependent on the mass ratio q (if a bin is dependent on q then r hier would be implicitly dependent on q, but this is not the case in our current setup). Therefore, changing q only changes the numerator of the fraction Li et al. (2015) included 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 a further explanation).

We also compare our results to those of the three-body scattering simulations of the EKL 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 the smaller SMBH on stars surrounding the 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 SMBHBs produces a burst of TDEs at rates much higher than those expected from two-body relaxation (see Figure 5). This phenomenon has been previously explored for TDEs, 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 means that the majority of disruptions do not actually come from EKL 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 EKL oscillations in particular are important, we follow Li et al. (2015)ʼs lead and explore 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 of 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 applicable only 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 SMBHBs 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 SMBHBs 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 e (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, Figure 7. Left: the number of TDEs as a function of m 1 , 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 fractions of TDEs for the α = 1 and 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 m 2 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 twobody relaxation over 10 8 yr for a single-SMBH system with mass m 1 (from Stone & Metzger 2016). Right: the number of TDEs as a function of the number of stars within the hierarchical radius. Each m 1 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 left panel.

2012 found a major-merger fraction of up to ≈25%-40% at z = 3, which is ≈10 10 yr of 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 this work we explore the rates around SMBHBs from minor mergers with mass ratios of q = 10 and q = 100. These are theoretically expected to be more common than major mergers (e.g., by factors of ∼2 and ∼6 at z = 0.1 for mergers of q = 10 and q = 100, respectively; Rodriguez-Gomez et al. 2015), although this is observationally difficult to test except for very massive galaxies (e.g., López-Sanjuan et al. 2012 found that for galaxies with M * > 10 11 M e , the minor-merger rate is ∼2-5× higher than the major-merger rate at low z, but appears relatively flat with increasing redshift). It is predicted that SMBHBs with less equal mass ratios take longer to harden than SMBHBs with more equal mass ratios; however the exact timescales vary dramatically with the model (e.g., with how dynamical friction is treated; Dosopoulou & Antonini 2017;Kelley et al. 2017).

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 short-lived. 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 but the most massive systems (see Figure 5). Therefore, it is possible that even in galaxies with close SMBHBs, 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 SMBHB, 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 highermass 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 Myr for a disrupting black hole mass of 10 6 (perturbing M h = 10 7 M e ). 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 them to the two-body relaxation rates. They found that the average rate from SMBHBs approaches the rate from two-body relaxation for black hole masses 10 7.5 M e . 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 they 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 approaches the rate from two-body relaxation for primary SMBH masses around 10 7 M e . 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 e ).

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 SMBHBs 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) pointed out, if the starburst is triggered by a merger, and the increased rates are from the resulting SMBHB, 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 SMBHBs to harden to distances of ≈1 pc are thought to be of the order of 10 9 yr (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 2019;Ramsden et al. 2022). Therefore, even if the galaxies host SMBHBs, 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 SMBHBs

Despite their rarity, bursts of TDEs produced by the EKL mechanism could help us find hidden SMBHBs. Tidal features from minor mergers dissipate quickly (Conselice et al. 2000), and the gravitational potential of a galactic nucleus hosting an SMBHB with a mass ratio 10 will be dominated by the mass of the larger black hole, making it difficult to find SMBHB candidates. However, as shown here, disruptions produced by the EKL mechanism will produce a significant number of TDEs around the smaller black hole in an SMBHB, comparable to the number of 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., µ t M h 1 2 ; Mockler 2019); therefore if a very short flare occurred around what seemed to be a much larger black hole, that would be an indication of a potential SMBHB. If the mass of the larger black hole is 10 8 M e , 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:

1. The EKL mechanism produces a burst of TDEs around the smaller black hole in an 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, 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. 2. 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 the parameter space where the EKL mechanism is the dominant dynamical mechanism; see Figures 3 and 7 and Table 1). The timescale is most strongly dependent on the black hole mass and mass ratio (see Figure 5). 3. The EKL mechanism generates a larger fraction of disruptions around the smaller black hole (∼5%-45% of stars within the hierarchical radius are disrupted for q = 10) 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). 4. Disruptions around the smaller black hole can provide evidence of the existence of the SMBHB as their lightcurve timescales will likely be shorter than what is expected from the larger black hole, which dominates the potential (Mockler 2019).

We have shown that tidal disruptions from the secondary black hole in an 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 for us to resolve distinct sources.

( ) ( ) = --t a c e a e e G m m 16 9 1 1 A4 GR,int bin 3 2 bin 2 3 2 2 3 2 1 1 2 2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ( ) ( ) ò p p yy = p t Gm a m e d M r 2 cos A5 NT 1 3 1 0 1 ( ) ( ) p = W t f m M r m 2 1 A6 RR,v vrr 1 ( ) s r = L t G m 0.34 ln A7 rel 3 2 ( ) ( ) = t a c e G m s 1 2 A8 LT 3 3 2 3 2 2 1 2 ( ) ( ) = + t a c G m m m m 4 5 64 A9 GW bin 4 5 3 1 2 1 2

The equations above represent timescales for the following physical processes: t quad is the quadrupole EKL timescale (Naoz et al. 2013;Antognini 2015); t GR,1 and t GR,2 are the firstorder post-Newtonian GR precession timescales for the inner and outer black holes; and t GR,int is the timescale of the firstorder 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 SMBHB 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 (2011). In Equation (A5), ψ is the true anomaly of the inner orbit, is the spin angular momentum of the inner SMBH (Peters 1964;Kocsis & Tremaine 2011;Naoz et al. 2013).