Ejected and Rejected: Strong Scatterings in Dense Star Clusters

Introduction

The evolution of many-body collisional star systems is a foundational problem in modern astrophysics: most stars are born in the collisional environment of a gas-embedded cluster [30], a substantial minority pass through an open cluster environment, and even denser systems like globular (GCs) and nuclear star clusters (NSCs) are responsible for the production of unique high-energy sources and transients [15, 23, 51, 49, 17, 39, 29] The dynamical evolution of dense star systems can be studied in a first-principles way using direct $N$-body simulations. However, because the cost of these simulations often becomes prohibitive for $N \gg1$, techniques from statistical mechanics can be more directly useful. In this picture, the background population of stars is idealized as a distribution function (DF), and its fluctuations and discreteness can be treated statistically.

One of the most fundamental tools for modeling the evolution of many-body systems is the calculation of diffusion coefficients (DCs); [14, 11], which most often represent the statistical outcomes of many weak, pair-wise two-body scatterings¹. These DCs are calculated from a rate-weighted average of the underlying distribution of velocity kicks, itself determined from the DF. Once computed, DCs can be used in different evolutionary frameworks: for example, the Fokker-Planck equation uses DCs to evolve the entire test star DF, while the Langevin equation merely evolves the orbit of a single test star. In these frameworks, the first and second-order DCs dictate deterministic drift and stochastic diffusion, respectively.

The two-body drift coefficient was first used in astrophysical environments to identify the dynamical friction force [14], which has since been understood to play a major role across many scales of astrophysical dynamics, from mass segregation in globulars [57] and the Milky Way’s Galactic Center [41], to the infall of satellite galaxies into their primary hosts [45]. Likewise, the second-order diffusion coefficients are the drivers of other important astrophysical phenomena, such as core collapse of globulars and the production of loss cone transients [36]. Taken as a set, the DCs are used to calculate the evolution of the DF through the Fokker-Planck equation, for example in the search for steady state, multi-mass solutions [6, 7, 1, 34]. Individual particle orbits can also be simulated through the use of the Langevin equation [53], which is particularly useful for investigating the properties of extreme mass ratio inspirals (EMRIs) in nuclear star clusters (NSCs; [24, 50, 35]).

Despite their evident scientific value, these approaches cannot capture strong gravitational encounters that result in non-diffusive motion through phase space. Although usually rare, a single strong encounter can qualitatively alter the evolution of a particle by sending it to parameter space not immediately accessible to the first and second-order DCs. Previous works focusing on strong scatterings have investigated properties of the underlying kick distributions [22, 19]. This kick distribution can yield further insights by conveniently providing higher-order DCs [9], estimating rates of cluster evaporation from ejections [33], and identifying novel phenomena, such as loss cone shielding [60, 26, 59].

Notably, previous works on strong scatterings have focused either on rates of specific phenomena such as ejections (e.g. [33, 60]) or on a spectrum of outcomes as is quantified by transition probability operators between different states (energy, for example). However, the calculation of this operator is typically complex, which both limits the direct intuition it provides and hampers its extension into higher dimensions (for example, we are unaware of the calculation of a strong scattering transition probability operator in energy and angular momentum space).

In this paper, we instead focus our effort on deriving the differential rate distribution of strong velocity kicks. This differential rate distribution contains the same information content (and quantifies the same dynamical physics) as a transition probability operator [19], but as we shall see, it can be written in closed form in two-dimensional velocity space. The differential rate distribution readily yields higher moments beyond the first and second-order DCs, and upon specifying the geometry of the stellar system, it can also be converted to obtain the rate distribution of strong kicks to energy and angular momentum (AM), yielding further insights into ejections and gravitational wave phenomena.

The work is organized as follows: § 2 outlines the analytic derivation of the two-dimensional differential rate distribution, along with two concrete calculations for globular clusters (§ 2.1) and NSCs (§ 2.2). § 2.3 outlines a discussion of the third moment’s calculation and approximations for the NSC case; § 3 outlines and approximates the recasting of the rate distribution into action space, and introduces applications to EMRI cancellation and generic ejections. § 4 summarizes this work and discusses further applications and extensions.

The Differential Rate Distribution

We consider the interaction between a test object of mass $M$ that undergoes a hyperbolic encounter with a field star of mass $M_{\mathrm{s}}$. The encounter possesses a relative velocity $V = \lvert\vec{V} \rvert= \lvert\vec{v} - \vec{v}_{\mathrm{f}} \rvert$ (with $\vec{v}_{\mathrm{f}}$ being the field star velocity and $\vec{v}$ being the test star velocity, both in the lab frame) and impact parameter $b = \lvert\vec{b} \rvert$. We treat this encounter as inherently local, unaffected by the presence of a broader GC or NSC potential during the scattering. The changes to the velocity of the test star after this encounter can be written [37]

with a hat denoting a unit vector $\hat{V} = \vec{V}/\lvert V \rvert$ and

Here, $\mu= M_{\mathrm{s}}/(M + M_{\mathrm{s}})$, $b_{0}$ is the impact parameter $b_{0} = G(M_{\mathrm{s}} + M)/V^{2}$ that results in a 90-degree deflection, and $G$ is Newton’s constant. The magnitude of $\Delta\vec{v}$ is invariant under rotational and translational coordinate changes (see Appendix § 6 for more explicit vector geometry), and we can write that the magnitude of change in velocity is

In the limit where $b/b_0$ is large, we can write that $b \approx2GM_s/(V\Delta v)$. To construct the differential rate distribution (DRD) of velocity kicks, we estimate the encounter rate $\Gamma$ using the $n-\Sigma-V$ approximation, which is generally valid in collisional star systems (E. Modelevsky et al. 2025). We will begin by calculating a very approximate form of the DRD to build intuition, and afterwards we will estimate the DRD more rigorously. First, we will assume a constant number density of field stars $n_*$ with a uniform relative velocity. The differential rate can then be written

and in the limit of large $b/b_0$,

While the overall scaling of $\mathrm{d}\Gamma/\mathrm{d}\Delta v \propto\Delta v^{-3}$ is evident, we are more specifically interested in the differential rate distribution as a function of changes to the test object’s original direction of motion ($\Delta v_\parallel$), and the motion perpendicular to this ($\Delta v_\perp$). Defining the lab frame such that the object’s original velocity is parallel to the $x$-axis, we can also identify a generic boundary to the kick distribution. From Eq. (1) and Eq. (2), a dot product of the velocity kick vector with the relative velocity vector yields

Upon expanding $\vec{V}=\vec{v}-\vec{v}_f$,

In systems such as globulars or NSCs, there is a maximum velocity a particle can possess before it escapes the system, defined as $v_{\mathrm{esc}}$. By the Cauchy-Schwarz inequality and bounding $|\vec{v}_f|\leq v_{\mathrm{esc}}$ for the generic escape velocity of a system, we have that

At the boundary, we have equality of $|\kappa|=v_{\mathrm{esc}}$; squaring both sides yields the following equation for the maximum kicks $\Delta v_\parallel$ and $\Delta v_\perp$ that can be received:

A figure displaying the changing kick boundary is shown in Figure 1 for a $M=M_s=1M_\odot$ object orbiting a $10^6M_\odot$ supermassive black hole (SMBH), assuming a Keplerian potential. We see that for $v=0$, the kick boundary is symmetric in $\{\Delta v_\parallel,\Delta v_\perp\}$; however, a strong asymmetry emerges as the test particle approaches the escape velocity. In this limit, a strong headwind of particles allows for strong scatterings that can remove a sizable fraction of the test particle’s velocity, or even reverse their instantaneous velocity, but that limits the opportunity for test particles to get “boosted” by strong scatterings along their current direction of motion.

To determine the maximum allowed kick accounting for any scattering geometry, we introduce the parameterization $\Delta v_\parallel=\Delta v\cos(\chi)$ and $\Delta v_\perp=\Delta v\sin(\chi)$ with $0\leq\chi\leq\pi$. Substituting into Eq. (9) yields $\Delta v=2\mu[v_{\mathrm{esc}}-v\cos(\chi)]$. To then find the corresponding maxima and minima for $\Delta v_\parallel$ and $\Delta v_\perp$, we can find extrema of $\Delta v_\parallel(\chi)$ and $\Delta v_\perp(\chi)$. We find that if $v<v_{\mathrm{esc}}/2$, then

Conversely, if $v>v_{\mathrm{esc}}/2$, then

The corresponding perpendicular maxima are

These bounds are generic, as they make no assumptions about the underlying field star velocity profile², number density, or (an)isotropy. Instead, the boundary is set by the maximum velocity a field star can possess, the velocity of our test object, and the masses of the two components of the scattering.

We now perform a more accurate calculation of the rate distribution by using a more generic version of the differential rate, i.e. by improving on Eq. (4). Assuming a spherical system and an isotropic DF, we have that

with $\psi$ a scattering angle (Appendix § 6) and $f(v_f)$ the DF. For the stellar systems we consider, $f$ is a function of the distance to the origin of the system, but we omit writing $f(v_{\rm f},r)$ for brevity. The differential distribution of velocity kicks can then be written

where the $\delta$-functions enforce the specific geometry yielding a given pair of observed kicks, $(\Delta v_{\parallel},\Delta v_{\perp})$. The overtilde labels the kicks $(\widetilde{\Delta v}_{\parallel},\widetilde{\Delta v}_{\perp})$ given by a specific scattering geometry $(b,V,\psi)$. The equations for $\widetilde{\Delta v}_{\parallel}(b,V,\psi)$ and $\widetilde{\Delta v}_{\perp}(b,V,\psi)$ can be found in Eqs. (74) and A6, respectively. Following substantial simplification (which we omit here but can be found in Appendix § 7), the DRD integral can be written as

Here, $\Delta\hat{\mathbf{v}}\cdot\mathbf{v}_{\rm f}=\kappa$ defines a plane in the space of field star velocities, with minimum distance $\kappa$ from the $\mathbf{v}_{\rm f}$ origin. For a DF that can be written solely as a function of the specific energy $E$, we will find it helpful to define

Then, the specific energy can be recast in terms of $\xi$ and $\kappa$

where $\xi=0$ is the minimum distance from the plane to the origin and $\Phi(r)$ is the gravitational potential. We note that we are using the positive-bound convention of e.g. J. Binney & S. Tremaine [11], such that $\Phi>0$ and $E>0$ for bound orbits. Correspondingly, we can define

The disk ends where $E(\xi)=0$, or where $\xi=\sqrt{2E_{*}}$, since we can write that

So then our rate distribution finally becomes

Further simplification of the DRD is not possible without specifying a particular DF. We will now consider the two most relevant DFs for collisional star clusters.

Globular Cluster - King Model

Globular clusters are composed of $\sim10^{4-6}$ stars and are abundant constituents of most galaxies. Observations of Milky Way GCs show that they are, to a reasonable approximation, both spherical and isotropic, an observation which may be explicable given that their relaxation times are typically much less than the Hubble time ([28, 55]).

As globulars are dominated by the self-potential of their own stars, collisional relaxation drives the DF towards an isothermal (Maxwellian) distribution. However, due to tidal effects within the host galaxy, a GC will often truncate at a tidally limited radius, allowing stars to evaporate ([56]). As a result, a lowered Maxwellian is a more realistic model for the DF of a globular cluster; a particular family which is the King model ([27]), with the DF given by

Here, $\sigma$ is the (constant) velocity dispersion and $D$ is a normalization constant related to the density³. Combining Eq. (20) with the DF of the King model, we find the globular DRD

This equation can be used to analytically predict the rates of interesting astrophysical phenomena in globular clusters such as runaway star ejections [13], and possibly also collisions [18]; it may be particularly relevant for the buildup of a central very massive star via collisional runaway in a newborn GC.

A substantial fraction of nuclear star clusters in dwarf galaxies seem to lack central SMBHs [42], and are thus also expected to relax into nearly isothermal spheres, with DRDs describable by Eq. (22). For the remainder of this work, however, we focus our attention on NSCs dominated by the potential of a central SMBH.

NSC - Power Law Density Profile

Consider a spherically symmetric and isotropic NSC with a central SMBH of mass $M_{\bullet}$ and a density profile $\rho=\rho_{\mathrm{inf}}(r/r_{\mathrm{inf}})^{-\gamma}$ up to the radius of influence $r_{\mathrm{inf}}$. At $r=r_{\mathrm{inf}}$, where the enclosed mass equals that of $M_{\bullet}$, the field star density is $\rho_{\mathrm{inf}}$. An Eddington integral yields the DF $f(E)=C_{\gamma}E^{\gamma-3/2}$, with $C_{\gamma}$ being a normalization constant that is a function of the power-law index $\gamma$, which itself must lie below $3$ (at which point the enclosed mass diverges) and above $1/2$ (at which point the velocity field can no longer be isotropic; [61] (J. H. An & N. W. Evans 2006)).

Many different density profiles for NSCs are plausible, as the complex processes of star formation and globular cluster infall (the two main formation pathways for NSCs) do not inherently need to produce power-laws in $\rho$. Nonetheless, power-law density profiles are the natural expectation for NSCs that have had time to collisionally relax since their formation. This result was first seen in the seminal work of J. N. Bahcall & R. A. Wolf [6], which found that $\gamma=7/4$ is a quasi-stationary⁴ state (“zero flux”) solution to the single-mass 1D Fokker-Planck equation in energy space. Subsequent studies that incorporated a multi-mass model (which is crucial for systems with a large spread in mass, such as mixtures of stellar-mass BHs and lower-mass field stars) find that rarer heavy objects tend to form a steeper power-law cusp with slope $\gamma_{\mathrm{H}}$, and light objects form a shallower cusp with $\gamma_{\mathrm{L}}$, leading to mass segregation. In the “weak segregation” regime, the heavy

species dominates energy diffusion rates; here $\gamma_{\mathrm{H}}=7/4$ and $\gamma_{\mathrm{L}}=3/2$ [7]. In contrast, the “strong segregation” regime, in which the heavy species is rarer and thus experiences a more pronounced form of dynamical friction, remains a subject of theoretical controversy. Early work found that the lighter stars maintain in this regime a shallow cusp of $\gamma_{\mathrm{L}}\approx3/2$–$7/4$, while the heavy objects attain extreme slopes, with $\gamma_{\mathrm{H}}\approx7/4$–$11/4$ [1]. The finite-flux strong segregation solutions of T. Alexander & C. Hopman [1] cannot be self-consistent across all radii, however [34], and more recent work has explored zero-flux alternatives [34].

For practical calculations, we take the one-dimensional galaxy-scale velocity dispersion $\sigma_{\mathrm{BH}}$ to be based on the empirical $M$–$\sigma$ relation (K. Gültekin et al. 2009)

The influence radius can then be defined such that $r_{\mathrm{inf}}=GM_{\bullet}(1+\gamma)^{-1}/\sigma_{\mathrm{BH}}^{2}$. The local velocity dispersion at a radius $r<r_{\mathrm{inf}}$ is then $\sigma(r)^{2}=GM_{\bullet}(1+\gamma)^{-1}/r$. Applying Eq. (20) to the NSC DF results in

A plot showing the differential rate distribution as a function of $\Delta v_{\parallel}$ and $\Delta v_{\perp}$ is shown in Figure 2. As the test object traverses its eccentric orbit, the allowed kick distribution changes shape and magnitude. Notably, some kick magnitudes are accessible from a single strong scattering only along individual subsections of the orbit; the most extreme kicks are heavily dependent on the mass ratio between the test object and the field object.

Beyond the bivariate DRD in Eq. (24), we are also interested in the DRD’s marginal distributions, wherein the differential rate can be written as a function of solely $\Delta v_{\parallel}$ or $\Delta v_{\perp}$ (after integrating over the other velocity kick component). Unfortunately, there is no general closed-form solution for these marginal distributions, even for the simple power-law DFs we consider for NSCs here. However, it is straightforward to compute these marginal distributions for the special case of $\gamma=3/2$ (where the math is dramatically simplified by the constancy of the DF), and they are

with the $\mathcal{F}$ and $\mathcal{G}$ functions written explicitly in Eq. (93). The minimum and maximum velocity kicks used as inputs for the $\mathcal{F}$ and $\mathcal{G}$ functions come from solving the corresponding extrema from the kick boundary.

To approximate the marginal distribution for arbitrary $\gamma$, we begin with the exact marginal distributions above, and then hypothesize that the leading-order difference between scattering rates in varied NSCs will not be the shape of the DF but rather its normalization. Since the scattering rate is always proportional to the number density, we now scale the $\gamma=3/2$ result by the ratio $\mathcal{N}(r)$ between the number density for arbitrary $\gamma$ and the number density for $\gamma=3/2$ profiles. As we shall see, this seemingly crude approach provides a surprisingly accurate approximation, which for the parallel kicks becomes

with the density correction factor $\mathcal{N}(r)=n_{\star}(r)/n_{\star}(r)|_{\gamma=3/2}$. At leading order (i.e. the weak scattering limit)

This approximation is notably similar to the even cruder estimate in Eq. (5). Here, we also encounter a similar problem faced by the DCs, namely the “Coulomb logarithm,” denoted $\ln(\Lambda)=\ln(b_{\max}/b_{0})$, with $b_{\max}$ a rough upper bound of the largest viable⁵ impact parameter. Upon taking the second moment in $\Delta v_{\parallel}$ using Eq. (28), we see that $\langle\Delta v_{\parallel}^{2}\rangle\propto\ln(\Delta v_{\parallel,\max}/\Delta v_{\parallel,\min})$. We can then identify that $\Delta v_{\parallel,\min}=2GM_{\mathrm{s}}/(b_{\max}V)$, recovering the ambiguity of the Coulomb logarithm even in this formalism.

To recover the classic form of dynamical friction, we notice that the next highest order term in Eq. (27) scales as $\Delta v_{\parallel}^{-2}$. Upon taking the first moment of Eq. (27), we identify that the leading-order term containing $\Delta v_{\parallel}^{-3}$ vanishes, while the term that scales as $\Delta v_{\parallel}^{-2}$ then itself acquires a Coulomb logarithm factor. Thus, we recover basic properties of the DCs through the DRD, where the first and second moments both contain the Coulomb logarithm.

A similar process can be done for the marginal distribution for perpendicular velocity kicks, with there again being a fully analytic solution for $\gamma= 3/2$ that can be exploited to approximate a general $\mathrm{d}\Gamma/\mathrm{d}\Delta v_{\perp}$:

Since the motion of our particle does not introduce headwinds in the $\Delta v_{\perp}$ direction, multiplying by $2/3$ for isotropy and substituting $V \sim\sqrt{3}\sigma(r)$ into Eq. (5) yields

Examples of the exact marginal distributions and their approximations are shown in Figure 3, where we notice that the heuristic scalings yield good results across a broad range of astrophysically plausible $\gamma$ values, from 1 to 2. This is because although $\gamma$ changes the velocity distribution, it does so relatively weakly. In contrast, the number density of field stars scales as $n_{*} \propto r^{-\gamma}$.

Moments for Nuclear Star Clusters

From our previous results, we identify a heavy tail allowed by the strong scatterings in $\Delta v_{\parallel}$. Part of this information can be captured in the higher-order moments of the DCs. The moment $\langle\Delta v_{\parallel}^{i}\Delta v_{\perp}^{j}\rangle$ can be calculated as

To demonstrate the similarity between this formalism and that of the standard DC calculations (Eq. (94)), we calculated the DCs from both Eq. (31) and Eq. (94) letting $v$ range from 5% to 99.9% of $v_{\mathrm{esc}}$, fixing $\gamma= 7/4$, $r = 0.1\,r_{\mathrm{inf}}$, $M = M_{\mathrm{s}} = 1\,M_{\odot}$, and $M_{\bullet} = 10^{6}\,M_{\odot}$. We find that the relative error between the DRD-derived and standard formalism has a median 0.8%, 1.8%, and 1.8% for $\langle\Delta v_{\parallel}\rangle$, $\langle\Delta v_{\parallel}^{2}\rangle$, and $\langle\Delta v_{\perp}^{2}\rangle$, respectively, across this velocity range. The Coulomb-logarithm ambiguity is dealt with by specifying a minimum velocity kick scale for integration to avoid the divergence at 0. We naturally take $\Delta v_{\min}=2GM_{\mathrm{s}}/(\sqrt{3}\sigma(r)b_{\max})$, with $b_{\max}$ taken to be the same as in the Coulomb logarithm.

Leveraging the same technique of fixing $\gamma= 3/2$ to find closed-form expressions for the marginal distributions, the third moment for parallel velocity kicks can be written

To extend to other values of $\gamma$, the same technique of scaling by the ratio of the number densities provides a sufficient approximation:

We can repeat this procedure for the other non-zero⁶ third moment in the special case of $\gamma= 3/2$:

And once again, performing a similar scaling procedure, we approximate for general $\gamma$ as:

Both of these third moments are presented in Figure 4, where we see that the approximations derived above (solid lines) provide an excellent scaling to the numerically calculated exact third moments (circles). The actual values between both third moments are similar:

For velocities $0 \leq v \leq v_{\mathrm{esc}}$, the ratio of the third moments is bounded

While third moments of the DRD are generally less important for overall rates of energy transport in collisional star systems than are low-order moments (J. Binney & S. Tremaine [11]), they can result in anomalous diffusion effects (B. Bar-Or et al. [9]). More speculatively, by introducing a skewness into the DRD, they may eventually be helpful at regulating numerical issues in studying diffusion in coordinate systems where particular variables have mathematically imposed boundaries, for example diffusion of the Carter constant for Kerr metric geodesics (Y. Alush et al. [2]).

Action Space and Astrophysical Applications

Although working in velocity space has inherent benefits due to the simplicity of the DRD, many types of astrophysical processes instead depend on integrals of motion of the test object, such as energy or AM. For an NSC, sufficiently large energy changes dictate ejections from the NSC; likewise, tidal disruption events and extreme mass ratio inspirals are driven (or throttled) by AM evolution.

To express the previously derived rate distribution in terms of the specific energy and AM, it would be ideal to use the Jacobian of the map $(\Delta v_{\parallel}, \Delta v_{\perp}) \rightarrow(\Delta E, \Delta J)$. In practice, this is substantially complicated by the vector nature of AM. However, in an NSC, the typical orbits relevant for ejections, tidal disruption, and EMRIs are all highly eccentric orbits, with specific angular momentum $J$ being a minute fraction of the circular angular momentum $J_{\mathrm{circ}}$. In this section, we will perform a series of approximations appropriate for this highly eccentric limit, deriving simplified but tractable rates of differential action space kicks that can then be straightforwardly applied to problems in NSC astrophysics where strong scatterings may be relevant.

For such an eccentric object with $1-e \ll1$ in an NSC, there are two regimes of motion: the pericentric regime, comprised of motion that is nearly radial and $v/v_{\mathrm{esc}} \approx1$; and the apocentric regime, where $v/v_{\mathrm{esc}} \ll1$, which is where the object spends most of its time. For the radial part of the orbit, perpendicular velocity kicks are nearly orthogonal to the radius vector, and so

Then, the mapping $(\Delta E,\Delta J)\rightarrow\Delta v_{\parallel}$ is characterized as

and thus under this approximation, the Jacobian map is straightforward,

Thus, given some $(\Delta E,\Delta J)$, we can map to velocity space and then compute the rate densities for all valid velocity coordinates. Explicitly, the transformation becomes

The subscript $i$ allows for the two different branches of $\Delta v_{\parallel}$ allowed by the quadratic nature of the variables. Due to the trivial mapping between $\Delta v_{\perp}$ and $\Delta J$ under our assumptions, the results from the boundaries of the velocity kick distribution dictate that

as $\lambda\approx-1/2$ for $v\approx v_{\mathrm{esc}}$. For the energy, a similar approach can be taken to the velocity kick boundary, although the maximum energy kick need not lie along the velocity boundary itself due to the more complex mapping. Instead, expressing $\Delta v_{\parallel}(\chi)=\Delta v_{\eta}\cos(\chi)$ and $\Delta v_{\perp}(\chi)=\Delta v_{\eta}\sin(\chi)$, where $0\leq\eta\leq1$ allows us to maximize over both $\chi$ and $\eta$. Doing so, we find that the maximum specific energy changes are (with $v\sim v_{\mathrm{esc}}$)

and that the corresponding minimum energy kicks are

The 2D DRD in action space and the corresponding marginal distributions are shown in Figure 5, where geometrically complex distributions emerge from the nonlinear Jacobian factor as $\mu$ varies. The specific results shown are for a particle traveling on an orbit with $a=0.1r_{\mathrm{inf}}$, with $e=0.999$ and the location chosen such that $v/v_{\mathrm{esc}}=0.9$ (i.e. a part of the orbit moderately into the “pericentric” regime of motion). For large test objects, there is a strong preference for strong encounters with a light field object to more tightly bind them. There is a clear emergence of approximate power law behavior for action space kicks in the marginal distributions, up to the cutoffs set by the aforementioned kick boundaries. This quasi-power-law behavior is expected to be inherited from the velocity space DRD, since the mapping from $\Delta v_{\perp}$ to $\Delta J$ is simple, and also (for the energy changes at rapid velocity and weak kicks) $\Delta E\approx-v\Delta v_{\parallel}$.

To understand the validity of the approximations we have introduced by assuming a radial orbit, we also show in Figure 5 data points with the exact marginalized distributions. These are computed by Monte Carlo sampling of Eq. (13), where we convert the drawn velocity kick variables to changes in energy and AM, with no assumptions made that the orbit is radial. We see that for the marginal distributions in energy, the agreement is excellent. Our approximations work less well for changes in AM. This can be attributed mainly to the fact that our prescription $\Delta J\approx r\Delta v_{\perp}$ implies

where $J_i$ and $J_f$ are the magnitudes of the AM before and after a scattering, respectively. As the product of the tangential velocity and the radius is conserved through an orbit, the above equation holds best at apocenter, when the tangential velocity is minimized, and vector addition can be approximated by scalar addition in the strong regime. Even in this case, we find that relative errors between the Monte Carlo draws for the marginalized AM and our radial approximation have a median of $\sim100\%$, with the error being dominated in the weak regime. For the same parameters but instead $1-e = 10^{-4}$, we find the median error drops to $\sim50\%$, and for $1-e = 10^{-5}$, the median relative error falls to be $\sim20\%$.

These errors are larger than one would ideally hope for, but are still modest enough to enable many astrophysical applications, particularly for the near-radial orbits characteristic of loss cone problems. We find the energy is well approximated in the radial regime, as we made no underlying approximations about its behavior.

Canceled EMRIs

EMRIs are a key science case for the upcoming Laser Interferometer Space Antenna (LISA; P. Amaro-Seoane et al. [3]), where they will provide a particularly stringent test of general relativity (L. Barack & C. Cutler [10]). An EMRI is expected to spend many orbits inside the LISA bandpass before plunging into the SMBH, building up incredible signal-to-noise [5]. However, a key uncertainty for LISA is the EMRI rate, as it is unknown whether there will be tens or thousands of EMRIs in the LISA band at any given moment. We aim below to understand the plausibility of derailing a proto-EMRI from ever reaching the LISA band, and the corresponding rate at which this “cancellation” may occur.

Typically, an EMRI’s evolution is governed by two competing regimes: the timescale for GW-inspiral ($t_{\mathrm{GW}}$) between the test object and the SMBH; and the AM diffusion timescale ($t_{\mathrm{AM}}$) caused by the sea of field particles around the test object. When $t_{\mathrm{GW}} \gg t_{\mathrm{AM}}$, the particle exhibits a random walk in $(1-e,a)$ space, potentially crossing into the GW-dominated regime, becoming circular, or plunging into the SMBH without an appreciable GW-inspiral [24, 50, 35].

If the particle crosses sufficiently into the GW-dominated regime, $t_{\mathrm{GW}} \ll t_{\mathrm{AM}}$ and the particle will first circularize from GW at nearly fixed pericenter, and then decrease its semi-major axis until it inspirals into the SMBH [47]. Late in this process, the particle will enter the LISA band and provide GW signals that will be key to LISA’s broader science goals.

In this standard picture, EMRIs totally decouple from stellar perturbations once $t_{\mathrm{GW}} \ll t_{\mathrm{AM}}$. However, one can hypothesize that although AM diffusion appears to be subdominant, there is still a rare chance of a close encounter at apocenter that might send the particle from the GW-dominated regime to the AM-dominated regime, canceling its subsequent inspiral.

Consider a proto-EMRI which is still highly eccentric and has yet to enter the LISA bandpass, but is nevertheless clearly on the EMRI track ($t_{\mathrm{GW}} \ll t_{\mathrm{AM}}$). At this stage, $v/v_{\mathrm{esc}} \ll1$ at apocenter, so the asymmetries in the kick distribution caused by the particle’s motion are minute. Thus, we can consider the two-dimensional velocity plane normal to the radius vector at apocenter. Although kicks in this plane have a geometry which is technically different from the parallel and perpendicular directions we modeled previously, in practice they will behave the same as $\Delta v_{\perp}$ (again, due to $v \ll v_{\mathrm{esc}}$). Other kicks out of this plane, i.e. parallel or antiparallel to the radius vector, do not contribute to the AM. In other words, at apocenter

From Eq. (30), the marginal distribution of AM kicks can be expressed as

The critical AM $J_c$ at the boundary between the GW and AM-dominated regimes is found by setting $t_{\mathrm{GW}} = t_{\mathrm{AM}}$ [24],

Here, $\langle M_s^2\rangle$ is the second moment of the mass distribution — taken to be $(1\,M_\odot)^2$ — and $k \approx0.34$ [11]. The cancellation of a proto-EMRI requires an encounter to deposit $\Delta J \approx J_c$ to the object to move it to the relaxation-dominated regime. The instantaneous canceled EMRI rate $\Gamma_{\mathrm{CE}}$ can be written using Eq. (47):

Here, $\Delta v_c = (J_c-J)/r$, and $\Delta v_{\perp,\max}$ is known from the kick boundary. Since $\Delta J \propto r$, cancellation will occur primarily at apocenter. Since the object spends the vast majority of its time near apocenter for an eccentric orbit, the orbit-averaged rate of the object is just its rate evaluated at apocenter.

Here, since $v\sim0$, $\Delta J_{\max}=2\mu r_a v_{\mathrm{esc}}(r_a)$, and so then our orbit-averaged apocentric cancellation rate is

Notably, all particles that begin on a track of $t_{\mathrm{GW}}=pt_{\mathrm{AM}}$, with $p$ being a dimensionless constant less than one, have the same expected number of cancellation encounters. To understand this relation, we first notice that the AM of the object is approximately constant during the initial stages of GW inspiral, and so the pericentric distance is also approximately constant. From Eq. (48), $J_c \propto a^{\frac{(\gamma-3)}{5}}$. Then, we observe that

Since the AM or pericenter is nearly constant along the high-$e$ phase of the inspiral, we can write that

We have from P. C. Peters [47] that the GW timescale scales as

and thus the total number is

Thus, for a given set of objects that all start along the same line of constant $t_{\mathrm{GW}} = p t_{\mathrm{AM}}$, we have an effectively constant canceled EMRI number. It is the choice of initial AM (and thus pericenter) that allows for such nice cancellation.

The calculation of $N_{\mathrm{CE}}$ across the relevant $(1-e,a)$ space is shown in Figure 6, where we see a large disparity in the expected number of accumulated strong kicks, depending on the location of the start of the GW inspiral. To prevent any potential cancellation being swamped by diffusion of AM, we exclude points where $J_{\mathrm{c}} - J < 5\sqrt{\langle\Delta J^{2}\rangle_{\mathrm{oP}}}$, or where the change in AM needed from one strong encounter is less than five typical diffusion step sizes. Notably, tracks even well within $t_{\mathrm{GW}} = t_{\mathrm{AM}}$ are subject to finite probabilities of strong encounters during their lifetime, although the $N_{\mathrm{CE}} = 1$ curves are quite close to the $t_{\mathrm{GW}} = t_{\mathrm{AM}}$ one.

In this calculation, we did not consider the resulting changes in energy from velocity kicks, which would result in a vertical displacement in $(1-e,a)$ space. However, to a first approximation, a kick at apocenter will likely not change the energy significantly. Consider an object undergoing a scattering at apocenter on an eccentric orbit, where after scattering, the final velocity is a factor $g$ times the initial velocity ($v_{\mathrm{F}} = gv_{\mathrm{I}}$). Then, the ratios of the initial to final energy can be written

with the final term truncated to first order. From the vis viva equation at apocenter,

Thus, for an orbit with $1-e_{\mathrm{I}} \sim10^{-2}-10^{-3}$, $g \lesssim10^{1}-10^{2}$ to stay within an order-unity correction. Thus, if some projection of this velocity kick is in the plane normal to the radius vector, then the AM can change by nearly an order of magnitude while marginally changing the energy.

The net effect of cancellations may not change overall EMRI rates in steady state, where the cancellation phenomenon can be viewed (averaging over many cycles of cancellation and relaxation) as a slowdown in the rate of AM diffusion towards the $t_{\mathrm{GW}} = t_{\mathrm{AM}}$ separatrix. However, a potential consequence of EMRI cancellations is a shift in the distribution of eccentricities observed by LISA. As an object undergoes the GW-dominated evolution, it will cross the regime where $N_{\mathrm{CE}} \sim1$, at which point it will either: (i) be moved back to the AM-dominated regime, potentially wandering back near the GW-dominated regime; or (ii) be placed onto an even more eccentric orbit, where $N_{\mathrm{CE}} \ll1$, and it can safely continue on the GW-inspiral track. Thus, strong scatterings will likely increase the number of eccentric EMRIs entering the LISA band, while suppressing the number of low-$e$ EMRIs, although direct/MC simulations would be needed to verify this hypothesis.

Ejection Rates

The rates of interesting phenomena like TDEs and EMRIs are governed by the distribution of objects near the SMBH and their diffusion into low-AM orbits vulnerable to tidal disruption or GW-driven inspiral. However, processes like ejections can act to throttle the supply of particles into this regime, a concept more generally referred to as “loss cone shielding” (O. Teboul et al. [60]). Ejections act as a leak in the TDE/EMRI supply line, and precise estimates of ejection rates are therefore needed to accurately predict the corresponding rates of loss cone phenomena.

A key advantage of the action space coordinates is their ability to predict outcomes of strong scatterings, such as ejections, which is straightforward once given the marginal distribution $\frac{\mathrm{d}\Gamma}{\mathrm{d}\Delta E}$. During the radial part of the eccentric orbit, the object moves at $v \approx v_{\mathrm{esc}}$ and only needs a slight velocity kick to become unbound from the SMBH. To approximate the rate of ejections, one can utilize the leading order terms in Eq. (28); taking $\Delta E \approx-v\Delta v_{\parallel}$,

An object becomes unbound when $\Delta E + E < 0$, and so the instantaneous local ejection rate is

where $\Delta E_{\min}$ is the local maximum allowed (negative) energy kick. Using Eq. (57) and under the assumption of $E \ll\lvert\Delta E_{\min}\rvert$ along with $v \approx v_{\mathrm{esc}}$, the ejection rate becomes

A more useful expression is the orbit-averaged ejection rate. Using Eq. (98) in the pericentric regime, we have

We display our calculations of the orbit-averaged ejection rates in Figure 7. We see excellent agreement comparing the ejection rates from Eq. (41) and the rates calculated from O. Teboul et al. [60], along with a close approximation from Eq. (60).

However, a key assumption in both the previous ejection rate and the calculation from O. Teboul et al. [60] is that kicks capable of ejecting a test star are rare, such that they will likely not happen more than once per orbit. To investigate when this assumption breaks down, we first define the “saturated strong kick regime”

with $t_c$ being the characteristic timescale over which our actual DRD changes. When this equation is satisfied for $-\Delta E > E$, there is at least one strong kick that can act to eject the particle. The characteristic radius at which the DRD changes is

where we approximated $v \approx v_{\mathrm{esc}}$ and have taken the weak limit using Eq. (57). The characteristic timescale is simply the crossing time of $r_c$

We can then solve for $\Delta E_c$, the critical energy for which we expect to have one encounter during the time $t_c$:

The value of $\Delta E_c$ will fluctuate throughout the orbit, but $E$ will not. The scaling implied by the above equation shows that for $\gamma> 1$, the maximum value of $\Delta E_c$ is at pericenter, and so enforcing $E = \Delta E_{c,\max}$ yields the characteristic semimajor axis beyond which we are in the saturated ejection regime:

Another assumption implicit in this analysis is that diffusion in energy from weak scatterings alone is much slower than the orbital period, namely that

where $\langle\Delta E^2\rangle_o$ is an orbit-averaged energy diffusion coefficient encapsulating the effects of weak scatterings. If energy relaxation happens at or faster than the orbital timescale, strong kicks that might have led to an ejection in isolation can ultimately be swamped by the rapid diffusion from weak encounters, leading to potentially enhanced loss cone shielding. During the radial part of the orbit, where most of the diffusion in energy happens, we can take $v \approx v_{\mathrm{esc}}$; then

the local second-order DC in energy is approximated by (Appendix D)

Here, $B_\gamma= B(5/2,\gamma- 1/2)$, with $B$ being the beta function. We take $\ln(\Lambda) = \ln(0.4M_\bullet/M_{\mathrm{s}})$. Performing an approximate orbit-average with $v_r \approx v_{\mathrm{esc}}$,

To find situations where diffusion can swamp singular ejections ($t_E \lesssim P$), setting $t_E = P$ yields the minimum semimajor axis for the “diffusively saturated regime,”

Because $\mathrm{d}\Gamma/\mathrm{d}\Delta E \propto|\Delta E|^{-3}$, strong scatterings contribute a fraction $\sim1/\ln\Lambda$ to the total time-averaged rate of energy evolution. Outside of the diffusively saturated regime, it is only rare strong scatterings that can eject stars from the system, but in this more exotic regime, even weak scatterings may in principle trigger an ejection over the course of a single orbit, and the greater number of such weak kicks means that they will (for $a > a_{\mathrm{diff}}$) dominate ejection rates.

Figure 8 displays both the saturated ejection regime, where $\Delta E_{c,\max} = E$ (solid lines), and the diffusively saturated regime, where $\langle\Delta E^2\rangle_oP = E^2$ (dashed lines). In these regimes, the single kick ejection approximation criterion likely misrepresents the true ejection rate, and regions possessing one or both of the aforementioned criteria become especially significant for a smaller-mass SMBH. We see in Figure 8 that the diffusively saturated regime is always broader than the saturated strong kick regime, implying that saturation of diffusive energy evolution will occur first as test stars random walk to extreme values of eccentricity.

The exact implications of energy kick saturation (in either form) for loss cone shielding are not trivial, and we defer a detailed examination of this for future work. In the saturated regime, test stars can move very quickly to either (i) lower $|E|$, enabling ejections, or to (ii) higher $|E|$, which may eventually freeze out rapid energy evolution (from both strong and weak scatterings). The first of these processes is itself subtle, as NSCs are embedded in a broader galactic potential, so that even if a test star is kicked to $v > v_{\mathrm{esc}} = \sqrt{2GM_\bullet/r}$, it is unlikely to be unbound from the galaxy itself; however, this crude ejection criterion may suffice in practical terms if the stellar potential outside the NSC is sufficiently aspherical [58]. The net result of the competition between these two processes is not obvious, and in principle could either greatly strengthen loss cone shielding, if process (i) predominates, or weaken it, if process (ii) does. Notably, however, Figure 8 shows that saturation of energy kicks occurs in precisely those parts of parameter space (large $\gamma$, large $a$) where loss cone shielding through rare strong scatterings is itself most effective [60, 59].

Saturation of energy-space evolution also calls into question the accuracy of orbit-averaged schemes to reliably estimate other loss-cone processes around a low-mass SMBH. “Cliffhanger” EMRIs are EMRIs that originate in the $t_{\rm GW} \gg t_{\rm AM}$ regime, far above the canonical “critical semimajor axis” that separates true EMRIs from failed plunges [24], but which eventually attain $t_{\rm GW} \lesssim P$. They are thus able to dissipate such a significant amount of energy in one orbit that they may move non-diffusively into the $t_{\rm GW} \lesssim t_{\rm AM}$ regime [50], hanging on and becoming true EMRIs. Interestingly, there appears to be a discrepancy in the rates of cliffhanger EMRI success between orbit-averaged and local Monte Carlo schemes for solving the Langevin equation, especially at lower SMBH masses [35]. The origin of this discrepancy, in which orbit-averaging reduces the cliffhanger rate (i.e. increases the plunge rate), remains unclear, but we speculate⁷ that this could be attributed to the incredibly rapid relaxation in energy that may often coincide with the cliffhanger regime, see Figure 8. As ejected objects cannot contribute to the EMRI or plunge rate, they are usually neglected in EMRI statistics for Langevin-type simulations [1, 50, 35], and so regardless of the net direction of energy evolution, the submersion of the $t_{\rm GW} < P$ zone inside the diffusively⁸ saturated kick regime is likely to bias EMRI statistics towards fewer plunges, and more apparent cliffhangers.

Overall, obtaining accurate ejection rates in the regimes where $\langle\Delta E^2\rangle_o P \gtrsim E^2$ or where $\Delta E_{\mathrm{c,max}} \gtrsim E$ using an approach that only considers rates of singular kicks is unlikely to be accurate. Instead, N-body simulations, local Langevin equation integrations, or schemes that directly sample the full DRD are likely the best way to obtain the ejection rate when energy diffusion is rapid or when the ejection kick regime is saturated.

Conclusion

In this work, we have derived the explicit differential rate distribution in velocity space for the two prototypical examples of collisional star systems: globular and nuclear star clusters. We have also computed the explicit boundaries of this DRD. For the NSC, we determined approximate formulae for the marginal distributions of kicks and also estimated the third moments of the overall kick rate distribution. We then recasted the velocity kicks into action space using a series of approximations appropriate for loss cone problems in NSCs, and applied this framework to study loss cone shielding via ejections, and the novel phenomenon of EMRI cancellation.

Past work on the collisional dynamics of dense star systems has almost always operated in the limit of weak encounters and diffusive evolution [57, 11]. The motivations for this choice are clear enough: weak encounters dominate overall rates of energy evolution by a factor of $\sim\ln\Lambda$. However, even if uncommon strong scatterings are subdominant drivers of the bulk evolution of their host cluster, they may unlock qualitatively novel phenomena forbidden to even a lengthy sequence of weak, diffusive encounters. The classic example of this kind of non-diffusive effect is stellar ejection [33], which has recently seen a revival of interest as a mechanism for shielding SMBH loss cones [60] and suppressing rates of transients such as TDEs [59] and EMRIs [26]. Other non-diffusive types of evolution, such as the cliffhanger EMRI phenomenon, have also been the topic of recent interest [50, 35].

Past efforts to systematically study the effects of strong scatterings have usually operated either with transition probability operators [19, 9] or by making precise rate estimates for a single process enabled by strong scatterings, such as ejection [60]. While transition probability operators are a completely general formalism for handling non-diffusive evolution, their mathematical complexity has so far restricted their application to 1D problems in dense star cluster dynamics (i.e. energy-only evolution). Our calculation of an explicit differential rate distribution encompasses the exact same physics as a transition probability operator, but is more analytically tractable by working entirely in the differential (rather than integral) mode. As a result, we are able to compute exact bivariate DRDs in velocity space; while we did not find closed-form translations of the DRD into action space, we present useful approximations, and note that this translation is straightforward to accomplish numerically.

For ejections of high-$e$ stars, we can use the DRD to recover the same orbit-averaged ejection rates as [60], but we also identify new regimes where the assumption of strong scatterings being rare breaks down. First, we identified a regime where strong scatterings can happen more than once per orbit, likely leading to a Poisson-like regime for determining ejection rates. Second, we found that this regime, in which strong kicks saturate, is likely encompassed and superseded by a wider regime in which even weak scatterings lead to non-diffusive evolution of energy. This broader regime of diffusive energy-kick saturation (where $t_{\mathrm{E}} < P$) may feature surprising implications for loss cone feeding rates, by either strengthening or attenuating loss cone shielding. We did not quantify these implications in this work, but we did map out the parameter space where they are likely relevant, and found that diffusive kick saturation applies to much of the parameter space ($a \sim r_{\mathrm{inf}}$, $\gamma\gtrsim2$) where standard loss cone shielding has been found to be important [60].

In applying the DRD formalism to EMRIs, we found that a sizable fraction of EMRIs with $t_{\mathrm{GW}} < t_{\mathrm{AM}}$ can undergo a scattering that sends them back into the AM-dominated regime. The cancellation of these EMRIs will preferentially select against the most circular EMRIs, biasing the loss cone channel for EMRI formation even further towards high-eccentricity signals, and increasing its distinguishability against other formation scenarios, such as Hills separations [38, 52] and wet EMRIs [32].

It is important to note the assumptions inherent in our analyses. The largest assumption was the choice of a simplistic DF, as most of this paper assumed a DF that was single-mass, spherically symmetric, and isotropic, none of which fully describe the Milky Way’s NSC [16, 48]. In principle, our formalism can be generalized straightforwardly to multi-mass or anisotropic systems; asphericity would introduce greater complications [62]. During the scattering process, we have also omitted tidal [31] and relativistic [43] effects, which can act to dissipate energy throughout particularly close encounters. Regarding the canceled EMRI number calculation, a more correct approach would be to fully integrate the rate throughout the expected GW-inspiral track [47], potentially even using an MC approach to include both DC and strong scattering effects.

There are likely many other aspects of collisional star clusters that could be studied with the DRD derived in this work. Particularly likely applications include rate calculations for interesting dynamical encounters, including the production of hypervelocity stars in pair-wise scatterings [44], close encounters producing LIGO-band gravitational waves through capture of compact objects [43], or depletion of dark matter particles in the center of galaxies [54]. While we have examined some implications of the DRD for TDEs and EMRIs, others likely exist as well, such as occasional refill of even deeply empty loss cones [64], which carries implications for both the formation of repeating partial TDEs [12] and the eccentricities of EMRIs.

Acknowledgments

AD and NCS are happy to acknowledge helpful discussions with Yael Alush during early stages of the project.

AD gratefully acknowledges support from the Peter Livingston Scholars Program at the University of Wisconsin–Madison, made possible by the generosity of Sharon Stark in memory of her husband Peter Livingston. AD acknowledges support from the Fannie and John Hertz Foundation. Support for NSC and AD was provided by the University of Wisconsin–Madison, Office of the Vice Chancellor for Research and Graduate Education, with funding from the Wisconsin Alumni Research Foundation. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-2140743. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

AI Disclosure Statement: The authors implemented generative AI (Claude Opus 5.5) in the use of text refinement suggestions and optimization of code written by the authors.

Software: matplotlib [25], numpy [21], scipy [63], pandas [46].

Scattering Geometry

Consider the lab frame (unprimed) and the encounter frame (primed). The encounter frame defines the velocity changes relative to the initial relative velocity vector, and one must then project the velocity changes in the encounter frame to that of a chosen lab frame. In the lab frame, we choose to orient our unit vector in the x-direction ($\hat{e}_x$) along the direction of our test object’s original velocity, so that $\vec{v} = v\hat{e}_x + 0\hat{e}_y + 0\hat{e}_z$.

In the encounter frame, a natural basis can be defined using the unit vector of the relative velocity $\hat{V}$ and impact parameter vector $\hat{b}$. Then, we can write $\Delta\vec{v} = \Delta v_V\hat{V} + \Delta v_b\hat{b}$. To write the relative velocity unit vector in the lab frame

With $V_x = v-v_{fx}$. Defining $\alpha= V_x/V$ and $\beta= -\sqrt{1-\alpha^2}$ and introducing a variable $\phi$ as the azimuthal angle around $\hat{e}_x$ such that $\hat{e}_\phi= \cos(\phi)\hat{e}_y + \sin(\phi)\hat{e}_z$, we see

To construct the basis in the primed frame, first define $\hat{e}'_x \equiv\hat{V}$. Then, define another basis vector in the plane spanned by $\hat{e}_x$ and $\hat{V}$. Guessing a form of $\hat{e}'_y = c\hat{e}_x + d\hat{e}_\phi$, one can determine $c=-\beta$ and $d=\alpha$ suffices as another basis vector. Thus,

To construct $\hat{e}'_z$, take the cross product of $\hat{e}'_x \times\hat{e}'_y$ using the lab basis components, finding $\hat{e}'_z \equiv\hat{e}'_x \times\hat{e}'_y = 0\hat{e}_x - \sin(\phi)\hat{e}_y + \cos(\phi)\hat{e}_z$. To then express Eq. (1) in terms of the primed basis, introduce an azimuthal angle in the scattering frame $\psi$ such that $\hat{b} = \cos(\psi)\hat{e}'_y + \sin(\psi)\hat{e}'_z$. The geometry of both primed and unprimed frames is shown in Figure 9. We can then simplify to write

One sees that

For the magnitude of the perpendicular component, we have

Which, as one would expect, is independent of $\phi$, although the quantities $\Delta v_{\parallel}$ and $\Delta v_{\perp}$ depend on $\psi$. We have that the magnitude of the velocity kick yields

Thus, the magnitude of the kick is invariant under our coordinate change.

Kick Distribution Derivation

From Eq. (13), first make the variable change from $\{v_{fx},v_{fy},v_{fz}\} \to\{V_x,V,\theta\}$. Then, a Jacobian transformation results in $d^3v_f = V\,dV\,dV_x\,d\theta$. Assuming $f$ is isotropic

We have then that

Where we have already integrated over $\theta\in[0,2\pi)$. Then, let $\psi_*$ to be where the left $\delta$ function is zero. This requires

and so defining $\Omega^{1/2} = \sqrt{\Delta v_b^2\beta^2-\left(\Delta v_V\alpha-\Delta v_{\parallel}\right)^2}$ we have

where $\Theta$ is the Heaviside function. We adopt the convention with $\Theta(0)=0$. To evaluate the remaining delta function, it is convenient to change variables. Let $u=b/b_0$, with $b_0=G(M+M_s)/V^2$. Then, performing this change of variables results in

To perform the integral involving the second $\delta$-function, define $u_*$ to be the value of $u$ where the remaining $\delta$-function argument is zero. Then $u=u_*$ implies

The argument of the root requires $4\mu^2V^2\ge\Delta v_{\parallel}^2+\Delta v_{\perp}^2$. For $u\in[u_{\min},u_{\max}]$,

Enforcing these restrictions on $\Omega$ we have

Then, we have

Without loss of generality, take the perpendicular kick direction to be aligned with the y-axis, so that $\Delta\vec{v}=\Delta v_{\parallel}\hat{e}_x+\Delta v_{\perp}\hat{e}_y+0\hat{e}_z$. Now, take a general relative velocity $\vec{V}=V_x\hat{e}_x+\sqrt{V^2-V_x^2}\cos(\theta)\hat{e}_y+\sqrt{V^2-V_x^2}\sin(\theta)\hat{e}_z$. Then,

Define a quantity $h=h(\theta)$ so that

Notice that the integration of $h(\theta)$ yields

But then, notice that

Here, the Heaviside functions have been omitted, which is permitted in the case of large $u_{\max}$ and $u_{\min}=0$. Further, $2\mu V>\Delta v$ from Eq. (6) and the Cauchy-Schwarz inequality. The definition of $\kappa$ allows us to rewrite the scattering condition as

Inserting this allows for the expression of the integral in Eq. (85) as

This is a two-dimensional integral, and can be interpreted as intersecting the three-dimensional field star velocity space with the plane where $\Delta\vec{v}\cdot\vec{v}_f/\Delta v=\kappa$ is enforced. This intersection then has minimum distance $|\kappa|$ from the origin, and conveniently traces out a disk to integrate over if the field star distribution is isotropic. A similar form and interpretation is derived in J. Goodman [19]. Generically, we can write Eq. (85) instead as

Where $\mathrm{d}A$ is the area element on the plane in the allowed region where $|\kappa|\leq v_{\mathrm{esc}}$. Then, one can introduce a radial coordinate in the disk to perform the area integral.

Marginal Distributions

The fully explicit forms of the $\mathcal{F}$ and $\mathcal{G}$ functions are

Diffusion Coefficients

For the velocity-space DCs, we have [37]

Here, $F_n(v)$ and $E_n(v)$ are defined as

To calculate the local second-order DC in action space, one has

In the NSC discussed in this work, when computing $\langle\Delta E^{2}\rangle$ for the object in the radial part of its orbit, $E_{1}\sim0$ as the bounds of integration are between $v$ and $v_{\mathrm{esc}}$. Correspondingly, the bounds on $F_{4}$ are approximately between 0 and $v_{\mathrm{esc}}$. Then, we may write

Where the last integral is the beta function $B(5/2,\gamma-1/2)\equiv B_{\gamma}$. Then, we can express $\langle\Delta E^{2}\rangle$ as in Eq. (67).

Orbit Averaging

To average a general quantity $X$ over a Keplerian orbit,

For the case of the highly eccentric orbit, one can approximate that $v_{r}\approx v_{\mathrm{esc}}$ on the radial portion of the orbit, and $v\sim0$ at apocenter.