Counting the Unseen. II. Tidal Disruption Event Rates in Nearby Galaxies with REPTiDE

When a star passes sufficiently close to a massive black hole (MBH), it will be torn apart by extreme tidal forces (J. G. Hills 1975). After such an event, conservation of angular momentum will lead to the eventual circularization of the bound stellar material (M. J. Rees 1988). Although the underlying hydrodynamics of the circularization process are complex (K. Hayasaki et al. 2013, 2016; J. Guillochon et al. 2014; H. Shiokawa et al. 2015; C. Bonnerot & W. Lu 2020; C. Bonnerot & N. C. Stone 2021; Z. L. Andalman et al. 2022; E. Steinberg & N. C. Stone 2024), observations indicate that the end result is a luminous, multiwavelength flare (N. Bade et al. 1996; S. Gezari et al. 2006; J. S. Bloom et al. 2011; S. van Velzen et al. 2011, 2016, 2020; S. Gezari et al. 2012; I. Arcavi et al. 2014; T. W. S. Holoien et al. 2016; S. Gezari 2021) powered by a combination of shock dissipation (T. Piran et al. 2015) and MBH accretion (M. J. Rees 1988; A. Loeb & A. Ulmer 1997; A. Ulmer 1999; B. D. Metzger & N. C. Stone 2016). The resulting electromagnetic transient is referred to as a tidal disruption event (TDE). Individually, TDEs are excellent probes of accretion physics as they present the rare opportunity to witness the formation and dissipation of an accretion disk around an MBH (e.g., S. Gezari 2021), sampling the poorly understood regime of super-Eddington accretion at early times (A. Ulmer 1999; G. Lodato & E. M. Rossi 2011; L. Dai et al. 2018) while exploring old questions in accretion disk stability at late times (R.-F. Shen & C. D. Matzner 2014; S. van Velzen et al. 2019; K. Kaur et al. 2023). The ability to “light up” previously invisible MBHs makes these events perfect for detecting new MBHs and potentially weighing them through modeling of their light curves (e.g., B. Mockler et al. 2019; T. Ryu et al. 2020a; S. Wen et al. 2020, 2022; A. Mummery et al. 2023).

Despite allowing for the direct detection of distant and dormant MBHs, observed rates suggest just one TDE every 10⁴⁻⁵ yr per galaxy on average (e.g., J. L. Donley et al. 2002; S. van Velzen 2018; S. Sazonov et al. 2021; Y. Yao et al. 2023). This scarcity prevents us from using individual TDE detections as probes of MBH masses in any particular galaxy of interest. However, the rates of TDEs across different galaxy and MBH masses are likely correlated with both the MBH occupation fraction and mass distribution (N. C. Stone & B. D. Metzger2016); understanding these correlations is critical for enabling future TDE science. For example, one might hope to use near-future samples of ∼10³⁻⁴ TDEs (K. Bricman & A. Gomboc 2020; Y. Shvartzvald et al. 2024) to measure the abundance of intermediate-mass black holes (IMBHs; with masses M_• ≲ 10⁶ M_⊙), which can distinguish between debated high-redshift MBH formation channels (M. Volonteri 2010; K. Inayoshi et al. 2020). A practical approach to this goal requires some understanding of how TDE rates depend on MBH mass and other parameters to distinguish between otherwise degenerate scenarios (e.g., an IMBH occupation fraction of unity and a low per-galaxy IMBH TDE rate versus a low IMBH occupation fraction and a high per-galaxy IMBH TDE rate). Therefore, the future of TDEs as MBH demographic probes necessitates the study of TDE rates and the comparison of theoretical predictions to observations.

Theoretical TDE rates are most commonly estimated by solving the orbit-averaged Fokker–Planck equation, which models the diffusion of stars through angular momentum space via collisional two-body relaxation (see e.g., D. Merritt 2013a; N. C. Stone et al. 2020, for reviews). In this approach, the TDE rate is defined by how often stars are scattered into the “loss cone,” which is the region of phase space (position and velocity) describing orbits that will undergo disruption. As TDE rates can vary dramatically with changes to the nuclear stellar density profile, realistic attempts to estimate TDE rates are usually semiempirical, calibrating stellar profiles from high-resolution Hubble Space Telescope (HST) observations of nearby galactic nuclei (J. Magorrian & S. Tremaine 1999; D. Syer & A. Ulmer 1999; J. Wang & D. Merritt 2004; E. Vasiliev & D. Merritt 2013; N. C. Stone & B. D. Metzger 2016; H. Pfister et al. 2020), although other, more theoretical, efforts exist too (M. Polkas et al. 2024).

As outlined in C. H. Hannah et al. (2024, hereafter Paper I), the overall goal of our multipaper project is to constrain MBH demographics in low-mass galaxies by comparing TDE rate predictions for model galaxy samples to observations. While those comparisons will be presented in our upcoming Paper III, we showcase here a range of astrophysically interesting intermediate results by dynamically modeling the loss cones of the sample galaxies presented in Paper I. Given the nuclear resolution requirements of this sample, these rates represent the most accurate set of TDE rates for a sample of nearby galaxies with resolved densities on scales below r = 5 pc. This galaxy sample also covers a wide range of galaxy stellar mass and includes both early- and late-type galaxies making it more representative of the overall galaxy population than previous empirical TDE rate predictions.

In the process of computing TDE rates for our galaxy sample, we have developed an observer-friendly Python implementation of the usual steady-state loss-cone formalism called Rate Estimation Program for Tidal Disruption Events (REPTiDE), which we present here (C. H. Hannah 2025). REPTiDE was designed as an easy-to-use Python module capable of computing the loss-cone TDE rate for a given stellar mass distribution and central MBH. A practical aspect of this code is the ability to input discrete density profiles produced directly from surface-brightness data rather than assuming some parametric model. Keeping with its observer-friendly theme, REPTiDE also contains a built-in function to measure the 3D stellar density profile from a supplied surface-brightness profile and mass-to-light ratio (M/L), following a procedure identical to Paper I.

The remainder of this paper is as follows. The theoretical loss-cone methodology we use is described in Section 2. Details on the functionality, accuracy, and use of the REPTiDE software are described in Section 3. In Section 4, we present the actual TDE rates for our galaxy sample and compare them with the most recent observed TDE rates, as well as other theoretical results. Lastly, our conclusions and plans for Paper III are discussed in Section 5.

Loss cone theory was first developed half a century ago to study TDE rates (J. Frank & M. J. Rees 1976; A. P. Lightman & S. L. Shapiro 1977; H. Cohn & R. M. Kulsrud 1978), and in more recent times has seen a number of applications, not just to TDEs (J. Magorrian & S. Tremaine 1999; J. Wang & D. Merritt 2004; N. C. Stone & B. D. Metzger 2016) but also to other kinds of loss-cone phenomena such as extreme mass ratio inspirals (C. Hopman & T. Alexander 2005; L. Broggi et al. 2022; I. Qunbar & N. C. Stone 2024), hypervelocity stars (H. B. Perets et al. 2007), and MBH binary coalescence (M. Milosavljević & D. Merritt 2003). Subsequent refinements to the theory have also incorporated the role of asphericity (J. Magorrian & S. Tremaine 1999; D. Merritt & M. Y. Poon 2004; E. Vasiliev & D. Merritt 2013; E. Vasiliev et al. 2015; K. Kaur & N. C. Stone 2025), sizeable velocity anisotropies (K. Lezhnin & E. Vasiliev 2015; N. C. Stone et al. 2018), strong scatterings (A. Weissbein & R. Sari 2017; O. Teboul & H. Perets 2025; O. Teboul et al. 2024), and repeated partial disruptions (L. Broggi et al. 2024), although for now we neglect these effects in our study, as their importance is challenging to quantify through photometric data alone. Below, we review the computation of theoretical TDE rates using the simplified, stationary state formalism appropriate for relaxed, quasi-spherical galactic nuclei (J. Wang & D. Merritt2004; N. C. Stone & B. D. Metzger 2016), conditions we expect to be reasonable approximations in the lower-mass galaxies of greatest interest for our study.

For a star to be tidally disrupted, it must reach a critical radius from the MBH during its orbit. This distance is referred to as the tidal radius (r_t), and it depends⁵ on the mass (M_⋆) and radius (R_⋆) of the star as well as the mass of the black hole (M_•):

For a TDE to occur, the pericenter of the star’s orbit (r_p) must be ≤r_t. The term “loss cone” comes from the geometric interpretation of this criterion, which corresponds to a cone in the velocity phase space of stellar orbits whose members will be tidally disrupted.⁶ When considering stars on highly eccentric orbits (i.e., the most likely stars to experience tidal disruption) in a spherical gravitational potential, this condition translates to a simple cut in action space. Stars having a specific angular momentum, J, less than the loss-cone angular momentum, , will suffer tidal disruption. Near the MBH, the portion of action space with J < J_LC is devoid of stars, and TDE rates are set by the rate at which two-body scatterings can diffuse stars into the loss cone.

While this condition defines which stars will be tidally disrupted, it is agnostic as to the observability of such events. For larger black holes (or more compact stars), the disruption itself can occur within the event horizon, leaving no observable signature. The maximum MBH mass at which a star will be disrupted outside of the event horizon is referred to as the Hills mass (J. G. Hills 1975):

Here we have set r_t = 2 r_g, the event horizon size for a nonspinning black hole (here r_g = GM_•/c² is the gravitational radius). This Newtonian definition of the Hills mass is clearly crude (for example, it should be computed by using the minimum parabolic pericenter) but agrees surprisingly well with general relativistic calculations in the Schwarzschild metric; on the other hand, rapid MBH spin can increase the true Hills mass over the Newtonian value by almost an order of magnitude (M. Kesden 2012).

Giant stars can be disrupted by these more massive MBHs, but their larger tidal radii and longer timescales lead to lower observability in time-domain surveys (M. MacLeod et al. 2012, 2013). While MBH spin can increase M_H to allow for the disruption of main-sequence stars by higher-mass MBHs, these objects are, in any case, greatly outnumbered by their low-mass counterparts, which minimizes their contributions to the total TDE rate. Therefore, we do not consider any spin corrections or giant disruptions in this work.

Under the assumption of spherical symmetry for the stellar density profile ρ(r), the total potential⁷ of the system is described as

where M_enc(r) is the stellar mass enclosed within a radius of r. Assuming an isotropic velocity distribution for the stars, the stellar distribution function (DF) is given by an Eddington integral:

where is the negative of the usual specific orbital energy and is the average stellar mass of the system.

The time evolution of stellar DFs in dense environments can be understood via the collisional Boltzmann equation, which can be well approximated⁸ by a Fokker–Planck equation in which uncorrelated two-body scatterings are accounted for as velocity-space diffusion coefficients. In the spherical Kepler potential of an MBH, the local (coordinate-space) Fokker–Planck equation can be transformed into an orbit-averaged action-space Fokker–Planck equation (A. P. Lightman & S. L. Shapiro 1977). In this picture, stellar orbits diffuse through the 2D space of actions {, J} until some kind of quasi-stationary state is achieved. In the presence of a loss cone, this quasi-stationary state involves a logarithmic profile in J (H. Cohn & R. M. Kulsrud 1978).

While the loss-cone problem is thus best studied with a 2D, time-dependent diffusion equation (L. Broggi et al. 2022), the lack of detailed phase space information of distant galactic nuclei makes it challenging to apply such an approach to observations without introducing further assumptions and uncertainties. Here, we will follow J. Magorrian & S. Tremaine(1999), J. Wang & D. Merritt (2004), and N. C. Stone & B. D. Metzger (2016) in applying a simplified loss-cone formalism to observations of galactic nuclei. This formalism exploits the timescale hierarchy that exists for the eccentric orbits relevant to the loss-cone problem. For such orbits, angular momentum diffusion is orders of magnitude faster than energy diffusion, so the latter process is neglected. Furthermore, for highly eccentric orbits, loss-cone refilling is dominated by a single orbit-averaged angular momentum diffusion coefficient:

Here, where J_c() is the angular momentum of a circular orbit and r_a specifies the apocenter radius for a star with energy . The orbital period, P(), is defined as .

In Equation (5), the orbit-averaged diffusion coefficient is computed from the limit of a local angular momentum diffusion coefficient, which can be expressed as

Here, the I_n terms encode combinations of local velocity-space diffusion coefficients (see e.g., D. Merritt 2013b for more details) and can be practically computed as moments of the DF:

Additionally, is the Coulomb logarithm and is the second moment of the stellar present-day mass function (PDMF):

The diffusion coefficient is ultimately a weighted average over scatterings from a population of stars, and the contribution of different stellar mass bins is captured by . With these definitions, the flux of stellar scattering into the loss cone (i.e., the number of stars experiencing tidal disruption) per unit of time and energy can be expressed as

where R₀() defines the minimum populated angular momentum, below which no stars exist due to disruption. This can be expressed (H. Cohn & R. M. Kulsrud 1978; D. Merritt 2013a) as

where the α_m are consecutive zeros of the Bessel function J₀(α). Here, q is a dimensionless diffusivity, the ratio of the per-orbit change in to the loss-cone value :

When q ≪ 1, we are in the “empty loss-cone” regime: , and stars slowly diffuse through J-space until they are destroyed in a grazing disruption (generally; see A. Weissbein & R. Sari 2017). In this regime, the TDE rate is determined by . Conversely, when q ≫ 1, we are in the “full loss-cone” regime (also referred to as the pinhole regime): , and stars may wander in and out of the loss cone many times in a particular orbit. In this regime, the TDE rate is independent of and is determined entirely by the DF. In realistic galaxies, small radii (high ) have empty loss cones and large radii (low ) have full loss cones.

REPTiDE provides a simple, user-friendly Python implementation of this theory, which we will now discuss further in Section 3.

We present REPTiDE, a new publicly available software package⁹ that implements the loss-cone physics from Section 2 as a series of functions written in Python and computes the expected TDE rate for a spherically symmetric galactic nucleus (C. H. Hannah 2025). In this section, we describe the functionality of REPTiDE along with its key features.

The primary inputs required for this calculation are the central black hole mass (M_•) and the radial 3D stellar density profile ρ(r). There are two versions of the code that depend on the definition of the stellar density profile; this can be either “Discrete” or “Analytic.” In the discrete version, the density profile is supplied as a discrete, nonparametric set of data points. The analytic version uses a density profile that follows a power law with Sérsic decay:

Here, ρ_5pc is the stellar density at a radius of 5 pc, γ is the power-law slope, and r_d is the radius where the exponential decay starts. For the Sérsic decay component, I_eff is defined such that the density matches at r_d given the user-defined Sérsic index (n) and effective radius (r_eff).

The analytic version is designed to take inputs from the galaxy stellar mass versus density and power-law slope relations presented in Paper I. This version of input was developed for our planned application of REPTiDE to model galaxy samples in Paper III (C. H. Hannah et al. 2025, in preparation).

3.1. One-dimensional Surface-brightness Profile to Three-dimensional Density

3.2. Procedure

The calculation begins by computing the gravitational potential defined in Equation (3), which is subsequently used to define the range of specific orbital energies to consider. The boundaries on orbital energies are defined as the potential at 100 r_t (for a solar-mass star) and the potential at 10⁶ pc. These boundaries define a wide range in energy space and ensure the peak of the loss-cone flux curve is captured. The default energy grid consists of 1000 energies, but this can be altered with the “n_energies” keyword, which specifies the desired number of logarithmically spaced orbital energies. However, we do not encourage using a small value for “n_energies,” for reasons quantified later in Section 3.5.

The code then computes the stellar DF (Equation (4)) using a double exponential transform to handle the integrable singularity in the DF integral. The next step is to compute the following quantities as functions of orbital energy: orbit-averaged angular momentum diffusion coefficient (), period (P), the squared ratio of the loss-cone angular momentum to the angular momentum of a circular orbit (), and the diffusivity parameter (q). As with the DF integral, we use a double exponential transform for the edge singularity in the period integral. These quantities are first calculated for a “monochromatic” population of 1 M_⊙ stars and are later adjusted when accounting for the present-day stellar mass function (see below). With all of the above functions defined, the loss-cone flux (Equation (10)) is tabulated and integrated over all orbital energies to define the “solar” TDE rate.

Monochromatic stellar populations do not exist in reality, and thus, the final step in this calculation involves accounting for the PDMF. Implementing a realistic PDMF has two competing effects on the TDE rate. The inclusion of more subsolar-mass stars increases the TDE rate due to their large numbers, but this can also change the angular momentum diffusion coefficient (Equation (5)), altering the TDE rate. N. C. Stone & B. D. Metzger (2016) tested two different PDMFs (Salpeter and Kroupa), and both were shown to have similar rate enhancements (≈1.63 for Salpeter and ≈1.53 for Kroupa with M_⊙ and M_⊙, respectively). In their investigation, had a larger effect on the level of rate enhancement (which is directly correlated with the age of the stellar system) because the diffusion coefficient is set by the most massive stars present. REPTiDE implements the Kroupa PDMF (P. Kroupa 2001) as shown below:

where and specify the minimum and maximum stellar masses to consider. These values can be set using the “M_min” and “M_max” parameters, respectively.

To determine the “full” TDE rate, the PDMF is first defined over the specified mass range and normalized. This PDMF is then used to compute and , which are set to M_⊙ and M in the monochromatic case, respectively. Each stellar mass considered possesses a different tidal radius (r_t), and thus, some quantities have to be adjusted to compute the TDE rate per stellar mass. First, must be adjusted given the new value for and the new integral bounds (as we set r_p to r_t in the orbit-averaging calculations). We also adjust for the new value of r_t. These quantities then define the adjusted q() values needed to compute the loss-cone flux curve for a given stellar mass. This process is performed for each of the 50 stellar masses, and the integrated loss-cone flux curves give the TDE rate per stellar mass. Lastly, this TDE rate per stellar mass is multiplied by the PDMF and integrated to give the “full” TDE rate.

When computing the “full” TDE rate, REPTiDE can account for TDE rate suppression due to direct capture, which is enabled through the “EHS” flag. When set, the final integration of the TDE rate per stellar mass is only computed for stellar masses with a Hills mass less than the MBH (as defined by Equation (2)). All TDE rates presented in this work are “full” rates with minimum and maximum stellar masses of 0.08 M_⊙ and 1 M_⊙, respectively.

3.3. Handling Small Radii

In these calculations, extrapolation of the density profile to small radii (well inside the resolution of the density measurements) is sometimes required. Extrapolating very steep power-law density profiles can cause a few issues: if γ ≥ 2, diffusion of stars through energy space becomes strongly nonlocal; if γ ≥ 9/4, the TDE rate actually diverges as one considers more tightly bound stars; if γ ≥ 3, the stellar mass enclosed diverges. However, it is generally unphysical for nuclear density profiles to be characterized by power laws with very large γ values over an arbitrary range of radii (N. C. Stone et al. 2018). To overcome these issues, we have implemented the ability to set a Bahcall–Wolf cusp density profile for the innermost radii.

The Bahcall–Wolf cusp is a theoretically expected steady-state distribution of stars around an MBH formed through collisional relaxation processes (J. N. Bahcall & R. A. Wolf 1976, 1977). When γ > 3/2, relaxation times are the shortest near the MBH, so this cusp is expected to form from the inside out, resulting in a broken power-law density profile with a break radius related to the dynamical age of the system (N. C. Stone et al. 2018). Main-sequence stars in relaxed stellar cusps usually follow a power law with γ_BW ≈ 1.5–1.75, where the exact value depends on the distribution of stellar masses (T. Alexander & C. Hopman 2009; M. Preto & P. Amaro-Seoane 2010).

An approximate relation between the relaxation radius (i.e., the radial extent of the Bahcall–Wolf cusp) and the dynamical age of the system can be derived by setting the relaxation time equal to the dynamical age. For a power-law density profile described as , this relation is

where κ = 0.34 (J. Binney & S. Tremaine 2008) and is the stellar mass contained within the relaxation radius.

When applying REPTiDE to the sample galaxies from Paper I (see Section 4), extrapolation of the density profiles inward of our resolution limits is required. Thus, for all galaxies with density profiles steeper than a Bahcall–Wolf cusp (γ ≥ 1.75), we use Equation (16) to estimate the relaxation radius under the conservative assumption that the system has relaxed over a Hubble time. For a few galaxies, this assumption is inconsistent as a Bahcall–Wolf cusp is expected to be resolved. In these cases, we apply the cusp at the resolution limit. We note that the TDE rate is quite sensitive to the value of r_relax with the effect becoming largest for the steepest density profiles. This uncertainty is discussed in more detail in Section 4.2.

3.4. Accuracy

To ensure the accuracy of REPTiDE’s calculations, we perform direct comparisons with the output from N. C. Stone & B. D. Metzger (2016) and the PHASEFLOW Fokker–Planck solver (E. Vasiliev 2017), which has been used to derive TDE rates (e.g., H. Pfister et al. 2020; E. Bortolas 2022; J. N. Y. Chang et al. 2025; M. Polkas et al. 2024; Figure 1). The top panel shows a one-to-one comparison of the TDE rates computed by REPTiDE using identical density profiles and black hole masses to those used in N. C. Stone & B. D. Metzger (2016).¹⁰ The median rate residual from this investigation is 0.019 dex and is believed to be caused by minor differences in orbital energy grid densities (described further in Section 3.5). Overall, the rates agree very well. The lower panel of Figure 1 shows the loss-cone flux curve for solar-mass stars derived by PHASEFLOW and REPTiDE for an identical density distribution and MBH mass (Dehnen profile with γ = 1.2, parameters “scaleRadius” = 10, “mass” = 10¹⁰, and “Mbh” = 10⁸; see E. Vasiliev 2017, for more details). In this comparison, the maximum difference in loss-cone flux near the peak (log) is 0.09 dex, resulting in a TDE rate difference of only 0.02 dex. These results together highlight the accuracy of the computations performed by this software.

Zoom In Zoom Out Reset image size

3.5. Effect of Orbital Energy Resolution

In Paper I, we constructed a sample of 91 galaxies where 3D stellar densities could be reliably measured on parsec scales. This sample is focused primarily on nucleated galaxies (i.e., ones hosting nuclear star clusters (NSCs) at their centers) extending down to stellar masses of 10^6.58 M_⊙ and was used to construct scaling relations between NSC density structure and galaxy mass. While this sample is not a volume-limited sample, it includes a range of Hubble types and includes galaxies selected from distance-limited samples (R. Pechetti et al. 2020; N. Hoyer et al. 2023), making these galaxies representative of those in the local Universe. In this work, we present the theoretically expected TDE rate for each galaxy in this sample computed with REPTiDE.

For these calculations, we use the best available central MBH mass measurement for each galaxy along with the discrete density profiles (i.e., using the discrete version of REPTiDE). The MBH masses were derived in one of three methods (in decreasing order of reliability): (1) dynamical mass measurements (available for 30 galaxies; A. E. Reines & M. Volonteri 2015; R. C. E. van den Bosch 2016; J. E. Greene et al. 2020), (2) a galaxy morphology dependent M_BH–σ relation from J. E. Greene et al. (2020), and (3) a galaxy morphology dependent M_BH–M_gal relation from J. E. Greene et al. (2020; see Paper I for more details). We fix the slope used in the inward extrapolation of the density profiles to the 3D power-law slopes derived in Paper I, unless the nucleus is steeper than a Bahcall–Wolf cusp (γ < −1.75). As discussed in Section 3.3, these systems are not stable and should relax into such a cusp. This expectation, paired with diverging TDE rates resulting from “ultrasteep” power-law density slopes (γ < −2.25), motivates our decision to apply a Bahcall–Wolf cusp slope of −7/4 inward of our resolution limit for these galaxies (47/91 galaxies). The radial extent of the cusps in these galaxies is defined through Equation (16), assuming relaxation over a Hubble time (unless this suggests a relaxation radius larger than the resolution limit, in which case we apply the cusp at the resolution limit; this occurs in 12/47 galaxies). For large radii beyond our measurements, we apply a simple exponential decay with a width of 1 kpc.

Figure 2 shows the TDE rates for these galaxies as functions of galaxy stellar mass (left panel) and black hole mass (right panel). These purely dynamical rate predictions do not account for event horizon rate suppression (which we will return to later) and assume a Kroupa PDMF with masses spanning from 0.08 M_⊙ to 2 M_⊙, which represents an ∼1 Gyr old stellar population. While there is some uncertainty introduced by unknown nuclear stellar ages, the PDMF enhancement factor ranges from 1.5 to 3.8 for stellar ages of 10 Gyr–30 Myr (N. C. Stone & B. D. Metzger 2016). Uncertainties in the TDE rates are indicated by symbol size, with the least certain being smaller and vice versa. The symbol shapes indicate the origin of the M_• estimate. Both of these plotting schemes are used for all results figures in this paper.

Zoom In Zoom Out Reset image size

Previous TDE rate measurements from N. C. Stone & B. D. Metzger (2016) are shown as gray crosses. Although the rates from N. C. Stone & B. D. Metzger (2016) are higher on average, there is considerable overlap in these results. For the 30 galaxies belonging to both samples, all have higher M/L ratios in N. C. Stone & B. D. Metzger (2016), which is the most likely reason for the higher rates, as increasing the M/L ratio will increase the overall stellar density at all radii. It is also worth noting here that none of the 150 density profiles used in N. C. Stone & B. D. Metzger (2016) require Bahcall–Wolf cusps at small radii (γ ≥ 2.25, of which we have 26), which further highlights the importance of resolving NSC scales when computing TDE rates. This is likely due to the use in N. C. Stone & B. D. Metzger (2016) of parameterized “Nuker-law” surface-brightness profiles (T. R. Lauer et al. 2005, 2007) that were fit over a relatively large range of spatial scales and were therefore less sensitive to steep slopes produced by NSCs.

To characterize these relations and aid future comparisons, we fit broken power laws to the data weighted by the TDE rate uncertainties, which are shown as solid black lines in Figure 2 and described as

Here, M_BH/gal represents the MBH mass or galaxy stellar mass (depending on the desired relation), M_b is the break mass, α and β are the slopes, and the normalization constants (A and B) are related by. Because the MBH and galaxy mass ranges differ, we normalize each power law at a different mass specified by M_norm; the galaxy mass fit uses M_norm = 10⁹ M_⊙ while the MBH mass fit uses M_norm = 10⁶ M_⊙. For TDE rates as a function of galaxy mass, we find , α = 1.19 ± 0.22, β = −0.88 ± 0.38, and . Our best-fit parameters for as a function of MBH mass are, α = 1.01 ± 0.25, β = −0.54 ± 0.23, and . Note that all logarithms in this paper are base 10.

The uncertainties on fit parameters were measured via 1000 iterations of bootstrap resampling (as are all fit parameter errors presented in this work). We estimate the scatter in these fits using Gaussian deconvolution of the fit residuals assuming the residuals follow a Gaussian distribution as in C. Pryor & G. Meylan (1993). We include the TDE rate uncertainties here as the heteroscedastic errors on the fit residuals. The scatter was measured both above and below the break mass and is shown as the gray shaded regions around the best-fit broken power law; at low masses, the scatter is 0.98 and 0.91 dex (M_gal/M_•), while at higher mass, it is 0.52 and 0.53 dex.

The black dashed line in Figure 2 shows the maximum black hole mass given a constant TDE rate over a Hubble time, assuming each TDE contributes . If a galaxy were to lie above this line, its position would indicate that it could not have maintained its present-day mass and TDE rate for a Hubble time (as otherwise stellar consumption would have dramatically increased the mass of the MBH, moving its position below/to the right of this line). This line can be viewed as a crude self-consistency check in our TDE rate calculation, at least at the ensemble level, however, one can always invoke an unusual evolutionary history for any individual outlier.

Almost all galaxies in our sample do indeed fall below this line, suggesting that the uncertainties in our modeling (most notably, the extrapolation of galaxy scaling relations into the IMBH regime to estimate black hole masses) are not causing dramatic errors. The nine galaxies above this line could arise from variable TDE rates over the life of these galaxies (i.e., with lower TDE rates at earlier times), a misestimation of stellar M/L ratios, or an underestimate of M_•. The lack of a galaxy pileup just below this line hints that for most galaxies in our sample, TDEs have not dominated MBH growth (N. C. Stone et al. 2017), though reaching a firm conclusion in this regard would require more robust IMBH mass estimates.

4.1. Comparisons with the Literature

Tension has existed between observed and theoretical TDE rates for years. Previous loss-cone-based TDE rates from J. Magorrian & S. Tremaine (1999), J. Wang & D. Merritt (2004), and N. C. Stone & B. D. Metzger (2016) have overpredicted the expected TDE rates for Milky Way-like galaxies by a factor of a few to an order of magnitude when compared to rate inferences from time-domain surveys. Most past loss-cone modeling predicted a rate of ∼10⁻⁴ galaxy⁻¹ yr⁻¹, while observed rates inferred from optical and X-ray surveys (e.g., J. L. Donley et al. 2002; S. van Velzen & G. R. Farrar 2014; S. van Velzen et al. 2020; Y. Yao et al. 2023) usually suggest a rate of ∼ few 10⁻⁵ galaxy⁻¹ yr⁻¹.

The resolution of this tension is not clear at present. While most assumptions of standard loss-cone modeling are conservative (N. C. Stone & B. D. Metzger 2016), there are dynamical mechanisms to reduce TDE rates. Most notably, a tangential velocity anisotropy, arising either from the aftermath of an MBH merger (K. Lezhnin & E. Vasiliev 2015) or from loss-cone “shielding” by a steep cusp of stars and/or stellar mass black holes (O. Teboul & H. Perets 2025; O. Teboul et al. 2024), will be capable of dramatically suppressing rates relative to the quasi-isotropic cusps we consider here. On the other side of the equation, observationally inferred rates are sensitive to other uncertainties: the poorly characterized bottom end of the TDE luminosity function (S. van Velzen 2018), nuclear dust obscuration (N. Roth et al. 2021), and uncertainties in TDE emission mechanisms (L. Dai et al. 2018; S. van Velzen et al. 2021). In this work, however, we reexamine this basic tension using our larger galaxy sample, with its relatively careful treatment of central light profiles.

Figure 3 compares our TDE rates to the most recent observed rates from Y. Yao et al. (2023, light blue lines) in addition to other recent theoretical rate results. Unlike Figure 2, the TDE rates presented here account for event horizon suppression as this dramatically affects the observability of TDEs. We use the same mass functions as Y. Yao et al. (2023), namely the I. K. Baldry et al. (2012) galaxy stellar mass function and the black hole mass function from E. Gallo & A. Sesana (2019), to convert the observed volumetric rates to per-galaxy rates.

Zoom In Zoom Out Reset image size

Our fits to the distributions of TDE rates from Figure 2 (black line and gray shaded band) are not ideal for comparison to observed rates for two reasons. First is the lack of event horizon suppression at the high-mass end. Second, at a given galaxy or MBH mass, galaxies will span a wide range of TDE rates and the observed rates will be dominated by the highest TDE rate galaxies, even when subdominant in the overall population. Thus, we perform additional fits to the binned means of our rate results including event horizon rate suppression via direct capture (yellow lines show the binned means and blue lines show the fits; Figure 3). Given that the Hills mass has a clear physical origin and is separate from a break in TDE rates near the IMBH–supermassive black hole (SMBH) transition, it is important to capture both of these effects separately. As such, we fit doubly broken power laws to this data described by

Similarly to the single broken power-law fits in Figure 2, M_BH/gal specifies the MBH or galaxy mass, and M_norm gives the normalization mass (same as previous fits). Here, the normalization constants are related by and .

To avoid introducing uncertainties from uneven numbers of galaxies in fixed bin sizes (e.g., J. Maíz Apellániz & L. Úbeda 2005), we implement variable bin sizes ensuring seven galaxies per bin. Uncertainties on the TDE rates were used to perform a weighted mean in each bin with the error on the mean defined by adding the individual rate uncertainties in quadrature. We find the following best-fit parameters for the binned mean TDE rates: : , α = 0.77 ± 0.27, β = −1.38 ± 1.01, ζ = −2.34 ± 2.13, , and ; : , α = 1.08 ± 0.31, β = −0.18 ± 0.79, ζ = −8.24 ± 3.64, , and . Across the full range of galaxy masses, the left panel of Figure 3 shows our average rates agree reasonably well with observed TDE rates from the Zwicky Transient Facility (ZTF; Y. Yao et al. 2023). For the first time, our theoretical rate predictions for Milky Way-mass galaxies are consistent with observational estimates of a few × 10⁻⁵ galaxy⁻¹ yr⁻¹. For rates as a function of MBH mass (right panel), we find remarkable agreement with observations extending well into the IMBH regime.

We are not the first of recent studies to bring theoretical and observed TDE rates into better agreement, and we include these in Figure 3 for comparison. M. Polkas et al. (2024) used semianalytic galaxy formation and evolution models to measure TDE rate evolution over cosmological timescales while considering galaxy nucleation. We show the z = 0 results from these simulations as bright green dotted–dashed lines in Figure 2 (converted from volumetric to per-galaxy rates using their fiducial black hole mass function and the I. K. Baldry et al. 2012 galaxy stellar mass function).

Another recent work, J. N. Y. Chang et al. (2025), presents a similar characterization of TDE rates as a function of galaxy mass, but using a very different sample of galaxies. Rather than focusing on high-quality nuclear stellar density measurements in the nearest galaxies as we do here, this study included two samples of galaxies, (1) 23 galaxies with potential IMBHs both in galaxy nuclei and in massive star clusters, and (2) 37 spiral galaxies with well-measured MBHs from B. L. Davis et al. (2019). Their IMBH candidates include both some with secure dynamical measurements, as well as IMBH candidates with just dynamical upper limits on their masses from N. Neumayer & C. J. Walcher (2012). Thus, the black hole masses in this sample are quite uncertain (as they are in our sample as well). In many cases, the light profiles for the IMBH hosts are created based solely on a fitted NSC light profile excluding the role of the galaxy. We note that although the TDE rate is dominated by the NSCs (H. Pfister et al. 2020), this does not mean the depth of the host galaxy potential does not impact the TDE rates; we find the inclusion/exclusion of a host galaxy can impact our TDE rates by up to 1 dex. The higher-mass MBH sample in J. N. Y. Chang et al. (2025) has well-determined black hole masses, but the mass profile information comes from low-resolution data relative to our sample and does not include NSC components (B. L. Davis et al. 2019).

A novel commonality in our results and those of J. N. Y. Chang et al. (2025) is the rollover in per-galaxy TDE rates for IMBHs located in dwarf galactic nuclei. In contrast to almost all past dynamical modeling of SMBH TDE rates (which decrease with increasing SMBH mass), both J. N. Y. Chang et al. (2025) and this paper find that IMBH TDE rates increase with increasing IMBH mass. This result is particularly notable insofar as our methodology and samples are quite different as noted above. We include the average TDE rate as a function of MBH mass from J. N. Y. Chang et al. (2025) as a green dashed line in the right panel of Figure 3. Compared to our binned mean fit (blue line) we find a similar shape, but with a significantly offset peak MBH mass, as this was not a fitted parameter in J. N. Y. Chang et al. (2025). Regardless of these differences, the qualitatively similar conclusions reached in both papers are likely robust to their largest uncertainty (choice of IMBH mass), and the reduced TDE rate in dwarf galactic nuclei likely represents lower density conditions in the NSCs there (R. Pechetti et al. 2020; Paper I). We note that low TDE rates in IMBH nuclei can also be understood from the dashed line in Figure 2: if dwarf galactic nuclei had TDE rates , smaller IMBHs would quickly grow into small SMBHs through star capture alone.

Overall, theoretical TDE rates appear to be getting closer to observations. This tentative agreement should be investigated further with new observational TDE rates based on larger samples of TDEs, as opposed to the sample of 33 TDEs used in Y. Yao et al. (2023). In addition, the previously unobserved class of infrared TDEs found by M. Masterson et al. (2024) shows a rate similar to optical and X-ray TDEs, suggesting that the existing rate estimates from a single wavelength band may underestimate the total TDE rate.

4.2. Rate Uncertainties

To determine the uncertainties on the TDE rates, we explore four primary sources of error: MBH mass, density normalization (largely set by M/L uncertainties), the power-law density slope used in extrapolation, and the radial extent of the Bahcall–Wolf cusp, if applied. For the first three sources of error, we perform three rate calculations for each galaxy with the following parameter values: the nominal value, nominal − 1σ error, and nominal + 1σ error. The errors on MBH masses come from their respective sources, while the uncertainties on the density normalization and slope were derived in Paper I. We then take the standard deviation of each set of three rates to measure the uncertainty originating from each source.

For steep galaxies with a Bahcall–Wolf cusp, we estimate the uncertainty through two rate calculations: one using our nominal definition of the Bahcall–Wolf radius based on relaxation over a Hubble time (as discussed in Section 3.3) and another using the observed resolution limit of each galaxy. We define the resulting rate uncertainty as half of the difference in TDE rates between the two cases.

Lastly, these individual errors were added in quadrature to define each galaxy’s TDE rate uncertainty (listed in Table 1). The individual uncertainties are shown in Figure 4 as a function of MBH mass where the symbol colors indicate the dominant source of error. Clearly, most TDE rate uncertainties are driven by the error in the MBH mass measurements. However, some of the largest uncertainties are driven by unknown dynamical ages, which define the radial extent of the Bahcall–Wolf cusps applied in galaxies with steep density profiles. The largest TDE rate uncertainty is related to this effect and belongs to NGC 4592. The resolution limit for this galaxy is 2.80 pc while relaxation over a Hubble time infers a Bahcall–Wolf cusp radius of 0.09 pc. Although this is not the largest difference in Bahcall–Wolf cusp radius between the definitions, it is the steepest galaxy (γ = 2.99) with such a difference, which leads to the high rate uncertainty.

Table 1. Tabulated Tidal Disruption Event Rate Results for all 91 Galaxies

Name	Type	log(Mgal)	log(MBH)	Source	γ	log(ρ5 pc)			BW Cusp	rBW
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
BTS 076	Late	7.51		5	−2.41	0.41	−7.66 ± 0.21	−7.66 ± 0.21	True	2.14e+00
BTS 109	Late	6.58		5	−2.97	0.15	−6.4 ± 0.57	−6.41 ± 0.57	True	6.92e–01
DDO 084	Late	8.53		5	−2.9	1.1	−6.14 ± 0.39	−6.14 ± 0.39	True	9.50e–01
ESO 274-1	Late	9.12		5	−2.71	2.3	−4.74 ± 0.29	−4.74 ± 0.29	True	4.73e–01
IC 5052	Late	9.18		5	−2.42	2.4	−5.16 ± 0.28	−5.16 ± 0.28	True	9.33e–01
LeG 09	Early	6.99		5	−2.15	1.24	−6.22 ± 0.23	−6.22 ± 0.23	True	1.73e+00
NGC 0584	Early	10.56		3	−1.38	3.74	⋯	−4.98 ± 0.07	False	⋯
NGC 0596	Early	10.46		4	−1.45	3.75	−5.83 ± 0.16	−4.89 ± 0.16	False	⋯
NGC 0821	Early	10.7		1	−1.36	4.03	⋯	−4.78 ± 0.10	False	⋯
NGC 1374	Early	10.3		3	−0.4	3.64	⋯	−5.27 ± 0.34	False	⋯
NGC 1399	Early	11.16		3	−1.12	2.63	⋯	−5.47 ± 0.08	False	⋯
NGC 1427	Early	10.39		4	−2.24	3.7	−6.0 ± 0.52	−4.71 ± 0.52	True	6.82e–03
NGC 1439	Early	10.49		4	−2.23	4.02	−4.99 ± 0.56	−4.16 ± 0.56	True	2.89e–02
NGC 2300	Early	10.81		4	−0.82	2.47	⋯	−5.68 ± 0.19	False	⋯
NGC 2434	Early	10.58		4	−1.15	3.98	⋯	−4.96 ± 0.22	False	⋯
NGC 2778	Early	9.94		3	−2.05	4.29	−3.33 ± 0.40	−3.33 ± 0.40	True	5.30e–01
NGC 2787	Late	9.95		3	−1.31	3.7	−5.24 ± 0.18	−4.9 ± 0.18	False	⋯
NGC 2903	Late	10.4		1	−1.65	3.45	−5.0 ± 0.11	−5.01 ± 0.11	False	⋯
NGC 3115B	Early	9.03		4	−1.85	2.9	−4.5 ± 0.20	−4.5 ± 0.20	True	1.65e+00
NGC 3274	Late	8.22		5	−2.84	1.92	−5.07 ± 0.24	−5.07 ± 0.24	True	1.18e+00
NGC 3344	Late	9.9		4	−1.62	3.49	−4.69 ± 0.12	−4.69 ± 0.12	False	⋯
NGC 3379	Early	10.63		3	−1.08	3.51	⋯	−4.9 ± 0.07	False	⋯
NGC 3384	Early	10.13		3	−1.98	3.75	−4.12 ± 0.29	−4.12 ± 0.29	True	1.12e–01
NGC 3412	Early	10.01		4	−2.23	2.87	−5.66 ± 0.65	−5.61 ± 0.65	True	9.28e–03
NGC 3522	Early	9.69		4	−1.94	3.73	−4.17 ± 0.40	−4.17 ± 0.40	True	4.41e–02
NGC 3585	Early	10.93		3	−1.23	3.8	⋯	−4.94 ± 0.18	False	⋯
NGC 3607	Early	10.87		3	−0.95	3.48	⋯	−5.0 ± 0.03	False	⋯
NGC 3608	Early	10.45		3	−0.98	3.74	⋯	−4.88 ± 0.20	False	⋯
NGC 3610	Early	10.8		4	−1.77	4.05	⋯	−4.72 ± 0.22	True	1.27e–07
NGC 3640	Early	10.9		3	−0.23	3.56	−5.65 ± 0.04	−4.76 ± 0.04	False	⋯
NGC 3945	Early	10.54		3	−2.4	4.11	−3.35 ± 0.33	−3.35 ± 0.33	True	7.80e–01
NGC 4026	Early	10.06		3	−1.67	3.93	⋯	−5.07 ± 0.35	False	⋯
NGC 4242	Late	9.1		5	−3.07	1.67	−5.51 ± 0.38	−5.51 ± 0.38	True	1.20e+00
NGC 4262	Early	9.97		4	−1.32	4.17	⋯	−4.82 ± 0.21	False	⋯
NGC 4291	Early	10.5		3	−0.56	2.91	⋯	−5.45 ± 0.26	False	⋯
NGC 4342	Early	9.93		3	−1.99	4.48	⋯	−4.35 ± 0.33	True	6.26e–05
NGC 4365	Early	11.09		4	−0.5	2.97	⋯	−5.28 ± 0.08	False	⋯
NGC 4377	Early	10.02		4	−2.29	3.16	−5.64 ± 0.77	−5.28 ± 0.77	True	9.28e–03
NGC 4379	Early	9.92		4	−2.32	3.96	−4.1 ± 0.50	−3.92 ± 0.50	True	1.85e–01
NGC 4382	Early	11.02		3	−0.57	2.99	−5.18 ± 0.25	−5.18 ± 0.25	False	⋯
NGC 4387	Early	9.84		4	−2.15	3.75	−4.28 ± 0.40	−4.28 ± 0.40	True	1.68e–01
NGC 4434	Early	10.03		3	−1.87	4.17	−5.11 ± 0.18	−4.31 ± 0.18	True	3.07e–04
NGC 4458	Early	9.73		4	−0.62	3.79	−4.53 ± 0.08	−4.53 ± 0.08	False	⋯
NGC 4464	Early	9.58		4	−1.53	3.7	−5.19 ± 0.21	−4.91 ± 0.21	False	⋯
NGC 4467	Early	9.04		4	−1.33	3.58	−4.73 ± 0.04	−4.73 ± 0.04	False	⋯
NGC 4472	Early	11.3		3	−1.46	2.76	⋯	−5.31 ± 0.05	False	⋯
NGC 4473	Early	10.49		3	−0.81	3.52	−6.16 ± 0.11	−4.96 ± 0.11	False	⋯
NGC 4474	Early	9.97		4	−2.18	3.61	−4.1 ± 0.42	−4.1 ± 0.42	True	3.08e–01
NGC 4478	Early	10.08		4	−1.79	3.61	−5.57 ± 0.13	−4.89 ± 0.13	True	2.67e–07
NGC 4483	Early	9.72		4	−1.66	3.87	−4.44 ± 0.19	−4.44 ± 0.19	False	⋯
NGC 4486B	Early	9.33		3	−1.0	3.46	⋯	−6.59 ± 0.55	False	⋯
NGC 4489	Early	9.65		4	−1.35	3.42	−4.77 ± 0.10	−4.77 ± 0.10	False	⋯
NGC 4494	Early	10.68		4	−1.44	4.24	−5.36 ± 0.17	−4.39 ± 0.17	False	⋯
NGC 4517	Late	9.86		4	−2.72	2.12	−5.33 ± 0.42	−5.33 ± 0.42	True	7.85e–01
NGC 4528	Early	9.85		4	−1.55	3.71	−4.69 ± 0.13	−4.69 ± 0.13	False	⋯
NGC 4551	Early	9.86		4	−1.88	3.62	−4.68 ± 0.28	−4.68 ± 0.28	True	4.02e–03
NGC 4552	Early	10.7		3	−2.15	3.52	⋯	−4.77 ± 0.04	True	3.42e–05
NGC 4589	Early	10.72		4	−0.79	3.33	⋯	−5.61 ± 0.25	False	⋯
NGC 4592	Late	9.45		5	−2.99	2.33	−2.89 ± 1.42	−2.9 ± 1.42	True	8.74e–02
NGC 4600	Early	9.01		4	−2.02	2.83	−5.31 ± 0.37	−5.31 ± 0.37	True	4.65e–02
NGC 4605	Late	9.17		4	−2.9	2.97	−3.0 ± 0.38	−3.0 ± 0.38	True	4.18e–01
NGC 4612	Early	9.97		4	−2.53	3.7	−3.81 ± 0.41	−3.81 ± 0.41	True	7.64e–01
NGC 4623	Early	9.71		4	−2.85	3.58	−3.08 ± 0.73	−3.08 ± 0.73	True	3.87e–01
NGC 4638	Early	10.14		4	−1.36	4.37	−4.3 ± 0.13	−4.09 ± 0.13	False	⋯
NGC 4649	Early	11.18		3	−1.47	3.02	⋯	−5.44 ± 0.12	False	⋯
NGC 4660	Early	9.93		4	−2.65	4.57	⋯	−3.15 ± 0.86	True	7.15e–02
NGC 4733	Early	9.83		4	−3.98	3.24	−4.24 ± 0.20	−4.24 ± 0.20	True	2.95e+00
NGC 5011C	Early	7.52		5	−2.63	0.98	−5.91 ± 0.24	−5.91 ± 0.24	True	6.33e–01
NGC 5055	Late	10.5		1	−1.66	3.58	⋯	−6.67 ± 0.38	False	⋯
NGC 5068	Late	9.39		5	−1.84	2.27	−5.93 ± 0.23	−5.93 ± 0.23	True	7.44e–02
NGC 5195	Early	10.1		4	−1.77	4.09	−4.65 ± 0.25	−4.36 ± 0.25	True	1.03e–04
NGC 5236	Late	10.44		5	−2.54	3.95	−3.28 ± 0.57	−3.28 ± 0.57	True	7.06e–01
NGC 5238	Late	7.75		5	−2.94	1.59	−4.69 ± 0.25	−4.69 ± 0.25	True	5.90e–01
NGC 5457	Late	10.21		1	−1.81	3.05	−5.06 ± 0.07	−5.06 ± 0.07	True	5.33e–02
NGC 5557	Early	10.88		4	−0.72	3.39	⋯	−5.05 ± 0.21	False	⋯
NGC 5576	Early	10.58		3	−1.11	4.0	⋯	−4.71 ± 0.14	False	⋯
NGC 5813	Early	11.08		3	−0.01	2.37	⋯	−5.13 ± 0.01	False	⋯
NGC 5982	Early	10.97		4	−0.53	2.68	⋯	−5.31 ± 0.13	False	⋯
NGC 6503	Late	9.64		1	−2.32	3.06	−4.36 ± 0.06	−4.36 ± 0.06	True	9.94e–01
NGC 7457	Early	9.88		3	−2.58	3.79	−3.52 ± 0.38	−3.52 ± 0.38	True	6.27e–01
NGC 7713	Late	8.97		5	−3.36	1.5	−5.65 ± 0.34	−5.65 ± 0.34	True	1.32e+00
NGC 7727	Late	10.67		4	−1.29	3.68	−5.1 ± 0.13	−4.79 ± 0.13	False	⋯
PGC 4310323	Late	6.61		5	−2.76	0.87	−5.98 ± 0.21	−5.98 ± 0.21	True	1.09e+00
UGC 07242	Late	7.75		5	−2.88	1.05	−6.01 ± 0.22	−6.01 ± 0.22	True	9.25e–01
VCC 1199	Early	8.6		5	−1.59	2.99	−5.12 ± 0.22	−5.12 ± 0.22	False	⋯
VCC 1440	Early	8.83		5	−1.61	2.85	−5.37 ± 0.15	−5.37 ± 0.15	False	⋯
VCC 1545	Early	8.85		5	−1.23	2.37	−6.27 ± 0.10	−6.28 ± 0.10	False	⋯
VCC 1627	Early	8.8		5	−1.27	3.08	−5.45 ± 0.12	−5.45 ± 0.12	False	⋯
[KK2000] 03	Early	8.16		5	−1.84	1.21	−7.2 ± 0.09	−7.2 ± 0.09	True	3.39e–01
[KK2000] 53	Early	6.85		5	−1.65	0.84	−7.62 ± 0.33	−7.62 ± 0.33	False	⋯
[KK98] 096	Early	7.07		5	−2.31	1.1	−6.37 ± 0.22	−6.37 ± 0.22	True	1.70e+00

Note. (1) Galaxy name; (2) galaxy type; (3) galaxy stellar mass [units of M_⊙]; (4) black hole mass [units of M_⊙]; (5) source of MBH mass: 1 = R. C. E. van den Bosch (2016), 2 = A. E. Reines & M. Volonteri (2015), 3 = J. E. Greene et al. (2020), 4 = J. E. Greene et al. (2020) M_BH–σ relations, 5 = J. E. Greene et al. (2020) M_BH–M_gal relations; (6) and (7) give the power-law fit parameters (γ = slope and ρ_{5 pc} = 3D stellar density at r = 5 pc [units of M_⊙ pc⁻³]) to the nuclear density profiles from Paper I; (8) full TDE rate with event horizon suppression; (9) full TDE rate without event horizon suppression; (11) flag to indicate if a Bachall–Wolf cusp was applied to the density profile at small radii; and (12) radius at which the Bachall–Wolf cusp begins [units of pc].

Download table as:  ASCIITypeset images: Typeset image Typeset image

Zoom In Zoom Out Reset image size

4.3. Pinhole Fraction

The pinhole fraction f_pinhole is a common metric presented for TDE rate estimates that describes the fraction of TDEs coming from the “full loss-cone” or “pinhole” regime (i.e., q > 1). These TDEs are interesting observationally because the pinhole scenario can produce TDEs with large penetration factors (β = r_t/r_p), while in the “empty loss-cone” regime, nearly all TDEs have β ≈ 1 (N. C. Stone & B. D. Metzger 2016; L. Broggi et al. 2022; although see also A. Weissbein & R. Sari 2017). This fraction is measured by integrating the loss-cone flux (Equation (10)) over orbital energies corresponding to q ≤ 1 and dividing this by the full integral over all energies.

Figure 5 shows f_pinhole for our sample galaxies as a function of MBH mass. Similarly to previous works, we find that f_pinhole increases with decreasing MBH mass. However, we do not observe a nearly 100% pinhole fraction for M_• ≲ 10⁶ M_⊙ as in the parametric fits provided by J. N. Y. Chang et al. (2025) and N. C. Stone & B. D. Metzger (2016), which are shown as the solid green and black dashed lines, respectively. Instead, we observe a convergence toward f_pinhole ≈ 0.75. Our sample contains mostly nucleated galaxies, with all our late-type and low-mass galaxies being nucleated (see Paper I for more details). This suggests that the high densities of the NSCs in these galaxies may place an upper limit on the pinhole fraction that has not been observed before.

Zoom In Zoom Out Reset image size

4.4. Angular Momentum Diffusion Distributions

While the pinhole fraction provides some insight into the dynamics of individual TDEs in a galaxy, we take this investigation a step further by deriving the normalized cumulative distributions of q values for each galaxy (see Figure 6). This is accomplished by multiplying the loss-cone flux () by ∣d/dq∣, which we compute with numerical derivatives. These results further highlight the preference for diffusive TDEs (the “empty loss-cone” regime) in more massive galaxies. Although most TDEs in a galaxy are generally sourced from stars with q ∼ 1, these distributions highlight subdominant contributions to the rates and how they differ with host galaxy mass. As galaxy mass decreases, the distributions tend to encompass larger q values, visibly highlighting the pinhole TDE preference for lower-mass galaxies.

Zoom In Zoom Out Reset image size

4.5. Massive Black Hole Spin Evolution

Assuming an isotropic distribution of impact trajectories, TDEs should serve to spin down an MBH over time (D. Merritt 2013b; B. D. Metzger & N. C. Stone 2016), in analogy to the “chaotic accretion” paradigm for MBH spin evolution through many active galactic nucleus (AGN) episodes (A. R. King & J. E. Pringle 2006). Based on this idea, we estimate the maximum dimensionless spin () allowed for each of our MBHs. While in principle, one should calculate by combining MBH growth histories (from e.g., mergers and AGN episodes) with a time-dependent TDE rate, here we present a simplified calculation that nonetheless captures basic trends. We imagine beginning with a maximally spinning (a_• = 0.997) black hole that experiences some mass growth over time due to a constant TDE rate. Under these assumptions, the maximum spin is given by (B. D. Metzger & N. C. Stone 2016)

where t_H is the age of the Universe (≈14 × 10⁹ yr) and we have assumed that half the star’s mass accretes during the TDE. Note that it is possible in extreme cases for , in which case we set it to zero.

Figure 7 shows the results of this calculation. It is clear that spin-down from TDEs is negligible for the highest-mass MBHs in our sample. However, TDEs appear to place nontrivial upper limits on the MBH spin for many lower-mass black holes. In particular, for M_• ≲ 10⁷ M_⊙, nuclear TDE rates are generally high enough to prevent time-averaged MBH spin values from reaching the Thorne limit of a = 0.997 (K. S. Thorne 1974). As MBH masses drop, so does (on average) , sometimes reaching values as low as ≈0.9, well below many astrophysical spin measurements from iron Kα spectroscopy (shown as gray crosses; C. S. Reynolds 2021).

Zoom In Zoom Out Reset image size

We note here that represents a temporally averaged limit on MBH spin. If indeed one finds (via X-ray observations or other techniques) a relatively small MBH (≲10⁵ M_⊙) with very large a, that could indicate (i) an abnormally low TDE rate in the galaxy, or, more likely (ii) that the galaxy has grown dramatically from an ongoing episode of coherent accretion, more than doubling its mass to spin up beyond the temporally averaged upper limit provided by TDEs. The existence of a preferred direction of TDE angular momentum (e.g., an NSC with net rotation) could also weaken this constraint, but except in the extreme case of a disk-like cluster, should only affect by a factor of order unity.

In this work, we have developed a new software package, REPTiDE, for computing two-body relaxation-based TDE rates given a galaxy density profile and MBH mass. The output from this code agrees well with previous implementations of the loss-cone procedure and will be a useful tool for future TDE rate studies enabled by the next generation of surveys. We demonstrate here that REPTiDE agrees well with the N. C. Stone & B. D. Metzger (2016) implementation of the standard loss-cone theory (Figure 1). While developed primarily for application to large samples of model galaxies, in this work we applied REPTiDE to the real, observed galaxies used to construct the new density scaling relations presented in Paper I.

Given our strict nuclear density requirements during sample construction (minimum radial resolution of 5 pc), the rate calculations presented in this paper represent the most accurate individual TDE rate estimations for a sample of nucleated galaxies. We found that our TDE rate distributions peak around log(M_gal/M_⊙) ≈ 9.4 and log(M_•/M_⊙) ≈ 6.5. We fit broken power laws directly to these distributions resulting in relations that describe the typical TDE rates as functions of galaxy and MBH mass.

The large variance in TDE rates between galaxies of similar properties (e.g., mass or MBH mass) means that the observed TDE rates will often be dominated by the galaxy subpopulations with the highest TDE rates. Additionally, TDE rate suppression due to direct capture heavily alters the observability of TDEs at the high-mass end. We, therefore, fitted a double broken power law to the binned mean of each distribution with rates that include event horizon suppression. We compared these with the most recent observed TDE rates from Y. Yao et al. (2023) and found excellent agreement. At galaxy masses ≲ 10^10.5 M_⊙, this sample is composed entirely of NSC-hosting galaxies (for both early and late types). Thus, achieving this agreement with a limited sample of nearby nucleated galaxies further highlights the importance of considering NSC density contributions when modeling TDE rates (e.g., H. Pfister et al. 2020).

We also presented the pinhole fraction for each of our sample galaxies, which describes the fraction of TDEs resulting from the “full” loss-cone (pinhole) regime (q > 1). Pinhole regime TDEs generally result in full disruptions while diffusive regime events (i.e., q < 1) can produce overwhelmingly partial TDEs (where the star is not fully destroyed on a single passage; L. Broggi et al. 2024). The distributions of individual TDE parameters in a galaxy could impact the observed TDE rate. In agreement with previous works (i.e., N. C. Stone & B. D. Metzger 2016; J. N. Y. Chang et al. 2025), we find that the pinhole fraction increases rapidly with decreasing MBH mass (Figure 5), possibly biasing larger SMBHs against producing (more easily observable) full disruptions. However, unlike in past works, we observed a convergence in the pinhole fraction to ≈0.75 at low MBH masses, which we attributed to the elevated densities in NSCs causing a nonnegligble fraction of TDEs to be sourced from very near the MBH. This has not been observed before, and should be tested against a sample of nonnucleated low-mass galaxies with similar density resolution.

The normalized cumulative q distributions for each galaxy (Figure 6) shows a more detailed view of the preference for pinhole TDEs in low-mass galaxies. They also highlight how the overall ranges of q values contributing to the TDE rates change with galaxy mass, where higher-mass galaxies produce TDEs down to much lower q values.

Lastly, we investigated the ability for TDEs to place upper limits on MBH spin (a_•), assuming that the quasi-isotropic nature of stellar orbital inclinations will serve to spin down an MBH. Assuming a near-maximal MBH spin (a_• = 0.997; K. S. Thorne 1974), we compute the maximum allowed spin () for each MBH (Equation (19)) assuming a constant TDE rate over a Hubble time (Figure 7). While TDE rates are expected to evolve with time, comparing these static calculations with other MBH spin measurements could place constraints on the time evolution of the TDE rate in a galaxy. AGN spin measurements that find a_• values far above can indicate ongoing episodes of coherent accretion.

In general, we do not find evidence for the “rate discrepancy” identified by some past semiempirical loss-cone modeling (J. Wang & D. Merritt 2004; N. C. Stone & B. D. Metzger 2016), i.e., a mismatch between (higher) theoretical rates and (lower) observed rates. Although it is possible that physical changes to simpler loss-cone models could reduce dynamical rates (K. Lezhnin & E. Vasiliev 2015; L. Broggi et al. 2024; O. Teboul & H. Perets 2025; O. Teboul et al. 2024), our underlying dynamical treatment is broadly the same as in past semiempirical work (N. C. Stone & B. D. Metzger 2016); the only thing that has changed is the underlying sample of galaxies we average over. We believe that this highlights the importance of constructing theoretical (per-galaxy or volumetric) TDE rate estimates from the most representative galaxy samples possible, something that is challenging to do using archival HST data alone because of the ad hoc selection criteria that have generally determined the availability of high-resolution nuclear HST photometry. One solution to this problem is to rely on a purely theoretical catalog of galaxies and nuclear dynamics (M. Polkas et al. 2024), but it is also possible to retain the direct link to observed nuclear dynamics.

In Paper III (C. H. Hannah et al. 2025, in preparation), REPTiDE will be applied to model galaxy samples constructed with existing galaxy scaling relations as well as the new NSC density scaling relations from Paper I. Each sample will be created with varying MBH occupation fractions and mass distributions based on expectations from different MBH formation scenarios. The intrinsic TDE rates computed for each model sample will be forward modeled into detection rates for ZTF, Rubin, and ULTRASAT. As ZTF is the current leader in TDE detections, immediate comparisons with the observed rates will be possible. However, ZTF samples are still small (≲100), so constraints may be limited in this case. Thus, we will forward model our predictions for Rubin and ULTRASAT as well, enabling future rate comparisons for expected near-future samples of thousands of observed TDEs.

C.H.H. and A.C.S. acknowledge support from National Science Foundation astronomy and astrophysics grant AST-2108180. N.C.S. acknowledges financial support from the Israel Science Foundation (Individual Research grant 2565/19) and the Binational Science Foundation (grant Nos. 2019772 and 2020397).

The number of orbital energies (n_ens) considered when calculating the TDE rate is essentially a free parameter that is directly linked to the runtime. Therefore, we performed a convergence test for TDE rates as a function of (n_ens) to minimize execution time while maintaining accuracy. Figure A1 shows the results of this investigation for NGC 4551.

Zoom In Zoom Out Reset image size