Examining the Effects of Dark Matter Spikes on Eccentric Intermediate Mass Ratio Inspirals Using 𝑁 -body Simulations

Recent studies have postulated that the presence of dark matter (DM) spikes around IMBHs could lead to observable dephasing effects in gravitational wave (GW) signals emitted by Intermediate Mass Ratio Inspirals (IMRIs). While prior investigations primarily relied on non-self-consistent analytic methods to estimate the influence of DM spikes on eccentric IMRIs, our work introduces the first self-consistent treatment of this phenomenon through 𝑁 -body simulations. Contrary to previous studies, which suggested that dynamical friction (DF), a cumulative effect of two-body encounters, is the primary mechanism responsible for energy dissipation, we reveal that the slingshot mechanism, a three-body effect, plays a more significant role in driving the binary system’s energy loss and consequent orbital shrinkage, similar to stellar loss cone scattering in Massive Black Hole (MBH) binaries. Furthermore, our work extends its analysis to include rotation in DM spikes, a factor often overlooked in previous studies. We find that binaries that counter-rotate with respect to the spike particles merge faster, while binaries that co-rotate merge slower, in opposition to the expectation from DF theory. While our models are idealistic, they offer findings that pave the way for a more comprehensive understanding of the complex interactions between DM spikes, IMRIs, GW emission, and the ability to constrain DM microphysics. Our work systematically includes Post-Newtonian (PN) effects until 2.5 order and our results are accurate and robust.


INTRODUCTION
The nature of dark matter (DM) remains one of the most pressing mysteries in modern astrophysics.While its existence has been inferred through a range of indirect observations such as galaxy rotational curves (e.g., Rubin & Ford 1970;Rubin et al. 1978Rubin et al. , 1980;;Persic et al. 1996), and gravitational lensing of galaxy clusters (e.g., Tyson et al. 1990;Hammer 1991;Fort et al. 1992;Le Fevre et al. 1994), its properties and interactions remain largely unknown.Cosmological simulations suggest that DM resides in galactic halos and is distributed according to the Navarro-Frenk-White profile (Navarro et al. 1996;Bertone 2010, see second reference for a review).Galaxies, and therefore galactic halos, often contain massive black holes (MBHs) at the center.Gondolo & Silk (1999) proposed that the adiabatic growth of these MBHs would modify the dark matter profile near them and create extremely dense density spikes.The density profile of these DM spikes  DM is parameterized as a power-law profile and is given as (e.g., Gondolo & Silk 1999;Eda et al. 2013;Kavanagh et al. 2020) where  denotes the distance from the central MBH,  sp is a normalization factor for the density profile,  sp characterizes the power-★ E-mail: diptajym@andrew.cmu.edulaw profile of the spike and  sp is a characteristic radius of the spike.Usually  sp is taken to be between 2 to 2.5.For a central MBH with a mass of 10 3  ⊙ , under these parameters,  sp ≈ 0.5pc and the density 10 −6 pc from the MBH is ∼ 10 15  ⊙ pc −3 .
Although the authors showed that such spikes would typically form around Supermassive Black Holes (SMBHs) with   ≥ 10 6  ⊙ , where   is the mass of the black hole (BH), later work (e.g., Ullio et al. 2001;Merritt et al. 2002;Merritt 2004;Bertone & Merritt 2005) found that spikes could be depleted via different astrophysical processes such as galaxy mergers, scattering with nearby stars, and off center seed BHs.While SMBHs are expected to undergo many mergers, the same cannot be said about Intermediate Mass Black Holes.DM spikes around Intermediate Mass Black Holes (IMBHs) with 10 3 ≤   ≤ 10 6 are expected to survive as they undergo fewer mergers (Zhao & Silk 2005;Bertone et al. 2005).Additionally, Ferrer et al. (2017) showed that spinning IMBHs could form denser spikes.Thus, IMBHs are expected to be the predominant source for black holes surrounded by DM spikes.
The advent of gravitational wave (GW) astronomy using LIGO-Virgo interferometers (e.g., Abbott et al. 2017Abbott et al. , 2020a,b,c) ,b,c) and pulsar timing array (PTA) (e.g., Mingarelli et al. 2017;Kelley et al. 2018;Agazie et al. 2023a,b) has opened up new avenues to detect and explore the properties of DM.In particular, it has been suggested that the effects of a DM spike could be imprinted on the GW signal of an Intermediate Mass Ratio Inspiral (IMRI) and has drawn a lot of interest (e.g., Eda et al. 2013;Eda et al. 2015;Yue & Han 2018;Yue et al. 2019;Kavanagh et al. 2020;Becker et al. 2022;Dai et al. 2022).A visual representation of such a system can be seen in Figure 1.IMRIs, formed from the inspiral of compact objects such as white dwarfs or neutron stars, or stellar mass black holes into IMBHs, will be detectable by future space based mHz GW detectors like LISA (Amaro-Seoane et al. 2017) and TianQin (Luo et al. 2016) or deciHz GW detectors like DECIGO (Kawamura et al. 2019) or MAGIS (Abe et al. 2021).GWs emanating from such systems are expected to be in the LISA and DECIGO band for periods of months to years allowing environmental effects, including DM spikes, to play a major role in modifying the signal.Since matched filtering relies on a careful determination of the simulated signal to a few cycles over hundreds to thousands of cycles, it is imperative to take into account the environmental effects while numerically calculating the simulated GW signal (e.g., Eda et al. 2013;Macedo et al. 2013;Zwick et al. 2022;Coogan et al. 2022;Baumann et al. 2022;Zwick et al. 2023;Cole et al. 2023a,b) Eda et al. (2013) and Eda et al. (2015) first proposed that the gravitational effects of the spike could leave an imprint on the GW signal of an IMRI.The spike particles would get scattered by the inspiraling object leading to an extra drag force experienced by the inspiraling object.This drag force has been attributed to dynamical friction (DF; Chandrasekhar 1943) on the inspiraling object.Although the DM mass contained in the spike is quite small relative to the mass of the central IMBH and the inspiraling object, it has been shown that the gravitational drag can lead to a phase shift in the GW signal over thousands of cycles.Recent studies have shown that dephasing due to DF can lead to anywhere between 10 3 − 10 7 fewer cycles which would be above the signal to noise detector theshold for LISA and could be detected and distinguished as an imprint of the DM spike surrounding the IMBH (e.g., Eda et al. 2015;Kavanagh et al. 2020;Becker et al. 2022;Dai et al. 2022).Different models of DM can lead to different spike parameters and as such change the amount of dephasing itself, which can be detected and distinguished.As such, GW astronomy provides a unique pathway to ascertain the presence and nature of DM itself.Detection of even one DM spike modified signal would be effective towards eliminating other theories of gravity such as modified Newtonian dynamics and can help constrain the particle nature of DM (Hannuksela et al. 2020).
A careful determination of the the expected dephasing, therefore, is vital.While initial studies (e.g., Eda et al. 2013;Eda et al. 2015) used a static DM background, Kavanagh et al. (2020) showed that the inspiral of the IMRI can inject energy into the spike itself, often comparable to the binding energy of the spike.As such the feedback of the IMRI onto the spike cannot be ignored.Through a semianalytic framework called HaloFeedback, Kavanagh et al. (2020) showed that the including the back reaction from the IMRI can lead upto 100× reduction in expected amount of dephasing.Therefore, it is necessary to self-consistently follow the dynamics of the binary in order to calculate the amount of dephasing.While HaloFeedback relies on the assumption that the binary is on a circular orbit at all times, studies (e.g., Yue & Han 2018;Yue et al. 2019;Cardoso et al. 2021;Dai et al. 2022) showed that such an assumption might not hold in reality.A self-consistent framework has only been developed for binaries on circular orbits and studies including eccentric binaries neglect the backreaction from the binary to the spike (e.g., Becker et al. 2022;Dai et al. 2022).The feedback can not only affect the evolution of the semi-major axis but also the eccentricity which has a major impact on inspiral times due to GW emission.The rate of circularization of the IMRI has also recently been suggested as a signature of the spike (Becker et al. 2022) which would require a self-consistent framework that takes into account the effects of the feedback from the binary onto the spike and vice versa to be determined accurately.
Additionally, recent work has shown that the analytic Chandrasekhar approximation may lead to inconsistent evolution of semimajor axis and eccentricity due to the lack of inclusion of drag force from fast moving particles, which can be significant in certain cases (Dosopoulou 2023).As it is difficult to create a self-consistent analytic/semi-analytic framework to include these effects, we have to resort to -body simulations.
-body simulations, although expensive, are considered to be the gold standard of dynamical modelling.Previous approaches have rarely relied on -body simulations as they are extremely computationally expensive.Kavanagh et al. (2020) report that a simulation of 150 orbits of a 1  ⊙ -100  ⊙ binary takes about 3 days using the body code GADGET-2 (Springel 2005).This is very computationally expensive to study long term effects of the IMRI on the spike over tens of thousands of orbits.Therefore, in order to study the secular effects of the IMRI on the spike, different strategies are required.
In our study we describe a novel -body code that is over 100 times faster than traditional -body codes for simulating IMRIs embedded in DM spikes.Using our code, we present self-consistent results from eccentric IMRIs embdedded in DM spikes.
In addition, we also study the effect of rotation in the DM spike which has not been done before.We argue that the inclusion of rotation is important since conservation of angular momentum would dictate that DM spikes should rotate upon formation from rotating galactic halos.Spinning IMBHs can also transfer angular momentum to the surrounding spike, making it rotate (Ferrer et al. 2017).We systematically study the effects of the different post-Newtonian (PN) terms including precession and radiation reaction terms upto 2.5.
We begin by describing our computational methods in Section 2, followed by the models of the IMRIs embedded in DM spikes in Section 3. The results are then described in Section 4, followed by discussions of our results and assumptions and conclusions in Sections 5 and 6 respectively.

COMPUTATIONAL METHODS
The numerical evolution of the IMRI using -body methods is an extremely computationally challenging problem owing to the enormous density of the spike near the central IMBH and the O ( 2 ) pairwise force calculations.An accurate evolution requires fine tuned timesteps for the particles near the central IMBH (primary) and inspiraling object (secondary) to resolve the scattering effects accurately.Previous studies have all relied on analytic/semi-analytic methods for long term simulations of the inspiraling IMRI.To the best of our knowledge, only Kavanagh et al. (2020) performed -body simulations of the binary embedded in a spike.However, the simulations were stopped after a period of ∼ 100 orbits.To understand realistic effects of the spike on the IMRI (and vice-versa) and to calibrate semi-analytic methods, we need long term -body simulations.Here we present a novel -body code1 that is specifically tuned for simulations of IMRIs in DM spikes.
A close analysis of the -body simulation reveals that the bulk of computational time is spent calculating the forces of DM particles near the primary.Since the mass of DM in the region of interest is very small compared to the mass of the primary and the secondary, the self-gravity of the spike can be neglected.This allows us to safely neglect the interactions between the DM particles themselves.Full force calculations are only required for the primary and the secondary.Neglecting the self interaction of the spike effectively reduces the number of force calculations to O (), speeding up the simulations massively.The mean relative force accuracy between a full O ( 2 ) calculation and our approximate method is ≤ 10 −5 in the region of interest as shown in Figure 2, on par with the force accuracies obtained by current Barnes-Hut tree (Barnes & Hut 1986) or Fast Multipole Method (FMM) (Greengard & Rokhlin 1987;Cheng et al. 1999;Dehnen 2002;Zhu 2021) based -body codes.To verify the accuracy of our method, we compare the evolution of a 100 ⊙ −1 ⊙ binary in a DM spike with the FMM based code Taichi (Zhu 2021;Mukherjee et al. 2021Mukherjee et al. , 2023) ) and this method and find negligible differences after ∼ 1500 orbits.
To evolve the particles, we use a 2nd order hierarchical Hamiltonian splitting based integration scheme HOLD in a Drift-Kick-Drift (DKD) fashion with symmetrized timesteps (Pelupessy et al. 2012).The usage of symmetrized timesteps ensures that there is no secular drift in energy and renders better energy conservation than other timestepping schemes (e.g., Makino et al. 2006;Pelupessy et al. 2012;Mukherjee et al. 2021).For more information on the integration scheme, we refer the reader to Pelupessy et al. (2012).The timesteps are controlled by a timestepping parameter  which is proportional to the analytically calculated timestep.In our simulations, we set  = 0.025.This results in a relative energy error of ≤ 10 −10 per orbit, sufficient to follow the dynamics over few hundreds of thousands orbits.
To account for relativistic effects, we add PN terms to the equations of motions for the primary and secondary.The PN equations of motion can be added to the standard Newtonian equation of motion as follows: where a   represents the  th PN term.We include up to 2.5 terms in our calculations following equation 203 from Blanchet (2014).While the 1 and 2 terms only lead to precession of the binary's orbit, the 2.5 term leads to GW radiation resulting We find that in the region of interest, the force accuracy is ≤ 10 −5 , which is comparable to the the force accuracy obtained in tree/FMM based codes.
in shrinkage of the orbit.The DM particles experience the Newtonian potential of the primary and the secondary.In principle this is not fully self-consistent since the DM particles will also experience PN effects in interactions with other DM particles.One can use the Einstein-Infield-Hoffman equations (e.g., Will 2014;Portegies Zwart et al. 2022) to self-consistently simulate the PN evolution of all the particles but this is beyond the scope of this study.
A problem arises when adding velocity dependent forces like PN precession and radiation reaction terms since a leapfrog-like explicit splitting of the Hamiltonian is not achievable with velocity dependent forces.In such a scenario, the simple DKD integration scheme cannot be utilized.An implicit method suggested by Mikkola & Aarseth (2002) is popular but is quite inefficient.The implicit scheme requires multiple iterations to solve for the force.For highly non-linear vector fields, it proves to be quite computationally expensive.Hellström & Mikkola (2010) show that using a clever mathematical trick one can create an explicit DKD-like scheme by extending the phase space of variables by introducing an auxiliary velocity w.An updated integration between step  and  + 1 using the auxiliary velocity w can now be written as: (5) where x is the position, v is the velocity, f is the force and ℎ is the timestep.We note that w 0 = v 0 .
To test the integration scheme, we simulate the evolution of a 1000 ⊙ − 1 ⊙ binary in vacuum and compare the evolution of the eccentricity and semi-major axis to those derived using Peters' equations (Peters 1964) and find that the results are consistent with one another.We also compare the evolution of the binary against that from the publicly available regularization code SpaceHub (Wang et al. 2021) which also includes PN terms upto 2.5PN and find no differences.
For three-body simulations, we utilize the 15th order Gauss-Radau integrator IAS15 (Rein & Spiegel 2015) from the rebound package (Rein & Liu 2012).IAS15 can handle close encounters and is extremely accurate allowing for a closer and in-depth analysis into the dynamics of the scattering process of the binary with a DM particle.The energy is conserved to machine precision.
All simulations are performed on the Vera computing cluster utilizing AMD Epyc 7742 nodes.The majority of the full -body simulations take about 100-120 core hours to finish to completion using a single core.The interactions between the secondary and DM particles are set to have zero softening while we set a relatively conservative softening value of 10 −10 pc between the primary and DM particles in our simulations.This is equal to the Schwartzchild radius for the 10 3  ⊙ IMBH and one-tenth that for the 10 4  ⊙ IMBH.For the three-body simulations, all interactions are unsoftened.A set of simulations were run with softening between the DM particles and the secondary and no major differences were noticed between the simulations that used softening and those that did not.A small value of softening was necessary to prevent some particles near the primary from taking extremely small timesteps.The interaction between the primary and the secondary object is not softened.

MODELS
We are interested in scenarios wherein the IMRI is visible in the LISA/DECIGO band for a duration of ∼ 5 years.To understand the effects of varying eccentricity, primary and secondary masses, and density, we generate physically motivated models corresponding to IMRIs in vacuum where the GW signal is in the LISA/DECIGO band for ∼ 5 years and above the LISA/DECIGO signal to noise threshold.We define the mass-ratio of the secondary to the primary  as  ≡  2  1 where  1 is the mass of the primary and  2 is the mass of the secondary.We choose  1 = 10 3  ⊙ , 10 4  ⊙ and  2 = 1 ⊙ , 10 ⊙ , representing three different scenarios with mass ratios  = 10 −2 , 10 −3 , 10 −4 as depicted in Figure 3.  2 = 1 ⊙ is representative of a scenario the inspiraling object is similar to a white dwarf of a neutron star whereas  2 = 10 ⊙ represents a scenario where the inspiraling object is a stellar mass black hole.The initial properties of the binary are characterized by its initial semi-major axis  0 and eccentricity  0 .We set  0 = 2 × 10 −8 pc for  = 10 −2 , 10 −3 models and  0 = 5 × 10 −8 for the  = 10 −4 model.We also set  0 = 0.4 for  = 10 −2 model,  0 = 0.75 for  = 10 −3 model, and  0 = 0.65 for  = 10 −4 model.Figure 3 provides a representation of the strain versus frequency curves of the models superimposed over the strain-frequency curves for LISA and DECIGO.

Non-rotating models
To generate the -body realizations of the DM spike, we utilize the galactic modeling toolkit Agama (Vasiliev 2019).Following previous studies (e.g., Eda et al. 2015;Kavanagh et al. 2020), we use equation 1 as our density profile.This profile is valid for all for all radii  >  ISCO where  ISCO is the innermost stable circular orbit (ISCO) for the central IMBH.The density profile is taken to be 0 where  <  ISCO .Furthermore,  sp is not considered to be a free parameter but is, instead, calculated as For all our simulations we set  sp = 226 ⊙ /pc 3 following Kavanagh et al. (2020).Using equation 8, we find  sp ≈ 0.54 pc when  1 = 10 3  ⊙ and  sp ≈ 1.17 pc when  1 = 10 4  ⊙ .Additionally, for most of our simulations we set  sp = 7/3.Such a spike profile can arise out of adiabatic growth of an IMBH in an NFW halo (Eda et al. 2015).However, this is not universal and  sp depends on how the DM spike originates.For example, primordial BHs have been shown to produce a  sp = 9/4 spike around them (Boudaud et al. 2021).
To understand how the dephasing changes as  sp changes, we run a set of simulations with  sp = 9/4.Furthermore, to understand how the dephasing changes at a fixed density profile for different , we perform a set of simulations with  1 = 10 4  ⊙ ,  sp = 0.54pc and  sp = 7/3.However, it is hard to robustly predict the existence of such dense spikes.In realistic environments where effects of stars cannot be neglected, DM "crests" are more likely.Such "crests" are formed due to the scattering of DM particles by the more massive more massive stellar particles and have been shown to produce a  sp = 1.5 profile (Merritt et al. 2007a).To understand if such profiles produce any observable dephasing effects, we perform a set of simulations with  sp = 1.5,  1 = 10 3 , 10 4  ⊙ and  2 = 1 ⊙ .The density spikes considered in our simulations have been plotted in Figure 4.
To generate models with finite total mass, we use an exponential truncation function with a truncation radius of 10 −5 pc for  = 10 −2 , 10 −3 models and 10 −6 pc for the  = 10 −4 model.We verify that our choice of truncation radius does not affect the mass profile of the spike in the region of interest.We also compare our initial -body profiles to those from Kavanagh et al. (2020) and find that the density profiles in the region of interest and the velocity profiles match.Additionally, we ensure that the density profile remains stable for the duration of the simulations.The mass of the DM particle,  DM , is chosen carefully.Since  1 ≫  2 , the dynamics of inspiral is dependent on the mass ratio of the secondary to the mass of the DM particle (e.g., Merritt 2013, chapter 8 and references therein).We define the mass-ratio of the secondary to the DM particle as  DM ≡  DM  2 .Three-body simulations show that using  DM ≤ 10 −2 , the results are convergent.The energy, angular momentum, and ejection time distributions are consistent between all of the models where  DM ≤ 10 −2 .Such a low mass ratio is desirable since it reduces any spurious three-body effects leading to sudden jumps in the evolution of the semi-major axis or eccentricity and ensures that the evolution is smooth.Accordingly, based on resolution and computational expenses, we choose  DM = 5 × 10 −5  ⊙ for the simulations using  sp = 7/3 and 10 −5  ⊙ when  sp = 9/4.This results in ≈ 8k − 10k particles in the simulation.Extra care is taken to ensure that the region surrounding the secondary have a sufficient number of particles to resolve the scattering process.Since the scattering process is quite stochastic, results from an individual -body simulation are unreliable.Accordingly, five independent realizations of every model are generated and simulated.
For most of the simulations, only the 2.5 term is added to the equation of motion of the binary.This is done in order to provide a direct comparison to previous studies such as Kavanagh et al. (2020) and Becker et al. (2022) where the effects of relativistic precession are neglected.In general, we do not expect the dephasing to be affected by any precession effects.However, the net precession can itself be affected by the Newtonian precession of the binary in the potential of the spike in addition to the relativistic precession, which can also be a The characteristic strain of the IMRIs used as our initial conditions as a function of the frequency  in vacuum for different mass-ratio  models.The central IMBH has a mass of  1 while the inspiraling object has a mass of  2 .The binary is defined by its initial semi-major axis  and eccentricity .The luminosity distance is denoted as   .All of the IMRIs presented here have a merger time of ∼ 5 years.Since eccentric binaries radiate in multiple harmonics, we have plotted the strain from the second and third harmonics.In all of the scenarios, we find that the second harmonic has an equivalent or higher strain than the third harmonic.Even higher harmonics have lower strains and are harder to detect using LISA but would be detectable using DECIGO.1.0 10 −4 5 × 10 −8 , 0.65 7/3, 1.17 5 × 10 −5 2.5PN 2.5 yr 5 Table 1.A summary of the initial conditions for our non-rotating models.The rotating models have similar initial conditions but are either in a net prograde motion or retrograde motion with respect to the motion of the IMRI.For more information on how to generate the initial conditions, we refer the reader to section 3.  signature of the spike (Dai et al. 2022).It can, also, potentially affect the exchange of angular momentum between the spike to the IMRI.In order to understand if relativistic precession affects our results, we run a set of simulations with all PN terms enabled up to 2.5 for  = 10 −3 models.We leave a more systematic study of precession effects to future work.A summary of the initial conditions for the non-rotating models can be found in Table 1.
We evolve the  = 10 −2 , and 10 −3 models to merger.This is de-fined as the time when  <  ISCO .The simulations where  = 10 −4 are extremely computationally expensive.We only achieve an evolution to 2.5 yr in 15-20 days of computing time.This results in about 2 × 10 6 orbits of the secondary around the primary.We calculate the dephasing cycles by extrapolating the rest of the evolution in vacuum to save on computational resources.To verify if the extrapolation produces valid results, we run a set of lower resolution simulations with 1024 DM particles to completion and compare our extrapolated results to the results from our lower resolution simulations.The extrapolated results can be seen as a lower limit on the number of dephasing cycles.We compare our results against the vacuum evolution computed using the Peters (1964) analytic formula and also against those calculated using the Chandrasekhar DF assuming there is no backreaction on to the spike from the IMRI.This is done using the publicly available code IMRIPy (Becker et al. 2022;Becker & Sagunski 2023;Becker 2023).The IMRIPy results are denoted as "Static DF" in the results section.Unfortunately, we could not use HaloFeedback (Kavanagh et al. 2020) for our purposes since it can only simulate circular binaries.

Rotating models
As mentioned in the introduction, we aim to understand the effect of rotation in the DM spike and its effect on the IMRI.If the DM halo surrounding the IMBH is rotating, angular momentum conservation will dictate that DM spike should rotate as well.Additionally, spinning IMBHs will transfer angular momentum from the IMBH to the spike, torquing the spike up or down depending on the direction of the primitive rotation of the spike (Ferrer et al. 2017).In a rotating spike, exchange of angular momentum between the IMRI and the spike will be enhanced compared to the non-rotating models and the eccentricity of the IMRI can be significantly affected.We aim to understand this effect in our study.To the best of our knowledge, there is no literature available quantifying the level of rotation in the spike and its correlation to the galactic halo rotation or the spin of the central IMBH.Therefore, we caution the reader that our rotating models are somewhat ad hoc and therefore exploratory.However, they are still useful in understanding the qualitative effects a spike that is in a counter-rotating motion (retrograde motion) or co-rotating motion (prograde motion) with respect to the binary.
To include rotation, we follow the Lynden-Bell trick (Lynden-Bell 1960).This has been motivated by several studies which use the method to generate rotating models of galactic nuclei (e.g., Holley-Bockelmann & Khan 2015; Khan et al. 2020;Khan & Holley-Bockelmann 2021;Varisco et al. 2021) to examine the effects of rotation upon the dynamics of MBH binaries.In this method, we flip the -component of the angular momentum of the DM particles,  ,  , to generate prograde or retrograde models.We note that the -component of the angular momentum of the IMRI,   , is always positive.Accordingly, to generate co-rotating or prograde models, we reverse the direction of the  and  components of the velocity for all DM particles which satisfy  ,  < 0. To generate counterrotating or retrograde models, we do the same, except with particles with  ,  > 0. This leaves the density and velocity profiles of the spike unchanged.In principle, this is not a fully self-consistent method of generating rotating models.Rotation is expected to flatten out the models to some extent so it is somewhat ad hoc to include rotation in spherical spikes.One can construct a distribution function  (,   ) as a function of the energy  and -component of the angular momentum   as is done in Wang et al. (2014) to generate self-consistent rotating models which are flattened.However, this is beyond the scope of this work.Irrespective of that, our models are able to qualitatively understand the effect of rotation on the IMRI.In future studies, we plan on investigating self-consistent models of rotating DM spikes to understand how the geometry of the spike affects the dynamics.
We denote the ratio of retrograde particles to the total number of particles as F .In our retrograde simulations, F = 1, indicating that all particles are moving in a retrograde motion with respect to the binary, and in our prograde simulation, F = 0, which is opposite of the previous scenario.Realistic spikes are expected to have 0 ≤ F ≤ 1.Since the simulations are quite expensive even with our fast -body code, we restrict our study only to the F = 0, 1 cases.This helps us put limits on the dephasing effect and compare the rotating models to fully isotropic non-rotating models.We note that the Lynden-Bell method is not the same as introducing rigid body rotation into the DM spike.The latter generates much stronger rotation whereas our method introduces a weaker level of rotation.
We generate prograde and retrograde rotating versions of all of the models presented in Table 1 and create five independent realizations for each configuration.All models were evolved to merger, except in the case where  1 = 10 4  ⊙ as explained in the section above.

RESULTS
Owing to the inherent stochasticity present due to the discrete nature of -body simulations, we note that the results presented here are calculated after taking the average of the quantities across all five independent realizations.In all of the calculated values, we noticed a maximum of 1 percent difference between results from individual simulations This does not affect our overall results and indicates the robustness of our simulations.We first present the results from the non-rotating models before moving on to the rotating models.

Non-rotating models
The mean orbital frequency of the binary  can be written as where  =  1 +  2 .We calculate the amount of dephasing by taking the difference between the number of GW cycles completed by the binary with and without the spike.Following Becker et al. (2022), we write the number of GW cycles  () ( 0 ,  final ) for the  th harmonic between some initial time  0 and final time  final as follows: We can, then, calculate the difference in the number of GW cycles Δ () ( 0 ,  final ), or dephasing, as: In our calculations,  final is taken to be the time of merger For comparison amongst different models, we choose  = 2, as has been the strategy in previous studies.Eccentric binaries radiate in multiple harmonics and the dephasing will be larger for larger .However, these higher harmonics are harder to detect.We present a comparison of the evolution of the mean orbital frequency of the binary  , the absolute value of number of dephasing cycles of the 2nd harmonic | (2) |, and the evolution of eccentricity of the binary  as a function of its semi-major axis  using -body models and IMRIPy for the non-rotating models in  sp = 7/3 spike with  sp calculated using equation 8 in Figure 5. Examining the  = 10 −2 model, our -body simulations predict that the binary merges 19 days earlier if it is embedded in a DM spike than in vacuum.It takes about ≈ 10 4 fewer GW cycles than in vacuum to merge.Comparing the -body model to the DF model with a static spike, we find a 100× reduction in the number of dephasing.The binary disrupts the spike almost completely within the first 0.1 years of inspiral leading to a drastic reduction in the amount of dephasing.This is unsurprising as the mass ratio of the binary is large which injects a substantial amount of energy ejecting DM particles around the secondary leading to the disruption of the spike.As such, the circularization rate of the binary is similar to that in vacuum as is evident from the  −  plot.
As we lower the mass ratio to  = 10 −3 , the binary merges earlier.It takes about 141 fewer days to merge in the  sp = 7/3 -body model than its vacuum counterpart.This creates a larger amount of dephasing.In a  sp = 7/3 spike, the dephasing with respect to vacuum is ≈ 2 × 10 5 , only 3× lower than what is predicted by the DF with static spike model.Although not directly comparable, we highlight the fact that this is much larger than the 100× reduction for the same mass ratio found by Kavanagh et al. (2020) in the semianalytic HaloFeedback models but with the secondary on circular orbit.Since the self gravity of the spike is minimal, we should expect a ∼ 3× reduction in the amount of dephasing compared to the static DF models irrespective of the density profile.To verify that, we compare the dephasing cycles in the  = 9/4 model with  = 10 −3 in Figure 6.We find that |Δ (2) | ≈ 1.8×10 5 for the static DF models whereas |Δ (2) | ≈ 5.2 × 10 4 , consistent with the ∼ 3× reduction we expected.The rate of circularization is quite similar between the A comparison of the binary parameters for different mass ratio models evolving in a  sp = 7/3 spike with  sp calculated using equation 8.The evolution in vacuum (purple line) is calculated using Peters (1964) analytic formula, while the evolution calculated using the Chandrasekhar DF formula assuming a static spike (orange line) is calculated using IMRIPy.They are compared to the evolution from our  -body simulations (green line).Left column: the mean orbital frequency  as a function of time  in years.Middle column: the estimated number of dephasing cycles of the second harmonic |Δ (2) | as a function of the frequency  in Hz.Right column: the eccentricity  as a function of the semi-major axis .We notice that in higher mass ratio models, the evolution of the binary is similar to that in vacuum.For  = 10 −2 model, there is a 100× reduction in the estimated number of dephasing cycles compared to the evolution calculated using Chandrasekhar DF formula as the spike has been disrupted in a very short time span.As we decrease the mass ratio, the disruption decreases.We notice that in case of the  = 10 −3 model, the dephasing is only reduced by 3× compared to the evolution calculated using IMRIPy.For  = 10 −4 model, we find that dephasing is a factor of 3 larger than what we obtain using the Chandrasekhar formula, with little to no disruption of the spike.This signals that the Chandrasekhar DF might be insufficient to explain the evolution of the binary in DM spikes.vacuum, static DF, and the -body models.Nevertheless, we find that the eccentricity of the binary as a function of its semi-major axis for the -body models lies approximately in between the static DF and the vacuum cases.
For the  sp = 1.5,  = 10 −3 case representing perturbed DM "crests", our findings are less optimistic.We note that such crests can also be formed due to adiabatic growth in an off-center IMBH (Zhao & Silk 2005).In Figure 7 we plot the evolution of the binary frequency  as a function of time for this model.We find that the frequency evolution in the -body models practically overlaps with that from the vacuum model.Unsurprisingly, the calculated dephasing is ≤ O (10) which is larger than the limits of uncertainty from our -body simulations.Such a small dephasing would also be much harder to detect than those from the denser spike models.Since the formation of the DM crests are more robust than that of DM spikes, our findings pose an important question of whether the spike models are overly optimistic.This is discussed further in the section 5.3.
We examine the evolution of the density profile in the  = 10 −3 ,  sp = 7/3 model over time due to the effect of the binary in Figure 8.The binary injects energy into the spike by scattering the DM  2) | as a function of the binary frequency  in Hz for a  = 10 −3 binary embedded in a  sp = 9/4 spike.The color scheme is the same as that used in Figure 5.We notice that similar to the  sp = 7/3 case, the amount of dephasing in our  -body models is reduced by ∼ 3× compared to the DF models. sp = 1.5 model.We find that the low density of the spike results in no discernable differences between the inspiral in vacuum and that in the spike.This results in ≤ O (10) cycles which might not be detectable.
particles, thereby reducing the density of the spike.It preferentially interacts and scatters particles near itself, so we expect a drop in the density profile near the binary over time.The local efficiency of scattering will depend on the orbital parameters of the binary and density profile near the secondary.During a 2 year time span, we notice that the binary has drastic effects on the density of the spike.There is up to a factor of 100 reduction in the density of DM particles near the initial semi-major axis of  = 2×10 −8 .The density of DM particles reduces from ∼ 10 20  ⊙ pc −3 to ∼ 10 18  ⊙ pc −3 .The binary effectively carves out a flat core near its semi-major axis.Since the secondary spends a larger time near its apoapsis than its periapsis, we expect the scattering to be most effective near the apoapsis carving out a core.The apoapsis of the secondary is 3.75×10 −8 pc.We notice that the flat core extends from ∼ 10 −8 pc to 4 × 10 −8 pc consistent with the hypothesis.The impact further away, and near the periapsis is smaller.This has important implications for the survival and detectability of DM spikes and is discussed at length in section 5.3.Similar findings were noted in Kavanagh et al.We notice that the inspiral of the IMRI leads to a significant disruption in the spike near the semi-major axis  of the IMRI.The density of DM in that region is reduced by a factor of 100 compared to the original density after an inspiral time of 2 years.This is qualitatively consistent with the findings from Kavanagh et al. (2020).However, we find that the disruption to the spike due to the IMRI in our  -body is much slower than what is found by Kavanagh et al. (2020) who find that the spike is significantly disrupted within 0.1 years.
(2020) in case of circular binaries.However, Kavanagh et al. (2020) find that the disruption happens quite quickly, effectively within 0.1 yr from the beginning of the simulation for  = 10 −3 , whereas in our case the disruption is much slower, by almost 20×.This leads to a much larger dephasing effect compared to the HaloFeedback models.Notably, the evolution of the density profile does not show the presence of a wake or a density "bump" behind the secondary, which is present in the self-consistent HaloFeedback models.According to the authors, the "bump" is caused due to scattering of DM particles near the secondary which give rise to the DF effect.The absence of the "bump" and the longer spike disruption time suggests that DF theory is inconsistent with the results from our -body simulations.Additionally, we find that the effect of softening is minimal on the feedback time and effects on the density profile.When a softening of ≈ 10 −10 pc is used between the secondary and DM particles, we find that the feedback at  < 10 −8 pc is decreased but the differences between the density profiles in that region between non-softened and the softened sims never vary substantially.
For our  = 10 −4 non-rotating model with  sp = 1.17 pc, examining Figure 5, we notice a very surprising result.While the binary in the static DF model merges 60.3 days faster then the vacuum case, the binary in the -body model merges 169 days faster.This leads to dephasing effects that are almost 3× larger than the DF model.Whereas in the static DF case we find that |Δ (2) | ≈ 6 × 10 4 , our -body simulation suggests that |Δ (2) | ≈ 1.5 × 10 5 .Since the results are extrapolated beyond 2.5 years using the Peters (1964) equation, we compared them to those obtained from set of lower resolution simulations which are run to completion.We find that the extrapolated results are consistent with those from the lower resolution simulations.The extrapolated results fall within a standard deviation of the results from the lower resolution simulations.For the models where we fix  sp = 0.54 pc, similar to the  = 10 −2 and 10 −3 models, we notice a decrease in the amount of dephasing   6 but for  = 10 −4 binary embedded in a  sp = 7/3,  sp = 0.54pc spike.We notice that the net dephasing predicted by the body simulation is comparable to that predicted by the IMRIPy simulation, although at higher frequencies the dephasing falls off faster in the  -body models.The density in the region of interest for this particular model is about 6 times lower than that for the  = 10 −4 ,  sp = 1.17pc model and we notice a similar decrease in the amount of dephasing in this scenario compared to that model.compared to the  sp = 1.17 pc model.We present the dephasing in  = 10 −4 ,  sp = 0.54 pc model in Figure 9 and find that the -body and the IMRIPy simulations predict a similar amount of dephasing, with |Δ (2) | = 3 × 10 4 .This is ≈ 6× lower than that obtained in the  sp = 1.17 pc models which is unsurprising given the lower density in the  sp = 0.54 pc models.The density in the region of interest is 6.1× lower in the model with  sp = 0.54 pc compared to the model with  sp = 1.17 pc, which explains the factor of 6 reduction in the dephasing.This, along with the results from the  = 10 −3 models, strongly suggests that the dephasing is proportional to the local density of the profile at the position of the secondary.Similar to our findings from the  = 10 −3 model, we find the spike plays a very little role in dephasing if  sp = 1.5 even when  1 = 10 4  ⊙ .The dephasing is again ≤ O (10), and below our uncertainty thresholds.
The results obtained for the  sp = 1.17 pc case for the  = 10 −4 binary is a factor of 15 greater than what was found by Kavanagh et al. (2020) in their HaloFeedback simulations for similar models on circular orbits.We also compare the spike density profiles at the beginning and end of our simulation and find that there is no degradation in the density profile of the spike, while the reduction of the DM density is observed in the HaloFeedback models even with such low mass ratios.This, along with our previous results, suggests that previous analytic/semi-analytic calculations using DF underestimated the total number of dephasing cycles by a factor of one to ten (or even larger) for lower-mass ratio binaries.This is especially evident while comparing the -body simulations to the HaloFeedback simulations, where we find a difference of ten to hundred times.We use three-body simulations in section 4.3 to show that three-body scattering provides a better description of the dynamics of the binary in the spike and could explain the inconsistency observed in this section, while showing that the impact of DF is subdominant.

Rotating models
We now move on to the analysis of the binary inspiral in the rotating models.From DF theory, we would expect the inspiral of the binary in counter-rotating or retrograde models to be slower than that in co-rotating or prograde models.This is caused due to the fact that in the prograde models, the relative velocity of the DM particles is smaller than that in the retrograde models.Since the force of DF is inversely proportional to the relative velocity of the interacting particles, we would expect the DF force to be larger in the prograde scenario, resulting in a faster inspiral.
Examining the properties of the binary over time in both the prograde models and retrograde models in Figure 10 we notice a very surprising result.We find that in all of the cases, the retrograde models merge faster than the prograde and even the non-rotating models.Even more astonishing is the fact that the prograde models are slower to inspiral than the vacuum models.This is the opposite of what we would have expected from DF. Examining the  = 10 −2 models, we find that the expected amount of dephasing increases by almost 4× compared to its non-rotating counterpart.Since we are plotting the abolute value of the number of dephasing cycles with respect to vacuum, we find that the prograde model has a larger amount of dephasing, almost 1.5× that of the non-rotating model.However, we want to stress that the prograde model actually take longer to inspiral than the vacuum model.So it actually takes 1.5 × 10 4 more GW cycles than the vaccum model to merge.
As we lower the mass ratio, we notice the same trend.In the  = 10 −3 models, we find that |Δ (2) | ≈ 7 × 10 5 in the retrograde model, whereas |Δ (2) | ≈ 2 × 10 5 in the non-rotating model, almost 3.5× lower.Similar to the 10 −2 model, we find that the prograde model takes 4.5 × 10 5 more GW cycles to merge than the vacuum model.Interestingly, we also notice that the retrograde model eccentrifies quickly in the beginning, before circularizing later due to GW effects.The opposite is observed in prograde models.There is an accelerated circularization in the beginning, followed by a period of circularization led by the emission of GW.
As we lower the mass-ratio to  = 10 −4 , we notice that the eccentrification and circularization effects increase.Examining the eccentricity as a function of the semi-major axis, we find that the evolution in the retrograde and prograde models diverge quickly.The binary eccentrifies quickly, reaching a maximum eccentricity of 0.67 before circularizing.The opposite happens in the prograde models where the binary circularizes quickly, within the first one year time span, reaching an eccentricity of 0.61.This eccentrification and circularization lead to accelerated or deccelerated inspiral due to GW emission since the GW emission is very sensitive to the binary eccentricity.This leads to the larger dephasing we see in retrograde models.We find that the retrograde model merges about 550 days earlier compared to the vacuum model, whereas the prograde model takes 484 days longer.In the retrograde models, this produces a dephasing effect that is 2.5 times as large as the non-rotating models, with the binary taking 5 × 10 5 fewer cycles to merge compared to the vacuum models.The opposite effect is noted in the prograde models with the binary taking 3.5 × 10 5 more cycles to merge.Interestingly, in all of our simulations we find that the ratio of dephasing cycles in retrograde rotating to non-rotating models is between 2.5 − 3. Similarly for prograde models, we find that the ratio of dephasing cycles in rotating to non-rotating models is about 1.5 − 2, independent of the density profile and mass ratio.This is likely to be an artifact of our chosen initial conditions rather than a fundamental property of rotation itself.
We also compare the density profiles of the non-rotating to the rotating models over time to understand if the feedback from the binary to the spike changes upon the inclusion of rotation.Surprisingly, we find that the density profiles are consistent among the rotating and non-rotating models.Although rotation changes the dynamics of the binary significantly, the effect of binary on the spike is consistent among non-rotating and rotating models.

Three-body simulations
The results obtained in the previous section suggest that DF theory is insufficient at explaining the long term dynamics of the binary.In fact, this is not a very surprising result as it has been known that upon the formation of a hard binary, three-body scattering becomes more effective than DF at dissipating energy, similar to SMBH binaries (e.g., Begelman et al. 1980;Merritt 2013).A binary is said to become a hard binary when the separation between the primary and secondary falls below the hard binary radius  h which is defined as where  infl is the influence radius of the primary and is defined as the radius that encloses twice the mass of the primary (e.g., Mukherjee et al. 2023).For our models, we find that  h ranges from 6.75×10 −3 pc in the  = 10 −2 model to 1.46 × 10 −4 pc in the  = 10 −4 model.Therefore, all of our binaries are in the hard binary limit at the beginning of our simulations.
In the three-body scattering scenario, the binary undergoes a complicated three-body interaction in which the incoming particle interacts with the binary multiple times.This results in the ejection of the particle eventually, leading to shrinkage of the binary's orbit.This process has been called the slingshot mechanism and is fundamentally different from Chandrasekhar DF which is a cumulative effect of two-body encounters.
To understand which of the two above-mentioned processes dominate, we can analyze the -body simulations and compare the efficiency of three-body scattering to that of dynamical friction.The DM particles are orbiting the primary and are bound to it initially.This allows us to obtain the semi-major axis  DM and eccentricity  DM of each particle.From our  = 10 −3 model with  sp = 7/3, we select a subset of DM particles which satisfy 0.5 0 ≤  DM ≤ 1.5 0 , where  0 is the initial semi-major axis of the binary.In principle, all particles interact with the binary, but particles with semi-major axis close to the binary's interact strongly.We, then, calculate the initial and final energies of the selected particles at the beginning and end of our simulations to calculate the change in energy of the DM particle (Δ DM ).We find that all of the selected particles are ejected from the spike by the end of the simulation.
The selected particles are evolved from their initial positions along with the binary individually using IAS15 for the same duration as the full -body simulations.The energy of each particle is recorded at the beginning and end of the simulation.The difference in energies provides an estimate of the change in energy due to three-body effects.On the other hand, calculating the energy dissipation due to dynamical friction from the -body simulations is somewhat nontrivial.According to Chandrasekhar (1943), as  2 moves through the medium of DM particles, it experiences a number of two-body encounters that change the velocity of the secondary in a direction parallel and opposite to the initial velocity of the secondary.Each particle contributes a net change in the velocity of the secondary giving rise to dynamical friction force over time.We can estimate the dynamical friction force due to each individual particle using the method described in Ma et al. (2023).The energy change due to dynamical friction attributed to the  th DM particle  i,df can be written as where a i,df can be calculated using equation 9 from Ma et al. (2023) and v 2 is the velocity of the secondary.We calculate  i,df using the saved snapshots and integrate over time to find the net energy change due to dynamical friction.The results from our simulations are stored at a fine enough time resolution that this method is able to provide a good approximation of the dynamical friction energy loss.We present the relative distribution of the energy change of the selected particles from the -body simulation, and compare it to the energy change due to three-body effects and dynamical friction in Figure 11.For clarity, we present the normalized change in energy Δ Ẽ where we normalize the change in energy with respect to the binding energy of the DM particle with  DM =  0 .We find that the distribution of the change in energy from the -body simulations matches that from the three-body simulations but is inconsistent with the energy change due to dynamical friction.The -body and threebody simulations show that each selected particle experiences, on average, a normalized energy change of 1.0-1.5.In the dynamical friction scenario, the average energy loss per particle is about 10× lower.Additionally, we find that, akin to the -body simulations, all particles are ejected from the system in the three-body simulations.From the estimates of the dynamical friction force, we find that only 10-20 per cent of the particles would be ejected from the system entirely, consistent with the findings from Kavanagh et al. (2020), but is in tension with the findings from the -body simulation.This confirms our hypothesis that three-body scattering, and not dynamical friction is the predominant method of energy dissipation in our simulations.Although not presented here, we find a similar story across all of the models used in this study.Interestingly, as the mass ratio is decreased, the fraction of energy loss due to the three-body scattering increases.This is consistent with Merritt (2013) where the three-body scattering efficiency is proportional to  −1 .We find a similar distribution in the change in energy of the particles from both the three-body and full -body simulations in our  = 10 −4 models as well.Our results also highlight that our -body simulations are robust, accurate, and consistent with results from the extremely accurate IAS15 integrator where the net energy is conserved to machine precision.
What causes the counter intuitive results from our rotating simulations?Merritt et al. (2009), Iwasawa et al. (2011), andSesana et al. (2011) hold clues that are able to shed some light on this mystery.During a three-body encounter, the binary exchanges energy and angular momentum with the particle in a complicated fashion.Iwasawa et al. (2011), in mergers of MBH binaries in galactic nuclei, noted that counter-rotating stars are much more effective in extracting angular momentum from the binary during the three-body scattering phase.As explained by Merritt et al. (2009) and Sesana et al. (2011), this is caused due to the torquing mechanism of the binary's potential on to the particle which leads to a secular evolution of the particle's eccentricity and inclination.This mechanism converts the counterrotating particles to co-rotating which are then preferentially ejected by the binary (Iwasawa et al. 2011).This results in a larger change in angular momentum of the particle and as such, counter-rotating particles are able to extract angular momentum from the binary more efficiently compared to co-rotating particles.This process becomes more efficient for more eccentric binaries.
To understand if the mechanism mentioned above is able to explain the counter-intuitive results from the previous sections and to provide a better description of the dynamics of the binary, we perform threebody simulations.This is done to understand the transfer of energy and angular momentum during the scattering of a DM particle by the binary.We follow the similar steps as in Sesana et al. (2011) to set up our simulation with some differences.A comparison of the binary evolution parameters for different mass ratio models inspiraling in a  sp = 7/3 spike similar to Figure 5 but also including the evolution from rotating models.We notice that the in all of the retrograde models, the binary merges faster than in the prograde and even the non-rotating models.This leads to a major enhancement in the number of dephasing cycles, by as much as 3× that of the non-rotating model in the case of  = 10 −4 .Since we are plotting the absolute value of the number of dephasing cycles, the dephasing in case of the prograde models appears to be positive.However, the prograde models actually merge slower than even the vacuum models and Δ (2) is actually negative (indicated by the dashed lines) .Thus the number of GW cycles in prograde models is larger than that in vacuum.We also find that in the retrograde models, the binary circularization rate is slower than that in the prograde models.In the latter scenario, the circularization is enhanced leading to a slower inspiral.As the mass ratio decreases, the effect of rotation becomes more prominent suggesting that transfer of angular momentum between the spike and the IMRI contributes majorly and cannot be ignored.
We first generate a  = 1 particle realization using Agama containg a primary with mass  1 = 10 3  ⊙ and  sp = 7/3 and then use the Lynden-Bell trick to generate prograde and retrograde models.We then place the secondary with mass  2 = 10 ⊙ at  0 = 2 × 10 −8 pc with varying eccentricity.We ensure that the binary lies in the  −  plane.From the generated initial model, we only select DM particles which have  DM ≈  0 .This results in ∼ 2500 DM particles being selected.We, then, evolve each particle with the binary individually, running 2500 three-body simulations using IAS15, until the DM particle is completely ejected from the system.We record the initial and final values of the energy and angular momentum of the particle and the binary.We note that the simulations are purely Newtonian.
The eccentricity of the binary can be written as follows: where  is the binding energy of the binary,  is the angular momentum of the binary,  =  1 +  2 , and  =  1  2  1 + 2 .The last equality follows from the fact that our binary is in the  −  interacting DM particles from our  -body simulations compared to that from three-body scattering simulations and dynamical friction approximation.We find that the distribution of energy change of strongly interacting particles is consistent between the  -body simulations and the three-body simulations but not with the calculated values using the dynamical friction approximation.This confirms our hypothesis that three-body scattering, and not dynamical friction is responsible for dissipating energy from the binary.
provides us with the change in eccentricity Δ, which we find to be: where Δ is the change in , Δ  is the change in   and , called the eccentrification parameter, has been defined using the last equality.We find that for  > 0, the binary becomes more eccentric after the encounter with the DM particle, and for  < 0, the binary becomes more circular.Since the energy and angular momentum in the simulation are conserved to machine precision, we expect that the change in energy and angular momentum of the binary are equal and opposite to the change in energy and angular momentum of the DM particle.Therefore, For the sake of clarity, we use the normalized change in energy Δ Ẽ (as defined previously), change in angular momentum Δ L , and eccentrification parameter χ in our calculations which are defined as follows: where   ,  is the angular momentum of a DM particle on a circular orbit with    =  0 , and    is the binding energy of the DM particle.
To understand how the retrograde and prograde families contribute differently to the change in energy and angular momentum of the binary, we plot the distribution of Δ Ẽ and Δ L for a binary with  = 0.1, 0.7 in Figure 12.Quite surprisingly, we find that the energy exchanged during the scattering process is similar between the retrograde and prograde models with relative differences of at most 1 − 2 percent between them.This is in tension with the DF theory where we would expect the retrograde family to have a much lower final energy than the prograde family owing to the larger relative velocity in the former case.Δ Ẽ is always positive signaling that all particles contribute to dissipating energy from the binary, even fast moving particles with velocities larger than the velocity of the secondary, in contrast with the assumptions made in previous studies (e.g., Kavanagh et al. 2020;Becker et al. 2022).Additionally, we note that the distribution of energy exchanged does not change substantially as the eccentricity is increased.Although not presented here, we find that the mean energy exchanged Δ Ẽmean remains almost constant across all values of , consistent with the findings from Sesana et al. (2011).
Contrary to the distribution of Δ Ẽ, we find major differences between the retrograde and prograde models while comparing the distribution of Δ L .We find that the distribution of angular momentum exchanged always skews towards lower values for prograde models compared to the retrograde models.The mean angular momentum exchanged is almost 37% lower for the prograde family than the retrograde family when  = 0.1.As the eccentricity increases, the differences become even more drastic.At  = 0.7, the mean angular momentum exchanged by the retgrograde family is almost 10× larger than the prograde family.Thus, we find that counter-rotating particles are more efficient at stealing angular momentum from the binary than co-rotating particles, especially for more eccentric binaries.
We calculate the mean values of Δ Ẽ and Δ L across different binary eccentricities to calculate the eccentrification parameter using equation 15.We present the normalized eccentrification parameter χ as a function of  in Figure 13 for both the retrograde and prograde models.Consistent with Sesana et al. (2011), we find that χ is always positive in retograde models while it is always negative in prograde models.Thus, in our retrograde models, the binary always undergoes eccentrification as a result of interaction with the spike, whereas the binary undergoes circularization in the prograde scenario.This is consistent with findings from Fokker-Planck models of binary evolution in rotating nuclei.According to Rasskazov & Merritt (2017) equation 84a, where ⟨Δ⟩ is the mean change in eccentricity,  is the dimensionless binary hardening rate, and  is the dimensionless eccentricity growth rate (Quinlan 1996;Sesana et al. 2006).Using Fokker-Planck methods, Rasskazov & Merritt (2017) find that, in rotating nuclei, where  is the orbital inclination of the binary and  is the fraction of prograde to retrograde particles (as opposed to F which is the ratio of retrograde to prograde particles).Since our binary inclination is 0, cos() = 1.For a retrograde rotation,  = 0 and we find,  > 0 so the binary eccentrifies as a result of interactions with DM particles.In the prograde case,  = 1 and  < 0 and the binary circularizes due to three body interactions.This explains why the binary eccentrifies or circularizes quickly at the beginning of our full -body simulations.
The rate of circularization is lower in the retrograde models than in the non-rotating and the prograde models, resulting in a faster inspiral as the effects of GW radiation are stronger at larger eccentricities.The inspiral is the slowest in prograde models as the binary circularizes faster than in non-rotating and even vacuum scenarios.While not presented here, we verified that our results are valid across the massratios considered in this study.One can also use the three-body simulations to understand how long it takes for the three-body interactions to disrupt the spike significantly.For a fixed mass ratio of  = 10 −2 , looking at Figure 13, we find that for lower eccentricity binaries, the mean ejection time  ej is larger for particles on retrograde orbits than prograde orbits.As the Top: probability distribution of the normalized energy change Δ Ẽ of 2500 DM particles from the three-body simulations performed until the particle is ejected by interactions with the binary for both prograde and retrograde rotation models with two different binary eccentricities .We notice that in both prograde and retrograde models, the distribution of the normalized energy is similar and does not change with the binary eccentricity.Bottom: similar to top but the probability distribution of the normalized angular momentum Δ L of 2500 DM particles after ejection.We notice that there are major differences in the distribution of Δ L between prograde and retrograde models.In particular, in the retrograde models, the angular momentum of the particle is larger than that in the prograde models, especially for more eccentric binaries.This suggests that retrograde particles are more efficient at "stealing" angular momentum from the binary, especially at larger .
eccentricity is increased, the ejection time becomes similar for both prograde and retrograde orbits.Since the particles are preferentially ejected when they are co-rotating with the binary, we postulate that difference between the prograde and the retrograde models across eccentricity is caused due to the fact that a less eccentric binary exterts weaker torques on the particle leading to a longer secular timescale over which the binary converts the particle from retrograde to prograde.For moderately or highly eccentric binaries, the conversion from retrograde to prograde rotation happens quickly, leading to a similar ejection time between the two models.
We can also study the ejection time as a function of the mass ratio.We take a non-rotating version of the model generated above and run multiple simulations with different mass ratio binaries with  0 = 2 × 10 −8 pc,  0 = 0.7 and plot the mean ejection time  ej as a function of the mass ratio  in Figure 14.We find that the relationship between the ejection time and the mass-ratio can be described by a power-law.We find that  ej ∝  −1.5 .The proportionality constant is a function of the central IMBH mass, semi-major axis, and eccentricity of the binary.
As pointed in section 3.2, we caution the reader, that our models are somewhat unphysical as F , the fraction of counter-rotating particles to the total number of particles, is 0 (prograde models) or 1 (retrograde models).The eccentrification or circularization sensitively depends on this fraction F .Sesana et al. (2011) find that when this fraction is greater than 0.5, the binary eccentrifies as a result of three-body interactions.This suggests that in our non-rotating models, the binary undergoes mild eccentrification as a result of interactions with the spike.A systematic study of the change in the rate of circularization as a function of F is left for a future study.
The three-body simulations are easier to parallelize, and faster to run than the -body simulations which can provide the foundation for a more extensive parameter space study in the future.We can Figure 13.Top: the normalized eccentrification parameter χ as a function of the eccentricity  of the binary.Bottom: the mean time of ejection of the DM particle  ej in units of the number of orbits of the binary as function of the eccentricity of the binary.We notice that for retrograde models χ is always positive and always negative for prograde models.This results in an eccentrification of the binary in case of retrograde rotation and circularization of the binary in case of prograde rotation.We also notice that the time of ejection of the DM particle is almost uniform for the prograde case across different values of .For retrograde scenario, for larger values of ,  ej is similar to that in the prograde case.However, it rapidly increases as  decreases.also use the three-body simulations to derive the distribution of the energy and angular momentum changes as a function of the binary semi-major axis and eccentricity which can later be used in semianalytic models like HaloFeedback.Future work might also involve considering the binary's angular momentum to be tilted relative to the angular momentum of the DM spike.

Precession effects
We present a brief analysis of the effects of relativistic precession on the results from both our rotating and non-rotating models.Relativistic precession is included by using the PN1 and PN2 terms in a set of rotating and non rotating  = 10 −3 models in  sp = 7/3 spikes.We plot the mean orbital frequency of the binary  as a function of the inspiral time  in Figure 15 for the rotating and the non rotating models.Comparing the evolution of the non rotating model in Figure 15 to that in Figure 5, we find that there are minimal differences in the evolution with and without relativistic precession.On the other hand, we immediately notice that the inclusion of precession severely dampens the effect of the rotational effects that were evident in the non precessing simulations.Although, we find the same qualitative 10 5 10 4 10 3 10 2 q 10 1 10 0 10 1 10 2 10 3 t ej [yr] t ej = 5 yr t ej q 1.5 q = 4.8 × 10 4 Figure 14.The mean ejection time  ej as a function of the binary mass ratio .The dots represent the values calculated from the three-body simulations whereas the dashed line represents the best fit to the datapoints.We find that  ej ∝  −1.5 .This implies that for any  < 4.8 × 10 −4 ,  ej > 5 yr, larger than the inspiral time of a LISA detectable IMRI.For  = 10 −4 or lower, the ejection time is much larger implying that the backreaction on to the halo is expected to be minimal.This is consistent with the findings from our full  -body simulations for  = 10 −4 where we did not find any substantial feedback on the DM density profile.

Model
|Δ (2) | Non rotating 2.4 × 10 5 Retrograde rotation 2.6 × 10 5 Prograde rotation 1.9 × 10 5 Table 2.An estimate of the dephasing in the second harmonic |Δ (2) | of the GW signal for  = 10 −3 models in  sp = 7/3 spike when relativistic precession effects are included.We find that the net dephasing for the non-rotating model is consistent among the non-precessing simulations and the precessing simulations.However, the effects of rotation are significantly dampened upon the inclusion of rotation leading to lower estimates of dephasing than before.Evolution of the mean orbital frequency of the binary  in Hz as a function of time  in years for  = 10 −3 in a  sp = 7/3 spike.We find that precession dampens the effects of rotating spikes substantially.In the nonprecessing scenario the difference in dephasing between the non-rotating and rotating models was O (10 5 ) GW cycles, whereas in the precessing scenario that difference drops to O (10 4 ) GW cycles.Nevertheless, we find that the prograde rotation model takes longer to merge than the non-rotating model while the retrograde rotation model merges faster.This indicates that the results of our rotating models are robust, at least qualitatively.effects as in the non precessing simulations, i.e., retrograde model merges faster than the prograde and non rotating model, the results indicate that precession reduces the effects of rotation.Whereas in the non-precessing retrograde model, the binary takes 433.4 fewer days to merge compared to the inspiral in vacuum, in the precessing retrograde model, the binary takes 145.4 fewer days to merge.In the non-precessing prograde models, we found that the binary actually took 227.8 days longer to merge than in vacuum but we find that upon the inclusion of precession it merges 105 days earlier than in vacuum.We note that this is still ∼ 38 days slower than the non rotating model.The changes in the inspiral time are also reflected in the amount of dephasing over the 5 year inspral timespan.We present the number of dephasing cycles of the second harmonic |Δ (2) | over the course of the full inspiral in Table 2.The amount of dephasing for the non rotating model is consistent among both sets of simulations, i.e. precessing and non precessing, but we notice differences in the rotating models.In the non-precessing retrograde model |Δ (2) | ≈ 7 × 10 5 but in the precessing retrograde model that drops to 2.4 × 10 5 .We note that in the non precessing prograde models, we required 5 × 10 5 more cycles compared to vacuum for the binary to merge but due to the decreased effect of rotation, the binary merges faster, taking about 1.9 × 10 5 fewer cycles.The differences in dephasing between the non rotating and rotating models in the non-precessing case was O (10 5 ) but upon the inclusion of relativistic precession, that difference drops to O (10 4 ).
The reason for the drastic changes in the rotating models upon the inclusion of precession is unclear.Sesana et al. (2011) noted that in self-consistent -body simulations where Newtonian precession of an extended gravitational system played a role, the effects of rotation were dampened.Here we observe a similar effect but with relativistic precession indicating the similarities between the two scenarios.Notably, the effect of the eccentric binary on the secular evolution of a particle can be described through the lens of the eccentric Kozai-Lidov effect (e.g., Merritt et al. 2009;Merritt 2013).Recent simulations of hierarchical MBH triplets have found that inclusion of PN effects can diminish or even extinguish the Kozai-Lidov evolution of the inner binary (Tanikawa & Umemura 2011;Bonetti et al. 2018;Mannerkoski et al. 2021;Koehn et al. 2023).These simulations may hold clues to understanding the effect of precession and how it affects our results.However, a more thorough analysis requires future work.
Although the effects of rotation are decreased in the  = 10 −3 models, we expect rotation to still play a significant role in lower mass ratio binaries where the mass enclosed within the orbit of the binary would be larger.In such a case, the Newtonian precession of the spike could also play a part and the net change in precession could be used as an indicator for the presence of DM.However verifying this requires simulations that are beyond the scope of this current study.

Accretion effects
The secondary is also expected to accrete from the spike during the inspiral.For the parameters considered in our simulations, the accretion effect is expected to be quite minimal, especially when  2 = 1 ⊙ .For the stellar mass BH scenario with  2 = 10 ⊙ , we can estimate the rate of accretion assuming that the secondary undergoes Bondi-Hoyle accretion (Bondi & Hoyle 1944;Edgar 2004;Macedo et al. 2013;Mach & Odrzywołek 2021).Following Yue & Han (2018), the change in mass of the secondary over time (   2  ) can be written as where  is the velocity of the secondary.Assuming that  DM ∼ 10 20  ⊙ pc −3 near the secondary, and  ∼ 2 × 10 4 kms −1 , we find that   2  ∼ 0.005 ⊙ yr −1 .Since the spike gets disrupted within the first 0.1 yr, we expect the accretion onto the secondary to be of the order of 10 −4  ⊙ having minimal effects on our results.One can imagine that in lower mass ratio scenarios where the spike is not disrupted as much, the accretion effects might be larger, but still subdominant compared to dephasing due to three-body scattering.
We note that   2  ∝ 1  .Since the relative velocity between the DM particles is lower in prograde rotating models than isotropic or retrograde models, the accretion effects could be larger.An estimate of the difference is beyond the scope of this current study but will be considered in the future.Nichols et al. (2023) recently explored a self-consistent treatment of accretion in IMRIs with stellar mass BHs as the secondary and found that inclusion of accretion can lead to difference of 100 − 1000 GW cycles compared to the models where accretion is not included.Although this represents about < 1 per cent difference in the number of dephasing cycles in our  = 10 −3 and 10 −4 models, future studies will need to account for accretion to generate LISA waveforms since matched filtering requires a proper waveform determination over a few cycles.More importantly, the change in eccentricity due to accretion would also need to be studied considering the large effect eccentricity plays in GW radiation dominated binary evolution, as demonstrated in this work.

Dark matter annihilation and EM signatures
In addition to modifying GW signals of IMRIs, DM spikes are considered to be a source of gamma radiation as a result of DM annihilation.Thus, they provide strong indirect signatures of weakly interacting DM particles and allow us to probe DM microphysics (e.g., Gondolo & Silk 1999;Ullio et al. 2001;Fields et al. 2014;Shapiro & Shelton 2016;Lacroix 2018).In our work we do not consider the effect of DM annihilation on the spike profile.Our models are, therefore, representative of the non-annihilating DM case.Nevertheless, we provide estimates on how the spike profile can be changed and discuss how that would affect our results in case of annihilating DM.
For self-annihilating DM, the spike profile is significantly depleted due to the interactions between DM particles.It has been suggested that near the central MBH, a flat plateau, or core forms as a result of DM annihilation (Gondolo & Silk 1999, but see also Vasiliev (2007); Shapiro & Shelton (2016)).The density of this plateau  pl is given as where   is the mass of the DM particle, ⟨⟩ is the interaction cross section, and  represents the age of the MBH.The interaction crosssection of the DM particles is considered to be a constant in standard thermal weakly interacting particle (WIMP) models.Assuming a scenario with   = 1TeV, and  ≳ 10 6 yr we can calculate the density of the plateau and the radius of the plateau ( pl ) by setting  DM ( pl ) =  pl ( pl ) for various values of ⟨⟩.We present this information in Table 3.We find that unless the DM annihilation cross section is very small, the density of the plateau is lower than  the density of the spike where the binary is situated.The lower density would lead to a smaller amount of dephasing.Since the dephasing is proportional to the density of DM near the binary, as argued before, when ⟨⟩ = 10 −27 and 10 −30 cm 3 s −1 , the dephasing would be reduced by factors of 50 − 5 × 10 5 .In such a case, the dephasing would be significantly reduced, and non-detectable in the ⟨⟩ = 10 −27 cm 3 s −1 case.However, it should be noted that, according to MAGIC Collaboration (2016) the upper limit on the cross-section is ∼ 10 −25 cm 3 s −1 and the annihilation cross-section in reality could be much lower.In such a scenario, one could observe GW dephasing in tandem with an EM signature.We note that the above mentioned case is an optimistic scenario.Under a different scenario, we consider  ∼ 10 10 yr, typical for the ages of nearby galaxies.Using more conservative values of   = 35GeV and  = 10 10 yr leads to a plateau density of 2.9×10 9  ⊙ pc −3 (assuming that ⟨⟩ = 10 −27 cm 3 s −1 ) which would leave no imprints on the GW signal.It should be noted, however, that over such a long duration, the system is not expected to be isolated.
Alternately, one could use the detection of GW signals from DM spikes to put upper-limits on the annihilation cross section of DM.Setting  pl =  DM , we get where we used equation 7 for  DM .For  1 = 10 3  ⊙ with a  sp = 7/3 spike,  sp ≈ 0.5 pc.As in our previous approximations, we take  sp = 226 ⊙ pc −3 .We assume that the annihilation radius  pl ∼  = 2 × 10 −8 pc, the semi-major axis of the binary.Since the dephasing due to three-body scattering depends on the local density as mentioned in previous sections, we should expect to obtain a similar amount of dephasing in the annihilation scenario as in the non-annihilating scenario.We, then, obtain the following upper-limit on the cross section: Similar analysis can be performed with different spike parameters and central IMBH masses.We point to the reader, however, that a lower amount of dephasing than expected can arise from different DM properties that lead to a larger annhilation, or from dynamical factors such as rotation, as pointed in this study.In such a scenario, a different signature, possibly electromagnetic would be needed to break the degeneracy.In any case, as Hannuksela et al. (2020) points out, any detection of DM spike using GWs will be in strong tension with current thermal WIMP models and place constraints on the mass of the particle.
Other models can suggest a lower cross section allowing for the co-existence of EM signals along with GW signals from the spike (e.g., Shelton et al. 2015).However, Hannuksela et al. (2020) report that in such a case, even in the most optimistic scenario with  1 = 10 6  ⊙ , the electromagnetic counterpart will be detectable only up to a distance of 90 Mpc by detectors such as ASTROGAM (de Angelis et al. 2018), FERMI, or CTA.Only a few IMRI events would happen in such a small volume and an EM counterpart would require a large fraction of IMRIs to be embedded in DM spikes.Since the luminosity of the gamma radiation  is proportional to the squared density of the DM spike profile  2 DM , the EM signals from IMBHs are expected to be weaker and the prospects for finding EM counterparts are more pessimistic.

Implications for binary inspiral in realistic environments
As we have seen in the previous sections, the amount of dephasing is sensitive to the density profile near the IMBH.Furthermore, the actions of the inspiraling compact object can have a drastic effect on the spike, completely unbinding it.This implores us to ask an important question: how realistic are the spikes considered in this study and previous studies?How do realistic spikes affect the dephasing of an inspiraling compact object?While a comprehensive study to understand spike profile in realistic environments is beyond the scope of this work, we present a brief analysis of the impact of surrounding environments and past inspirals on the spike density profile.
In a non-isolated environment interactions between the spike and the surrounding material can reduce the density of the spike compared to isolated adiabatic growth models.Previous works have only considered the effect of stars surrounding the spike (e.g., Merritt et al. 2007a).Using the Fokker-Planck code Phaseflow (Vasiliev 2017) we study, for the first time, the effect of stellar mass BHs along with stars on the final DM spike profile.Here we only consider the case where the central SMBH has a mass of 10 6  ⊙ and the total mass of the stellar mass BH particles is 1% of that of the total mass of stars in the galactic nucleus.The total stellar mass in the nucleus is set to be 10 7  ⊙ .We note that this excludes the bulge mass.Our two-component model is consistent with the expectation from the Kroupa initial mass function (IMF) (Kroupa 2001) The stellar mass BH particles are given masses 10 times larger than that of the star particles.The DM halo has a Hernquist profile initially representing a  = 1 slope in the center and the stars and BHs in the nuclei are given a shallow  = 0.5 cusp.
In Figure 16, we plot the slope of the DM component at different times and find that in an equilibrium state the spike reaches a  = 1.5 profile near the MBH.This is similar to what is observed in a stars only model.However, due to mass segregation in the presence of a two component mass spectrum, the rate of growth of the spike is enhanced.Comparing our two component mass model to that from a single component model, we find that presence of two mass species accelerates the growth by a factor of 4. In fact, this is quite sensitive to the IMF and the initial density profile of the nuclei.A top-heavy IMF (e.g., Chabrier 2005) results in an even faster growth.This can have major implications for the time taken by the spike to regrow after it has been disrupted.A detailed study of regrowth time after disruption will be considered in a future work.
A major caveat of the above result is that we do not consider the effect of stellar evolution on relaxation.Stellar evolution can lead to mass loss and kicks to compact objects over time which reduces their population and therefore affects the relaxation timescale.We consider the effect of stellar evolution in an approximate manner by evolving a population of stars drawn from a Kroupa IMF with different metallicities and in different environments using SSE (Hur- The DM profile is embedded in a nuclei consisting of stars and stellar mass BHs, drawn from a Kroupa IMF, and the subsequent evolution is performed using a Fokker-Planck model.We find that, similar to the stars only model, the spike reaches a  = 1.5 profile at the center, near the MBH.However, unlike the stars only model, two body relaxation is enhanced due to the presence of a two component mass function, leading to an accelerated growth.ley et al. 2000).In a dense environment such as a galactic nucleus where the escape velocity is higher, we find the mass-ratio of BH particles approaches  BH / * ∼ 10 −3 of the mass of stellar particles after 5 Gyr of evolution when the metallicity  = 0.1 ⊙ where  ⊙ is the solar metallicity.This drops to  BH / * ∼ 10 −4 when  =  ⊙ .The decrease in mass is due to stellar evolution mass loss and random kicks imparted by SSE onto compact objects during their formation.In a globular cluster like environment with a total mass of 10 6  ⊙ , the escape velocity is lower and we find heavier compact objects are hardly retained ( BH / * ∼ 0).This affects the timescale of growth of the DM spike.Using Fokker-Planck models with the evolved profile, we find that when  BH / * ∼ 10 −3 , the timescale of spike growth increases by 40% compared to the twocomponent model without stellar evolution where  BH / * ∼ 10 −2 .When  BH / * ∼ 10 −4 , the timescales increases by a factor of 2 compared to the non-stellar evolved model.This highlights the importance of inclusion of a mass-species which accelerates the growth of a spike even when the population of the heavier mass species is much smaller than lower mass species.
Once a spike has been disrupted, the regrowth happens on timescales that are on the order of collisional relaxation time within the sphere of influence of the MBH.The relaxation time is affected by the mass species surrounding the MBH.In a single component model, relaxation takes longer than when a mass spectrum is present, as evident from our results above.We can estimate the relaxation time  relax in a stellar only environment as follows (e.g., Babak et al. 2017;Becker 2024) : where ln(Λ) is the Coulomb logarithm, is the velocity dispersion, and  infl is the influence radius of the MBH. can be estimated from the well known  −  relationship (e.g., Gültekin et al. 2009;Kormendy & Ho 2013) where  is the ratio of plunges to inspirals, and  is the characteristic mass of the compact object.Taking  = 0 to find the upper limit and  = 10 ⊙ representing BHs , we find that   ≲ 0.1 Gyr.Thus, the time between inspirals is comparable to the two body relaxation time in high density environments and is shorter in the low density environments.As such, we expect the spike to not exist in its equilibrium state in low density environments.However, under the presence of a mass spectrum the relaxation time decreases in which case   ∼  relax which can lead to the spike existing in the equilibrium  = 1.5 state, even in low density environments.This highlights the importance of examining the spike growth embedded in a realistic mass spectrum.We note one caveat of this calculation is that relaxation times after depletion of a spike a larger than the ones calculated from the above equation.Therefore, our relaxation timescale should be considered as a lower limit.
A big question also lingers regarding the starting point of the compact object in the simulations.In a realistic scenario, a compact object would be unbound initially and become bound over time due to the effect of dynamical friction, two-body relaxation, and hardening in the stellar and DM environment.In a stellar dominated environment the compact object is going to be driven to inspiral orbits mainly due to the effects of two-body relaxation.In such a case, we would expect the initial inspiral to cause minimal disruption to the spike.However, in such a case, the spike profile would have a  sp = 1.5 slope which we showed produces minimal dephasing effects for the systems considered.This was also highlighted in Becker (2024) where the author found that GW signals from the spikes in this case would be hardly detectable.On the other hand, in an isolated environment with no erosion, the compact object is going to inspiral due to the effect of dynamical friction and three body hardening from the DM environ-ment.In such a case, the formation and inspiral itself would produce a flat core that can massively reduce dephasing effects.As an example, we produce the impact on the density profile from a simulation with a  sp = 7/3 spike where the compact object now starts at an initial semi-major axis of  0 = 2 × 10 −3 pc in Figure 17.Within 0.1 Myr, we find a core has formed with a density of ∼ 10 7  ⊙ pc −3 .According to Merritt et al. (2007b), the formation of a hard binary is accompanied by an ejection of mass comparable to the mass of the binary.Therefore in DM only environments, we expect the density to be even lower and the core to be larger.The relaxation time of the spike in this isolated environment is so long that we do not expect a regrowth within a Hubble time.This indicates that the simulations with the  sp = 7/3 models represent optimistic scenarios where the spike exists in isolation and has not undergone a previous inspiral.This is a major caveat of our work and that of previous works.Additionally, on account of their lower masses, IMBHs can be off-center where the adiabatic growth of the spike can be diminished.In such a scenario, a  = 1.5 is typically formed, as noted by Zhao & Silk (2005).All of the above noted issues are pertinent and we emphasize the need to understand spike growth in realistic environments.We leave this as a future study.

Prospects for GW detection
As shown in Figure 3, IMRI sources are promising for multiband GW astronomy with LISA and potential decihertz detectors.Several formation scenarios have been proposed for the formation of IMRIs in DM spikes, including host sites such as the nuclear star clusters of dwarf galaxies (Yue et al. 2019) and merger remnants of elliptical galaxies (Vázquez-Aceves et al. 2023).The merger rate of such systems with and without DM spikes, however, is still uncertain.More detailed population synthesis of these sources with either semi-analytic models, cosmological simulations, or some combination thereof would provide more insight into what merger rates might be possible.Conversely, the GW detection of IMRI sources and constraints on potential DM environments would provide tests of astrophysical population models in addition to DM physics.
In light of the critical examination of realistic spike profiles from the previous section and Becker (2024), the prospects for dephasing due to DM spikes on IMRIs appear to be somewhat not optimistic.However, we emphasize that a full parameter space study is required before we can conclusively determine whether spikes would have tangible dephasing impact on GWs from IMRIs.On the other hand, Dai et al. (2022Dai et al. ( , 2023) ) and Becker (2024) noted that inspirals often happen on highly eccentric orbits with larger semi-major axes.In such a case, the DM spike can lead to a periapse precession in addition to the relativistic precession.The Newtonian precession of the spike is opposite in direction to the relativistic precession and the net effect can be quantified as a deshifting index (Becker 2024), which is a measure of the change in GW cycles due to change in pericenter precession from the DM spike.Dai et al. (2023) noted that deshifting can be present even in low density environments.Our preliminary analysis suggests that deshifting may be a more optimistic signature than dephasing.Unfortunately, for the set of models considered in this work, we were not able to quantify deshifting.In the future, we plan on studying an extended parameter space where we study dephasing as well as deshifting across spike models.
A potential space-based decihertz GW detector provides a significant advantage for boosting the IMRI detection rate as well as breaking parameter degeneracies that might be encountered with just LISA data.In addtion, if the inspiral phase of IMRIs is at a sub-threshold level in LISA data, the resolved merger phase in a space-based decihertz GW detector can provide priors to search for the subthreshold inspiral phase in archival LISA data.This is similar to how LIGO and 3G/XG ground-based detectors like Einstein Telescope or Cosmic Explorer may provide prior information on stellar-mass BH binary mergers to search for potential inspiral counterparts in archival LISA data (Ewing et al. 2021).
In addition to the proposed tests of dynamical friction and accretion on gravitational waveforms, one can also test for the presence of three-body/loss-cone scattering for a potential IMRI event.One caveat for testing the presence of such effects might be that nonlinear feedback and DM-spike relaxation may complicate waveform parameterization.Though this certainly motivates more numerical simulations of IMRIs in DM spikes in order to explore the parameter space more comprehensively and to inform waveform parameterizations that are interpretable.

CONCLUSIONS
IMRIs are considered among the most crucial sources of lowfrequency gravitational waves (GWs) that future space-based GW detectors like LISA and DECIGO can potentially detect.Recent studies have suggested that if the intermediate-mass black hole (IMBH) in these systems is surrounded by a dark matter (DM) spike, the gravitational effects of the spike may induce dephasing in the observed GW signal, serving as a distinctive signature of the DM spike.Prior investigations have predominantly employed analytic methods to model the interaction between the spike and the binary, treating the gravitational impact of the spike on the binary through dynamical friction (DF).However, these approaches neglect the influence of the binary on the spike itself, a factor that can significantly affect the dynamics and, consequently, the degree of dephasing.
In our study, we employ -body methods to delve into the complex dynamics involved in eccentric IMRIs embedded in both non-rotating and rotating DM spikes, evolving the system self-consistently.Our work represents the first attempt at modeling such systems using -body simulations.The simulations reveal that, contrary to previous assumptions, the primary mode of dissipating energy from the system is not DF, a cumulative effect of two-body encounters, but rather three-body scattering, similar to stellar loss-cone scattering in SMBH/MBH binaries.We note that during the submission of this work, Kavanagh et al. (2024) also finalized their work on -body simulations employing methods similar to ours.They qualitatively find similar behavior of the IMRI but find that the interactions can be described using a semi-analytic method combining dynamical friction loss with the time dependent potential of the binary.While their simulations are run for a shorter duration than ours, the combined effect which they propose is of three-body nature and similar to what we describe in this work.
The main conclusions of our work can be summarized as follows: • The evolution of the non-rotating models reveals that in larger mass ratio cases, the binary is highly efficient at rapidly degrading the spike, resulting in minimal dephasing effects.Depending on the mass-ratio of the binary, our dense spike models with  sp = 7/3, 9/4 show that dephasing of O (10 4 ) − O (10 5 ) can be expected.However, in lower density  sp = 3/2 spikes produced in a collisional environment, we find minimal dephasing of ≤ O (10).
• We find that all of our simulations with dense spikes ( sp = 7/3, 9/4) predict that the dephasing is much larger, by factors of 10-100, than that predicted by the self-consistent semi-analytic method HaloFeedback for similar binaries in circular orbits.This indicates that DF theory is unable to fully capture the dynamics of the binary.
• In our rotating models with dense  sp = 7/3 spikes, we find that, spikes that counter-rotate with the binary lead to faster inspirals compared to spikes that co-rotate with the binary, which lead to slower inspirals, even slower than inspirals in vacuum.This leads to dephasing effects that are 2.5 − 3.5 times higher in the retrograde/counterrotating cases than their non-rotating counterparts.
• We use three-body simulations to investigate the nature of the interaction of the binary with spike particles and find that the binary primarily dissipates energy via three-body interactions rather than DF, a two-body effect.The impact of DF in all of our simulations is sub-dominant shedding some light on the discrepancy between our results and those obtained in previous studies.Additionally, threebody simulations reveal that in rotating spikes, particles on retrograde motion eccentrify the binary whereas particles on prograde motion circularize the binary.This indicates that rotation in spikes cannot be neglected and should be the subject of further investigations in the future.
• We conduct simulations to investigate the impact of relativistic precession on dephasing in our  = 10 −3 models.Our findings indicate that precession has minimal to no effect on dephasing in the non-rotating models but significantly diminishes the effects of rotation.
• By examining our initial models with Fokker-Planck methods, we assess the presence and detectability of spikes in realistic environments.Our results suggest that non-isolated environments have DM spikes with shallower slopes than previously considered, leading to smaller signals and lower detection prospects via dephasing.We show that the densest DM spike models with  sp = 7/3 represent scenarios where the spike is isolated and has not undergone any previous inspirals, as any inspiral leads to the formation of a lowdensity core where the regeneration of the spike takes much longer than the Hubble time.In the presence of lower density  sp = 3/2 spikes, present in non-isolated environments, dephasing is minimal and non-detectable.High stellar and stellar-mass black hole densities can accelerate the growth of  sp = 3/2 spikes, but a more extensive parameter space study is needed to determine which E/IMRIs embedded in such spikes can produce detectable dephasing signatures.Our preliminary analysis, coupled with recent studies, suggests that "deshifting" rather than dephasing might be a more optimistic signature, as it is more robust even in low-density environments (Dai et al. 2022;Becker 2024).
Although our work uses idealized models, it provides a foundational exploration that paves the way for more comprehensive analyses in the future.Despite critical uncertainties, particularly regarding the growth, retention, and population of DM spikes around IMBHs, our research highlights the promise and potential benefits of delving deeper into this intricate problem.We demonstrate the need to consider more realistic environments, as detection prospects depend sensitively on spike density, which is influenced by the surrounding environment.Our work also underscores the importance of multiband GW astronomy and the need for a dedicated decihertz detector.By enhancing our approach to include accretion onto the inspiraling object and incorporating higher-order post-Newtonian terms, our -body code emerges as an effective method for simulating the evolution of IMRIs embedded in DM spikes. .The evolution of the relative change in the semi-major axis Δ/ 0 as a function of the number of orbits of the binary.The solid lines indicate the average of five independent simulations whereas the shaded region indicates the standard deviation.We find that the results from method presented in this work (purple line) agree with that from Taichi (orange line), ph4 (green line),  -body codes where the self-gravity of the spike is not neglected.This indicates that the effect of the DM-DM interactions is minimal and can be safely ignored, validating our method.Also presented is the evolution of the same system using HaloFeedback (black line).While the results from HaloFeedback and the  -body codes agree over the first 100-150 orbits, they diverge after that.The spike is able to dissipate more energy from the binary in the  -body models compared to the HaloFeedback models since the impact of DF is subdominant compared to three-body effects.
from the five simulations whereas the shaded regions represent the standard deviation.We find that the evolution using our approximate force -body method is consistent with that obtained from Taichi and ph4 where the DM-DM interactions are not neglected.This indicates the effect of the self gravity of the spike is minimal and can be safely neglected.However, we find major differences between the results from the -body codes and HaloFeedback.The three methods produce consistent results for the first 150 orbits after which the HaloFeedback results deviate because the binary enters the hard binary phase in which three-body interactions play a dominant role.The binary stalls at Δ/ 0 ≈ −5 × 10 −4 in the HaloFeedback models whereas a longer evolution indicates that the binary actually stalls at Δ/ 0 ≳ −2×10 −3 in the -body simulations.This is almost 4× larger than what is predicted by HaloFeedback.As a result of this discrepancy, the dephasing obtained from the HaloFeedback models should be smaller than those obtained from the -body models, a fact discussed in the results section.The discrepancy arises from the fact that DF is subdominant to three-body scattering as we showed in section 4.3.Moreover, the three-body effects occur over longer timescales, indicating that semi-analytic methods should be validated against long term, secular -body simulations.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.A visual representation of an eccentric IMRI embedded in a DM spike.The mass of the central BH is  1 and that of the inspiraling object is  2 .This figure has been inspired by Figure 1 of Kavanagh et al. (2020).

Figure 2 .
Figure 2. The relative force accuracy | f |/|f | of DM particles as a function of their separation  from the central IMBH with mass  1 = 10 3  ⊙ .The spike follows a  sp = 7/3 density profile.The inspiraling object, whose mass is  2 = 1 ⊙ , is situated at an initial semi-major axis of  0 = 2 × 10 −8 pc.We find that in the region of interest, the force accuracy is ≤ 10 −5 , which is comparable to the the force accuracy obtained in tree/FMM based codes.
Figure3.The characteristic strain of the IMRIs used as our initial conditions as a function of the frequency  in vacuum for different mass-ratio  models.The central IMBH has a mass of  1 while the inspiraling object has a mass of  2 .The binary is defined by its initial semi-major axis  and eccentricity .The luminosity distance is denoted as   .All of the IMRIs presented here have a merger time of ∼ 5 years.Since eccentric binaries radiate in multiple harmonics, we have plotted the strain from the second and third harmonics.In all of the scenarios, we find that the second harmonic has an equivalent or higher strain than the third harmonic.Even higher harmonics have lower strains and are harder to detect using LISA but would be detectable using DECIGO.

Figure 4 .
Figure 4.The density of the DM spike  DM as a function of the distance from the central IMBH .We use these density profiles to generate our  -body initial conditions.
Figure5.A comparison of the binary parameters for different mass ratio models evolving in a  sp = 7/3 spike with  sp calculated using equation 8.The evolution in vacuum (purple line) is calculated usingPeters (1964) analytic formula, while the evolution calculated using the Chandrasekhar DF formula assuming a static spike (orange line) is calculated using IMRIPy.They are compared to the evolution from our  -body simulations (green line).Left column: the mean orbital frequency  as a function of time  in years.Middle column: the estimated number of dephasing cycles of the second harmonic |Δ (2) | as a function of the frequency  in Hz.Right column: the eccentricity  as a function of the semi-major axis .We notice that in higher mass ratio models, the evolution of the binary is similar to that in vacuum.For  = 10 −2 model, there is a 100× reduction in the estimated number of dephasing cycles compared to the evolution calculated using Chandrasekhar DF formula as the spike has been disrupted in a very short time span.As we decrease the mass ratio, the disruption decreases.We notice that in case of the  = 10 −3 model, the dephasing is only reduced by 3× compared to the evolution calculated using IMRIPy.For  = 10 −4 model, we find that dephasing is a factor of 3 larger than what we obtain using the Chandrasekhar formula, with little to no disruption of the spike.This signals that the Chandrasekhar DF might be insufficient to explain the evolution of the binary in DM spikes.

Figure 6 .
Figure6.The dephasing of the second harmonic |Δ(2) | as a function of the binary frequency  in Hz for a  = 10 −3 binary embedded in a  sp = 9/4 spike.The color scheme is the same as that used in Figure5.We notice that similar to the  sp = 7/3 case, the amount of dephasing in our  -body models is reduced by ∼ 3× compared to the DF models.

Figure 7 .
Figure7.Binary frequency  as a function of time  for the  = 10 −3 ,  sp = 1.5 model.We find that the low density of the spike results in no discernable differences between the inspiral in vacuum and that in the spike.This results in ≤ O (10) cycles which might not be detectable.

Figure 8 .
Figure8.The density of the DM spike  DM as a function of the distance  from the primary in pc for the  = 10 −3 model in a  sp = 7/3 spike.The different colors represent the evolution of the density profile over time with darker colors representing earlier times and lighter colors representing later times.We notice that the inspiral of the IMRI leads to a significant disruption in the spike near the semi-major axis  of the IMRI.The density of DM in that region is reduced by a factor of 100 compared to the original density after an inspiral time of 2 years.This is qualitatively consistent with the findings fromKavanagh et al. (2020).However, we find that the disruption to the spike due to the IMRI in our  -body is much slower than what is found byKavanagh et al. (2020) who find that the spike is significantly disrupted within 0.1 years.

Figure 9 .
Figure 9. Similar to Figure6but for  = 10 −4 binary embedded in a  sp = 7/3,  sp = 0.54pc spike.We notice that the net dephasing predicted by the body simulation is comparable to that predicted by the IMRIPy simulation, although at higher frequencies the dephasing falls off faster in the  -body models.The density in the region of interest for this particular model is about 6 times lower than that for the  = 10 −4 ,  sp = 1.17pc model and we notice a similar decrease in the amount of dephasing in this scenario compared to that model.
Figure10.A comparison of the binary evolution parameters for different mass ratio models inspiraling in a  sp = 7/3 spike similar to Figure5but also including the evolution from rotating models.We notice that the in all of the retrograde models, the binary merges faster than in the prograde and even the non-rotating models.This leads to a major enhancement in the number of dephasing cycles, by as much as 3× that of the non-rotating model in the case of  = 10 −4 .Since we are plotting the absolute value of the number of dephasing cycles, the dephasing in case of the prograde models appears to be positive.However, the prograde models actually merge slower than even the vacuum models and Δ (2) is actually negative (indicated by the dashed lines) .Thus the number of GW cycles in prograde models is larger than that in vacuum.We also find that in the retrograde models, the binary circularization rate is slower than that in the prograde models.In the latter scenario, the circularization is enhanced leading to a slower inspiral.As the mass ratio decreases, the effect of rotation becomes more prominent suggesting that transfer of angular momentum between the spike and the IMRI contributes majorly and cannot be ignored.

Figure 11 .
Figure11.The distribution of the normalized energy change Δ Ẽ of strongly interacting DM particles from our  -body simulations compared to that from three-body scattering simulations and dynamical friction approximation.We find that the distribution of energy change of strongly interacting particles is consistent between the  -body simulations and the three-body simulations but not with the calculated values using the dynamical friction approximation.This confirms our hypothesis that three-body scattering, and not dynamical friction is responsible for dissipating energy from the binary.
Figure12.Top: probability distribution of the normalized energy change Δ Ẽ of 2500 DM particles from the three-body simulations performed until the particle is ejected by interactions with the binary for both prograde and retrograde rotation models with two different binary eccentricities .We notice that in both prograde and retrograde models, the distribution of the normalized energy is similar and does not change with the binary eccentricity.Bottom: similar to top but the probability distribution of the normalized angular momentum Δ L of 2500 DM particles after ejection.We notice that there are major differences in the distribution of Δ L between prograde and retrograde models.In particular, in the retrograde models, the angular momentum of the particle is larger than that in the prograde models, especially for more eccentric binaries.This suggests that retrograde particles are more efficient at "stealing" angular momentum from the binary, especially at larger .

Figure
Figure15.Evolution of the mean orbital frequency of the binary  in Hz as a function of time  in years for  = 10 −3 in a  sp = 7/3 spike.We find that precession dampens the effects of rotating spikes substantially.In the nonprecessing scenario the difference in dephasing between the non-rotating and rotating models was O (10 5 ) GW cycles, whereas in the precessing scenario that difference drops to O (10 4 ) GW cycles.Nevertheless, we find that the prograde rotation model takes longer to merge than the non-rotating model while the retrograde rotation model merges faster.This indicates that the results of our rotating models are robust, at least qualitatively.

Figure 16 .
Figure16.Slope of the DM profile  as a function of distance  from the central MBH at various times.The time is presented in  -body units with 1 time unit corresponding to ≈ 14.9 Myr.The initial profile of the DM surrounding the MBH is a  = 1 model, representing a Hernquist type halo.The DM profile is embedded in a nuclei consisting of stars and stellar mass BHs, drawn from a Kroupa IMF, and the subsequent evolution is performed using a Fokker-Planck model.We find that, similar to the stars only model, the spike reaches a  = 1.5 profile at the center, near the MBH.However, unlike the stars only model, two body relaxation is enhanced due to the presence of a two component mass function, leading to an accelerated growth.
Figure A1.The evolution of the relative change in the semi-major axis Δ/ 0 as a function of the number of orbits of the binary.The solid lines indicate the average of five independent simulations whereas the shaded region indicates the standard deviation.We find that the results from method presented in this work (purple line) agree with that from Taichi (orange line), ph4 (green line),  -body codes where the self-gravity of the spike is not neglected.This indicates that the effect of the DM-DM interactions is minimal and can be safely ignored, validating our method.Also presented is the evolution of the same system using HaloFeedback (black line).While the results from HaloFeedback and the  -body codes agree over the first 100-150 orbits, they diverge after that.The spike is able to dissipate more energy from the binary in the  -body models compared to the HaloFeedback models since the impact of DF is subdominant compared to three-body effects.

Table 3 .
The annihilation radius  pl and annihilation plateau density  pl under different DM cross sections ⟨ ⟩ for a  sp = 7/3 spike with a central BH with mass  1 .
Babak et al. 2017;Becker 2024/3 DM spike  as a 10 ⊙ BH is in.Even at larger separations, we find that the inspiraling object forms a core extremely rapidly.This suggests that dense spikes of the form of  sp = 7/3 can only be present if there were no prior inspirals.infl.Taking ln(Λ) ∼ 10,  ∼ 50kms −1 , and  infl ∼ 0.1 pc ,we find  relax ≲ 0.025 Gyr for a 10 4  ⊙ IMBH.This is representative of relaxation time in a high density environment.On the other hand, in a low density environment,  infl is going to be larger.Assuming  infl ∼ 0.5 pc in that case, we find  relax ≲ 0.6 Gyr.The inspirals happen on timescales over which two-body relaxation depletes the cusp by driving the compact objects on loss cone orbits.The depletion timescale  d is estimated as follows (e.g.,Babak et al. 2017;Becker 2024): Nguyen et al. (2017Nguyen et al. ( , 2018) )Since we focus on IMBHs in this work, we refer to the density profiles of known dwarfs fromNguyen et al. (2017Nguyen et al. ( , 2018) )to calculate  and