# Cosmological and idealised simulations of dark matter haloes with velocity-
dependent, rare and frequent self-interactions

Moritz S. Fischer,1,2,3 Lenard Kasselmann,3 Marcus Brüggen,3 Klaus Dolag,1,4
Felix Kahlhoefer,5 Antonio Ragagnin,6,7,8 Andrew Robertson,9 Kai Schmidt-
Hoberg10
1Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität
München, Scheinerstr. 1, D-81679 München, Germany
2Excellence Cluster ORIGINS, Boltzmannstrasse 2, D-85748 Garching, Germany
3Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029
Hamburg, Germany
4Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748
Garching, Germany
5Institute for Theoretical Particle Physics (TTP), Karlsruhe Institute of
Technology (KIT), D-76128 Karlsruhe, Germany
6Dipartimento di Fisica e Astronomia "Augusto Righi", Alma Mater Studiorum
Università di Bologna, via Gobetti 93/2, I-40129 Bologna, Italy
7INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143
Trieste, Italy
8IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2,
I-34014 Trieste, Italy
9Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA
10Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
E-mail<EMAIL_ADDRESS>(LMU)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Dark matter self-interactions may have the capability to solve or at least
mitigate small-scale problems of the cosmological standard model,
$\Lambda$CDM. There are a variety of self-interacting dark matter (SIDM)
models that lead to distinguishable astrophysical predictions and hence
varying success in explaining observations. Studies of dark matter (DM)
density cores on various mass scales suggest a velocity-dependent scattering
cross-section. In this work we investigate how a velocity dependence alters
the evolution of the DM distribution for _frequent_ DM scatterings and compare
to the velocity-independent case. We demonstrate that these cases are
qualitatively different using a test problem. Moreover, we study the evolution
of the density profile of idealised DM haloes and find that a velocity
dependence can lead to larger core sizes and different time scales of core
formation and core collapse. In cosmological simulations, we investigate the
effect of velocity-dependent self-interaction on haloes and satellites in the
mass range of $\approx 10^{11}$–$10^{14}\,\mathrm{M_{\odot}}$. We study the
abundance of satellites, density and shape profiles and try to infer
qualitative differences between velocity-dependent and velocity-independent
scatterings as well as between frequent and rare self-interactions. We find
that a strongly velocity-dependent cross-section can significantly amplify the
diversity of rotation curves, independent of the angular dependence of the
differential cross-section. We further find that the abundance of satellites
in general depends on both the velocity dependence and the scattering angle,
although the latter is less important for strongly velocity-dependent cross-
sections.

###### keywords:

astroparticle physics – methods: numerical – galaxies: haloes – dark matter

††pubyear: 2023††pagerange: Cosmological and idealised simulations of dark
matter haloes with velocity-dependent, rare and frequent self-interactions–F

##  1 Introduction

Historically, dark matter (DM) self-interactions have been motivated to solve
problems on small, i.e. galactic scales. It was found that cosmological DM-
only simulations can explain the large-scale structure of the universe quite
well. But on smaller scales, deviations between the observations and
simulations were encountered (e.g. Moore et al., 1998). Spergel & Steinhardt
(2000) proposed self-interacting dark matter (SIDM) as a solution to two
problems on small scales. Namely, SIDM can reduce the abundance of satellites
and the central density of haloes. As the self-interactions lead to heat flow
into the central region of a Navarro-Frenk-White (NFW, Navarro et al., 1996)
halo, they reduce the central density and can form density cores. The first
$N$-body simulation using a Monte Carlo scheme of this core formation has been
performed by Burkert (2000). Since then SIDM has been found to be capable of
solving or at least mitigating further small-scale problems of Cold Dark
Matter (CDM) (for a review see Tulin & Yu, 2018; Adhikari et al., 2022). This
does not only include the core-cusp problem (e.g. Davé et al., 2001), but also
diverse rotation curves (e.g. Creasey et al., 2017; Kamada et al., 2017;
Robertson et al., 2018; Correa et al., 2022) and the too-big-to-fail problem
(e.g. Zavala et al., 2013; Elbert et al., 2015; Kaplinghat et al., 2019). For
a review of small-scale problems in Lambda cold dark matter ($\Lambda$CDM), we
refer the reader to Bullock & Boylan-Kolchin (2017).

Meanwhile, it has also emerged that there are other avenues to solve these
small-scale problems. On the one hand, it was found that including the
baryonic physics, in particular, feedback mechanisms from supernovae (e.g.
Read & Gilmore, 2005; Governato et al., 2012; Pontzen & Governato, 2012) and
black holes can form density cores (e.g. Martizzi et al., 2013; Silk, 2017;
Peirani et al., 2017). On the other hand, researchers have become more
cautious about inferring density profiles from rotation curves (e.g. Pineda et
al., 2016; Read et al., 2016b; Genina et al., 2018; Oman et al., 2019; Roper
et al., 2023; Downing & Oman, 2023). Beyond SIDM, other DM models have been
investigated, including warm DM (Dodelson & Widrow, 1994) and fuzzy DM (Hu et
al., 2000).

Although SIDM has initially been mainly motivated by small-scale issues, it
provides DM candidates worth investigating, independent of the state of the
small-scale crisis. The nature of DM is still unknown and could have
properties which we can only infer indirectly via astronomical observations.
This is true for models of SIDM, and studying them is essentially constraining
particle physics properties of DM. Particle candidates that fall into the
class of SIDM can have various characteristics. The scattering may be elastic
or inelastic, it may involve multiple states and can feature different angular
dependencies. Another aspect is how the cross-section depends on the relative
velocity of the scattering particles.

Velocity-dependent self-interactions have been recently studied by various
authors (e.g. Colin et al., 2002; Nadler et al., 2020; Yang & Yu, 2022;
Outmezguine et al., 2023; Yang et al., 2023b). Such studies were performed not
only with DM-only (DMO) simulations but also within hydrodynamical
cosmological simulations (e.g. Vogelsberger et al., 2014; Robertson et al.,
2019, 2020; Rose et al., 2022; Mastromarino et al., 2023; Rahimi et al.,
2023). They are well-motivated for different angular dependencies, including
forward-enhanced cross-sections from light mediator models (e.g. Buckley &
Fox, 2010; Loeb & Weiner, 2011; Bringmann et al., 2017). But also models of
resonant scattering (e.g. Chu et al., 2019; Tsai et al., 2022) can explain a
velocity dependence while featuring an isotropic cross-section.

From an astronomical perspective, velocity-dependent self-interactions are
well motivated (e.g. Kaplinghat et al., 2016; Correa, 2021; Gilman et al.,
2021; Sagunski et al., 2021; Silverman et al., 2022; Lovell & Zavala, 2023).
They would allow fulfilling stringent constraints from galaxy clusters while
having a fairly large effect on low-mass haloes. When the self-interaction
cross-section decreases with velocity, it has a weaker effect in galaxy
clusters because their typical relative DM velocities are larger than in
galaxies. Furthermore, they can lead to a qualitative different evolution of
systems that involve multiple velocity scales. For instance, this is true for
the evolution of the satellite distribution (e.g. Zeng et al., 2022) and could
lead to an increase in the diversity of density profiles and rotation curves
(e.g. Nadler et al., 2023; Yang et al., 2023c).

The aim of this study is to explore qualitative differences arising from the
velocity dependence of the self-interactions and to understand their
implications on constraining the angular dependence of the cross-section. In
this paper, we consider two different angular dependencies: Firstly, isotropic
scattering, to which we refer as rare self-interactions (rSIDM). Secondly, a
cross-section with typical scattering angles that are very small. In
consequence, frequent interactions are needed to significantly alter the DM
distribution. Hence, we call it frequent self-interactions (fSIDM). We explore
rSIDM and fSIDM models with several velocity dependencies to study qualitative
differences arising from the velocity and angular dependence. The scattering
of all SIDM models we consider is elastic. For our study, we employ idealised
N-body simulations of a test problem and DM haloes as well as cosmological
simulations. Unlike velocity-independent models (Fischer et al., 2022), fSIDM
with a velocity-dependent interaction has not been studied in a cosmological
context. Finally, all our simulations are DM-only, i.e. we ignore the effects
of baryons. In a companion paper (Sabarish et al., prep) velocity-dependent
fSIDM is studied in the context of merging galaxy clusters.

This paper is structured as follows: In Section 2 we describe the numerical
setup of our simulations including a novel time-stepping criterion. A
presentation of the simulations and our results follows for the idealised
setups in Section 3 and the cosmological simulations in Section 4.
Shortcomings and directions for further research are discussed in Section 5.
Finally, in Section 6 we conclude. Additional information can be found in the
appendices.

##  2 Numerical setup

In this section, we describe our numerical setup. First, we begin by
describing the simulation code and the SIDM implementation. We continue with
the parametrisation for the velocity-dependent cross-section. Next, we
introduce a novel time-step criterion for the velocity-dependent self-
interaction. Lastly, the simulations with their initial conditions and the
identification of the substructure are described. In addition, a description
of our improved parallelisation scheme for SIDM can be found in Appendix A.

###  2.1 SIDM implementation and simulations

For our simulations, we use the cosmological hydrodynamical $N$-body code
gadget-3. The predecessor gadget-2 has been described in Springel (2005).
Various additional modules have been developed for the OpenGadget3 version
that we are using. The implementation of DM self-interactions has been
described by Fischer et al. (2021a, b, 2022). To date, this is the only
implementation for simulating frequently self-interacting DM. The code is also
able to simulate isotropic cross-sections. Here, the employed scheme is very
similar to the one introduced by Rocha et al. (2013), except that we use an
adaptive kernel size set by the 64 next neighbours and a different time step
criterion.

We have run several simulations of CDM, rSIDM and fSIDM for idealised setups
with individual haloes as well as cosmological simulations. For all
simulations we used the cosmological $N$-body code gadget-3. The details of
the simulations can be found in the corresponding Sections 3 and 4. In
addition, we ran simulations to test the code, they can be found in the
Appendices C and D.

###  2.2 Velocity-dependent cross-section

There are numerous studies in the literature considering a cross-section,
$\sigma$, that depends on the scattering velocity, $v$. A typical choice –
that we employ as well – is a cross-section that scales as $\sigma\propto
v^{-4}$ in the limit of high $v$. This dependence may be motivated by particle
physics (e.g. Ibe & Yu, 2010; Tulin et al., 2013) and has been employed in
numerous studies (e.g. Kaplinghat et al., 2016; Robertson et al., 2017b).

Following Robertson et al. (2017b); Kahlhoefer et al. (2017), we consider the
momentum transfer cross-section

$\sigma_{\mathrm{T}}=2\pi\int_{-1}^{1}\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega_{\mathrm{cms}}}\left(1-|\cos\theta_{\mathrm{cms}}|\right)\mathrm{d}\cos\theta_{\mathrm{cms}}\,.$
(1)

We parameterise the velocity dependence of the momentum transfer cross section
as

$\frac{\sigma_{\mathrm{T}}}{m}=\frac{\sigma_{0}}{m}\left(1+\left(\frac{v}{w}\right)^{\beta}\right)^{\alpha/\beta}\,.$
(2)

Here, $\sigma_{0}$ corresponds to the cross-section in the velocity-
independent regime, $w$ denotes the velocity cutoff, $\alpha$ sets the decline
at high velocities and $\beta$ describes the transition from the constant
cross-section at low velocities to the decreasing cross-section at high
velocities. In this study, we always set $\alpha=-4$ and $\beta=2$. This
choice is motivated by the fact that in the limit of the Born-approximation,
the velocity dependence of the total and the transfer cross-section are very
similar (Ibe & Yu, 2010). More details on the transfer cross-section and the
possible connections to the underlying particle physics can be found in the
companion paper (Sabarish et al., prep).

In most physically motivated cases, a velocity dependence also implies an
angular dependence of the differential scattering cross-section. $N$-body
simulations had been limited in simulating frequent scatterings about small
angles until the work by Fischer et al. (2021a). Here, we go beyond the common
large-angle scattering and investigate small-angle as well as isotropic
scattering combined with a velocity dependence.

In order to probe different velocity regimes, we use several combinations of
$\sigma_{0}$ and $w$. These are described together with the details of the
simulations in Section 3 and Section 4. Each parameter set is simulated with
fSIDM and rSIDM, the latter corresponding to isotropic scattering. Note that
we use the momentum transfer cross-section (Eq. 1) to match fSIDM and rSIDM.
In the case of isotropic scattering, the total cross-section is twice as large
as the momentum transfer cross-section.

###  2.3 Time step criterion

For velocity-dependent self-interactions, a separate time-step criterion can
become more important than for velocity-independent scatterings because cross-
sections can become large at low velocities. Depending on the cross-section
this can give more stringent limitations on the time step than imposed by the
gravity scheme. We found that the time-step criterion introduced by Fischer et
al. (2021b) for velocity-independent self-interactions is not always well-
suited for a velocity-dependent cross-section (this has been previously
described by Kasselmann, 2021). Here we introduce a new time step criterion
for velocity-dependent scattering that has a velocity-dependence as described
by Eq. 2. In more general terms, our time step criterion requires that there
is a finite velocity for which the fractional velocity change due to the drag
force becomes maximal and finite. We remind the reader that the effective drag
force for fSIDM was introduced by Kahlhoefer et al. (2014) and employed to
develop a numerical scheme by Fischer et al. (2021a). It is given as

$F_{\mathrm{drag}}=\frac{1}{2}\,\frac{\sigma_{\mathrm{T}}(v)}{m}\,v^{2}\,m_{\mathrm{n}}^{2}\,\Lambda\,.$
(3)

The relative particle velocity is denoted by $v$, $m_{\mathrm{n}}$ is the
numerical particle mass and $\Lambda$ is the kernel overlap, a geometrical
factor (for details see Fischer et al., 2021a).

Assuming the parametrisation according to Eq. 2 the fractional velocity change
($\Delta v/v$) due to the drag force becomes maximal for the velocity

$v_{e}=\frac{w}{\left(-1-\alpha\right)^{1/\beta}}\,.$ (4)

Note that this is only applicable if $\alpha<-1$ and $\beta>0$. For our choice
of $\alpha=-4$ and $\beta=2$ this implies $v_{e}=w/\sqrt{3}$.

Using the maximum allowed fractional velocity change $\tau$, we can express
the time-step criterion for particle $i$ as

$\Delta
t_{i}<\tau\,\frac{2}{v_{e}}\frac{1}{m_{\mathrm{n}}\,\Lambda_{ii}}\left(\frac{\sigma_{\mathrm{T}}(v_{e})}{m}\right)^{-1}\,.$
(5)

Here, $\Lambda_{ii}$ gives the maximal possible kernel overlap by calculating
it with the particle itself.

It is worth pointing out that this time step depends on the chosen number of
neighbours, $N_{\mathrm{ngb}}$. With a larger number of neighbours
$\Lambda_{ii}$ becomes smaller and thus the time step is larger and vice
versa. Finally, we note that this time step criterion also applies to rSIDM
when using the total cross-section, $\sigma$, instead of
$\sigma_{\mathrm{T}}$. For rSIDM, the scattering probability reaches a maximum
at $v_{e}$ (see Eq. 4) too.

In Appendix B, we provide further discussion on issues related to the
formulation of a time step criterion.

##  3 Idealised simulations

In this section, we present and analyse our idealised simulations and show the
results we obtain. First, we start with a simple test problem in Sec. 3.1.
Secondly, the evolution of the core size for isolated haloes is shown (Sec.
3.2) for both initial Hernquist and NFW profiles.

###  3.1 Thermalisation problem

name | type | $\sigma_{0}/m$ | $w$
---|---|---|---
|  | $[\mathrm{cm}^{2}\,\mathrm{g}^{-1}]$ | $[\mathrm{km}\,\mathrm{s}^{-1}]$
f10 | frequent | $10$ | –
r10 | rare | $10$ | –
f4.5e6w0.1 | frequent | $4.5\times 10^{6}$ | $0.1$
r4.5e6w0.1 | rare | $4.5\times 10^{6}$ | $0.1$
Table 1:  The table shows the different cross-sections that we used for the
thermalisation problem. The first column gives the name that we use in the
paper to abbreviate the cross-section. It follows the type of self-
interaction. Here, “rare” corresponds to isotropic scattering. The third
column gives $\sigma_{0}/m$ and the last one $w$ (see also Eq. 2).

To learn about the differences between a constant and a velocity-dependent
cross-section, we first consider the thermalisation problem previously studied
by Fischer et al. (2021a). This has the advantage that we study the pure
effect of DM self-interactions without the influence of gravity.

The numerical setup consists of a periodic box with a constant density of
$10^{7}\,\mathrm{M_{\odot}}\,\mathrm{kpc}^{-3}$ sampled by $10^{4}$ particles.
The cubic box has a side length of 10 kpc and its particles have a velocity of
$2\,\mathrm{km}\,\mathrm{s}^{-1}$ which points into a random direction. In
Tab. 1, we describe the employed cross-sections. The scattering broadens the
velocity distribution such that it evolves towards a Maxwell–Boltzmann
distribution. We can characterise the width of the distribution of the
absolute velocities by computing its variance.

Figure 1: The variance for the distribution of absolute velocities of the
thermalisation problem introduced by Fischer et al. (2021a) is shown. We
display the results for different SIDM models as a function of time. In black
we indicated the variance of the final Maxwell–Boltzmann distribution.

In Fig. 1 we show the results as a function of time. For frequent self-
interactions, this has been previously studied by Kasselmann (2021). In line
with his results, we find that the evolution of the thermalisation rate
evolves qualitatively differently for velocity-dependent self-interactions
compared to a constant cross-section. The thermalisation process evolves
faster at early times and slower at late times for the velocity-dependent
self-interactions. For the isotropic cross-section, we find the same.
Initially, the system evolves faster for the velocity-dependent cross-
sections, because the cross-section evaluated at the typically relative
velocity of the particles is larger compared to the velocity-independent
cross-section. The lower thermalisation rate at late times, i.e. when the
velocity distribution is already close to the Maxwell–Boltzmann distribution,
stems mainly from a deviation at the high-velocity tail. The decrease of the
cross-section with velocity makes velocity-dependent self-interactions less
efficient in scattering particles to high velocities. In consequence, the
thermalisation rate in a late stage is reduced.

###  3.2 Isolated haloes

Here, we study the evolution of isolated haloes subject to velocity-dependent
self-interactions. Firstly, we investigate the density profile of an isolated
halo with a density following a Hernquist profile (Hernquist, 1990) and
secondly, we do the same for a halo with an NFW profile (Navarro et al.,
1996). For the two haloes, we also compare rare and frequent self-
interactions.

####  3.2.1 Hernquist Halo

name | type | $\sigma_{0}/m$ | $w$
---|---|---|---
|  | $[\mathrm{cm}^{2}\,\mathrm{g}^{-1}]$ | $[\mathrm{km}\,\mathrm{s}^{-1}]$
c0 | collisionless | $0.0$ | –
f0.8 | frequent | $0.8$ | –
r0.8 | rare | $0.8$ | –
f1e5w100 | frequent | $10^{5}$ | $100$
r1e5w100 | rare | $10^{5}$ | $100$
Table 2:  The cross-sections that we employed for simulating a Hernquist halo
are shown. The columns are the same as in Tab. 1.  Figure 2: The size of the
density core for a Hernquist halo as a function of time is shown when evolved
with different DM models. We indicate the cross-section in the legend. The
first number refers to $\sigma_{0}/m$ in units of
$\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ and the second one to $w$ in units of
$\mathrm{km}\,\mathrm{s}^{-1}$ (see Tab. 2). The first two SIDM simulations
are for a velocity-independent cross-section and the number gives
$\sigma_{\mathrm{T}}/m$.

We simulate the same Hernquist halo as first described by Robertson et al.
(2017b). It has a mass of $M=2.46\times 10^{14}\,\mathrm{M_{\odot}}$ and a
scale radius of $r_{s}=279\,\mathrm{kpc}$. We generate the initial conditions
by sampling the halo up to $r=400\,r_{s}$ using $N=10^{7}$ particles. For the
gravitational softening length we employ $\epsilon=0.56\,\mathrm{kpc}$. The
simulations include velocity-independent and velocity-dependent cross-sections
both for fSIDM and rSIDM. In detail, the cross-sections are shown in Tab. 2.
It is worth noting that for the velocity-dependent simulations, the SIDM time-
step constraint was tighter than the one from gravity, at least for a fraction
of the particles. This led to a significant increase in computational costs.
We determine the core size, $r_{\mathrm{core}}$, as previously done by
Robertson et al. (2017a) and Fischer et al. (2021a) by fitting a cored
Hernquist profile. It is given as

$\rho(r)=\frac{M}{2\pi}\frac{r_{s}}{(r^{4}+r^{4}_{\mathrm{core}})^{1/4}}\frac{1}{(r+r_{s})^{3}}\,.$
(6)

As in the original Hernquist profile, $M$ denotes the halo mass and
$r_{\mathrm{s}}$ the scale radius. The evolution of the core size is shown in
Fig. 2 for the different DM models.111We found the exact core size to be
sensitive to details of the optimisation procedure, which might be caused by a
noisy likelihood. This might be the main source of different core sizes for
the same halo in the literature (Robertson et al., 2017b; Fischer et al.,
2021a; Correa et al., 2022).

In the early stages, the density core grows due to self-interactions whose
effect can be described as heat transfer (e.g. Lynden-Bell & Eggleton, 1980;
Balberg et al., 2002) that follows the gradient of the velocity dispersion. As
a result, the central region of the halo heats up and its density is
decreasing. For the collisionless DM, we find a small core caused by
gravitational two-body interactions, a process known as numerical core
formation (e.g. Dehnen, 2001). At later stages, the core size is decreasing
and the halo enters the collapse phase. In this phase, heat is only
transported outward, as the central region cools it also contracts.
Gravitational bound systems are characterised by a negative heat capacity.
This is for example well known from star clusters but also applies to the
haloes we study here. In consequence, the velocity dispersion at the central
region of the halo is increasing and leads to a runaway process called the
gravothermal catastrophe.

In previous studies it was found that the maximum core size that is reached
during the haloes evolution is roughly independent of the strength of the
cross-section (e.g. Kochanek & White, 2000), but also its angular dependence
(e.g. Robertson et al., 2017a; Fischer et al., 2021a). In contrast, we find
that the velocity-dependent cross-sections give a larger maximum core size.
However, we have to note that this only occurs for sufficiently small values
of $w$. For the initial Hernquist halo heat is flowing inwards for radii
smaller than the radius of the maximal velocity dispersion,
$r({\nu^{2}_{\mathrm{max}}})$, this should set the core formation time. In
contrast, for radii larger than $r({\nu^{2}_{\mathrm{max}}})$, heat is flowing
outward, determining the core collapse time. The maximum core size should be a
result of the ratio of the total heat in and outflow. In consequence, a DM
candidate that is more efficient in transporting heat inward than outward
compared to other DM models would produce a larger maximum core size. We
discuss this further in Sec. 3.2.3, after we have shown the results for the
isolated NFW halo.

However, to gain further insights into the halo following initially a
Hernquist profile, we first plot various quantities at the time of maximum
core expansion in Fig. 3. The upper panel shows the density and velocity
dispersion profile, and the bottom panel displays quantities related to heat
conductivity.

In the following, we describe how we compute the quantities of the bottom
panel. Assuming identical particles the viscosity cross-section is given by

$\sigma_{v}=4\pi\int_{0}^{1}\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}\sin^{2}\theta\,\mathrm{d}\cos\theta\,.$
(7)

Based on this we can express the effective cross-section of Yang & Yu (2022)
as

$\sigma_{\mathrm{eff}}=\frac{3}{2}\frac{\langle
v^{5}\sigma_{v}(v)\rangle}{\langle v^{5}\rangle}\,.$ (8)

They introduced the effective cross-section with the aim of matching
differential cross-sections with various angular and velocity dependencies. It
thus allows transferring constraints on the strength of self-interactions to
various SIDM models. Here, the average is computed assuming the velocities are
well described by a Maxwell–Boltzmann distribution. Next, we give the heat
conductivity using $\sigma_{\mathrm{eff}}$. Strictly speaking, we do not
specify the heat conductivity $\kappa$, but use
$\kappa^{\prime}=m/\mathrm{k_{B}}\kappa$, with $m$ the DM particle mass and
$\mathrm{k_{B}}$ the Boltzmann constant. This is commonly used in the
gravothermal fluid model (e.g. Koda & Shapiro, 2011). Note, Kummer et al.
(2019) took the angular dependence into account by expressing the heat
conductivity in terms of the viscosity cross-section. Here we go further and
use the effective cross-section for $\kappa^{\prime}$. For the short-mean-
free-path (smfp) regime it is given as

$\kappa^{\prime}_{\mathrm{smfp}}=\frac{9\,b\,\nu}{4}\left(\frac{\sigma_{\mathrm{eff}}}{m}\right)^{-1}\quad\textnormal{with}\quad
b=\frac{25\sqrt{\pi}}{32}\,.$ (9)

The one-dimensional velocity dispersion is expressed by $\nu^{2}$. In the
long-mean-free-path (lmfp) regime the heat conductivity can be expressed as

$\kappa^{\prime}_{\mathrm{lmfp}}=\hat{a}\,C\,\frac{\nu^{3}\rho}{4\pi\mathrm{G}}\left(\frac{\sigma_{\mathrm{eff}}}{m}\right)\quad\textnormal{with}\quad\hat{a}=\sqrt{\frac{16}{\pi}},\,C\approx
0.75\,.$ (10)

Here, $\rho$ denotes the density and $\mathrm{G}$ is the gravitational
constant.

The Knudsen number, $\mathrm{Kn}$, is usually used to distinguish between the
lmfp and smfp regime and is defined as

$\mathrm{Kn}=\frac{3}{2}\,\sqrt{\frac{4\pi\mathrm{G}}{\rho\nu^{2}}}\left(\frac{\sigma_{\mathrm{eff}}}{m}\right)^{-1}\;.$
(11)

Numerically $\mathrm{Kn}>1$ corresponds to the lmfp regime and $\mathrm{Kn}<1$
to the smfp regime.

Figure 3: Various properties of the halo following initially a Hernquist
profile are shown at the evolution stage when its density core is the largest.
In the upper panel, we show the density (black) and the velocity dispersion
(blue) as a function of radius. Moreover, the scale radius, $r_{s}$, and the
radius at which the velocity dispersion of the initial profile reaches its
maximum, $r(\nu^{2}_{\mathrm{max}})$, are indicated. The lower panel gives
$\kappa^{\prime}$ for the smfp (grey) and lmfp (black) regime (see Eq. 9 and
Eq. 10) as well as the Knudsen number (see Eq. 11). These quantities are
computed based on the effective cross-section, $\sigma_{\mathrm{eff}}/m$. In
addition, the maximum core sizes are shown for the runs with the frequent
self-interactions, i.e. for the velocity-independent and velocity-dependent
cross-sections. To compute the quantities that are shown as a function of
radius we used the simulation with frequent self-interactions and without
velocity-dependence.

From the upper panel of Fig. 3, we can see that the velocity dispersion at the
time of maximum core expansion is roughly constant for radii smaller than
$r(\nu^{2}_{\mathrm{max}})$. However, the density core itself is much smaller,
resulting in a steep density gradient at $r(\nu^{2}_{\mathrm{max}})$.

The bottom panel shows the maximum core sizes for the velocity-dependent and
velocity-independent cross-sections. It is visible that the maximum core size
is smaller than $r(v_{\mathrm{max}}^{2})$. Moreover, we can see that the
Knudsen number is increasing as a function of radius and is always much larger
than unity. Implying that the halo is always in the lmfp regime. For radii
smaller than $r_{\mathrm{s}}$, the corresponding heat conductivity
($\kappa^{\prime}_{\mathrm{lmfp}}$) is larger for the velocity-independent
cross-section. In contrast, $\kappa^{\prime}_{\mathrm{smfp}}$ is larger for
the velocity-dependent cross-section. If the cross-section is decreasing as a
function of velocity, smaller scattering velocities may play a more important
role compared to large velocities in the heat conduction than for velocity-
independent cross-sections (see also Sec. 3.1).

However, using the effective cross-section may eventually be problematic for
extreme velocity dependencies. Depending on the velocity of a DM particle, it
sees different distributions of relative velocities and thus has a mean free
path that depends on its velocity. Unfortunately, it is not understood how the
evolution in the lmfp regime could be derived from first principles. This
complicates a precise description of the heat conduction in the halo.

####  3.2.2 NFW Halo

name | type | $\sigma_{0}/m$ | $w$
---|---|---|---
|  | $[\mathrm{cm}^{2}\,\mathrm{g}^{-1}]$ | $[\mathrm{km}\,\mathrm{s}^{-1}]$
c0 | collisionless | $0$ | –
f10 | frequent | $10$ | –
r10 | rare | $10$ | –
f5e3w720 | frequent | $5\times 10^{3}$ | $720$
r5e3w720 | rare | $5\times 10^{3}$ | $720$
f2.5e5w180 | frequent | $2.5\times 10^{5}$ | $180$
r2.5e5w180 | rare | $2.5\times 10^{5}$ | $180$
Table 3: The cross-sections that we employed for simulating an NFW halo are
shown. The columns are the same as in Tab. 1.

We studied the core formation in an isolated NFW halo using various DM models.
These include velocity-independent cross-sections for fSIDM and rSIDM each
with $\sigma_{\mathrm{T}}/m=10.0\,\mathrm{cm}^{2}\mathrm{g}^{-1}$ and
velocity-dependent fSIDM and rSIDM cross-section with
$\sigma/m=5000.0\,\mathrm{cm}^{2}\mathrm{g}^{-1}$,
$w=720\,\mathrm{km}\,\mathrm{s}^{-1}$ and $\sigma/m=2.5\times
10^{5}\,\mathrm{cm}^{2}\mathrm{g}^{-1}$,
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$. The cross-sections and the
abbreviations we use for them are also shown in Tab. 3

For the NFW halo we use the same initial conditions as used by Fischer et al.
(2021a) for fig. 5. Our halo has a virial mass of
$10^{15}\,\mathrm{M_{\odot}}$, a scale radius of $300\,\mathrm{kpc}$ and a
density parameter of $\rho_{\mathrm{0}}\equiv 4\rho(r_{\mathrm{s}})=2.9\times
10^{6}\,\mathrm{M_{\odot}\,kpc}^{-3}$. The halo is sampled up to the virial
radius ($r_{\mathrm{vir}}=1626\,\mathrm{kpc}$) and resolved by $N=10^{6}$
particles. For the simulations we employ a gravitational softening length of
$\epsilon=0.56\,\mathrm{kpc}$.

Figure 4: We display the core size for an NFW halo, which we simulated with
different DM models. The abbreviations for the cross-sections are explained in
Tab. 3.  Figure 5: We show the same as in Fig. 3, but for the NFW simulations.
The the maximum core size of the velocity-dependent model refers to the run
with f2.5e5w180.

We measure the core size by fitting a cored NFW profile222In the literature
also other descriptions of a cored NFW profile exist (e.g. Read et al., 2016a,
b; Ray et al., 2022). The one we use corresponds to the one employed by Yang
et al. (2023a) for their parametric model of the evolution of a halo following
initially an NFW profile.. It is given by,

$\rho(r)=\frac{\rho_{0}}{(r^{4}+r^{4}_{\mathrm{core}})^{1/4}}\,\frac{r_{\mathrm{s}}}{\left(1+r/r_{\mathrm{s}}\right)^{2}}\,.$
(12)

For the fitting procedure we have $\rho_{0}$, $r_{\mathrm{s}}$ and
$r_{\mathrm{core}}$ as free parameters. We maximise a likelihood based on
Poisson statistics as described in sec. 4 of Fischer et al. (2021a).

The core sizes for different DM models are shown in Fig. 4. First, we consider
the cross-sections, f10, r10, f5e3w720 and r5e3w720. For the phase of the core
formation and the onset of core collapse up to $\approx 4\,\mathrm{Gyr}$ the
core sizes are very similar. Only the velocity-dependent rSIDM cross-section
yields slightly larger core sizes. Hence the momentum transfer cross-section
provides a good match between fSIDM and rSIDM in the given case. Only at later
stages of the halo evolution do differences between the models occur. When the
core size is almost zero, it seems that small-angle scattering slows down the
core collapse compared to isotropic scattering. These results are partially in
line with previous work. Yang & Yu (2022) found that a constant and velocity-
dependent cross-section behave qualitatively very similarly for most of the
halo evolution but differ at the late stages of the collapse phase. They also
found that the viscosity cross-section provides a better match between
different angular dependencies than the momentum transfer cross-section. In
the companion paper (Sabarish et al., prep), it is found that the viscosity
cross-section can indeed provide a reasonable, but not perfect match between
isotropic scattering and a very anisotropic cross-section in the fSIDM limit.
In contrast, for our setup with a much stronger cross-section the momentum
transfer cross-section provides a very good match regardless of the velocity
dependence. However, we should point out that the quality of the match depends
on the halo properties and the strength of the self-interactions, see, e.g.
fig. 9 from Fischer et al. (2022) (we show this result again in Sec. 4.2.2).
Here, one can see for the larger cross-section that the momentum transfer
cross-section match yields a larger effect of fSIDM on the central densities
of DM haloes at the high-mass end compared to rSIDM. For lower-mass haloes, it
changes and rSIDM has a stronger effect on the central halo density. As Yang &
Yu (2022) simulated NFW haloes with a mass of $M_{\mathrm{200}}\approx
10^{7}\,\mathrm{M_{\odot}}$ and a concentration of $c_{\mathrm{200}}\approx
20$ (for details see their tab. 1) they probed a different regime than we do
here. Hence, the quality of a matching procedure for the angular dependence
could depend on the halo properties and the strength of the self-interactions.
It is also important to note that the inner regions of our NFW halo are in the
smfp regime or close to it ($\mathrm{Kn}<1$) and not in the lmfp regime for
the velocity-independent cross-sections.

For the strongly velocity-dependent cross-section, i.e. the one with
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$, we find that the evolution differs
qualitatively from the ones with a weaker velocity-dependence. The results are
somewhat similar to the results for the Hernquist halo, the maximum core size
becomes larger and the collapse time longer compared to the core formation
time. However, the increase in the maximum core size is weaker compared to the
Hernquist halo. This could be because the cross-section we have simulated is
not as extremely velocity-dependent as for the Hernquist halo
($w=100\,\mathrm{km}\,\mathrm{s}^{-1}$ for the Hernquist halo and
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$ for the NFW halo). Note that the NFW
halo has a larger total mass and hence a larger velocity dispersion than the
Hernquist halo, such that the two simulations cannot be directly compared. But
when $w$ is compared to the typical scattering velocity of the halo, the
velocity dependence appears to be similar. In consequence, it is plausible
that the difference in maximum core size stems primarily from a different
reason such as the details of the density profile.

Analogous to the Hernquist halo we have computed the same quantities as in
Fig. 3, but for the NFW halo and show them in Fig. 5. In contrast to the
Hernqusit halo, we find that the central region of the halo has a Knudsen
number smaller than unity when simulated with the velocity-independent cross-
section and thus would be considered to be in the smfp regime. In addition,
the heat conductivity in the two regimes is more similar. But the Knudsen
number varies strongly with velocity dependence. As for the Hernquist halo
$\kappa^{\prime}_{\mathrm{smfp}}$ has a larger value in the case of the
velocity-dependent cross-section and $\kappa^{\prime}_{\mathrm{lmfp}}$ is
larger for the velocity-independent cross-section.

####  3.2.3 Discussion of isolated halo evolution

In this last part on isolated haloes, we discuss the physics driving their
evolution. During the evolution of the halo, the central velocity dispersion
is increasing and the effective strength of the self-interactions may change
according to the velocity dependence of the cross-section. An increasing
velocity dispersion implies higher relative velocities of the DM particles and
for a cross-section that decreases with velocity this leads to fewer
scatterings.

The halo may reach its maximum core size when the gradient of the velocity
dispersion has become zero. Afterwards, heat is only flowing outwards, which
leads to a shrinking density core and the gravothermal collapse of the halo.
While the density core is shrinking the central velocity dispersion is
increasing. Given this increase in velocity dispersion, one would expect that
the collapse is slowing down for a velocity-dependent cross-section compared
to a velocity-independent one. However, in our simulations, we do not find an
indication that the rate at which the density core is shrinking changes due to
the velocity dependence (see Fig. 4). Instead, we only found that the core
collapse time scale relative to the core formation time scale changes.

The evolution of the halo may not only be determined by the central region but
also by larger radii, at least radii up to $r(\nu^{2}_{\mathrm{max}})$ and a
bit beyond may play a crucial role. A core-collapse rate that is insensitive
to the velocity dependence might be caused by the relevant velocity dispersion
staying roughly constant. Indeed the velocity dispersion at larger radii is
less affected during the evolution and may play a crucial role in the core
collapse. Right from the beginning of the simulation, during core formation,
heat flows outward at radii larger than $r(\nu^{2}_{\mathrm{max}})$. This heat
flow takes place at velocities that are larger than in the central region of
the halo. In consequence, the ratio of heat inflow and outflow depends on the
velocity dependence of the scattering. For example, this is visible in the
core formation and core collapse times. They are set by heat inflow and
outflow.

The cross-sections we have simulated lead to roughly the same core formation
time. For strongly velocity-dependent cross-sections, less heat outflow takes
place during that time. This can result in a larger maximum core size as we
found for the Hernquist halo (see Fig. 2). The maximum core size depends on
the transition radius between heat inflow and outflow. Initially, this radius
is set by $r(\nu^{2}_{\mathrm{max}})$ but evolves according to the ratio of
heat in and outflow. As we found, this evolution is only significantly
affected by strongly velocity-dependent cross-sections.

Overall it becomes clear that if the scattering is velocity dependent, the
evolution of an isolated halo can change qualitatively. However, we do not
have a precise understanding of the physical mechanisms driving this
difference. How effective the heat outflow taking place in the lmfp regime
could depend on the gradient of the gravitational potential and the ability to
scatter particles to large velocities. It could be mainly the high-velocity
particles exceeding the escape velocity and carrying energy away that drive
the core collapse. In this context, the exact density profile may eventually
matter. For example the Hernquist and NFW profiles that we have investigated,
have a different slope in the outskirts. Implying a different gradient of the
gravitational potential. Further investigation is needed to fully understand
the evolution of isolated haloes.

##  4 Cosmological Simulations

We present our cosmological simulations in this section and show the results
we obtain. First, we describe the simulations, followed by the analysis of the
data. This includes many aspects such as the density and shape profiles of the
DM haloes and the abundance of satellites.

###  4.1 Simulations

We have run several simulations of CDM, rSIDM and fSIDM. For the SIDM models
we use two different velocity dependencies, namely
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$ and
$w=560\,\mathrm{km}\,\mathrm{s}^{-1}$. For each of them we have models that
differ in $\sigma_{0}$ by one order of magnitude. Our simulations are run with
fSIDM and a momentum transfer matched isotropic cross-section. The details of
the DM models are given in Tab. 4 and their velocity-dependence is plotted in
Fig. 6. Here, we also show the scattering velocities inside the centres of
haloes from three different mass bins, which we use in Section 4.2. The
velocities are indicated with a Maxwell–Boltzmann distribution that runs
logarithmically in velocity.

$f_{\mathrm{log}}(v_{\mathrm{scat}})=\sqrt{\frac{2}{\pi}}\frac{v_{\mathrm{scat}}^{3}}{a^{3}}\,e^{-\frac{v_{\mathrm{scat}}^{2}}{2\,a^{2}}}\quad\textnormal{with}\quad
a=\sqrt{2\,\nu^{2}}\,.$ (13)

The distribution of scattering velocities, $v_{\mathrm{scat}}$, depends on the
one-dimensional velocity dispersion, $\nu^{2}$, of the halo. In appendix F, we
put those DM models in the context of current observational constraints on the
strength of DM self-interactions.

name | type | $\sigma_{0}/m$ | $w$
---|---|---|---
|  | $[\mathrm{cm}^{2}\,\mathrm{g}^{-1}]$ | $[\mathrm{km}\,\mathrm{s}^{-1}]$
c0 | collisonless | $0.0$ | –
f0.1 | frequent | $0.1$ | –
r0.1 | rare | $0.1$ | –
f1 | frequent | $1.0$ | –
r1 | rare | $1.0$ | –
f10w180 | frequent | $10.0$ | $180$
r10w180 | rare | $10.0$ | $180$
f100w180 | frequent | $100.0$ | $180$
r100w180 | rare | $100.0$ | $180$
f0.3w560 | frequent | $0.3$ | $560$
r0.3w560 | rare | $0.3$ | $560$
f3w560 | frequent | $3.0$ | $560$
r3w560 | rare | $3.0$ | $560$
Table 4: The table shows the different cross-sections that we used for the
cosmological simulations. Analogously to Tab. 1, we use the same columns.
Note, that the simulations of the first five DM models have been presented by
Fischer et al. (2022).  Figure 6: We illustrate the cross-sections used for
our cosmological simulations. In blue, we show the velocity-independent cross-
sections from Fischer et al. (2022). The velocity-dependent cross-sections are
displayed in orange ($w=180\,\mathrm{km}\,\mathrm{s}^{-1}$) and purple
($w=560\,\mathrm{km}\,\mathrm{s}^{-1}$). In green, we indicate typical
scattering velocities. The Maxwell-Boltzmann distributions (see Eq. 13)
correspond to the scattering velocities in the centres of the haloes from the
three halo mass bins that we use in Sec. 4.2.

For the full box cosmological simulations, we use the same ICs as by Fischer
et al. (2022). They are similar to box 4 of the Magneticum simulations
333Magneticum: http://www.magneticum.org and have a comoving side length of
$48\,\mathrm{Mpc}h^{-1}$. The employed cosmological model is described by the
following parameters: $\Omega_{\mathrm{M}}=0.272$,
$\Omega_{\mathrm{\Lambda}}=0.728$, $h=0.704$, $n_{s}=0.963$ and
$\sigma_{8}=0.809$ (WMAP7, Komatsu et al., 2011). Further properties can be
found in Tab. 5.

Name | $l_{\mathrm{box}}$ | $N_{\mathrm{DM}}$ | $m_{\mathrm{DM}}$
---|---|---|---
| $(\mathrm{cMpc}\,h^{-1})$ |  | $(\mathrm{M_{\odot}}\,h^{-1})$
hr | 48 | $216^{3}$ | $8.28\times 10^{8}$
uhr | 48 | $576^{3}$ | $4.37\times 10^{7}$
Table 5: The table gives the different simulations we run. The first column
denotes the name, the second one the box size, the third one the number of
numerical DM particles and the last one the mass of the numerical DM
particles. Each setup was run with 8 different velocity-dependent cross-
sections as described in Tab. 4.

The DM haloes are identified using the friends-of-friends algorithm444A
description of the friends-of-friends algorithm can, for example, be found in
the work by More et al. (2011)., which is implemented in gadget-3. The mass of
a halo, $M$, is computed as the sum of the gravitationally bound particles.
The virial radius, $r_{\mathrm{vir}}$, and the virial mass,
$M_{\mathrm{vir}}$, are measured with the spherical-overdensity approach based
on the over-density predicted by the generalised spherical top-hat collapse
model (e.g. Eke et al., 1996). Here, $r_{\mathrm{vir}}$ is defined as the
radius at which the mean density becomes larger than the one of the top-hat
collapse model and $M_{\mathrm{vir}}$ is the mass inside $r_{\mathrm{vir}}$.

We use SubFind (Springel et al., 2001; Dolag et al., 2009), which is
implemented as part of gadget-3, to identify the substructure in the
simulation. Every halo contains at least one subhalo, which is the primary
subhalo located at the same position as the halo (determined by the location
of the most gravitationally bound particle of the halo). The primary subhalo
typically contains most of the particles that belong to the halo, but this is
not necessarily the case.

###  4.2 Results

In the following, we show the results of our cosmological simulations. The
simulation setup we used is described in Sec. 4.1. We begin with the surface
density of a massive halo (Section 4.2.1). Subsequently, we discuss the
density profiles of the haloes in Section 4.2.2 and continue with their shapes
(Section 4.2.3). We investigate the abundance of satellites (Section 4.2.4) as
well as their diversity in terms of the circular velocity (Sec. 4.2.5).
Finally, in Sec. 4.2.6, we study differences between frequent and rare self-
interactions in the context of velocity-dependent scattering.

Figure 7: The surface density of the fourth most massive system in our
simulation is shown. We cross-identified it among all simulations and show it
from the same perspective. We rotate the system such that for CDM the semi-
major axis is parallel to the $x$-axis and the semi-minor axis parallel to the
$y$-axis. We scale the axes in terms of $r_{1/2}$, the half mass radius of the
primary subhalo in the CDM simulation. The surface density is indicated with a
logarithmic colour scaling. We use the same for each panel. The abbreviation
of the cross-section is given in the lower left corner of each panel and the
detailed parameters can be looked up in Tab. 4.

####  4.2.1 Surface Density

In Fig. 7 we show the surface density of the same halo but in different DM
models. It is the fourth most massive halo ($M=9.3\times
10^{13}\,\mathrm{M_{\odot}}h^{-1}$) in our simulation and nicely illustrates
the effects of SIDM. They are most pronounced when comparing the two panels on
the left-hand side, as the fSIDM simulation of the two has relatively strong
self-interactions
($\sigma_{\mathrm{T}}/m=1.0\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$). Typical
effects of SIDM that can be seen here are the formation of a density core, the
rounder shape of haloes and the suppression of substructure. Many of the
satellites visible in the CDM run do not exist in the fSIDM run. However, in
the other SIDM runs shown here the suppression of the satellite abundance is
weaker. There exist even objects for which no counterpart in the CDM
simulation can be identified by eye. This is in particular the case for the
velocity-dependent cross-section shown in the right-hand side panels. In the
following sections, we quantify these self-interaction-induced changes in the
DM distribution.

####  4.2.2 Density Profiles

Figure 8: We show the median density profile for haloes from three different
mass bins. The results for the velocity-independent and velocity-dependent
cross-sections are displayed together. However, we show the results only for
fSIDM as the rSIDM results are similar. The density is plotted as a function
of the radius in units of the virial radius. The shaded regions indicate the
scatter among the haloes, and the range between the 25th and 75th percentiles
is displayed. The virial mass and the virial radius given in the panels
indicate the median of the corresponding mass bin from the CDM simulation. All
plots show the profiles for a redshift of $z=0$ and are produced from the full
cosmological box with the highest resolution. Note, we have used all
particles, not only those that belong to the halo as identified by SubFind.
Figure 9: The central density of the DM haloes is shown as a function of their
virial mass. We measure the central density as the mean density within a
radius of $0.01r_{\mathrm{vir}}$. In the left-hand side panel the simulations
with a velocity-independent cross-section are shown (reprint of fig. 9 of
Fischer et al., 2022). The middle panel gives the velocity-dependent
scattering with $w=560\,\mathrm{km}\,\mathrm{s}^{-1}$ and the right-hand side
panel displays the self-interactions with
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$. Individual systems are indicated by “+”
when evolved with the smaller cross-section. For the larger cross-section, we
use “$\times$” and the CDM case is marked by “$\downY$”. In addition, we
computed the mean of the distribution as a function of virial mass, shown by
the lines. The shaded regions give the corresponding standard deviation.

A quantity commonly measured for SIDM is the density profile of haloes. In
particular the formation of a central density core that is characterised by a
shallow gradient and a lower density compared to CDM (except O’Neil et al.,
2023). We have studied this in an idealised setup in Section 3. Within the
cosmological context this has been measured by various authors (e.g. Stafford
et al., 2020; Eckert, D. et al., 2022; Mastromarino et al., 2023) and used to
constrain the strength of DM self-interactions (see Appendix F).

We investigate the DM density profile for the haloes of our cosmological
simulations. In particular, we study the median density profile within three
halo mass bins. This is shown in Fig. 8, where we indicated the median virial
mass and virial radius of the haloes contained in the three mass bins. We show
all cross-sections we have simulated the ones with
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$ are shown in orange and the ones with
$w=560\,\mathrm{km}\,\mathrm{s}^{-1}$ are shown in purple. The small cross-
sections, i.e. the one with the smaller $\sigma_{0}/m$ for each $w$ show
hardly any core formation for the most massive haloes (left-hand panel). But
for the less massive haloes, the core size is increasing in terms of the
virial radius, $r_{\mathrm{vir}}$. This is a consequence of the relative
velocities between the DM particles being smaller for less massive systems. As
a result, the particles typically scatter at smaller relative velocities for
which the interaction strength is larger compared to high velocities (see also
Fig. 6).

In Fig. 9, we show the central density of the DM haloes as a function of their
virial mass. For the velocity-independent cross-section (left-hand panel), we
find that it is decreasing as a function of halo mass when self-interactions
are present. When considering the velocity-dependent runs it becomes clear
that the gradient with halo mass depends on the velocity-dependence of the
self-interactions. For $w=560\,\mathrm{km}\,\mathrm{s}^{-1}$ there is no or
only a weak trend with halo mass (middle panel). But for the
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$ cross-section (right-hand panel), the
central density is increasing with halo mass and thus the trend is opposite to
the simulations with a constant cross-section.

Note that we used the momentum transfer cross-section to match rSIDM and
fSIDM. If we would have used the viscosity cross-section, the fSIDM cross-
section would only have 2/3 of its value to correspond to the simulated rSIDM
cross-section. A detailed derivation of this factor has been presented by
(Sabarish et al., prep). This would imply larger central densities for the
fSIDM cross-sections. In consequence, it probably would often provide a better
matching. Except for haloes with masses lower than $M_{\mathrm{vir}}\approx
10^{13}\mathrm{M_{\odot}}$ and simulated with the strong and velocity-
independent scattering. Here the matching would become worse. It should be
noted, that not all haloes used in Fig. 9 are relaxed which makes the picture
more complicated.

####  4.2.3 Shapes

Figure 10: We show the median shape, $s=c/a$, of the DM haloes within three
mass bins as a function of the major semi-axis, $a$. Each panel displays a
different mass bin with its median mass being indicated. This figure is build
analogously to the density profiles in Fig. 8. The shaded regions indicate the
scatter among the haloes, and the range between the 25th and 75th percentiles
is displayed. We show it only for the collisionless DM and the strongest fSIDM
model of each velocity-dependence. In addition, we indicate at which radius
the shape sensitivity (Fischer & Valenzuela, 2023) for the 25th percentile
drops below a value of 25. This is indicative of a radius above which the
shape measurements are reliable. Note, in particular for CDM the presence of
satellites reduces the shape sensitivity.  Figure 11: The shape of the DM
haloes is shown as a function of their virial mass. The left-hand side panel
gives the results for the velocity-independent cross-sections (previously
shown in fig. 14 by Fischer et al., 2022). In the middle panel, we display the
results for the velocity-dependent scattering with
$w=560\,\mathrm{km}\,\mathrm{s}^{-1}$ and in the right-hand panel for
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$. This figure is built analogously to
Fig. 9.

A commonly studied property of DM haloes is their shape. This has for SIDM
been investigated by several authors (e.g. Peter et al., 2013; Sameie et al.,
2018; Robertson et al., 2019; Banerjee et al., 2020; Chua et al., 2020; Harvey
et al., 2021; Despali et al., 2022; Shen et al., 2022). DM self-interactions
significantly affect the shape of the haloes up to larger radii than the
density profile (Fischer et al., 2022). Furthermore, how large the affected
radii are depends on the strength of the self-interactions (Vargya et al.,
2022).

To compute the shapes of our simulated DM haloes we proceed as previously
described by Fischer et al. (2022). We compute the mass tensor of particles
within an ellipsoidal selection volume using their mass, $m$, and position,
$r$:

$\mathbf{M}_{ij}=\sum_{k}m_{k}r_{k,i}r_{k,j}$ (14)

Here, $k$ denotes a particle and $i$, $j$ are the coordinate indices. The
selection volume for the next iteration is determined by the eigenvalues and
eigenvectors of the mass tensor. We iterate until the shape of the selection
volume converges against the one inferred from the mass tensor. It is
important to note that shapes close to the centre of the haloes cannot be
measured accurately. The vanishing density gradient within the density core of
SIDM haloes renders the shape undefined (Fischer & Valenzuela, 2023).

In Fig. 10 we plot $s=c/a$ as a function of the semi-major axis, $a$, in units
of the virial radius. The semi-minor axis is denoted by $c$. In general, we
find that SIDM makes the haloes more round, as one would expect, and that
fSIDM and rSIDM are qualitatively very similar.

Moreover, we show the shape of the haloes as a function of mass in Fig. 11.
Here, we compute the shape from the innermost particles within a volume equal
to a sphere of radius $0.078r_{\mathrm{vir}}$. For CDM we find that haloes
become more ellipsoidal with increasing mass. This trend is well known in the
literature (e.g. Jing & Suto, 2002; Allgood et al., 2006; Muñoz-Cuartas et
al., 2011; Despali et al., 2013, 2014). This can change when including self-
interactions, especially for a velocity-independent cross-section. Here, the
effect of the self-interactions is increasing with halo mass (see left panel
of Fig. 11). However, for the most massive systems in our simulation we find
the haloes to become more elliptical even with SIDM. This might be due to few
objects which on average might be less relaxed than the ones at lower masses.
Given a velocity dependent cross-section haloes become more elliptical with
mass at the high-mass end. But the gradient is steeper compared to CDM as
self-interactions lead to rounder haloes at lower masses and at the high mass
end the shape becomes similar to CDM (middle and right panel).

####  4.2.4 Satellites

Figure 12: We show the number of satellites per logarithmic mass as a function
of their total mass relative to the virial mass of their host (upper panels).
In the lower panels we dispaly the ratio of the DM models to CDM. All panels
give the result of the 100 most massive groups in our full cosmological box.
The left-hand side panels shows the results for the velocity-independent
cross-sections (previously shown in fig. 6 of Fischer et al., 2022). The
middle panels gives the velocity-dependent self-interactions with
$w=560\,\mathrm{km}\,\mathrm{s}^{-1}$ and the right-hand side panels for
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$. All subhaloes, except for the primary
one, within a radius of $5\,r_{\mathrm{vir}}$ were considered. The results are
for a redshift of $z=0$. Note that the least resolved satellites used here
contain about 100 particles. Figure 13: For the 100 most massive haloes of our
simulations, we show the cumulative number of satellites per halo as a
function of radius (upper panels). We also give the ratio of the DM models to
CDM (lower panels). The left-hand side panel shows the results for the
velocity-independent cross-sections (previously shown in fig. 7 of Fischer et
al., 2022). The middle panel gives the velocity-dependent self-interactions
with $w=560\,\mathrm{km}\,\mathrm{s}^{-1}$ and the right-hand side panel for
$w=180\,\mathrm{km}\,\mathrm{s}^{-1}$. The results are shown for $z=0$ and
subhaloes were only considered if they are less massive than the primary
subhalo and more massive than $M>9.6\times
10^{10}\,\mathrm{M_{\odot}}\,h^{-1}$.

The properties of satellite systems are a promising probe for studies of DM.
Depending on the DM model, fewer or more satellites are predicted, and they
may differ in their density profiles. This has been studied in the context of
multiple DM models, including SIDM (e.g. Banerjee et al., 2020; Nadler et al.,
2020, 2021; Bhattacharyya et al., 2022).

In Fig. 12 we show the number of satellites per logarithmic mass as a function
of their mass in units of the virial mass of their host system. We find that
DM self-interactions can reduce the abundance of satellites, and the number of
less massive subhaloes is stronger affected than the more massive satellites.
Moreover, the momentum-transfer-matched frequent self-interactions lead to a
stronger suppression than the isotropic scattering (as previously described
for a constant cross-section in Fischer et al., 2022). All this seems to be
independent of the velocity dependence. Interestingly, the difference between
fSIDM and rSIDM is shrinking for the strong velocity dependence. For the
velocity-independent simulations (left-hand side panel) and the mildly
velocity-dependent runs ($w=560\,\mathrm{km}\,\mathrm{s}^{-1}$, middle panel),
the stronger rSIDM cross-section has a similar effect to the weak fSIDM cross-
section. But for the strongly velocity-dependent run
($w=180\,\mathrm{km}\,\mathrm{s}^{-1}$, right-hand panel), the strong rSIDM
cross-section is no longer similar to the weak fSIDM one but closer to the
strong fSIDM one. Hence, we find a strong velocity dependence to reduce the
differences between cross-sections with different angular dependencies.

The difference between rSIDM and fSIDM may mainly arise from host-satellite
scattering as those interactions take place with a preferred direction and
thus are far from an equilibrium state. Also, these interactions significantly
contribute to the suppression of the satellite abundance (e.g. Zeng et al.,
2022). To understand the reduced difference between rSIDM and fSIDM it is
important to note, that the host-satellite interactions take place at higher
velocities than the scatterings within the satellite between its particles.
Consequently, a velocity-dependent cross-section can reduce the host-satellite
scattering compared to the satellite-satellite interactions and thus reduce
the difference between rSIDM and fSIDM.

Figure 14: We show the circular velocity at $3.5\,\mathrm{kpc}$ for satellites
with a mass of at least $\approx 4.9\times
10^{10}\,\mathrm{M_{\odot}}\,h^{-1}$. We consider all satellites that are not
the primary subhalo. The lines indicate the mean and the shaded regions the
standard deviation for the corresponding DM models. This is analogous to Fig.
9, as well as the markers.

In addition, we find that the suppression of the satellite abundance for the
mildly velocity-dependent cross-sections (middle panel) is less strong than
for the other two velocity dependencies. We would not have expected this
difference in strength from the density profiles that we show in Sec. 4.2.2.
Though, there is a velocity scale at which the mildly velocity-dependent
cross-sections are weaker than the corresponding ones with a different
velocity-dependence (see Fig. 6). Interestingly, this becomes even more
pronounced when computing the effective cross-section introduced by Yang & Yu
(2022), see Appendix F. Given that the host-satellite scattering, which drives
the suppression of the satellite abundance, takes preferentially place in this
velocity regime, it could explain the different strengths of the satellite
suppression.

In Fig. 13 we display the number of satellites as a function of the distance
to their host in units of the host’s virial radius. The upper panels show the
cumulative number of satellites and the lower panels display the ratio to CDM.
We note that the ratios at small distances are subject to a considerable
amount of noise as they are computed from a small number of satellites. Here,
we find again that self-interactions can suppress the number of satellites.
The inner ones are more affected than the distant ones and frequent self-
interactions lead to a stronger suppression than rare scattering if the same
momentum transfer cross-sections are compared. This is well visible for the
velocity-independent cross-sections in the left-hand panel. The simulations
with frequent self-interactions show roughly a reduction in the number of
satellites twice as large as for the corresponding simulations with rare self-
interactions. As in Fig. 12 we find that the difference between rSIDM and
fSIDM becomes less for the strongest velocity-dependence
($w=180\,\mathrm{km}\,\mathrm{s}^{-1}$).

####  4.2.5 Diversity of satellites

One of the small-scale issues is the diversity problem. It usually refers to
the variation between the rotation curves of galaxies (e.g. Kamada et al.,
2017; Ren et al., 2019; Zentner et al., 2022). To study their diversity, we
focus on the circular velocity at a radius of $3.5\,\mathrm{kpc}$ instead of
looking at the full profile. The velocity at $3.5\,\mathrm{kpc}$ is sensitive
to the core formation or core collapse. In Fig. 14 we show the circular
velocity at that radius for satellites more massive than $\approx 4.9\times
10^{10}\,\mathrm{M_{\odot}}\,h^{-1}$ as a function of their mass. Note that we
consider all subhaloes identified by SubFind satellites if they are not a
primary subhalo (see Sec. 2.1).

For the velocity-independent cross-sections (left-hand panel of Fig. 14) we
find that self-interactions decrease the circular velocity at
$3.5\,\mathrm{kpc}$. This corresponds to the formation of a density core. For
the larger cross-sections the circular velocity is lower, i.e. the density
core is larger. Basically, the same applies for the cross-sections with
$w=560\,\mathrm{km}\,\mathrm{s}^{-1}$ (middle panel of Fig. 14). But it is
noticeable that the most massive subhaloes experience less suppression of
$v_{\mathrm{circ}}$ in the inner region. This is simply a consequence of the
velocity dependence, as the DM particles in the more massive subhaloes have
higher typical relative velocities. For the cross-section with the strong
velocity dependence ($w=180\,\mathrm{km}\,\mathrm{s}^{-1}$), we find
qualitatively different results. For the more massive subhaloes, we find the
suppression of the circular velocity as for the other simulations. But on
average, the least massive objects show an increase in circular velocity for
the stronger cross-sections compared to CDM. Moreover, the distribution of
values for the circular velocity is broader compared to CDM. The other cross-
sections do not show such a significant increase in diversity.

When comparing the results for rSIDM and fSIDM, we do not find a clear
qualitative difference arising from the typical scattering angle of the self-
interactions. In contrast, the momentum transfer cross-section provides a
matching that is not far off but surprisingly accurate for the velocity-
dependent cross-sections.

The diversity of rotation curves has been studied a lot with SIDM, and it has
been shown that self-interactions can create more diverse density profiles. In
particular, low-mass objects have been studied. There are several papers that
studied MW-like satellites and dwarf galaxies (e.g. Creasey et al., 2017;
Zavala et al., 2019; Correa et al., 2022; Lovell & Zavala, 2023). It has been
found in DMO simulations that cross-sections with a strong velocity dependence
can even trigger core collapse within satellites (e.g. Turner et al., 2021;
Yang et al., 2023c; Nadler et al., 2023). Especially for satellites the core
collapse can be enhanced by tidal stripping (e.g. Kahlhoefer et al., 2019;
Nishikawa et al., 2020). This is in line with our finding of more compact
objects at low masses for our strongly velocity-dependent cross-sections.

####  4.2.6 Frequent vs. rare self-interactions

Figure 15: We show how the DM model affects the maximum circular velocity in
the satellites and the host’s central density. We have cross-identified the
haloes in the different DM runs. The lines connect the same halo, i.e.
indicated how the properties of a halo change when varying the cross-section.
The shown haloes are among the most massive ones, the details of the selection
criterion are explained in the text.

Finally, we want to investigate how the different DM models affect the
satellites of our most massive haloes. Previously, we found that fSIDM can
lead to a stronger suppression of the number of satellites than rSIDM does
(Fischer et al., 2022). Identifying such differences is crucial to constrain
the angular dependence of dark matter self-interactions. In contrast to our
previous work, we investigate the maximum circular velocity in the satellites
here. But show the number of satellites in Appendix E.

We cross-identify the haloes and their satellites among the simulations based
on their particles. As we start from the same initial conditions, we can match
the haloes with the same particles identified based on their unique
identification numbers. To evaluate how well two haloes match we make use of
the gravitational potential at the particle’s location. Particles at a lower
gravitational potential are stronger weighted to find the best matching
analogue. Given a list of the halo particles sorted according to how deep they
sit in the gravitational potential, starting with the one at the lowest
potential, we compute weights for them. These weights are given as

$w_{i}=\left(\frac{1}{i+1}\right)^{\alpha}\,.$ (15)

Note, here we assume the first list index to be $i=0$. The parameter $\alpha$
allows for different weightings, we use $\alpha=0.8$. In practice, we compute
the weight for the CDM run only. This is because we use the CDM haloes as a
benchmark and ask how well the SIDM haloes match them. The quality of a
potential match is given by the sum of the weights $w_{i}$ for the particles
that the CDM halo and the SIDM halo have in common.

For the analysis, we do not consider all haloes but apply different selection
criteria. Firstly, the hosts and their satellites should be well-resolved. We
consider only the 13th most massive haloes and limit the selection further by
requiring that we are able to match at least five satellites with a minimum
mass of $9.6\times 10^{10}\,\mathrm{M_{\odot}}\,h^{-1}$ (2200 particles).
Furthermore, we require the haloes to be relaxed. Here, we assume a halo to be
relaxed if the centre of mass and the most bound particle of the primary
subhalo are separated by not more than $10\%$ of the virial radius. In
addition, we tested a further limitation by excluding haloes based on the
ratio of the halo and primary subhalo mass. However, in practice, this did not
exclude any halo. At least when we have required that the primary subhalo does
not contain less than $75\%$ of the halo mass.

In Fig. 15, we display our results for how the central halo densities
correlate with the relative change of the maximum circular velocity in the
satellites. We show the average relative change multiplied by the average
maximum circular velocity in the CDM satellites. Here, we use the maximum
velocity as computed by SubFind. It is given by the maximum of the circular
velocity, $v_{\mathrm{circ}}=\sqrt{\mathrm{G}\,M(<r)/r}$, in radial distance,
$r$, from the centre of the subhalo.

We find the maximum circular velocity in the satellites altered by the DM
self-interactions. For the velocity-independent scattering it typically
decreases with increasing cross-section. This implies that the satellites are
less concentrated. In contrast, a velocity-dependent cross-section can also
lead to a larger value for the maximum circular velocity. Whether this is the
case or not depends in our model on the parameter $w$, i.e. how strongly
velocity-dependent the scattering is. It is worth pointing out that our
selection criterion of subhaloes above a mass threshold that we can match
might favourably pick subhaloes that have become more concentrated due to the
velocity-dependent self-interactions. Thus the increase in maximum circular
velocity may not be representative of all the subhaloes.

We find that frequent self-interactions tend to lead to a smaller maximum
circular velocity than rare scattering. For the larger cross-sections we have
simulated, we find that the maximum circular velocity for rare self-
interactions compared to frequent ones is increased for the typical system
(median) by $\approx 8\%$ (velocity-independent), $\approx 2\%$
($w=560\,\mathrm{km}\,\mathrm{s}^{-1}$), and $\approx 1\%$
($w=180\,\mathrm{km}\,\mathrm{s}^{-1}$). This means that the difference
between fSIDM and rSIDM decreases for our simulations with stronger velocity
dependence. Hence, this is in line with our finding of a qualitative
difference for the abundance of satellites in Section 4.2.4. However, the
difference we find here might also largely be due to the fact that the
stronger velocity-dependent cross-section we study has a weaker effect on
massive haloes. For example, this becomes visible when comparing the central
densities. In consequence, the reduced qualitative difference between large-
and small-angle scattering might be better visible from Fig. 12. But here we
can see that not only for a constant cross-section the angular dependence
matters but also for strongly velocity-dependent self-interactions even if the
subhaloes are becoming more compact on average.

We note that the analysis above is not based on a larger statistical sample
and thus the exact numbers may change. But we expect the qualitative trend to
be the same. It is also worth pointing out that the less massive satellites
might be affected more strongly by the self-interactions (see Fig. 12) and
thus differences between models are larger for them. Hence, this should be
followed up with simulations with a much higher spatial resolution.

##  5 Discussion

In this section, we discuss the assumptions and limitations of our simulations
as well as the implications of our results. We begin with technical
considerations and end by discussing what the next steps for a follow-up study
may look like.

In contrast to our previous work (Fischer et al., 2022) we explored velocity-
dependent cross-sections. We found that simulating those interactions requires
a separate time step criterion (i.e. different from the one of Fischer et al.,
2021b). Especially cross-section with a strong velocity-dependence, i.e. a
small value for $w$ (see Eq. 2), can be computationally very expensive
compared to a velocity-independent cross-section with a similar effective
cross-section. A more detailed discussion of building a time step criterion
can be found in Appendix B.

When measuring the core sizes in Sec. 3.2, we found that the resulting fit is
surprisingly sensitive to the optimisation method. This may limit the
comparability of core sizes inferred by different authors. In particular,
Correa et al. (2022) describe in their appendix B, that results on the
evolution of the core size differ in terms of the maximum cores size in the
literature.

The results of our cosmological simulations depend on the algorithms employed
to identify haloes and their substructure. For this task, we used the build-in
module SubFind (Springel et al., 2001; Dolag et al., 2009). There exist a
number of codes that are capable of identifying substructure (e.g. Knollmann &
Knebe, 2009; Maciejewski et al., 2009; Tweed, D. et al., 2009; Behroozi et
al., 2012; Han et al., 2017; Elahi et al., 2019). These codes use different
algorithms and are known to give somewhat different results Knebe et al.
(2013). In consequence, our results could change a bit when employing a
different substructure finder.

In this paper, we aimed to understand how a velocity dependence of the self-
interactions affects differences arising from the angular dependence of the
cross-section. Very anisotropic cross-sections are typically expected to be
velocity-dependent (e.g. Buckley & Fox, 2010; Loeb & Weiner, 2011; Bringmann
et al., 2017). It is known that fSIDM and rSIDM differ mainly in systems that
are far from equilibrium, such as mergers (Fischer et al., 2021a) and the
abundance of satellites Fischer et al. (2022). The evolution of those systems
is governed by multiple velocity scales, where typically the larger velocity
scale is the one that is mainly responsible for differences arising from the
angular dependence of the self-interactions. Consequently, the difference
becomes less when the self-interactions at large velocities are suppressed due
to velocity-dependent scattering. We found this for the abundance of
satellites. In consequence, it could be interesting to probe less massive
systems for distinguishing rSIDM and fSIDM as the velocity dependence could be
weaker. At least in the model employed in our study, a system with typical
velocities smaller than $w$ would only experience a weak velocity dependence
(see Eq. 2). The relevant mass scales for the cross-sections we simulated are
visible from the effective cross-section as a function of mass shown in
Appendix F.

Despite our studies of satellites, it is worth mentioning that very
anisotropic cross-sections have been mainly studied in the context of merging
galaxy clusters (e.g. Kahlhoefer et al., 2014; Harvey et al., 2015; Fischer et
al., 2023; Wittman et al., 2023). At about the pericentre passage, such cross-
sections can give rise to an effective drag force decelerating the DM
component and creating an offset between the galaxies and the DM. Cross-
sections that are velocity dependent and strongly anisotropic, have not been
studied in the context of such mergers yet. Only a Bullet Cluster-like system
has been simulated by Robertson et al. (2017b) using a velocity-dependent
anisotropic cross-section, but it does not fall within the limit of fSIDM.
Studying merging systems with velocity-dependent fSIDM is crucial to
understand their power to constrain such models and is the subject of a
companion paper (Sabarish et al., prep).

Our simulations are all DM-only. On the one hand, it allows us to understand
the qualitative differences between DM models better compared to simulations
including further physical processes. But on the other hand, it limits the
possibility to compare the results to observations and derive constraints on
the cross-section. Consequently, the next step would be to include baryonic
physics, i.e. run hydrodynamical simulations. Several authors have found that
taking baryons into account can reduce the differences between collisonless
and self-interacting DM and thus would mitigates constraints derived from DM-
only studies (e.g. Fry et al., 2015; Despali et al., 2022; Sirks et al., 2022;
Mastromarino et al., 2023). SIDM can be more responsive to the baryon
distribution than CDM in Milky Way-mass galaxies (e.g. Sameie et al., 2018,
2021). In the presence of baryons, effects from SIDM can even be reversed – at
least for a fraction of the haloes. It has been shown that for galaxies with
Milky Way-like masses and above the interplay of baryons and self-interactions
can lead to cuspier density profiles than in CDM (e.g. Despali et al., 2019;
Rose et al., 2022). In principle, baryons could also affect the ability to
constrain the angular dependence with the abundance of satellites.

Aside from constraining the angular dependence, one would like to have a
procedure to compare the effect of SIDM with different angular dependencies.
This would allow to transfer constraints between models that differ in their
typical scattering angle. Yang & Yu (2022) introduced the effective cross-
section for this purpose, where the angular matching is based on the viscosity
cross-section. However, the quality of the matching may depend on the physical
system, i.e. how relaxed the system is. But not only on this, as we found the
momentum transfer cross-section can at least for some setups provide an
excellent match (see Fig. 4) excluding that the viscosity cross-section does
as well. However, this does not contradict the viscosity cross-section
providing a better match usually. But it implies that the matching is more
complicated and may depend on the properties of the astrophysical system. It
may matter how strong the self-interactions are and whether the system evolves
in the smfp or lmfp regime. In the latter one gravity plays an important role
between two consecutive scattering events (assuming an isotropic cross-
section) and thus may make the evolution of the halo and the matching of
different angular dependencies sensitive to the details of the density
profile.

##  6 Conclusions

In this paper, we have studied SIDM with velocity-dependent scattering,
considering isotropic cross-sections and strongly forward-enhanced ones. For
accurate modelling of velocity-dependent self-interactions, we introduced a
new time-step criterion and enhanced the performance with an improved
parallelisation scheme. To learn about qualitative differences arising from
the velocity dependence, we first simulated the thermalization problem, a
simple test problem without gravity. Secondly, we studied the evolution of the
density profile of isolated haloes including Hernquist and NFW profiles. For
the remainder of the paper, we focused on cosmological simulations and
investigated the qualitative differences between the DM models concerning the
velocity and angular dependence of the self-interactions. Our most important
results can be summarized as follows:

  * •

We found that velocity-dependent self-interactions lead to a slower population
of the high-velocity tail of the Maxwell-Boltzmann distribution during
thermalisation due to the suppressed cross-section at high velocities.

  * •

The evolution of the density profile of isolated haloes is qualitatively
affected by the velocity dependence, i.e. it is not self-similar. This can
lead to a longer collapse time relative to the core formation time and a
larger maximum core size. However, we found a significant difference between
velocity-independent and velocity-dependent cross-sections only for strong
velocity dependencies, i.e. when $w$ is much smaller than the typical
scattering velocity.

  * •

The velocity dependence of the self-interactions controls whether the central
density of haloes is increasing or decreasing as a function of halo mass.

  * •

Given a strong velocity dependence (small $w$), frequent self-interactions can
diversify the density profile similar to an isotropic cross-section. We found
that the two angular dependencies can create haloes that are less compact as
well as haloes that are more compact at the same subhalo mass.

  * •

A strong velocity dependence of the cross-section, i.e. a small value of $w$
can reduce the differences between fSIDM and rSIDM regarding the abundance of
satellites.

The simulations we conducted were DM-only and allowed us to understand
phenomenological differences arising from the velocity dependence of DM
scattering. Our results may be instructive for more detailed studies of
qualitative differences between SIDM models and helpful in designing more
sophisticated simulations that include baryonic matter and additional physics
such as cooling, star formation, AGN, and associated feedback mechanisms.
Undertaking such a study to learn about the chances to discriminate between
rSIDM and fSIDM when baryonic physics is taken into account, is the subject of
forthcoming work.

## Acknowledgements

We thank all participants of the Darkium SIDM Journal Club for helpful
discussions. MSF thanks Lucas Kimmig for a fruitful discussion. This work is
funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe”
– 390833306, Germany’s Excellence Strategy – EXC-2094 “Origins” – 390783311
and the Emmy Noether Grant No. KA 4662/1-2. Antonio Ragagnin acknowledge
support from the grant PRIN-MIUR 2017 WSCC32. KD acknowledges support by the
COMPLEX project from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme grant agreement
ERC-2019-AdG 882679. The simulations have been carried out on the computing
facility Hummel (HPC-Cluster 2015) of the University of Hamburg and the
computing facilities of the Computational Center for Particle and Astrophysics
(C2PAP). Preprint numbers: DESY-23-154, TTP23-041.

Software: NumPy (Harris et al., 2020), Matplotlib (Hunter, 2007), SciPy
(Virtanen et al., 2020)

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author. In addition, some of the data can be retrieved from our
webpage: https://www.darkium.org.

## References

  * Adhikari et al. (2022) Adhikari S., et al., 2022, Astrophysical Tests of Dark Matter Self-Interactions (arXiv:2207.10638)
  * Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
  * Andrade et al. (2021) Andrade K. E., Fuson J., Gad-Nasr S., Kong D., Minor Q., Roberts M. G., Kaplinghat M., 2021, MNRAS, 510, 54
  * Balberg et al. (2002) Balberg S., Shapiro S. L., Inagaki S., 2002, ApJ, 568, 475
  * Banerjee et al. (2020) Banerjee A., Adhikari S., Dalal N., More S., Kravtsov A., 2020, J. Cosmology Astropart. Phys., 2020, 024
  * Behroozi et al. (2012) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2012, ApJ, 762, 109
  * Bhattacharyya et al. (2022) Bhattacharyya S., Adhikari S., Banerjee A., More S., Kumar A., Nadler E. O., Chatterjee S., 2022, ApJ, 932, 30
  * Bonafede et al. (2011) Bonafede A., Dolag K., Stasyszyn F., Murante G., Borgani S., 2011, MNRAS, 418, 2234
  * Bringmann et al. (2017) Bringmann T., Kahlhoefer F., Schmidt-Hoberg K., Walia P., 2017, Phys. Rev. Lett., 118, 141802
  * Buckley & Fox (2010) Buckley M. R., Fox P. J., 2010, Phys. Rev. D, 81, 083522
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Burkert (2000) Burkert A., 2000, ApJ, 534, L143
  * Chu et al. (2019) Chu X., Garcia-Cely C., Murayama H., 2019, Phys. Rev. Lett., 122, 071103
  * Chua et al. (2020) Chua K. T. E., Dibert K., Vogelsberger M., Zavala J., 2020, MNRAS, 500, 1531
  * Colin et al. (2002) Colin P., Avila-Reese V., Valenzuela O., Firmani C., 2002, ApJ, 581, 777
  * Correa (2021) Correa C. A., 2021, MNRAS, 503, 920
  * Correa et al. (2022) Correa C. A., Schaller M., Ploeckinger S., Anau Montel N., Weniger C., Ando S., 2022, MNRAS, 517, 3045
  * Creasey et al. (2017) Creasey P., Sameie O., Sales L. V., Yu H.-B., Vogelsberger M., Zavala J., 2017, MNRAS, 468, 2283
  * Davé et al. (2001) Davé R., Spergel D. N., Steinhardt P. J., Wandelt B. D., 2001, ApJ, 547, 574
  * Dehnen (2001) Dehnen W., 2001, MNRAS, 324, 273
  * Despali et al. (2013) Despali G., Tormen G., Sheth R. K., 2013, MNRAS, 431, 1143
  * Despali et al. (2014) Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208
  * Despali et al. (2019) Despali G., Sparre M., Vegetti S., Vogelsberger M., Zavala J., Marinacci F., 2019, MNRAS, 484, 4563–4573
  * Despali et al. (2022) Despali G., Walls L. G., Vegetti S., Sparre M., Vogelsberger M., Zavala J., 2022, MNRAS, 516, 4543
  * Dodelson & Widrow (1994) Dodelson S., Widrow L. M., 1994, Phys. Rev. Lett., 72, 17
  * Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
  * Downing & Oman (2023) Downing E. R., Oman K. A., 2023, MNRAS, 522, 3318
  * Dutton & Macciò (2014) Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
  * Eckert, D. et al. (2022) Eckert, D. Ettori, S. Robertson, A. Massey, R. Pointecouteau, E. Harvey, D. McCarthy, I. G. 2022, A&A, 666, A41
  * Eke et al. (1996) Eke V. R., Cole S., Frenk C. S., 1996, MNRAS, 282, 263
  * Elahi et al. (2019) Elahi P. J., Cañas R., Poulton R. J. J., Tobar R. J., Willis J. S., Lagos C. d. P., Power C., Robotham A. S. G., 2019, Publications of the Astronomical Society of Australia, 36
  * Elbert et al. (2015) Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., Oñorbe J., Peter A. H. G., 2015, MNRAS, 453, 29
  * Fischer & Valenzuela (2023) Fischer M. S., Valenzuela L. M., 2023, A&A, 670, A120
  * Fischer et al. (2021a) Fischer M. S., Brüggen M., Schmidt-Hoberg K., Dolag K., Kahlhoefer F., Ragagnin A., Robertson A., 2021a, MNRAS, 505, 851
  * Fischer et al. (2021b) Fischer M. S., Brüggen M., Schmidt-Hoberg K., Dolag K., Ragagnin A., Robertson A., 2021b, MNRAS, 510, 4080
  * Fischer et al. (2022) Fischer M. S., Brüggen M., Schmidt-Hoberg K., Dolag K., Kahlhoefer F., Ragagnin A., Robertson A., 2022, MNRAS, 516, 1923
  * Fischer et al. (2023) Fischer M. S., Durke N.-H., Hollingshausen K., Hammer C., Brüggen M., Dolag K., 2023, MNRAS
  * Fry et al. (2015) Fry A. B., et al., 2015, MNRAS, 452, 1468–1479
  * Genina et al. (2018) Genina A., et al., 2018, MNRAS, 474, 1398
  * Gilman et al. (2021) Gilman D., Bovy J., Treu T., Nierenberg A., Birrer S., Benson A., Sameie O., 2021, MNRAS, 507, 2432
  * Gopika & Desai (2023) Gopika K., Desai S., 2023, Physics of the Dark Universe, 42, 101291
  * Governato et al. (2012) Governato F., et al., 2012, MNRAS, 422, 1231
  * Han et al. (2017) Han J., Cole S., Frenk C. S., Benitez-Llambay A., Helly J., 2017, MNRAS, 474, 604–617
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Harvey et al. (2015) Harvey D., Massey R., Kitching T., Taylor A., Tittley E., 2015, Science, 347, 1462
  * Harvey et al. (2019) Harvey D., Robertson A., Massey R., McCarthy I. G., 2019, MNRAS, 488, 1572
  * Harvey et al. (2021) Harvey D., Chisari N. E., Robertson A., McCarthy I. G., 2021, MNRAS, 506, 441–451
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hu et al. (2000) Hu W., Barkana R., Gruzinov A., 2000, Phys. Rev. Lett., 85, 1158
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Ibe & Yu (2010) Ibe M., Yu H.-B., 2010, Physics Letters B, 692, 70
  * Jing & Suto (2002) Jing Y. P., Suto Y., 2002, ApJ, 574, 538
  * Kahlhoefer et al. (2014) Kahlhoefer F., Schmidt-Hoberg K., Frandsen M. T., Sarkar S., 2014, MNRAS, 437, 2865
  * Kahlhoefer et al. (2017) Kahlhoefer F., Schmidt-Hoberg K., Wild S., 2017, J. Cosmology Astropart. Phys., 08, 003
  * Kahlhoefer et al. (2019) Kahlhoefer F., Kaplinghat M., Slatyer T. R., Wu C.-L., 2019, J. Cosmology Astropart. Phys., 2019, 010
  * Kamada et al. (2017) Kamada A., Kaplinghat M., Pace A. B., Yu H.-B., 2017, Phys. Rev. Lett., 119, 111102
  * Kaplinghat et al. (2016) Kaplinghat M., Tulin S., Yu H.-B., 2016, Phys. Rev. Lett., 116, 041302
  * Kaplinghat et al. (2019) Kaplinghat M., Valli M., Yu H.-B., 2019, MNRAS, 490, 231
  * Kasselmann (2021) Kasselmann L., 2021, N-body simulations of self-interacting dark matter with a velocity-dependent cross section
  * Knebe et al. (2013) Knebe A., et al., 2013, MNRAS, 435, 1618
  * Knollmann & Knebe (2009) Knollmann S. R., Knebe A., 2009, ApJS, 182, 608–624
  * Kochanek & White (2000) Kochanek C. S., White M., 2000, ApJ, 543, 514
  * Koda & Shapiro (2011) Koda J., Shapiro P. R., 2011, MNRAS, 415, 1125
  * Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
  * Kummer et al. (2019) Kummer J., Brüggen M., Dolag K., Kahlhoefer F., Schmidt-Hoberg K., 2019, MNRAS, 487, 354
  * Loeb & Weiner (2011) Loeb A., Weiner N., 2011, Phys. Rev. Lett., 106
  * Lovell & Zavala (2023) Lovell M. R., Zavala J., 2023, MNRAS, 520, 1567
  * Lynden-Bell & Eggleton (1980) Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483
  * Maciejewski et al. (2009) Maciejewski M., Colombi S., Springel V., Alard C., Bouchet F. R., 2009, MNRAS, 396, 1329–1348
  * Martizzi et al. (2013) Martizzi D., Teyssier R., Moore B., 2013, MNRAS, 432, 1947
  * Mastromarino et al. (2023) Mastromarino C., Despali G., Moscardini L., Robertson A., Meneghetti M., Maturi M., 2023, MNRAS, 524, 1515
  * Moore et al. (1998) Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ, 499, L5
  * More et al. (2011) More S., Kravtsov A. V., Dalal N., Gottlöber S., 2011, ApJS, 195, 4
  * Muñoz-Cuartas et al. (2011) Muñoz-Cuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584
  * Nadler et al. (2020) Nadler E. O., Banerjee A., Adhikari S., Mao Y.-Y., Wechsler R. H., 2020, ApJ, 896, 112
  * Nadler et al. (2021) Nadler E. O., Banerjee A., Adhikari S., Mao Y.-Y., Wechsler R. H., 2021, ApJ, 920, L11
  * Nadler et al. (2023) Nadler E. O., Yang D., Yu H.-B., 2023, An SIDM Solution to the Extreme Diversity of Low-mass Halo Properties (arXiv:2306.01830)
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Nishikawa et al. (2020) Nishikawa H., Boddy K. K., Kaplinghat M., 2020, Phys. Rev. D, 101, 063009
  * O’Neil et al. (2023) O’Neil S., et al., 2023, MNRAS, 524, 288
  * Oman et al. (2019) Oman K. A., Marasco A., Navarro J. F., Frenk C. S., Schaye J., Benítez-Llambay A. r., 2019, MNRAS, 482, 821
  * Outmezguine et al. (2023) Outmezguine N. J., Boddy K. K., Gad-Nasr S., Kaplinghat M., Sagunski L., 2023, MNRAS, 523, 4786
  * Peirani et al. (2017) Peirani S., et al., 2017, MNRAS, 472, 2153
  * Peter et al. (2013) Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
  * Pineda et al. (2016) Pineda J. C. B., Hayward C. C., Springel V., Mendes de Oliveira C., 2016, MNRAS, 466, 63
  * Planelles et al. (2013) Planelles S., Borgani S., Fabjan D., Killedar M., Murante G., Granato G. L., Ragone-Figueroa C., Dolag K., 2013, MNRAS, 438, 195
  * Pontzen & Governato (2012) Pontzen A., Governato F., 2012, MNRAS, 421, 3464
  * Rahimi et al. (2023) Rahimi E., Vienneau E., Bozorgnia N., Robertson A., 2023, J. Cosmology Astropart. Phys., 2023, 040
  * Rasia et al. (2015) Rasia E., et al., 2015, ApJ, 813, L17
  * Ray et al. (2022) Ray T. S., Sarkar S., Shaw A. K., 2022, J. Cosmology Astropart. Phys., 2022, 011
  * Read & Gilmore (2005) Read J. I., Gilmore G., 2005, MNRAS, 356, 107
  * Read et al. (2016a) Read J. I., Agertz O., Collins M. L. M., 2016a, MNRAS, 459, 2573
  * Read et al. (2016b) Read J. I., Iorio G., Agertz O., Fraternali F., 2016b, MNRAS, 462, 3628
  * Ren et al. (2019) Ren T., Kwa A., Kaplinghat M., Yu H.-B., 2019, Phys. Rev. X, 9, 031020
  * Robertson et al. (2017a) Robertson A., Massey R., Eke V., 2017a, MNRAS, 465, 569
  * Robertson et al. (2017b) Robertson A., Massey R., Eke V., 2017b, MNRAS, 467, 4719
  * Robertson et al. (2018) Robertson A., et al., 2018, MNRAS, 476, L20–L24
  * Robertson et al. (2019) Robertson A., Harvey D., Massey R., Eke V., McCarthy I. G., Jauzac M., Li B., Schaye J., 2019, MNRAS, 488, 3646–3662
  * Robertson et al. (2020) Robertson A., Massey R., Eke V., Schaye J., Theuns T., 2020, MNRAS, 501, 4610–4634
  * Rocha et al. (2013) Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81
  * Roper et al. (2023) Roper F. A., Oman K. A., Frenk C. S., Benítez-Llambay A., Navarro J. F., Santos-Santos I. M. E., 2023, MNRAS, 521, 1316
  * Rose et al. (2022) Rose J. C., Torrey P., Vogelsberger M., O’Neil S., 2022, MNRAS, 519, 5623
  * Sabarish et al. (prep) Sabarish V. M., Brüggen M., Schmidt-Hoberg K., Fischer M. S., Kahlhoefer F., in prep., arXiv e-prints
  * Sagunski et al. (2021) Sagunski L., Gad-Nasr S., Colquhoun B., Robertson A., Tulin S., 2021, J. Cosmology Astropart. Phys., 2021, 024
  * Sameie et al. (2018) Sameie O., Creasey P., Yu H.-B., Sales L. V., Vogelsberger M., Zavala J., 2018, MNRAS, 479, 359
  * Sameie et al. (2021) Sameie O., et al., 2021, MNRAS, 507, 720
  * Shen et al. (2022) Shen X., Brinckmann T., Rapetti D., Vogelsberger M., Mantz A., Zavala J., Allen S. W., 2022, MNRAS, 516, 1302
  * Shi et al. (2021) Shi Y., Zhang Z.-Y., Wang J., Chen J., Gu Q., Yu X., Li S., 2021, ApJ, 909, 20
  * Silk (2017) Silk J., 2017, ApJ, 839, L13
  * Silverman et al. (2022) Silverman M., Bullock J. S., Kaplinghat M., Robles V. H., Valli M., 2022, MNRAS, 518, 2418
  * Sirks et al. (2022) Sirks E. L., Oman K. A., Robertson A., Massey R., Frenk C., 2022, MNRAS, 511, 5927
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Stafford et al. (2020) Stafford S. G., Brown S. T., McCarthy I. G., Font A. S., Robertson A., Poole-McKenzie R., 2020, MNRAS, 497, 3809
  * Tsai et al. (2022) Tsai Y.-D., McGehee R., Murayama H., 2022, Phys. Rev. Lett., 128, 172001
  * Tulin & Yu (2018) Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1
  * Tulin et al. (2013) Tulin S., Yu H.-B., Zurek K. M., 2013, Phys. Rev. D, 87, 115007
  * Turner et al. (2021) Turner H. C., Lovell M. R., Zavala J., Vogelsberger M., 2021, MNRAS, 505, 5327
  * Tweed, D. et al. (2009) Tweed, D. Devriendt, J. Blaizot, J. Colombi, S. Slyz, A. 2009, A&A, 506, 647
  * Vargya et al. (2022) Vargya D., Sanderson R., Sameie O., Boylan-Kolchin M., Hopkins P. F., Wetzel A., Graus A., 2022, MNRAS, 516, 2389
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2012) Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740
  * Vogelsberger et al. (2014) Vogelsberger M., Zavala J., Simpson C., Jenkins A., 2014, MNRAS, 444, 3684
  * Wittman et al. (2023) Wittman D., Stancioli R., Finner K., Bouhrik F., van Weeren R., Botteon A., 2023, ApJ, 954, 36
  * Yang & Yu (2022) Yang D., Yu H.-B., 2022, J. Cosmology Astropart. Phys., 2022, 077
  * Yang et al. (2023a) Yang D., Nadler E. O., Yu H.-B., Zhong Y.-M., 2023a, A Parametric Model for Self-Interacting Dark Matter Halos (arXiv:2305.16176)
  * Yang et al. (2023b) Yang S., Du X., Zeng Z. C., Benson A., Jiang F., Nadler E. O., Peter A. H. G., 2023b, ApJ, 946, 47
  * Yang et al. (2023c) Yang D., Nadler E. O., Yu H.-B., 2023c, ApJ, 949, 67
  * Zavala et al. (2013) Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20
  * Zavala et al. (2019) Zavala J., Lovell M. R., Vogelsberger M., Burger J. D., 2019, Phys. Rev. D, 100, 063007
  * Zeng et al. (2022) Zeng Z. C., Peter A. H. G., Du X., Benson A., Kim S., Jiang F., Cyr-Racine F.-Y., Vogelsberger M., 2022, MNRAS, 513, 4845
  * Zentner et al. (2022) Zentner A., Dandavate S., Slone O., Lisanti M., 2022, J. Cosmology Astropart. Phys., 2022

##  Appendix A Parallelisation

The improved MPI parallelization orders communication with the aim to reduce
waiting time. Therefore, a separate communication list is computed at the cost
of additional overhead.

The communication list is computed in seven steps:

  1. 1.

Every MPI rank already knows how many particles it will receive from and send
to other processes. It creates a list with all its communications to other
processes, here the sum of the particles that are exchanged with each process
(sum of send and receive) is stored. This number will later be used to assign
individual communications a priority.

  2. 2.

Each process sorts the previously created list according to priority.

  3. 3.

All processes exchange their communications list. In turn, every process has
all communication lists. We note that the length of these lists does scale
with less than $N^{2}$ ($N$ being the number of MPI ranks) if no communication
between all processes is required (this is typically the case).

  4. 4.

Every process builds a single list with all communications. Starting with the
highest priority elements of all individual lists.

  5. 5.

The list from the previous step is sorted according to priority.

  6. 6.

The overall communication plan is built by every process to avoid additional
communication. It contains several communication steps. At each step, several
pairs do their communication and also the scattering computations of the
exchanged particles. One process can be assigned only one pairwise
communication per step. The plan is built by trying to fulfil the highest
priorities first. We start with the first communication step and try to fill
in the communications sorted after priority. If a communication does not fit
in a step because a corresponding process is already busy, it is queued and
retried for the next step. We turn to the next step when the queue for the
communication that should be assigned for the current step is empty or when
all processes are already assigned a communication for the current step.

  7. 7.

Each process extracts its communication schedule from the overall
communication plan.

For each pair of processes we consider, the first sends the particles and the
second receives them and does the computation. When the computation has
finished, the particles are sent back. Subsequently, the communication is done
in the other direction, i.e. the other task will compute the scattering of the
particles.

Figure 16: We show the execution time as a function of simulated time for the
Magneticum Box 4 with high resolution (hr) and ultra-high resolution (uhr).
The simulations are DMO. We also show the time spend on the fSIDM-related
computations, which make up the majority of the computational costs. Note,
that we do not use adaptive gravitational softening. Consequently, the
computation of the kernel sizes is counted as fSIDM-related. In grey, we show
the ratio of execution time between the old and new parallelisation.

We have run a performance test for Box 4 of the Magneticum simulations with DM
only using the high-resolution (hr) and ultra-high-resolution (uhr) initial
conditions (the same as in Fischer et al., 2022). The results using 32 and 128
MPI ranks are displayed in Fig. 16. For this test, we used only the MPI
parallelisation and did not make use of an OpenMP parallelisation that exists
for other parts of the code. Clearly, the improved parallelisation leads to a
significant speed-up. It is visible that we save with the improvements more
than $15\%$ of the computational costs This parallelisation has been used for
some of the simulations described in this paper.

##  Appendix B Time step criterion discussion

In the following, we discuss the thoughts behind the construction of a time
step criterion for SIDM in further detail. Here, we aim to build a time step
criterion that ensures that all or almost all interactions take place at a
sufficiently small time step. This differs from the approach taken by
Vogelsberger et al. (2012), which considered the local velocity dispersion.
Instead, we are concerned about the full velocity distribution, i.e. which
relative velocities a particle actually sees.

Previously we introduced a time step criterion for velocity-independent self-
interactions (Fischer et al., 2021b). That time step criterion estimates the
time step based on the maximum velocity that a particle experienced in the
previous time step555Strictly speaking, we use $\sigma_{\mathrm{T}}/m\,\Delta
v$, but for a velocity-independent cross-section $\sigma_{\mathrm{T}}/m$ is a
constant and thus the same for all interactions. and the maximal possible
kernel overlap, $\Lambda_{ii}$666Actually, we did not directly compute
$\Lambda_{ii}$, but used $h_{i}^{-3}$ as an estimate ($h$ denotes the kernel
size).. Given that the neighbour number, $N_{\mathrm{ngb}}$, is sufficiently
large, the velocities that a particle has seen in the previous time step allow
us to roughly estimate the maximum velocity it may experience in the next time
step. In contrast, for a velocity-dependent cross-section, the relevant
velocity is not the maximal velocity that a particle may see but how close it
gets to the velocity for which the interaction probability, or generally
speaking, $\sigma_{\mathrm{T}}/m\,\Delta v$, becomes maximal. Actually, this
is also what the time step criterion by (Fischer et al., 2021b) tries to
estimate when used for a velocity-dependent cross-section. But in contrast to
a constant cross-section, it is much harder to estimate this for a velocity-
dependent cross-section based on the velocities a particle has seen, as the
probability of seeing a relative velocity close to $v_{e}$ (see Eq. 4) might
be small. Our tests showed that we would often overestimate the time step.
This problem can be circumvented by directly using $v_{e}$, instead of making
estimates based on what a particle has seen in the previous time step as
described in Sec. 2.3. Hence, it is possible to build a timestep criterion
that guarantees for each particle pair the interaction probability or drag
force is sufficiently small.

Lastly, we want to explain why one should not directly use the interaction
probabilities a particle has encountered in the previous time step. The
disadvantage is that large interaction probabilities are less likely to be
seen by a numerical particle than the relevant velocities (for a constant
cross-section this would be the maximum velocity). This is because in most
cases the kernel overlap, $\Lambda_{ij}$, is small. Or in other words, the
probability of having a numerical particle pair with a relative velocity close
to the relevant velocity and a large kernel overlap is smaller than having
only a relative velocity close to the relevant velocity.

##  Appendix C Comoving Integration Test

Figure 17: The cosmic deceleration problem in terms of canonical momentum is
shown. It is simulated from $a=0.5$ to $a=1.0$ with $122500$ particles in a
cubic box with a comoving side length of $1400\,\mathrm{kpc}\,h^{-1}$. The
total mass is $22.8465\times 10^{10}\,\mathrm{M_{\odot}}\,h^{-1}$,
corresponding to a comoving density of
$83.26\,\mathrm{M_{\odot}}\,\mathrm{kpc}^{-3}\,h^{2}$. The initial snapshot
velocity of the test particle is $100\,\mathrm{kpc}\mathrm{Gyr}^{-1}$, which
corresponds to an initial canonical momentum of
$35.35534\,\mathrm{kpc}\,\mathrm{Gyr}^{-1}$. The particles are evolved with a
cross-section of $\sigma_{0}/m_{\chi}=7\times
10^{7}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$,
$w=10.0\,\mathrm{km}\,\mathrm{s}^{-1}$ and the SIDM kernel sizes are computed
using $N_{\mathrm{ngb}}=64$.

For testing the implementation of velocity-dependent self-interactions, we
introduce and use a new test problem in this Appendix. The test problem is
very similar to the one we have used by Fischer et al. (2022). A single
particle is travelling within an expanding space through a background density.
This background is at rest (zero canonical momentum) and has no density
gradient. For this test problem, we only consider the drag force which
decelerates the particle but do not re-add the energy as described in section
2.2 of Fischer et al. (2021a). We do not take any further physics into
account, i.e. run the test problem without gravity. Hence, we expect the test
particle to be decelerated over time. We calculate semi-analytically how the
canonical momentum of the test particle evolves over time and compare the
results from the simulation to it. This is shown in Fig. 17. Note, in the
absence of self-interactions the canonical momentum would stay constant over
the cosmic expansion.

Further, we want to point out that this test problem is more susceptible to
numerical errors than a typically fSIDM simulation. The interaction between a
pair of numerical particles does not change their relative velocity. This
makes the pairwise interaction in some sense time-implicit and more stable.
However, if we compute the drag force only and do not re-add the energy, as
done for the test problem, we break this. This is true for the test problem
when conducted with a velocity-independent cross-section (see appendix A by
Fischer et al., 2022) too. But in contrast, the velocity dependence makes it
even more unstable. Assuming that the test particle is slightly faster than it
is supposed to be, one would expect the drag force to be stronger (velocity-
independent cross-section) or weaker (velocity-dependent cross-section) than
it is supposed to experience. The first case would suppress the deviation but
the second enhances it. In the opposite case where the particle is slower than
supposed, one finds again that the velocity-dependent cross-section tends to
increase the deviation. In consequence, the test problem we show is quite
unstable. However, in general, this depends on how strong the velocity
dependency is, in our model specified by $\alpha$.

Overall, we find that the test simulation agrees sufficiently enough with the
prediction and we can conclude the implementation of velocity-dependent self-
interactions works as supposed.

##  Appendix D Convergence of density profile

name | $N_{\mathrm{high\,res}}$ | $m_{\mathrm{DM}}$
---|---|---
|  | $[\mathrm{M_{\odot}}\,h^{-1}]$
1x | $\sim 4.51\times 10^{4}$ | $8.3\times 10^{8}$
10x | $\sim 4.52\times 10^{5}$ | $8.3\times 10^{7}$
25x | $\sim 1.13\times 10^{6}$ | $3.3\times 10^{7}$
250x | $\sim 1.13\times 10^{7}$ | $3.3\times 10^{6}$
2500x | $\sim 1.13\times 10^{8}$ | $3.3\times 10^{5}$
Table 6: The properties of the zoom-in simulations we use for the convergence
test are given. We provide the name of the simulation, the number of particles
in the highly resolved region ($N_{\mathrm{high\,res}}$) and the mass of the
high-resolution particles ($m_{\mathrm{DM}}$). All simulations share the same
initial conditions but with different resolutions. Figure 18: We show the
density profile of the most massive subhalo of our zoom-in simulations. We
give the profile for a CDM simulation (left-hand panel), which is the same as
given by Fischer et al. (2022). The other two panels show the result of
velocity-dependent fSIDM simulations. The different colours indicate different
resolutions. This allows us to see that the profiles are converging for
increasing resolution. In the highest-resolved run, the subhalo is represented
by $\approx 2.3\times 10^{6}$ particles.

In this Appendix, we study the convergence of our simulations with velocity-
dependent self-interactions. To do so we run zoom-in simulations of the same
object but with varying resolutions. The zoom-in region is selected from a
large box with a comoving side length of $1\,\mathrm{Gpc}\,h^{-1}$ and its
most massive halo has a virial mass of $\sim 8.8\times
10^{11}\,\mathrm{M_{\odot}}\,h^{-1}$. Several publications (e.g. Planelles et
al., 2013; Rasia et al., 2015) made use of this box for zoom-in initial
conditions and it was first described by Bonafede et al. (2011). We run the
simulation with two cross-sections,
$\sigma_{0}/m=1.0\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$,
$w=360.0\,\mathrm{km}\,\mathrm{s}^{-1}$ and
$\sigma_{0}/m=5.0\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$,
$w=360.0\,\mathrm{km}\,\mathrm{s}^{-1}$. The different resolutions we
simulated are described in Tab. 6.

In Fig. 18, we show the density profile for the most massive halo of the zoom-
in simulations. We can see that the density profile converges for
collisionless DM (upper panel), but also when velocity-dependent self-
interactions are present (middle and lower panel).

##  Appendix E Frequent vs. rare self-interactions

Figure 19: We show the central density (upper panels) and the shape (lower
panels) of the host halo as a function of the number of their satellites for
different DM models. The velocity-independent cross-sections are shown in the
left-hand side panels (this has previously been shown in fig. 16 by Fischer et
al., 2022). The middle panels give the results for the models with
$w=560\,\mathrm{km}\,\mathrm{s}^{-1}$. And cross-sections with the strongest
velocity dependence ($w=180\,\mathrm{km}\,\mathrm{s}^{-1}$) are displayed in
the right-hand side panels.

Here we show the central density of the host halo as a function of the number
of satellites (Fig. 19). The three most massive haloes are displayed as
previously done in fig. 16 by Fischer et al. (2022). We find that the frequent
self-interactions independent of $w$ reduce the number of satellites stronger
than rare scattering when comparing them at levels of the same central host
density (upper panels) or the roundness of the host’s shape (lower panels).
However, one halo for the $w=180\,\mathrm{km}\,\mathrm{s}^{-1}$ poses an
exception.

##  Appendix F SIDM constraints

Figure 20: We show constraints for a velocity-independent cross-section
together with the fSIDM models that we simulated. This is given in terms of
the effective cross-section, $\sigma_{\mathrm{eff}}$, (Yang & Yu, 2022) as a
function of the virial DM halo mass. Constraints on the self-interaction
strength obtained by various authors are shown. The colours of our SIDM models
correspond to the ones shown in Fig. 6. Note, our rSIDM models are $1/3$
weaker than the fSIDM cross-sections when compared in terms of
$\sigma_{\mathrm{eff}}$.

In Fig. 20, we show constraints on the strength of DM self-interactions
together with our SIDM models. Here, we compute the effective cross-section as
introduced by Yang & Yu (2022) (see also Eq. 8). This requires an estimate of
an effective velocity dispersion, which we compute from a given virial mass,
$M_{\mathrm{vir}}$. To do so, we use the halo mass-concentration relation
given by Dutton & Macciò (2014). With the obtained concentration parameter,
$c$, we infer the maximum velocity dispersion, $\nu^{2}_{\mathrm{max}}$. For
the effective velocity dispersion we employ
$\sigma^{\mathrm{eff}}_{\mathrm{1D}}=0.9\times\nu_{\mathrm{max}}$. We choose
the factor of $0.9$ as it provides a good match for our isolated NFW
simulation shown in Section 3.2.2. This concerns the match of the velocity-
independent cross-sections with the ones that are described by
$w=720\,\mathrm{km}\,\mathrm{s}^{-1}$ and $\sigma_{0}/m=5\times
10^{3}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$. However, we have to note that the
viscosity cross-section-like matching for the angular dependence in
$\sigma_{\mathrm{eff}}$ does not provide a match as good as the one from the
momentum transfer cross-section in this particular case. If we would have used
the viscosity cross-section for the matching the isotropic cross-section would
have $3/2$ of the strength we obtained from the momentum transfer matching
while leaving the fSIDM cross-section unchanged.

The constraints shown in Fig. 20 stem from measures of different effects that
SIDM has on the distribution of DM. This includes the formation of a density
core (Sagunski et al., 2021; Andrade et al., 2021; Correa, 2021; Shi et al.,
2021; Eckert, D. et al., 2022; Gopika & Desai, 2023), oscillations of the
brightest cluster galaxy (Harvey et al., 2019), and the shapes of the haloes
(Peter et al., 2013; Despali et al., 2022).