© The Authors 2025

1. Introduction

Most galaxies host a massive black hole (MBH) at their core, with masses in the range of 10⁵–10¹⁰ M_⊙ (Kormendy & Richstone 1995; Magorrian et al. 1998; Kormendy & Ho 2013). Often these MBHs are found at the centre of nuclear clusters (Carollo et al. 1997; Böker et al. 2002; Côté et al. 2006; Schödel et al. 2007, 2009; Turner et al. 2012): highly dense clusters of stars and compact objects (COs) observed in the nuclei of many galaxies. In the central region of nuclear clusters, densities can reach up to 10⁷ M_⊙pc⁻³ (Neumayer et al. 2020), facilitating the scattering of COs onto tight and eccentric orbits around the central MBH, thus forming binaries with extreme mass ratios (Binney & Tremaine 2008; Merritt 2013a).

General relativity (GR) predicts that these binary systems emit gravitational waves (GWs), mostly during each CO pericentre passage (Misner et al. 1973). Over time the CO orbit shrinks and circularises (Peters 1964) as a consequence of the loss of energy and angular momentum of the binary system, which are carried away by GWs (Maggiore 2007). Accordingly, the CO slowly inspirals towards the MBH, completing thousands of orbital cycles, eventually merging with it (Amaro-Seoane 2018). These GW sources are known as extreme mass ratio inspirals (EMRIs) due to the large mass imbalance between the pair of objects that make up the binary system: the MBH and the orbiting CO (Amaro-Seoane et al. 2007; Amaro Seoane 2022). Alternatively, if the merger occurs after few orbits or in a head-on collision, the event is referred to as a direct plunge (DP).

Extreme mass ratio inspirals that form around MBHs with masses of the order of 10⁴ − 10⁷ M_⊙ are expected to emit GWs in the millihertz frequency range (Barack & Cutler 2004; Amaro-Seoane 2018). The numerous bursts of GWs emitted during an EMRI are anticipated to build up a significant signal-to-noise ratio (Barack & Cutler 2004; Babak et al. 2017), enabling the detection of such events by the upcoming Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017, 2023; Colpi et al. 2024). Detecting EMRIs will enable us to test GR (Ryan 1997; Barack & Cutler 2007; Ghosh & Chakravarti 2024), estimate cosmological parameters (MacLeod & Hogan 2008; Laghi et al. 2021), and investigate dormant MBHs (Gair et al. 2010; Gallo & Sesana 2019).

Recent works (e.g. Metzger et al. 2022; Franchini et al. 2023; Linial & Metzger 2023) have also proposed that EMRIs might be responsible for the novel, currently unexplained, X-ray variability phenomenon of quasi-periodic eruptions (Miniutti et al. 2019; Giustini et al. 2020; Arcodia et al. 2021, 2024), and recent observations agree with such an interpretation (Nicholl et al. 2024). This would offer a way to detect and study EMRIs within the electromagnetic spectrum, enabling multi-messenger observations that would boost their potential as astrophysical and cosmological probes (e.g. Toubiana et al. 2021).

Extreme mass ratio inspirals can form in many different ways, including (i) the tidal separation of a stellar-mass black hole (BH) binary system by the central MBH followed by the capture of one of the binary components (Hills 1988; Miller et al. 2005; ii) the formation or capture and subsequent migration of BHs within accretion discs surrounding MBHs (Levin 2007; Pan & Yang 2021); (iii) the capture of a BH remnant by an MBH as a consequence of supernova kicks (Bortolas et al. 2017; Bortolas & Mapelli 2019); and (iv) the scatter of BHs onto tight orbits due to the gravitational influence of a secondary MBH in the aftermath of galaxy mergers (Chen et al. 2011; Bode & Wegg 2014; Naoz et al. 2022; Mazzolari et al. 2022).

The standard and ubiquitous formation channel for EMRI formation, which we explore in this work, is the scatter of a CO onto a sufficiently tight and eccentric orbit as a consequence of two-body relaxation (Alexander 2017; Amaro-Seoane 2018). The estimated EMRI formation rate through this channel is in the range of 10⁻⁸ − 10⁻⁶ EMRIs per year for a central MBH of mass M_• = 10⁶ M_⊙ (Hils & Bender 1995; Gair et al. 2004; Rubbo et al. 2006). Considering also the uncertainty in the number density of MBHs inside the observable range of LISA and the fraction of these that are hosted in nuclear clusters, the predicted number of EMRIs that will be detected by LISA during its lifetime spans three orders of magnitude, from a few events per year to a few thousand (Babak et al. 2017).

Simulating the formation of EMRIs in this standard channel poses considerable challenges. The most accurate description of the whole nuclear cluster would involve a direct N-body computation of the orbits of all stars and COs within it. However, the need to follow these objects for thousands of orbits, and the high mass ratio between them and the central MBH pose an insurmountable challenge for N-body approaches, which show problems even in the simplest Newtonian regime (Amaro-Seoane 2018). Alternative methods, such as solving the steady-state Fokker-Planck equation for the distribution function of the system (e.g. Amaro-Seoane & Preto 2011; Bar-Or & Alexander 2016; Vasiliev 2017; Broggi et al. 2022; Kaur et al. 2024) or employing Monte Carlo codes (e.g. Hénon 1971; Giersz 1998; Joshi et al. 2000; Freitag & Benz 2001; Pattabiraman et al. 2013; Fragione & Sari 2018; Zhang & Amaro Seoane 2024), offer faster approximated solutions, but carry inherent limitations and assumptions.

From the aforementioned body of work a picture has emerged in which there exists a critical initial semi-major axis a_c for the CO orbit around the MBH determining its eventual fate (Hopman & Alexander 2005; Bar-Or & Alexander 2016; Raveh & Perets 2021). In short, COs starting their inspiral from orbits with a₀ < a_c would only result in the formation of EMRIs, while only DPs would form when a₀ > a_c, with a sharp transition between the two regimes around a_c. Notably, this finding has been derived from detailed simulations of nuclear clusters hosting an MBH with M_• ≥ 10⁶ M_⊙ and then extrapolated at lower masses. Recently, however, this dichotomy has been questioned in Qunbar & Stone (2024) (hereafter QS24) for the case of low-mass MBHs and intermediate-mass BHs, opening to the possibility of EMRIs forming even from initially wide orbits, dubbed cliffhanger EMRIs. This is an intriguing claim that might bear interesting consequences for estimating LISA EMRI detection rates.

As mentioned above, the techniques and codes employed to verify both the existence of a_c and the appearance of cliffhanger EMRIs rely on a number of assumptions and approximations. A notable approximation is that of being allowed to average the effects on the orbit of two-body relaxation over an entire orbital period (e.g. Cohn & Kulsrud 1978; Hopman & Alexander 2005, QS24). However, this simplification might not accurately describe the system, as two-body encounters strongly depend on the instantaneous velocity of the orbiting CO (Binney & Tremaine 2008; Merritt 2013a), which significantly changes during a highly eccentric orbit. Another is to use Newtonian dynamics to describe the system. GR effects are usually included only by adding extra orbit-averaged loss terms due to GW emission and by adjusting the CO capture condition by a suitable factor.

In this work we present a novel post-Newtonian (PN) Monte Carlo code that evolves the PN dynamics of a CO orbiting an MBH. The code treats two-body relaxation locally, by applying multiple small kicks to the CO orbit within a single orbit, and can account for an arbitrary number of stellar components contributing to relaxation and to the potential in which the CO is evolving. As a first application, here we apply this code to study the fraction of EMRIs relative to DPs that form around MBHs of different masses as a function of the initial semi-major axis of the orbit a₀. This allows us to verify the existence of a_c and its dependence on the system parameters, and the emergence of the novel population of cliffhanger EMRIs found in QS24.

The study presented in QS24 has some limitations which our code can improve upon. It employs a Monte Carlo approach to evolve in time the Newtonian energy and angular momentum of an orbit subjected to two-body relaxation, which is orbit-averaged over more than a full orbit at times. The loss of energy and angular momentum due to GW radiation is accounted for through equations detailed in Peters (1964), which are orbit-averaged as well. Moreover, only two-body encounters between a stellar-mass BH of mass m_• = 10 M_⊙ and a single-mass population of m_⋆ = 1 M_⊙ stars surrounding the central MBH are considered. The stars are assumed to follow a simple isotropic power-law numeric density, and their global gravitational potential is not considered. In comparison, in this work we track the orbits of a stellar-mass BH in the total potential well of an MBH and of two populations of objects surrounding it: stars with masses m_⋆ and stellar-mass BHs with masses m_•. We describe both populations with Dehnen numeric density profiles (Dehnen 1993; Tremaine et al. 1994). We also solve the equations of motion up to the 2.5PN term¹, thus naturally accounting for GW radiation, and we locally correct the shape of the orbit for two-body relaxation multiple times within each orbital period.

The paper is organised as follows. In Sect. 2 we review the aspects of the theory relevant to our implementation and predict some behaviours of cliffhanger EMRIs with analytical means. In Sect. 3 we delve into the specifics of our methods, to then present the results of our simulations in Sect. 4. Finally, in Sect. 5 we discuss our findings and draw our conclusions.

2. Theory of relaxation-driven EMRI formation

In this section we first introduce some notions on the treatment of two-body interactions in the loss cone framework. Then, we give operative definitions for EMRIs and DPs. Finally, we briefly review the literature on the study of the EMRI-to-DP ratio and draw some first theoretical results on cliffhanger EMRIs.

2.1. Loss cone and diffusion theory

The main actor in a nuclear cluster is the central MBH, as it strongly affects the motion of stars and COs within its influence radius, which we define in this work as (Peebles 1972)

$$\begin{array}{c}\hfill R_{\inf }=\frac{G{M}_{\bullet }}{{\sigma }^2}.\end{array}$$(1)

Here G is the gravitational constant, M_• is the mass of the MBH, and σ is the velocity dispersion of stars in the bulge of the host galaxy. The velocity dispersion can be related to the mass of the central MBH through the M-σ relation (Ferrarese & Merritt 2000; Gebhardt et al. 2000), that we employ using the best fit from Gültekin et al. (2009):

$$\begin{array}{c}\hfill \sigma =70{\left(\frac{M_{\bullet }}{1.53\times {10}^6{\mathrm{M}}_{\odot }}\right)}^{1/4.24}\mathrm{km}/\mathrm{s}.\end{array}$$(2)

The nuclear star cluster is then composed by many other celestial bodies, which can be grouped into populations based on their nature and mass (e.g. Sun-like stars, stellar-mass BHs, neutron stars, and so on). We start by considering an arbitrary number of populations, later specialising to a system with only two (which can be identified generically as solar-mass stars and stellar-mass BHs).

We label each population with an α index² and model their number distribution with Dehnen profiles (Dehnen 1993; Tremaine et al. 1994):

$$\begin{array}{c}\hfill n_{\alpha }(r)=\frac{3-{\gamma }_{\alpha }}{4\pi }\frac{M_{\mathrm{tot},\alpha }}{m_{\alpha }}\frac{r_{\mathrm{a},\alpha }}{r^{{\gamma }_{\alpha }}{(r+r_{\mathrm{a},\alpha })}^{4-{\gamma }_{\alpha }}}.\end{array}$$(3)

Here m_α is the individual mass of the objects in the α population, while M_tot, α is their total mass. The exponent γ_α, comprised in the range [0, 3), is the slope of the profile for r → 0 and r_a, α is the scale radius at which the distribution smoothly switches to a slope of −4 for r → ∞. Finally, r measures the distance from the centre of the distribution. We consider each profile in a steady-state configuration, with their centres fixed at the origin of our coordinate system, just like the MBH.

The gravitational potential produced by population α is (Dehnen 1993; Tremaine et al. 1994)

$$\begin{array}{c}\hfill {\varphi }_{\alpha }(r)=\{\begin{array}{cc}{\displaystyle -\frac{1}{2-{\gamma }_{\alpha }}(1-{\left(\frac{r}{r+r_{\mathrm{a},\alpha }}\right)}^{2-{\gamma }_{\alpha }})\frac{G{M}_{\mathrm{tot},\alpha }}{r_{\mathrm{a},\alpha }}}\hfill & \text{if}{\gamma }_{\alpha }\ne 2\hfill \\ {\displaystyle ln\left(\frac{r}{r+r_{\mathrm{a},\alpha }}\right)\frac{G{M}_{\mathrm{tot},\alpha }}{r_{\mathrm{a},\alpha }}}\hfill & \text{if}{\gamma }_{\alpha }=2\hfill \end{array},\end{array}$$(4)

so the total potential at a given position is given by

$$\begin{array}{c}\hfill \psi (r)=-(-\frac{G{M}_{\bullet }}{r}+{\displaystyle \sum _{\alpha }}{\varphi }_{\alpha }(r)),\end{array}$$(5)

where we flipped the overall sign adopting the convention that the potential experienced by any object in the distribution is positive definite.

The central MBH affects the distribution of surrounding objects also by removing bodies that pass too close to it at orbital pericentre. For example, stars can be tidally disrupted and subsequently accreted or directly swallowed, whereas COs can be gravitationally captured resulting in EMRIs and DPs. This happens if the velocity vector v of a bound object points towards the central MBH within a small solid angle, inside the so-called loss cone: a cone-shape region in velocity space containing objects which will be lost inside the MBH at the pericentre of their current orbit (see Alexander 2017; Amaro-Seoane 2018). In nuclear clusters, objects can be brought closer to or directly inside the loss cone as a consequence of two-body encounters within each other.

Let us now consider a particular stellar-mass BH of mass m_• = 10 M_⊙ and velocity magnitude v having a sequence of two-body encounters with the background bodies; we decompose its velocity change as the sum of Δv_∥, parallel to the initial line of motion, and Δv_⊥, in the plane orthogonal to the original direction of motion. At any point, given the position r and velocity v of the stellar-mass BH, its orbital binding energy per unit mass can be computed as

$$\begin{array}{c}\hfill \mathcal{E}(r,v)=-(\frac{H(r,v)}{m_{\bullet }}+{\displaystyle \sum _{\alpha }}{\varphi }_{\alpha }(r)),\end{array}$$(6)

where H is the Hamiltonian of the binary comprised of the MBH and the orbiting stellar-mass BH, and we define the energy so that ℰ < 0 for bound orbits. We report the expression of H we used in this work in Eq. (41).

Given the random nature of the scattering process, different encounters will result in different values for Δv_∥ and Δv_⊥. It is for this reason that the problem is better approached statistically, by looking at the diffusion coefficients (i.e. the expectation values per unit of time) of Δv_∥, Δv_⊥, and related quantities. We use the symbol D[X] for the diffusion coefficient of the generic quantity X, and we explicitly write those relevant to our study in Appendix A. Crucially, these are local diffusion coefficients, meaning they depend on the exact position and velocity of the orbiting stellar-mass BH. The usual approach in literature is instead that of employing orbit-averaged diffusion coefficients D_o[X] to account for two-body encounters, which follow from their local versions as (Binney & Tremaine 2008; Merritt 2013a):

$$\begin{array}{c}\hfill {\mathrm{D}}_{\mathrm{o}}[X]=\frac{2}{P}{\displaystyle {\int }_{r_{-}}^{r_{+}}}\frac{\mathrm{d}r}{v_{\mathrm{r}}}\mathrm{D}[X].\end{array}$$(7)

Here r₊ and r₋ respectively indicate the apocentre and pericentre of the orbit, P its radial period, and v_r the radial component of the velocity. If we indicate with a and e the semi-major axis and the eccentricity of an orbit, their Keplerian estimates are

$$\begin{array}{c}\hfill r_{+}=a(1+e),r_{-}=a(1-e)\end{array}$$(8)

and

$$\begin{array}{c}\hfill P=2\pi {\left(\frac{a^3}{G(M_{\bullet }+m_{\bullet })}\right)}^{1/2}.\end{array}$$(9)

We note that the orbit is an ellipse with semi-major axis a and eccentricity e only in the Keplerian approximation, that holds where the stellar potential is negligible (r₊ ≲ R_inf), and the orbit is non-relativistic (r₋ ≫ R_g). Here

$$\begin{array}{c}\hfill R_{\mathrm{g}}=\frac{G{M}_{\bullet }}{c^2}\end{array}$$(10)

is the gravitational radius of the MBH, where c is the speed of light in vacuum. For orbits where the Keplerian approximation does not hold, our definitions of a and e do not have a simple, direct interpretation (see Sect. 2.2 below).

In Appendix A we show that a key ingredient to compute the local diffusion coefficients is f_α(ℰ): the ergodic distribution function of objects in population α. In particular, ℰ is the orbital binding energy per unit mass of the generic α particle. The ergodic distribution function corresponding to a spherical density profile can be computed using Eddington’s formula (Eddington 1916; Merritt 2013a, see also our Appendix A), which simplifies to the expression

$$\begin{array}{c}\hfill f_{\alpha }(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}{\displaystyle {\int }_0^{\mathcal{E}}}\frac{\mathrm{d}\psi }{\sqrt{\mathcal{E}-\psi }}\frac{{\mathrm{d}}^2{n}_{\alpha }}{\mathrm{d}{\psi }^2}\end{array}$$(11)

given our assumptions. For a Dehnen density profile, n_α ∝ ψ⁴ at infinite distance, as n_α ∝ r⁻⁴ and ψ ∝ r⁻¹. Thus,

$$\begin{array}{c}\hfill {\frac{\mathrm{d}{n}_{\alpha }}{\mathrm{d}\psi }|}_{\psi =0}={4{\psi }^3|}_{\psi =0}=0.\end{array}$$(12)

In general, the spherically symmetric distribution function f can depend also on the quantity

$$\begin{array}{c}\hfill \mathcal{R}(\mathcal{E},J)=\frac{J^2}{J_{\mathrm{c}}^2(\mathcal{E})},\end{array}$$(13)

where J is the magnitude of the specific angular momentum of the orbiting body, and J_c is its value for a circular orbit of energy ℰ. Cohn & Kulsrud (1978) showed that

$$\begin{array}{c}\hfill f(\mathcal{E},\mathcal{R})=f_{\mathrm{N}}(\mathcal{E})\frac{ln\mathcal{R}-ln{\mathcal{R}}_{\mathrm{q}}}{{\mathcal{R}}_{\mathrm{q}}-1-ln{\mathcal{R}}_{\mathrm{q}}}\end{array}$$(14)

is a steady-state expression for f near the loss cone. Here f_N is a normalisation factor, while ℛ_q is defined as

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{q}}(\mathcal{E})={\mathcal{R}}_{\mathrm{lc}}(\mathcal{E})exp(-\sqrt[4]{q^2+q^4}).\end{array}$$(15)

In this definition, ℛ_lc denotes the value of ℛ at the edge of the loss cone, while q is the loss cone diffusivity, defined as the ratio between the expected relative change in ℛ² due to two-body relaxation at the edge of the loss cone:

$$\begin{array}{c}\hfill q(\mathcal{E})={\frac{⟨\mathrm{\Delta }{\mathcal{R}}^2⟩}{{\mathcal{R}}^2}|}_{\mathcal{R}={\mathcal{R}}_{\mathrm{lc}}}.\end{array}$$(16)

The q parameter distinguishes between the full (q ≫ 1) and the empty (q ≪ 1) loss cone regimes. In the full loss cone regime, two-body relaxation is particularly efficient in scattering angular momentum, with variations as large as the angular momentum itself. Therefore, objects instantaneously located on an orbit penetrating the loss cone are likely to be scattered away from it before reaching the pericentre, and escape the capture by the MBH. On the other hand, in the empty loss cone regime, once a body enters the loss cone it has very little chance of escaping it, as two-body relaxation induces very small kicks in angular momentum over a single orbit (Cohn & Kulsrud 1978; Merritt 2013a).

2.2. Gravitational wave emission and loss cone captures

The dynamics of a stellar-mass BHs orbiting an MBH in a nuclear cluster is defined by two fundamental timescales:

1. t_rlx, the average time needed for two-body encounters to significantly change the orbital parameters;

2. t_GW, the average time needed for the emission of GWs to significantly change the orbital parameters.

When the stellar-mass BH and the MBH are far enough, the emission of GWs is negligible. In this configuration, the orbit evolution is stochastic and results only from the effect of two-body encounters. These encounters shuffle the position of the orbit on the (1 − e, a) plane and are particularly efficient in changing 1 − e. Observing that ℛ = 1 − e² for Keplerian orbits, we can then define the average time needed for two-body relaxation to significantly change 1 − e as

$$\begin{array}{c}\hfill t_{\mathrm{rlx}}=\frac{1-e}{{\mathrm{D}}_{\mathrm{o}}[\mathrm{\Delta }(1-\mathrm{e})]}\simeq \frac{1-e^2}{{\mathrm{D}}_{\mathrm{o}}[\mathrm{\Delta }(1-{\mathrm{e}}^2)]}=\frac{\mathcal{R}}{{\mathrm{D}}_{\mathrm{o}}[\mathrm{\Delta }\mathcal{R}]},\end{array}$$(17)

where we use the fact that 1 − e² goes like 2(1 − e) for e → 1. The expression we used for D_o[Δℛ] can be found in Appendix A and in Broggi et al. (2022).

During this random motion on the (1 − e, a) plane, the orbit might end up in a situation where the stellar-mass BH passes so close to the MBH that the emission of GWs begins to contribute to the evolution of the orbit, pushing it to become tighter and less eccentric, as the system loses energy and angular momentum. We can measure the strength of this effect by looking at the t_GW timescale, which we define as

$$\begin{array}{c}\hfill t_{\mathrm{GW}}=(1-e){\langle \frac{\mathrm{d}}{\mathrm{d}t}(1-e)\rangle }_{\mathrm{GW}}^{-1}=-(1-e){\langle \frac{\mathrm{d}e}{\mathrm{d}t}\rangle }_{\mathrm{GW}}^{-1},\end{array}$$(18)

where angle brackets indicate orbit-averages. Peters (1964) showed the following:

$$\begin{array}{c}\hfill {\langle \frac{\mathrm{d}e}{\mathrm{d}t}\rangle }_{\mathrm{GW}}=-\frac{304}{15}\frac{G^3{M}_{\bullet }{m}_{\bullet }(M_{\bullet }+m_{\bullet })}{c^5{a}^4}\frac{e}{{(1-e^2)}^{5/2}}(1+\frac{121}{304}{e}^2).\end{array}$$(19)

Thus³,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill t_{\mathrm{GW}}& =\frac{15}{304}\frac{c^5{a}^4}{G^3{M}_{\bullet }{m}_{\bullet }(M_{\bullet }+m_{\bullet })}\frac{(1-e){(1-e^2)}^{5/2}}{e}{(1+\frac{121}{304}{e}^2)}^{-1}\hfill \\ \hfill & \stackrel{e\to 1}{\to }\frac{12\sqrt{2}}{85}\frac{c^5{a}^4}{G^3{M}_{\bullet }{m}_{\bullet }(M_{\bullet }+m_{\bullet })}{(1-e)}^{7/2}.\hfill \end{array}\end{array}$$(20)

At this point, we can define two regions in the (1 − e, a) plane (see Fig. 1):

1. the kick-dominated region, where t_rlx < t_GW and the orbit evolution is stochastic and dominated by two-body encounters;

2. the GW-dominated region, where t_GW < t_rlx and the orbit evolution is quasi-deterministic and dominated by the emission of GWs.

Since strong velocity kicks might still push the orbit back into the kick-dominated region, or directly inside the loss cone after few orbital periods⁴, we caution against identifying the EMRI formation region with the GW-dominated region. Following Hopman & Alexander (2005), a safest choice is that of labelling as EMRIs the orbits that reach the condition t_GW = 10⁻³t_rlx. Here the orbital elements are so deep into the GW-dominated region that we can be confident that a merger will eventually take place after several thousands of periods. We instead define DPs as those events where the orbit enters the loss cone before reaching the t_GW = 10⁻³t_rlx condition.

2.3. Loss cone definition in post-Newtonian gravity

From GR, we know that a test particle on a stable orbit of eccentricity e must always keep a distance from an MBH of at least (6 + 2e)/(1 + e) gravitational radii in order to avoid plunging into it (Cutler et al. 1994). Consistently, the innermost stable circular (e = 0) orbit around an MBH has radius 6R_g, while a body on a parabolic (e = 1) trajectory can approach the MBH up to 4R_g, still managing to escape it (Misner et al. 1973). Therefore, 4R_g can be taken as the loss cone radius for very eccentric orbits in GR.

Relaxation in dense systems, however, is usually modelled with Newtonian gravity, which poses a problem for the loss cone definition. Crucially, there is no single method to translate Keplerian orbits of Newtonian gravity into geodesics of the Schwarzschild metric in GR (Servin & Kesden 2017). In particular, two different mappings can both recover the same Newtonian ellipse for large pericentres. One mapping is that where both the general relativistic and Newtonian orbits share the same pericentre. This mapping has the property of conserving the circumference of circular orbits between the two theories (Servin & Kesden 2017), which is a property that does not appear to be much of use for our study, as we are interested in highly eccentric orbits. Instead, a more useful mapping is the one that preserves the angular momentum between the theories. In this case, we can identify plunging orbits around Schwarzschild MBHs with the same condition used in GR (Merritt et al. 2011; Sadeghian & Will 2011; Will 2012; Merritt 2013a):

$$\begin{array}{c}\hfill J<\frac{4G{M}_{\bullet }}{c}.\end{array}$$(21)

From the Newtonian relation J² = GM_•a(1 − e²), then follows the equivalent condition for the pericentre when e → 1:

$$\begin{array}{c}\hfill r_{-}<8R_{\mathrm{g}}.\end{array}$$(22)

To get our point across let us consider an example. Take a body moving in Newtonian gravity on a highly eccentric orbit, passing at slightly less than 8R_g from the MBH at pericentre and evolve it back in time to the apocentre of its orbit. Given the high eccentricity we assume, the apocentre is large enough that GR simplifies to Newtonian gravity. If we now turn on GR and let the system evolve forward in time, the same body will reach a minimum distance from the MBH slightly lower than 4R_g, plunging into its horizon. Thus, when simulating similar orbits in Newtonian dynamics, we have to consider as plunges all bodies that get to less than 8R_g from the MBH.

However, for this work we neither solved the Einstein field equations to evolve our system nor relied on a Newtonian description of the dynamics. We evolved our systems using PN dynamics, which operates in an intermediate regime by definition. We are not aware of studies on the capture of eccentric orbit at the PN level, unlike in the case of circular orbits (e.g. Schäfer & Jaranowski 2024). Thus, we first test how the PN dynamics affects the position of the pericentre for a highly eccentric and deeply penetrating orbit within our code.

We took a test orbit initialised such that the stellar-mass BH reaches a distance from the MBH of 8R_g at pericentre in Newtonian dynamics. Then, we ran the same initial conditions for three more times, adding a new PN term each time. We repeated this test for different combinations of masses and eccentricities, and found that only the latter have some small influence on the result. In particular, we tested for all combinations of e = { 0.9, 0.999, 0.99999 }, m_• = { 10, 100 } M_⊙, and M_•/m_• = { 10³, 10⁴, 10⁵, 4 × 10⁵ }.

We find that, regardless of the choice of parameters, a highly eccentric orbit with r₋ = 8R_g in Newtonian dynamics reaches a pericentre of ∼6.45R_g at 1PN order and a pericentre of ∼5.6R_g at 2PN and 2.5PN order. The results of this investigation are shown in Fig. 2.

Thus, we define the loss cone radius as

$$\begin{array}{c}\hfill R_{\mathrm{lc}}=\zeta R_{\mathrm{g}},\end{array}$$(23)

with

$$\begin{array}{c}\hfill \zeta =\{\begin{array}{cc}8\hfill & \text{if Newtonian dynamics}\hfill \\ 6.45\hfill & \text{if 1PN dynamics}\hfill \\ 5.6\hfill & \text{if 2PN or 2.5PN dynamics}\hfill \end{array}.\end{array}$$(24)

2.4. EMRI-to-DP ratio and cliffhanger EMRIs

Through Monte Carlo simulations, Hopman & Alexander (2005) were the first to explore how the probability of EMRI versus DP depends on the initial semi-major axis of the orbit a₀.

The code they presented works in (ℰ, J) space and follows an object on a relativistic orbit, with small perturbations due to GW emission and stochastic two-body scattering. It accounts for GW radiation using (Peters 1964)

$$\begin{array}{c}\hfill {\langle \frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}\rangle }_{\mathrm{GW}}=\frac{32}{5}\frac{G^4{M}_{\bullet }^2{m}_{\bullet }(M_{\bullet }+m_{\bullet })}{c^5{a}^5{(1-e^2)}^{7/2}}(1+\frac{73}{24}{e}^2+\frac{37}{96}{e}^4)\end{array}$$(25)

and

$$\begin{array}{c}\hfill {\langle \frac{\mathrm{d}J}{\mathrm{d}t}\rangle }_{\mathrm{GW}}=-\frac{32}{5}\frac{G^{7/2}{M}_{\bullet }{m}_{\bullet }{(M_{\bullet }+m_{\bullet })}^{3/2}}{c^5{a}^{7/2}{(1-e^2)}^2}(1+\frac{7}{8}{e}^2),\end{array}$$(26)

which are multiplied by P to compute an orbit-averaged estimate of the energy and angular momentum carried away by GWs during a period. Similarly, the code also updates the angular momentum once per orbit to account for two-body relaxation, by employing the orbit-averaged diffusion coefficients D_o[ΔJ] and D_o[(ΔJ)²]. It instead neglects scattering in ℰ, as two-body relaxation is much more efficient in shuffling J rather than ℰ (Binney & Tremaine 2008; Merritt 2013a).

After running many simulations for a central MBH of mass M_• = 3 × 10⁶ M_⊙ surrounded by stellar-mass BHs of mass m_• = 10 M_⊙, Hopman & Alexander (2005) focused on the function

$$\begin{array}{c}\hfill S(a_0)=\frac{N_{\mathrm{EMRI}}(a_0)}{N_{\mathrm{EMRI}}(a_0)+N_{\mathrm{DP}}(a_0)},\end{array}$$(27)

where N_EMRI and N_DP are respectively the counts of EMRIs and DPs formed given an initial a₀, and discovered the existence of a critical value a_c^HA ≃ 5 × 10⁻³R_inf such that S(a₀) = 1 for a₀ ≪ a_c^HA, and S(a₀) = 0 for a₀ ≫ a_c^HA. Raveh & Perets (2021) fixed a sign error in the original implementation of Hopman & Alexander (2005), updating the result to a_c^RP ≃ 2 × 10⁻²R_inf.

Recently, QS24 studied again the problem, considering a more realistic set up in which the stellar-mass BH is scattered by stars of mass m_⋆ = 1 M_⊙ and scattering in energy is not neglected. Most notably, they extended the study for low-mass MBHs, finding that it is possible for EMRIs to form even from large a₀, thus breaking the classical notion of a somewhat sharp transition between a₀ values from which only EMRIs form to those that result in DPs only.

Qunbar & Stone (2024) called EMRIs that form from initially wide orbits cliffhanger EMRIs. The name has to do with the suspenseful ending of the orbital evolution, leading to an EMRI in a region of the parameter space where only DPs were originally expected. In the classical picture, stellar-mass BHs on wide orbits generally end up plunging into the loss cone following a sequence of two-body encounters, which culminate in an almost radial orbit towards the MBH. However, a fraction of these orbits will not enter the capture sphere of the MBH, but graze its surface at a slightly larger distance. In this case, the burst of GW radiation emitted is so strong to easily both reduce a and increase 1 − e by a factor of ten during a single passage (see Fig. 3). If the mass of the MBH is small and the semi-major axis large enough, following the large loss of energy and angular momentum, the stellar-mass BH will end up on a much tighter orbit, closer to or inside the GW-dominated region, with a high probability of later becoming an EMRI.

To better quantify how small the MBH mass should be and to better characterise the phenomenon, we here introduce the ‘cliffhanger region’: a region on the (1 − e, a) plane from which cliffhanger EMRIs come from. Given Eqs. (9) and (20), the semi-major axis a_cliff where t_GW = P as a function of the eccentricity is

$$\begin{array}{c}\hfill a_{\mathrm{cliff}}={\left(\frac{85\pi }{6\sqrt{2}}\right)}^{2/5}{\left(\frac{m_{\bullet }^2(M_{\bullet }+m_{\bullet })}{M_{\bullet }^3}\right)}^{1/5}{(1-e)}^{-7/5}{R}_{\mathrm{g}},\end{array}$$(28)

where we assume e → 1. The loss cone boundary, from Eqs. (8) and (23), is instead identified by the equation r₋ = R_lc, that is

$$\begin{array}{c}\hfill a_{\mathrm{lc}}=\frac{\zeta }{1-e}{R}_{\mathrm{g}}.\end{array}$$(29)

The region in the (1 − e, a) plane identified by the condition

$$\begin{array}{c}\hfill a_{\mathrm{lc}}(e)<a<a_{\mathrm{cliff}}(e)\end{array}$$(30)

is what we refer to as the cliffhanger region (see Fig. 3). As the curves cross at eccentricity

$$\begin{array}{c}\hfill \stackrel{\sim }{e}=1-\frac{85\pi }{6\sqrt{2}}{\zeta }^{-5/2}{\left(\frac{m_{\bullet }^2(M_{\bullet }+m_{\bullet })}{M_{\bullet }^3}\right)}^{1/2},\end{array}$$(31)

the cliffhanger region is defined only for large eccentricities ($e>\stackrel{\sim }{e})$.

Whenever a stellar-mass BH is scattered by two-body encounters onto an orbit which lays inside the cliffhanger region, the orbit will be drastically modified by GW radiation during a single pericentre passage, meaning the variation happens in less than an orbital period. The stellar-mass BH will not plunge onto the MBH after this strong emission. During this process the orbital parameters are related through (Peters 1964)

$$\begin{array}{c}\hfill a(e)=c_0\frac{e^{12/19}}{1-e^2}{(1+\frac{121}{304}{e}^2)}^{870/2299},\end{array}$$(32)

where c₀ is determined by the initial conditions. For e → 1, this implies that a(e)∝(1 − e)⁻¹. Therefore, GWs emission conserves the pericentre⁵ for e → 1, and in the (1 − e, a) plane the orbit will evolve parallel to the loss cone edge.⁶

The scaling of the two curves with the mass of the central MBH is different: Eq. (29) implies that a_lc ∝ M_•, while Eq. (28) shows that a_cliff ∝ M_•^3/5. As we consider larger MBHs, the loss cone edge will cross the t_GW = P curve at larger eccentricities and semi-major axes, eventually to a > R_inf (see Fig. 4). As a result, it gets increasingly harder with larger values of M_• for two-body relaxation to scatter orbits inside the cliffhanger region while avoiding DPs, as the corresponding range in e becomes extremely small. We can therefore predict the maximum MBH mass which allows the formation of cliffhanger EMRIs as in QS24. From Eqs. (9) and (25), follows that the change in energy in a period is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Delta }{\mathcal{E}}_{\mathrm{GW}}(a,e)& =P{\langle \frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}\rangle }_{\mathrm{GW}}\hfill \\ \hfill & \stackrel{e\to 1}{\to }\frac{85\pi }{12\sqrt{2}}\frac{G^{7/2}{M}_{\bullet }^2{m}_{\bullet }{(M_{\bullet }+m_{\bullet })}^{1/2}}{c^5}{a}^{-7/2}{(1-e)}^{-7/2},\hfill \end{array}\end{array}$$(33)

which depends only on the pericentre in the e → 1 limit. Plugging Eq. (28) into Eq. (33) gives

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathcal{E}}_{\mathrm{GW}}^{\mathrm{cliff}}(a_{\mathrm{cliff}},e(a_{\mathrm{cliff}}))=\frac{G{M}_{\bullet }}{2a_{\mathrm{cliff}}}\end{array}$$(34)

which is the energy change in a period due to GWs emission on the t_GW = P curve. On this curve, GW radiation doubles ℰ during a single period,

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{f}}^{\mathrm{cliff}}={\mathcal{E}}_{\mathrm{i}}^{\mathrm{cliff}}+\mathrm{\Delta }{\mathcal{E}}_{\mathrm{GW}}^{\mathrm{cliff}}(a_{\mathrm{cliff},\mathrm{i}},e(a_{\mathrm{cliff},\mathrm{i}}))=\frac{G{M}_{\bullet }}{a_{\mathrm{cliff},\mathrm{i}}}=2{\mathcal{E}}_{\mathrm{i}}^{\mathrm{cliff}},\end{array}$$(35)

meaning that a is halved in the process.

Going deeper into the cliffhanger region allows larger energy jumps. At the limit of the loss cone edge, r₋ = R_lc and Eq. (33) gives

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathcal{E}}_{\mathrm{GW}}^{\mathrm{lc}}(a_{\mathrm{lc}},e(a_{\mathrm{lc}}))=\frac{85\pi }{12\sqrt{2}}{\zeta }^{-7/2}{c}^2{m}_{\bullet }{(M_{\bullet }+m_{\bullet })}^{1/2}{M}_{\bullet }^{-3/2},\end{array}$$(36)

which does not depend on a_lc as we are setting the value of the pericentre. On the edge of the loss cone, the energy jump will be large, thus we can approximate

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{f}}^{\mathrm{lc}}={\mathcal{E}}_{\mathrm{i}}^{\mathrm{lc}}+\mathrm{\Delta }{\mathcal{E}}_{\mathrm{GW}}^{\mathrm{lc}}\simeq \mathrm{\Delta }{\mathcal{E}}_{\mathrm{GW}}^{\mathrm{lc}},\end{array}$$(37)

from which follows that all the cliffhanger jumps on the loss cone edge will end up with the same semi-major axis:

$$\begin{array}{c}\hfill a_{\mathrm{f}}^{\mathrm{lc}}=\frac{6\sqrt{2}}{85\pi }{\zeta }^{7/2}\frac{G}{c^2}\frac{M_{\bullet }^{5/2}}{m_{\bullet }{(M_{\bullet }+m_{\bullet })}^{1/2}}\simeq \frac{6\sqrt{2}}{85\pi }{\zeta }^{7/2}\frac{M_{\bullet }}{m_{\bullet }}{R}_{\mathrm{g}}.\end{array}$$(38)

Finally, we can get a conservative estimate for the maximum MBH mass expected to produce cliffhanger EMRIs. The burst of GWs will end up forming an EMRI as long as the jump in a is inside or close to the EMRI formation region. Using the classical estimate consistent with our results (see Sect. 4.1), a_f^lc ≲ 10⁻²R_inf, we obtain the following:

$$\begin{array}{c}\hfill M_{\bullet }\lesssim \frac{17\pi }{120\sqrt{2}}{\zeta }^{-7/2}\frac{c^2}{{\sigma }^2}{m}_{\bullet }.\end{array}$$(39)

Given our choice of parameters for the M–σ relation (Gültekin et al. 2009), we finally get

$$\begin{array}{c}\hfill M_{\bullet }\lesssim 1.8\times {10}^7{\zeta }^{-2.38}{\left(\frac{m_{\bullet }}{10{\mathrm{M}}_{\odot }}\right)}^{0.68}{\mathrm{M}}_{\odot }.\end{array}$$(40)

It should be noted that this derivation takes into account only the proprieties of the MBH and the stellar-mass BHs, while the distributions of background bodies only marginally enters in the effective critical radius 10⁻² R_inf. In particular, it does not factor in how likely is for two-body relaxation to push the orbit inside the cliffhanger region, nor the fact that a cliffhanger EMRIs can be the result of more than one energy jump, as it is possible for the orbit to enter the cliffhanger region multiple times.

Plugging ζ = 5.6 and m_• = 10 M_⊙ into Eq. (40), we get M_• ≲ 3 × 10⁵ M_⊙, which is indeed roughly the MBH mass for which we start seeing cliffhanger EMRIs forming (see Sect. 4.1). For ζ = 8 and m_• = 10 M_⊙, we find M_• ≲ 1.3 × 10⁵ M_⊙, which is remarkably close to the 10⁵ M_⊙ value obtained by QS24, despite our different assumptions.

3. Numerical simulations of EMRI and DP formation

In this section we illustrate the features of our numerical implementation. We begin with an overview of the integration of the Hamiltonian evolution, then we explain how our treatment of local and orbit-averaged two-body relaxation works. Finally, we detail our procedure to generate the initial conditions for the runs we perform, together with the stopping criteria we use to identify the formation of EMRIs and DPs.

3.1. Numerical overview and setup

In this work we employ an extended version of the code presented in Bonetti et al. (2016), a PN few-body integrator with adaptive step size, which we use to study the evolution of a binary system comprised of an MBH and a stellar-mass BH. The integrator leverages on a C++ implementation of the Bulirsch-Stoer algorithm (Bulirsch & Stoer 1966) coupled with Richardson extrapolation (Richardson 1911) to increase the accuracy of the numerical solution. The adaptive time step allows us to get smaller steps when the dynamics is faster, such as at the pericentre of the orbit.

The code integrates the equation of motion of two-body systems in the Hamiltonian form:

$$\begin{array}{c}\hfill H=H_0+\frac{1}{c^2}{H}_1+\frac{1}{c^4}{H}_2+\frac{1}{c^5}{H}_{2.5}+\mathcal{O}\left(\frac{1}{c^6}\right).\end{array}$$(41)

Here, H₀ is the non-relativistic Newtonian Hamiltonian, while H₁, H₂ and H_2.5 are corrective terms introducing purely relativistic effects. Complete expressions for these terms, matching the form implemented numerically, can be found for example in Bonetti et al. (2016)⁷. H_2.5 is the first non-conservative term and accounts for GW radiation (Jaranowski & Schäfer 1997; Blanchet 2014). We only consider Schwarzschild BHs, thus we ignore the 0.5PN and 1.5PN terms of the expansion, as they account for spin and spin-orbit effects (Schäfer & Jaranowski 2024).

We use the Hamiltonian H to evolve the position x and the momentum p of the stellar-mass BH by integrating its Hamilton’s equations in time t:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle \frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}={\displaystyle {\sum }_n}\frac{1}{c^{2n}}\frac{\partial H_n}{\partial \mathit{p}}}\hfill \\ {\displaystyle \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}=-{\displaystyle {\sum }_n}\frac{1}{c^{2n}}\frac{\partial H_n}{\partial \mathit{x}}}\hfill \end{array}.\end{array}$$(42)

At the same time we forcefully keep the MBH fixed and still at the origin of our Cartesian coordinate system. This choice is motivated by the fact that in a nuclear cluster the MBH will simultaneously interact with virtually all the objects in the stellar distribution, which will only cause a small Brownian motion-like displacement from the centre.

We ignore such a stochastic displacement of the MBH, even if it might be relevant as the number of stellar objects and the mass of the MBH gets lower (Merritt 2013a). Together with the evolution of the binary system, the code can account for the presence of additional external gravitational potentials, possibly representing a distribution of background bodies. We note again that only a single stellar-mass BH is evolved by directly solving Hamilton’s equations.

In this work, we extend the code to include the effects of two-body relaxation on the orbit of the stellar-mass BH. These are due to the gravitational interaction of the stellar-mass BH with the rest of the stellar distribution. These stochastic interactions are fully characterised by the diffusion coefficients of the velocity component computed from their distribution function (Shapiro & Marchant 1978; Cohn & Kulsrud 1978; Merritt 2013a). We include the effects of two-body relaxation via a Monte Carlo approach by repeating the following procedure at the end of each time step of the integration:

1. we compute the mean and the standard deviation of the expected variation ⟨Δv⟩ of the velocity v of the stellar-mass BH, considering the two-body encounters with the scattering objects which happened during the last time step;

2. we randomly draw a velocity change Δv from the appropriate Gaussian distribution, which we build according to the previous step;

3. before moving to the next time step, we update the momentum of the stellar-mass BH to p′=p + m_•Δv, where m_• is the mass of the stellar-mass BH⁸.

Following this procedure, we can statistically account for two-body encounters locally and multiple times within a single orbit, without leveraging on orbit-by-orbit approaches (QS24) or orbit-averages of the diffusion coefficients (Zhang & Amaro Seoane 2024).

For our simulations, we consider a stellar-mass BH of mass m_• = 10 M_⊙ and explore five possible masses for the MBH: M_• = { 10⁴, 10⁵, 3 × 10⁵, 10⁶, 4 × 10⁶ } M_⊙. We also group the background objects into two populations:

1. stars with m_⋆ = 1 M_⊙, M_tot, ⋆ = 20M_•, r_a, ⋆ = 4R_inf and γ_⋆ = 1.5;

2. stellar-mass BHs with m_• = 10 M_⊙, M_tot, • = 0.2M_•, r_a, • = 4R_inf and γ_• = 1.7.

This choice corresponds to the two-components Bahcall-Wolf solution commonly used in literature (Bahcall & Wolf 1977, γ ≃ 1.5 for the lightest component, γ ≃ 1.75 for the heaviest component), where we slightly reduced the slope of the heaviest component to better represent their distribution in the inner regions of the stellar cluster, where they dominate relaxation (see Fig. 4 of Broggi et al. 2022, for the quasi-relaxed distribution of the stellar-mass BHs). The scale radius of the distribution r_a is set so that the stellar cluster in which our 4 × 10⁶ M_⊙ MBH is immersed is consistent with observational estimates of the density at the influence radius of SgrA^* (e.g. Schödel et al. 2007).

3.2. Local treatment of two-body relaxation

As explained in Appendix A, in order to compute local diffusion coefficients we have to first find f_α(ℰ) given n_α and then integrate it over ℰ. This process involves nested integrals and derivatives, making it impractical to numerically perform such computations on the fly at each time step of our simulations. Thus, we precompute the auxiliary functions ℐ_k, α and 𝒥_k, α, that are directly related to the diffusion coefficients, for some logarithmically spaced values of ℰ and ψ, and then, during the execution, we bilinearly interpolate these values of ℐ_k, α and 𝒥_k, α. In this way, we get approximate local values for those coefficients given the exact position and velocity of the stellar-mass BH at a particular time step.

The diffusion coefficients in velocity are then directly related to the parameters of the Gaussian distributions characterising Δv_∥ and a Δv_⊥ at each time step of the integration:

$$\begin{array}{c}\hfill \mathrm{\Delta }{v}_{\Vert }={\mu }_{\Vert }+s_{\Vert }{\sigma }_{\Vert },\mathrm{\Delta }{v}_{\perp }={\mu }_{\perp }+s_{\perp }{\sigma }_{\perp }.\end{array}$$(43)

Here, s_∥ and s_⊥ are two uncorrelated random values, each extracted from a Normal distribution, while

$$\begin{array}{c}\hfill \begin{array}{cccc}{\mu }_{\Vert }\hfill & =\mathrm{D}\left[\mathrm{\Delta }{v}_{\Vert }\right]\mathrm{\Delta }t,\hfill & \hfill {\sigma }_{\Vert }^2& =\mathrm{D}\left[{\left(\mathrm{\Delta }{v}_{\Vert }\right)}^2\right]\mathrm{\Delta }t-{\left(\mathrm{D}\left[\mathrm{\Delta }{v}_{\Vert }\right]\mathrm{\Delta }t\right)}^2,\\ {\mu }_{\perp }\hfill & =0,\hfill & \hfill {\sigma }_{\perp }^2& =\mathrm{D}\left[{\left(\mathrm{\Delta }{v}_{\perp }\right)}^2\right]\mathrm{\Delta }t.\end{array}\end{array}$$(44)

In these equations, Δt is the length of the particular time step we are considering. Thus, the values of Δv_∥ and Δv_⊥ we draw represent a statistical realisation of the change in velocity due to all of the encounters that have occurred during that time step.

Once Δv_⊥ is drawn, we set its direction randomly, as it is distributed uniformly in the plane perpendicular to the vector v. Finally, we vectorially sum together Δv_∥ and Δv_⊥ to get the total variation Δv of the velocity of the stellar-mass BH, and update its momentum to p′=p + m_•Δv before the next step of the integration of Hamilton’s equations begins. Figure 5 shows how this procedure affects a few consecutive apocentre passages. Deviations from regular elliptical orbits are evident, albeit small. Local velocity kicks are more frequent near pericentre due to the adaptive nature of the time step used by the main integrator.

3.3. Orbit-averaged treatment of two-body relaxation

The usual approach to implement two-body relaxation into the orbit evolution differs from what we have just described. Codes usually work in the (ℰ, J) space and thus employ orbit-averaged diffusion coefficients of ℰ and J, which can be computed from those of velocity which we present in Appendix A (see Binney & Tremaine 2008; Merritt 2013a). In this work, we are interested in testing the reliability of the orbit-averaging approximation in opposition to our local approach. To make a meaningful comparison between the two methods, we present some simulations using the orbit-averaged approximation in Sect. 4.3.

Including orbit-averaging in our code is not straightforward, as we do not work with ℰ and J. Here we summarise our implementation:

1. at each time step of the integration, we draw δv_∥ and δv_⊥⁹ as explained in Sect. 3.2;

2. from δv_∥ and δv_⊥, we compute the corresponding change in energy δℰ and angular momentum δJ;

3. we move to the next time step without applying the kick and let the integration proceed, keeping track of the cumulative change in energy Δℰ = ∑_iδℰ_i and angular momentum ΔJ = ∑_iδJ_i;

4. once per orbit, from Δℰ and ΔJ we compute a compatible kick Δv and update p to p′ = p + m_•Δv;

5. once the velocity of the stellar-mass BH is changed, we reset Δℰ and ΔJ to zero, and repeat steps 1–4.

Following this procedure, we include the effects of stochastic kicks on the orbit once per period. The changes in ℰ and J are naturally orbit-integrated as the result of this procedure: the δℰ_i and δJ_i terms we sum together are indeed intrinsically weighted by the length of the time steps δt_i, which contribute through the δv_∥, i and δv_⊥, i terms.

In Appendix B we show that, once δv_∥ and δv_⊥ are drawn, then

$$\begin{array}{cc}\hfill \delta \mathcal{E}& =-v\delta v_{\Vert }-\frac{1}{2}{\left(\delta v_{\Vert }\right)}^2-\frac{1}{2}{\left(\delta v_{\perp }\right)}^2,\hfill \end{array}$$(45)

$$\begin{array}{cc}\hfill \frac{\delta J}{J}& =\begin{array}{cc}& \frac{\delta v_{\Vert }}{v}+\delta v_{\perp }\frac{v_{\mathrm{r}}}{v{v}_{\mathrm{t}}}sin\theta +\frac{1}{2}{\left(\frac{\delta v_{\perp }cos\theta }{v_{\mathrm{t}}}\right)}^2+\mathcal{O}\left({(\delta v)}^3\right).\hfill \end{array}\hfill \end{array}$$(46)

Here v_t = J/r is the tangential velocity, $v_{\mathrm{r}}=\sqrt{v^2-v_{\mathrm{t}}^2}$ is the radial velocity, and θ is an angle uniformly drawn from the interval [0, 2π). We compute J as

$$\begin{array}{c}\hfill J=|\mathit{J}|=\frac{\mathit{x}\times \mathit{p}}{m_{\bullet }}.\end{array}$$(47)

In Appendix C we show that a suitable orbit-integrated kick in velocity can be computed from ΔE and ΔJ as

$$\begin{array}{cc}\hfill \mathrm{\Delta }{v}_{\mathrm{r}}& =\sqrt{v^2-2\mathrm{\Delta }\mathcal{E}-{\left(\frac{J+\mathrm{\Delta }J}{r}\right)}^2}-v_{\mathrm{r}},\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{v}_{\mathrm{t}}& =\frac{(J+\mathrm{\Delta }J)cos\beta }{r}-v_{\mathrm{t}},\hfill \end{array}$$(49)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{v}_3& =\frac{(J+\mathrm{\Delta }J)sin\beta }{r}.\hfill \end{array}$$(50)

Here the subscript 3 indicates the direction perpendicular to the orbital plane (thus, v₃ = 0 before each kick), while β is an angle uniformly drawn from the interval [0, 2π).

When updating the orbit, we change the velocity of the stellar-mass BH without displacing it. This is possible only if the old and the new orbit intersect at the location where Δv is applied, a fact that we ensure by applying the kick when r ≃ a (Fig. 6 illustrates the situation with a sketch). It is safe to assume that the orbits intersect there, as long as the energy of the particle (and thus a, since ℰ ≃ GM_•/2a) does not change drastically along a single period¹⁰. In case this assumption is not verified, the orbit-averaging approximation is not valid any more, as the time needed for two-body relaxation to significantly change ℰ is less than the orbital period (Binney & Tremaine 2008; Merritt 2013a).

This motivates us to apply the kick Δv to the stellar-mass BH once per orbit at the time step when r becomes smaller than a, coming from the apocentre. This choice reduces the number of cases where Eq. (48) has no real solutions to a negligible fraction, where orbit-averaging would fail because of the significant change received in one single orbit. In this case, the code simply halts, and we remove the run from the analysis.

The code computes a and e as

$$\begin{array}{c}\hfill a=\frac{r_{+}+r_{-}}{2},e=\frac{r_{+}-r_{-}}{r_{+}+r_{-}},\end{array}$$(51)

each time the pericentre or apocentre of an orbit is reached, where these are identified as the distance from the centre when the radial velocity flips sign. If a and e are computed at apocentre, then r₋ is taken as the last pericentre passed. Instead, if they are computed at pericentre, then r₊ is taken as the last apocentre passed. Figure 7 shows how this procedure affects a few consecutive orbits. In the figure, orbits are almost perfect ellipses, with small deviations due to PN effects.

3.4. Stopping criteria and initial conditions

The main aim of the code is that of identifying the formation of EMRIs and DPs. To this scope, we employ the definitions presented in Sect. 2.2, with some corrections needed given the computational cost of our hybrid few-body/Monte Carlo code.

As our aim was to gather a statistically significant number of EMRIs and DPs, we chose to limit the parameter space that we explored to a sub-portion close to the loss cone edge. Moreover, as we assume a steady-state stellar profile, we are interested in exploring orbital parameters for which the relaxation time t_rlx is not too large, as the initial conditions and the diffusion coefficients would be affected by the overall relaxation of the system. In particular, we consider as our timescale of investigation ∼1% of the time-to-peak of (classical) EMRI rates in the systems we are considering, which estimates the timescale of evolution of the distribution function of scatterers (Broggi et al. 2022)

$$\begin{array}{c}\hfill t_{\max }=10{\left(\frac{M_{\bullet }}{4\times {10}^6{\mathrm{M}}_{\odot }}\right)}^{1.29}\mathrm{Myr},\end{array}$$(52)

and exclude from our investigations the region in the (1 − e, a) phase space where t_rlx > t_max. The scaling with M_• for t_max comes from the fact that (Merritt 2013a)

$$\begin{array}{c}\hfill t_{\mathrm{rlx}}\propto M_{\bullet }^2{\sigma }^{-2}\propto M_{\bullet }^{1.29},\end{array}$$(53)

where we leverage on the M–σ relation (Gültekin et al. 2009) to obtain the final scaling with the MBH mass. We note that t_max is not to be intended as a dynamical stopping condition (meaning the simulations do not halt at the time t = t_max), rather the run is ended only when the orbital elements cross over the t_rlx = t_max curve on the (1 − e, a) plane.

We initialise all our simulations inside the region of interest where t_rlx < t_max following a procedure described at the end of this section and we stop and discard any run where the orbit leaves the region. This choice does not affect the EMRI-to-DP ratio. In fact, t_max is so long with respect to the orbital period that any plunge or EMRI must enter this region before the capture. As long as the procedure is not biased towards EMRI or DP candidates, once a particle is discarded we are safe to reinitialise it as an uncorrelated particle.

The reinitialisation procedure becomes biased when the restriction to t_rlx < t_max excludes a significant part of the t_GW = 10⁻³t_rlx curve. In this case, it can happen that some orbits well inside the GW emission dominated region but still far from the loss cone, can reach the t_rlx = t_max boundary before crossing the t_GW = 10⁻³t_rlx separatrix, thus being erroneously stopped and reinitialised instead of being counted as EMRIs. A clear example of such an orbit is shown by the blue line in Fig. 8: the event will clearly produce an EMRI, but according to our reinitialisation criterion is effectively discarded and resimulated.

As DP candidates are not affected by this issue, the EMRI-to-DP ratio is influenced by these discarded EMRIs. This problem mostly affects the EMRI-to-DP ratio at small initial semi-major axes for large M_•, as we show in Sect. 4.1.

For this reason, we correct our counts of EMRIs to also include orbits which cross into the excluded region while inside the GW-dominated region. Before accepting the orbit as an EMRI, we also verify that the final value of a is below the evolution line of an orbit which evolves as Eq. (32), where c₀ is chosen so that the line crosses the minimum of the t_GW = 10⁻³t_rlx curve (again, see Fig. 8). This means that, ignoring stochastic kicks due to two-body encounters, the orbit would eventually reach our condition for EMRI formation, while inside the excluded region. We also point out that near the t_rlx = t_max curve, the condition t_GW < t_rlx is likely sufficient to identify EMRIs anyway, as the stricter condition t_GW < 10⁻³ t_rlx is mainly introduced to avoid misclassifying DPs as EMRIs, and such a misclassification can only occur for 1 − e ≲ 10⁻⁴, where the t_GW = t_rlx curve closely approaches the loss cone.

To summarise, the possible outcomes of a simulation are:

1. a DP, if the separation r between MBH and stellar-mass BH satisfies the condition r < R_lc;

2. an EMRI, either if (i) the condition t_GW = 10⁻³t_rlx is satisfied, or (ii) the curve t_rlx = t_max is reached with a final a such that the curve t_GW = 10⁻³t_rlx is eventually reached at a later time assuming only GW emission;

3. a reinitialisation, in all other cases in which the condition t_rlx = t_max is reached.

Condition (ii) for EMRI formation is the correction introduced above, and its importance will be shown in Sect. 4.1. Examples of each of these outcomes are shown in Figs. 1, 3, and 8.

Finally, the code halts if Eq. (48) has no real solutions when looking for an orbit-averaged kick, or if ℰ < 0 at some point during the simulation (meaning the stellar-mass BH has become unbound). In practice, the former situation excluded a negligible number of simulations, and the latter never occurred.

For each nuclear cluster configuration, we select a set of possible initial semi-major axes a₀ uniformly distributed in logarithm. Then, we adopt the steps described below to randomly draw a different realistic value for the initial eccentricity e₀ of each simulation. Given a₀, we compute

$$\begin{array}{c}\hfill {\mathcal{R}}_{\max }={\mathcal{R}|}_{t_{\mathrm{rlx}}=t_{\max }}\end{array}$$(54)

and

$$\begin{array}{c}\hfill {\mathcal{R}}_{\min }=min\{{\mathcal{R}}_{\mathrm{lc}},{\mathcal{R}|}_{t_{\mathrm{GW}}={10}^{-3}{t}_{\mathrm{rlx}}}\}\end{array}$$(55)

at the corresponding energy ℰ₀ = GM_•/2a₀. Then, we randomly draw an initial value for ℛ₀ from the probability distribution function

$$\begin{array}{c}\hfill p(\mathcal{R})=\frac{log(\mathcal{R}/{\mathcal{R}}_{\mathrm{q}})}{{\mathcal{R}}_{\max }(ln({\mathcal{R}}_{\max }/{\mathcal{R}}_{\mathrm{q}})-1)-{\mathcal{R}}_{\min }(ln({\mathcal{R}}_{\min }/{\mathcal{R}}_{\mathrm{q}})-1)}\end{array}$$(56)

to reproduce the Cohn-Kulsrud distribution (Cohn & Kulsrud 1978, also see Sect. 2.1), where the normalisation factor is such that ${\displaystyle {\int }_{{\mathcal{R}}_{\min }}^{{\mathcal{R}}_{\max }}}\mathrm{d}\mathcal{R}p(\mathcal{R})=1$. Finally, from ℛ₀, we compute

$$\begin{array}{c}\hfill e_0=\sqrt{1-{\mathcal{R}}_0}\end{array}$$(57)

and we initialise the orbit with semi-major axis a₀ and eccentricity e₀ with Newtonian initial conditions, placing the stellar-mass BH at the apocentre.

In our procedure, we initialise all particles at the apocentre of their initial orbits, rather than assigning a random location along the orbit. We select this point because PN corrections are least significant there and this choice does not bias the likelihood of a particle becoming a DP or an EMRI. Before reaching a disruptive pericentre, the particle will either quickly modify its orbital parameters through two-body scattering (when t_rlx ≲ P), effectively erasing any memory of the initial condition, or complete more than one period on the same orbit (when t_rlx ≳ P), rendering the specific initial location along the orbit irrelevant.

4. Results

In this section we present the results we obtained running a large number of simulations with the code described in Sect. 3.

4.1. EMRI-to-DP ratio and cliffhanger EMRIs

For each simulated central MBH mass, we build the S(a₀) function given the counts of EMRIs and DPs. Given the success/failure nature of the problem (we can see EMRI formation as a success, and DP formation as a failure), we employ the binomial statistic to characterise the uncertainties on our findings. Specifically, we compute the corresponding binomial proportion confidence interval at 1σ confidence level using the Wilson score interval (Wilson 1927), an improved version of the commonly used standard Wald interval. The Wilson score interval avoids issues when S(a₀) approaches 0 or 1, such as falsely implying certainty or error bars overshooting outside the [0, 1] range (Brown et al. 2001).

We report our main findings in the left panel of Fig. 9. For M_• = 10⁶ M_⊙ and M_• = 4 × 10⁶ M_⊙, we recover the expected dichotomy between EMRIs and DPs. In this case, EMRIs can only form from orbits with a₀ ≲ a_c, while those with a₀ ≳ a_c exclusively result in DPs. We find that the transition between the two regimes occurs over about an order of magnitude in a₀ at a_c ≃ 10⁻²R_inf, which is consistent within a factor of two with previous works (Hopman & Alexander 2005; Raveh & Perets 2021, QS24). We believe that the quantitative difference is driven by differences in the nuclear cluster models, in the accuracy of the dynamical model (e.g. PN order) and in the treatment of two-body relaxation. We better explore the impact that some of these factors have on the shape of S(a₀) in Sects. 4.2 and 4.3.

Conversely, for M_• < 10⁶ M_⊙, we observe a relevant number of EMRIs formed from initially wide orbits, with a₀ ≫ a_c. These events are consistent with the cliffhanger EMRI population first seen in Monte Carlo simulations by QS24 and, to the best of our knowledge, never found in previous work.

We observe cliffhanger EMRIs forming for M_• ≲ 3 × 10⁵ M_⊙ and they appear to be more common as M_• decreases, consistently with the argument based on the (1 − e, a) phase space presented in Sect. 2.4.

Our results for S(a₀) are compared to those obtained by QS24 and other authors in the right panel of Fig. 9. First, we can notice the large scatter in a_c/R_inf that spans almost an order of magnitude, from ≈3 × 10⁻³ (QS24) to ≈2 × 10⁻² (Raveh & Perets 2021). This is likely due to the differences in the parameters of the simulated systems as well as in the treatment of the dynamics and loss cone definition, as we discuss in Sect. 4.3. Focusing on cliffhanger EMRIs, QS24 observe their formation only for MBHs with M_• ≲ 5 × 10⁴ M_⊙, a threshold set almost an order of magnitude below our limiting value. Moreover, we see a much larger fraction of EMRIs over DPs at large a₀, with a maximum ratio of S(a₀)≃0.65 at a₀ ≫ 10⁻²R_inf for M_• = 10⁴ M_⊙. As a further difference, our S(a₀) functions hint at a downward trend for a₀ ≫ 10⁻²R_inf, as opposed to tending to a constant at a₀ → ∞. This could be due to the fact that QS24 stopped their investigation at a_max ≃ R_inf/3, as their implementation does not account for the stellar potential.

Consider now a_α and a_β as the values of the semi-major axis of the orbit at two subsequent generic apocentres. For each run that resulted in the formation of an EMRI, we consider the largest ratio a_α/a_β to set a criterion to distinguish cliffhanger and traditional EMRIs. Consistently with the reasoning that led us to Eq. (35), we determine the passage across the cliffhanger region if at some point in the orbital evolution the semi-major axis is cut in half or more during a single period, thus for all the runs where (a_α/a_β)_max^EMRI > 2. This threshold is somehow arbitrary, but offers a quantitative way to identify cliffhanger EMRIs. Figure 10 shows the result of this investigation for M_• = 10⁵ M_⊙ and M_• = 10⁶ M_⊙. For both MBH masses, all EMRIs forming for a₀ ≪ 10⁻²R_inf are classical EMRIs. Both cases show an increase in the average value of (a_α/a_β)_max^EMRI for larger a₀, but only for M_• = 10⁵ M_⊙ the threshold of 2 is crossed. Moreover, we verify that the prediction for the largest possible (a_α/a_β)_max^EMRI ratio (which occurs at the loss cone edge) given in Eq. (38) is verified for M_• = 10⁵ M_⊙ for large a₀, with some deviations for smaller a₀. This is likely because the approximation given in Eq. (37) is less precise for smaller semi-major axes.

Finally, in Fig. 11, we show the impact on the function S(a₀) of the correction to the stopping criterion for EMRIs detailed in Sect. 3.4. Largest masses are more affected by the issue, in particular for low a₀ and near the inflection point at a ≃ 10⁻²R_inf. The figure clearly shows that this correction is needed to correctly recover the limit to 1 of S(a₀) at a₀ ≪ 10⁻²R_inf, while it does not affect the curve at a₀ ≫ R_inf region, where we observe cliffhanger EMRIs.

4.2. Local versus orbit-averaged approach

One of the main objectives of our investigation is that of testing the reliability of the orbit-averaged approximation to describe the complex interactions in nuclear clusters. As explained in Sect. 2.1, this approach is based on the idea of orbit-averaging the local diffusion coefficients, which in reality strongly depend on the exact position r and velocity v of the orbiting body in consideration. In highly eccentric orbits, such as the one needed to form EMRIs, both r and v oscillate by orders of magnitude in a single period, and orbit-averaging diffusion coefficients might lead to a mismodelling of the very different features displayed by relaxation at r₋ or r₊. Another common practice in the orbit-averaging approach that is likely to lead to inaccuracies is that of accounting for two-body relaxation only once per period to update the orbital elements. This effectively means that once the energy and angular momentum of the particle are updated (usually at apocentre) the orbit is frozen for a full period. This treatment is, however, at odd with the concept of full loss cone itself, as we come across the contradiction of being in the regime where escaping from the loss cone should be easy while being effectively locked onto an orbit, that thus cannot escape.

The relaxation timescale for a particle with energy ℰ and pericentre r₋ estimates the average time that the particle passes on the orbit before being scattered away. When r₋ is much smaller than the semi-major axis of the orbit, it estimates the time before r₋ is changed by the order of itself, as energy remains almost constant (Cohn & Kulsrud 1978). In this small r₋ regime the relaxation timescale for the orbit can be estimated as t_rlx ≃ P r₋/(q r_lc) (Merritt 2013a; Bortolas et al. 2023; Broggi et al. 2024), where P is the orbital period. Consider now an orbit with energy ℰ and r₋ < r_lc: in a time P it will be reached by a number dN of stellar objects due to two-body relaxation. In the orbit-averaged approach, they will all be DPs produced in the subsequent orbital period, with a rate from that orbit of $\mathrm{d}\dot{N}=\mathrm{d}N/P$. In the local approach, if t_rlx < P, only the stellar objects that will take less than t_rlx to reach the pericentre will be DPs, and the corresponding rate will be $\mathrm{d}\dot{N}\simeq \mathrm{d}N{t}_{\mathrm{rlx}}/P^2=r_{-}/(q{r}_{\mathrm{lc}})\mathrm{d}N/P$. Therefore, the orbit-averaged approach will overestimate the rate of DPs for q ≳ 1 due to the ‘frozen orbit’ approach, as they have r₋ < r_lc by definition. The same reasoning applies to cliffhangers, that are produced at very small pericentres, too. However, being their pericentre larger than r_lc by definition, the overestimation of cliffhangers in the orbit-averaged assumption will be lower. Consequently, we expect that the orbit-averaged approach will underestimate S(a).

In order to test this expectation we re-ran the three sets of simulations with M_• = 10⁴, 10⁵, 3 × 10⁵ M_⊙ with the exact same system parameters but employing the orbit-averaged procedure described in Sec. 3.3. Results of these test runs are compared to our default simulations in Fig. 12, and appear to confirm our expectation. The S(a₀) functions computed in the local vs orbit-averaged simulations nicely match up to the value of a₀/R_inf for which q = 1, which marks the transition from empty to full loss cone (vertical dotted lines in the figure). At larger a₀, the orbit-averaged approach tends to produce less cliffhanger EMRIs and more plunges, thus resulting in a lower S(a₀). We note that the trend becomes more severe by lowering the central MBH mass, consistent with the fact that the cliffhanger region becomes larger and larger. In practice, the smaller the central MBH, the bigger gets the cliffhanger region compared to the loss cone, and it becomes increasingly more difficult for the orbit-averaged approximation to capture the detailed dynamics of particle scattered within the cliffhanger region in real time along the orbit close to the pericentre, before getting to the loss cone, thus leading to an over-estimation of DPs.

4.3. Detailed investigation of the differences with QS24

Our main results display some deviations with respect to other works in the literature. In particular the emergence of cliffhanger EMRIs and their abundance does not match findings by QS24. This is not surprising, given the many differences between our simulations and theirs. In particular, we simulate a two-population system, including the effect of solar-mass stars and stellar-mass BHs in the computation of the diffusion coefficient, we adopt a local treatment of two-body relaxation and we evolve the stellar-mass BH using PN equations of motion. Conversely, in QS24 scattering comes from a single population of solar-mass object, is treated as orbit-averaged and the dynamics of the stellar-mass BH is Newtonian.

In the previous section, we already highlighted the importance of the local treatment of relaxation, which leads to a much larger relative fraction of EMRIs in the full loss cone regime. This can certainly account for some of the differences with QS24. For example in the case of M_• = 10⁴ M_⊙, the orbit-averaged tests shown in Fig. 12 result in S(a₀)≈0.4 at large a₀, which is much closer to the value of ≈0.25 found in QS24 compared to the ≳0.6 found in our default simulations employing a local treatment of relaxation.

In order to further test how the differences between the QS24 and our studies affect the results, we perform further targeted simulations, focusing on the M_• = 3 × 10⁵ M_⊙ case, as it is close to the MBH mass value at which the majority of LISA EMRI detections is expected (Babak et al. 2017). We try to reproduce the results by QS24¹¹ by running our simulations with a setup which is as close as possible to theirs. We use a single population of scattering objects of m_⋆ = 1 M_⊙, distributed according to a Denhen density profile characterised by M_tot, ⋆ = 1425M_•, r_a, ⋆ = 223.37R_inf and γ_⋆ = 1.75. With these parameters, our density profile departs from the single r^−1.75 power law simulated by QS24 by less than 1% within the MBH influence radius.

Qunbar & Stone (2024) use Newtonian dynamics, and take into account emission of GWs through the equations detailed in Peters (1964). We do the same by setting H₁ = H₂ = H_2.5 = 0 and adding the contributions

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathcal{E}}_{\mathrm{GW}}=P{\langle \frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}\rangle }_{\mathrm{GW}},\mathrm{\Delta }{J}_{\mathrm{GW}}=P{\langle \frac{\mathrm{d}J}{\mathrm{d}t}\rangle }_{\mathrm{GW}}\end{array}$$(58)

to the orbital sums of the stochastic kicks Δℰ = ∑_iδℰ_i and ΔJ = ∑_iδJ_i once per orbit, before applying the velocity kick following the orbit-averaged procedure described in Sect. 3.3.

We report the results of these tests in Fig. 13. Despite our attempt to reproduce the system employed in QS24 as close as possible, we cannot replicate their results. Although the shape of the S(a₀) function is very similar, our a_c is approximately a factor of two smaller than what they found. It is likely that this discrepancy stems from the neglect of the stellar potential in QS24, or in differences in the specific implementation of orbit-averaging. Further investigation is still needed to get a definitive answer. Moreover, we observe that S(a₀) is not modified by using the 2.5PN term to account for GW radiation instead of the equations detailed in Peters (1964). As discussed in the previous section, for this central MBH mass, the local treatment of two-body relaxation does not have a strong impact on the S(a₀) function. We find that in Newtonian dynamics the local treatment results in no discernible deviation of S(a₀) compared to the orbit-averaged approximation. At second PN order instead, we observe some more significant deviations. In particular, as noted in the previous section, our local approach predicts only slightly more DPs for a₀ ≃ 3 × 10⁻³R_inf and more EMRIs for a₀ ≳ 2 × 10⁻²R_inf. Finally, a large difference comes from the choice of dynamics. When going from Newtonian to PN dynamics, S(a₀) is shifted to larger semi-major axes and cliffhanger EMRIs become more likely. This is because PN orbits can get closer to the MBH without being accreted as DPs (R_lc = 5.6R_g instead of 8R_g), thus having the chance of experiencing much larger jumps in the (a, 1 − e) plane along the loss cone line. This, in turn, means that cliffhanger EMRIs can form from larger a₀, consistent with Fig. 13.

These results show that a_c (i.e. the EMRI-DP separatrix) as well as the abundance of cliffhanger EMRIs depends on the treatment of the stellar-mass BH orbit dynamics. Since PN orbits are a better approximation to GR than Newtonian gravity, we speculate that our PN results are closer to reality than those obtained by employing Newtonian dynamics. However, a full geodesic orbital description in GR would be needed to provide an accurate characterisation of the S(a₀) function.

4.4. Effect on EMRI rates

We now apply our results on S(a₀) to assess the impact of cliffhanger EMRIs on the rates of EMRI predicted by Fokker-Planck codes. Classically, the rate of EMRIs is computed from the differential rate in energy of loss cone captures

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{lc}}(\mathcal{E})\equiv \frac{\mathrm{d}\dot{N}}{\mathrm{d}\mathcal{E}}\end{array}$$(59)

that can be directly computed during the evolution (Cohn & Kulsrud 1978). The total rate of EMRIs and DPs per unit time can be computed as follows:

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{tot}}={\displaystyle {\int }_0^{\infty }}\mathrm{d}\mathcal{E}{\mathcal{F}}_{\mathrm{lc}}(\mathcal{E}).\end{array}$$(60)

In general, S(a₀) acts on the differential rate as a transfer function, giving the fraction of loss cone events happening at semimajor axis a₀ that will result in EMRIs. The function 1 − S(a₀) gives, symmetrically, the fraction of DPs. In order to apply the transfer function S(a₀), one needs to map the orbit in a generic potential to a typical scale radius, interpreted as the semi-major axis when the stellar object enters the loss cone or the GW-dominated region, meaning it is considered to be captured. A reliable choice for objects captured at a given energy is the radius r_c of the circular orbit with the same energy, that reduces to the semi-major axis in the Keplerian regime:

$$\begin{array}{c}\hfill S(\mathcal{E})=S(r_{\mathrm{c}}(\mathcal{E})).\end{array}$$(61)

Another possibility is a₀ = (r₋ + r₊)/2, but it is affected by the stellar potential at smaller ℰ. The two mappings for typical loss cone orbits depend on the total potential evaluated at r_c = a₀ and r₊ ≃ 2 a₀ respectively. The rates of DPs and EMRIs can be computed as

$$\begin{array}{cc}\hfill {\dot{N}}_{\mathrm{DP}}& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}\mathcal{E}[1-S(\mathcal{E})]{\mathcal{F}}_{\mathrm{lc}}(\mathcal{E}),\hfill \end{array}$$(62)

$$\begin{array}{cc}\hfill {\dot{N}}_{\mathrm{EMRI}}& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}\mathcal{E}S(\mathcal{E}){\mathcal{F}}_{\mathrm{lc}}(\mathcal{E}).\hfill \end{array}$$(63)

When using the classical discriminant between EMRIs and DPs, one sets S(ℰ) = Θ(ℰ − ℰ_c), where Θ(x) is the Heaviside step function, that is zero when its argument is negative and one otherwise, and ℰ_c is the energy such that r_c(ℰ_c) = a_c. Accordingly, the classical rate of the two classes of captures, which we denote with the upper index ‘cl’, is

$$\begin{array}{cc}\hfill {\dot{N}}_{\mathrm{DP}}^{\mathrm{cl}}& ={\displaystyle {\int }_0^{{\mathcal{E}}_{\mathrm{c}}}}\mathrm{d}\mathcal{E}{\mathcal{F}}_{\mathrm{lc}}(\mathcal{E}),\hfill \end{array}$$(64)

$$\begin{array}{cc}\hfill {\dot{N}}_{\mathrm{EMRI}}^{\mathrm{cl}}& ={\displaystyle {\int }_{{\mathcal{E}}_{\mathrm{c}}}^{\infty }}\mathrm{d}\mathcal{E}{\mathcal{F}}_{\mathrm{lc}}(\mathcal{E}).\hfill \end{array}$$(65)

4.4.1. The case M_• = 3 × [sans]10⁵ M_⊙: implications for the most relevant LISA systems

We first specialise to the case M_• = 3 × 10⁵ M_⊙, which is where the bulk of LISA detections are expected, and use the code presented in Broggi et al. (2022) to simulate the relaxation of the nuclear cluster surrounding it. The initial Dehnen profile of stars and stellar-mass BHs has γ_⋆ = γ_• = 1.5, while the other properties match those of Sect. 3.1. At later times, the system relaxes and two-body scattering in the inner region are dominated by stellar-mass BHs, whose distribution reaches γ_• = 1.7, as in Sect. 3.1. The code accounts for a non-evolving stellar potential set by the initial conditions and includes the effects of two-body relaxation, but does not account for GWs emission. Since the GWs driven region of phase space encloses the loss cone for a < a_c, we interpret ℱ_lc for ℰ > ℰ_c as the differential rate entering the GW dominated region¹². For this system R_inf = 0.57 pc and ℰ_c = 63 σ_inf², where σ_inf = σ(R_inf). We build S(ℰ) by interpolating the function shown in Fig. 9 linearly in the log-log space. We extrapolate as S(ℰ) = 0 for ℰ → 0 (a₀ → ∞), and S(ℰ) = 1 for ℰ → ∞ (a₀ → 0).

In Fig. 14, we show the differential rate ℱ_lc computed at t = 25 Myr and t = 45 Myr, both after the initial numerical transient (coming from the perfectly isotropic initial conditions), but before the peak of the classical EMRI rates, which happens at t = 75 Myr. We also show the differential rate of EMRIs

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{EMRI}}(\mathcal{E})\equiv {\mathcal{F}}_{\mathrm{lc}}(\mathcal{E})S(\mathcal{E})\end{array}$$(66)

for a direct comparison with its classical estimate

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{EMRI}}^{\mathrm{cl}}(\mathcal{E})\equiv {\mathcal{F}}_{\mathrm{lc}}(\mathcal{E})\mathrm{\Theta }(\mathcal{E}-{\mathcal{E}}_{\mathrm{c}}).\end{array}$$(67)

The instantaneous rate of EMRIs and DPs according to our and the classical prescriptions are reported in Table 1. At t = 25 Myr, when EMRIs classically account for 15% of the total loss cone event rate, Eq. (62) gives a total rate of EMRIs larger by ∼35%. The rate of DPs is only reduced by ∼5%, as they are overall more numerous. Therefore, with our new estimates EMRIs account for 20% of the total loss cone event rate. At t = 45 Myr, EMRIs classically account for ∼30% of the total loss cone event rate, with our new estimate of the EMRI and DP rates comparable to the classical estimates.

Our model predicts a fraction of EMRIs produced in the classical region ℰ > ℰ_c of 45% at t = 25 Myr, and 60% at t = 45 Myr. Even with a small value S(a₀)≃0.1 at large semi-major axes, cliffhanger EMRIs are a relevant contribution to the overall count of EMRIs due to the shape of the loss cone flux.

Aside from the rates, the inclusion of S(a₀) has three direct consequences on the energy distribution of EMRIs and DPs:

• about 10% of DPs are produced for a₀ < a_c, mainly between ℰ_c (r_c = 0.01 R_inf) and 5 ℰ_c (r_c = 0.002 R_inf);

• there is a significant deviation from the classical expectation that most EMRIs are produced with a₀ ≃ a_c, even restricting to a₀ ≤ a_c;

• the differential rate of EMRIs has a local maximum where the overall loss cone rate has a peak, in this case at r_c ≃ 0.07R_inf, corresponding to cliffhangers.

The first two points depend mainly on the fact that S(a₀) takes approximately an order of magnitude in semi-major axis (or, equivalently, in energy) to smoothly decrease from 1 to its minimum value S_min, with S_min < S(a_c) < 1 at any value of M_•. This shows the limitations of approximating S(ℰ) as a step function. The last point, on the other hand, is exclusively caused by the novelty that S(a₀) > 0 even for a₀ ≫ a_c introduced by QS24.

4.4.2. Effects on EMRIs from lower central MBH masses

Having explored in detail the most relevant LISA case, it is also interesting to have a look at the EMRI rate implications for lower MBH masses, where cliffhanger EMRIs have been shown to dominate the S(a₀) function. Results are shown in Fig. 15 and quantified in Table 1 for the cases M_• = 10⁴, 10⁵ M_⊙. Rates are computed at t/t_p = 0.6, where t_p is the peaking time of the classical EMRI rate for the system under scrutiny (Broggi et al. 2022).

Despite cliffhangers becoming increasingly more common at lower MBH masses, their effect on the overall EMRI rates gets progressively less important. This is mainly because the peak energy of the loss cone flux, ℰ_lc, p, becomes comparable (and even bigger) then ℰ_c at low MBH masses. In the M_• = 10⁵ M_⊙ case, ℰ_c ≳ ℰ_lc, p, therefore adding a significant population of cliffhanger EMRIs at lower energies has a relevant impact on the overall EMRI population. Despite the overall EMRI rate being substantially unaffected (Table 1) we find that ∼37% of those EMRIs are cliffhanger and that the peak of the EMRI flux occurs at a > a_c (ℰ < ℰ_c) which might result in extremely eccentric systems. Conversely, the overall EMRI rate decreases in the M_• = 10⁴ M_⊙ case, with cliffhangers accounting only to ∼10% of the overall EMRI population. Notably, for such a low MBH mass, ℰ_c < ℰ_lc, p meaning that most EMRIs form well inside a_c potentially resulting in lower eccentric systems.

This analysis shows that while the impact of the detailed treatment of nuclear dynamics and the appearance of cliffhangers on the EMRI rate might be only moderate, the energy distribution of EMRIs (and hence their detailed properties) can be highly affected. In particular, a simplistic treatment of EMRIs as systems mostly forming at a = a_c fails in capturing the rich and MBH mass dependent features of this class of GW sources.

5. Discussion and conclusions

In this work we developed a formalism to locally account for the effects of two-body relaxation on the formation of EMRIs and DPs around an MBH at the core of a nuclear cluster, for which we considered a two-population model. We examined a binary system comprised of an MBH and a stellar-mass BH orbiting around it, evolving the system by using PN dynamics up to the 2.5PN term. We ran a number of simulations and studied the EMRI-to-DP ratio as a function of the initial semi-major axis a₀ of the orbit through the examination of the function S(a₀), for five different central MBH masses. Our main findings can be summarised as follows:

1. For low MBH masses the classical picture, for which there is a critical value a_c such that EMRIs only form at a₀ ≪ a_c while no EMRIs occur at a₀ ≫ a_c, breaks down. Instead, we find a significant population of EMRIs forming for large values of a₀. Such a population is consistent with the novel phenomenon of cliffhanger EMRIs, first identified in QS24.

2. Cliffhanger EMRIs are more common for low-mass MBHs because there is a larger region in the (1 − e, a) phase space which allows for a drastic reduction of the orbit in a single pericentre passage due to GW radiation. Following this analytical argument, we predict that the maximum MBH mass expected to produce cliffhanger EMRIs, given our assumption, is roughly 3 × 10⁵ M_⊙, which was verified in our numerical investigation.

3. Employing PN dynamics to evolve the system leads to more EMRIs, strongly modifying the shape of the function S(a₀) by shifting a_c to larger values and boosting the tail of cliffhanger EMRIs produced at a₀ ≫ a_c.

4. The choice of locally accounting for two-body relaxation does significantly influence the EMRI-to-DP ratio for M_• ≲ 10⁵ M_⊙, leading to a much larger fraction of cliffhanger EMRIs compared to that found in the orbit-averaged approximation.

5. Cliffhanger EMRIs account for a significant fraction of the total EMRI rate for the MBH masses where the loss cone capture differential rate peaks at a_p > a_c (up to 55% in the snapshots we tested). This is due to a large fraction of EMRIs being produced near this peak.

These results might have a significant impact on the estimates of EMRI detection rates by LISA. Currently, these estimates are still highly uncertain, ranging from a few to a few thousand EMRIs detected per year (Babak et al. 2017), and do not account for the novel population of cliffhanger EMRIs. Furthermore, we note that many cliffhanger EMRIs get to the GW emission region while grazing the loss cone edge, which might imply a significant sub-population of extremely eccentric EMRIs in the LISA band. We defer the investigation of observable cliffhanger EMRI properties by LISA to future work.

Interestingly, we observed that our results are quite sensitive to the choice of the loss cone radius R_lc. As explained in Sect. 4.3, part of the reason for the differences between the S(a₀) obtained with Newtonian and PN dynamics likely stems from switching from R_lc = 8R_g to R_lc = 5.6R_g, as plunges are more difficult to form and the cliffhanger region grows. For this reason, we believe that the closest approach to a BH on an eccentric orbit should be better explored at PN order, similarly to how it has already been done for circular orbits (e.g. Schäfer & Jaranowski 2024). Compared to Schwarzschild MBHs, spinning MBHs present innermost stable orbits which are closer or further away, respectively, depending on whether the orbit is prograde or retrograde. This again influences the loss cone radius, and consequently the EMRI-to-DP ratio. For this reason, we plan on updating our code to introduce spin effects.

Despite improving on realism by not employing orbit-averages in accounting for two-body relaxation, our approach is still limited on some fronts. Our analysis currently addresses spherical systems with a central Schwarzshild MBH. However, these results might not hold for axisymmetric or triaxial systems due to aspherical two-body relaxation and orbits naturally penetrating the loss cone (Vasiliev & Merritt 2013; Kaur & Stone 2024), or for spinning MBHs (Amaro-Seoane et al. 2013). Furthermore, the two-component Bahcall-Wolf profile we assumed might underestimate the number of stellar-mass BHs for r < R_inf, as previous works (Hopman & Alexander 2006; Broggi et al. 2022; Rom et al. 2024) have found more segregated profiles. As a result, in our simulations stars dominate two-body relaxation over the entire r range, while for a completely relaxed system stellar-mass BHs should dominate for r ≲ 0.1R_inf. While this could affect the shape of S(a₀), and in particular the value of a_c, we note that for the EMRI rate computations of Sect. 4.4 the density profiles were consistently relaxed using the code presented in Broggi et al. (2022), thus resulting in a more segregated distribution of stellar-mass BHs.

Finally, even restricting the discussion to the class of systems we are investigating, we explored a limited region of the (1 − e, a) plane. Since this choice was not biased (at least for low MBH masses) towards EMRIs or DPs, this did not influence our results on the EMRI-to-DP ratio, but it required the introduction of corrections on the counts of EMRIs for high MBH masses and limits the maximum a₀ that we can explore. This is particularly true for low MBH masses, for which we are not able to reach a₀ ≈ R_inf, as most simulated orbits starting from large a₀ result in reinitialisations, thus requiring a very large number of runs to find the EMRI-to-DP ratio within a reasonably small error. We plan on solving this issue in a follow-up work, as pushing the exploration to larger a₀ values is fundamental in order to understand whether cliffhanger EMRIs allow S(a₀) to reach a plateau S(∞)≠0 or not.

---

Acknowledgments

We thank Ismail Qunbar and Nicholas C. Stone for useful discussions and willingness to share data with us. DM acknowledges that this publication was produced while attending the PhD program in Space Science and Technology at the University of Trento, Cycle XXXIX, with the support of a scholarship financed by the Ministerial Decree no. 118 of 2nd March 2023, based on the NRRP – funded by the European Union – NextGenerationEU – Mission 4 “Education and Research”, Component 1 “Enhancement of the offer of educational services: from nurseries to universities” – Investment 4.1 “Extension of the number of research doctorates and innovative doctorates for public administration and cultural heritage” - CUP E66E23000110001. MB acknowledges support provided by MUR under grant “PNRR - Missione 4 Istruzione e Ricerca – Componente 2 Dalla Ricerca all’Impresa – Investimento 1.2 Finanziamento di progetti presentati da giovani ricercatori ID:SOE_0163” and by University of Milano-Bicocca under grant “2022-NAZ-0482/B”. AS acknowledges financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691).

References

Appendix A: Diffusion coefficients

In this work, we consider a stellar-mass BH of mass m_• and velocity magnitude v having a sequence of two-body encounters with background bodies, grouped into populations (identified by the index α) based on their masses m_α. We decompose its velocity change as the sum of Δv_∥, parallel to the initial line of motion, and Δv_⊥, in the plane orthogonal to the original direction of motion.

We use the notation D[X] to indicate the expectation value of the generic quantity X per unit of time. The diffusion coefficients we use in this work are (Chandrasekhar 1943; Binney & Tremaine 2008; Merritt 2013a):

$$\begin{array}{cc}\hfill \mathrm{D}\left[\mathrm{\Delta }{v}_{\Vert }\right]& =-16{\pi }^2{G}^{2}ln\mathrm{\Lambda }{\displaystyle \sum _{\alpha }}{m}_{\alpha }(m_{\bullet }+m_{\alpha }){\mathcal{I}}_{2,\alpha }(r,v),\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \mathrm{D}\left[\mathrm{\Delta }{v}_{\perp }\right]& =0,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill \mathrm{D}\left[{\left(\mathrm{\Delta }{v}_{\Vert }\right)}^2\right]& =\frac{32{\pi }^2}{3}{G}^{2}ln\mathrm{\Lambda }v{\displaystyle \sum _{\alpha }}{m}_{\alpha }^2({\mathcal{I}}_{4,\alpha }(r,v)+{\mathcal{J}}_{1,\alpha }(r,v))\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill \mathrm{D}\left[{\left(\mathrm{\Delta }{v}_{\perp }\right)}^2\right]& =\begin{array}{cc}& \frac{32{\pi }^2}{3}{G}^{2}ln\mathrm{\Lambda }v\times \hfill \\ \hfill & \times {\displaystyle \sum _{\alpha }}{m}_{\alpha }^2(3{\mathcal{I}}_{2,\alpha }(r,v)-{\mathcal{I}}_{4,\alpha }(r,v)+2{\mathcal{J}}_{1,\alpha }(r,v)).\hfill \end{array}\hfill \end{array}$$(A.4)

Here lnΛ is the Coulomb logarithm, which we set to 10 as an order of magnitude approximation (Spinnato et al. 2003; Merritt 2013b), while the functions ℐ_k, α and 𝒥_k, α are defined as

$$\begin{array}{cc}\hfill {\mathcal{I}}_{k,\alpha }(r,v)& ={\displaystyle {\int }_0^v}\mathrm{d}{v}^{\prime }{\left(\frac{v^{\prime }}{v}\right)}^k{f}_{\alpha }(r,v^{\prime }),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\mathcal{J}}_{k,\alpha }(r,v)& ={\displaystyle {\int }_v^{\infty }}\mathrm{d}{v}^{\prime }{\left(\frac{v^{\prime }}{v}\right)}^k{f}_{\alpha }(r,v^{\prime }),\hfill \end{array}$$(A.6)

with f_α being the distribution function of objects in population α on the six-dimensional phase-space of position and velocity. Assuming that the distribution of scatterers is spherical, and approximating it as locally isotropic in velocity, f_α depends only on the magnitudes of its arguments r and v.

In practice, the diffusion coefficients tell us how the velocity of the stellar-mass BH changes per unit time as it moves with velocity v through the distributions of background bodies at distance r from their centres. Notably, the perpendicular variation of the velocity D[Δv_⊥] vanishes since the stellar-mass BH is equally likely to be deflected along any of the infinite directions perpendicular to its motion, while the parallel variation D[Δv_∥] < 0 describes the effect of dynamical friction (Chandrasekhar 1943).

At Newtonian order, $v=\sqrt{2(\psi -\mathcal{E})}$ for bound orbits. Under this approximation, we can express ℐ_k, α and 𝒥_k, α as functions of ψ(r) and ℰ(r, v):

$$\begin{array}{cc}\hfill {\mathcal{I}}_{k,\alpha }(\psi ,\mathcal{E})& =\frac{1}{\sqrt{2(\psi -\mathcal{E})}}{\displaystyle {\int }_{\mathcal{E}}^{\psi }}\mathrm{d}{\mathcal{E}}^{\prime }{\left(\frac{\psi -{\mathcal{E}}^{\prime }}{\psi -\mathcal{E}}\right)}^{(k-1)/2}{f}_{\alpha }({\mathcal{E}}^{\prime }),\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill {\mathcal{J}}_{k,\alpha }(\psi ,\mathcal{E})& =\frac{1}{\sqrt{2(\psi -\mathcal{E})}}{\displaystyle {\int }_0^{\mathcal{E}}}\mathrm{d}{\mathcal{E}}^{\prime }{\left(\frac{\psi -{\mathcal{E}}^{\prime }}{\psi -\mathcal{E}}\right)}^{(k-1)/2}{f}_{\alpha }({\mathcal{E}}^{\prime }).\hfill \end{array}$$(A.8)

Here f depends only on ℰ as we assume that the distribution function is ergodic, meaning it depends on r and v only through the Hamiltonian: f(r, v) = f(H(r, v)).

The ergodic distribution function corresponding to a spherical density profile can be computed using Eddington’s formula (Eddington 1916; Merritt 2013a):

$$\begin{array}{c}\hfill f_{\alpha }(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}(\frac{1}{\sqrt{\mathcal{E}}}{\frac{\mathrm{d}{n}_{\alpha }}{\mathrm{d}\psi }|}_{\psi =0}+{\displaystyle {\int }_0^{\mathcal{E}}}\frac{\mathrm{d}\psi }{\sqrt{\mathcal{E}-\psi }}\frac{{\mathrm{d}}^2{n}_{\alpha }}{\mathrm{d}{\psi }^2}).\end{array}$$(A.9)

We can simplify this expression by observing that ψ = 0 when r → ∞. For a Dehnen density profile, at infinite distance we have that n_α ∝ ψ⁴, since n_α ∝ r⁻⁴ and ψ ∝ r⁻¹. Thus,

$$\begin{array}{c}\hfill {\frac{\mathrm{d}{n}_{\alpha }}{\mathrm{d}\psi }|}_{\psi =0}={4{\psi }^3|}_{\psi =0}=0,\end{array}$$(A.10)

which gives us the simpler expression:

$$\begin{array}{c}\hfill f_{\alpha }(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}{\displaystyle {\int }_0^{\mathcal{E}}}\mathrm{d}\psi g(\psi ),\end{array}$$(A.11)

where

$$\begin{array}{c}\hfill g(\psi )=\frac{1}{\sqrt{\mathcal{E}-\psi }}\frac{{\mathrm{d}}^2{n}_{\alpha }}{\mathrm{d}{\psi }^2}.\end{array}$$(A.12)

Finally, we can observe that:

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{n}_{\alpha }}{\mathrm{d}{\psi }^2}=\frac{\mathrm{d}r}{\mathrm{d}\psi }\frac{\mathrm{d}}{\mathrm{d}r}\left(\frac{\mathrm{d}r}{\mathrm{d}\psi }\frac{\mathrm{d}}{\mathrm{d}r}{n}_{\alpha }\right)=\frac{1}{{\left({\psi }^{\prime }\right)}^2}{n}_{\alpha }^{″}-\frac{{\psi }^{″}}{{\left({\psi }^{\prime }\right)}^3}{n}_{\alpha }^{\prime },\end{array}$$(A.13)

where the prime symbol indicates derivation with respect to r.

In practice, to find f_α we have to integrate g(ψ) over ψ between 0 and ℰ. We perform this numerically by evaluating the function g(ψ) at discrete values ψ_i inside the range [0, ℰ]. Thus, to use the expression for the second derivative of n_α that we just found, for each of these values we first numerically find the r_i such that ψ(r_i) = ψ_i by solving Eq. (5). Then, we can compute ψ′(r_i), ψ″(r_i), n′(r_i), and n″(r_i), as their analytic expressions can be easily obtained from Eqs. (3), (4), and (5).

In order to estimate the orbit-averaged relaxation timescale, we compute the orbit averaged diffusion coefficient D_o[ΔR] as

$$\begin{array}{c}\hfill {\mathrm{D}}_{\mathrm{o}}[\mathrm{\Delta }\mathrm{R}]={\displaystyle \underset{R\to 0}{lim}}\frac{D_{\mathit{RR}}(\mathcal{E},R)}{4{\pi }^2{J}_{\mathrm{c}}^2(\mathcal{E})P(\mathcal{E},R)},\end{array}$$(A.14)

where D_RR is one of the coefficients of the orbit-averaged Fokker-Planck equation in the flux conservation form. It can be computed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_{\mathit{RR}}& =\frac{256}{3}{\pi }^4{G}^{2}log\mathrm{\Lambda }R{J}^2\times \hfill \\ \hfill & {\displaystyle \sum _{\alpha }}{m}_{\alpha }^2{\displaystyle {\int }_{r_{-}}^{r_{+}}}\frac{v\mathrm{d}r}{v_{\mathrm{r}}}[(\frac{2}{v_{\mathrm{t}}^2}+\frac{2v^2}{v_{\mathrm{c}}^4}-\frac{4}{v_{\mathrm{c}}^2}){\mathcal{J}}_{1,\alpha }+(\frac{3}{v_{\mathrm{t}}^2}-\frac{3}{v^2}){\mathcal{I}}_{2,\alpha }\hfill \\ \hfill & +(\frac{3}{v^2}-\frac{1}{v_{\mathrm{t}}^2}+\frac{2v^2}{v_{\mathrm{c}}^4}-\frac{4}{v_{\mathrm{c}}^2}){\mathcal{I}}_{4,\alpha }],\hfill \end{array}\end{array}$$(A.15)

where v_c = J_c/r is the circular velocity at a given energy.

Appendix B: From δv_∥ and δv_⊥ to δℰ and δJ

The aim of this appendix is to show how we can express the change in energy and angular momentum of the orbiting stellar-mass BH as a function of its velocity change in the parallel and perpendicular directions relative to its original motion.

At first, the velocity of the stellar-mass BH is

$$\begin{array}{c}\hfill \mathit{v}=v_{\mathrm{r}}{\widehat{\mathit{e}}}_{\mathbf{r}}+v_{\mathrm{t}}{\widehat{\mathit{e}}}_{\mathbf{t}},\end{array}$$(B.1)

where ${\widehat{\mathit{e}}}_{\mathbf{r}}$ and ${\widehat{\mathit{e}}}_{\mathbf{t}}$ respectfully represent the radial and tangential unit vectors of the orbit.

The variation of velocity δv due to a stochastic velocity kick is

$$\begin{array}{c}\hfill \delta \mathit{v}=(\delta v_{\Vert }\frac{v_{\mathrm{r}}}{v}-\delta v_{\perp }sin\theta \frac{v_{\mathrm{t}}}{v}){\widehat{\mathit{e}}}_{\mathbf{r}}+(\delta v_{\Vert }\frac{v_{\mathrm{t}}}{v}+\delta v_{\perp }sin\theta \frac{v_{\mathrm{r}}}{v}){\widehat{\mathit{e}}}_{\mathbf{t}}+\delta v_{\perp }cos\theta {\widehat{\mathit{e}}}_{\mathbf{3}},\end{array}$$(B.2)

where ${\widehat{\mathit{e}}}_{\mathbf{3}}$ is the unit vector in the direction perpendicular to the orbital plane and θ is the angle between a direction perpendicular to the motion of the stellar-mass BH and ${\widehat{\mathit{e}}}_{\mathbf{3}}$. In practice, in our implementation this angle is always uniformly selected in the range [0, 2π), thus we do not properly define the direction of δv_⊥ here.

We can now compute the change in energy given δv_∥ and δv_⊥:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \delta \mathcal{E}& =-\frac{{(\mathit{v}+\delta \mathit{v})}^2}{2}+\frac{{\mathit{v}}^2}{2}=-\frac{{(\delta \mathit{v})}^2}{2}+\mathit{v}·\delta \mathit{v}\hfill \\ \hfill & =-v\delta v_{\Vert }-\frac{1}{2}{(\delta v_{\Vert })}^2-\frac{1}{2}{(\delta v_{\perp })}^2.\hfill \end{array}\end{array}$$(B.3)

The vector change of J is instead

$$\begin{array}{c}\hfill \delta \mathit{J}=r{\widehat{\mathit{e}}}_{\mathbf{r}}\times \delta \mathit{v}=-r\delta v_{\perp }sin\theta {\widehat{\mathit{e}}}_{\mathbf{t}}+r(\frac{v_{\mathrm{t}}}{v}\delta v_{\Vert }+\frac{v_{\mathrm{r}}}{v}\delta v_{\perp }cos\theta ){\widehat{\mathit{e}}}_{\mathbf{3}},\end{array}$$(B.4)

from which follows that J = |J| changes as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \delta J& =|\mathit{J}+\delta \mathit{J}|-|\mathit{J}|\hfill \\ \hfill & =J(\frac{\delta v_{\Vert }}{v}+\delta v_{\perp }\frac{v_{\mathrm{r}}}{v{v}_{\mathrm{t}}}sin\theta +\frac{1}{2}{\left(\frac{\delta v_{\perp }cos\theta }{v_{\mathrm{t}}}\right)}^2)+\mathcal{O}\left({(\delta v)}^3\right),\hfill \end{array}\end{array}$$(B.5)

where we used J = rv_t.

More detailed steps can be found in Broggi (2024).¹³

Appendix C: From Δℰ and ΔJ to Δv

The aim of this appendix is to show how we can kick the stellar-mass BH to recover a given change in energy and angular momentum. Crucially, we wish to achieve this only by changing the velocity of the stellar-mass BH, while keeping its position the same.

At first, the velocity of the stellar-mass BH is

$$\begin{array}{c}\hfill \mathit{v}=v_{\mathrm{r}}{\widehat{\mathit{e}}}_{\mathbf{r}}+v_{\mathrm{t}}{\widehat{\mathit{e}}}_{\mathbf{t}},\end{array}$$(C.1)

where ${\widehat{\mathit{e}}}_{\mathbf{r}}$ and ${\widehat{\mathit{e}}}_{\mathbf{t}}$ respectfully represent the radial and tangential unit vectors of the orbit. The most generic velocity kick we can give to the stellar-mass BH is

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathit{v}=\mathrm{\Delta }{v}_{\mathrm{r}}{\widehat{\mathit{e}}}_{\mathbf{r}}+\mathrm{\Delta }{v}_{\mathrm{t}}{\widehat{\mathit{e}}}_{\mathbf{t}}+\mathrm{\Delta }{v}_3{\widehat{\mathit{e}}}_{\mathbf{3}},\end{array}$$(C.2)

where ${\widehat{\mathit{e}}}_{\mathbf{3}}$ is the unit vector in the direction perpendicular to the orbital plane. Thus, the new velocity of the stellar-mass BH v₁ is

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathbf{1}}=(v_{\mathrm{r}}+\mathrm{\Delta }{v}_{\mathrm{r}}){\widehat{\mathit{e}}}_{\mathbf{r}}+(v_{\mathrm{t}}+\mathrm{\Delta }{v}_{\mathrm{t}}){\widehat{\mathit{e}}}_{\mathbf{t}}+\mathrm{\Delta }{v}_3{\widehat{\mathit{e}}}_{\mathbf{3}}.\end{array}$$(C.3)

In Newtonian dynamics, J = r × v. Thus,

$$\begin{array}{c}\hfill \mathit{J}=r{v}_{\mathrm{t}}{\widehat{\mathit{e}}}_{\mathbf{3}}\end{array}$$(C.4)

and

$$\begin{array}{c}\hfill {\mathit{J}}_{\mathbf{1}}=\mathit{r}\times {\mathit{v}}_{\mathbf{1}}=-r\mathrm{\Delta }{v}_3{\widehat{\mathit{e}}}_{\mathbf{t}}+r(v_{\mathrm{t}}+\mathrm{\Delta }{v}_{\mathrm{t}}){\widehat{\mathit{e}}}_{\mathbf{3}},\end{array}$$(C.5)

where J₁ represents the angular momentum after the change in velocity. Moreover, the energy will change from

$$\begin{array}{c}\hfill \mathcal{E}=-\frac{v^2}{2}+\psi (r)\end{array}$$(C.6)

to

$$\begin{array}{c}\hfill {\mathcal{E}}_1=-\frac{v_1^2}{2}+\psi (r).\end{array}$$(C.7)

Now, we write the following system:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle {\mathcal{E}}_1=\mathcal{E}+\mathrm{\Delta }\mathcal{E}}\hfill \\ {\displaystyle J_1=J+\mathrm{\Delta }J}\hfill \\ {\displaystyle \frac{{\mathit{J}}_{\mathbf{1}}·\mathit{J}}{J_{1}J}=cos\beta }\hfill \end{array},\end{array}$$(C.8)

where β is the angle between vectors J₁ and J. Using the relations found above, we can rewrite the system as

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle -\frac{1}{2}\sqrt{{(v_{\mathrm{r}}+\mathrm{\Delta }{v}_{\mathrm{r}})}^2+{(v_{\mathrm{t}}+\mathrm{\Delta }{v}_{\mathrm{t}})}^2+{(\mathrm{\Delta }{v}_3)}^2}=-\frac{1}{2}{v}^2+\mathrm{\Delta }\mathcal{E}}\hfill \\ {\displaystyle r\sqrt{{(\mathrm{\Delta }{v}_3)}^2+{(v_{\mathrm{t}}+\mathrm{\Delta }{v}_{\mathrm{t}})}^2}=J+\mathrm{\Delta }J}\hfill \\ {\displaystyle \frac{r(v_{\mathrm{t}}+\mathrm{\Delta }{v}_{\mathrm{t}})}{J+\mathrm{\Delta }J}=cos\beta }\hfill \end{array}.\end{array}$$(C.9)

In order to solve this system, we need to know the value of β. In practice, during our simulations we randomly drew it from the interval [0, 2π), thus assuming an isotropic distribution.

Finally, after a sequence of substitutions, the system above can be written as

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle \mathrm{\Delta }{v}_{\mathrm{r}}=\sqrt{v^2-2\mathrm{\Delta }\mathcal{E}-{\left(\frac{J+\mathrm{\Delta }J}{r}\right)}^2}-v_{\mathrm{r}}}\hfill \\ {\displaystyle \mathrm{\Delta }{v}_{\mathrm{t}}=\frac{(J+\mathrm{\Delta }J)cos\beta }{r}-v_{\mathrm{t}}}\hfill \\ {\displaystyle \mathrm{\Delta }{v}_3=\frac{(J+\mathrm{\Delta }J)sin\beta }{r}}\hfill \end{array}.\end{array}$$(C.10)

All Tables

All Figures