2. Models

Following the approach of Dai et al. (Reference Dai, Gong, Jiang and Liang2022) and Trani and Di Cintio (Reference Trani and Di Cintio2025), we integrate the orbit of the smaller companion $m_{\rm BH}$ in the central field of the massive BH $M_{\rm BH}$ , accounting for the effect of DF as parametrized by the friction coefficient $\eta$ given by

(1) \begin{equation} \frac{{\rm d}^2\mathbf{r}}{{\rm d}t^2}=\mathbf{F}_{\rm tot}-\eta\mathbf{v},\end{equation}

using the Mikkola and Merritt (Reference Mikkola and Merritt2006) modified midpoint leapfrog scheme with adaptive timestep. In Eq. (1) $\mathbf{F}_{\rm tot}=\mathbf{F}_{\rm BH}-\nabla\Phi_{\rm DM}$ is the total force per unit mass. The contribution of the central BH is

(2) \begin{equation} \mathbf{F}_{\rm BH}=-\frac{GM_{\rm BH}}{r^2}\bigg[(1+\mathscr{A})\frac{\mathbf{r}}{r}+\mathscr{B}\mathbf{v}\bigg],\end{equation}

with the coefficients of the post-Newtonian (PN) expansion in powers of the speed of light c, $\mathscr{A}$ and $\mathscr{B}$ , summed up to the order 2.5 (for the explicit expression see e.g. Mora and Will Reference Mora and Will2004), corresponding to the lowest order accounting for energy loss in GW emission. The DM cusp generates the potential $\Phi_{\rm DF}$ via the usual (classical) Poisson equation associated to the density

(3) \begin{equation}\rho_{\mathrm{DM}}(r) = \rho_6 \left(\frac{r}{r_6}\right)^{-\beta} \left( 1 + \frac{r}{r_t}\right)^{-\alpha},\end{equation}

where $r_6 = 10^{-6}\,\mathrm{pc}$ is a reference radius; $\rho_6$ specifies the density normalization (at said radius); $\alpha$ and $r_t$ are the slope and the cutoff radius of the truncation, and $\beta$ is the inner density slope.

The DM density also produces DF on the companion $m_{\rm BH}$ , here we use the approximate form of the relativistic coefficient evaluated in Chiari and Di Cintio (Reference Chiari and Di Cintio2023)

(4) \begin{equation} \eta\approx 4\pi G^2m_{\rm BH}\rho_{\mathrm{DM}}(r)\frac{\log\Lambda}{\sigma^3\gamma^3},\end{equation}

where $\gamma=(1+v^2/c^2)^{-1/2}$ and we have assumed that the velocity dispersion $\sigma$ of the DM particles in the cusp locally equals the the circular velocity in the (classical) Keplerian potential of $M_{\rm BH}$ . The argument of the Coulomb logarithm $\Lambda=b_{\rm max}/b_{\rm min}$ is constrained as in Kavanagh et al. (Reference Kavanagh, Karydas, Bertone, Di Cintio and Pasquato2025) with $b_{\max}\approx r/3$ and $b_{\rm min}\approx Gm_{\rm BH}/v^2$ .

3. Preliminary results and implications

In the same spirit of our previous work with $N-$ body simulations, we integrated the orbit of the companion around an intermediate mass BH with $M_{\rm BH}=10^3{\rm M}_\odot$ , with $m_{\rm BH}$ ranging from one to $30{\rm M}_\odot$ . We fixed the cusp parameters in Eq. (3) to $\alpha=3$ ; $\beta=2.33$ and $\rho_6=10^{16}{\rm M_\odot {\rm pc}^{-3}}$ . For this specific choice of cusp parameters, the associated gravitational potential is less than one part in $10^2$ at a separation of $10^{-3}$ astronomical units, so that its contribution to $\mathbf{F}_{tot}$ could be in principle neglected. We considered different values of the orbital eccentricity e, from a moderately eccentric orbits with $e=0.5$ to circular orbits. All integrations of Eq. (1) where extended up to the time at which the binary separation reaches ten times the Schwartzschild radius of the central BH $r_S=2GM_{\rm BH}/c^2$ .

Figure 1 shows the evolution of e and the semimajor axis a for two sample initial conditions and three different simulation set-ups. Namely the purely classical case (i.e. $\mathscr{A}=\mathscr{B}=0$ ; $\gamma=1$ ), in a way corresponding to the simulations in NbodyIMRI; the full relativistic case, and an additional relativistic case where the contribution of the DM cusp is neglected for both DF and potential.

Figure 1.

Evolution of the orbital eccentricity e (top panels) and the semimajor axis a (bottom panels) in the last phase of a IMRI for a 20 $M_\odot$ black hole orbiting an IMBH of $10^3M_\odot$ on a nearly circular orbit (left panels) and a moderately eccentric orbit (right panels).

We observe that for non vanishing e both DF and GW emission-induced dissipation tend to circularize the orbit, however at about $r\approx20r_S$ the PN corrections to the gravitational acceleration result in a drastic increase of the eccentricity. This latter behavior is expected because of the Newtonian definition of eccentricity that we employ (see Will Reference Will2019, and reference therein). If on one hand in classical DF-driven inspirals the decay of the semimajor axis a is always nearly linear, once the PN terms are dominant the decay rate increases with a steep power-law resulting in a reduction of the infall time $t_f$ of a factor of about $1/2$ , as shown in Figure 2 for a circular BH of mass 1, 3, 10 and 30 ${\rm M}_\odot$ on an initially circular orbit of semimajor axis $6\times 10^{-3}$ astronomical units. Larger values of the companion’s mass result in stronger relativistic corrections (since $m_{\rm BH}$ enters explicitly in the definitions of $\mathscr{A}$ and $\mathscr{B}$ ) and hence $t_f$ is further reduced with respect to the classical case for analogous initial orbital parameters.

Figure 2.

Evolution of the semimajor axis a in the last phase of an IMRI for different black hole masses ( $m_{BH}=1,$ 3, 10 and $30{\rm M}_\odot$ ) orbiting an IMBH of $10^3{\rm M}_\odot$ on a (initially) nearly circular orbit (left panel). The solid lines refer to the relativistic simulations with up to 2.5PN terms while the dashed lines mark the parent purely classical integration. Ratio of the inspiral time down to $10r_S$ in relativistic and classical simulations as function of $m_{BH}$ (right panel).

As a general trend, we confirm the result of our $N-$ body simulations (and those of similar works such as Mukherjee et al. Reference Mukherjee, Holgado, Ogiya and Trac2024) that the presence of a DM spike surrounding the primary component of the binary significantly enhances the inspiral (cfr cyan and orange curves in Fig. 1), even in the cases when its gravitational field is negligible. We stress the fact that our simple model could be also applied to a broad range of mass ratios $m_{\rm BH}/M_{\rm BH}$ , in principle allowing to simulate EMRI in galactic cores. This established, the observability of dynamical friction effects in the gravitational wave band remains to be determined, which will be addressed in our follow-up work. Moreover, as the BH mass ratios explored here are likely to be associated to a non-negligible DM profile evolution (see e.g Coogan et al. Reference Coogan, Bertone, Gaggero, Kavanagh and Nichols2022), an extension of the numerical method that incorporates this feature is currently in the works. In particular, this approach would become relevant to probe alternative models of DM, such as for example the superfluid DM, where DF departs significantly from the Chandrasekhar formulation (Berezhiani et al. Reference Berezhiani, Cintia, De Luca and Khoury2024).