Loss cone shielding

ABSTRACT

A star wandering close enough to a massive black hole can be ripped apart by the tidal forces of the black hole. The advent of wide-field surveys at many wavelengths has quickly increased the number of tidal disruption events (TDEs) observed, and has revealed that (i) observed TDE rates are lower than theoretical predictions and (ii) E+A galaxies are significantly overrepresented. This overrepresentation further worsens the tension between observed and theoretically predicted TDEs for non-E+A galaxies. Classical loss cone theory focuses on the cumulative effect of many weak scatterings. However, a strong scattering can remove a star from the distribution before it can get tidally disrupted. Most stars undergoing TDEs come from within the radius of influence, the densest environments of the Universe. In such environments, close encounters rare elsewhere become non-negligible. We revise the standard loss cone theory to take into account classical two-body interactions as well as strong scattering, collisions, tidal captures, and study under which conditions close encounters can shield the loss cone. We (i) analytically derive the impact of strong scattering and other close encounters, (ii) compute time-dependent loss cone dynamics including both weak and strong encounters, and (iii) derive analytical solutions to the Fokker–Planck equation with strong scattering. We find that (i) TDE rates can be reduced to up to an order of magnitude and (ii) strong shielding preferentially reduces deeply plunging stars. We also show that stellar overdensities, one possible explanation for the E + A preference, can fail to increase TDE rates when taking into account strong scattering.

1 INTRODUCTION

A star wandering too close to a massive black hole (MBH) can be ripped apart by the MBH’s tidal forces. Roughly half of the gaseous debris from the disrupted star falls back on to the MBH, eventually circularizing into an accretion disc and powering a luminous flare (Rees 1988; Evans & Kochanek 1989).

Tidal disruption events (TDEs) have been theoretically predicted since the 1970s (e.g. Hills 1975) but were not observed until the advent of all-sky surveys (beginning with the ROSAT X-ray all-sky survey; Bade, Komossa & Dahlem 1996; Komossa & Bade 1999). The sample of observed TDEs has rapidly grown in the last decades and now spans all wavelengths, from the radio to gamma-rays. The last decade has seen a high rate of TDE discovery in wide-field optical surveys such as Sloan Digital Sky Survey (van Velzen et al. 2011), Pan-STARRS (Gezari et al. 2012; Chornock et al. 2014; Nicholl et al. 2019), PTF (Arcavi et al. 2014), iPTF (Hung et al. 2017; Blagorodnova et al. 2019), ASASSN (Holoien et al. 2016; Hinkle et al. 2021), and ZTF (van Velzen et al. 2019; Hammerstein et al. 2023). X-ray surveys are also returning to prominence as a driver of TDE discovery; at present, the SRG/eROSITA all-sky survey (Sazonov et al. 2021) seems to be finding ∼20 TDEs per year, comparable to the optical TDE discovery rate of ZTF. As time domain capabilities continue to expand, the observed sample of TDEs is expected to grow by orders of magnitude, with hundreds to thousands of new TDEs expected to be found in optical wavelengths by the LSST/Vera Rubin Observatory (Ivezić et al. 2019; Bricman & Gomboc 2020), and in the near-ultraviolet by ULTRASAT (Ben-Ami et al. 2022).

The opening of this new era for TDEs has allowed observers to empirically constrain the per-galaxy TDE rate to a range |$\dot{N}_{\rm TDE} \sim 10^{-5}-10^{-4}$| yr⁻¹ gal⁻¹ (Holoien et al. 2016; van Velzen 2018), and revealed a surprising over-representation of these events in post-starburst (e.g. ‘E + A’) galaxies (Arcavi et al. 2014; French, Arcavi & Zabludoff 2016, 2017; Law-Smith et al. 2017; Graur et al. 2018; French et al. 2020; Hammerstein et al. 2021). The latest estimations suggest that E+A galaxies are overrepresented by a factor of ≈22 (Hammerstein et al. 2021). This further reduces the observationally inferred TDE rate in non-E+A galaxies, which is already in tension with theoretical (but empirically calibrated) rate estimates of |$\dot{N} \sim 10^{-4}-10^{-3}$| yr⁻¹ gal⁻¹ (e.g. Wang & Merritt 2004; Stone & Metzger 2016).

TDE rates are traditionally estimated from classical loss cone theory, which focuses on the cumulative effect of many weak scatterings. These scatterings are modelled as local and uncorrelated, and lead to effective diffusion coefficients describing the drift and diffusion of stellar populations through phase space (see Merritt 2013; Stone et al. 2020, for reviews). Strong or small impact parameter scatterings are generally neglected, as they are largely outnumbered by weak scatterings.

However, strong scatterings can have a qualitatively different effect than weak ones on stellar populations. In particular, a strong scattering can eject a star from the distribution before it gets tidally disrupted. Most stars undergoing TDEs come from within the radius of influence, the densest environments in the Universe. Such an environment can facilitate ejection in strong scatterings, but it also increases the rate of close encounters that are rare elsewhere. Such close encounters can be divided into star–star interactions, including collisions and tidal captures, and star–stellar mass black hole (sBH) interactions, such as μ-TDEs and tidal captures. In this paper, we will revise the standard loss cone theory to take into account not simply weak two-body scatterings, but also destructive processes such as strong scattering ejections, collisions, and tidal captures, with the ultimate goal of determining under which conditions close encounters can shield the loss cone.

Close encounters are highly dependent on both the stellar and sBH slopes. An important prediction of stellar dynamics is the redistribution of arbitrary initial distributions of stars around MBHs over time into quasi-stationary states (QSSs). In the presence of one population of stars, a Bahcall–Wolf density cusp with stellar density ρ ∝ r^−γ, will develop, with γ = 7/4 (Bahcall & Wolf 1976). If the system is composed of both stars and heavier objects, the heavier objects are expected to segregate towards the centre of the galactic nucleus and settle on a steeper cusp γ = 7/4 − 2, while the light objects will have a weaker cusp γ ≈ 1.3−1.5 (Bahcall & Wolf 1977; Preto & Amaro-Seoane 2010; Amaro-Seoane & Preto 2011; Broggi et al. 2022). In the strong segregation limit, when the heavier objects are relatively rare, the heavy objects have been predicted to settle to even steeper power-law slopes of γ = 2 − 11/4 (Alexander & Hopman 2009). Moreover for an adiabatically growing MBH, the stellar slope may assume steeper values of γ ≳ 2 (Young 1980).

We introduce the different kinds of close encounters that are relevant for our problem in Section 2, and derive their analytical rates in Section 2.2. In Section 3, we solve the Fokker–Planck equation in angular momentum space in the presence of close encounters. In Section 3.3, we present our analytical time-dependent solutions of the Fokker–Planck equation with strong scattering. In Section 4, we compute the impact of close encounters on TDE rates and discuss our findings for the E+A preference. In Section 5, we discuss the implications of our work and briefly summarize it.

2 CLOSE ENCOUNTERS

In a spherical galactic nucleus, one can characterize the population of stars with an orbit-averaged distribution function f(ϵ, J), where ϵ is the orbital energy and J the orbital angular momentum of stars. A variety of processes can remove stars from the distribution f. Most notably, a star of mass M_⋆ and radius R_⋆ will be destroyed by the central MBH (of mass M_•) if it completes an orbit with galactocentric pericenter R_p less than the tidal radius,

A finite R_t (i.e. the presence of a central MBH) carves a conical gap in velocity space (at fixed radius), or equivalently places a non-zero lower limit on the angular momentum that particle orbits may have (Frank & Rees 1976; Lightman & Shapiro 1977; Cohn & Kulsrud 1978). Particles with less angular momentum suffer destructive interactions at pericenter.

However, before the typical star can get to an orbit with R_p ∼ R_t, its pericenter must endure a long and potentially dangerous random walk through regions of high stellar and sBH density. During this random walk, sufficiently close encounters with the stellar/sBH population can eliminate our test star through a variety of different processes. Here, we overview the relevant ones, roughly in order of increasing cross-section:

• Direct physical collisions: Star–star (or star–sBH) collisions with a pairwise pericenter p < R_⋆ will often destroy the test star. At large distances r from the MBH, where the velocity dispersion |$\sigma \ll V_\star =\sqrt{2GM_\star /R_\star }$|⁠, the stellar surface escape speed, collisions with other stars can be coagulative rather than destructive.¹

• μ-TDEs: Slightly more distant encounters with a sBH [p < r_t, where r_t is the same as R_t from equation (1) but with the sBH mass m_bh substituted into equation (1) for M] can destroy the test star through the action of the sBH’s tides in a scaled down version of a standard TDE. The observational signatures of these ‘μ-TDEs’ are poorly understood (Perets et al. 2016). However, the different dynamics and greatly reduced Eddington limit of an sBH (compared to an MBH) means that the resulting transient is unlikely to be mistaken for a standard TDE, and a star experiencing this fate can be safely removed from the distribution function.

• Tidal captures: At slightly larger pericenter distances than the μ-TDE radius, tides can excite mechanical oscillation modes in the test star. If this process converts enough orbital energy into mechanical energy, then an unbound parabolic/hyperbolic two-body orbit will capture into a bound binary (Fabian, Pringle & Rees 1975; Press & Teukolsky 1977; Lee & Ostriker 1986). A tidal capture can only occur when the relative kinetic energy between the two objects is less than the binding energy of the test star. The long-term evolution of a tidal capture system is not entirely clear (Generozov et al. 2018). Some recent hydrodynamical simulations indicate that the captured star will be destroyed in a brief sequence of runaway partial (micro)-TDEs (Kremer et al. 2022), but even if this is not the case, a star that has become tightly bound to a sBH is unlikely to ever experience a TDE. If a tight binary’s centre-of-mass angular momentum starts falling to levels approaching the loss cone, then direct merger (i.e. μ-TDE) becomes the most probable outcome (Bradnick, Mandel & Levin 2017), followed by tidal separation via the Hills mechanism. Star-star tidal captures are also possible, but again the most likely final result (if the binary is excited towards small p values) is a merger, which would effectively remove one star from the distribution function.

• Ejections: A strong (large-angle) scattering between the test star and another stellar-mass object can accelerate it above the local escape speed, v_esc. When this happens, the star is removed from the system (Lin & Tremaine 1980) and can, under certain conditions, become a hypervelocity star (Yu & Tremaine 2003; O’Leary & Loeb 2008).

We will now calculate the different rates of the processes that can remove stars from the distribution. We will begin with an order of magnitude estimate in Section 2.1 and then present our analytical derivations for strong scattering in Section 2.2, star–star interactions in Section 2.3, and star–stellar mass black hole (BH) interactions in Section 2.4.

2.1 Approximate description

The rates of processes (i–iv) are governed by basically the same underlying physics, and can be calculated in the usual ‘n − Σ − v’ way as

where |$\dot{N}$| is the per-target rate of close encounters, p_c is the critical pericenter for the close encounter (e.g. direct collision, μ-TDE, or tidal capture) to occur, v_∞ is the relative velocity at infinity, M_⋆ is the mass of the target star, M₂ is the mass of the secondary (another star for star–star encounters or a sBH for TCs/μ-TDEs), and n is the number density of secondaries.

Far from the central MBH, the rates of these close encounters will be set by gravitational focusing, but close to the MBH, encounter rates will be determined by straight-line flyby trajectories. We will now approximate the rates of close encounters in these two regimes by making the simplified assumption that all stars have the same mass M_⋆ and are distributed in space with a power-law |$n_\star (r) \propto r^{-\gamma _\star }$|⁠, and likewise that all sBHs have a common mass m_bh and a power-law density profile |$n_{\rm bh}(r) \propto r^{-\gamma _{\rm bh}}$|⁠. We further approximate the situation by assuming that all relative velocities are equal to the stellar velocity dispersion,

sparing us the need to integrate over a relative velocity distribution. Under these simplifying assumptions, the (per-target) rates of star–star collisions, star–sBH μ-TDEs, and star–sBH tidal captures scale, respectively, as

Here, the transition radii between the focused and unfocused interaction regimes are

and we have further defined α = M_⋆/m_bh as the sBH–star mass ratio, and λ ≈ 2 as the rough enhancement to the μ-TDE cross-section for tidal captures.² In the approximate equalities above, we have taken γ_⋆ = γ_bh = 7/4, λ = 2, and α = 0.1 for concreteness, and used the shorthand variables M₆ = M_•/(10⁶ M_⊙), m_⋆ = M_⋆/M_⊙, and r_⋆ = R_⋆/R_⊙.

It is important to note that tidal capture becomes impossible when the relative kinetic energy of the target and secondary becomes larger than the maximum amount of mechanical energy that can be injected into the target star’s normal modes in a single pericenter passage. Crudely approximating this condition as |$M_\star m_{\rm bh} (M_\star + m_{\rm bh})^{-1} \sigma ^2/2 \gt u_\star GM_\star ^2/R_\star$|⁠, where u_⋆ < 1 is a fudge factor representing the fraction of the star’s binding energy that can be absorbed by mechanical oscillations, we find that tidal capture deactivates at radii below

Because |$r_{\rm TC} \sim r_{\rm TC}^{\rm tr}\alpha ^{-2/3}u_\star ^{-1}$|⁠, and most likely |$u_\star \sim \mathcal {O}(0.1)$|⁠, r_TC is generally greater than |$r_{\rm TC}^{\rm tr}$|⁠, and tidal captures are always gravitationally focused.

We also want to understand the basic scaling laws describing rates of stellar ejection through strong scatterings. We approximate the specific impulse delivered to the target star in a strong scattering as Δv ∼ (Gm₂/p²)(p/v_p), with v_p the velocity at pericenter. Angular momentum conservation lets us rewrite this in terms of the velocity (v_∞) and impact parameter (b) at infinity: Δv ∼ Gm₂/(bv_∞). Thus, the critical impact parameter that allows strong scatterings to eject stars will be b_ej = Gm₂/(vΔv_ej), where Δv_ej is the critical specific impulse to eject a star. While it is tempting to equate |$\Delta v_{\rm ej} = v_{\rm esc}(r) = \sqrt{2GM/r}$|⁠, we must remember that for stars on loss cone orbits, 1 − e ≪ 1 and most of the orbit is spent traveling at speeds v(r) ≈ v_esc(r). Therefore, we approximate δv_ej ≈ v_esc(r) − v(r), or

for a star on a highly eccentric orbit with r ≪ a. We now see that the critical impact parameter for ejection is |$b_{\rm ej} \sim 2\sqrt{2} a m_2/M$|⁠. The local ejection rate from scatterers with number density n will thus be |$\dot{N}_{\rm ej} = n \pi b_{\rm ej}^2 v_\infty$|⁠, so that

Here, we have left m₂ and γ general, as they may represent either star–star or star–sBH strong scatterings. In contrast to equations (4)–(6), there is no leading-order distinction between gravitationally focused or unfocused regimes for strong scatterings. Furthermore, there is no leading-order dependence on the mass m_⋆ of the test star experiencing the scattering (although note that this equation implicitly assumes m_⋆ ≲ m₂).

While equations (4)–(6) and (12) give local encounter rates |$\dot{N}(r)$|⁠, we also want to compute their orbit-averaged versions, |$\langle \dot{N} \rangle$|⁠. Because of our ultimate interest in the loss cone, we consider highly eccentric orbits with pericenter p ≪ a, the stellar semimajor axis. At the order of magnitude level, we can approximate the orbit-averaged rate of some process at a given radius as |$\langle \dot{N}(r) \rangle \sim \dot{N}(r) \times (r/a)^{3/2}$|⁠, and determine whether the rates are pericenter- or apocenter-dominated to get |$\langle \dot{N} \rangle$|⁠.

For ejection in strong scatterings,

so ejection rates will be pericenter-dominated as long as γ > 1, which we will assume for the remainder of this paper, as the shallowest profiles are γ ∼ 1.3–1.5 (e.g. Bahcall & Wolf 1977). For other types of close encounters, encounter rates will likewise be pericenter-dominated at small radii r < r^tr so long as γ > 1, but the situation at large radii r > r^tr is more complicated. From equations (4)–(6), we see that close encounter rates will always be pericenter-dominated when γ > 2, but when 1 < γ < 2, the situation depends on whether p > p^tr, a transitional pericenter.³ More specifically, we find that when 1 < γ < 2, close encounter rates scale as

Here, we have approximated the transitional pericenter as

Note that it is possible that p^tr > a; in this situation, rates are always pericenter-dominated, so even though p^tr no longer has a physical meaning, the inequality conditions in equation (14) are still accurate.

Alternatively, when γ > 2,

Because of the many permutations of possibilities, we have not disaggregated collision, μ-TDE, and tidal capture rates in the above formulae, instead referring to a ‘close encounter radius’ R_ce which will be equal to R_⋆, r_t, or r_TC for these three processes, respectively. Likewise, we have denoted a general transition radius r^tr differentiating between focused and unfocused local encounter rates; the appropriate transitional radius should be chosen from among equations (7)–(9).

With these scaling relations in hand to provide physical intuition, we now proceed to more rigorously derive ejection and close encounter rates.

2.2 Strong scatterings

In this section, we present our analytical derivations of orbit-averaged rates for strong scattering. We employ a stellar distribution function calculated under the assumptions that (i) the stellar density profile |$\rho (r) = \rho _{\rm infl} (r/r_{\rm infl})^{-\gamma _\star }$|⁠, and (ii) the potential ψ = GM_•/r is everywhere Keplerian.⁴ Under the assumption of isotropy, we apply an Eddington integral to ρ to obtain

where the (positive-definite) specific orbital energy is, for a given star at radius r and velocity v, ϵ = ψ(r) − v²/2; 〈m_⋆〉 is the average mass in the stellar population; and r_infl is the radius that encloses a mass in stars equal to the MBH mass.

2.2.1 Equal mass scatterer

We first consider strong scattering from an equal mass scatterer, e.g. from other stars. We will assume that all stars have the same mass m_⋆ and that the distribution of velocities is isotropic at every point. Let us consider a test star whose velocity is |${\boldsymbol V}$|⁠. The probability for such a star to have an encounter during a time dt that will increase its velocity to |${{\boldsymbol V} + \delta {\boldsymbol v}}$| is (Hénon 1960a):

where f_⋆(v) is the distribution function of stars and v₀ is given by

Let θ be the angle between |$\boldsymbol {V}$| and |${\delta \boldsymbol v}$|⁠, with v_esc the escape velocity at this point. The star will be ejected if

Replacing Cartesian coordinates by spherical coordinates, the local ejection rate, i.e. the probability for one star to be ejected during a time dt, is given by (Hénon 1960b)

Assuming Keplerian motion and an escape velocity |$v_{\rm esc}=\sqrt{2GM_\bullet /r}$|⁠, the local ejection rate becomes

with V the local Keplerian velocity of the test star and a its semimajor axis. As the relaxation time is much longer than an orbital period, the local ejection rate per star can then be orbit-averaged following equation (A2) as

with ν the true anomaly, and e the eccentricity of the star.

The orbit-averaged ejection rate does not have a general closed form. However, for relevant values of γ_⋆ which are integers or half-integers, an analytical solution exists. Analytic solutions for some physically motivated values of γ_⋆ are presented in Appendix  A.

2.2.2 Unequal mass scatterer

Let us now consider strong scattering from an unequal mass scatterer, i.e. stellar mass BHs with α = m_⋆/m_bh. The probability for a star with velocity |${\boldsymbol V}$| to have an encounter during a time dt that will increase its velocity to |${{\boldsymbol V} + \delta {\boldsymbol v}}$| becomes

where f_bh(v) is the distribution function of stellar mass BHs (also assumed to reflect a power-law density profile, with |$\rho _{\rm bh} \propto r^{-\gamma _{\rm bh}}$|⁠), and v₀ is given by

The condition for the star to escape remains unchanged. Therefore, after some algebraic manipulations, the local ejection rate becomes

where

and the integration limits are given by

The local ejection rate per star can then be orbit averaged as in equation (A2). This orbit-averaged ejection rate does not have a general closed form, although the local ejection rate (equation 26) can be written in closed form for physically motivated values of γ_bh.

2.3 Star–star encounters: physical collisions and tidal captures

Tidal captures and physical collisions between stars, although rare in most astrophysical environments, are non-negligible in the extreme densities of galactic nuclei. Assuming that all stars have the same mass m_⋆, and that the test star has a velocity |$\boldsymbol {V}$|⁠, the local collision rate is given by

with the gravitationally focused impact parameter |$b^2 = r_{\rm coll}^2 + 4 G m r_{\rm coll}/V_0^2$| and |$V_0=| \boldsymbol {V}- \boldsymbol {v}|$|⁠. To take into account both collisions and tidal captures from stars, we can rewrite the formula as

with |$b_{\rm max}^2 = r_{\rm max}^2 + 4 G m r_{\rm max}/V_0^2$|⁠, where r_max = max[r_coll, r_capt(V₀, m_⋆)] and r_capt = 2R_{t, ⋆} ≈ 2R_⋆ is the approximate tidal capture radius we employ for a high-σ environment (Stone, Küpper & Ostriker 2017). This formula is only an approximation to a more rigorous calculation of the tidal capture radius (see e.g. Lee & Ostriker 1986), but it will generally be correct to within |$\sim 10~{{\ \rm per\ cent}}$| in a high-σ galactic nucleus (Generozov et al. 2018).

Note that tidal capture cannot occur if the relative velocity |$V_0 \gtrsim \eta v_\star = \eta \sqrt{Gm_\star /r_\star }$|⁠, where η < 1 is a dimensionless number of uncertain magnitude. When V₀ ≪ v_⋆, tidal capture can proceed cleanly through the linear excitation of normal modes in the star. However, when V₀ ∼ v_⋆, the oscillation modes reach non-linear amplitudes and may begin to dissipate their energy rapidly through non-linear mode–mode couplings or by steepening into shocks. Rapid dissipation of the mode mechanical energy can inflate the star and result in its destruction via a sequence of partial μ-TDEs (Kremer et al. 2022). The fudge factor η encodes the boundary between successful and unsuccessful captures; for simplicity, in the remainder of this paper we will take η = 0.1. Regardless of the true value of η, it is clear that at small galactocentric radii, stars are tidally disrupted rather than captured.

The orbit-averaged star–star collision/tidal capture rate |$\langle \dot{N}_{\star , \star } \rangle$| are computed following standard orbit-averaging procedures (equation A2).

2.4 Star–sBH encounters: μ-TDEs and tidal captures

A star wandering too close to a stellar mass BH can be tidally captured or tidally disrupted. The characteristic tidal radius is r_t = r_⋆(m_bh/m_⋆)^1/3, while the maximum initial pericentre distance resulting in tidal capture, r_capt, is approximated as before by r_capt = 2r_t, with the same caveat that large V₀ are treated as μ-TDEs rather than tidal captures.

Considering a test star with velocity |$\boldsymbol {V}$|⁠, its probability to be tidally captured or disrupted by one of the stellar mass BHs is given by

where

For practical evaluation, the tidal capture rate is orbit averaged following equation (A2).

2.5 Comparison between the different close encounters

Given the many close encounters at play in a galactic nucleus, it is useful to compare their relative rates. We perform this comparison in Fig. 1, which illustrates how strong scattering rates and other close encounter rates vary as one changes stellar eccentricity e, and various properties of the central cusp (e.g. γ_⋆, γ_bh). The rate of angular momentum relaxation,

is also plotted as a useful reference rate [here μ(ϵ) is the diffusion coefficient as computed in equation (35)]: loss cone shielding is likely to only become important in so far as the rate of a destructive close encounter can become comparable to |$\dot{N}_{\rm AM}$|⁠.

We see that in general, |$\dot{N}_{\rm AM} \gg \dot{N}_{\rm ej}$| so long as the relevant power-law indices γ ≲ 2. In contrast, γ = 5/2 leads to ejection rates at large a that can exceed |$\dot{N}_{\rm AM}$|⁠. At small a, ejection rates decline and can be exceeded by close encounter rates, but regardless of γ and a, |$\dot{N}_{\rm AM} \gg \dot{N}_{\rm ce}$| always.

To summarize, these simple rate comparisons show that (i) destructive close encounters will never effectively shield the loss cone; (ii) ejections in strong scatterings are likely to shield the loss cone for a ∼ r_infl and γ > 2; (iii) at a ≪ r_infl, strong scatterings become subdominant to close encounters as a way to remove stars from the distribution function. Ultimately, however, these rate comparisons only provide an order of magnitude check of the competition between different physical processes, so we now move on to computing time-dependent solutions to the loss cone problem in the presence of strong scattering and close encounters.

3 THE SHIELDED LOSS CONE

3.1 Loss cone theory

The radial distribution of stars and compact objects evolves over time due to two-body relaxation. In the continuum limit, and assuming spherical symmetry, the stellar population is represented with a distribution function f(ϵ, J), where J is the specific angular momentum of a stellar orbit. For the near-radial orbits relevant for TDEs, angular momentum relaxation is much faster than energy relaxation and we will use a two-time-scale argument: we separate f(ϵ, J) = f_ϵ(ϵ)f_j(J) and assume a frozen distribution of energy. By assumption, stars are fixed in bins of orbital energy, but allowed to diffuse through angular momentum space in a random walk. This process is captured by the orbit-averaged Fokker–Planck equation (Merritt & Wang 2005),

where |$j \equiv J/J_{c}(\epsilon)=\mathcal {R}^{1/2}$| is a dimensionless angular momentum variable (normalized by the angular momentum of a circular orbit, J_c) and τ ≡ μ(ϵ)t ≈ t/t_r is a dimensionless version of time t, with μ(ϵ) the orbit-averaged diffusion coefficient at specific energy ϵ:

Here, P(ϵ) is the orbital period of a radial orbit of energy ϵ, v_r is the star’s radial velocity, and the local diffusion coefficient |$\langle (\Delta \mathcal {R})^2\rangle$| is presented in Appendix  B. Hence, τ ≈ t/t_r, with t_r, the energy relaxation time.

In the presence of strong scatterings and destructive close encounters, stars can be ejected or otherwise removed from the distribution, which we model by adding a sink term to the Fokker–Planck equation:

Here, |$\langle \dot{N}_{\rm ej}\rangle$| is the orbit-averaged rate of close encounters, as derived in Section 2.2 for strong scattering and Section 2.3 for star–star encounters and Section 2.4 for star–BH encounters.

The intial and inner boundary conditions depend on a dimensionless diffusivity parameter |$q(\epsilon)=\mu (\epsilon)P(\epsilon)/j^2_{\rm lc}(\epsilon)$|⁠, which is defined in terms of the size of the loss cone in dimensionless angular momentum space (⁠|$j_{\rm lc} \approx \sqrt{ 2GM R_{\rm t}}$|⁠). The value of q determines whether the loss cone evolves in the ‘empty’ (q ≪ 1; stars immediately destroyed once j ≤ j_lc) or ‘full’ (q ≫ 1; stars may move in and out of the loss cone multiple times per orbit) limits.

When q(ϵ) ≪ 1, one can assume an absorbing boundary condition at the loss cone, and a zero-flux boundary condition at j = 1:

Conversely, when q(ϵ) ≫ 1, the distribution function does not go to zero until a much smaller value of dimensionless angular momentum, j₀ = j_lc(ϵ)exp (− α/2), where

is an approximate flux variable that smoothly bridges the empty and full loss cone limits, which has been derived by returning to the local (non-orbit-averaged) Fokker–Planck equation and determining how f varies with radial phase assuming f = 0 at periapsis (Cohn & Kulsrud 1978; Merritt 2013).⁵ We thus have a more general set of boundary conditions that we use for all numerical solutions:

In both cases, we take as an initial condition an isotropic distribution, f(j) = 1, for j₀ ≤ j ≤ 1, and set f(j) = 0 elsewhere.

The flux of stars that scatter into the loss cone per unit time and energy is given by

Then, the TDE rate is obtained by integrating |$\mathcal {F}(\epsilon)$| across many bins of energy ϵ, such as

3.2 Numerical results

We begin by solving the modified Fokker–Planck equation (equation 36) numerically. An example of these time-dependent solutions is presented in Fig. 2, where we show the evolution of the angular momentum distribution (in a bin of fixed ϵ), both with and without strong scatterings, in the diffusive regime (q = 10⁻²). In the absence of strong scatterings, we see that we quickly achieve the logarithmic Cohn–Kulsrud profile (Cohn & Kulsrud 1978). However, when strong scatterings are turned on, a range of outcomes are possible. For (γ_⋆ = 3/2 and γ_bh = 2), both stars and sBHs have a very small impact at early times and a moderate one at later times. The distribution function still has a very similar form to the logarithmic Cohn–Kulsrud solution, albeit with a small depletion almost evenly distributed across angular momentum. For γ_bh = 5/2 (or an ultra-steep density star profiles γ_⋆ = 5/2), the distribution function can be radically modified relative to the standard Cohn–Kulsrud solution, with a major depletion of low-j orbits. The ‘cavity’ carved out of the distribution function at low j in this case reflects the impact of strong scattering, which preferentially ejects stars on the most radial orbits (as it is generally the case that ejection rates for an individual star are pericenter-dominated, see e.g. Equation 23).

When strong scatterings are capable of strongly depleting the low-j end of the distribution function, the result is a major reduction in loss cone flux (as |$\mathcal {F} \propto (\partial f_j/\partial j)|_{\rm j=j_{\rm lc}}$|⁠). We further explore this effect by computing flux reduction for different BH slopes and BH mass, applying equation (40) to our numerical solutions, and plot the results in Fig. 3. In this figure, we see that in all cases, there is relatively little evolution of the flux suppression factor |$\mathcal {F}/\mathcal {F}_{\rm w}$| (here |$\mathcal {F}_{\rm w}$| is the flux into the loss cone that would be achieved without the presence of strong scatterings). Suppression factors |$\mathcal {F}/\mathcal {F}_{\rm w}$| are modest ≈1 for γ_bh = 2, reach the factor of 2 level for γ_bh = 9/4, and can produce multiple-order-of-magnitude drops in the TDE rate for steeper γ_bh values. As shown in equation (26), we see stronger levels of flux suppression for larger values of m_bh.

3.3 Analytical solutions of Fokker–Planck with strong scattering

The flux of stars inside the loss cone |$\mathcal {F} (\epsilon)$| is a sharply peaked function, and most of the stars come from energies close to the influence radius. As we have shown, in Section 2, the dominant close-encounter mechanism at influence radius is strong scattering. Therefore, we searched for an analytical solution to the Fokker–Planck equation with a sink term corresponding to strong scattering (equation 36).

3.3.1 Equal mass scatterer

For an equal mass scatterer, we found an exact analytical form for the orbit-averaged ejection rate for all relevant integer or half-integer values of γ_⋆. The formulas for γ_⋆ = 3/2, 5/2 can be found in Appendix  A.⁶ We assume a slow variation with respect to τ such as df/dτ ≈ 0, and use the method of Frobenius to derive an analytic solution in angular momentum to the modified Fokker–Planck equation for these two power-law distributions. For the shallower slope, γ_⋆ = 3/2, the solution is very close to the Cohn–Kulsrud profile, whereas for the steeper slope of γ_⋆ = 5/2, it has an exponential behaviour. The solutions for γ_⋆ = 3/2 and γ_⋆ = 5/2 can be written explicitly as

where |$A = \langle \dot{N}_{ej} \rangle / \mu (\epsilon)$| is a per-star sink term, |$\langle \dot{N}_{ej} \rangle$| can be found in Appendix  A, |$g_{\pm } = 1 \pm \frac{\sqrt{j}}{32 \sqrt{ 2 A} } + 0.0022 \frac{j}{ A}$| and two undetermined constants – b and a – need to be found. These solutions are valid at small j, which is the relevant parameter space for stars that could undergo a disruption. Their explicit time dependence comes from the undetermined constants a and b. In practice, b is deterministically set as a function of a using the absorbing boundary condition at j = j_o, so there is only one true time-dependent free parameter. We find that a depends on the dimensionless time τ as a ∼ τ^−1/2. Since the approximate time-scale for angular momentum relaxation to occur is t_j(j) ∼ j²t_r, our solution is accurate for |$j \lesssim j_{\rm CK} = \sqrt{\tau }$|⁠. In Fig. 4, the numerical solutions of the Fokker–Planck equation with strong scatterings (equation 36) are compared to our analytical solution (γ_⋆ = 5/2) at different dimensionless times. The agreement is excellent for values of j ≲ j_CK, the ‘Cohn–Kulsrud’ angular momentum below which the distribution function has had time to relax into a QSS.

3.3.2 Unequal mass scatterer

For an unequal mass scatterer, we found close form solutions for the local rates and derived approximate solutions for the orbit-averaged rates. We then repeat our application of the Frobenius method to find time-dependent solutions for physically motivated slopes. Specifically, we find solutions for the time evolution of the distribution function of stars scattering off a population of sBHs with power-law slopes of indices γ_bh = 3/2, 7/4, 2, 9/4, and 5/2. The resulting closed form solutions are as follows:

where γ_E is Euler’s constant and |$h_{\pm } = 1 \pm \frac{\sqrt{j}}{32 \sqrt{A} } + 0.0044 \frac{j}{ A}$|⁠. As before, explicit time dependence enters through two undetermined coefficients, a and b, although application of the absorbing inner boundary condition allows us to eliminate one of these. We find that the remaining coefficient can be determined quite accurately with a matching procedure: we set the value of our solution, |$f_{\gamma _\star , u}(x j_{\rm cK})$| equal to the value of the Cohn–Kulsrud distribution at j = xj_CK. Here, x is a dimensionless number of order unity; we find that x ∼ 0.1 yields a good match.

In Fig. 5, we compare the analytical solutions we obtained for unequal-mass strong scatterings at various γ_bh with the numerical solutions to equation (36). In all cases (analytic and numerical), the comparison is done at a fixed dimensionless time τ = 0.1. As can be seen, the agreement between analytical and numerical solutions is excellent for |$j \lesssim \sqrt{\tau }$|⁠.

4 IMPACT ON TDE RATES

4.1 Rates

As we have seen in Fig. 3, the flux inside the loss cone reaches a quasi-stationary state after a time t ∼ 0.1 t_r, with t_r, the relaxation time. Therefore, in this section we explore the impact on TDE rates of strong scattering after a time where the quasi-stationary state has been reached. To compute numerical fluxes we compute the diffusion coefficients μ(ϵ) for a grid of energies. Equation (36) is then integrated forward for the appropriate dimensionless time at each ϵ-value and the flux computed from equation (40). Finally, TDE rates are obtained by integrating across the energy bins (equation 41). To compute analytical fluxes, we use our solutions presented in Section 3.3 for the appropriate dimensionless time at each ϵ-value, then compute the fluxes with equation (40) and finally integrate across multiple energy bins to obtain TDE rates (equation 41).

In Fig. 6, we show energy-integrated TDE rates as a function of MBH mass both with and without the effects of strong scattering for cusp galaxies with radius of influence |$r_{\rm infl} = 11~{\rm pc} \left(\frac{M_\bullet }{10^8 {\rm M}_\odot } \right)^{0.58}$| (Stone & Metzger 2016).

We can note that our TDE rates without strong scattering shown in Fig. 6 are a factor of a few higher than in previous loss cone calculations (e.g. Wang & Merritt 2004; Stone & Metzger 2016). This effect is due to the inclusion of stellar mass BHs. Indeed, stellar mass BHs dominate the total relaxation rate and increase diffusion coefficients hence increasing the TDE rates by a factor of a few. This effect is further explored in Fig. 7, where we show the influence of different sBH masses on TDE rates both with and without strong scattering.

We see here that the effects of strong scattering are quite modest for γ_bh = 2 (reducing TDE rates by |$\lesssim 50~{{\ \rm per\ cent}}$|⁠), rise to a factor of 2 reduction for γ = 9/4, and then become very large (at least an order of magnitude reduction) for γ = 5/2. It is important to note that the TDE rates we obtain with γ_bh = 5/2 are in agreement with the observed TDE rates. In this example, we have taken m_bh = 30 M_⊙. The impact of different SBH masses is explored in Fig. 7.

In Fig. 7, we show that in the absence of strong scattering, TDE rates increase by a factor ≈3 as m_bh increases from 10 to 50 M_⊙. In the presence of strong scattering, however, the TDE rate is substantially lower at all values of sBH mass and is relatively independent of m_bh. Although larger m_bh values still increase diffusion coefficients, this is offset by the greater rate at which stars are ejected in close encounters (note that both effects scale |$\propto m_{\rm bh}^2$|⁠).

In these calculations, we considered a broad range of sBH masses, from m_bh = 10−50 M_⊙. Here, we are motivated partially by the ‘initial mass function’ for sBHs formed via the deaths of massive stars. Depending on the metallicity of the stellar population, the maximum masses for sBHs formed in this way is likely to be ≈50 M_⊙ (Spera & Mapelli 2017), although even with a favourable metallicity, the large majority of sBHs will have substantially lower masses.⁷ However, more recent simulations have shown that 30 M_⊙ can be formed at stellar metallicity (Bavera et al. 2022). In addition, once a population of sBHs is situated in a galactic centre environment, it may be able to gain mass through mergers and accretion. Most notably, any episode of large-scale MBH accretion will grind down and capture many of the pre-existing sBHs into the MBH accretion disc (Syer, Clarke & Rees 1991); once embedded in this disc, they can be processed to larger masses either via mergers (Bellovary et al. 2016; Tagawa, Haiman & Kocsis 2020) or from direct gas accretion (Gilbaum & Stone 2022). After the AGN episode ends, the population may relax into a more spherical configuration of the type we assume here (although see Szölgyén & Kocsis 2018). Efficient sBH growth through merger may also occur during periods of time when a dense nuclear star cluster lacks a MBH (Antonini & Rasio 2016), either early during its evolution or at later points if an MBH is temporarily ejected via gravitational wave recoil (Gualandris & Merritt 2008).

The rate suppression due to strong scatterings can depend strongly on the semimajor axis, or equivalently the energy ϵ, of the stars in question. At low values of ϵ, TDE rates can be suppressed by one or more orders of magnitude in the presence of high-mass and steeply cusped sBHs, but the suppression is often milder (especially at high γ_bh) for the subset of stars on tightly bound, high-ϵ orbits. We have seen this numerically, and established it analytically in equations (22) and (23). A consequence of this is that galactic nuclei can see a relative suppression of high-β TDEs, which are overwhelmingly produced in the full loss cone (or ‘pinhole’) regime where rates are suppressed most heavily by strong scatterings. Although the observational implications of high- versus low-β TDEs are not yet clear, high-β events may exhibit faster circularization (Hayasaki, Stone & Loeb 2013; Bonnerot et al. 2016), produce stronger thermal soft X-ray emission (Dai, McKinney & Miller 2015), produce stronger sub-relativistic outflows from the circularization process (Metzger & Stone 2016; Lu & Bonnerot 2020). These observational signatures could be less frequently produced in galaxies with large and steeply cusped populations of sBHs.

We quantify this effect in Fig. 8, which shows f_pinhole, the fraction of TDEs occurring in the pinhole regime, as a function of MBH mass M_• (with and without strong scattering). Without strong scattering, the value of f_pinhole decreases with increasing M_• and drops steeply for M_• ∼ 10⁸ M_⊙, as expected for cusp galaxies. Interestingly, the influence of strong scattering modifies the pinhole fraction in two opposite directions, with the net effect depending on γ_bh. The first effect is already explained in the preceding paragraph: at fixed τ, loss cone flux is more greatly reduced (by strong scatterings) when ϵ is low than when ϵ is high. However, a countervailing force is visible in Fig. 2, where we see the mild but relevant time evolution of the suppression due to strong scatterings. The suppression factor |$\mathcal {F}(\epsilon)/\mathcal {F}_{\rm w}(\epsilon)$| becomes smaller and smaller as τ grows, so when computing an energy-integrated |$\dot{N}_{\rm TDE}$| rate, the bins of large ϵ (with larger τ) will be further along the time evolution of their suppression curves. In Fig. 8, we see that for large γ_bh, the first effect wins out and the pinhole fraction decreases due to ejections. But for γ_bh ≲ 9/4, the second effect wins out and the pinhole fraction sees little net change from strong scatterings (the net change that does exist is a very mild increase in f_pinhole).

4.2 E+A preference

Observations have revealed that E+A and other post-starburst galaxies are significantly overrepresented among TDE hosts (Arcavi et al. 2014; French, Arcavi & Zabludoff 2016, 2017; Law-Smith et al. 2017; Graur et al. 2018; French et al. 2020; Hammerstein et al. 2021). A few solutions have been proposed to explain this preference: stellar overdensities (Stone et al. 2018), radial anisotropies (Stone et al. 2018), SMBH binaries (e.g. Arcavi et al. 2014), and accounting for a complete stellar mass function (Bortolas 2022).

In the framework of our revised loss cone theory, we have shown that steep slopes induce a strong reduction of TDE rates. Motivated by our findings, here we study the impact of stellar overdensities with strong scattering.

E + A galaxies are post-starburst galaxies and hence, in this section, we consider only stars and no sBHs. If we assume that most of the central stars formed impulsively in a starburst, these ultra-steep stellar central density cusps will relax over time towards a steady-state configuration: the Bahcall–Wolf cusp, with γ_⋆ = 7/4 (Bahcall & Wolf 1976). An ultra-steep stellar density profile will result in much shorter relaxation times, especially at small energies, and hence the relaxation to a BW cusp is quicker for steeper star profiles. Following Stone et al. (2018), we define the Bahcall–Wolf radius r_BW as the radius where the post-starburst age t equals the local energy relaxation time:

In Fig. 9, we plot the Bahcall–Wolf radius r_BW divided by the influence radius for different MBH masses. As r_BW expands, energy relaxation erodes the population of the most tightly bound stars, decreasing their density relative to the initial conditions. The density decrease at high ϵ reduces the efficiency of loss cone shielding, as it is the most tightly bound stars which are responsible for the strongest shielding effects. From Fig. 9, we see that this cusp erosion is fastest for small MBHs and slowest for large MBHs. As we do not self-consistently model this cusp erosion in our calculations, we only compute TDE rate reductions at times where this effect would be minor.

In Fig. 10, we present TDE rates that consider the influence of strong scattering, normalized by TDE rates that do not account for strong scattering. These rates are depicted at two different times where the cusp erosion can be disregarded: |$\frac{r_{\rm BW}}{r_{\rm inf}} \le 10^{-3}$|⁸ (see Fig. 9). We find that, for the steep stellar profile: γ ≥ 9/4, accounting for strong scattering can reduce TDE rates by up to an order of magnitude. Additionally, it is worth noting that the effect is more pronounced for the dashed lines corresponding to times t = 100 Myr and t = 150 Myr (the cusp erosion is still negligible for the higher MBH masses). Our results suggest that ultra-steep stellar profiles can effectively self-shield, which makes them a less plausible explanation for the preference of TDEs for E+A galaxies.

5 DISCUSSION AND CONCLUSIONS

In this paper, we have proposed a revised loss cone theory that takes into account both classical weak scattering and close encounters: strong scattering, collisions, tidal captures, and μ-TDEs. These close encounters can occur at non-negligible rates inside the radius of influence, where stellar densities can be higher than anywhere else in the Universe. Below we summarize our key findings:

• For stars with semimajor axis a ∼ r_inf, i.e. the primary sources of TDEs, strong scattering rates are significantly larger than other close encounter rates. For very high orbital energies, collisions become the dominant close encounter. However, those energies have a very small contribution to TDE rates.

• We derive time-dependent analytical solutions of the Fokker–Planck equation with strong scattering for a broad range of physically motivated slopes. This family of analytic solutions, which we can write in closed form for all relevant integer, half-integer, and quarter-integer γ values, is a generalization of the standard Cohn & Kulsrud (1978) f(j) distribution function in the presence of strong scatterings.

• We find that for BH slopes γ_bh ≥ 9/4, the loss cone will be effectively shielded by strong scattering. TDE rates are reduced by a factor ∼2 for γ_bh = 9/4 and by a factor ≥10 for γ_bh = 5/2.

• Strong scattering is more effective at removing stars with a ∼ r_inf compared to stars with a ≪ r_inf. Thus, in addition to shielding the loss cone, strong scattering can also reduce the fraction of TDEs in the pinhole regime, and therefore the fraction of deeply plunging TDEs for γ_bh = 5/2.

• Stellar overdensities have been proposed as one possible explanation for the observed preference of tidal disruption events (TDEs) in E+A galaxies. However, our revised loss cone theory challenges this idea. In classical loss cone theory, ultra-steep stellar slopes increase the diffusion coefficients and, thus, the TDE rates. However, we have shown that they also increase the rates of strong scattering, which can reduce TDE rates. As a result, the net enhancement of TDEs induced by stellar overdensities is milder than previously assumed. This finding makes other possible explanations, such as radial anisotropies (e.g. Stone et al. 2018) or a more complete mass function (e.g. Bortolas 2022), more likely to explain the observed preference of TDEs for E + A galaxies.

In this paper, we have focused on the most dramatic qualitative difference between weak and strong scatterings: ejection. However, stars also have mildler but still strong encounters that can non-diffusively move them to new – but still bound – energies. This effect will be the subject of a future work. The shielding effect we have studied depends on BH slopes that have not been constrained observationally. N-body simulations have shown that sBHs can settle to a slope of γ_bh ≥ 2 (Preto & Amaro-Seoane 2010; Amaro-Seoane & Preto 2011), in agreement with some predictions of strong mass segregation (Alexander & Hopman 2009). It can also be noted that, for an adiabatically growing MBH, the stellar slope may assume steeper values of γ ≳ 2 (Young 1980). Moreover, Generozov et al. (2018) showed that, in the presence of a source term that is continuously depositing new stars or sBHs at small radii, BHs will settle into a density profile with γ_bh = 5/2. Motivated by these studies, we considered different BH slopes with γ_bh ≤ 2.5.

If steep profiles for BHs are common in galactic nuclei, this may explain persistent discrepancies between the high predicted rates of TDEs (Wang & Merritt 2004) and the (relatively) low observed rates of TDEs in low-mass galactic nuclei. While the steep TDE luminosity function (van Velzen 2018) helps to resolve this rate discrepancy, it may be unable to resolve the entire breadth of the discrepancy if the occupation fraction of small MBHs (∼10⁵ M_⊙) is high in dwarf galactic nuclei. Conversely, loss cone shielding may pose a challenge for the ‘overdensity’ solution (Stone & Metzger 2016; Stone & van Velzen 2016; Stone et al. 2018) to the observed post-starburst preference (Arcavi et al. 2014; French, Arcavi & Zabludoff 2016, 2017; Hammerstein et al. 2021) among the host galaxies of TDEs. While an ultra-steep stellar distribution can boost TDE rates up to a point, our work here shows that the gains from increasing γ_⋆ beyond 9/4 may be self-limiting, and that TDE rates may even drop as γ_⋆ increases to very large values.

It is worth mentioning that we considered a monochromatic distribution of stars throughout the paper. However, we did also compute the impact of a Kroupa IMF and found that it increased all TDE rates with and without strong scattering by a common factor of 1.6. It is also possible that shielding effects could be important for the production of hypervelocity stars in binary break-ups (Hills 1988), as these are generally sourced from binaries with galactocentric semimajor axes a ≫ r_infl (Perets, Hopman & Alexander 2007).

To conclude, we emphasize that we considered the 1D Fokker–Planck equation in angular momentum, neglecting energy evolution of the stellar (and sBH) distribution function. This may be important for steep density profiles that are not quasi-steady states. This, together with the effect of strong scattering that can non-diffusively move stars to higher energies will be the subject of a future work.

ACKNOWLEDGEMENTS

The authors would like to thank Cole Miller, Eliot Quataert, and Luca Broggi for helpful discussions. We also warmly thank the anonymous referee for their very useful comments and suggestions. OT gratefully acknowledges the support of the Einstein-Kaye scholarship and of the Israel Ministry of Science and Technology. This work was partially funded by the Israel Science Foundation (Individual Research Grant 2565/19) and the United States-Israel Binational Science Foundation (grant no. 2019772).

DATA AVAILABILITY STATEMENT

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: ORBIT-AVERAGED EJECTION RATE

Throughout this paper, we performed orbit averages of various interaction rates X in the following way:

In the special case of a Kepler potential, this simplifies to

We found that the orbit averaged ejection rate has an analytical form for γ_⋆ values that are integers or half-integers. Hereafter are the analytic forms we found for interesting values of γ. For an equal mass scatterer and γ_⋆ = 3/2, we obtained

For an equal mass scatterer, with γ_⋆ = 5/2, we obtained

APPENDIX B: LOCAL DIFFUSION COEFFICIENT

The local diffusion coefficient we evaluate in equation (35) is given by Magorrian & Tremaine (1999) and Wang & Merritt (2004):

with

and