GW emission and relativistic dynamical friction in intermediate mass ratio inspirals

TL;DR

This paper presents preliminary simulations of intermediate mass ratio inspirals within dark matter spikes, incorporating relativistic effects and dynamical friction corrections, revealing significant impacts on inspiral times.

Contribution

It introduces the first simulations including relativistic and dynamical friction corrections for IMRIs in dark matter spikes.

Findings

Relativistic corrections halve the inspiral time compared to classical estimates.

Neglecting dynamical friction results in longer inspiral times.

Relativistic effects are crucial for accurate modeling of IMRIs in dark matter environments.

Abstract

We present a set of preliminary simulations of intermediate mass ratio inspirals (IMRIs) inside dark matter (DM) spikes accounting for post-Newtonian corrections the interaction between the two black holes up to the order 2.5 in $c^2$, as well as relativistic corrections to the dynamical friction (DF) force exerted by the DM distribution. We find that, incorporating relativity reduces of a factor $1/2$ the inspiral time, for equivalent initial orbital parameters, with respect to the purely classical estimates. Vice versa, neglecting the DF of the spike systematically yields longer inspiral times.

Equations

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}={\mathbf{F}}_{\mathrm{tot}}-\eta \mathbf{v},$$

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

$${\rho }_{\mathrm{DM}}(r)={\rho }_6{\left(\frac{r}{r_6}\right)}^{-\beta }{\left(1+\frac{r}{r_t}\right)}^{-\alpha },$$

$$\eta \approx 4\pi G^2{m}_{\mathrm{BH}}{\rho }_{\mathrm{DM}}(r)\frac{\log \Lambda }{{\sigma }^3{\gamma }^3},$$

Full text

††editors: Hyung Mok Lee, Rainer Spurzem and Jongsuk Hong, eds.

1 Introduction

Massive black holes (BHs) are likely to be sitting at the centre of dark matter (DM) structures such as spikes or mounds that may or may not evaporate across cosmic time according to the specific nature of DM itself or the black hole mass.

It has been speculated that future gravitational waves (GW) observatories can probe the DM distribution around intermediate or supermassive BHs with masses $M_{\rm BH}$ ranging from $10^{2}$ to $10^{9}{\rm M}_{\odot}$, and possibly obtain constraints on the nature of DM, by detecting the dephasing in the waveform of binary BH mergers induced by said DM overdensities (see Kavanagh et al., 2020). The accurate modelling of the dynamical friction (Chandrasekhar, 1943, DF), induced by the DM on the secondary compact object of mass $m_{\rm BH}\approx 1{\rm M}_{\odot}$ in intermediate (IMRI) and extreme mass ratio inspirals (EMRI) is therefore of pivotal importance.

In Karydas et al., (2025) and Kavanagh et al., (2025) we introduced a simple semi-analytical scheme to account for the effect of DF and the back reaction of the inspiralling compact object on the DM distribution and demonstrate that including stirring (i.e. cusp density fluctuations induced by the time-dependent potential of the binary BH) tends to slow the rate of dark matter depletion and therefore enhances the impact of the drag exerted by DM on the dynamics of the binary. Additionally, we validated our models using the numerical code NbodyIMRI, a new publicly available code designed for simulating binary systems within DM spikes (Kavanagh et al., 2024) that models the DM distribution as a set of $N$ non-interacting simulation particles with tunable mass and force acting on $m_{\rm BH}$.

So far, our treatment was purely classical, while it is of course natural that, in particular in the last phases of the DF driven inspirals, relativistic effects do play a non-negligible role. Here we introduce another simple numerical code incorporating corrections to both the DF induced by the DM spike and the force field of the central massive BH.

2 Models

Following the approach of Dai et al., (2022) and Trani and Di Cintio, (2025), 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

[TABLE]

using the Mikkola and Merritt, (2006) 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

[TABLE]

with the coefficients of the post-Newtonian (PN) expansion in powers of the speed of light $c$, $\mathcal{A}$ and $\mathcal{B}$, summed up to the order 2.5 (for the explicit expression see e.g. Mora and Will, 2004), 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

[TABLE]

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, (2023)

[TABLE]

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., (2025) 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. $\mathcal{A}=\mathcal{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.

We observe that for non vanishing $e$ both DF and GW emission-induced dissipation tend to circularize the orbit, however at about $r\approx 20r_{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,, 2019, 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 $\mathcal{A}$ and $\mathcal{B}$) and hence $t_{f}$ is further reduced with respect to the classical case for analogous initial orbital parameters.

As a general trend, we confirm the result of our $N-$body simulations (and those of similar works such as Mukherjee et al., 2024) 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., 2022), 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., 2024).