Darkness visible: N-body simulations of dark matter spikes in Hernquist haloes

ABSTRACT

Dark matter is theorized to form massive haloes, which could be further condensed into so-called spikes when a black hole grows at the centre of such a halo. The existence of these spikes is instrumental for several dark matter detection schemes such as indirect detection and imprints on gravitational wave inspirals, but all previous work on their formation has been (semi-)analytical. We present fully numerically simulated cold dark matter spikes using the swift code. Based on these results, we propose a simple empirical density profile – dependent on only a single mass-ratio parameter between the black hole and total mass – for dark matter spikes grown in Hernquist profiles. We find that the radius of the spike scales differently compared to theoretical predictions, and show a depletion of the outer halo that is significant for high mass-ratio systems. We critically assess approximations of the spike as used in the field, show that our profile significantly deviates, and contextualize the potential influence for future DM detections by simulating binary black hole inspirals embedded in our profile.

1 INTRODUCTION

One of the biggest mysteries in modern cosmology is the origin of dark matter (DM). While never directly detected, it has been indirectly observed in systems ranging from small galaxies to clusters of galaxies, with more work underway focused on direct detection or solving theoretical challenges (for an overview, see Navas et al. 2024). The most recent Planck satellite survey found, assuming standard Lambda cold dark matter (⁠|$\Lambda$|CDM) cosmology, that 26.1 per cent of the total mass-energy of the universe consists of this DM, which is 83.9 per cent of the total matter density (Aghanim et al. 2020). DM is theorized to form small density perturbations in the early universe, which stop expanding after sufficient growth and collapse into a spherical halo (Navarro, Frenk & White 1996; Bullock & Boylan-Kolchin 2017).¹ While this work focuses on collisionless cold DM, as it is part of |$\Lambda$|CDM cosmology, many other forms of DM have been proposed in order to solve known conflicts with this cosmology (see e.g. Bullock & Boylan-Kolchin 2017).

Several methods of detecting DM have been proposed and tried out: from direct detection in a laboratory to indirect detection using signals from astrophysical systems (Navas et al. 2024). One of these systems is a DM halo with an embedded black hole (BH) at the centre. For potential detections, the density distributions of these systems must be well understood. This article will numerically simulate DM haloes with a central BH, and contextualize the results using the possible detection of DM through gravitational waves (GWs). There are more detection methods where the halo profile is of great importance, such as the detection of gamma-rays from DM self-annihilation (e.g. Gondolo & Silk 1999; Bertone, Zentner & Silk 2005; Aschersleben et al. 2024), but this work will only consider the impact on the GW interpretation.

Eda et al. (2013) suggested that DM haloes slow down the inspiral of a compact object into an embedded BH. This leads to a dephasing of the generated GWs compared to the vacuum solution, which could be measurable with the next generation of GW detectors (e.g. Eda et al. 2015; Yue, Han & Chen 2019; Kavanagh et al. 2020; Coogan et al. 2022). Recent numerical studies found that this dephasing is both larger than previously predicted, and the DM halo is depleted at a slower rate than anticipated, reinforcing the potential of this method (Kavanagh et al. 2024; Mukherjee et al. 2024). Future GW detectors such as the Laser Interferometer Space Antenna (LISA) GW detector should not only be able to measure this dephasing, but also differentiate between DM haloes and other environments such as accretion discs (Cole et al. 2023; Colpi et al. 2024).

LISA will be able to detect extreme and intermediate mass-ratio inspirals where the central BH has a mass between |$10^3$| and |$10^6$| |$\text{M}_\odot$| (Colpi et al. 2024). These intermediate mass BHs (IMBHs) have been detected in our universe (e.g. Farrell et al. 2009; Pasham, Strohmayer & Mushotzky 2014; Häberle et al. 2024), and there exists a catalogue of candidates from low-luminosity active galactic nuclei (Barrows, Mezcua & Comerford 2019). The formation of these IMBHs at galactic centres likely mirrors that of supermassive BHs (SMBHs), though having experienced very few major galactic mergers, leaving the surrounding DM halo relatively unperturbed.

It has been known for some time that the growth of a BH influences the distribution of matter around it (e.g. Quinlan, Hernquist & Sigurdsson 1995). Gondolo & Silk (1999) (From here on: G&S) predicted that a growing BH in a DM halo creates a DM spike, where the DM is redistributed into a steep cusp of uniform slope in the central regions. They analytically showed that if this growth is adiabatic, an initial distribution of |${\rho \propto r^{-\gamma }}$| transforms into a spike of slope |${\gamma _\text{GS} = (9-2\gamma)/(4-\gamma)}$|⁠. These spikes are predicted to start at one-fifth of the radius of gravitational influence of the central BH (Merritt 2004, 2004). While DM spikes have been of interest for the detection of DM, previous studies into their shape has predominately been (semi-)analytical. To the best of the authors’ knowledge, fully numerical N-body simulations of the formation of these DM spikes have never been published. These allow for important checks of (semi-)analytical formalisms, and should lead to more realistic DM halo density profiles.

In this article, we present the first numerical step towards realistic spike density distribution functions, by demonstrating the formation of a DM spike in N-body simulations using realistic Hernquist halo profiles. In Section 2, we give the used DM halo distribution functions and analytical predictions of the expected spikes. In Section 3, we describe the used code and numerical schemes (Section 3.1), how the initial haloes are constructed (Section 3.2), and a description of the simulated systems and their processing (Section 3.3). In Section 4, we propose a new profile for the spike (Section 4.1) and contextualize the results by estimating their impact on the detectability of dephasing in GWs (Section 4.2). Finally, in Section 5, we draw conclusions and some concluding remarks are given. Appendices with additional material regarding numerical solutions and validation of the simulations are also supplied after the bibliography.

2 THEORETICAL BACKGROUND

Many different DM density profiles exist and are actively compared to observational data (e.g. Li et al. 2020). Arguably the most prominent spherical profiles are the Navarro, Frenk & White (NFW, 1996) and the Hernquist (1990) profiles. Both feature an inner cusp of |$\rho \propto r^{-1}$|⁠, with the Hernquist having a more steep outer density falloff of |$\rho \propto r^{-4}$| compared to |$\rho \propto r^{-3}$| of the NFW.² We will use the Hernquist profile, because it does not require the use of artificial cut-offs and is easier to work with analytically. Due to the identical inner cusp, we expect the results presented in this article to be equivalent for both profiles. The density profile of a Hernquist halo is given by

where |$M_\text{h}$| is the total mass of the halo and a is the characteristic scale radius. Its gravitational potential |$\psi _\text{Hernq}$| is then given by

Simulations by Correa et al. (2015) have shown that the two free parameters a and |$M_\text{h}$| are correlated for realistic Hernquist halos, and we thus define our haloes using only |$M_\text{h}$|⁠. This article still includes a in equations for ease of reading and consistency with other works. Assuming zero redshift, this relation becomes:

where |$r_\text{vir}$| is the virial radius, defined as the radius where the encapsulated halo has a density of 200 times the critical density of the universe, |$\rho _\mathrm{ c} = H^2/8\pi G^2$|⁠, with H the Hubble parameter.³

In this study, we introduce a BH to the centre of Hernquist haloes and grow it adiabatically. Specific growth rates have been proposed for BHs, such as the Bondi–Hoyle–Lyttleton accretion rate (Hoyle & Lyttleton 1941; Bondi & Hoyle 1944), however the final state of the system is independent of the actual rate of growth within the assumption of adiabatical growth, due to the invariance of the actions of the DM halo (Binney & Tremaine 2008). Following G&S, a spike of slope |$\gamma _\text{sp} = 7/3$| is predicted for the inner regions of the halo. This spike is predicted to start at

for a Hernquist profile, where |$\alpha _\gamma$| is a factor obtained from numerical integration, being 0.122 for |$\gamma = 1$|⁠. If the growth is non-adiabatic, a milder slope is predicted (Ullio, Zhao & Kamionkowski 2001). The BHs in this study are treated as Newtonian point-masses, as relativistic effects on the haloes only become significant at radii of the order of 10 Schwarzschild radii |$R_\mathrm{ s}$| (Sadeghian, Ferrer & Will 2013). However, the exact shape of a G&S spike close to the BH could be less steep based on the BH formation history (Bertone et al. 2024). This effect is again below our resolution, such that we predict our simulated spikes to follow a power law. The growth of the BH would be fueled by a baryonic component present in the total astrophysical system. However, it can be excluded from the spike calculations as long as the total baryonic mass is much smaller than the mass of the combined BH + DM halo system, which is often the case in realistic systems and is also assumed by G&S.

Numerical studies are often performed with either the NFW, Hernquist, or toy profiles of a constant slope of |$\rho \propto r^{-1}$|⁠. A spiked DM halo profile is then often approximated as an initial profile with a spike added after a certain spike radius |$r_\text{sp}$| (e.g. Bertone et al. 2005; Eda et al. 2013, 2015; Kavanagh et al. 2020; Aschersleben et al. 2024; Mukherjee et al. 2024):

where |$r_\text{sp}$| is not taken to be the G&S prediction, but instead as one-fifth of the radius of gravitational influence |$r_\mathrm{ h}$|⁠, defined as the radius where the enclosed halo mass is equal to twice the mass of the central BH (Merritt 2004):

When solved for the Hernquist halo, the radius of the spike becomes

where |$\mu = M_\text{BH}/M_\text{tot}$| is the mass ratio. We note that equations (4) and (7) are equivalent for |$\mu \lessapprox 0.01$|⁠, but quickly deviate after. We refer to the profile of G&S, expanded using equations (5) and (6), as the Modified G&S profile. We note that this profile was originally derived for an isothermal density profile with |$\rho \propto r^{-\gamma }$| for |$0.5 \le \gamma \le 2$|⁠, and that deviations are thus expected as the Hernquist starts to deviate from the isothermal assumption.

All particles in our simulated haloes are gravitationally bound and must thus have energies below the binding energy |$\mathcal {E}$|⁠. The distribution function of these energies is the Eddington formula (Binney & Tremaine 2008):

In appendix  A, we give the solutions for both an isolated Hernquist halo, and one where a central BH is already present.

3 NUMERICAL METHODS

3.1 Code and numerical schemes

The simulations presented here are performed with a modified version of the swift code (Schaller et al. 2024).⁴ We introduced a new BH model to the code that can grow its mass over time at a uniform pace until a final mass is reached. This BH model is called the ‘Dark Matter Attracting Black Hole’ or DAB for short. The BHs in this model have a set rate of growth and a final mass, after which the growth is stopped. Expanding DAB to grow following accretion formalisms where baryonic matter is deleted is easily done, and fully supported by swift, however we have chosen not to do so in order to reduce computational complexity.

swift’s gravity solver employs the Wendland-C2 kernel, named after Wendland (1995), where the gravitational interactions between two particles are smoothed below a certain smoothing radius  |$\epsilon$|⁠. This value cannot be taken too small, as to prevent two-body interactions, nor too large, since the resolution of our simulations is determined by this softening. We follow the scheme as developed by Power et al. (2003),

where |$\alpha$| is a numerical constant set empirically. While different values exist in the literature, we find |$\alpha = 6$| to be stable in all our simulations.⁵ The smallest radius yielding physical results is called the convergence radius  |$r_\text{conv}$|⁠, the radius below which particles start to show two-body behaviour (Duffy et al. 2008). We found that |$r_\text{conv} = 2\epsilon$| for our simulations after 4 Gyrs, where we refer to appendix  B for the full analysis.

Having a large amount of particles in the simulation will thus lower |$r_\text{conv}$| and improve our resolution. However, the complexity of the calculations performed by swift scale as |$\mathcal {O}(N\log N)$| compared to the resolution scaling as |$\mathcal {O}(1/\sqrt{N})$|⁠, and the balance between a manageable computation time and a small resolution is a fine one.⁶ As the mass of the central BH grows, and larger forces are to be computed, the timestepping of the simulation gets smaller. swift allows for dynamical timestepping |$\Delta t$|⁠, where the time-steps are scaled to the forces acting upon the particle. The time-stepping is defined as (Power et al. 2003)

where we found that |$\eta = 0.005$| in our simulations. Due to this complexity and ever decreasing timestepping as the mass of the central BH grows, individual simulations took up to two weeks to finish.

Finally, this article employs a numerical implementation of the G&S formalism to compare our results to previous studies, allowing us to compute results for the exact same Hernquist haloes as the main N-body simulations. We refer to this as the Numerical G&S formalism, and it is based on code developed for Bertone et al. (2024).

3.2 Initialization of the haloes

Every DM particle in our simulations has seven properties: positions |$(x,y,z)$|⁠, velocities |$(v_x,v_y, v_z)$|⁠, and a mass set equal across all DM particles. The initial conditions are chosen so that an isolated halo is stable over the time span of the simulation, in the following way:

• N DM particles are created and radially distributed using inverse transform sampling. The azimuthal and polar angles are sampled uniformly.

• A BH particle of desired mass is placed in the centre. In practice, this was set to be the DM mass as to not have a large impact on the distribution function.

• An array of binding energies |$\mathcal {E}_x$| ranging from 0 to |$\psi$| – equation (2) – is created for each DM particle. The distribution function of the energies in this array is numerically calculated,

  $$\begin{eqnarray} F(\mathcal {E}_x) = \frac{\int _0^{\mathcal {E}_x} f(\mathcal {E})\sqrt{\psi -\mathcal {E}} \mathrm{d}\mathcal {E}}{\int _0^{\psi } f(\mathcal {E})\sqrt{\psi -\mathcal {E}} \mathrm{d}\mathcal {E}}. \end{eqnarray}$$

  (11)

  where |$f({\mathcal {E}})$| is the solution of the Eddington formula of equation (8) (see also Appendix  A).

• |$\mathcal {E}_x$| is sampled from |$F(\mathcal {E}_x)$|⁠, resulting in the velocity of the particle using |$\left|{\vec{v}}\right| = \sqrt{2(\psi -{\mathcal {E}_x})}$|⁠. Its direction is taken at random.

Using the above, all particles in the simulation are both bound and isotropically distributed as long as |$N \gg 1$|⁠,⁷ but have otherwise completely free orbits. The central BH is not bound in place and will move due to its gravitational interactions. This movement is completely random, but increases as the spike forms and the central densities around the BH increase. The halo is bound to the BH and it follows it around, making it so the final results are not significantly influenced by this movement.⁸

These initial conditions produce a very small initial shock wave when the simulation is started with a BH present. This effect is thoroughly treated in Appendices  A and C3, where we conclude that this of no significant effect to our final results.

3.3 Simulated systems and processing

Simulations and detections show that Super Massive BHs in the centres of galaxies are orders of magnitude lighter compared to their surrounding haloes (e.g. Ferrarese 2002; Bandara, Crampton & Simard 2009; Booth & Schaye 2010). However, simulating these kinds of systems would place the spike fully within |$r_\text{conv}$| and thus unresolvable. We furthermore use IMBHs due to the computational difficulty of a system with a SMBH and halo of the same order of magnitude. Therefore, this paper simulates systems where the final BH mass is between |$10^3 \,\text{M}_\odot$| and |$5\times 10^3\, \text{M}_\odot$|⁠, and the mass of the corresponding halo is of equal order of magnitude or one higher. As a result, the mass ratios |${\mathit{ M}}_\text{BH}/{\mathit{ M}}_\text{h}$| in our study are higher than those treated by previous studies (e.g. Bandara et al. 2009; Booth & Schaye 2010).

We deem that these systems will yield results that can also be applied to larger haloes, as long as the mass ratio stays the same. For |$r\ll a$|⁠, the density slope is equal, and the density |$\rho \approx M_\text{h}/2\pi a^3\tilde{r}$|⁠, with |$\tilde{r} = r/a$|⁠. The radial velocity dispersion |$\bar{v_\mathrm{ r}^2}$| is also equivalent in the same regime up to a scaling factor with |$\bar{v_\mathrm{ r}^2} \approx GM_\text{h}/\tilde{r}\ln ({1/\tilde{r}})$| (Hernquist 1990). The internal dynamics of larger and smaller haloes are thus equivalent up to scaling factors.

We validate that the rate of growth of the central BH is adiabatic by growing a |$10^3\,$||$\text{M}_\odot$| BH using different rates-of-growth in a |$10^4\,\text{M}_\odot$| halo, varying from 250 to 4000 |$\text{M}_\odot$| Gyr⁻¹. The resulting spikes are indistinguishable above the convergence radius for rates of growth of 2000 |$\text{M}_\odot$| Gyr⁻¹ and below. The highest absolute and relative rate-of-growths of the simulations are below this reference value. We also compute the orbital times of particles in a Hernquist halo in Appendix C1, and verify that those within the G&S spike radius have orbital times below the maximum simulation time of 4 Gyr. We thus conclude that the assumption of adiabatic growth is valid.

The simulations were performed in seven different runs, each with a different halo mass and final BH mass, grown over a period of 4 Gyr. Every 0.1 Gyr, the system is recorded, and as the growth is adiabatic, each of these snapshots is a result. This creates a large parameter space of 235 different values of |$\mu$|⁠, see also Fig. 1. The simulated systems were chosen such that the largest possible range of |$\mu$| is probed and are given in Table 1. Some lower mass ratios did not produce measurable spikes, as these manifest below |$r_\text{conv}$|⁠. Due to the BHs growing from near-zero mass, low values of |$\mu$| are overrepresented in the final data set.

The final particle locations are radially binned from the location of the BH and are assumed to be Poisson-distributed in said bins. The error of such a distribution are given by |$\sqrt{N_\text{bin}}$| and as more than |$10^4$| particles are present in the least populated bins, these errors are relatively very small. These histograms are then fitted with a least-squares procedure using Minuit through the pyroot package (James & Roos 1975; Brun et al. 2019). The lower fitting boundary was taken to be the radius of convergence |$r_\text{conv}$| and the upper was chosen such that the fitting stops below the small numerical shock waves discussed above.

4 RESULTS

Our simulations produce measurable DM spikes. An example of these spikes can be found in Fig. 2a, where the the normalized density is fitted as a function of radius. Also visible is a numerically calculated G&S spike, which shows a clear deviation. In total, 235 unique systems with spikes were recorded, demonstrating the existence of this phenomenon in numerical simulations. This large data set allows us to fit a DM spike profile, where normalized parameters were found to be only dependent on the mass ratio |$\mu$|⁠. We first discuss our proposed profile, compare it to used approximations in the field, and then make an estimation of the effect of these results on the dephasing GWs due to DM spikes.

4.1 Spike profile

We empirically propose the final profile of a Hernquist halo with a spike after central adiabatic BH growth to be

where |$\tilde{\rho }_\text{final}$| is the final density normalized by |$\rho _\text{Hernq}$|⁠, |$\tilde{r}$| is radius normalized by a, |$r_\text{sp}$| is the spike radius, |$\beta$| governs the depletion of the original Hernquist profile as particles are pulled inward, and |$\gamma _\text{sp}$| is the slope of the spike. This profile is equivalent to (modified) G&S when |$\beta$| is set to one.

We check whether this profile can be fitted within our limited fitting range in Appendix C2, where we fit artificial data with Gaussian noise, both over a large and a limited range. While the fitting uncertainties do increase as the fitting range decreases, we observe no systematic differences between the two fits. The fitting uncertainties for |$\beta$| and |$r_\text{sp}$| are stable within an order of magnitude but the uncertainty on |$\gamma _\text{sp}$| grows significantly. We thus expect the fits of the N-body data to yield precise values for |$\beta$| and |$r_\text{sp}$|⁠, while |$\gamma _\text{sp}$| will be harder to determine.

We present the fits of all systems in Figs 2b, 2c, and 2d as a function of the mass-ratio |$\mu$|⁠. Both |$\beta$| and |$r_\text{sp}$| show a clear dependence on |$\mu$|⁠, and we thus perform further fits to determine their functional forms. The halo depletion parameter |$\beta$| is fitted to a simple power law

with the errors given in 1|$\sigma$|⁠, and a goodness-of-fit of |${\chi ^2_\text{red} = 1.81\times 10^{-5}}$|⁠. The data-points and resulting fit are given in Fig. 2b. For |$\mu = 0$|⁠, |$\beta$| is imposed to be 1. We note that for small |$\mu$| the relative depletion is approximately linear and independent of the halo mass.

The spike radius |$\tilde{r}_{\text{sp}}$| is fitted to a more complex function

with the errors given in 1|$\sigma$|⁠, and a goodness-of-fit of |${\chi ^2_\text{red} = 3.36 \times 10^{-6}}$|⁠. For systems where |$\mu \gtrapprox 0.1$|⁠, this relation simplifies to |$\tilde{r}_\text{sp} = 0.8\sqrt{\mu }$|⁠. The data-points and resulting fit are given in Fig. 2c. Also included are the spike radius as found by G&S (equation 4) and the Modified G&S profile (equation 7). The root mean squared error (RMSE) between every profile and the data is calculated and given in Table 2. The RMSE is lowest for our profile by two orders of magnitude, showing a very clear improvement over previous models for the spike radius.

The disagreement between the data and the theoretical profiles appears both at high values of |$\mu$| and at low values of |$\mu$|⁠. At high values of |$\mu$|⁠, we find |$\tilde{r}_\text{sp}$| to be larger than predicted by G&S. We suspect that this is due to the (lack of) self-consistency in the G&S formalism compared to the original formalism derived by Young (1980), which included multiple iterations to allow for changes in the density distribution (and thus the gravitational potential) of the stellar distribution during adiabatic growth. We suspect that this feedback effect is largest at large values of |$\mu$|⁠, as the pull of the DM halo is significant. At low values of |$\mu$|⁠, the data shows smaller spike radii than predicted by G&S and Modified G&S. However, we refrain from hypothesising about the cause as this could be the result of low resolutions of these small spikes.

Our simulations did not find a singular value for the spike slope as predicted, with the values for large |$\mu$| agreeing with the G&S value of |$\gamma _\text{sp} = 7/3$|⁠. For lower values of |$\mu$|⁠, this spike slope becomes smaller. We suspect that this is due to the relatively large fitting uncertainty (as discussed in Appendix C2). The data points are given in Fig. 2d. As |$\mu$| increases, the error of the data decreases and their values converges to the value as predicted by G&S. This is most likely due to the spike growing, and as a result being more visible in the data, resulting in a better fit. Concluding, for |$\mu \gtrapprox 0.06$|⁠, we find that |$\gamma _\text{sp}$| is in agreement with the value found by G&S. For values below this, we find lower values, likely due to limited fitting range as discussed in Appendix C2. More work needs to be done to show whether the low |$\mu$| values will also converge to |$\gamma _\text{sp} = 7/3$|⁠, or if the true value of |$\gamma _\text{sp}$| is indeed smaller.

Substituting equations (14) and (15) into our spike profile (equation 12) and integrating this profile over space for |$\gamma _\text{sp} = -7/3$| demonstrates mass conservation with deviations of less than 5 per cent for values of |$\mu \le 0.05$| and less than 10 per cent for values of |$\mu \le 0.1$|⁠.

4.2 Implications for future detections

We demonstrate the impact of our proposed spike profile (equation 12) compared to currently used approximations by simulating BH inspirals for both the proposed profile and Modified G&S. The feedback on the inspiralling object differs between different halo profiles, generating different dephasings in the resulting GWs. This has been done using the HaloFeedback code (Kavanagh et al. 2020).⁹ As |$\gamma _\text{sp}$| is not uniquely determined, we compute inspirals in our proposed profile for both |$\gamma _\text{sp} = 2$| and |$7/3$|⁠.

We can quantify the dephasing by calculating the number of cycles that the secondary object travels between two times as

where we fix |$t_f = 0$| as the time that the secondary object reaches |$r = 4 R_\mathrm{ S}$|⁠. We then fix |$t_i$| as five years before |$t_f$|⁠, regardless of environment. This allows us to make a comparison between a vacuum inspiral and an embedded inspiral.

In Table 3, the difference in the number of cycles |$\Delta N_\text{cycles}$| is given for a toy system with a central BH |$M_\text{1} = 10^3 \,\text{M}_\odot$| and an inspiralling object of |$M_\text{2} = \text{M}_\odot$|⁠. For |$M_\text{h} = 10^4 \,\text{M}_\odot$| we see an increase of |$\Delta N_\text{cycles}$| of a factor 2 when our proposed profile is used compared to Modified G&S, with a |$\gamma _\text{sp} = 7/3$|⁠. However, when the slope is dulled to |$\gamma _\text{sp} = 2$|⁠, |$\Delta N_\text{cycles}$| decreases until it is practically naught. We note that for the corresponding value of |$\mu$|⁠, our simulation found |$\gamma _\text{sp} \approx 7/3$|⁠.

If the halo mass is increased to |$M_\text{h} = 10^5 \text{M}_\odot$|⁠, our profile yields slightly lower results of about thirty per cent compared to Modified G&S for |$\gamma _\text{sp} = 7/3$|⁠. This is caused by the slightly smaller value of |$r_\text{sp}$| at low |$\mu$|⁠, which in turn lowers the normalization of the spike density. This is also the |$\mu$| regime where our simulations indicate a potentially lower value for the slope: |$\gamma _\text{sp} = 2$| again shows |$\Delta N_\text{cycles}$| nearly disappear. More shallow spikes thus yield a significant decrease of the dephasing of GWs compared to often used approximations. If |$\gamma _\text{sp}$| is instead |$7/3$| as is predicted by theory, then the dephasing effect is decreased by thirty per cent, a smaller yet still significant effect.

5 CONCLUSIONS

The mass distribution of DM around BHs is of vital importance to some astrophysical DM detection methods. We demonstrated the formation of DM spikes in fully numerical N-body simulations using realistic Hernquist halo profiles. The constant inner slope of |$\rho \propto r^{-1}$| means our results should be equivalent for other profiles such as the NFW, especially at low mass ratio |$\mu$|⁠. We proposed a new profile for these DM spikes, which deviates from theoretical predictions and common approximations. This profile only depends on a single new parameter compared to the initial halo profile, the mass ratio of BH and the total system |$\mu$|⁠. We observe a depletion of the original Hernquist of up to 20 per cent, and show the spike radius follows a different dependence on |$\mu$|⁠. This impacts the overall normalization of the spike compared to previous assumptions. We found a value of |$\gamma _\text{sp} = 7/3$| for large values of |$\mu$|⁠, but were unable to confidently determine it for smaller mass ratios due to the nature of the simulations. The impact of these results on future DM detections has been demonstrated using the expected dephasing of GWs due to the presence of a DM halo. This dephasing was shown to deviate between our profile and previous approximations, increasingly so as the spike slope and mass ratio |$\mu$| increases.

Our study unlocks the potential for future investigations into the exact shape of DM spikes and predictions for future DM detections. An improved resolution would yield smaller errors on the fitted parameters of the proposed profile, especially the slope, and allow for lower |$\mu$| systems to be probed. Simply increasing the number of particles is not feasible however, due to the scaling nature of N-body simulations. If achieved however, an exact value for the spike slope can be determined and larger mass differences in the system can be probed, bringing us closer to realistic DM spike profiles.

ACKNOWLEDGEMENTS

JLK offers their sincerest gratitude to Camila Correa for the plentiful and insightful discussions on the code that was developed for this work. The insightful comments of the anonymous peer reviewer were instrumental to solving numerical issues and improving the quality of the manuscript, and both authors wish to express their sincerest gratitude for them. The authors furthermore express their sincerest gratitude to Tim Linden for insightful feedback on an early draft of this manuscript, and to Gerben Wierda for crucial technological support. JLK acknowledges SURF and the Snellius supercomputer on which part of the research was conducted. The authors acknowledge the work of the swift – and other collaborations, providing us with open-source software that is vital to scientific research, specifically the MatPlotLib, Numpy, and SciPy python packages.

JLK performed the N-body simulations underlying this article and the resulting data was jointly processed by both authors. JLK analysed the results with critical input from RW. RW performed the analysis of the dynamical time-scales of the systems, the cross-reference of the data using the numerical G&S formalism, and the inspiral analysis. JLK wrote the manuscript with critical input from RW. No AI was used during any part of the research underlying, and writing of, this work.

DATA AVAILABILITY

The code underlying this article is a modified version of swift and can be found at https://github.com/JLKamermans/SWIFT-DAB. The data is reproducible using this code and will be shared upon reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: SOLUTIONS TO THE EDDINGTON EQUATION AND INITIAL CONDITIONS

The solution to the Eddington Equation for an isolated Hernquist halo is given by (Binney & Tremaine 2008)

where |$q = \sqrt{a\mathcal {E}/GM_\text{h}}$|⁠. It is also possible to determine |$f(\mathcal {E})$| with a massive BH already present. After the extra term |$-GM_\text{BH}/r$| is added to |$\psi$|⁠, the analytical solution involves a collection of elliptic integrals (Baes, Dejonghe & Buyle 2005), and a more manageable numerical result is preferred. Calculating the double derivative of |$\rho$| w.r.t. |$\psi$| yields

Rewriting r in terms of |$\psi$| yields

Equation (A4) is to be substituted into equation (A3), which is then inserted into equation (8) for numerical integration. In the limit with no BH, |$\mu = 0$|⁠, and the result of equation (A1) is obtained. The haloes generated using these initial conditions are found to be largely stable over a time span of 4 Gyr, with only a small numerical shockwave present that quickly travels radially outward. This shockwave is discussed further in Appendix C3. In Fig. A1, two Hernquist haloes are shown after 4 Gyr: one without a central BH and one with a central BH of constant mass equal to the initial mass as utilised in this work. Both show full equivalence to the initial Hernquist profile until the radius of convergence, below which a core forms due to numerical effects. The isolated Hernquist halo has also been simulated using the unmodified swift code upon which our modifications were built and an identical result was obtained, which is also included in Fig. A1.

APPENDIX B: DETERMINATION OF THE CONVERGENCE RADIUS

Determining the convergence radius |$r_\text{conv}$| is delicate due to the nature of our simulations. Our results are of an exponential shape and if |$r_\text{conv}$| is chosen too large, no significant result is fitted at all. In order to determine |$r_\text{conv}$|⁠, we fit data for |$r_\text{conv} = \delta \epsilon$| where |$1 \le \delta \le 3$|⁠, in steps of 0.5, and determine at which values of |$\delta$| the results remain consistent. A simulation of a |$10^4 \,\text{M}_\odot$| halo with a |$10^3\, \text{M}_\odot$| BH grown over 4 Gyr has been processed for these values of |$\delta$| and all found values for the proposed profile of equation (12) are compared, together with the fits’ values of |$\chi ^2_{\mathrm{ red}}$|⁠. See Figs B1, B2, and B3. The results for |$\beta$| are excluded, as this parameter remains consistent over all values of |$\delta$|⁠. We note that for some values, a spike was only found at later values of |$\mu$|⁠. The upper boundary of the fitting was set the same for every run and masked away the initial numerical shock wave (see Section C3).

We find that the results for |$\delta = 2\ \mathrm{ and}\ 2.5$| are consistent with each other, with best found values of |$r_\text{sp}$| and |$\gamma _\text{sp}$| falling comfortably between each others |$2\sigma$|-error margins. As we note the similar values of |$\chi ^2_\text{red}$| across these runs, we deem this to be a significant result. We have chosen to set |$\delta = 2.5$| for the remainder of this article.

APPENDIX C: VALIDATION OF RESULTS

C1 Dynamical time-scales

The adiabatic assumption is easily checked by comparing the orbital times of the particles in the halo with the total growth time of the central BH. By initializing a halo in a similar way to Bertone et al. (2024), the orbital times of a Hernquist profile has been computed as a function of initial radius |$r_i$|⁠, see also Fig. C1. This clearly shows that all particles with |$r_i < a/2$| have orbital times much smaller than 4 Gyr, indicating adiabatic growth for the initial Hernquist. As the spike grows and central densities increase, these orbital times are expected to only decrease, and adiabatic growth thus holds for the entire simulation within the region of spike growth.

C2 Influence of resolution on fitted values

In order to gauge the influence of the resolution of our simulations on the fitting of the parameters, two artificial systems were generated using similar binning size and Gaussian noise compared to the simulated data.¹⁰ The two artificial systems are chosen such that they represent a system with both a low |$(\mu \approx 0.02)$| and high (⁠|$\mu \approx 0.1$|⁠) mass central BH. These systems were fitted twice. One using a large range of |$10^{-9} \le r \le 10^3$|⁠, as often used in analytical studies (hereafter referred to as the full fit), and second on a small range of |$10^{-3} \le r \le 10^0$| (hereafter referred to as the zoomed fit), comparable to our fitting range. Unique gaussian noise of 1  per cent is applied and fits are performed a thousand times. The average fitted values are given in Tables C1 and C2, where we can see that there is no significant deviation between the zoomed and the full fits.

The average standard deviations are given in Tables C3 and C4. The size of the standard deviations is smaller than that of the N-body data, as the artificial data follows the fitted profile perfectly except for noise. Comparing the different results for the low mass BH in Table C3, we can see that the uncertainty on |$\gamma _\text{sp}$| improves by two orders of magnitude when going from the zoomed to the full fit. The improvements for |$r_\text{sp}$| and |$\beta$| are more modest, indicating that these two parameters are already quite accurately fitted with a more limited range. Similar observations can be made for the high mass BH in Table C4, where the uncertainties on |$r_\text{sp}$| and |$\beta$| are similar for the zoomed and full fits. The uncertainty on |$\gamma _\text{sp}$| improves by one order of magnitude, indicating that the zoomed fit has a less dramatic impact on the fitting uncertainty than for low-mass systems.

Concluding, we expect the fits of our systems to be very accurate for the |$r_\text{sp}$| and |$\beta$| parameters, with the fitting uncertainty being smaller as the mass ratio |$\mu$| of the system increases. However, the resolution of these simulations is too low to determine |$\gamma _\text{sp}$| with the same level of accuracy, and the effects of the fitting should be taken into account when interpreting the results.

One interesting quantity would be the necessary improvements to our N-body simulations to improve the determination of |$\gamma _\text{sp}$| for the low |$\mu$| system. If said system is fitted from |$10^{-4.5} \le r \le 10^0$|⁠, an improvement factor of 11 is reached for |$\gamma _\text{sp}$|⁠. Using equation (9) to calculate the particles necessary for such a convergence radius in a |$M_\text{h} = 10^{4}\, \text{M}_\odot$| system, we find approximately |$5\times 10^{10}$| particles are needed. Increasing the number of particles also decreases the error on the data points and might increase prefactor |$\alpha$| of equation (9), and this should thus only be seen as a ways of decreasing the fitting uncertainties. Nevertheless, these are 25 000 times more particles than currently present in the simulations of said systems, and since the calculations performed by swift scale as |$\mathcal {O}(N\log N)$|⁠, increasing the number of particles by such an amount is fully unfeasible without significant improvements to the efficiency of the code.

C3 Numerical shock waves

An initial shock wave of particles is launched radially outward in all simulations. It is created immediately at the start of the simulation. As shown in Appendix  A, this initial shock wave is present even if no BH is present or if the unmodified swift is used. This effect is also tested to be independent of the rate-of-growth of the central BH, the timestepping and the choice of gravitational softening. The velocity distributions of the particles are checked and behave as analytically expected. We therefore deem this effect to be due to slight imperfection in the initial conditions. We do not deem it likely this effect is due to the central BHs close range influence, as it is independent of gravitational softening. While the now missing particles will slightly change the potential of the halo, their absence is quickly overshadowed by the increasingly massive BH. Furthermore, after this shockwave has passed, the particle distribution remains unchanged, as demonstrated by the 4 Gyr old isolated haloes fully agreeing with the initial density distribution in Fig. A1. This shock wave is quickly traveling, therefore easily masked away during fitting, and only disruptive for the determination of the depletion parameter |$\beta$| (see Section 4.1) for the first few time-steps. We thus deem this numerical effect to not be of significant negative influence to our results.