# Gravothermal Catastrophe in Resonant Self-interacting Dark Matter Models

Vinh Tran  0009-0003-6068-6921 Department of Physics and Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA<EMAIL_ADDRESS>Daniel Gilman  0000-0002-5116-7287
Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL
60637, USA  Brinson Prize Fellow  Mark Vogelsberger  0000-0001-8593-7692
Department of Physics and Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Xuejian Shen
0000-0002-6196-823X Department of Physics and Kavli Institute for Astrophysics
and Space Research, Massachusetts Institute of Technology, Cambridge, MA
02139, USA  Stephanie O’Neil  0000-0002-7968-2088 Department of Physics and
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, MA 02139, USA  Xinyue Zhang  0000-0002-9502-5600
Department of Industrial Engineering and Operations Research, Columbia
University, New York, NY 10027, USA

###### Abstract

We investigate a self-interacting dark matter (SIDM) model featuring a
velocity-dependent cross section with an order-of-magnitude resonant
enhancement of the cross section at $\sim 16\,{\rm km}\,{\rm s}^{-1}$. To
understand the implications for the structure of dark matter halos, we perform
N-body simulations of isolated dark matter halos of mass $\sim 10^{8}{\,\rm
M_{\odot}}$, a halo mass selected to have a maximum response to the resonance.
We track the core formation and the gravothermal collapse phases of the dark
matter halo in this model and compare the halo evolving with the resonant
cross section with halos evolving with velocity-independent cross sections. We
show that dark matter halo evolution with the resonant cross section exhibits
a deviation from universality that characterizes halo evolution with velocity-
independent cross sections. The halo evolving under the influence of the
resonance reaches a lower minimum central density during core formation. It
subsequently takes about $20\%$ longer to reach its initial central density
during the collapse phase. These results motivate a more detailed exploration
of halo evolution in models with pronounced resonances.

keywords

††preprint: APS/123-QED

##  I Introduction

To date, dark matter (DM) remains a purely astrophysical phenomenon (see [1]
for a review), with the only positive evidence for its existence coming from
cosmological observations [2]. In particular, some of the strongest bounds on
the nature of DM from astrophysics derive from inferences of the properties of
DM halos, the concentrations of DM that host galaxy clusters and galaxies in
their centers. In many classes of DM theory, the details of the particle
physics model alter the abundance and internal structure of halos, leading to
constraints on DM particle physics from inferences of the abundance and
internal structure of halos [3, 4].

The concordance cosmological model includes cold dark matter (CDM) which, by
definition, behaves as a collisionless fluid. In the last two decades, self-
interacting dark matter (SIDM) has emerged as a viable alternative to CDM
[e.g. 5, 6, 7, 8, 9, 10]. SIDM refers to a class of theory in which the DM has
a non-zero self-interaction cross section (see [11] for a review). The
exchange of energy and momentum through scattering in SIDM causes a temporal
evolution in the internal structure of halos, whereas CDM predicts a static
internal structure well-approximated by a Navarro-Frenk-White profile [12,
hereafter NFW profile] for halos at all times. The evolution of halo internal
structure begins with a period of core formation [e.g. 10, 13, 14, 15], during
which time heat flows into the center of the halo. Following core formation,
the direction of heat flow eventually reverses, causing a runaway contraction
of the density profile referred to as ‘gravothermal catastrophe’, or core
collapse. This phenomenon was initially studied in the context of dense star
clusters [16]. Several decades later, Balberg  _et al._ [17] (see also Koda
and Shapiro [18]) pointed out that DM self-interactions provide a more
efficient mechanism to transfer heat throughout halo profiles than purely
gravitational scattering. Core collapse represents the final state of isolated
self-interacting DM halos [e.g. 19, 20, 21, 22]. Inside galactic halos,
various processes such as tidal stripping and scattering between subhalo and
host halo particles can delay or accelerate the onset of collapse in subhalos
[23, 24, 25, 26, 27, 28]. More complicated core collapse behaviors are found
under the influence of baryons [e.g. 29] or with dissipative DM self-
interactions [e.g. 30, 31, 32, 33].

Core formation and collapse in SIDM theories result in a variety of observable
consequences for the properties of dwarf galaxies [e.g. 34, 35, 36, 21, 37].
In addition, cored and collapsed halos produce distinct gravitational lenses
signatures that enable studies of SIDM from observations of strongly-lensed
quasars [38, 39], as well as galaxy-galaxy strong lenses [40]. The most
definitive existing constraints on SIDM theories come from observations of
galaxy clusters, which place stringent upper limits on the self-interaction
cross section $\sigma<0.1\rm{cm^{2}}\ \rm{g^{-1}}$ [41, 42, 43], effectively
requiring that the DM behave as a collisionless fluid on halo mass scales
$\sim 10^{14}M_{\odot}$. The total mass of a collapsed halo sets the velocity
scale for the scattering process, as the typical speed of a DM particle in a
virialized structure increases with the mass of the halo. Thus, the
constraints on the self-interaction cross section from galaxy clusters place
stringent upper limits on the self-interaction cross section at relative
velocities $v\sim 1000\ \rm{km}\ \rm{s^{-1}}$. For a velocity-independent
scattering cross section, this upper limit would rule an SIDM model capable of
altering the observable properties of DM halos on galactic and sub-galactic
scales. However, if the SIDM cross section has a velocity dependence that
increases the cross section at low relative velocities, heat transfer in low-
mass halos can become significantly enhanced, and core formation and collapse
occur within the age of the Universe [e.g. 44, 45]. Velocity-dependent cross
sections that become amplified at low speeds are a common feature of
scattering mediated by light force carriers with Yukawa interaction potentials
(for a detailed review of this topic, see Tulin  _et al._ [9]).

The current modeling of core collapse assigns a DM halo a single effective
cross section as a function of the velocity-dependent cross section and the
halo’s central velocity dispersion. From this effective cross section, one
obtains a characteristic timescale for which the evolution of SIDM halos is
self-similar, or universal [20, 22, 46]. This approach enables a predictive
framework for modeling core collapse, enabling the determination of the time
of core collapse as a multiple of the collapse timescale given the properties
of the SIDM cross section and structural parameters of a halo [e.g. 38]. One
potential complication for this relatively simple picture of halo evolution in
SIDM is that for a DM halo of a fixed total mass, the typical relative
velocity between particles changes as a function of radius, and thus the
scattering inside a given halo occurs over a range of velocities. While some
numerical simulations have verified the universality of halo evolution for
some velocity-dependent cross sections [20], in this work we consider a class
of SIDM models with velocity-dependent cross sections that change
significantly across the range of relative velocities typical of particles
inside a halo.

We consider a class of SIDM models in which a resonance in the cross section
causes an order-of-magnitude enhancement to the scattering cross section near
a particle relative velocity. Resonances refer to an enhancement or
suppression of the cross section that results from the formation of bound
states between DM particles and attractive potential. Resonances can manifest
as single pronounced peaks, multiple peaks, or as suppression or enhancement
of the cross section strength across a range of speeds, depending on the
details of the particle physics model [9, 47, 39]. For this class of model,
the cross section amplitude can signifcantly change with radius inside a halo,
complicating the relatively simple physical picture in which a halo evolves
under the influence of a single effective cross section. This class model of
model has received relatively limited study. Ref. [48] examined the evolution
of a halo analytically using the gravothermal fluid formalism, while Ref. [28]
studied a resonant model in a numerical simulation, but did not distinguish
its structural evolution from models with simpler velocity-dependence.

To gain a detailed physical understanding of the evolution of self-interacting
DM halos with pronounced resonances in the cross section, we conduct idealized
N-body simulations of isolated halos with a resonant cross section. We pick a
halo mass expected to maximize the effect of the resonance on the halo
evolution, and compare the evolution of the halo with the resonant cross
section with the evolution of halos evolving under the influence of velocity-
independent cross sections. We examine whether existing treatments of SIDM
halo evolution and universality hold for this class of model.

This paper is organized as follows: Section II gives an overview of the
particle physics model and the calculation of the resonant cross section
considered in this work. Section III discusses the calculation of two
effective scattering cross sections we will consider when interpreting the
evolution of the halo with a resonant cross section. In Section IV, we review
the details of our N-body simulations. Section V presents the results of the
simulations. We summarize our findings and make concluding remarks in Section
VI.

##  II Particle physics model

We consider an SIDM cross section that results from interactions between DM
particles of mass $m_{\chi}$ through a force carrier of mass $m_{\phi}$. We
consider the case of an attractive Yukawa potential $\left(c=\hbar=1\right)$

$V\left(r\right)=\frac{-\alpha}{r}\exp{\left(-rm_{\phi}\right)}$ (1)

where $\alpha$ determines the strength of the potential. While no exact
closed-form solutions exist for the differential scattering cross section
corresponding to $V\left(r\right)$, in certain regimes of parameter space
various approximations give analytic and closed-form solutions that match the
exact cross section to high precision. For example, for weak potentials
$\left(\alpha m_{\chi}/m_{\phi}<<1\right)$, one obtains an approximate
analytic expression for the differential scattering cross section using the
Born approximation

$\frac{{\rm d}\sigma}{{\rm
d}\Omega}=\frac{\sigma_{0}}{\left(w^{2}+v^{2}\sin^{2}\left(\theta/2\right)\right)^{2}}$
(2)

where $v$ represents relative velocity,
$\sigma_{0}\equiv\alpha^{2}m_{\chi}^{2}/m_{\phi}^{4}$ and $w\equiv
m_{\phi}/m_{\chi}$.

Outside the range of validity of the Born approximation, one must compute the
differential scattering cross section using partial wave analysis. In this
case, the differential scattering cross section is given by

$\frac{{\rm d}\sigma}{{\rm
d}\Omega}=\frac{1}{k^{2}}\Big{|}\sum_{\ell=0}^{\infty}\left(2\ell+1\right)e^{i\delta_{l}}\sin\delta_{l}P_{\ell}\left(\cos\theta\right)\Big{|}^{2}$
(3)

where $\delta_{\ell}$ represents the phase shifts of partial waves labeled by
their angular momentum quantum number $\ell$, $k=m_{\chi}v/2$ is the momentum
of the DM particle, and $P_{\ell}\left(x\right)$ are Legendre polynomials of
order $\ell$.

It has become common practice to average over the angular dependence of the
differential scattering cross section when predicting the internal structure
of DM halos in SIDM. This simplification can yield accurate predictions for
the evolution of halo profiles provided one chooses a suitable proxy for ${\rm
d}\sigma/{\rm d}\Omega$. As shown by Yang and Yu [20], the viscosity transfer
cross section $\sigma_{V}=\int\frac{{\rm d}\sigma}{{\rm
d}\Omega}\sin^{2}\theta{\rm d}\theta$, for which $\sigma_{V}$ is given by [49]

$\sigma_{V}=\frac{4\pi}{k^{2}}\sum_{\ell=0}^{\ell_{\rm{max}}}\frac{\left(\ell+1\right)\left(\ell+2\right)}{2\ell+3}\sin^{2}\left(\delta_{\ell+2}-\delta_{\ell}\right),$
(4)

allows one to predict the structural evolution of a halo without explicitly
accounting for the angular dependence of the scattering cross section in
N-body simulations. This represents a considerable simplification for the
numerical scattering implementation, and for obtaining analytical predictions
for the structural evolution of the halo. We truncate the summation at a value
of $\ell_{\rm{max}}$ for which $\sigma_{V}$ converges to within $1\%$

For an attractive potential given by Equation 1 in the quantum regime
$\left(k/m_{\phi}<1\right)$, the formation of bound states with the potential
gives rise to resonances in the cross section. A resonance refers to the
enhancement or suppression of the scattering cross section across a range of
relative velocities. The velocities where resonances in the cross section
occur depend on the values of $\alpha$ and $m_{\phi}/m_{\chi}$ for the
potential in Equation 1. Due to the computational cost of solving the full
scattering problem and then running N-body simulations that implement the
resulting cross sections, SIDM models with significant resonances have not
been extensively studied in the literature. As a result, we lack a detailed
understanding of structure formation in this class of SIDM theory, and it is
unclear whether analytical treatments of halo evolution in SIDM calibrated for
other velocity-dependent cross sections [e.g. 20, 46] apply to this class of
model.

To evaluate the cross section in the resonant regime, we solve the Schrödinger
equation to compute the phase shifts $\delta_{\rm{\ell}}$, and use Equation 4
to compute the cross section. To calculate $\delta_{\ell}$ we use an auxiliary
function for the Schrödinger equation [47]

$\dfrac{\partial\delta_{\ell}\left(r\right)}{\partial
r}=-k\,m_{\chi}\,r^{2}\,V\left(r\right)\,\left[\cos\left(\delta_{\ell}\left(r\right)\right)\,j_{\ell}\left(kr\right)-\sin\left(\delta_{\ell}\left(r\right)\right)\,n_{\ell}\left(kr\right)\right]^{2},$
(5)

where $j_{\ell}\left(x\right)$ and $n_{\ell}\left(x\right)$ are spherical
Bessel functions of the first and second kind of order $\ell$. We obtain a
numerical solution for $\delta_{\ell}\left(r\right)$ using Mathematica with
the initial condition $\delta_{\ell}\left(0\right)=0$, and obtain the phase
shift by taking the limit
$\delta_{\ell}=\lim_{r\rightarrow\infty}\delta_{\ell}\left(r\right)$.

For the analysis presented in this work, we consider the resonant cross
section represented by the shaded region in the left-hand panel Figure 1. This
cross section was also considered by Ref. [39] in a strong-lensing analysis of
multiply-imaged quasars, and has $m_{\rm{\chi}}=38.1\ \rm{GeV}$,
$m_{\rm{\phi}}=6.8\ \rm{MeV}$, and $\alpha=1.6\times 10^{-3}$. The amplitude
of the cross section, as indicated by the axis labels on the right-hand side
of the figure, reaches $\sim 100\ \rm{cm^{2}}\ \rm{g^{-1}}$ at $v\sim 16\
\rm{km}\ \rm{s^{-1}}$, a feature that arises primarily from constructive
interference between the $\ell=3$ and $\ell=1$ partial waves. To explore the
consequences of this model for structure formation, we run a N-body simulation
of a DM halo evolving under the influence of this resonant cross section, and
track its evolution from core formation until core collapse.

##  III Analytic predictions for the evolution of halo structure

Figure 1: Particle collisions distribution as functions of relative velocity
(left) and cross section per DM mass (right) for halos of masses $10^{7}$
(cyan), $10^{8}$ (red), $10^{9}$ (purple), and $10^{10}$ (blue) $\rm
M_{\odot}$. The resonant cross section model is shown by the yellow profile in
the left-hand panel, with the right-hand y-axis labels showing the amplitude
of the cross section in $cm^{2}\ \rm{g^{-1}}$. The effect of the resonant
cross section can be observed clearly in the $10^{8}$ and $10^{9}{\,\rm
M_{\odot}}$ halos with strong peaks at the resonant relative velocity in the
particle collision distribution. In the right figure, the grey vertical dashed
line shows the upper limit of the cross section model. The vertical dash
dotted lines indicate the collision rate effective cross section of the
respective halo masses. The delta function-like features are due to the cross
section of the plateau region at low relative velocity. From both figures,
with the particle collision distribution and the effective cross sections, it
is clear that the halo mass that maximizes the number of collisions happening
around the resonant peak is around $10^{8}{\,\rm M_{\odot}}$. Figure 2: The
left y-axis with the yellow and orange solid lines show the average and
conductivity effective cross section (calculated from Equation 17 and 18,
respectively) as functions of halo mass. The right y-axis and the two dashed
lines show the portion of each halo evolving under a specific cross section
range as a function of halo mass. The red line accounts for cross section
within the FWHM of the resonant peak compared to the cross section of the
plateau region at low relative velocity, while the blue line accounts for
cross section near the resonant peak ($\sigma_{V}/m>90\,{\rm cm}^{2}\,{\rm
g}^{-1}$). These profiles peak at the halo of mass $10^{7.9}{\,\rm M_{\odot}}$
and $10^{7.8}{\,\rm M_{\odot}}$, respectively, although it must also be noted
that the values for these two halo masses are very similar.

In the current picture of SIDM halo physics, the change in a halo density
profile over time follows a self-similar or universal trajectory when the
central density is scaled by $\rho_{\rm s}$, or a characteristic central
density, and time is scaled by a characteristic timescale, $t_{\rm col}$,
often referred to as the collapse timescale [22, 20, 50]. This timescale is
given by [51]

$t_{\rm col}(C,\sigma_{\rm{eff}}/m)=\frac{150}{C}\frac{1}{\left(\sigma_{\rm
eff}/m\right)\rho_{\rm s}r_{\rm s}}\frac{1}{\sqrt{4\pi G\rho_{\rm s}}}.$ (6)

By convention, the numerical prefactor is expressed as $150/C$ with a value
set to match predictions from SIDM fluid models [e.g. 30, 46, 50] with the
results from N-body simulations [e.g. 18, 30, 32, 20]. In the current picture
of SIDM halo evolution, $C$ is an $\mathcal{O}\left(1\right)$ constant, and
all SIDM and thermodynamic considerations enter in the derivation of the
effective cross section $\sigma_{\rm{eff}}$.

To test the hypothesis that SIDM halo evolution is universal when expressed in
terms of $t_{\rm{col}}$, we perform N-body simulations of isolated halos with
different cross section models, and compare their evolution using different
definitions of $\sigma_{\rm{eff}}$. We evolve one halo with the resonant cross
section shown in the left-hand panel of Figure 1. In addition, we simulate two
other halos, the first evolving with a velocity-independent cross section with
an amplitude $\sigma_{\rm{eff}}$ calculated by considering the scattering rate
in a halo (see Section III.1), and the second halo with $\sigma_{\rm{eff}}$
calculated based on the heat conduction cross section $\langle\sigma
v^{5}\rangle$ (see Section III.2). If the universality of halo evolution holds
for resonant models, the relation $\rho/\rho_{s}$ versus
$t/t_{\rm{col}}\left(C,\sigma_{\rm{eff}}\right)$ should appear almost
identical among the various simulations, modulo modest numerical systematic
effects that affect the precise onset of core collapse [e.g. 52].

###  III.1 The collision rate effective cross section

In contrast to a velocity-dependent cross section, the resonant cross section
shown on the left-hand side of Figure 1, exhibits a pronounced peak that
increases the cross section amplitude by a factor of 10 at $\sim 16\rm{km}\
\rm{s^{-1}}$. Therefore, we expect the most interesting structure formation
consequences in halos that have the portion of scattering occurring near the
peak of the resonance reach the maximum. To find the halo mass that maximizes
scattering near the peak, we derive the particle collision distribution
$p_{\rm collision}^{\rm core}\left(v_{\rm{rel}}\right)$, which quantifies the
fraction of scattering events in a halo occurring in a range of relative
velocities $\left[v_{\rm rel},v_{\rm rel}+{\rm d}v_{\rm rel}\right]$. From
this distribution, we will also calculate $p_{\rm collision}^{\rm
core}\left(\sigma_{V}/m\right)$, which represents the probability that a given
particle in a halo sees a cross section amplitude $\sigma_{V}$.

We assume an isotropic NFW DM density profile [12]

$\rho\left(r\right)=\frac{\rho_{\rm s}}{\left(r/r_{\rm
s}\right)\left(1+\left(r/r_{\rm s}\right)\right)^{2}},$ (7)

where $\rho_{\rm s}$ is the scale density, and $r_{\rm s}$ is the scale
radius. These parameters are calculated from the virial radius $r_{200}$,
defined as the radius at which the average density of the halo is 200 times
the critical density of the universe $\rho_{\rm crit}\,=\,3H^{2}/8\pi G$, and
the concentration parameter $c$, with $c=r_{200}/r_{\rm s}$. Both of these
parameters can be obtained with the virial mass $M_{200}$, the mass enclosed
within $r_{200}$ \- the former straightforwardly from its definition and the
latter from the low-mass concentration-mass (c-M) relation at redshift $z=0$
following [53]. Thus, we only need the virial mass $M_{200}$ to parameterize
our halos.

For dynamical equilibrium and the velocity distribution, we use Eddington’s
inversion formula [54]

$\rho\left(r\right)=\int_{0}^{v_{\rm esp}}4\pi v^{2}{\rm
d}v\,f\left(\xi\right),$ (8)

where $v_{\rm esp}$ is the escape velocity at radius $r$,
$\xi=\psi\left(r\right)-v^{2}/2$ is the binding energy per unit of mass,
$\psi$ is the negative gravitational potential, and $f\left(\xi\right)$ is a
phase-space distribution function. Here, we assume that the system is bounded
and isotropic, thus $f$ only depends on $\xi$, and $f\left(<0\right)=0$. The
Eddington distribution function satisfying Equation 8 takes the form

$f\left(\xi\right)=\frac{1}{\sqrt{8}\pi^{2}}\int_{0}^{\xi}{\rm
d}\psi\left[\frac{1}{\sqrt{\xi-\psi}}\frac{{\rm d}^{2}\rho}{{\rm
d}\psi^{2}}+\frac{1}{\sqrt{\xi}}\frac{{\rm d}\rho}{{\rm
d}\psi}\Bigr{\rvert}_{\psi=0}\right].$ (9)

This can be evaluated from the analytical forms of $\rho\left(r\right)$ and
$\psi\left(r\right)$ with $r$ as an intermediate variable. Integrating, we
obtain the semi-analytical form of the Eddington distribution function.

From the Eddington distribution function (Equation 9), we get the particle
velocity distribution function at radius $r$ following

$p_{v}^{\rm shell}\left(v\mid r\right)=4\pi
v^{2}\frac{f\left(\psi\left(r\right)-\frac{v^{2}}{2}\right)}{\rho\left(r\right)}.$
(10)

The particle relative velocity distribution function of particles within a
thin shell at radius $r$ and width ${\rm d}r$ can then be retrieved using the
equation

$p_{v_{\rm rel}}^{\rm shell}\left(v_{\rm rel}\mid r\right)=\int p_{v}^{\rm
shell}\left(v_{1}\mid r\right){\rm d}v_{1}\,p_{v}^{\rm shell}\left(v_{2}\mid
r\right){\rm d}v_{2}\,p_{\theta}\left(\theta\right)\,{\rm
d}\theta\,\delta\left(v_{\rm
rel}-\sqrt{v_{1}^{2}+v_{2}^{2}-2v_{1}v_{2}\cos{\theta}}\right),$ (11)

with $v_{1}$, $v_{2}$ as the velocities of the particle pair, and $\theta$ as
the relative angle between them. We assume the particle direction of movement
as isotropic, which results in the distribution function of $\theta$ as
follows

$p_{\theta}\left(\theta\right)=\sin{\theta}/2.$ (12)

To simplify calculations, we make use of the distribution of $\cos{\theta}$
rather than $\theta$. This distribution function takes the form of

$p_{\cos{\theta}}\left(\cos{\theta}\right)=1/2$ (13)

and is a uniform distribution with the range of $\left[-1,1\right]$.

Simply integrating a complex triple integral such as Equation 11 is
unnecessarily complicated analytically and undesirably costly numerically.
Therefore, we make use of a numerical trick involving the probability point
functions (PPF) - the inverse function of the cumulative distribution function
(CDF) - of $p_{v}^{\rm shell}\left(v\mid r\right)$ and
$p_{\cos{\theta}}\left(\cos{\theta}\right)$. We create a 1500x1500x150 grid,
and within each cell of the grid, we assign a
$\left(v_{1},v_{2},\cos{\theta}\right)$ set of values which are decided by
applying the appropriated PPF to a linearly-spaced array of the range
$\left[0,1\right]$. What we do here is essentially similar to combining the
three distribution functions, $p_{v}^{\rm shell}\left(v_{1}\mid r\right)$,
$p_{v}^{\rm shell}\left(v_{2}\mid r\right)$, and
$p_{\cos{\theta}}\left(\cos{\theta}\right)$, into a simultaneous distribution
of all three values and discretize it. We then calculate the relative
velocity, $v_{\rm rel}$, in each cell and bin the grid’s data into a
histogram. This functions similarly to integrating the distribution with the
delta function, and by normalizing the histogram, we acquire the distribution
function of the relative velocity inside a thin shell situated at radius $r$.

From the relative velocity distribution (Equation 11) we obtain the particle
collisions distribution as a function of relative velocity $v_{\rm rel}$ using
the following formula:

$p_{\rm collision}^{\rm shell}\left(v_{\rm rel}\mid r\right)=Cp_{v_{\rm
rel}}^{\rm shell}\left(v_{\rm rel}\mid r\right)R\left(v_{\rm rel},r\right),$
(14)

where C is the normalizing factor; $R\left(v_{\rm rel},r\right)$ is the
particle rate of collision, which follows

$R\left(v_{\rm rel},r\right)=v_{\rm rel}\,\sigma_{V}\left(v_{\rm
rel}\right)\rho\left(r\right),$ (15)

where $\sigma_{V}\left(v_{\rm rel}\right)$ is the cross section at the
relative velocity $v_{\rm rel}$.

The particle collision distribution for the halo core can then be obtained by
taking the weighted average of the thin shells’ particle collision
distributions with the shell’s mass

$p_{\rm collision}^{\rm core}\left(v_{\rm
rel}\right)=C^{\prime}\int_{0}^{3r_{\rm s}}4\pi r^{2}{\rm
d}r\rho\left(r\right)p_{\rm collision}^{\rm shell}\left(v_{\rm rel}\mid
r\right),$ (16)

where $C^{\prime}$ is another normalizing factor. Here we take the halo core
to be the region within 3 $r_{\rm s}$. Outside of this region, the number
density is low enough that collision rarely happens. We also compute the
collision distribution as a function of cross section $p_{\rm collision}^{\rm
core}\left(\sigma_{V}\right)$, which one can interpret as the probability that
a particle inside the halo sees a cross section amplitude $\sigma_{V}$. We
compute this distribution by sampling from $p_{\rm collision}^{\rm
core}\left(v_{\rm rel}\right)$, and evaluating $\sigma_{V}$ at the resulting
relative velocities. The particle collision distribution as a function of
relative velocity, $v_{\rm rel}$, and as a function of cross section per mass,
$\sigma_{V}/m$, for halos of mass $10^{7}$, $10^{8}$, $10^{9}$, and
$10^{10}{\,\rm M_{\odot}}$ are shown in Figure 1 with our resonant cross
section model.

In Figure 1, we notice the distinguished peaks in the collision distribution
at the resonant relative velocity for the higher mass halos and a fainter sign
for the lower mass halo in the left figure. In the right figure, we observe a
delta function-like feature at the cross section of the plateau region at low
relative velocity for all halo masses, though the feature weakens as the mass
increases and almost disappears as the halo mass reaches $10^{10}{\,\rm
M_{\odot}}$. At higher cross sections, the probability density increases with
the halo mass reaching its maximum around the $10^{8}{\,\rm M_{\odot}}$ halo.
The probability density then decreases as the mass increases further. All of
these features are expected from the features of our resonant cross section
model.

From the core’s particle collision distribution (Equation 16), we acquire the
average cross section seen by the halo following the formula

$\sigma_{\rm avg}=\int{\rm d}v_{\rm rel}\,p_{\rm collision}^{\rm
core}\left(v_{\rm rel}\right)\sigma_{V}\left(v_{\rm rel}\right).$ (17)

Onward we called this cross section the collision rate effective cross
section.

###  III.2 The heat conduction cross section

Additionally, existing treatments of gravothermal evolution with certain
classes of velocity-dependent cross section evolve in a self-similar way when
computed with a heat conductivity motivated effective cross section
$\sigma_{\kappa}$, assuming a Maxwell-Boltzmann distribution. We also consider
this approach and inspect the effective cross section [20, 46]

$\sigma_{\kappa}=\frac{\langle\sigma_{V}v^{5}\rangle}{\langle
v^{5}\rangle}=\frac{\int\sigma_{V}\left(v\right)v^{5}p_{v}^{\rm
MB}\left(v\right){\rm d}v}{\int v^{5}p_{v}^{\rm MB}\left(v\right){\rm d}v}.$
(18)

Here, $p_{v}^{\rm MB}\left(v\right)$ is the Maxwell-Boltzmann distribution

$p_{v}^{\rm MB}\left(v\right)\propto
v^{2}\exp\left(-\frac{v^{2}}{4v_{0}^{2}}\right),$ (19)

where $v_{0}=0.64\,v_{\rm max}$ with $v_{\rm max}=1.65\sqrt{G\rho_{\rm
s}r_{\rm s}^{2}}$ as the maximum circular velocity. The average and
conductivity effective cross sections of halos with masses in the range of
$10^{7}$ to $10^{10}{\,\rm M_{\odot}}$ are displayed in Figure 2. In the
figure are also two lines indicating the portion of each halo that will evolve
under cross sections within the full width at half maximum (FWHM) of the
resonant peak (compared to the cross section of the plateau region at low
relative velocity) $P_{\sigma_{V}/m>\sigma_{\rm FWHM}/m}$, and the portion of
each halo that will evolve under cross sections near the resonant peak
$P_{\sigma_{V}/m>90\,{\rm cm}^{2}\,{\rm g}^{-1}}$. We observe that halo with
mass of $10^{7.9}{\,\rm M_{\odot}}$ have the largest effective cross section
and the largest $P_{\sigma_{V}/m>\sigma_{\rm FWHM}/m}$, while halo with mass
of $10^{7.8}{\,\rm M_{\odot}}$ have the largest $P_{\sigma_{V}/m>90\,{\rm
cm}^{2}\,{\rm g}^{-1}}$. Onward, we will keep our focus on the evolution of
the $10^{7.9}{\,\rm M_{\odot}}$ halo, which has an effective cross sections of
$\sigma_{\rm avg}/m=53.91\,{\rm cm}^{2}\,{\rm g}^{-1}$ and
$\sigma_{\kappa}/m=31.98\,{\rm cm}^{2}\,{\rm g}^{-1}$. Other parameters of the
halo are shown in Table 1.

##  IV Simulations

Figure 3: Density (blue, left) and velocity dispersion (orange, right) profiles of the $10^{7.9}{\,\rm M_{\odot}}$ CDM halo. The black dashed lines show the halo analytical profiles following Equation 7, 20, and 22. The colored solid lines and the colored dotted lines show the profiles at 0 Gyr and 15 Gyr, respectively. The formers are generated from the initial conditions sampling code and used as the initial conditions of the CDM simulation. The laters show signs of small fluctuations in the innermost region of the halo due to the statistical error arising from the limited number of particles within the region, as well as fluctuations in the outer edge of the halo during the halo’s evolution to stability. There also exists a small core due to the effect of the softening length. However, the deviation from the analytical profiles resulting from this small core is of the same order of magnitude as the statistical fluctuation. The evolved 15 Gyr CDM halo is then used as the initial conditions for the SIDM runs. $\log M_{200}$ | $N_{200}$ | $m_{\rm DM}$ | $r_{200}$ | $c$ | $r_{\rm s}$ | $\log\rho_{\rm s}$ | $\epsilon$ | $\sigma\,m_{\rm\chi}^{-1}$ | $t_{\rm col}$
---|---|---|---|---|---|---|---|---|---
[${\rm M}_{\odot}$] |  | [${\rm M}_{\odot}$] | [$\rm kpc$] |  | [$\rm kpc$] | [${\rm M}_{\odot}\,{\rm kpc}^{-3}$ ] | [$\rm pc$] | [${\rm cm}^{2}\,{\rm g}^{-1}$] | [Gyr]
7.9 | $3\times 10^{7}$ | 2.64 | 9.08 | 18.7 | 0.49 | 7.43 | 1.00 | Resonant | Vary
7.9 | $3\times 10^{7}$ | 2.64 | 9.08 | 18.7 | 0.49 | 7.43 | 1.00 | $\sigma_{\rm{avg}}=53.91$ | 22.94
7.9 | $3\times 10^{7}$ | 2.64 | 9.08 | 18.7 | 0.49 | 7.43 | 1.00 | $\sigma_{\kappa}=31.98$ | 38.67
Table 1: Simulations configuration. (1) $M_{200}$ is the virial mass of the
halo. (2) $N_{200}$ is the number of DM particle within the virial radius (4)
$r_{200}$. (3) $m_{\rm DM}$ is the mass of a DM particle. (5) $c$ is the halo
concentration parameter. (6) $\rho_{\rm s}$ and (7) $r_{\rm s}$ are the scale
density and radius of the NFW profile. (8) $\epsilon$ is the gravitational
softening length. (9) $\sigma\,m_{\rm\chi}^{-1}$ is the cross section per unit
of mass of the SIDM collisions. (10) $t_{\rm col}$ is the collapse time scale,
defined in Equation 6.

We perform dark-matter-only simulations of isolated halos using the multi-
physics, massively parallel, moving-mesh magneto-hydrodynamic (MHD) simulation
code Arepo, which is publicly available [55]. Gravity is solved using the
Tree-Particle Mesh (Tree-PM) method. Numeric details of the simulations are
summarized in Table 1. Below we will briefly describe the construction of
initial conditions and the modeling of DM self-interactions.

The isolated halos are initialized with $3\times 10^{7}$ DM particles within
$r_{200}$ with distribution of particles extending to $3\times r_{200}$. The
density profile within $r_{200}$ follows the NFW profile [12]. Beyond
$r_{200}$, we use an exponential cut-off following [56]

$\rho\left(r\right)=\frac{\rho_{\rm
s}}{c\left(1+c\right)^{2}}\left(\frac{r}{r_{200}}\right)^{\epsilon_{\rm
decay}}\exp\left(-\frac{r-r_{200}}{r_{\rm decay}}\right),$ (20)

where $r_{\rm decay}$ is the decay scale chosen to be $r_{200}$ and
$\epsilon_{\rm decay}$ is chosen so that the continuity of the logarithmic
slope of the density profile is preserved

$\epsilon_{\rm decay}=\frac{r_{200}}{r_{\rm decay}}-\frac{1+3c}{1+c}.$ (21)

We assume the dynamical equilibrium of the halo and again use the Eddington
distribution function (Equation 9) to calculate the phase-space distribution
function. To explicitly sample the positions and velocities of particles, we
use a combination of analytical formulas and numerical calculations with cubic
spline interpolations. The integrated mass $M\left(r\right)$ and negative
gravitational potential $\psi\left(r\right)$ are calculated analytically for
$r\leq r_{200}$ and numerically for $r>r_{200}$. Following the procedure in
[57], we obtain the density profile as a function of negative gravitational
potential, $\rho\left(\psi\right)$, from the interpolation of tabulated values
of $\rho\left(r\right)$, $M\left(r\right)$, and $\psi\left(r\right)$ from
$10^{-5}\,r_{200}$ to $10^{3}\,r_{200}$. We validate our calculation by
generating the halo density profile using the acquired Eddington’s
distribution function. The generated density profile converges almost
perfectly with the analytical formulas in equation 7 and 20. The relative
error is typically orders of magnitude below 10-2 for $r>10^{-4}\,r_{200}$ but
is significantly higher for $r<10^{-4}\,r_{200}$ which can be ignored as the
region is sufficiently small.

We sample the positions of particles with the inversion sampling method using
the density profile as the distribution function. From the isotropic
assumption, we then sample a random direction for each particle and calculate
the Cartesian coordinates. For the velocities, we sample a random direction
for each particle. In each radial bin, we sample the velocity magnitudes of
particles via rejection sampling with the Eddington distribution function at
the bin centers. We adopt stretched Gaussian proposal functions, which are
fitted to minimize the rejection rate in each bin.

The initial conditions of a $10^{7.9}{\,\rm M_{\odot}}$ halo derived through
this approach is shown in Figure 3. The stability of the ICs is tested by
running a CDM simulation for $15\,{\rm Gyr}$ and inspecting the evolution of
the density and velocity dispersion profiles. We compare the profiles with
their analytical formulas. The analytical density profile is described by
Equation 7 and Equation 20. We obtain the velocity dispersion profile by
solving the Jeans equation [54]

$\frac{1}{\rho\left(r\right)}\frac{\rm d}{{\rm
d}r}\left(\rho\left(r\right)\sigma_{\rm r}^{2}\left(r\right)\right)=\frac{{\rm
d}\psi\left(r\right)}{{\rm d}r}.$ (22)

Figure 3 shows the stability of the $10^{7.9}{\,\rm M_{\odot}}$ halo in CDM.
We find that sampling particles to 3 $r_{200}$ provide sufficient supporting
structures for the halo. There are small fluctuations in the density and
velocity dispersion profiles at the outer edge of the halo, which takes a few
Gyr to evolve into its equilibrium state. Meanwhile, the central density and
velocity dispersion also fluctuate, and eventually a small core forms due to
two-body relaxation. Nevertheless, the deviation from the analytical profiles
created by this core is of the same order of magnitude as the fluctuation due
to statistical noises from the limited number of particles inside the region.
After 15 Gyr, our halo becomes totally stable (as shown in Figure 3). Its
density and velocity dispersion profiles follow closely the analytical
profiles except for the small statistical fluctuation.

We also compared the particle relative velocity distribution of the sampled
halo with the analytical result from Section III. Using the k-d tree
algorithm, we find all particle pairs with distances smaller than
$0.01\,r_{\rm s}$ in the central region of the halo ($r_{\rm core}=r_{\rm s}$
or $3\,r_{\rm s}$). This distance is chosen empirically to ensure a
statistically significant number of particle pairs while keeping the mixing of
particles in different shells at a minimum. The relative velocities of the
pairs are then calculated and used to evaluate the relative velocity
distribution $p_{v_{\rm rel}}^{\rm core}\left(v_{\rm rel}\right)$. This
distribution is equivalent to

$p_{v_{\rm rel}}^{\rm core}\left(v_{\rm rel}\right)=D\int_{0}^{r_{\rm
core}}4\pi r^{2}{\rm d}r\rho\left(r\right)^{2}p_{v_{\rm rel}}^{\rm
shell}\left(v_{\rm rel}\mid r\right),$ (23)

where $D$ is a normalizing factor. This relative velocity distribution takes
into account the shell’s densities which influence the number of particle
pairs within each shell. We find that throughout the CDM halo evolution, its
relative velocity distribution for both choices of core radius ($r_{\rm s}$
and $3\,r_{\rm s}$) follows strictly (with only some minor fluctuations and
deviations) the analytical distributions.

The CDM halo evolved for 15 Gyr, at which point it settles into dynamic
equilibrium, is used as the initial condition for all of our SIDM simulations.
DM self-interactions are simulated in a Monte-Carlo fashion as implemented in
[10, 58]. Scatterings of DM particles are assumed to be elastic and isotropic.
The implementation ensures explicit conservation of energy and momentum. The
time steps of simulations are regulated by the estimated SIDM scattering rates
seen by each particle and are short enough to avoid multiple scatterings per
time step. We refer readers to [10] for detailed descriptions of the numeric
implementation.

We simulate one SIDM halo with the resonant velocity-dependent cross section
and two SIDM halos with velocity-independent cross sections. The velocity-
independent cross sections are set at the values of $\sigma_{\rm avg}/m$ and
$\sigma_{\kappa}/m$, calculated from Equation 17 and Equation 18,
respectively. The exact values for these effective cross section amplitudes
are shown in Table 1.

##  V Results

Figure 4: Time evolution of the particle relative velocity distribution inside
the halo core $p_{v_{\rm rel}}^{\rm core}\left(v_{\rm rel}\right)$ (Equation
23) measured directly from snapshots of the resonant SIDM simulation. Here,
distributions from two choices of halo core $r_{\rm core}=r_{\rm s}$ and
$r_{\rm core}=3\,r_{\rm s}$ are shown in solid and dashed lines, respectively.
The analytical distributions for CDM halo with NFW profile are also shown in
dark blue. As the SIDM halo goes into the core-formation phase, the relative
velocity distribution approaches the Maxwell-Boltzmann distribution (shown by
squares) calculated from Equation 19 with the central radial velocity
dispersion $\sigma^{2}_{r,{\rm c}}$ as $v_{0}^{2}$. The resonant model cross
section profile is shown in yellow and scaled logarithmically in both axices.

Figure 5: Time evolution of DM central (at $r\leq 0.1\,r_{\rm s}$) densities
in the VICS and the RCS SIDM models. The times are normalized by the collapse
time $t_{\rm col}$ (as defined in Equation 6) of the corresponding model. For
the halo in the RCS model, we also normalize the collapse time with an
adaptive timescale $\tilde{t}_{\rm col}$ described in Equation 24 that take
into account of the change in the relative velocity distribution and,
consequently, the conductivity cross section $\sigma_{\kappa}/m$ over time.
The top panel shows the results when the same prefactor of conduction
coefficient $C=0.75$ is assumed for both VICS and RCS models. The bottom panel
shows the results when the prefactor for the RCS model is changed to $C_{\rm
res}\sim 0.6$. In general, halos in SIDM undergo a first phase of core
formation until thermal equilibrium is reached by DM self-interactions. The
collapse time based on $\sigma_{\rm avg}/m$ is better in describing the core
formation phase in the RCS model. On the contrary, the later gravothermal
collapse phase is better described by $\sigma_{\kappa}/m$ despite about $20\%$
delay in collapse time, which we attempt to capture by adjusting $C_{\rm res}$
as indicated. This adjustment also better matches the normalized central
density evolution described by $\sigma_{\rm avg}/m$ to the core formation
phase.

Figure 6: Time evolution of DM density profiles (top), velocity dispersion
profiles (middle), and logarithmic profile slopes (bottom) of the simulated
halo in the VICS and RCS SIDM models. The mean relative velocity $\langle
v_{\rm rel}\rangle\equiv\sigma_{\rm r}\times 4/\sqrt{\pi}$, assuming a
Maxwell-Boltzmann distribution, is shown in the bottom panel as a reference.
The qualitative evolution of halo structure in the VICS and RCS models is
similar. The rapid core-formation phase ($t\lesssim 0.2\,t_{\rm col}$) is
followed by the gravothermal collapse phase, where the halo centers maintain
isothermal and the density profiles are self-similar. The velocity dispersion
profile of the halo increases mildly through the collapse phase.
Quantitatively, such a small ($\lesssim 1\,{\rm km}\,{\rm s}^{-1}$) increase
in $\langle v_{\rm rel}\rangle$ does not change the effective cross-section
and conduction efficiency substantially. Therefore, a fixed $t_{\rm col}$ or
equivalently $C_{\rm res}$ can approximate the entire collapse phase.

In this section, we inspect the evolution of the halo in the velocity-
independent cross section (VICS) and resonant cross section (RCS) SIDM models.

Figure 4 shows the relative velocity distributions $p_{v_{\rm rel}}^{\rm
core}\left(v_{\rm rel}\right)$ of the halo in the RCS model at 0, 0.1, and 1.1
$\tilde{t}_{\rm col}$ with $\tilde{t}_{\rm col}$ as the adaptive time scale
defined below in Equation 24. As mentioned in Section IV, this distribution is
read directly from the simulation snapshots and takes into account the
weighting factor of the local particle density. We observe that the relative
velocity distribution rapidly transforms into Maxwell-Boltzmann distributions
(Equation 19) with the central radial velocity dispersion $\sigma_{r,{\rm c}}$
as as the characteristic velocity $v_{0}$ in the core formation phase and
remain so during the core collapse. These Maxwell-Boltzmann distributions are
shown by squares in Figure 4.

In Figure 5, we show the evolution of halo central density versus time in the
RCS and two VICS models. The time is normalized by the collapse timescale
$t_{\rm{col}}$, as defined in Equation 6. The central density of each halo is
calculated as the average DM density within $r_{\rm c}=0.1\,r_{\rm s}$ 111The
value of $r_{\rm c}$ is chosen to reduce statistical uncertainties while
keeping the region small enough to represent the halo central density.. The
top panel shows the evolution of halos with collapse timescales calculated
using two different definitions of the effective cross section ($\sigma_{\rm
avg}/m$ and $\sigma_{\kappa}/m$) and $C=0.75$. The green and pink curves
represent halos evolving with VICS of amplitudes $\sigma_{\rm{avg}}/m$ and
$\sigma_{\kappa}/m$, respectively. The orange (blue) curve shows the evolution
of the halo in the RCS model when using $\sigma_{\rm{avg}}/m$
($\sigma_{\rm{\kappa}}/m$) to compute $t_{\rm{col}}$.

In addition, we define an adaptive time scale $\tilde{t}_{\rm{col}}$ to
account for the change of $\sigma_{\rm{\kappa}}/m$ over time due to the shift
of $p_{v_{\rm rel}}^{\rm core}\left(v_{\rm rel}\right)$ as seen in Figure 4.
The phase of collapse, $t/\tilde{t}_{\rm{col}}$, is defined following

${\rm d}\left(\frac{t}{\tilde{t}_{\rm{col}}}\right)=\frac{{\rm d}t}{t_{\rm
col}\left(\sigma_{\kappa}\left(t\right)/m\right)}.$ (24)

where $t$ is the physical time, and $\sigma_{\kappa}\left(t\right)/m$ is
calculated from Equation 18 using $p_{v}^{\rm MB}\left(v\right)=p_{v_{\rm
rel}}^{\rm core}\left(v_{\rm rel}\right)$ as measured according to the
procedural detailed in Section IV. The insert in the bottom panel of Figure 5
shows the evolution of the conductivity effective cross section
$\sigma_{\kappa}\left(t\right)/m$ over physical time. Generally,
$\sigma_{\kappa}\left(t\right)/m$ evolve proportional to $\sim t^{-2}$,
decreasing from $\simeq 33.5\,{\rm cm}^{2}\,{\rm g}^{-1}$ at $t=0$ Gyr to
$\simeq 26\,{\rm cm}^{2}\,{\rm g}^{-1}$ at the time the central density reach
its initial value $t\sim 55$ Gyr. The dashed blue line in Figure 5 shows the
evolution of the central density in the halo in the RCS model rescaled with
$\tilde{t}_{\rm col}$.

As expected, the halos evolving with velocity-independent cross sections
exhibit nearly universal evolution when time is scaled by $t_{\rm{col}}$, with
a small deviation in the core-collapse phase likely due to numerical effects
[e.g. 52, 60, 61]. For the resonant model, we find that computing
$t_{\rm{col}}$ with $\sigma_{\rm avg}/m$ (orange curve) provides a better
description of the initial core-formation phase until $t\sim 0.2\
t_{\rm{col}}$, but the agreement breaks down for $t\gtrsim 0.2\ t_{\rm{col}}$.
In the later stages of evolution, after the core-formation phase, calculating
$t_{\rm{col}}$ using $\sigma_{\kappa}/m$ (blue curve) more closely tracks the
evolution of the velocity-independent models. However, the halo evolving with
a resonant cross section still takes $\sim 20\%$ longer in rescaled time units
$t/t_{\rm{col}}$ to reach its initial central density during the collapse
phase than the halos evolving with VICS. This is also the case for the
evolution in the RCS model scaled by the adaptive time unit
$t/\tilde{t}_{\rm{col}}$.

In the bottom panel of Figure 5 we show the blue and orange curves for the
resonant model after calculating $t_{\rm{col}}$ with $C_{\rm res}=0.62$, which
one could also interpret a modification of $\sigma_{\rm{eff}}$ at the level of
$\sim 20\%$ to match the evolution in the core collapse phase. Using $C_{\rm
res}=0.62$ with the adaptive time scale $\tilde{t}_{\rm{col}}$ as shown by the
dashed blue line in the bottom panel of Figure 5 maps the central density
evolution in the RCS model exactly along the rescaled collapse track of the
halo with VICS of amplitude $\sigma_{\kappa}/m$ for $t\gtrsim
0.4\,t_{\rm{col}}$. In addition, the choice of $C_{\rm res}=0.62$ also better
matches the central density evolution of the halo in the RCS model when
rescaled with $t_{\rm{col}}\left(\sigma_{\rm avg}/m\right)$ to the universal
track in the core formation phase.

However, even after introducing this ad-hoc rescaling of the collapse
timescale, as well as utilizing different time scales to describe different
phases, the behavior during the late stage of core formation and early stage
of core collapse ($0.2\,t_{\rm col}\lesssim t\lesssim 0.4\,t_{\rm col}$) still
differ from those of the halos we evolve with velocity-independent cross
sections. The minimum core density of the halo in the RCS model is smaller
than that in the VICS model by a factor of $\simeq 10\%$. We interpret these
deviations as evidence of a violation of universality at the level of $\sim
10-20\%$ for halos evolving with resonant cross sections.

The evolution of the density and velocity dispersion profiles are shown in
Figure 6. We show the DM density and relative velocities at several snapshots
in one velocity-independent model (the halo with a velocity-independent cross
section $\sigma_{\kappa}/m$) and the resonant model. Based on Figure 5, we use
$\tilde{t}_{\rm col}$ and $C_{\rm res}=0.62$ to define the collapse timescale
for the resonant model to match the evolution during the collapse phase seen
in the velocity-independent models.

##  VI Discussion and conclusions

We have performed idealized N-body simulations of DM halos to examine halo
evolution in a SIDM model with a resonant cross section. We compare the
evolution of the halo with a resonant cross section with the evolution of
halos evolving with velocity-independent cross sections, and examine the
change in their internal structure over time using two definitions of an
effective cross section and core collapse timescale. Our main findings are
summarized as follows:

  * •

The resonant model exhibits deviation from universality in the evolution of
its density profile from halos evolving with velocity-independent cross
sections. The halo evolving with a resonant cross section achieves a minimum
core density $10\%$ lower than the halos with velocity-independent scattering
cross sections, and reaches its initial density $20\%$ later when time is
scaled by the thermal conduction cross section
$\sigma_{\kappa}\sim\langle\sigma v^{5}\rangle/\langle v^{5}\rangle$.

  * •

In the early thermalization and core-formation phase in the resonant model,
the halo evolution is better described by the averaged cross-section
($\sigma_{\rm avg}/m$) “seen” by particles having self-interactions. After
core formation, the late gravothermal collapse phase in the resonant model is
better described by $\sigma_{\kappa}/m$, although the time evolution resulting
from this definition of the effective cross section still deviates
significantly from the velocity-independent models.

  * •

For our specific resonant SIDM model and halo mass choice, using $C_{\rm
res}=0.62$ matches the central density evolution when rescaled with $t_{\rm
col}\left(C_{\rm res},\sigma_{\rm avg}\right)$ exactly with the self-similar
track in the core formation phase. Additionally, when rescaled with the
adaptive time scale $\tilde{t}_{\rm col}\left(C_{\rm
res},\sigma_{\kappa}\right)$, the evolution of the central density follows
precisely along that of the halo evolved with a velocity-independent cross
section of amplitude $\sigma_{\kappa}$ during most of the core collapse.

The utility of a universal evolution of SIDM halos lies in its predictive
power; given the physical properties of a CDM halo ($\rho_{s}$ and $r_{s}$)
and the form of the SIDM cross section, the universal treatment of halo
evolution based on the collapse timescale enables the prediction of minimum
core density and the structure of the halo during the collapse phase [20].
While we observe a deviation from universality, it is only at the level of
$10-20\%$. Other physical processes can impact SIDM halo density profiles to a
similar or even larger degree, particularly tidal stripping and tidal heating
effects relevant for galactic subhalos [e.g. 23, 26, 27].

A modified calculation of the core collapse timescale based on a time-
dependent definition of the thermal conduction cross section $\sigma_{\kappa}$
provides an improved description of the collapse phase for the halo evolving
with a resonant cross section. This suggests that a modified calculation of
the thermal conduction cross section could restore universality for resonant
cross sections, provided this modified treatment reverts to the definition of
$\sigma_{\kappa}$ for velocity-independent cross sections. Exploring this
topic further, and simulating the evolution of halos with different resonant
cross sections, will be the topic of future study.

###### Acknowledgements.

We thank Annika Peter, Antti Asikainen, and Yi-Ming Zhong for useful
discussions. DG acknowledges support for this work provided by the Brinson
Foundation through a Brinson Prize Fellowship grant. MV acknowledges support
through NASA ATP 19-ATP19-0019, 19-ATP19-0020, 19-ATP19-0167, and NSF grants
AST-1814053, AST-1814259, AST-1909831, AST-2007355 and AST-2107724.

## References

  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Particle dark matter: evidence, candidates and constraints, Phys. Rep. 405, 279 (2005), arXiv:hep-ph/0404175 [hep-ph] .
  * Buckley and Peter [2018] M. R. Buckley and A. H. G. Peter, Gravitational probes of dark matter physics, Phys. Rep. 761, 1 (2018), arXiv:1712.06615 [astro-ph.CO] .
  * Drlica-Wagner  _et al._ [2022] A. Drlica-Wagner, C. Prescod-Weinstein, H.-B. Yu, A. Albert, M. Amin, A. Banerjee, M. Baryakhtar, K. Bechtol, S. Bird, S. Birrer, T. Bringmann, R. Caputo, S. Chakrabarti, T. Y. Chen, D. Croon, F.-Y. Cyr-Racine, W. A. Dawson, C. Dvorkin, V. Gluscevic, D. Gilman, D. Grin, R. Hložek, R. K. Leane, T. S. Li, Y.-Y. Mao, J. Meyers, S. Mishra-Sharma, J. B. Muñoz, F. Munshi, E. O. Nadler, A. Parikh, K. Perez, A. H. G. Peter, S. Profumo, K. Schutz, N. Sehgal, J. D. Simon, K. Sinha, M. Valluri, and R. H. Wechsler, Report of the Topical Group on Cosmic Probes of Dark Matter for Snowmass 2021, arXiv e-prints , arXiv:2209.08215 (2022), arXiv:2209.08215 [hep-ph] .
  * Bechtol  _et al._ [2022] K. Bechtol, S. Birrer, F.-Y. Cyr-Racine, K. Schutz, S. Adhikari, M. Amin, A. Banerjee, S. Bird, N. Blinov, K. K. Boddy, C. Boehm, K. Bundy, M. Buschmann, S. Chakrabarti, D. Curtin, L. Dai, A. Drlica-Wagner, C. Dvorkin, A. L. Erickcek, D. Gilman, S. Heeba, S. Kim, V. Iršič, A. Leauthaud, M. Lovell, Z. Lukić, Y.-Y. Mao, S. Mau, A. Mitridate, P. Mocz, J. B. Muñoz, E. O. Nadler, A. H. G. Peter, A. Price-Whelan, A. Robertson, N. Sabti, N. Sehgal, N. Shipp, J. D. Simon, R. Singh, K. Van Tilburg, R. H. Wechsler, A. Widmark, and H.-B. Yu, Snowmass2021 Cosmic Frontier White Paper: Dark Matter Physics from Halo Measurements, arXiv e-prints , arXiv:2203.07354 (2022), arXiv:2203.07354 [hep-ph] .
  * Carlson  _et al._ [1992] E. D. Carlson, M. E. Machacek, and L. J. Hall, Self-interacting Dark Matter, ApJ 398, 43 (1992).
  * de Laix  _et al._ [1995] A. A. de Laix, R. J. Scherrer, and R. K. Schaefer, Constraints on Self-interacting Dark Matter, ApJ 452, 495 (1995), arXiv:astro-ph/9502087 [astro-ph] .
  * Firmani  _et al._ [2000] C. Firmani, E. D’Onghia, V. Avila-Reese, G. Chincarini, and X. Hernández, Evidence of self-interacting cold dark matter from galactic to galaxy cluster scales, MNRAS 315, L29 (2000), arXiv:astro-ph/0002376 [astro-ph] .
  * Spergel and Steinhardt [2000] D. N. Spergel and P. J. Steinhardt, Observational Evidence for Self-Interacting Cold Dark Matter, Phys. Rev. Lett. 84, 3760 (2000), arXiv:astro-ph/9909386 [astro-ph] .
  * Tulin  _et al._ [2013] S. Tulin, H.-B. Yu, and K. M. Zurek, Beyond collisionless dark matter: Particle physics dynamics for dark matter halo structure, Phys. Rev. D 87, 115007 (2013), arXiv:1302.3898 [hep-ph] .
  * Vogelsberger  _et al._ [2012] M. Vogelsberger, J. Zavala, and A. Loeb, Subhaloes in self-interacting galactic dark matter haloes, MNRAS 423, 3740 (2012), arXiv:1201.5892 [astro-ph.CO] .
  * Adhikari  _et al._ [2022] S. Adhikari, A. Banerjee, K. K. Boddy, F.-Y. Cyr-Racine, H. Desmond, C. Dvorkin, B. Jain, F. Kahlhoefer, M. Kaplinghat, A. Nierenberg, A. H. G. Peter, A. Robertson, J. Sakstein, and J. Zavala, Astrophysical Tests of Dark Matter Self-Interactions, arXiv e-prints , arXiv:2207.10638 (2022), arXiv:2207.10638 [astro-ph.CO] .
  * Navarro  _et al._ [1997] J. F. Navarro, C. S. Frenk, and S. D. M. White, A Universal Density Profile from Hierarchical Clustering, ApJ 490, 493 (1997), arXiv:astro-ph/9611107 [astro-ph] .
  * Rocha  _et al._ [2013] M. Rocha, A. H. G. Peter, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, J. Oñorbe, and L. A. Moustakas, Cosmological simulations with self-interacting dark matter - I. Constant-density cores and substructure, MNRAS 430, 81 (2013), arXiv:1208.3025 [astro-ph.CO] .
  * Zavala  _et al._ [2013] J. Zavala, M. Vogelsberger, and M. G. Walker, Constraining self-interacting dark matter with the Milky way’s dwarf spheroidals., MNRAS 431, L20 (2013), arXiv:1211.6426 [astro-ph.CO] .
  * Elbert  _et al._ [2015] O. D. Elbert, J. S. Bullock, S. Garrison-Kimmel, M. Rocha, J. Oñorbe, and A. H. G. Peter, Core formation in dwarf haloes with self-interacting dark matter: no fine-tuning necessary, MNRAS 453, 29 (2015), arXiv:1412.1477 [astro-ph.GA] .
  * Lynden-Bell  _et al._ [1968] D. Lynden-Bell, R. Wood, and A. Royal, The Gravo-Thermal Catastrophe in Isothermal Spheres and the Onset of Red-Giant Structure for Stellar Systems, Monthly Notices of the Royal Astronomical Society 138, 495 (1968).
  * Balberg  _et al._ [2002] S. Balberg, S. L. Shapiro, and S. Inagaki, Self-Interacting Dark Matter Halos and the Gravothermal Catastrophe, ApJ 568, 475 (2002), arXiv:astro-ph/0110561 [astro-ph] .
  * Koda and Shapiro [2011] J. Koda and P. R. Shapiro, Gravothermal collapse of isolated self-interacting dark matter haloes: N-body simulation versus the fluid model, MNRAS 415, 1125 (2011), arXiv:1101.3097 [astro-ph.CO] .
  * Yang and Yu [2021] D. Yang and H.-B. Yu, Self-interacting dark matter and small-scale gravitational lenses in galaxy clusters, Physical Review D 104, 103031 (2021).
  * Yang and Yu [2022] D. Yang and H.-B. Yu, Gravothermal evolution of dark matter halos with differential elastic scattering, Journal of Cosmology and Astroparticle Physics 2022 (09), 077.
  * Nadler  _et al._ [2023] E. O. Nadler, D. Yang, and H.-B. Yu, A Self-interacting Dark Matter Solution to the Extreme Diversity of Low-mass Halo Properties, ApJ 958, L39 (2023), arXiv:2306.01830 [astro-ph.GA] .
  * Outmezguine  _et al._ [2023] N. J. Outmezguine, K. K. Boddy, S. Gad-Nasr, M. Kaplinghat, and L. Sagunski, Universal gravothermal evolution of isolated self-interacting dark matter halos for velocity-dependent cross-sections, MNRAS 523, 4786 (2023), arXiv:2204.06568 [astro-ph.GA] .
  * Nishikawa  _et al._ [2020] H. Nishikawa, K. K. Boddy, and M. Kaplinghat, Accelerated core collapse in tidally stripped self-interacting dark matter halos, Phys. Rev. D 101, 063009 (2020), arXiv:1901.00499 [astro-ph.GA] .
  * Nadler  _et al._ [2020] E. O. Nadler, A. Banerjee, S. Adhikari, Y.-Y. Mao, and R. H. Wechsler, Signatures of Velocity-dependent Dark Matter Self-interactions in Milky Way-mass Halos, ApJ 896, 112 (2020), arXiv:2001.08754 [astro-ph.CO] .
  * Sameie  _et al._ [2020] O. Sameie, H.-B. Yu, L. V. Sales, M. Vogelsberger, and J. Zavala, Self-Interacting Dark Matter Subhalos in the Milky Way’s Tides, Phys. Rev. Lett. 124, 141102 (2020), arXiv:1904.07872 [astro-ph.GA] .
  * Zeng  _et al._ [2022] Z. C. Zeng, A. H. G. Peter, X. Du, A. Benson, S. Kim, F. Jiang, F.-Y. Cyr-Racine, and M. Vogelsberger, Core-collapse, evaporation, and tidal effects: the life story of a self-interacting dark matter subhalo, MNRAS 513, 4845 (2022), arXiv:2110.00259 [astro-ph.CO] .
  * Slone  _et al._ [2023] O. Slone, F. Jiang, M. Lisanti, and M. Kaplinghat, Orbital evolution of satellite galaxies in self-interacting dark matter models, Phys. Rev. D 107, 043014 (2023).
  * Carton Zeng  _et al._ [2023] Z. Carton Zeng, A. H. G. Peter, X. Du, S. Yang, A. Benson, F.-Y. Cyr-Racine, F. Jiang, C. Mace, and R. B. Metcalf, Till the core collapses: the evolution and properties of self-interacting dark matter subhalos, arXiv e-prints , arXiv:2310.09910 (2023), arXiv:2310.09910 [astro-ph.GA] .
  * Feng  _et al._ [2021] W.-X. Feng, H.-B. Yu, and Y.-M. Zhong, Seeding Supermassive Black Holes with Self-interacting Dark Matter: A Unified Scenario with Baryons, ApJ 914, L26 (2021), arXiv:2010.15132 [astro-ph.CO] .
  * Essig  _et al._ [2019] R. Essig, S. D. McDermott, H.-B. Yu, and Y.-M. Zhong, Constraining Dissipative Dark Matter Self-Interactions, Phys. Rev. Lett. 123, 121102 (2019), arXiv:1809.01144 [hep-ph] .
  * Vogelsberger  _et al._ [2019] M. Vogelsberger, J. Zavala, K. Schutz, and T. R. Slatyer, Evaporating the Milky Way halo and its satellites with inelastic self-interacting dark matter, MNRAS 484, 5437 (2019), arXiv:1805.03203 [astro-ph.GA] .
  * Shen  _et al._ [2021] X. Shen, P. F. Hopkins, L. Necib, F. Jiang, M. Boylan-Kolchin, and A. Wetzel, Dissipative dark matter on FIRE - I. Structural and kinematic properties of dwarf galaxies, MNRAS 506, 4421 (2021), arXiv:2102.09580 [astro-ph.GA] .
  * O’Neil  _et al._ [2023] S. O’Neil, M. Vogelsberger, S. Heeba, K. Schutz, J. C. Rose, P. Torrey, J. Borrow, R. Low, R. Adhikari, M. V. Medvedev, T. R. Slatyer, and J. Zavala, Endothermic self-interacting dark matter in Milky Way-like dark matter haloes, MNRAS 524, 288 (2023), arXiv:2210.16328 [astro-ph.GA] .
  * Kaplinghat  _et al._ [2016] M. Kaplinghat, S. Tulin, and H.-B. Yu, Dark Matter Halos as Particle Colliders: Unified Solution to Small-Scale Structure Puzzles from Dwarfs to Clusters, Phys. Rev. Lett. 116, 041302 (2016), arXiv:1508.03339 [astro-ph.CO] .
  * Tulin and Yu [2018] S. Tulin and H.-B. Yu, Dark matter self-interactions and small scale structure, Phys. Rep. 730, 1 (2018), arXiv:1705.02358 [hep-ph] .
  * Silverman  _et al._ [2023] M. Silverman, J. S. Bullock, M. Kaplinghat, V. H. Robles, and M. Valli, Motivations for a large self-interacting dark matter cross-section from Milky Way satellites, MNRAS 518, 2418 (2023), arXiv:2203.10104 [astro-ph.GA] .
  * Shah and Adhikari [2024] N. Shah and S. Adhikari, The abundance of core-collapsed subhaloes in SIDM: insights from structure formation in $\Lambda$CDM, MNRAS 529, 4611 (2024), arXiv:2308.16342 [astro-ph.CO] .
  * Gilman  _et al._ [2021] D. Gilman, J. Bovy, T. Treu, A. Nierenberg, S. Birrer, A. Benson, and O. Sameie, Strong lensing signatures of self-interacting dark matter in low-mass haloes, MNRAS 507, 2432 (2021), arXiv:2105.05259 [astro-ph.CO] .
  * Gilman  _et al._ [2023] D. Gilman, Y.-M. Zhong, and J. Bovy, Constraining resonant dark matter self-interactions with strong gravitational lenses, Phys. Rev. D 107, 103008 (2023), arXiv:2207.13111 [astro-ph.CO] .
  * Minor  _et al._ [2021] Q. Minor, S. Gad-Nasr, M. Kaplinghat, and S. Vegetti, An unexpected high concentration for the dark substructure in the gravitational lens SDSSJ0946+1006, MNRAS 507, 1662 (2021), arXiv:2011.10627 [astro-ph.GA] .
  * Peter  _et al._ [2013] A. H. G. Peter, M. Rocha, J. S. Bullock, and M. Kaplinghat, Cosmological simulations with self-interacting dark matter – II. Halo shapes versus observations, Monthly Notices of the Royal Astronomical Society 430, 105 (2013).
  * Harvey  _et al._ [2015] D. Harvey, R. Massey, T. Kitching, A. Taylor, and E. Tittley, The nongravitational interactions of dark matter in colliding galaxy clusters, Science 347, 1462 (2015), arXiv:1503.07675 [astro-ph.CO] .
  * Sagunski  _et al._ [2021] L. Sagunski, S. Gad-Nasr, B. Colquhoun, A. Robertson, and S. Tulin, Velocity-dependent self-interacting dark matter from groups and clusters of galaxies, Journal of Cosmology and Astroparticle Physics 2021 (01), 024.
  * Correa [2021] C. A. Correa, Constraining velocity-dependent self-interacting dark matter with the Milky Way’s dwarf spheroidal galaxies, MNRAS 503, 920 (2021), arXiv:2007.02958 [astro-ph.GA] .
  * Turner  _et al._ [2021] H. C. Turner, M. R. Lovell, J. Zavala, and M. Vogelsberger, The onset of gravothermal core collapse in velocity-dependent self-interacting dark matter subhaloes, MNRAS 505, 5327 (2021), arXiv:2010.02924 [astro-ph.GA] .
  * Yang  _et al._ [2023a] S. Yang, X. Du, Z. C. Zeng, A. Benson, F. Jiang, E. O. Nadler, and A. H. G. Peter, Gravothermal Solutions of SIDM Halos: Mapping from Constant to Velocity-dependent Cross Section, ApJ 946, 47 (2023a), arXiv:2205.02957 [astro-ph.CO] .
  * Chu  _et al._ [2020] X. Chu, C. Garcia-Cely, and H. Murayama, A practical and consistent parametrization of dark matter self-interactions, J. Cosmol. Astropart. Phys. 2020, 043 (2020), arXiv:1908.06067 [hep-ph] .
  * Kamada and Kim [2024] A. Kamada and H. J. Kim, Evolution of resonant self-interacting dark matter halos, Phys. Rev. D 109, 063535 (2024), arXiv:2304.12621 [astro-ph.CO] .
  * Colquhoun  _et al._ [2021] B. Colquhoun, S. Heeba, F. Kahlhoefer, L. Sagunski, and S. Tulin, Semiclassical regime for dark matter self-interactions, Phys. Rev. D 103, 035006 (2021), arXiv:2011.04679 [hep-ph] .
  * Zhong  _et al._ [2023] Y.-M. Zhong, D. Yang, and H.-B. Yu, The impact of baryonic potentials on the gravothermal evolution of self-interacting dark matter haloes, MNRAS 526, 758 (2023), arXiv:2306.08028 [astro-ph.CO] .
  * Yang  _et al._ [2024] D. Yang, E. O. Nadler, H.-B. Yu, and Y.-M. Zhong, A parametric model for self-interacting dark matter halos, J. Cosmol. Astropart. Phys. 2024, 032 (2024), arXiv:2305.16176 [astro-ph.CO] .
  * Mace  _et al._ [2024] C. Mace, Z. Carton Zeng, A. H. G. Peter, X. Du, S. Yang, A. Benson, and M. Vogelsberger, Convergence Tests of Self-Interacting Dark Matter Simulations, arXiv e-prints , arXiv:2402.01604 (2024), arXiv:2402.01604 [astro-ph.GA] .
  * Diemer and Joyce [2019] B. Diemer and M. Joyce, An Accurate Physical Model for Halo Concentrations, ApJ 871, 168 (2019), arXiv:1809.07326 [astro-ph.CO] .
  * Binney and Tremaine [1987] J. Binney and S. Tremaine,  _Galactic dynamics_ (1987).
  * Weinberger  _et al._ [2020] R. Weinberger, V. Springel, and R. Pakmor, The AREPO Public Code Release, ApJS 248, 32 (2020), arXiv:1909.04667 [astro-ph.IM] .
  * Springel and White [1999] V. Springel and S. D. M. White, Tidal tails in cold dark matter cosmologies, MNRAS 307, 162 (1999), arXiv:astro-ph/9807320 [astro-ph] .
  * Kazantzidis  _et al._ [2004] S. Kazantzidis, J. Magorrian, and B. Moore, Generating Equilibrium Dark Matter Halos: Inadequacies of the Local Maxwellian Approximation, ApJ 601, 37 (2004), arXiv:astro-ph/0309517 [astro-ph] .
  * Vogelsberger  _et al._ [2014] M. Vogelsberger, J. Zavala, C. Simpson, and A. Jenkins, Dwarf galaxies in CDM and SIDM with baryons: observational probes of the nature of dark matter, MNRAS 444, 3684 (2014), arXiv:1405.5216 [astro-ph.CO] .
  * Note [1] The value of $r_{\rm c}$ is chosen to reduce statistical uncertainties while keeping the region small enough to represent the halo central density.
  * Palubski  _et al._ [2024] I. Palubski, O. Slone, M. Kaplinghat, M. Lisanti, and F. Jiang, Numerical Challenges in Modeling Gravothermal Collapse in Self-Interacting Dark Matter Halos, arXiv e-prints , arXiv:2402.12452 (2024), arXiv:2402.12452 [astro-ph.CO] .
  * Shen  _et al._ [2022] X. Shen, P. F. Hopkins, L. Necib, F. Jiang, M. Boylan-Kolchin, and A. Wetzel, Dissipative Dark Matter on FIRE: II. Observational signatures and constraints from local dwarf galaxies, arXiv e-prints , arXiv:2206.05327 (2022), arXiv:2206.05327 [astro-ph.GA] .
  * Akerib  _et al._ [2017] D. S. Akerib, S. Alsum, H. M. Araújo, X. Bai, A. J. Bailey, J. Balajthy, P. Beltrame, E. P. Bernard, A. Bernstein, T. P. Biesiadzinski, E. M. Boulton, R. Bramante, P. Brás, D. Byram, S. B. Cahn, M. C. Carmona-Benitez, C. Chan, A. A. Chiller, C. Chiller, A. Currie, J. E. Cutter, T. J. R. Davison, A. Dobi, J. E. Y. Dobson, E. Druszkiewicz, B. N. Edwards, C. H. Faham, S. Fiorucci, R. J. Gaitskell, V. M. Gehman, C. Ghag, K. R. Gibson, M. G. D. Gilchriese, C. R. Hall, M. Hanhardt, S. J. Haselschwardt, S. A. Hertel, D. P. Hogan, M. Horn, D. Q. Huang, C. M. Ignarra, M. Ihm, R. G. Jacobsen, W. Ji, K. Kamdin, K. Kazkaz, D. Khaitan, R. Knoche, N. A. Larsen, C. Lee, B. G. Lenardo, K. T. Lesko, A. Lindote, M. I. Lopes, A. Manalaysay, R. L. Mannino, M. F. Marzioni, D. N. McKinsey, D. M. Mei, J. Mock, M. Moongweluwan, J. A. Morad, A. S. J. Murphy, C. Nehrkorn, H. N. Nelson, F. Neves, K. O’Sullivan, K. C. Oliver-Mallory, K. J. Palladino, E. K. Pease, P. Phelps, L. Reichhart, C. Rhyne, S. Shaw, T. A. Shutt, C. Silva, M. Solmaz, V. N. Solovov, P. Sorensen, S. Stephenson, T. J. Sumner, M. Szydagis, D. J. Taylor, W. C. Taylor, B. P. Tennyson, P. A. Terman, D. R. Tiedt, W. H. To, M. Tripathi, L. Tvrznikova, S. Uvarov, J. R. Verbus, R. C. Webb, J. T. White, T. J. Whitis, M. S. Witherell, F. L. H. Wolfs, J. Xu, K. Yazdani, S. K. Young, C. Zhang, and LUX Collaboration, Results from a Search for Dark Matter in the Complete LUX Exposure, Phys. Rev. Lett. 118, 021303 (2017), arXiv:1608.07648 [astro-ph.CO] .
  * Aprile  _et al._ [2018] E. Aprile, J. Aalbers, F. Agostini, M. Alfonsi, L. Althueser, F. D. Amaro, M. Anthony, F. Arneodo, L. Baudis, B. Bauermeister, M. L. Benabderrahmane, T. Berger, P. A. Breur, A. Brown, A. Brown, E. Brown, S. Bruenner, G. Bruno, R. Budnik, C. Capelli, J. M. R. Cardoso, D. Cichon, D. Coderre, A. P. Colijn, J. Conrad, J. P. Cussonneau, M. P. Decowski, P. de Perio, P. di Gangi, A. di Giovanni, S. Diglio, A. Elykov, G. Eurin, J. Fei, A. D. Ferella, A. Fieguth, W. Fulgione, A. Gallo Rosso, M. Galloway, F. Gao, M. Garbini, C. Geis, L. Grandi, Z. Greene, H. Qiu, C. Hasterok, E. Hogenbirk, J. Howlett, R. Itay, F. Joerg, B. Kaminsky, S. Kazama, A. Kish, G. Koltman, H. Landsman, R. F. Lang, L. Levinson, Q. Lin, S. Lindemann, M. Lindner, F. Lombardi, J. A. M. Lopes, J. Mahlstedt, A. Manfredini, T. Marrodán Undagoitia, J. Masbou, D. Masson, M. Messina, K. Micheneau, K. Miller, A. Molinario, K. Morâ, M. Murra, J. Naganoma, K. Ni, U. Oberlack, B. Pelssers, F. Piastra, J. Pienaar, V. Pizzella, G. Plante, R. Podviianiuk, N. Priel, D. Ramírez García, L. Rauch, S. Reichard, C. Reuter, B. Riedel, A. Rizzo, A. Rocchetti, N. Rupp, J. M. F. Dos Santos, G. Sartorelli, M. Scheibelhut, S. Schindler, J. Schreiner, D. Schulte, M. Schumann, L. Scotto Lavina, M. Selvi, P. Shagin, E. Shockley, M. Silva, H. Simgen, D. Thers, F. Toschi, G. Trinchero, C. Tunnell, N. Upole, M. Vargas, O. Wack, H. Wang, Z. Wang, Y. Wei, C. Weinheimer, C. Wittweg, J. Wulf, J. Ye, Y. Zhang, T. Zhu, and Xenon Collaboration, Dark Matter Search Results from a One Ton-Year Exposure of XENON1T, Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Kahlhoefer  _et al._ [2015] F. Kahlhoefer, K. Schmidt-Hoberg, J. Kummer, and S. Sarkar, On the interpretation of dark matter self-interactions in Abell 3827, MNRAS 452, L54 (2015), arXiv:1504.06576 [astro-ph.CO] .
  * Springel [2010] V. Springel, E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh, MNRAS 401, 791 (2010), arXiv:0901.4107 [astro-ph.CO] .
  * Eddington [1916] A. S. Eddington, The distribution of stars in globular clusters, MNRAS 76, 572 (1916).
  * Łokas and Mamon [2001] E. L. Łokas and G. A. Mamon, Properties of spherical galaxies and clusters with an NFW density profile, MNRAS 321, 155 (2001), arXiv:astro-ph/0002395 [astro-ph] .
  * Prada  _et al._ [2012] F. Prada, A. A. Klypin, A. J. Cuesta, J. E. Betancort-Rijo, and J. Primack, Halo concentrations in the standard $\Lambda$ cold dark matter cosmology, MNRAS 423, 3018 (2012), arXiv:1104.5130 [astro-ph.CO] .
  * Yang  _et al._ [2023b] S. Yang, X. Du, Z. C. Zeng, A. Benson, F. Jiang, E. O. Nadler, and A. H. G. Peter, Gravothermal Solutions of SIDM Halos: Mapping from Constant to Velocity-dependent Cross Section, ApJ 946, 47 (2023b), arXiv:2205.02957 [astro-ph.CO] .
  * Yang and Yu [2022] D. Yang and H.-B. Yu, Gravothermal evolution of dark matter halos with differential elastic scattering, J. Cosmol. Astropart. Phys. 2022, 077 (2022), arXiv:2205.03392 [astro-ph.CO] .