Issue 
A&A
Volume 690, October 2024



Article Number  A299  
Number of page(s)  13  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202451304  
Published online  18 October 2024 
Dynamical friction from selfinteracting dark matter
^{1}
UniversitätsSternwarte, Fakultät für Physik, LudwigMaximiliansUniversität München, Scheinerstr. 1, D81679 München, Germany
^{2}
Excellence Cluster ORIGINS, Boltzmannstrasse 2, D85748 Garching, Germany
^{3}
Institute for Theoretical Physics, Goethe University, 60438 Frankfurt am Main, Germany
Received:
28
June
2024
Accepted:
15
July
2024
Context. Merging compact objects such as binary black holes provide a promising probe for the physics of dark matter (DM). The gravitational waves emitted during inspiral potentially allow one to detect DM spikes around black holes. This is because the dynamical friction force experienced by the inspiralling black hole alters the orbital period and thus the gravitational wave signal.
Aims. The dynamical friction arising from DM can potentially differ from the collisionless case when DM is subject to selfinteractions. This paper aims to understand how selfinteractions impact dynamical friction.
Methods. To study the dynamical friction force, we use idealised Nbody simulations, where we include selfinteracting dark matter.
Results. We find that the dynamical friction force for inspiralling black holes would be typically enhanced by DM selfinteractions compared to a collisionless medium (ignoring differences in the DM density). At lower velocities below the sound speed, we find that the dynamical friction force can be reduced by the presence of selfinteractions.
Conclusions. DM selfinteractions have a significant effect on the dynamical friction for black hole mergers. Assuming the Chandrasekhar formula may underpredict the deceleration due to dynamical friction.
Key words: black hole physics / gravitational waves / methods: numerical / dark matter
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Since the first direct detection of gravitational waves (GWs) from merging stellarmass black holes (BHs) and neutron stars with the Laser Interferometer GravitationalWave Observatory (LIGO; Abbott et al. 2016, 2017a,b, 2019), 100 years after the formulation of Einstein’s theory of general relativity (Einstein 1916) and about 40 years of the first indirect evidence (Hulse & Taylor 1975; Taylor et al. 1976; Weisberg et al. 1981), GW astronomy has become an important and very promising field to probe fundamental physics. The GWs emitted by binary mergers do not only allow for testing general relativity at unprecedented levels, but also to probe other physical processes. This in particular includes the physics of neutron stars. Binary mergers of those objects could allow probing the equation of state of ultrahighdensity matter (e.g. Baiotti 2022). In the future, GW detectors will become sensitive to subHertz GW frequencies, thanks to spacebased GW observatories such as the Laser Interferometer Space Antenna (LISA; AmaroSeoane et al. 2017) and Taiji (Hu & Wu 2017). This will allow us to observe the evolution of intermediate and extrememass ratio inspirals (I/EMRIs), in which a small stellarmass compact object inspirals into a much larger intermediate/supermassive black hole, over long times and provide a sensitive probe of BH environments.
BHs might be ‘dressed’ with DM, in other words, surrounded by a dense DM spike (e.g. Gondolo & Silk 1999; Bertone et al. 2005; Bertone 2024). The first robust evidence for an observation of a DM spike around a BH is claimed to be found by Chan & Lee (2024), see also the work by Alachkar et al. (2023). They studied the supermassive black hole binary system OJ 287 with outburst observations covering a time span of 104 years. The presence of a DM spike can impact the orbit of an inspiralling BH and make the GW signal sensitive to DM physics (e.g. Barausse et al. 2014). This could in principle be detectable, for example for M87* (Daghigh & Kunstatter 2024). The inspiralling BH is decelerated by the dynamical friction force, which depends on the density and the velocity dispersion of the DM spike (Chandrasekhar 1943). In consequence, the orbital decay of an inspiralling BH allows probing properties of the DM spike and thus the particle nature of DM.
Potentially, such DM spikes could lead to a faster decay of the orbit at small separations (≈1 pc). This would provide an explanation for the final parsec problem, that is, how binary BHs can lose enough energy and angular momentum to reach separations small enough that the energy loss due to the emission of GWs is sufficient to make them merge within a sufficiently small timescale (e.g. Milosavljević & Merritt 2003). In this context, AlonsoÁlvarez et al. (2024) claim that collisionless DM is not sufficient to make the BH orbit decay fast enough because the DM spike would be destroyed by the friction energy injected from the BHs. Instead, to solve the final parsec problem, DM selfinteractions would be needed to recreate the spike on a short enough timescale.
Selfinteracting dark matter (SIDM) is a class of particle physics models assuming that DM has selfinteractions beyond gravity. Initially, SIDM was proposed to solve problems on galactic scales, namely reducing the abundance of Milky Way satellites and the density in the centre of DM haloes (Spergel & Steinhardt 2000). DM selfinteractions have been studied in systems covering many orders of magnitude in mass, ranging from dwarf galaxies to galaxy clusters (e.g. Vogelsberger et al. 2014; Fitts et al. 2019; Robertson et al. 2019; Zeng et al. 2022; Ragagnin et al. 2024; Sabarish et al. 2024). Observations of those systems have been used to constrain the strength of the selfinteractions (e.g. Sagunski et al. 2021; Despali et al. 2022; Kong et al. 2024; Yang et al. 2024; Zhang et al. 2024). For a review on SIDM, we refer the reader to Tulin & Yu (2018) and Adhikari et al. (2022).
Potentially, the range of systems used to constrain SIDM can be extended by merging binary BHs. Based on pulsar timing array observations (Agazie et al. 2023; Reardon et al. 2023; Antoniadis et al. 2023), it was claimed that a selfinteraction crosssection of order 0.5–5.0 cm^{2} g^{−1} would be needed to solve the final parsec problem (AlonsoÁlvarez et al. 2024). We note that this claim is made based on the assumption of the density spike being in the longmeanfreepath regime (see Sect. 2.3). Interestingly, those binary systems could be a very sensitive probe for DM scattering, as we discuss in this paper.
When DM is selfinteracting, the density profile of the DM spike may change compared to collisionless matter. For SIDM, the slope of the density profile of the DM spike could be smaller than in the case of collisionless DM (Shapiro & Paschalidis 2014). The lower density would allow relaxing constraints from DM annihilation (Alvarez & Yu 2021). Importantly, this also affects the dynamical friction, as it depends on the density of the DM spike, and thus also the GW signal (e.g. Macedo et al. 2013; Yue & Han 2018; Becker et al. 2022). However, given the properties of the density spike, one may still need to calibrate the strength of the dynamical friction with simulations. Kavanagh et al. (2020, 2024) run Nbody simulations of collisionless DM to measure the Coulomb logarithm for a BH orbiting in a DM spike.
Several studies of dynamical friction for various DM models have been conducted. This includes ultralight DM (Wang & Easther 2022; Boey et al. 2024), also with selfinteractions (Glennon et al. 2024), selfinteracting scalar DM (Boudon et al. 2022, 2023, 2024; Kadota et al. 2024), superfluid DM (Berezhiani et al. 2024), and fuzzy DM (Lancaster et al. 2020). In this paper, we investigate the strength of the dynamical friction for a medium consisting of SIDM.
Often the binary BHs have been described in the Newtonian limit, however, depending on their evolution stage this may not be valid. Studies beyond Newtonian gravity have been undertaken by various authors (e.g. Traykova et al. 2023; Montalvo et al. 2024; Mukherjee et al. 2024). For this study, we focus on the Newtonian regime, that is, when the separation of the two BHs is sufficiently large. This allows us to build upon earlier work and simplifies the numerical description compared to the relativistic case. Extending our study to the relativistic case is left for future work.
The rest of this paper is structured as follows: Sect. 2 deals with the analytic description of dynamical friction and discusses it in the context of a BH moving through a DM spike. In Sect. 3, we describe the numerical setup that we use to study dynamical friction in the context of DM selfinteractions. This includes the simulation code and the initial conditions (ICs) and parameters for the simulations. It follows the analysis of our simulations and the results (Sect. 4). In Sect. 5, we discuss the limitations of this work as well as directions for further research. Finally, we summarise and conclude in Sect. 6. Additional information is provided in the appendix.
2. Analytic estimates
In this section, we first describe the properties of DM spikes around BHs relevant to our study. Secondly, we review analytic estimates of the dynamical friction force for collisionless and gaseous media. Last, we discuss how to describe the strength of selfinteractions and their impact on the dynamical friction arising from DM spikes acting on inspiralling BHs.
2.1. DM spike
In the vicinity of a BH the DM density can potentially reach high values due to the gravitational influence of the BH. Such a DM spike is commonly parameterised by a power law,
$$\begin{array}{c}\hfill {\rho}_{\mathrm{DM}}(r)={\rho}_{6}{\left(\frac{r}{{r}_{6}}\right)}^{{\alpha}_{\mathrm{sp}}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(1)
The radius r_{6} = 10^{−6} pc is a reference radius at which the density is given by ρ_{6}. We note that with this choice, we follow Coogan et al. (2022).
When taking relativistic effects into account (Sadeghian et al. 2013) the spike profile becomes,
$$\begin{array}{c}\hfill {\rho}_{\mathrm{DM}}(r)=\{\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\begin{array}{cc}\phantom{\rule{0.166667em}{0ex}}0\hfill & for\phantom{\rule{1em}{0ex}}r\le 2\phantom{\rule{0.166667em}{0ex}}{R}_{s}\hfill \\ \phantom{\rule{0.166667em}{0ex}}{\rho}_{\mathrm{sp}}{(1\frac{2{R}_{s}}{r})}^{3}{\left(\frac{r}{{r}_{\mathrm{sp}}}\right)}^{{\alpha}_{\mathrm{sp}}}\hfill & for\phantom{\rule{1em}{0ex}}2\phantom{\rule{0.166667em}{0ex}}{R}_{s}<r\le {r}_{\mathrm{sp}}\hfill \\ \phantom{\rule{0.166667em}{0ex}}\frac{{\rho}_{s}\phantom{\rule{0.166667em}{0ex}}{r}_{s}}{r}\hfill & for\phantom{\rule{1em}{0ex}}{r}_{\mathrm{sp}}<r\ll {r}_{s}\hfill \end{array}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(2)
For the relativistic corrections on the properties of the DM spike, we refer the reader to Eda et al. (2015), Tang et al. (2021), Capozziello et al. (2023), and John et al. (2024). In Fig. 1 we show the corresponding spike density and velocity dispersion as a function of radius. We note that the density profile of the DM around the BH could be even more complicated when considering a realistic formation scenario of the BH (Bertone et al. 2024).
Fig. 1. Density and onedimensional velocity dispersion for a DM density spike as a function of radius according to Eq. (2). We assume r_{sp} = 1 μpc, ρ_{sp} = 5.45 × 10^{−3} M_{⊙} μpc^{−3}, α_{sp} = 7/3, and M_{BH} = 10^{3} M_{⊙}. In addition, we indicate the Schwarzschild radius R_{s} = 2GM/c^{2}. For computing the velocity dispersion, we followed Kremer (2022). Here, we assumed an isotropic velocity distribution. In contrast to Kremer (2022), we solved their Eq. (27) including higherorder terms instead of Eq. (28). Moreover, we show the results for m/(k_{B}T)≈6.53 × 10^{3} s^{2} km^{−2}, with m being the rest mass of a physical DM particle, k_{B} the Boltzmann constant and T the absolute temperature of the system. 
As we see in the next subsection the Mach number of the secondary BH moving through the density spike plays a crucial role in determining the dynamical friction force acting on it. The Mach number, ℳ = v_{BH}/c_{s}, is given by the velocity, v_{BH}, of the BH and the sound speed, c_{s}. The local speed of sound depends on the adiabatic index γ, and is given by
$$\begin{array}{c}\hfill {c}_{\mathrm{s}}=\sqrt{\gamma}\phantom{\rule{0.166667em}{0ex}}{\sigma}_{\mathrm{DM}}.\end{array}$$(3)
To estimate c_{s}, we assume for simplicity that the local velocity distribution obeys a MaxwellBoltzmann distribution with a onedimensional velocity dispersion σ_{DM}. For the Mach number, ℳ, we make the assumption that the DM spike at sufficiently large radii is well described by Eq. (1).
Furthermore, we assume that the gravitational potential is dominated by the mass, M_{BH} of the primary (more massive) BH, and that the contribution of the DM spike and of the less massive inspiralling BH are negligible. Under these assumptions, the onedimensional velocity dispersion of the DM is
$$\begin{array}{c}\hfill {\sigma}_{\mathrm{DM}}=\sqrt{\frac{\mathrm{G}\phantom{\rule{0.166667em}{0ex}}{M}_{\mathrm{BH}}}{1+{\alpha}_{\mathrm{sp}}}\frac{1}{r}}\xb7\end{array}$$(4)
Assuming that the BH is moving on a circular orbit (${\mathit{v}}_{\mathrm{circ}}=\sqrt{\mathrm{G}\phantom{\rule{0.166667em}{0ex}}M(<r)/r}$), it follows that the Mach number is
$$\begin{array}{c}\hfill \mathcal{M}=\sqrt{\frac{1+{\alpha}_{\mathrm{sp}}}{\gamma}}\xb7\end{array}$$(5)
The adiabatic growth model for a DM spike assuming collisionless DM suggests, for an NFW halo, a spike index of α_{sp} = 7/3 (Gondolo & Silk 1999; Fields et al. 2014; Eda et al. 2015). For SIDM, Shapiro & Paschalidis (2014) derived the DM spike to be α_{sp} = 7/4. We note that this is based on the assumption that the selfinteraction crosssection has the velocity dependence σ ∝ v^{−4}. Moreover, Shapiro & Paschalidis (2014) assumed that the spike resides in the longmeanfreepath regime, in Sect. 2.3, we explain that this may eventually not always be the case. This implies that the motion of the inspiralling BH is supersonic when being on a circular orbit (see Eq. 5) given a spike index of α_{sp} = 7/3 or α_{sp} = 7/4. Even for a smaller value of α_{sp}, it is not feasible to reach the subsonic regime on a circular orbit within the range typically considered for the spike index (1 < α_{sp} < 3). However, a lower value for α_{sp} might be possible depending on the formation scenario (Ullio et al. 2001).
2.2. Dynamical friction
The dynamical friction force for a perturber moving through a collisionless medium of constant density was derived by Chandrasekhar (1943). The strongly collisional case with an object moving through a gaseous medium was studied analytically by several authors (e.g. Rephaeli & Salpeter 1980; Just et al. 1986; Just & Kegel 1990). Ostriker (1999) derived an analytic description for this regime. Following her, the dynamical friction force can be expressed as
$$\begin{array}{c}\hfill {F}_{\mathrm{df}}=4\pi \phantom{\rule{0.166667em}{0ex}}{\rho}_{\mathrm{DM}}{\left(\frac{\mathrm{G}\phantom{\rule{0.166667em}{0ex}}{m}_{\mathrm{BH}}}{{v}_{\mathrm{BH}}}\right)}^{2}\phantom{\rule{0.166667em}{0ex}}I.\end{array}$$(6)
Here, m_{BH}, is the mass of the perturber, in our case a BH. Its velocity is given by v_{BH} and the density of the perturbed medium, in our case DM, is denoted by ρ_{DM}. The factor I depends besides other things such as the Mach number, ℳ, on the collisional nature of the perturbed medium and can be specified for the different cases. In the collisionless case, I, is given by
$$\begin{array}{cc}\hfill {I}_{\mathrm{collisionless}}=& ln\left(\frac{{r}_{\mathrm{max}}}{{r}_{\mathrm{min}}}\right)[\mathrm{erf}(X)\frac{2X}{\sqrt{\pi}}\phantom{\rule{0.166667em}{0ex}}{e}^{{X}^{2}}]\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}with\phantom{\rule{1em}{0ex}}X\equiv \frac{{v}_{\mathrm{BH}}}{{\sigma}_{\mathrm{DM}}\sqrt{2}}\xb7\hfill \end{array}$$(7)
The effective size of the BH is indicated by r_{min}, and r_{max} denotes the effective size of the surrounding medium. Where, r_{min}, can be understood as the 90° deflection radius and r_{max} as the maximum impact parameter (Binney & Tremaine 2008). We note that while Chandrasekhar’s formula given by Eq. (7) approximates the dynamical friction well, corrections to the Coulomb logarithm, ln(r_{max}/r_{min}), can be necessary (Just & Peñarrubia 2005; Just et al. 2010), for example for velocity distributions deviating from a Maxwell–Boltzmann distribution. For corrections arising from fast moving DM particles on the Chandrasekhar formula, we refer the reader to the work by Dosopoulou (2023) on the dynamical friction in DM spikes.
For gaseous media, one distinguishes the subsonic case where the BH is moving with a velocity v that is smaller than the speed of sound c_{s} and the supersonic case where it is faster than c_{s}. In the first one, the Mach number, ℳ = v/c_{s} is smaller than unity and I is given as
$$\begin{array}{c}\hfill {I}_{\mathrm{subsonic}}=\frac{1}{2}ln\left(\frac{1+\mathcal{M}}{1\mathcal{M}}\right)\mathcal{M}.\end{array}$$(8)
We note that this equation is based on the assumption that (c_{s} − v_{BH}) t exceeds the effective size, r_{min}, of the perturber and that (c_{s} + v_{BH}) t is smaller than the effective size, r_{max}, of the surrounding medium. The time, t, indicates for how long the BH is perturbing the DM distribution. In the supersonic case with ℳ > 1, we have
$$\begin{array}{c}\hfill {I}_{\mathrm{supersonic}}=\frac{1}{2}\phantom{\rule{0.166667em}{0ex}}ln(1\frac{1}{{\mathcal{M}}^{2}})+ln\left(\frac{{v}_{\mathrm{BH}}\phantom{\rule{0.166667em}{0ex}}t}{{r}_{\mathrm{min}}}\right)\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(9)
This implicitly assumes (v_{BH} − c_{s}) t > r_{min} and (v_{BH} + c_{s}) t < r_{max}.
In Fig. 2, we plot the drag force computed according to Ostriker (1999) for our test setup used in Sect. 4. We can see that the dynamical friction for velocities slightly above the sound speed is becoming fairly strong in the gaseous case and sharply drops when the velocity falls below the speed of sound.
Fig. 2. Dynamical friction force as a function of the Mach number, i.e. the velocity in units of the sound speed. The dashed lines indicated the dynamical friction for a collisionless medium according to the Chandrasekhar formula (Chandrasekhar 1943) and the solid lines are for a gaseous medium as given by Ostriker (1999). We note that for the plot, we assume r_{max} = ℳ c_{s} t and ln(r_{max}/r_{min}) = ln(ℳ)+C. The line colours correspond to different values of C ∈ {4.0, 6.0, 8.0, 10.0, 12.0}. In addition, the red dotted lines indicate the velocity of a secondary BH moving on a circular orbit for different adiabatic indices. 
In the previous subsection, we have seen that the BH would typically move with a supersonic velocity. In consequence, the inspiralling BH experiences enhanced dynamical friction if the DM is not collisionless. However, if the BH is on an eccentric orbit, it might partially, when it is close to its apocentre, move with a subsonic velocity. We also note that dynamical friction can contribute to a circularisation of the orbit (Becker et al. 2022; Karydas et al. 2024).
In general, the dynamical friction acting on inspiralling BHs may not only arise from the DM spike but also from the accretion disk. The latter can also be described with the Ostriker formula (e.g. Szölgyén et al. 2022; Becker & Sagunski 2023).
2.3. The strength of the selfinteractions
A remaining question is how close a given crosssection would be to the collisionless or gaseous regime. This means we need to understand how much the trajectory of the DM particles is impacted by gravity compared to the scattering. A common measure for this in the context of SIDM haloes is the Knudsen number (e.g. Koda & Shapiro 2011). It is given by
$$\begin{array}{c}\hfill \mathrm{Kn}=\sqrt{\frac{4\pi \mathrm{G}}{\rho {\sigma}_{\mathrm{DM}}^{2}}}{\left(\frac{{\sigma}_{\mathrm{eff}}}{m}\right)}^{1}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(10)
This is comparing the ratio of the mean free path to the Jeans length. The latter is a measure for how susceptible the DM medium is to gravitational collapse. We note that the Knudsen number, Kn, is independent of the mass of the perturber, i.e. the BH mass. However, the strength of the dynamical friction depends on it. For Eq. (10) we used the effective crosssection per mass, σ_{eff}/m, for defining the Knudsen number. The effective crosssection (Yang & Yu 2022) allows one to map crosssections with different velocity dependencies onto each other and is given by
$$\begin{array}{c}\hfill {\sigma}_{\mathrm{eff}}=\frac{\u27e8{v}^{5}{\sigma}_{\mathrm{V}}(v)\u27e9}{\u27e8{v}^{5}\u27e9}\xb7\end{array}$$(11)
Where an average integrating over the distribution of relative velocities, v, and weighting the crosssection with v^{5} is computed assuming that the relative velocities follow a Maxwell–Boltzmann distribution. Furthermore, σ_{V} is the viscosity crosssection,
$$\begin{array}{c}\hfill {\sigma}_{\mathrm{V}}=3\pi {\displaystyle {\int}_{1}^{1}}\frac{\mathrm{d}\sigma}{\mathrm{d}\mathrm{\Omega}}{sin}^{2}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}cos\theta .\end{array}$$(12)
Here, we follow the normalisation of Yang & Yu (2022) implying that viscosity and total crosssection have equal values for isotropic scattering (σ_{V}_{iso} = σ_{tot}_{iso}). The viscosity crosssection has the advantage that it often provides a good matching between different angular dependencies (e.g. Yang & Yu 2022; Colquhoun et al. 2021) and is thus a good characteristic for the strength of the scattering. If the Knudsen number is large, i.e. Kn ≫ 1, a system is considered to be in the longmeanfreepath regime, i.e. gravity dominates over the scattering. In contrast, if Kn ≪ 1, the system would be in the shortmeanfreepath regime and scattering dominates over gravity. So we expect that the system is close to the collisionless case when the Knudsen number is larger, and for small numbers falls into the gaseous regime.
For the DM spike, the Knudsen number becomes smaller when getting closer to the primary BH as density and velocity dispersion increase. The Knudsen number scales as Kn ∝ r^{(1 + αsp)/2} when assuming a density spike profile as given by Eq. (1) and the gravitational potential being dominated by the BH. This implies that the behaviour of the dynamical friction could transition from the collisionless regime to the gaseous one while the secondary BH is inspiralling. However, it depends on the strength of the DM selfinteractions.
Last, we note that also in the case of a primary BH hosting a DM spike, the heat conduction due to selfinteractions being in the longmeanfreepath or shortmeanfreepath regime is independent of the mass of the primary BH if Eq. (10) holds. In other words, the Knudsen number does not depend on the mass of the primary BH. We note that this is in contradiction to the assumption by Shapiro & Paschalidis (2014) and Shapiro (2018). Instead of using the Jeans length for the Knudsen number, they approximated the gravitational scale height as the distance, r, to the primary BH. In fact, their Knudsen number is given as the mean free path divided by r. The influence of the BH changing the gravitational scale height can be motivated by studies of star clusters (e.g. Bettwieser & Spurzem 1986). It is unclear if this is in general the same for the SIDM case. If instead Eq. (10) holds, a DM spike given a crosssection within the typically studied range of ≈0.1–100 cm^{2}g^{−1} (e.g. Tulin & Yu 2018; Adhikari et al. 2022) could be to a large extent rather in the shortmeanfreepath regime than in the longmeanfreepath regime. This would be in line with the work by Pollack et al. (2015), where they studied a DM model with a fraction of the matter being ultrastrongly selfinteracting. Their heat conductivity does only depend on the density of the ultrastrongly interacting matter component.
3. Numerical setup
In this section, we first describe our simulation code with its implementation of DM selfinteractions. Next, it follows the setup for our simulations to understand the role of selfinteractions for dynamical friction.
3.1. Simulation code
We used the cosmological Nbody code OPENGADGET3 (see Groth et al. 2023, and the references therein). The SIDM module of the code has been described by Fischer et al. (2021a,b, 2022, 2024b). We employed it to model DM interactions beyond gravity. The code can simulate both large and smallangle scatterings. To the first one, we also refer to as rare scattering since a single scattering event on average leads to a larger energy and momentum transfer between the particles compared to smallangle scattering, which must be much more frequent to have a similar effect. For the first one, the angular dependence is modelled explicitly by drawing random angles with respect to the differential crosssection, for the second one an effective description is used based on a draglike behaviour (see also Kahlhoefer et al. 2014). The interactions between the numerical particles are modelled on a pairwise basis. We have computed the scattering probability for the rarely selfinteracting dark matter (rSIDM) and the drag force for the frequently selfinteracting dark matter (fSIDM) using a spline kernel (Monaghan & Lattanzio 1985). Its size is determined with the next neighbour approach using N_{ngb} = 64 neighbours. This choice is known to provide accurate results (Fischer et al. 2021a). We employed a separate time step criterion for SIDM to keep the interaction probability per interaction sufficiently low. The SIDM module is designed to explicitly conserve linear momentum and energy, even when executed in parallel.
To model the gaseous regime when selfinteractions are very strong, we used smoothedparticle hydrodynamics (SPH, Gingold & Monaghan 1977). Compared to GADGET2 (Springel 2005), an updated SPH formulation (Beck et al. 2016) with an updated viscosity scheme (Dolag et al. 2005; Beck et al. 2016) was employed. We chose a Wendland C^{6} kernel (Dehnen & Aly 2012) with its size being determined by the, N_{ngb} = 230, next neighbours.
3.2. Simulation setup
To test the dynamical friction, we simulated a BH travelling through a constant density consisting of DM. The mass density is set to ρ = 4.815 × 10^{−4} M_{⊙} μpc^{−3}. The velocities of the DM particles follow a MaxwellBoltzmann distribution, and their onedimensional velocity dispersion is chosen to be σ_{1D} = 625.5 km s^{−1}. The BH is represented by a single particle with a mass of m_{BH} = 7 × 10^{3} M_{⊙}. Initially, it is moving with a velocity of v_{ini} = 1037 km s^{−1}. These values are selected because they lead to a strong dynamical friction force, and we expect to find a significant impact in our simulations. Furthermore, we chose a cubic simulation volume with a side length of l_{box} = 829.5 μpc for most of our simulations but also used a box length of l_{box} = 663.6 μpc. We varied the box length to study its impact on the dynamical friction force. In addition, periodic boundary conditions were employed. We note that this also implies a periodic computation of the gravitational forces. Otherwise, the distribution of DM would easily collapse. The DM is resolved by N_{DM} = 5 × 10^{6} particles, which leads to a numerical particle mass of m_{DM} = 5.5 × 10^{−2} M_{⊙}. For testing the accuracy of our simulations, we also performed runs with N_{DM} = 1 × 10^{6} DM particles, i.e. a five times higher particle mass. We set the Plummerequivalent gravitational softening length to ϵ = 24 μpc. The gravitational forces are computed between all particles, i.e. we also take the selfgravity of the perturbed medium into account. In principle, the selfgravity can be relevant for the resulting dynamical friction (see Just & Kegel 1990). However, we do not separately study its contribution.
We note that the Nbody representation of the constant density is subject to fluctuation due to the sampling noise. The fluctuations can collapse under gravitational influence. We have checked that, for our choice of the box length and resolution, this does not occur within the time we simulate.
We ran simulations with a very anisotropic crosssection by employing the numerical scheme for frequent selfinteractions developed by Fischer et al. (2021a). Here the strength of the selfinteractions is specified in terms of the momentum transfer crosssection,
$$\begin{array}{c}\hfill {\sigma}_{\mathrm{T}}=2\pi {\displaystyle {\int}_{1}^{1}}\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathrm{\Omega}}_{\mathrm{cms}}}\left(1cos{\theta}_{\mathrm{cms}}\right)\mathrm{d}cos{\theta}_{\mathrm{cms}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(13)
We simulated velocityindependent elastic scattering with crosssections of σ_{V}/m ∈ {3 × 10^{−9}, 3 × 10^{−8}, 3 × 10^{−7}, 9 × 10^{−7}, 3 × 10^{−6}} (cm^{2} g^{−1})^{1}. In addition, we simulated velocityindependent isotropic crosssections with a total crosssection of σ_{V}/m ∈ {2 × 10^{−8}, 3 × 10^{−8}, 2 × 10^{−7}, 3 × 10^{−7}} (cm^{2} g^{−1}) to investigate the role of the angular dependence. We summarise all the crosssections used in Table 1.
Simulated crosssections.
In addition, we also investigated velocitydependent selfinteractions. Here, we used the implementation by Fischer et al. (2024b). We parameterise the velocity dependence of the viscosity crosssection as
$$\begin{array}{c}\hfill \frac{{\sigma}_{\mathrm{V}}}{m}=\frac{{\sigma}_{0}}{m}{(1+{\left(\frac{v}{w}\right)}^{2})}^{2}.\end{array}$$(14)
For the velocity parameter, w, we choose w ∈ {300, 461.1,683.2, 1483.5} (km s^{−1}). We note that the mean relative velocity of the DM particles in the ICs is ⟨v⟩ = 1411.6 km s^{−1}. For all simulations with a velocitydependent crosssection, we chose the value of σ_{0} to result in an effective crosssection (see Yang & Yu 2022) of σ_{eff}/m = 3 × 10^{−8} cm^{2} g^{−1}. A summary of all crosssections is given in Table 1.
We note that all the crosssections mentioned above are many orders of magnitude smaller than the tightest constraints on SIDM in the current literature (e.g. Harvey et al. 2019; Andrade et al. 2021; Sagunski et al. 2021; Eckert et al. 2022; Gopika & Desai 2023). Crosssections as large as the ones typically discussed in the astrophysical context imply a mean free path substantially shorter than the kernel size, and may no longer be faithfully simulated with the standard SIDM schemes (see our Appendix A as well as Fischer et al. 2024a). In addition, the simulations become more costly with increasing crosssection as the required timestep decreases. To model the regime of such a large crosssection, we use SPH. For the ICs, we set the initial velocity of the SPH particles to zero and chose their internal energy to correspond to the same velocity dispersion as in the SIDM setup.
4. Test problems
In this section, we present the results from our simulations following the setup described in Sect. 3. We begin with the collisionless case and investigate the impact of variations in the numerical setup. Subsequently, we study the role of the strength of the selfinteractions and its impact on the density wake. Furthermore, we investigate the matching of different angular and velocity dependencies onto each other in the case of dynamical friction.
4.1. The collisionless case
In our test setup as described in Sect. 3, the BH is travelling through an initially constant density. The BH perturbs the DM and creates a density wake, leading effectively to a force decelerating the BH. This dynamical friction force depends on the velocity of the BH relative to the velocity dispersion of the perturbed DM density (see Fig. 2). In our setup, we first expect the dynamical friction to become stronger while the BH is slowing down. After the BH velocity has decreased enough, it enters the phase where the dynamical friction becomes weaker with decreasing velocity.
In Fig. 3, we show the position and velocity of the BH in our test setup for collisionless matter. As expected, the velocity of the BH decreases over the course of the simulation. The results for two different box sizes are shown. We find that the dynamical friction is reduced for a smaller box size. This is not surprising, as r_{max} in the Coulomb logarithm of the Chandrasekhar formula, Eq. (7), depends for our setup effectively on the box size because the box size limits the maximal possible impact parameter. The results agree qualitatively with the Chandrasekhar formula, after an initial phase during which the density wake forms, implying that Eq. (7) provides a good description when the Coulomb logarithm is taken as a free parameter. For the large box, we show the results using two different resolutions. In both cases, the trajectory of the BH is very small, such that we can safely assume that the resolution for our study is large enough.
Fig. 3. Velocity and position of the BH in our test setup as a function of time. All simulations are for collisionless DM. The solid lines are from the setup with a box length of l = 829.5 μpc. In contrast, the dashed lines correspond to a setup with a smaller box size of l = 663.6 μpc. The highresolution results correspond to N_{DM} = 5 × 10^{6} DM particles, and the lowresolution ones to N_{DM} = 10^{6}. We note that in the smaller box setup, the low resolution implies a higher mass resolution than in the large box setup. 
We also quantitatively determined the Coulomb logarithm, ln(Λ) = ln(r_{max}/r_{min}), by fitting Λ to the simulation data for the run with the high resolution (N_{DM} = 5 × 10^{6}) and the large box size (l = 829.5 μpc). Here we did not use all of the simulation data but exclude the first phase of the simulation where the density wake forms and used the BH velocities for the fit only. The simulation results together with the fitted curves are displayed in Fig. 4. Our fit yielded a value of Λ_{fit} ≈ 6.4. Besides we can analytically estimate the expected value for Λ. It can be expressed as (Binney & Tremaine 2008)
Fig. 4. As in Fig. 3, position and velocity of the BH as a function of time. The highresolution setup with the large box size of l = 829.5 μpc is displayed. We fitted Chandrasekhar’s formula (Eqs. 6 and 7) to the simulation results to determine the Coulomb logarithm, ln(Λ), with Λ = r_{max}/r_{min}. Besides the simulation result (dotted lines), we also show the fitted curves (solid lines). 
$$\begin{array}{c}\hfill \mathrm{\Lambda}=\frac{{b}_{\mathrm{max}}\phantom{\rule{0.166667em}{0ex}}{v}_{\mathrm{typ}}^{2}}{\mathrm{G}\phantom{\rule{0.166667em}{0ex}}M}\xb7\end{array}$$(15)
We approximated the maximum impact parameter, b_{max}, with half the size of the box length, l. The typical velocity, v_{typ}, between the BH and the DM particles is approximately ≈10^{3} km s^{−1}. The mass, M, of the perturber might effectively be larger than the BH mass as DM particles are bound to the BH. We discuss this in Sect. 4.3. Here we can assume, M ≈ 2 × m_{BH}. Last, G, denotes the gravitational constant. Based on these assumptions we obtained a value of Λ ≈ 6.9, which agrees well with our simulation result (Λ_{fit} ≈ 6.4), especially given that the logarithm of that value matters for the dynamical friction.
4.2. Dependence on the strength of the selfinteractions
Given that the dynamical friction force differs between the collisionless and the gaseous case as described in Sect. 2, we expect that selfinteractions can have an impact on the strength of the dynamical friction force. To study this in detail, we run simulations with SIDM. In Fig. 5, we show the results for fSIDM with different values for the momentum transfer crosssection. We also include the collisionless and gaseous case. It is visible that the selfinteractions can enhance or reduce the dynamical friction, depending on the strength of the crosssection. For the gaseous regime, we expect the dynamical friction to be enhanced above the sound speed and reduced below. Accordingly, we find a larger velocity of the BH at the end of the simulation, when the BH velocity has fallen well below the sound speed of the DM.
Fig. 5. Similar to Fig. 3, position and velocity of the BH particle as a function of time. The results are shown for several interaction strengths of the DM, indicated by the opacity of the lines. All simulations are for the large box with the high resolution of N_{DM} = 5 × 10^{6} particles. For the simulation of gaseous matter, we used the same number of SPH particles. 
A detailed investigation of the dynamic friction force follows. But before, we take a look at the DM density perturbed by the BH. For the same simulations as shown in Fig. 5, we show the surface density of the DM in the top row of Fig. 6. Here, we can see that the shape and density of the wake depend on the strength of the selfinteractions. We find that for σ_{V}/m ∈ {3 × 10^{−8}, 3 × 10^{−7}} (cm^{2} g^{−1}) the densities in the wake are higher and at the same time the deceleration is stronger. For these small crosssections, the scattering leads to more lowvelocity particles in the density wake compared to CDM effectively enhancing the deceleration of the BH.
Fig. 6. Surface density (top row) and the velocity dispersion (bottom row) for our test simulations at a time of t = 0.49 yr. The panels display the different types of matter we investigate, with increasing collisionality from the left to the right. The density and one dimensional velocity dispersion are indicated following a logarithmic colour scale spanning a factor of two. For all panels of the same quantity, the same colour scale is used. In addition, we indicate the position of the BH particle with a red dot. We provide the time evolution of the surface density as a video available online. 
In contrast, for the cases of even stronger selfinteractions, we find that the densities in the wake are decreasing and the deceleration of the BH is weaker. In principle, this can be seen from Fig. 6, but we study this quantitatively in Sect. 4.3. We also note that the shape of the density wake changes for an increasing selfinteraction strength, impacting the dynamical friction force as well. The reduced dynamical friction for the strong crosssections is in line with the expectation of a weaker deceleration in the gaseous case for the subsonic regime.
In addition to the surface density, we show the velocity dispersion in the bottom row of Fig. 6. To compute the velocity dispersion, we first subtracted the bulk motion. In the gaseous case, we computed the velocity dispersion directly from the temperature of the SPH particles. Similar to the density plots, it is visible how the wake changes when selfinteractions are present. For example, the velocity dispersion in the vicinity of the perturber increases with crosssection and the shape of the wake changes.
For a more detailed analysis of the dynamical friction force and its dependence on the crosssection, we computed the force and show it as a function of the velocity (Fig. 7). The force is computed from the momentum of the BH and its change between consecutive snapshots. We note that the strongest SIDM crosssection shown is σ_{V}/m = 3 × 10^{−6} cm^{2} g^{−1}. Given our resolution, we are not able to faithfully simulate a crosssection as large as σ_{V}/m = 3 × 10^{−5}cm^{2}g^{−1} because the mean free path becomes much smaller than the size of the kernel used for the selfinteractions (see Appendix A and Fischer et al. 2024a).
Fig. 7. Dynamical friction force as a function of the Mach number for our simulations. This plot is analogue to Fig. 2, which gives the analytic description of the dynamical friction force. Different crosssections are shown, indicated by the Knudsen number. Here, Kn = ∞ corresponds to a collisionless medium and Kn = 0.0 to a gaseous medium. We note that in the simulation the BH initially starts at a high velocity, i.e. at the right side of the plot and slows down throughout the simulation, i.e. moves towards the left side of the plot. 
We specify the strength of the crosssection in terms of the Knudsen number (Eq. (10)). The collisionless case is indicated by Kn = ∞ and the gaseous case by Kn = 0.0. In Fig. 7 we can see that the dynamical friction is initially increasing (starting at the highest velocity) and decreases for later times of the simulation. At the beginning of the simulation, where the BH moves with a supersonic velocity, the dynamical friction is the smallest in the collisionless case. At later stages, it depends on the crosssection whether the selfinteractions enhance or reduce the deceleration of the BH. In Fig. 7, we plot the dynamical friction as a function of the velocity. This shows that the decelerating force reaches its maximum roughly at the sound speed, as expected. Interestingly, the selfinteracting run with an intermediate crosssection (Kn = 1.35 × 10^{−1}) shows the strongest dynamical friction, stronger than in the collisionless or gaseous regime. In principle, the figure is built analogously to Fig. 2, but a direct comparison is limited by the fact that our simulation setup may not fulfil the assumptions of the Ostriker formula for the dynamical friction in the gaseous case (see Sect. 2.2). This is for example illustrated by Fig. 3, where we can see that the exact results depend on the choice of the setup, in particular on the box size.
Overall, we find that scatterings of the DM particles can alter the strength of the dynamical friction in a nontrivial way. Importantly, we have shown that the Knudsen number (Eq. (10)) provides a good measure for the medium being collisionless or collisional and its impact on the dynamical friction arising from the selfinteractions.
4.3. Properties of the wake
To provide more insight into the increased dynamical friction for mildly collisional media, we computed further quantities such as the density of the wake, its velocity dispersion, and bulk motion as well as the DM mass gravitational bound to the BH.
In Fig. 8, we show the mean density within a sphere of radius r = 43.09 μpc about the point of the maximum DM density in the wake as a function of the velocity of the perturber. The simulations showing a strong dynamical friction force also show a larger density compared to the ones with relatively weak dynamical friction. This is not surprising, as a more massive wake can exert a stronger gravitational pull on the BH and thus faster decelerate it. The presence of selfinteractions allows for the growth of a larger wake density, as the scattering produces more lowvelocity particles compared to the collisionless case. However, for sufficiently large selfinteraction when entering the shortmeanfreepath regime, the wake density starts to decrease with increasing crosssection.
Fig. 8. Density of the wake, i.e. the mean density about the point with the highest density, as a function of the velocity of the perturbing BH. The simulation results of the different crosssections are indicated by their Knudsen number. 
Furthermore, we show the velocity dispersion in the density wake in Fig. 9. We use the same volume as for the density computation. The velocity dispersion was computed after subtracting the mean velocity of the particles. The velocity dispersion in the wake increases compared to the ICs and is larger for stronger DM selfinteractions. Interestingly, the crosssection with Kn = 1.35 × 10^{−1}, which gives the strongest dynamical friction, shows a substantially lower velocity dispersion compared to the larger crosssections.
Fig. 9. Onedimensional velocity dispersion of the wake, i.e. the velocity dispersion in the vicinity of the point with the highest density, as a function of the velocity of the perturber. We note that the velocity dispersion is computed after the bulk motion has been subtracted. The simulation results for different crosssections are shown, indicated by the Knudsen number. 
The bulk motion of the density wake is displayed in Fig. 10. The stronger the selfinteractions the more efficiently the BH drags the DM along its trajectory, except for the largest crosssection for which the bulk velocity is again smaller. This is in line with the strength of the dynamical friction force found in Fig. 7, implying that when the DM is dragged along by the BH, it can more efficiently decelerate the BH.
Fig. 10. Mean velocity of the wake in the direction of the motion of the BH, i.e. the mean velocity in the vicinity of the point with the highest density. The results for different crosssections are shown as a function of the velocity of the perturber. They are indicated by their corresponding Knudsen number. 
Last, we computed the DM mass that is gravitationally bound to the BH. To do so, we compared the relative velocity of a DM particle and the BH to the escape velocity,
$$\begin{array}{c}\hfill {v}_{\mathrm{esc}}=\sqrt{\frac{2\phantom{\rule{0.166667em}{0ex}}\mathrm{G}\phantom{\rule{0.166667em}{0ex}}{m}_{\mathrm{BH}}}{r}},\end{array}$$(16)
of the BH. Here, the perturbing BH has a mass of m_{BH} and the gravitational constant is given by G. The distance between the DM particle and the BH is denoted by r. This may underestimate the bound mass to some extent as the DM bound to the BH increases the effective mass and thus may lead to further particles being gravitationally bound.
Our estimate for the gravitational bound mass is shown in Fig. 11. It becomes visible that the bound mass depends on the crosssection. In particular, strong selfinteractions can effectively increase the bound mass. This can be understood by the effect that the scattering has on the velocity distribution. We note that the velocities in the wake do not follow a Maxwell–Boltzmann distribution, but the scattering evolves the distribution towards it. This implies that the selfinteractions create more lowvelocity particles compared to the collisionless case. In consequence, more DM particles are gravitationally bound to the BH. This effectively enhances the dynamical friction acting on the perturber given the crosssection is small enough such that the behaviour known from the gaseous regime (Ostriker 1999) does not become dominant. Although the dynamical friction decreases with increasing crosssection for strong selfinteractions, the gaseous regime nevertheless yields a stronger dynamical friction force in the supersonic regime compared to the collisionless case.
Fig. 11. Dark matter mass gravitational bound to the perturbing BH as a function of the velocity of the perturber. We indicate the different crosssections by their Knudsen number. 
4.4. Angular dependence
The angular dependence of the selfinteraction crosssection can in principle lead to qualitative differences in the evolution of astrophysical systems. This is for example well known for merging DM haloes (e.g. Fischer et al. 2023; Sabarish et al. 2024) or the abundance of satellites (Fischer et al. 2022). However, for systems that are rather relaxed, the angular dependence plays a minor role only and the effect of selfinteractions with different angular dependencies can be mapped onto each other using the viscosity crosssection σ_{V} (Eq. 12, see also Yang & Yu 2022).
To test whether this mapping also works for our setup to study dynamical friction, we ran simulations with a very anisotropic crosssection (using the fSIDM framework) and an isotropic crosssection. The crosssections are chosen to have the same value in terms of σ_{V}/m or σ_{T}/m to test the matching with the momentum transfer crosssection (see Eq. 13) as well.
In Fig. 12, we compare the velocity of the BH for the different angular dependencies. We computed the ratio of the change in BH velocities for simulations with a matched isotropic crosssection to the fiducial fSIDM simulations. The two panels in Fig. 12 display the ratio of the velocities as a function of time for a weaker fSIDM crosssection of σ_{V}/m = 3 × 10^{−8} cm^{2} g^{−1} (top panel) and a stronger one of σ_{V}/m = 3 × 10^{−7} cm^{2} g^{−1} (bottom panel). For a perfect matching, we would expect the ratio to be unity. As we can see, the viscosity crosssection (Eq. 12) provides a better match than the transfer crosssection (Eq. 13). This is the case for the weaker and stronger crosssection that we have simulated to test this.
Fig. 12. Ratio between the velocity change, Δv = v(0)−v(t), due to dynamical friction of an fSIDM crosssection and the transfer (blue) / viscosity (orange) matched isotropic crosssections. The upper panel is for an fSIDM viscosity crosssection of σ_{V}/m = 3 × 10^{−8} cm^{2} g^{−1} and the lower panel is for σ_{V}/m = 3 × 10^{−7} cm^{2} g^{−1}. 
4.5. Velocity dependence
The effective crosssection σ_{eff} (introduced by Yang & Yu 2022) allows one to map crosssections with different velocity dependencies onto each other. It integrates the crosssection with its velocity dependence over the assumed MaxwellBoltzmann distribution of the scattering velocities, applying a weighting factor of v^{5}. The effective crosssection is given by
$$\begin{array}{c}\hfill {\sigma}_{\mathrm{eff}}=\frac{\u27e8{v}^{5}{\sigma}_{\mathrm{V}}(v)\u27e9}{\u27e8{v}^{5}\u27e9}\xb7\end{array}$$(17)
We note that this is very similar to the matching procedure proposed by Yang et al. (2023); see also the work by Outmezguine et al. (2023). It has been shown that this works fairly well for isolated haloes and, in turn, velocitydependent crosssections can be approximately described with constant crosssections. However, it is known that this does not work for systems that involve multiple velocity scales, such as merging galaxy clusters (Sabarish et al. 2024) or haloes containing substructure (Fischer et al. 2024b).
Here, we investigate how well the matching works for our setup to study dynamical friction. In Fig. 13, we show simulations involving different velocity dependencies specified in terms of w, but they are all matched using Eq. (11). We have simulated very anisotropic crosssections (top panel) and isotropic ones (bottom panel). For times up to t ≈ 0.5 yr, the difference between the runs is fairly small and increases at later stages when the BH has slowed down. We find that the matching of the velocity dependence using the effective crosssection does not work as well as for the angular dependence using the viscosity crosssection. However, in our case where the velocity of the BH is falling below half its initial value over the course of the simulation, the deviation between the runs hardly exceeds 2.5%. This might be enough for a rough estimate, and separate simulations for velocitydependent cases may only be needed when a higher accuracy is required.
Fig. 13. Ratio of the velocity change, Δv = v(0)−v(t), due to dynamical friction acting on the BH for a velocitydependent and a velocityindependent (constant) crosssection. All crosssections were matched using the effective crosssection (Eq. 11), i.e. all crosssections correspond to σ_{eff}/m = 3 × 10^{−8} cm^{2} g^{−1}. The upper panel displays runs with frequent selfinteractions, and the lower panel is for isotropic scattering. 
5. Discussion
In this section, we discuss the limitations of our study, mainly the simulation setup. Moreover, we elaborate on the implications of our results and further directions to provide a better understanding of the dynamical friction in merging binary BHs.
Currently, there does not exist a theory for the nonequilibrium dynamics of SIDM to describe dynamical friction. This is why we have to use simulations to understand the impact of a mildly collisional medium on the dynamical friction force. Ideally, one would develop a more general formalism that contains both the Chandrasekhar and the Ostriker formula in the limits of a collisionless and a gaseous medium but also describes the intermediate regime as a function of the selfinteraction crosssection. But even then one may still depend on simulations to calibrate such a formalism and make it applicable to problems such as merging binary BHs. Moreover, we want to note that the Ostriker formula assumes the trajectory of the perturber to be a straight line. Kim & Kim (2007) improved on this and provided a generalisation to circular orbits.
For binary BHs, we find that the motion of the secondary BH through the DM spike would typically happen at Mach numbers larger than unity, in other words, in the supersonic regime. Moreover, the Mach number does only very slightly depend on the radius, as the mass of the DM spike in the relevant radial range is typically much smaller than the mass of the primary BH. Also, the spike index does not play a significant role in the motion of the secondary BH being supersonic. But more importantly, a nonzero eccentricity could lead to a subsonic phase about the apocentre of the orbit of the secondary BH. This implies that the supersonic regime where the dynamical friction is enhanced by selfinteractions plays an important role and would lead to a faster decay of the orbit compared to a collisionless medium with the same density.
To study the dynamical friction with selfinteractions being present, simulations are required. We used an idealised setup with a periodic box hosting initially a constant density. We have to note that this setup comes with drawbacks that limit the conclusion we can draw. In particular, it is difficult to create a stable setup where the box size is small enough such that the velocity dispersion of the DM counteracts the gravitational collapse. This effectively limits the maximal box size that can be studied, which also has an impact on the dynamical friction. This is because the box size is related to the maximum impact parameter r_{max} (see Eq. 7). Moreover, the setup does not always fulfil the assumptions of analytic formulas discussed in Sect. 2.2. Also, we cannot exclude that perturbations arising from the supersonic phase still affect the trajectory of the BH in the subsonic phase. Nevertheless, we are still able to make qualitative statements and gain new insights into the dynamical friction. Especially for this purpose, we tuned the setup to give a relatively larger impact on the BH velocity and make it easier to obtain qualitative differences being much larger than the numerical errors.
To obtain accurate quantitative predictions on the strength of the dynamical friction in BH mergers, one would need to simulate a setup much closer to the real system. Kavanagh et al. (2020, 2024) studied intermediate massratio inspirals resolving the DM spike around the primary BH and measured the deceleration due to the dynamical friction of the secondary. In principle, such a setup could be used for studies of SIDM as well. Potentially this does not only allow us to calibrate the Chandrasekhar formula, but also the Ostrikers formula for the gaseous regime and with this make it applicable for merging BHs. Ideally, one would have a more general formalism for arbitrary crosssections that can be calibrated as well.
In our simulations, we find a regime of intermediate crosssections where the dynamical friction is stronger than in the collisionless case or the gaseous case. Given the drawback our setup suffers from, it would be great to study this in a setup closer to the physical systems as described above. Simulations of BHs moving through a DM spike could potentially tell whether the intermediate regime of enhanced dynamical friction is only due to our setup or would also occur in real BH mergers.
The prediction of the exact strength of the dynamical friction force and its impact on the GW signal is difficult. We find that not only does the DM density matter but its collisionless or collisional nature, too. Moreover, the latter plays a crucial role in the spike index and the DM density (Shapiro & Paschalidis 2014). But it must also be noted that the energy that is transferred due to dynamical friction from the secondary BH into the DM spike affects the DM spike as well. AlonsoÁlvarez et al. (2024) argues that only SIDM can effectively transport this heat outwards and replenish the DM spike. In contrast to collisionless matter, this effectively allows a SIDM spike to make the BH orbit decay and solve the final parsec problem. Overall, it remains challenging to predict the exact GW signal, especially beyond the Newtonian regime.
6. Conclusions
In this paper, we have reviewed analytic estimates for the dynamical friction force in the collisionless and gaseous regime, and discussed the two regimes in the context of binary BH mergers. Furthermore, we have performed idealised simulations to study the dynamical friction as a function of the selfinteraction crosssection. In addition, we also studied how well crosssections with different angular and velocity dependencies can be matched onto each other. Our main results can be summarised as follows:

DM selfinteractions can enhance but also lower the dynamical friction force, depending on the velocity of the perturber relative to the local velocity dispersion of the DM.

For merging binary BHs, we expect the perturber to be in the supersonic regime all the time when being in a circular orbit and partially when following an eccentric orbit. In this regime, selfinteractions enhance the dynamical friction force.

The Knudsen number, Kn, is a good proxy for how much the selfinteractions impact the dynamical friction, in other words, how close the system is to the collisionless or gaseous regime.

When the SIDM crosssection is large enough, Kn ≲ 2, we expect the dynamical friction in the supersonic regime to be significantly larger than in the collisionless case. We note that for typical binary BHs such a crosssection could be many orders of magnitude smaller than current upper bounds on the selfinteraction crosssection (Adhikari et al. 2022).

This implies that observed BH mergers, when being consistent with collisionless matter, can potentially imply extremely tight upper bounds on the SIDM crosssection.

We find that the viscosity crosssection provides a good matching procedure between isotropic and very anisotropic crosssections when computing the effect of selfinteractions on the dynamical friction.

The matching of different velocity dependencies onto each other works less well than the angular matching. However, depending on the application, the deviation could still be within an acceptable range.
The results we obtained are solely based on simulations of a constant density perturbed by a BH in a periodic box. This allowed us to gain insight into the qualitative behaviour but does not allow for accurately quantifying how strong the effect of selfinteractions is on the dynamical friction force. In consequence, the next step would be to simulate a SIDM spike with an orbiting BH, similar to the work done by Kavanagh et al. (2020, 2024). Such a setup would allow us to precisely determine the strength of the dynamical friction force acting onto BHs orbiting in DM spikes when selfinteractions are present. This is the subject of a forthcoming paper (Sabarish et al., in prep.).
Data availability
Movie is available at https://www.aanda.org
To make it easier to compare crosssections of different angular dependencies, we do not specify the crosssection in terms of the momentum transfer crosssection but use the viscosity crosssection (see Eq. 12).
Acknowledgments
The authors thank the organisers of the Pollica 2023 SIDM Workshop, where part of this work was done. In addition, the authors thank all participants of the Darkium SIDM Journal Club for discussions. MSF and LS thank Gonzalo AlonsoÁlvarez, Niklas Becker, Andrew Benson, Akaxia Cruz, Xiaolong Du, Frederick Groth, Charlie Mace, Annika H. G. Peter, Rainer Spurzem, Sabarish Venkataramani, Daneng Yang, Shengqi Yang, HaiBo Yu, Xingyu Zhang, and Zhichao Carton Zeng as well as the members of the CAST group at the University Observatory Munich for helpful discussions. They are also grateful to Bradley Kavanagh for helpful comments on the draft. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC2094 ‘Origins’ – 390783311. Software: NumPy (Harris et al. 2020), Matplotlib (Hunter 2007).
References
 Abbott, B., Abbott, R., & Abbott, T., et al. 2016, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B., Abbott, R., Abbott, T., et al. 2017a, Phys. Rev. Lett., 119, 161101 [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJ, 848, L12 [Google Scholar]
 Abbott, B., Abbott, R., Abbott, T., et al. 2019, Phys. Rev. X, 9, 031040 [NASA ADS] [Google Scholar]
 Adhikari, S., Banerjee, A., Boddy, K. K., et al. 2022, ArXiv eprints [arXiv:2207.10638] [Google Scholar]
 Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023, ApJ, 951, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Alachkar, A., Ellis, J., & Fairbairn, M., 2023, Phys. Rev. D, 107, 103033 [NASA ADS] [CrossRef] [Google Scholar]
 AlonsoÁlvarez, G., Cline, J. M., & Dewar, C. 2024, Phys. Rev. Lett., 133, 021401 [CrossRef] [Google Scholar]
 Alvarez, G., & Yu, H. B. 2021, Phys. Rev. D, 104, 043013 [NASA ADS] [CrossRef] [Google Scholar]
 AmaroSeoane, P., Audley, H., Babak, S., et al. 2017, ArXiv eprints [arXiv:1702.00786] [Google Scholar]
 Andrade, K. E., Fuson, J., GadNasr, S., et al. 2021, MNRAS, 510, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadis, J., Arumugam, P., Arumugam, S., et al. 2023, A&A, 678, A50 [CrossRef] [EDP Sciences] [Google Scholar]
 Baiotti, L. 2022, Arab. J. Math., 11, 105 [CrossRef] [Google Scholar]
 Barausse, E., Cardoso, V., & Pani, P., 2014, Phys. Rev. D, 89, 104059 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110 [Google Scholar]
 Becker, N., & Sagunski, L. 2023, Phys. Rev. D, 107, 083003 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, N., Sagunski, L., Prinz, L., & Rastgoo, S. 2022, Phys. Rev. D, 105, 063029 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhiani, L., Cintia, G., De Luca, V., & Khoury, J. 2024, JCAP, 2024, 024 [CrossRef] [Google Scholar]
 Bertone, G. 2024, ArXiv eprints [arXiv:2404.11513] [Google Scholar]
 Bertone, G., Zentner, A. R., & Silk, J. 2005, Phys. Rev. D, 72, 103517 [NASA ADS] [CrossRef] [Google Scholar]
 Bertone, G., Wierda, A. R. A. C., Gaggero, D., et al. 2024, ArXiv eprints [arXiv:2404.08731] [Google Scholar]
 Bettwieser, E., & Spurzem, R. 1986, A&A, 161, 102 [NASA ADS] [Google Scholar]
 Binney, J., & Tremaine, S., 2008, Galactic Dynamics (Princeton: Princeton University Press) [Google Scholar]
 Boey, R., Wang, Y., Kendall, E., & Easther, R. 2024, Phys. Rev. D, 109, 103526 [NASA ADS] [CrossRef] [Google Scholar]
 Boudon, A., Brax, P., & Valageas, P. 2022, Phys. Rev. D, 106, 043507 [NASA ADS] [CrossRef] [Google Scholar]
 Boudon, A., Brax, P., & Valageas, P. 2023, Phys. Rev. D, 108, 103517 [NASA ADS] [CrossRef] [Google Scholar]
 Boudon, A., Brax, P., Valageas, P., & Wong, L. K. 2024, Phys. Rev. D, 109, 043504 [NASA ADS] [CrossRef] [Google Scholar]
 Capozziello, S., Zare, S., & Hassanabadi, H. 2023, ArXiv eprints [arXiv:2311.12896] [Google Scholar]
 Chan, M. H., & Lee, C. M. 2024, ApJ, 962, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
 Colquhoun, B., Heeba, S., Kahlhoefer, F., Sagunski, L., & Tulin, S. 2021, Phys. Rev. D, 103, 035006 [NASA ADS] [CrossRef] [Google Scholar]
 Coogan, A., Bertone, G., Gaggero, D., Kavanagh, B. J., & Nichols, D. A. 2022, Phys. Rev. D, 105, 043009 [NASA ADS] [CrossRef] [Google Scholar]
 Daghigh, R. G., & Kunstatter, G. 2024, ArXiv eprints [arXiv:2308.15682] [Google Scholar]
 Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068 [NASA ADS] [CrossRef] [Google Scholar]
 Despali, G., Walls, L. G., Vegetti, S., et al. 2022, MNRAS, 516, 4543 [NASA ADS] [CrossRef] [Google Scholar]
 Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Dosopoulou, F. 2023, ArXiv eprints [arXiv:2305.17281] [Google Scholar]
 Eckert, D., Ettori, S., Robertson, A., et al. 2022, A&A, 666, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eda, K., Itoh, Y., Kuroyanagi, S., & Silk, J. 2015, Phys. Rev. D, 91, 044045 [NASA ADS] [CrossRef] [Google Scholar]
 Einstein, A. 1916, Ann. Phys., 354, 769 [Google Scholar]
 Fields, B. D., Shapiro, S. L., & Shelton, J. 2014, Phys. Rev. Lett., 113, 151302 [NASA ADS] [CrossRef] [Google Scholar]
 Fischer, M. S., Brüggen, M., SchmidtHoberg, K., et al. 2021a, MNRAS, 505, 851 [NASA ADS] [CrossRef] [Google Scholar]
 Fischer, M. S., Brüggen, M., SchmidtHoberg, K., et al. 2021b, MNRAS, 510, 4080 [Google Scholar]
 Fischer, M. S., Brüggen, M., SchmidtHoberg, K., et al. 2022, MNRAS, 516, 1923 [NASA ADS] [CrossRef] [Google Scholar]
 Fischer, M. S., Durke, N.H., Hollingshausen, K., et al. 2023, MNRAS, 523, 5915 [CrossRef] [Google Scholar]
 Fischer, M. S., Dolag, K., & Yu, H.B. 2024a, A&A, 689, A300 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fischer, M. S., Kasselmann, L., Brüggen, M., et al. 2024b, MNRAS, 529, 2327 [NASA ADS] [CrossRef] [Google Scholar]
 Fitts, A., BoylanKolchin, M., Bozek, B., et al. 2019, MNRAS, 490, 962 [NASA ADS] [CrossRef] [Google Scholar]
 Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Glennon, N., Musoke, N., Nadler, E. O., PrescodWeinstein, C., & Wechsler, R. H. 2024, Phys. Rev. D, 109, 063501 [NASA ADS] [CrossRef] [Google Scholar]
 Gondolo, P., & Silk, J. 1999, Phys. Rev. Lett., 83, 1719 [Google Scholar]
 Gopika, K., & Desai, S. 2023, Phys. Dark Univ., 42, 101291 [NASA ADS] [CrossRef] [Google Scholar]
 Groth, F., Steinwandel, U. P., Valentini, M., & Dolag, K. 2023, MNRAS, 526, 616 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Harvey, D., Robertson, A., Massey, R., & McCarthy, I. G. 2019, MNRAS, 488, 1572 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W.R., & Wu, Y.L. 2017, Natl. Sci. Rev., 4, 685 [CrossRef] [Google Scholar]
 Hulse, R. A., & Taylor, J. H. 1975, ApJ, 195, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 John, I., Leane, R. K., & Linden, T. 2024, Phys. Rev. D, 109, 123041 [NASA ADS] [CrossRef] [Google Scholar]
 Just, A., & Kegel, W. H. 1990, A&A, 232, 447 [NASA ADS] [Google Scholar]
 Just, A., & Peñarrubia, J. 2005, A&A, 431, 861 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Just, A., Kegel, W. H., & Deiss, B. M. 1986, A&A, 164, 337 [NASA ADS] [Google Scholar]
 Just, A., Khan, F. M., Berczik, P., Ernst, A., & Spurzem, R. 2010, MNRAS, 411, 653 [Google Scholar]
 Kadota, K., Kim, J. H., Ko, P., & Yang, X.Y. 2024, Phys. Rev. D, 109, 015022 [NASA ADS] [CrossRef] [Google Scholar]
 Kahlhoefer, F., SchmidtHoberg, K., Frandsen, M. T., & Sarkar, S. 2014, MNRAS, 437, 2865 [CrossRef] [Google Scholar]
 Karydas, T. K., Kavanagh, B. J., & Bertone, G. 2024, ArXiv eprints [arXiv:2402.13053] [Google Scholar]
 Kavanagh, B. J., Nichols, D. A., Bertone, G., & Gaggero, D. 2020, Phys. Rev. D, 102, 083006 [NASA ADS] [CrossRef] [Google Scholar]
 Kavanagh, B. J., Karydas, T. K., Bertone, G., Cintio, P. D., & Pasquato, M. 2024, ArXiv eprints [arXiv:2402.13762] [Google Scholar]
 Kim, H., & Kim, W. 2007, ApJ, 665, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Koda, J., & Shapiro, P. R. 2011, MNRAS, 415, 1125 [NASA ADS] [CrossRef] [Google Scholar]
 Kong, D., Yang, D., & Yu, H.B. 2024, ApJ, 965, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Kremer, G. M. 2022, Universe, 8, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Lancaster, L., Giovanetti, C., Mocz, P., et al. 2020, JCAP, 2020, 001 [Google Scholar]
 Macedo, C. F. B., Pani, P., Cardoso, V., & Crispino, L. C. B. 2013, ApJ, 774, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Milosavljević, M., & Merritt, D. 2003, ApJ, 596, 860 [Google Scholar]
 Monaghan, J. J., & Lattanzio, J. C. 1985, A&A, 149, 135 [NASA ADS] [Google Scholar]
 Montalvo, D., SmithOrlik, A., & Rastgoo, S., et al. 2024, ArXiv eprints [arXiv:2401.06084] [Google Scholar]
 Mukherjee, D., Holgado, A. M., Ogiya, G., & Trac, H. 2024, MNRAS, 533, 2335 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, E. C. 1999, ApJ, 513, 252 [Google Scholar]
 Outmezguine, N. J., Boddy, K. K., GadNasr, S., Kaplinghat, M., & Sagunski, L. 2023, MNRAS, 523, 4786 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J., Spergel, D. N., & Steinhardt, P. J. 2015, ApJ, 804, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Ragagnin, A., Meneghetti, M., Calura, F., et al. 2024, A&A, 687, A270 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Rephaeli, Y., & Salpeter, E. E. 1980, ApJ, 240, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, A., Harvey, D., Massey, R., et al. 2019, MNRAS, 488, 3646 [NASA ADS] [CrossRef] [Google Scholar]
 Sabarish, V. M., Brüggen, M., SchmidtHoberg, K., Fischer, M. S., & Kahlhoefer, F. 2024, MNRAS, 529, 2032 [NASA ADS] [CrossRef] [Google Scholar]
 Sadeghian, L., Ferrer, F., & Will, C. M. 2013, Phys. Rev. D, 88, 063522 [NASA ADS] [CrossRef] [Google Scholar]
 Sagunski, L., GadNasr, S., Colquhoun, B., Robertson, A., & Tulin, S. 2021, JCAP, 2021, 024 [CrossRef] [Google Scholar]
 Shapiro, S. L. 2018, Phys. Rev. D, 98, 023021 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, S. L., & Paschalidis, V. 2014, Phys. Rev. D, 89, 023506 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
 Szölgyén, Á., MacLeod, M., & Loeb, A. 2022, MNRAS, 513, 5465 [CrossRef] [Google Scholar]
 Tang, M., Xu, Z., & Wang, J. 2021, Chin. Phys. C, 45, 015110 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, J. H., Hulse, R. A., Fowler, L. A., Gullahorn, G. E., & Rankin, J. M. 1976, ApJ, 206, L53 [NASA ADS] [CrossRef] [Google Scholar]
 Traykova, D., Vicente, R., Clough, K., et al. 2023, Phys. Rev. D, 108, L121502 [NASA ADS] [CrossRef] [Google Scholar]
 Tulin, S., & Yu, H.B. 2018, Phys. Rep., 730, 1 [Google Scholar]
 Ullio, P., Zhao, H., & Kamionkowski, M. 2001, Phys. Rev. D, 64, 043504 [NASA ADS] [CrossRef] [Google Scholar]
 Vogelsberger, M., Zavala, J., Simpson, C., & Jenkins, A. 2014, MNRAS, 444, 3684 [CrossRef] [Google Scholar]
 Wang, Y., & Easther, R. 2022, Phys. Rev. D, 105, 063523 [NASA ADS] [CrossRef] [Google Scholar]
 Weisberg, J. M., Taylor, J. H., & Fowler, L. A. 1981, Sci. Am., 245, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, D., & Yu, H.B. 2022, JCAP, 2022, 077 [CrossRef] [Google Scholar]
 Yang, S., Du, X., Zeng, Z. C., et al. 2023, ApJ, 946, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, S., Jiang, F., & Benson, A., et al. 2024, MNRAS, 533, 4007 [NASA ADS] [CrossRef] [Google Scholar]
 Yue, X.J., & Han, W.B. 2018, Phys. Rev. D, 97, 064003 [NASA ADS] [CrossRef] [Google Scholar]
 Zeng, Z. C., Peter, A. H. G., Du, X., et al. 2022, MNRAS, 513, 4845 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, X., Yu, H.B., Yang, D., & An, H. 2024, ApJ, 968, L13 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Unresolved mean free path
In this appendix, we show the results of two simulations with a larger crosssection than the ones we have shown in the main text. These additional simulations are for fSIDM and employ a crosssection of σ_{V}/m = 3 × 10^{−5} cm^{2} g^{−1}. They differ in resolution, we use N = 1 × 10^{6} and N = 5 × 10^{6} particles. In Fig. A.1, we can see, in contrast to the previous simulations, that the results for the two resolutions differ significantly from each other.
Fig. A.1. Position and velocity of the BH in our test setup for a simulation with a fSIDM crosssection of σ_{V}/m = 3 × 10^{−5} cm^{2} g^{−1}. The high resolution corresponds to N = 5 × 10^{6} particles and the low resolution to N = 1 × 10^{6} particles. The larger box size of l_{box} = 829.5 μpc was used. 
With increasing crosssection, the mean free path is shrinking and may fall below the kernel size used for SIDM (see Fischer et al. 2024a). The ratio between the kernel size, h, and the mean free path l, is fairly large for these simulations, h/l ≈ 172 for the lower resolution run and h/l ≈ 100 for the higher resolution. Thus, it is not surprising that the simulations deviate from each other. The value for h/l at which the simulations start to become inaccurate may in general depend on the problem at hand. The values of our simulations here may not be indicative of other setups.
All Tables
All Figures
Fig. 1. Density and onedimensional velocity dispersion for a DM density spike as a function of radius according to Eq. (2). We assume r_{sp} = 1 μpc, ρ_{sp} = 5.45 × 10^{−3} M_{⊙} μpc^{−3}, α_{sp} = 7/3, and M_{BH} = 10^{3} M_{⊙}. In addition, we indicate the Schwarzschild radius R_{s} = 2GM/c^{2}. For computing the velocity dispersion, we followed Kremer (2022). Here, we assumed an isotropic velocity distribution. In contrast to Kremer (2022), we solved their Eq. (27) including higherorder terms instead of Eq. (28). Moreover, we show the results for m/(k_{B}T)≈6.53 × 10^{3} s^{2} km^{−2}, with m being the rest mass of a physical DM particle, k_{B} the Boltzmann constant and T the absolute temperature of the system. 

In the text 
Fig. 2. Dynamical friction force as a function of the Mach number, i.e. the velocity in units of the sound speed. The dashed lines indicated the dynamical friction for a collisionless medium according to the Chandrasekhar formula (Chandrasekhar 1943) and the solid lines are for a gaseous medium as given by Ostriker (1999). We note that for the plot, we assume r_{max} = ℳ c_{s} t and ln(r_{max}/r_{min}) = ln(ℳ)+C. The line colours correspond to different values of C ∈ {4.0, 6.0, 8.0, 10.0, 12.0}. In addition, the red dotted lines indicate the velocity of a secondary BH moving on a circular orbit for different adiabatic indices. 

In the text 
Fig. 3. Velocity and position of the BH in our test setup as a function of time. All simulations are for collisionless DM. The solid lines are from the setup with a box length of l = 829.5 μpc. In contrast, the dashed lines correspond to a setup with a smaller box size of l = 663.6 μpc. The highresolution results correspond to N_{DM} = 5 × 10^{6} DM particles, and the lowresolution ones to N_{DM} = 10^{6}. We note that in the smaller box setup, the low resolution implies a higher mass resolution than in the large box setup. 

In the text 
Fig. 4. As in Fig. 3, position and velocity of the BH as a function of time. The highresolution setup with the large box size of l = 829.5 μpc is displayed. We fitted Chandrasekhar’s formula (Eqs. 6 and 7) to the simulation results to determine the Coulomb logarithm, ln(Λ), with Λ = r_{max}/r_{min}. Besides the simulation result (dotted lines), we also show the fitted curves (solid lines). 

In the text 
Fig. 5. Similar to Fig. 3, position and velocity of the BH particle as a function of time. The results are shown for several interaction strengths of the DM, indicated by the opacity of the lines. All simulations are for the large box with the high resolution of N_{DM} = 5 × 10^{6} particles. For the simulation of gaseous matter, we used the same number of SPH particles. 

In the text 
Fig. 6. Surface density (top row) and the velocity dispersion (bottom row) for our test simulations at a time of t = 0.49 yr. The panels display the different types of matter we investigate, with increasing collisionality from the left to the right. The density and one dimensional velocity dispersion are indicated following a logarithmic colour scale spanning a factor of two. For all panels of the same quantity, the same colour scale is used. In addition, we indicate the position of the BH particle with a red dot. We provide the time evolution of the surface density as a video available online. 

In the text 
Fig. 7. Dynamical friction force as a function of the Mach number for our simulations. This plot is analogue to Fig. 2, which gives the analytic description of the dynamical friction force. Different crosssections are shown, indicated by the Knudsen number. Here, Kn = ∞ corresponds to a collisionless medium and Kn = 0.0 to a gaseous medium. We note that in the simulation the BH initially starts at a high velocity, i.e. at the right side of the plot and slows down throughout the simulation, i.e. moves towards the left side of the plot. 

In the text 
Fig. 8. Density of the wake, i.e. the mean density about the point with the highest density, as a function of the velocity of the perturbing BH. The simulation results of the different crosssections are indicated by their Knudsen number. 

In the text 
Fig. 9. Onedimensional velocity dispersion of the wake, i.e. the velocity dispersion in the vicinity of the point with the highest density, as a function of the velocity of the perturber. We note that the velocity dispersion is computed after the bulk motion has been subtracted. The simulation results for different crosssections are shown, indicated by the Knudsen number. 

In the text 
Fig. 10. Mean velocity of the wake in the direction of the motion of the BH, i.e. the mean velocity in the vicinity of the point with the highest density. The results for different crosssections are shown as a function of the velocity of the perturber. They are indicated by their corresponding Knudsen number. 

In the text 
Fig. 11. Dark matter mass gravitational bound to the perturbing BH as a function of the velocity of the perturber. We indicate the different crosssections by their Knudsen number. 

In the text 
Fig. 12. Ratio between the velocity change, Δv = v(0)−v(t), due to dynamical friction of an fSIDM crosssection and the transfer (blue) / viscosity (orange) matched isotropic crosssections. The upper panel is for an fSIDM viscosity crosssection of σ_{V}/m = 3 × 10^{−8} cm^{2} g^{−1} and the lower panel is for σ_{V}/m = 3 × 10^{−7} cm^{2} g^{−1}. 

In the text 
Fig. 13. Ratio of the velocity change, Δv = v(0)−v(t), due to dynamical friction acting on the BH for a velocitydependent and a velocityindependent (constant) crosssection. All crosssections were matched using the effective crosssection (Eq. 11), i.e. all crosssections correspond to σ_{eff}/m = 3 × 10^{−8} cm^{2} g^{−1}. The upper panel displays runs with frequent selfinteractions, and the lower panel is for isotropic scattering. 

In the text 
Fig. A.1. Position and velocity of the BH in our test setup for a simulation with a fSIDM crosssection of σ_{V}/m = 3 × 10^{−5} cm^{2} g^{−1}. The high resolution corresponds to N = 5 × 10^{6} particles and the low resolution to N = 1 × 10^{6} particles. The larger box size of l_{box} = 829.5 μpc was used. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.