# Astrophysical Tests of Dark Matter Self-Interactions

Susmita Adhikari1, 2, Arka Banerjee1, Kimberly K. Boddy3, Francis-Yan Cyr-
Racine4,
Harry Desmond5, 6, 7, Cora Dvorkin8, Bhuvnesh Jain9, Felix Kahlhoefer10,
Manoj Kaplinghat11, Anna Nierenberg12, Annika H. G. Peter13,
Andrew Robertson14, Jeremy Sakstein8, 15, Jesús Zavala16
susmita<EMAIL_ADDRESS>
( 1 Department of Physics, Indian Institute of Science Education and Research,
Homi Bhaba Road, Pashan, Pune 411008, India
2 Department of Astronomy and Astrophysics, The University of Chicago,
Chicago, IL, 60637
3 Department of Physics, University of Texas at Austin, Austin, TX, 78712, USA
4 Department of Physics and Astronomy, University of New Mexico, Albuquerque,
New Mexico 87106, USA
5 Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road,
Oxford OX1 3RH, UK
6 McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, 5000 Forbes Ave, Pittsburgh, PA 15213
7 Institute of Cosmology $\&$ Gravitation, University of Portsmouth, Dennis
Sciama Building, Burnaby Road, Portsmouth, PO1 3FX, UK
8 Harvard University, Department of Physics, Cambridge, Massachusetts, 02138,
U.S.A.
9 Center for Particle Cosmology, Department of Physics and Astronomy,
University of Pennsylvania, 209 S. 33rd St., Philadelphia, PA 19104, USA
${}^{1}0$ Institute for Theoretical Particle Physics (TTP), Karlsruhe
Institute of Technology (KIT), 76128 Karlsruhe, Germany
11 Department of Physics and Astronomy, University of California, Irvine,
California 92697-4575, USA
12 Department of Physics, University of California, Merced, 5200 North Lake
Rd., Merced CA 95343
13 CCAPP, Department of Physics, and Department of Astronomy, The Ohio State
University, 191 W. Woodruff Ave., Columbus, OH 43210
14 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Drive, Pasadena, CA 91109, USA
15 Department of Physics & Astronomy, University of Hawaii, Watanabe Hall,
2505 Correa Road, Honolulu, HI 96822, USA
16 Center for Astrophysics and Cosmology, Science Institute, University of
Iceland, Dunhagi 5, 107 Reykjavik, Iceland )

###### Abstract

Self-interacting dark matter (SIDM) arises generically in scenarios for
physics beyond the Standard Model that have dark sectors with light mediators
or strong dynamics. The self-interactions allow energy and momentum transport
through halos, altering their structure and dynamics relative to those
produced by collisionless dark matter. SIDM models provide a promising way to
explain the diversity of galactic rotation curves, and they form a predictive
and versatile framework for interpreting astrophysical phenomena related to
dark matter.

This review provides a comprehensive explanation of the physical effects of
dark matter self-interactions in objects ranging from galactic satellites
(dark and luminous) to clusters of galaxies and the large-scale structure. The
second major part describes the methods used to constrain SIDM models
including current constraints, with the aim of advancing tests with upcoming
galaxy surveys. This part also provides a detailed review of the unresolved
small-scale structure formation issues and concrete ways to test simple SIDM
models. The review is rounded off by a discussion of the theoretical
motivation for self-interactions, degeneracies with baryonic and gravitational
effects, extensions to the single-component elastic-interactions SIDM
framework, and future observational and theoretical prospects.

###### Contents

  1. 1 Introduction
  2. 2 Theoretical considerations: from micro to macro
    1. 2.1 What do we mean by “cross section”?
    2. 2.2 Modeling the Effective Macroscopic Phenomena Arising from Microscopic Interactions
  3. 3 Physical Effects in Galaxies and Clusters
    1. 3.1 Core formation with SIDM
    2. 3.2 Gravothermal collapse of SIDM halos
    3. 3.3 Lack of Dynamical Friction in Cored Halos
    4. 3.4 Enhanced Gravitational Tidal Stripping of Cored Subhalos
    5. 3.5 Halo shapes
    6. 3.6 Drag Force and Mass Loss due to Self-Interactions
    7. 3.7 Large-scale structure of SIDM models
    8. 3.8 Are these physical effects unique?
      1. 3.8.1 Degeneracies with baryonic physics
      2. 3.8.2 Degeneracies with novel dark matter or gravitational physics
  4. 4 Observations of SIDM phenomenology
    1. 4.1 Strong Gravitational Lensing in Clusters, Groups and Large Ellipticals
      1. 4.1.1 The total deflector mass distribution
      2. 4.1.2 Constraints from resolved subhalos and line-of-sight halos
      3. 4.1.3 Constraints from the collective effect of unresolved subhalos and line-of-sight halos
    2. 4.2 Stellar streams in the Milky Way
    3. 4.3 X-ray and Weak Lensing Observations of Clusters, Groups and Large Ellipticals
    4. 4.4 Major and Minor Mergers in Groups and Clusters (and associated observables)
    5. 4.5 Incidence of Multiple BCGs, their Stellar Kinematics and their Spatial Separations
    6. 4.6 Dwarf Galaxies in the Local Group and beyond
    7. 4.7 Rotation curves of spiral galaxies
    8. 4.8 Ultra-diffuse galaxies
    9. 4.9 Warping of stellar disks
    10. 4.10 Large-scale structure constraints
    11. 4.11 Summary of constraints on the self-interaction cross section
  5. 5 Beyond the simple SIDM framework
    1. 5.1 SIDM with dissipation
    2. 5.2 Subcomponent SIDM
    3. 5.3 Phenomenology of light mediators
    4. 5.4 Connections to the Standard Model
  6. 6 Future Prospects
    1. 6.1 Observational probes
    2. 6.2 Simulations and Modeling
    3. 6.3 Promising tests of simple SIDM models

##  1 Introduction

Ideas for physics beyond the Standard Model (SM) often predict new stable
particles, which could be part or all of the dark matter (DM) in the Universe.
Two well-known examples of such ideas introduce new symmetries, Supersymmetry
and the Peccei-Quinn symmetry, to alleviate problems with the SM. The
predicted DM particles, neutralinos and axions, are canonical examples of the
Cold Dark Matter (CDM) model wherein DM is born non-relativistic and only its
gravitational interaction is relevant for structure formation. CDM models
yield predictions for the evolution of the homogeneous and inhomogeneous
evolution of the Universe that agree very well with observations on large
scales [1, 2].

Collider and deep underground searches have steadily ruled out large swaths of
the parameter space of these new symmetry models. As this search progresses,
another front has opened up. In addition to new symmetries, physics beyond the
SM could include new forces among DM particles. This is a generic consequence
of a dark sector. Dark sectors play an important role in current model
building efforts in particle physics. The new particles in the dark sector do
not have to be heavy, but they are very weakly coupled to the SM. Dark sectors
also allow for qualitatively new production mechanisms for DM, such as freeze-
in [3]. The SIDM production mechanisms offer concrete benchmark models that
may be tested through laboratory experiments or astrophysical probes (e.g. [4,
5, 6, 7]).

A generic consequence of a new force is the ability for DM particles to
scatter off of each other. If the force mediator is light, then the resulting
dark sector interactions could leave observable signatures on galactic and
sub-galactic scales with a cross section that is generically velocity-
dependent (e.g., Refs. [8, 9, 10, 11]). However, a light mediator is not a
necessary ingredient of DM models with a large self-interaction cross section
that is velocity-dependent (e.g., Refs. [12, 13, 14]). In this review, we
discuss what happens if DM is assumed to have interactions with itself (self-
interactions) beyond gravitational interactions. If dark-matter particles have
a non-trivial probability of interacting on $\sim$ Gyr timescales, this will
allow energy and momentum to flow from one part of the dark matter halo to
another beyond what is enabled by gravity. As we will highlight in this
review, the introduction of DM scattering has profound implications for the DM
distribution within individual halos and in the hierarchical assembly of
structure on non-linear scales. Furthermore, the types of particle physics
models that admit strong DM self-scattering could also lead to imprints on the
DM power spectrum. Thus, self-interacting dark matter (SIDM) phenomenology
includes deviations from CDM on scales of individual DM halos and subhalos, as
well as their population statistics.

The initial interest in SIDM models began more than twenty years ago as a
response to two puzzles in observational astronomy: the shape of the rising
part of the rotation curves of spiral galaxies (the “cusp/core problem” [15,
16, 17, 18]) and the number of observed satellites of the Milky Way (MW; the
“missing satellites problem” [19, 20]). The term “self-interacting dark
matter" had been coined almost a decade earlier [21], but it took these
astrophysical puzzles to put dark-matter models with strong self-scattering
firmly in the view of particle physicists and astronomers.

The cusp/core problem refers to the tension between the predicted density
profiles of DM halos and those inferred from 21 cm observations. In the
absence of strong feedback from star formation physics (e.g., [22]), DM halos
are expected to be _cusped_ , with central density profiles $\rho\propto
r^{-1}$ [23]. However, the observationally inferred rotation curves of low-
surface-brightness (LSB) spiral galaxies almost universally indicate shallower
density profiles for their galaxies’ halos (e.g., [18, 24]). SIDM was invoked
as a way to produce isothermal, flat cores ($\rho\sim\hbox{ const}$ at small
radii) in the halos hosting spiral and dwarf spheroidal galaxies [25].

The missing satellites problem pertains to the most striking prediction of the
CDM paradigm: a hierarchy of dark-matter halos down to a minimum mass that
could be about Earth-mass [26, 27, 28, 29]. The halo mass function scales
approximately as $dN/dM\propto M^{-1.9}$ [30, 31, 32], which means there are
many more low-mass halos than higher mass ones. We can use galaxies to trace
the halo mass function [33] down to a mass scale at which gas cannot cool and
form stars [34, 35, 36, 37, 38, 39, 40]. A comparison of $N$-body simulation
predictions for the CDM subhalo mass function to the observed MW satellites
circa 2000 revealed a large discrepancy, which was dubbed the missing
satellites problem. This provided motivation for SIDM models since self-
interactions can suppress the abundance of subhalos [25]. Interactions between
DM particles in the host halo and a subhalo can kick particles out of the
subhalo, which, combined with the increased tidal disruption of cored subhalos
compared to cuspy subhalos, leads to fewer MW satellites compared with CDM
cosmologies.

The original motivations for SIDM—the cusp/core problem and the missing
satellites problem—have undergone significant evolution in the past twenty
years because of both new observations and simulations. The cusp/core problem
still exists—the inferred dark matter densities of many low surface brightness
and dwarf spheroidal galaxies are intriguingly low. Additionally, there is a
high degree of galaxy-to-galaxy scatter in the shapes of rotation curves
(e.g., [41, 42, 43, 44]). It is unclear if the physics of star formation can
create the diversity of density profiles observed, and it appears that SIDM
could accommodate the range of observed rotation curves [45, 46, 47, 48, 49,
50, 51]. This diversity is most apparent in galaxies with rotation velocities
of about $70-100\,\mathrm{km\,s^{-1}}$, but there is some evidence for
deviation from the standard CDM cuspy halo profile even on cluster scales [52,
53, 54]. Thus, instead of narrowly focusing on a narrow range of galaxy
masses, SIDM models should be challenged with explaining halo density profiles
spanning many decades in halo mass from those of dwarf galaxies to clusters of
galaxies [55], as we will discuss in this review.

The missing satellites problem has largely gone away, and the increasing
numbers of confirmed satellite galaxies of the MW imply that this galaxy
population can be used to place tight constraints on novel DM physics [56,
57]. Since the advent of wide-field digital sky surveys, dozens of faint dwarf
galaxies have been identified as satellites of the MW [58, 59, 60, 61, 62, 63,
64]. The abundance of satellites is consistent with CDM predictions under
reasonable and simulation-tested models for star formation before and after
reionization [35, 36, 37, 65, 66, 67, 68, 69, 70, 71, 72, 73]. In the context
of constant cross section SIDM models, evaporation (mass loss) due to
scattering events has been found to be inefficient at altering the subhalo
mass function unless the SIDM cross section is so high as to have already been
ruled out [74, 75]. Thus, to zeroth order, SIDM and CDM models share the same
issues or successes when it comes to the census of satellites.

In this review, we focus on the astrophysical phenomenology of SIDM across the
range of scales from unresolved subhalos to clusters of galaxies. In Sec. 2,
we provide the theoretical background to SIDM and the ways in which it is
modelled. Sec. 3 provides an exhaustive account of the physical effects of
self-interactions in galaxies, clusters and large-scale structure. Sec. 4
describes the ways in which these effects may be searched for, including hints
of the existence of self-interactions and the discovery potential of near-term
astronomical facilities. We focus predominantly on elastic interactions
because these typically dominate over inelastic interactions in particle
models, and most simulation work has been done in that context. We provide a
summary of the current constraints on the elastic self-interaction cross
section at the end of this section. In Sec. 5 we connect to specific particle
models and describe ways in which more complex interactions and effects may
extend the SIDM paradigm. We close with our view on the theoretical and
observational path forward to characterizing DM self-scattering in the next
decade, highlighting important issues that need more work and promising future
tests (Sec. 6).

The present review is complementary to a recent treatise by Tulin & Yu [76] on
DM self-interactions. Together, these reviews provide a more complete picture
of the theory, phenomenology, and current and future observational tests of
SIDM models. The present review also provides a self-contained summary of the
small-scale structure challenges, but readers will benefit from referring to
the influential review on this topic in Ref. [77].

This review forms the second arm of the Novel Probes Project, an initiative to
nurture the burgeoning field of astrophysical tests of the dark sector. The
first describes astrophysical tests of extended theories of gravity (Baker et
al. [78]). The website of the Novel Probes Project is
https://www.novelprobes.org.

##  2 Theoretical considerations: from micro to macro

We will be interested in signatures of the microphysical DM scattering process
on macrophysical scales, typical ${\cal O}(100)\leavevmode\nobreak\ \rm pc$ or
larger. In this section, we define various quantitative measures of the
interaction strength on the microscopic scale—the differential cross section,
the total cross section, and the momentum and viscosity cross section—and
outline the typical range of values of the cross section for which the effects
of these interactions become interesting from the perspective of structure
formation. In the second part, we focus on how the effects of the
microphysical interactions are modeled in simulations of structure formation.
Typically, these $N$-body simulations follow the evolution of particles with
masses many times $M_{\odot}$, so it is important to clarify the connections
between the micophysical cross sections, and the modification of trajectories
of the macrophysical particles in simulations. We present various methods that
have been used in the literature to perform this mapping, as well as a range
of tests which have been performed to validate these implementations.

###  2.1 What do we mean by “cross section”?

In this section, we tackle several issues related to the definition of a cross
section and show that the microscopic definition is not the most relevant one
to describe macrophysical phenomena. We show what range of cross sections is
astrophysically interesting and what types of microphysical models give cross
sections in this range.

Historically, most of the discussion of SIDM was framed in the context of
velocity-independent “hard-sphere" scattering [25, 76] where the outgoing
momentum direction is random in the center-of-mass frame. This simple setup
can be parameterized by one parameter, the total self-scattering cross section
$\sigma$.

We typically, however, discuss the quantity $\sigma/m_{\text{DM}}$, the DM–DM
cross section per unit DM particle mass $m_{\text{DM}}$. For a DM particle
moving at velocity $v_{0}$ through a background of stationary DM particles
with a number density $n$, the rate at which that particle scatters with
background particles is

$R=\sigma\,n\,v_{0}=\frac{\sigma}{m_{\text{DM}}}\rho\,v_{0}\;,$ (1)

where $\rho$ is the mass density (henceforth “density”) of the background
particles, such that $n=\rho/m_{\text{DM}}$. The total probability for the
moving particle to scatter is therefore given by

$p=1-\exp\left(-\frac{\sigma}{m_{\text{DM}}}\Sigma\right)\;,$ (2)

where $\Sigma\equiv\int v_{0}\,\rho\,\mathrm{d}t=\int\rho\,\mathrm{d}x$ is the
integrated background density through which the particle passes. Placing this
scenario in the context of virialized halos, the integrated background density
is given by the column density of the halo multiplied with the number of times
the DM particle passes through the halo. The former quantity is $\Sigma\sim
1\,\mathrm{g\,cm^{-2}}$ on cluster scales (halo mass $\sim 10^{15}M_{\odot}$,
virial velocity $\sim 10^{3}\hbox{ km}/\hbox{s}$) and somewhat smaller for
smaller halos [79], whereas the latter quantity can be estimated by the ratio
of the Hubble time and the dynamical time at the virial radius and is found to
be of order one. Thus, we must have $\sigma/m_{\text{DM}}\sim
1\,\mathrm{cm^{2}\,g^{-1}}$ in order for scatters to be frequent enough to
induce large astrophysical effects.

In realistic SIDM models, the cross section will generically change with the
relative velocity $v$ of the colliding DM particles. However, constant cross
section models are also possible in the case of DM with a finite size (the
“black disk” scattering limit), which can arise if DM is a composite particle
formed by a QCD-like confining force, like neutron–neutron scattering [80].
More prosaically, we expect velocity dependence if two mass scales are
relevant for the scattering process, for example, a DM and a mediator mass, or
the dark proton and electron masses in the case of atomic DM.

The best-known example of velocity-dependence is a long-range Yukawa
interaction, arising from the exchange of a light scalar or vector mediator.
In the massless limit, we recover the standard Rutherford scattering cross
section that scales as $v^{-4}$. If the mediating particle is light but not
massless, the differential scattering cross section with respect to the
scattering angle $\theta$ is

$\frac{\text{d}\sigma}{\text{d}\cos\theta}=\frac{\sigma_{0}}{2\left(1+\frac{v^{2}}{w^{2}}\sin^{2}\frac{\theta}{2}\right)^{2}}\;,$
(3)

where $\sigma_{0}$ denotes the total cross section for $v\to 0$ and
$w=m_{\text{med}}/m_{\text{DM}}$ denotes the ratio of mediator mass to DM
mass. Scattering is velocity-independent for $v\ll w$ and is suppressed as
$v^{-4}$ for large velocities. Equation (3) holds under the Born
approximation, in which the strength of the Yukawa potential
$\alpha_{\textrm{DM}}$ is small. If, however, we are in the regime
$w\ll\alpha_{\textrm{DM}}$, the computation of the differential cross section
is no longer perturbative in $\alpha_{\textrm{DM}}$ and must be calculated by
numerically solving the nonrelativistic Schödinger equation for the Yukawa
potential [9, 10, 11, 81]. In the low-velocity regime of $v\lesssim w$, the
presence of a bound-state spectrum permits a resonant enhancement of the DM
self-scattering cross section, making it possible to obtain SIDM from a
weakly-coupled theory. If, on the other hand, DM particles are strongly-
coupled, they may also form permanent bound states [82, 83]. While more
complicated to study, these models are also expected to exhibit velocity-
dependent self-interactions [84, 85, 86].

If the scattering cross section exhibits a velocity dependence, it is useful
to define the velocity-averaged cross section

$\langle\sigma\rangle=\frac{\int f_{\rm rel}(v_{\rm rel})\,\sigma(v_{\rm
rel})\,v_{\rm rel}^{n}\mathrm{d}v_{\rm rel}}{\int f_{\rm rel}(v_{\rm
rel})\,v_{\rm rel}^{n}\mathrm{d}v_{\rm rel}}\;,$ (4)

where $f_{\rm rel}(v_{\rm rel})$ denotes the one-dimensional distribution of
relative velocities in the system under consideration, and the integral is
weighted by a factor of $v^{n}$. The effect of a velocity dependence is then
approximately captured by replacing $\sigma\to\langle\sigma\rangle$ in the
discussion above for an appropriate choice of $n$ (discussed below). The
single-particle velocity distribution $f(v)$ is typically taken to be a
Maxwell–Boltzmann distribution, which is a good approximation within the
thermalized core but may not work as well in the outskirts of the halo:

$f(v)\propto v^{2}\exp\left(-\frac{v^{2}}{2\sigma_{v}^{2}}\right)$ (5)

with the one-dimensional velocity dispersion $\sigma_{v}$ of the order of the
virial velocity. The distribution of relative velocities is then given by a
Maxwell–Boltzmann distribution with velocity dispersion $\sqrt{2}\sigma_{v}$.
For an even rougher estimate, it is often sufficient to approximate
$\langle\sigma\rangle\approx\sigma(v=\langle v_{\rm rel}\rangle)$ where
$\langle v_{\rm rel}\rangle=4/\sqrt{\pi}$ for a Maxwell–Boltzmann
distribution.

For a velocity weighting of $n=1$ in Eq. (4), the substitution
$\sigma\to\langle\sigma\rangle$ appropriately captures the average rate that
self interactions occur within the halo. One of the main features of SIDM,
however, is that the self-scattering of DM particles permits heat flow from
one region of the halo to another. In order to capture this effect, a
different value of $n$ is needed. Under the assumption of a Maxwell–Boltzmann
equilibrium distribution, the amount of momentum and energy transferred in a
single collision can be used to obtain momentum- and energy-transfer rates,
which are proportional to Eq. (4) for $n=2$ and $n=3$, respectively [87],
reflecting the fact that the transferred momentum is proportional to $v$ and
the transferred energy to $v^{2}$. A more robust treatment of the velocity
dependence exists in the fluid regime, in which the DM distribution function
experiences a small departure from equilibrium. Integrating the non-
relativistic collisional Boltzmann equation at linear order [88, 89] yields an
energy flux that is proportional to Eq. (4) for $n=5$ [90].

From a phenomenological perspective, velocity-dependent self-interactions are
interesting because they strongly modify the relative importance of different
types of astrophysical systems. If, as in the example above, scattering is
suppressed for large velocities, the sensitivity of very massive systems such
as galaxy clusters (where relative velocities are typically large) is reduced
compared to less massive systems like galaxies, because $\langle\sigma\rangle$
is predicted to be much smaller for the former than for the latter. Velocity-
dependent interactions have, therefore, frequently been considered to evade
the strong constraints on $\sigma$ that arise from observations of galaxy
clusters [9, 10, 55], which will be discussed in more detail below.

Equation (3) exhibits another interesting property: the velocity dependence is
accompanied by an angular dependence. While scattering is isotropic for $v\ll
w$, the differential cross section is strongly peaked towards $\theta\to 0$
for larger velocities. The combination of a velocity dependence and an angular
dependence is in fact very generic, as they often arise simultaneously from a
fundamental dependence on the momentum transfer $q$, which is given by
$q^{2}=2\,m_{\textrm{DM}}^{2}\,v^{2}(1-\cos\theta)$. In other words, if
scattering with large velocity is suppressed, there is usually also a
suppression of large scattering angles.

If DM self-scattering is not isotropic, it is no longer clear that the total
cross section $\sigma=\int(\mathrm{d}\sigma/\mathrm{d}\theta)\mathrm{d}\theta$
is the relevant quantity for calculating effects in astrophysical systems. In
order for the phase space density of DM to change on account of scattering, we
must consider how scattering affects the exchange of momentum and energy
between particles. In fact, the more traditional quantity to use is the
momentum transfer cross section

$\sigma_{\mathrm{T}}=\int\frac{\mathrm{d}\sigma}{\mathrm{d}\theta}(1-\cos\theta)\mathrm{d}\theta\;,$
(6)

which down-weights forward-angle scattering (which does not appreciably change
the momentum of the scattering particles) and up-weights backward-angle
scattering (which yields the largest momentum transfer). In analogy to
$\langle\sigma\rangle$, one can then define the velocity-averaged momentum
transfer cross section $\langle\sigma_{\mathrm{T}}\rangle$ [11, 81].

This definition of the momentum transfer cross section is, however, not self-
consistent for the scattering of identical particles, as there should be no
physical difference between $\theta=0$ and $\theta=\pi$. To address this
issue, a modified definition is [91].

$\tilde{\sigma}_{\mathrm{T}}=\int\frac{\mathrm{d}\sigma}{\mathrm{d}\theta}(1-|\cos\theta|)\mathrm{d}\theta\;.$
(7)

Indeed, explicit numerical simulations show this latter definition is more
appropriate for a differential cross section of the form of Eq. (3)
[92].111This conclusion may be different for the scattering of non-identical
particles, e.g. for particle-antiparticle scattering [93]. Alternatively,
Refs. [81, 86] suggest considering the viscosity cross section

$\sigma_{\mathrm{V}}=\int\frac{\mathrm{d}\sigma}{\mathrm{d}\theta}\sin^{2}\theta\mathrm{d}\theta\;,$
(8)

which weighs the cross section by the fractional transverse energy transfer.
Simulations show that $\sigma_{V}$ is also more appropriate than $\sigma_{T}$
for Rutherford scattering in Eq. (3), as well as for Møller scattering [94].
Both the modified momentum transfer cross section and the viscosity cross
section exhibit the appropriate symmetry: they suppress forward- and backward-
angle elastic scattering, which are identical processes for identical
particles.

For numerical simulations of isotropic systems, such as isolated DM halos, it
has been shown that it is not necessary to explicitly include the angular
dependence of the scattering in the simulation, as long as the velocity
dependence is correctly included. This is also true when considering average
properties of entire populations, such as stacked weak lensing profiles or
stacked satellite counts around galaxy clusters [95]. For systems with a
preferred direction, such as merging galaxy clusters, however, the detailed
angular dependence does become relevant [92]. In other words, it is no longer
possible to capture all features of DM self-interactions from the momentum
transfer cross sections (or other variants of the cross section), since the
angular information has been lost in the integration over angles. We will
return to the effects of an angular dependence in Sec. 3.6.

In summary, what we mean by “the cross section” depends on the combination of
the type of microphysical model and the astrophysical application. Our guiding
principle is capturing the flow of particles in momentum and energy space
arising from interactions. For all astrophysical applications, though, we
expect the relevant cross section to be of order $\sigma/m_{\text{DM}}\sim
1\leavevmode\nobreak\ \,\mathrm{cm^{2}\,g^{-1}}$ in order for there to be a
significant effect for halos, which typically have DM column densities of
order $\Sigma\sim 1\,\mathrm{g\,cm^{-2}}$.

###  2.2 Modeling the Effective Macroscopic Phenomena Arising from Microscopic
Interactions

For most systems of interest in an SIDM cosmology, cross sections
$\mathcal{O}(1)\,\mathrm{cm^{2}\,g^{-1}}$ imply that they fall in an
interesting regime that is an intermediary between collisionless
dynamics—simulated using $N$-body methods—and fully collisional hydrodynamics
typically simulated with either smoothed particle hydrodynamics (e.g. Ref.
[96]) or grid-based methods (e.g. Ref. [97]). To see this, we note that the
mean free path —the average distance a particle travels between scattering—is

$\lambda=\frac{v_{0}}{R}=\frac{1}{\frac{\sigma}{m_{\text{DM}}}\rho}\;,$ (9)

which, evaluated for the DM-density in the solar-neighborhood ($\sim
0.4\,\mathrm{GeV\,cm^{-3}}\approx 10^{7}\,\mathrm{M_{\odot}\,kpc^{-3}}$ [98]),
is $\sim 0.5\,\mathrm{Mpc}$ for
$\sigma/m_{\text{DM}}=1\,\mathrm{cm^{2}\,g^{-1}}$. Thus, SIDM scattering is
neither frequent enough that the distribution of particles can be described by
fluid elements with isotropic pressure, nor infrequent enough to be ignored.
For spherically symmetric systems, the effects of these infrequent
interactions can be modelled by considering energy transport, where the
thermal conductivity depends on $\lambda$ (cf. Refs. [99, 100, 101, 102, 103,
90, 94, 104]). However, most work on astrophysical SIDM phenomena is performed
in the context of numerical simulations. The solution has been to use standard
$N$-body methods, with the addition of Monte Carlo scattering to account for
self-interactions, which we now describe further.

The generalization of equation (1) for the interaction rate to the case of a
non-stationary background is

$R(\vec{r},\vec{v})=\int
f_{v}(\vec{r},\vec{v}^{\prime})\frac{\sigma}{m_{\text{DM}}}\rho(\vec{r})\,|\vec{v}-\vec{v}^{\prime}|\,\text{d}^{3}\vec{v}^{\prime}\;,$
(10)

where $f_{v}$ is the velocity distribution function, here normalized such that
$\int f_{v}(\vec{r},\vec{v})\,\text{d}^{3}\vec{v}=1$. In order to determine
whether a particle in an $N$-body simulation should scatter, the quantities
$f_{v}$ and $\rho$ must first be estimated from the distribution of nearby
particles, and the quantity $R$ multiplied by the length of the particle’s
time-step to produce a probability for scattering. In most implementations,
the value of the interaction probability computed in this way is then compared
to a random number generated from a uniform distribution between $0$ and $1$,
and if the computed probability exceeds the random number, a scatter takes
place.

For velocity-independent elastic scattering, especially, treating the
scattering of simulation particles (which may represent $\sim 10^{60}$ or more
“real” DM particles) in the same fashion as microphysical scattering is a
sensible assumption. In other words, one can visualize the scattering of
simulation particles as a scaled up version of the scattering of real
particles. One reason for this is that the astrophysically relevant quantity
is the cross section per unit mass, rather than the cross section itself.
Second, it can be shown that the actual evolution of the fine-grained phase
space density of microphysical DM particles is mimicked by that of the coarse-
grained phase space density (the phase-space density averaged over many
particles) for macroscopic collections of particles. This suggests that the
overall flow of energy and momentum is correctly captured by this procedure.
It is important to remember that the transfer of energy and momentum among
particles motivated the discussion of the different types of cross sections
($\sigma_{\rm T}$ and $\sigma_{\rm V}$) in the previous section. The choices
presented in that section were also a result of the underlying assumption that
the interactions of simulation particles are completely analogous to the
scattering of microscopic particles. More complicated self-interaction models,
e.g. those with non-trivial angular or velocity dependence in the differential
cross section, have also been simulated in an exactly analogous manner, but
with appropriate use of the momentum transfer cross section or the viscosity
cross section.

In the context of velocity-independent elastic scattering, a number of early
SIDM simulations first calculated the probability for a simulation particle to
scatter, and then chose a neighboring partner with which to scatter [105].
Most recent simulations have implemented scattering on a pair-by-pair basis
where the probability for individual pairs of nearby particles to scatter is
calculated, and a random number drawn for each pair to see if they do [74,
106, 107, 95].222A notable exception are Refs. [108, 109, 110], which consider
the case of very frequent self-interactions leading to an effective drag
force. However, also in this case it is necessary to consider pairs of
simulation particles in order to re-add the energy lost due to the drag force
and ensure overall energy conservation. This ensures not just that particles
scatter at the correct rate, but that the particles they scatter with
correctly sample the local velocity distribution. The probability of two
nearby particles, $i$ and $j$, scattering within the next time step, $\Delta
t$, is given by

$P_{ij}=\frac{\sigma}{m_{\text{DM}}}\rho_{ij}|\vec{v}_{i}-\vec{v}_{j}|\Delta
t\;,$ (11)

where $\rho_{ij}$ is the contribution of particle $j$ to the density estimate
at the location of particle $i$. Modern implementations of SIDM within
$N$-body simulations have differed primarily in the form of $\rho_{ij}$. For
example:

$\displaystyle\mathrm{Vogelsberger}\leavevmode\nobreak\
\cite[cite]{[\@@bibref{}{Vogelsberger:2012ku}{}{}]}$ $\displaystyle:$
$\displaystyle\rho_{ij}=m_{\mathrm{p}}W(r_{ij},h_{\mathrm{SI}})$ (12)
$\displaystyle\mathrm{Rocha}\leavevmode\nobreak\
\cite[cite]{[\@@bibref{}{Rocha13}{}{}]}$ $\displaystyle:$
$\displaystyle\rho_{ij}=m_{\mathrm{p}}\int
W(|\vec{x}|,h_{\mathrm{SI}})W(|\vec{x}+\vec{r}_{ij}|,h_{\mathrm{SI}})\,\text{d}^{3}\vec{x}$
(13) $\displaystyle\mathrm{Robertson}\leavevmode\nobreak\
\cite[cite]{[\@@bibref{}{Robertson17BC}{}{}]}$ $\displaystyle:$
$\displaystyle\rho_{ij}=\begin{cases}m_{\mathrm{p}}/\frac{4}{3}\pi
h_{\mathrm{SI}}^{3},&r_{ij}\leq h_{\mathrm{SI}}\\\
0,&r_{ij}>h_{\mathrm{SI}}\end{cases}$ (14)

where

$W(r,h)=\frac{8}{\pi
h^{3}}\begin{cases}1-6(\frac{r}{h})^{2}+6(\frac{r}{h})^{3},&0\leq\frac{r}{h}\leq
1/2\\\ 2(1-\frac{r}{h})^{3},&1/2<\frac{r}{h}\leq 1\\\
0,&1<\frac{r}{h}\>.\end{cases}$ (15)

The other main difference between recent SIDM implementations has been the
choice of $h_{\mathrm{SI}}$. Ref. [74] used a variable $h_{\mathrm{SI}}$ that
adapted to the local density to keep the number of neighbours $\sim 38$, while
both Refs. [106] and [107] used $h_{\mathrm{SI}}$ that was fixed for all
particles. Ref. [106] found that $h_{\mathrm{SI}}\gtrsim
0.2\,(\rho/m_{\mathrm{p}})^{-1/3}$ was required for scattering to be correctly
implemented (i.e. $h_{\mathrm{SI}}$ must be larger than 20% of the mean inter-
particle separation), and set $h_{\mathrm{SI}}$ such that this was the case at
densities a few times lower than the lowest densities for which self-
interactions should be significant. However, Ref. [107] (who used a smaller
$h_{\mathrm{SI}}$, similar to the Plummer-equivalent gravitational softening
length) demonstrated that scattering can be correctly implemented when
$h_{\mathrm{SI}}\ll 0.2\,(\rho/m_{\mathrm{p}})^{-1/3}$, and that the result in
Ref. [106] was likely a case of pairwise scattering probabilities saturating
above 1.333For more details see Appendix A in Ref. [107]. Ref. [107] also
compared their fixed $h_{\mathrm{SI}}$ approach to that adopted by Ref. [74]
and found the results to be very similar.

While a detailed study comparing the differences between SIDM implementations
has not been performed, large differences are not expected. Given the
inherently stochastic nature of SIDM interactions, whether the interaction
probability is smooth (equations (12) and (13)) or discontinuous (equation
(14)) with respect to the particle separation, and exactly how that
probability varies with separation should be unimportant. Also, large
$h_{\mathrm{SI}}$ results in many neighbors being found, but the probability
of scattering from a particular one being small (in all cases
$\rho_{ij}\propto 1/h_{\mathrm{SI}}^{3}$), while reducing $h_{\mathrm{SI}}$
essentially moves the stochasticity from the drawing of an unlikely random
number, to the unlikely event of finding a particle within the search region.

A number of authors have performed tests of their simulation codes, both to
test the results against analytical predictions, and to check for numerical
convergence. One test for which the expected outcome can be calculated
analytically is the so-called _wind tunnel_ test, in which particles are sent
through a uniform field of stationary background particles. The rate at which
particles scatter from this background can be calculated from equation (1),
and the distribution of speeds and directions of scattered particles can be
calculated by converting the differential scattering cross section from the
center-of-mass frame of each collision, to the rest-frame of the simulation
(typically the rest-frame of the background particles). Refs. [106] and [107]
both presented this test, with results that matched the analytical predictions
(see e.g. Fig. 1).

Figure 1: This set of figures from [106] shows the results from a test without
gravity where a DM sphere is moving with uniform velocity with respect to a
uniform background using the formalism described in Eq. 13. The top row shows
the distribution in $\theta$ after scattering compared to the theory
prediction for single elastic scatterings. The right panel shows better
agreement when multiple scattering events are removed. The bottom row shows
the distribution after scattering in speed and azimuthal angle. Note the
excess in the speed distribution due to multiple scattering events.

Another test with an analytical prediction is the rate of DM scattering in a
DM halo with a known distribution function. A good example that has been used
in the literature is a Hernquist profile [111], which has a similar density
profile to the NFW profile found in CDM simulations [23], except for a sharper
fall off in density at large radius. This profile has an analytical
distribution function, which allows for an analytical calculation of the
expected scattering rate as a function of radius and also means that
simulation initial conditions (which in the absence of DM self-interactions
would have a stable density profile as a function of time) can be easily
generated. A test with the Hernquist profile was proposed and used by
Vogelsberger et al. [74] to validate their choice of scattering kernel (Eq.
12).

Figure 2: Density profile of a MW-like DM halo (Aquarius Aq-A) at three
resolution levels for both CDM (black) and isotropic SIDM with a cross section
of 10 ${\rm cm}^{2}/{\rm g}$ (red). The resolution levels correspond to
particle masses of $4.9\times 10^{4}\,\mathrm{M}_{\odot}$ (solid), $3.9\times
10^{5}\,\mathrm{M}_{\odot}$ (dashed), and $3.1\times
10^{6}\,\mathrm{M}_{\odot}$ (dotted). The vertical lines mark the
gravitational softening length (2.8 times the Plummer-equivalent softening
length). The figure is adapted from Ref. [74].

Next, we discuss the question of numerical convergence, i.e. whether the
profiles are converged with respect to a change in the mass resolution of the
simulations. Given the stochastic nature of SIDM scattering, it would be
natural to think that a large number of particles is required in SIDM
simulations, so that this stochastic interaction is adequately sampled. Given
much of the interest around SIDM focuses on the effects it has on particle
density profiles of halos, it is important to determine the minimum radius to
which the density profiles from simulations can be trusted. Interestingly, the
density profiles of simulated SIDM halos are converged with respect to
simulation particle mass down to smaller radii than their CDM counterparts, as
illustrated in Fig. 2. The reason for this is that numerical effects in CDM
simulations typically form cores, so when the physical model itself is one
where a core is expected to form, these numerical effects are less
significant. In fact, one of the drivers of non-convergence in the inner
regions of CDM halos is the spurious two-body gravitational scattering between
particles.444These interactions are unphysical, because the importance of
these gravitational scattering events decreases with decreasing particle mass,
such that they should be irrelevant for the Universe’s DM fluid which is
presumed to be made up of a very large number of particles with very low
masses compared with the mass of simulation particles. The cross section for
these gravitational scattering events is a dramatic departure from the correct
behaviour of CDM, but adds only a small perturbation to SIDM which has a high
physical cross section for DM–DM scattering.

Finally, we briefly discuss the limitations of the Monte Carlo implementation
of self-interactions that has been described above. As mentioned at the
beginning of this section, this implementation is tailored to work in the
intermediate regime between collisionless dynamics of CDM, and the fully
collisional hydrodynamic limit. Therefore, for those models of self-
interactions where the interaction rate can be much higher (typically in cases
where the average momentum transfer is low), a different implementation is
needed to correctly describe the evolution of the system. As a concrete
example, consider a differential cross section such that DM particles
experience a very large number of collisions, each with tiny scattering
angles. The resulting effect is that DM particles lose momentum when moving
through a sea of other DM particles; the momentum loss of the infalling
particle is converted into an increase of energy of the surrounding particles,
in the sense that the bulk motion of particles (i.e., the bulk kinetic energy)
is turned into heat. In spirit, this is similar to the effect of dynamical
friction, but can have a different velocity dependence depending on the
details of the interaction cross section (e.g. [91, 112]). The implementation
of such interactions in simulations depend on the range of the interactions.
For long range interactions, which introduce an additional $1/r^{2}$ force on
DM particles, the interactions can be implemented using the same framework
used for computing the gravitational forces between simulation particles
[113], in a completely deterministic manner. For short range interactions,
different implementations are being explored. For example, [114] have
attempted to use smoothed particle hydrodynamics methods to capture SIDM
effects in models with frequent short range interactions, while [115] used an
‘effective drag force’ approach. This area remains a field of active research,
and is important to extend the reach of using simulations to explore a wider
range of self-interaction models.

##  3 Physical Effects in Galaxies and Clusters

We begin this section with a discussion of the mechanism by which DM self-
interactions alter the density profile of a DM halo: the transport of heat
from hot to cold regions of the halo. In §3.1 we show that this typically
reduces the central density of DM halos, although things become more
complicated when halos are baryon-dominated in their inner regions. This
reduction in central density is only transitory, and in §3.2 we discuss the
ultimate fate of SIDM halos, a process known as “gravothermal collapse”. We
then go on to discuss some consequences of these altered DM density profiles,
including the lack of dynamical friction for objects orbiting near the centre
of cored SIDM halos (§3.3), and the enhanced tidal stripping of SIDM subhalos
due to their lower binding energies (§3.4).

We then look beyond the radial density profiles, and discuss how SIDM affects
the shapes of halos (i.e. how they depart from being spherically symmetric) in
§3.5, with the primary effect being that halos become rounder with SIDM,
although this is also complicated when baryons are considered. Then in §3.6 we
look at the effects that inter-halo DM self-interactions have in merging
systems, focusing on the drag force induced by the exchange of momentum
through self-interactions as well as the mass loss that happens when high-
velocity collisions unbind particles from their halos.

Particle physics models that give rise to significant self-interactions can
also affect the evolution of DM density fluctuations in the early universe.
These changes to the matter power spectrum can affect the Cosmic Microwave
Background (CMB) as well as the number density of low-mass halos, and we
present current constraints in §3.7. Finally in §3.8 we discuss degeneracies
that exist between DM self-interactions and other physical processes that
could have similar effects. These include the baryonic physics relevant for
galaxy formation, departures from cold and collisionless DM other than SIDM,
and modifications to General Relativity.

###  3.1 Core formation with SIDM

Halos grow through hierarchical mergers in SIDM models (as in CDM models). The
interactions between DM particles allows for energy transfer from one region
of the halo to another. For moderate cross sections
$\sigma/m_{\text{DM}}\lesssim 10\rm cm^{2}/g$, the predictions of $\Lambda$CDM
for the large scale structure correlations and the growth of halos over time
carries over unchanged to $\Lambda$SIDM models. For $\sigma/m_{\text{DM}}\gg
10\rm cm^{2}/g$, the time scale for thermalization becomes shorter than the
dynamical time scale for galactic halos and the predictions for the large-
scale structure should be re-evaluated.

As a halo grows, self-interactions transfer heat from the outer (hotter) parts
to the inner (colder) parts of the halo (see the bottom panels of Fig. 3 for
example velocity dispersion profiles of CDM and SIDM halos). Over time, the
inner region of the halo becomes isothermal. This drive towards a constant
temperature is a complicated process but the end result can be understood by
focusing attention on the region of the halo where DM particles have had only
one or fewer interactions over the age of the halo. Simulations show that in
this region, the density profile of the halo is very similar to what would
have resulted had the interaction strength been dialed to zero. Thus, the
outer profile is well-described by the Navarro–Frenk–White (NFW) profile [23],
while SIDM in the inner regions is approximately isothermal. The radius at
which the behaviour transitions is approximately given by the radius where the
scattering rate multiplied by the age of the halo, $t_{\text{age}}$, is unity.
This radius is known as $r_{1}$, which is defined through:

$\displaystyle\rho_{s}\tilde{\rho}(r_{1}/r_{s})\langle v_{\rm rel}\rangle
t_{\rm age}\sigma/m_{\text{DM}}=1\;,$ $\displaystyle\implies 2.5\frac{V_{\rm
max}^{3}}{G_{\rm N}R_{\rm max}^{2}}\tilde{\rho}(r_{1}/r_{s})t_{\rm
age}\sigma/m=1$ (16)

where $\langle v_{\rm rel}\rangle$ is the average relative velocity between DM
particles and $\tilde{\rho}(x)=x^{-1}(1+x)^{-2}$ is the functional form of the
NFW profile $\rho_{\rm NFW}(r)=\rho_{s}\tilde{\rho}(r/r_{s})$. In the above
equation, $V_{\rm max}$ is the maximum circular velocity within the halo (i.e.
the maximum value of $\sqrt{G_{\rm{N}}M(<r)/r}$) and $R_{\rm max}$ is the
radius at which this maximum occurs. We have implicitly assumed a spherical
halo in making these calculations, but that is sufficient to understand the
behavior of $r_{1}$. The concentration-mass relation in standard $\Lambda$CDM
dictates that $R_{\rm max}\propto V_{\rm max}^{1.3-1.5}$ and hence (3.1)
implies that $r_{1}$ is an approximately fixed fraction of $r_{s}$ [106]. For
a cross section of $1{\rm cm}^{2}/{\rm g}$ and $t_{\rm age}=10\,{\rm Gyr}$,
$r_{1}\simeq r_{s}$, such that the region inside $r_{s}$ is changed by self-
interactions across a wide range of mass scales.

For cross sections much larger than $1\rm cm^{2}/g$ the core size does not
continue to grow. This is because the CDM halo is already isothermal around
$r=r_{s}$. This implies that there is a lower limit to the central density in
SIDM halos, which seems to be around 2-3 times $\rho_{s}$ [116, 117, 118,
103]. As $t_{\rm age}\times\sigma/m_{\text{DM}}$ is dialed up, all SIDM halos
(evolving in the field without baryons) will hit this floor in the density and
then enter a phase where the core starts to shrink, with the core density
starting to increase [99, 116, 117]. This process is analogous to the
gravothermal core collapse process that is already well-understood from
studies of globular clusters and that we discuss further in §3.2.

While the arguments above indicate where interactions are important in the
halo, they do not tell us about the density profile of the halo interior to
$r_{1}$. Since the density profile is isothermal, we can write

$\rho_{\rm SIDM}=\rho_{0}\exp\left(-\frac{\Phi_{\rm tot}(r)-\Phi_{\rm
tot}(0)}{v_{\rm{1D}}^{2}}\right)$ (17)

where $v_{\rm{1D}}$ is the one dimensional velocity dispersion of DM particles
inside $r_{1}$ and can be related to $\langle v_{\rm rel}\rangle$. The
gravitational potential $\Phi_{\rm tot}$ is due to all of the mass in the halo
including the baryons. Thus, the density profile of an SIDM halo is tied to
the gravitational potential of the baryons, unless the baryons in the halo are
dynamically unimportant (as in low surface brightness galaxies; LSBs). For
cross sections $\sigma/m_{\text{DM}}\gtrsim 1\rm cm^{2}/g$, the core size is
close to $r_{1}$ for the case where the baryons are dynamically unimportant.
As the stellar density increases at fixed halo mass (that is, the stellar
distribution becomes more compact), the core size shrinks in response to the
increased gravitational potential of the stars. In the limit where the stars
dominate the central potential, the core size is set entirely by the half-
light radius of the stars. Thus, at fixed halo mass, SIDM predicts that the
core size is between (roughly) $r_{\rm half}$ of the stars and $r_{1}$ [119].
The central density $\rho_{0}$ increases as the core size shrinks due to the
gravitational influence of the baryons.

These effects can be seen in Fig. 3, which shows the density profiles and
velocity dispersion profiles from $N$-body simulations of SIDM, in the
presence of a baryonic potential. The halo in question has a decreasing
central density with increasing cross section in the DM-only case. However,
with the addition of the gravitational potential due to a compact disc, the
central densities increase above that expected with CDM, and increase with
increasing cross section. The lower-panels show the velocity dispersion
profiles ($\sigma_{\text{v}}$ in the figure is $v_{\rm{1D}}$ in our notation),
which are isothermal in the inner regions as expected from the discussion
above.

Returning our attention to the isothermal density profile inside of $r_{1}$,
we note that it has two parameters that have yet to be determined: the central
density, $\rho_{0}$, and temperature / velocity dispersion, $v_{\text{1D}}$. A
very simple model where the isothermal mass profile is matched continuously
and smoothly (both density and mass are continuous) at $r_{1}$ to the NFW
profile has been shown to be an excellent fit to a wide range of SIDM
simulations, both with and without baryons [120]. One may consider many other
matching procedures, for example matching at a radius where one has more or
fewer interactions than 1. However, this simple procedure of matching at
$r_{1}$ does a good job. Investigations into obtaining a better fitting
function to the density profile would be useful, for example, matching at $N$
interactions where $N\simeq 1$ but is a mild function of halo mass or $V_{\rm
max}$ [50].

Figure 3: Figure from [121] showing the DM density (top) and velocity
dispersions (bottom) in the central regions of halos in $N$-body simulation
with self-interacting DM and a central disk. The left panel shows the density
and dispersion in the presence of a compact baryonic disk and the right panel
corresponds to an extended disk. In the absence of a baryonic potential the
main effect of introducing self-interactions is thermalization and the
formation of a central core. A baryonic disk in the centre can shorten the
period of core expansion and trigger core contraction.

###  3.2 Gravothermal collapse of SIDM halos

The transformation of the cusp at the halo center into a core due to elastic
self-interactions is a transitory phase that leads to a quasi-equilibrium
configuration once the core has achieved its maximum size, which is about the
size of the radius where the velocity dispersion profile peaks. Prior to this
stage, the transfer of energy due to elastic collisions occurs from the
outside in, since the velocity dispersion profile has a positive gradient in
the inner region (see lower panel of Fig. 3). There is a net flux of heat from
the regions close to the maximum of the velocity dispersion to the halo centre
(e.g. [122]). Once the core reaches its maximum size, however, subsequent
collisions cause a heat flux from the inside out since the velocity dispersion
profile has a negative slope in the outer regions. This condition triggers the
gravothermal collapse phase of the inner region of the SIDM halo.

Gravothermal collapse [123] is a well known process in globular clusters,
where the inner regions have a negative specific heat that is smaller in
magnitude than the positive specific heat of the outer regions when they are
constrained to remain within a certain volume. In the case of an SIDM halo
both the inner and outer regions have negative heat capacity, so that the flow
of energy from the centre outwards heats the centre and cools the outskirts.
It is this runaway process that is responsible for the collapse. In the case
of globular clusters, the collapse is prevented by the formation of binary
stars. In the case of an SIDM halo, since interactions are purely elastic, the
core contracts to form a cusp and ultimately collapses to form a black hole.
This phase can be followed using the gravothermal fluid approximation [124,
100, 102, 125, 126], which combines the equations of hydrostatic equilibrium,
the heat conduction law (Fourier’s law) and the first law of thermodynamics.
These equations contain three characteristic scales in the center of the
halo—a time scale given by the relaxation time, and two length scales, the
mean free path and the Jeans length. There is no concrete derivation for the
effective conductivity in the long mean free path (LMFP) regime, and a semi-
empirical relation must be used to interpolate between the LMFP regime and the
short mean free path (SMFP) regime where the fluid approximation can be
applied. This relation to interpolate the conductivity between the LMFP and
SMFP regimes introduces a constant that must be calibrated to SIDM N-body
simulations [100].

Recently, investigations of the gravothermal fluid model have been extended to
velocity-dependent cross sections [90, 94, 104]. In the LMFP regime, the
gravothermal equations admit an approximate universality that allow, with an
appropriate averaging, a velocity-dependent model to be mapped on a velocity-
independent model [90]. This seems to be reflected in the evolution of halo
profiles in idealized N-body simulations [94]. Although further work with
idealized and cosmological N-body simulations is required, there is promise
here of a great simplification of the SIDM model space in as far as
predictions for halo profiles of field galaxies is concerned.

The gravothermal fluid approximation provides a time scale for the onset of
the gravothermal collapse phase in terms of the relaxation time scale at the
characteristic radius of the halo $r_{s}$ [99, 100, 102]: $t_{\rm gc}/t_{\rm
0}\sim 400$, where $t_{\rm 0}=1/(a\rho_{0}v_{\rm 1D}\sigma/m_{\text{DM}})$ is
the interaction time (i.e. $1/R$, where $R$ the scattering rate defined in
equation 1) in the core with central density $\rho_{0}=2.4\rho_{s}$,
$v_{1D}=0.65V_{\rm max}$ and $a=4/\sqrt{\pi}$ [100, 90]. For the median
densities and velocities of galaxy-sized halos, the collapse time exceeds a
Hubble time as long as the cross section per unit mass is $\lesssim$
$\mathcal{O}(50\,{\rm cm^{2}/g})$. Note that there are no constraints that
prohibit such large cross sections for velocities of order
$100\leavevmode\nobreak\ \,\mathrm{km\,s^{-1}}$ or smaller. However, the
collapse time varies inversely with the cube of the concentration parameter
[118, 103], which implies that a halo with concentration about 0.25 dex
(roughly 2-$\sigma$) higher would collapse within a Hubble time even for a
cross section of about $10\,{\rm cm^{2}/g}$.

An additional mechanism to speed up core collapse is tidal truncation of
satellite galaxies. An investigation using the gravothermal equations suggests
that the temporal evolution can be quickened sufficiently such that present-
day subhalos could be in the core collapse phase (density increasing with
time) even for moderate values of $\sigma/m_{\text{DM}}\sim
10\leavevmode\nobreak\ \mathrm{cm}^{2}/\mathrm{g}$ [103, 127]. The effect of
tidal truncation on halo collapse is demonstrated in Fig. 4 for an initially
NFW halo. The disruption is modeled as an abrupt truncation, so that beyond
the truncation radius $r_{t}$ (taken to be the NFW scale radius $r_{s}$), the
density is $\rho(r)=\rho_{\mathrm{NFW}}(r_{t})\times(r_{t}/r)^{5}$. For an
isolated halo, the core is gradually formed and persists through today.
However, if the halo is initially truncated, the steepened outer slope creates
a temperature gradient within the halo that allows for more rapid heat
transfer and faster collapse.

Figure 4: _Left panel_ : Gravothermal evolution for the central density of a
halo with an initial NFW profile (solid) and with truncated NFW profiles, with
truncation radii at $r_{s}$ (dotted) and $3r_{s}$ (dashed). The evolution of
the central density is normalized to the NFW scale density $\rho_{s}$ and
shown as a function of time, normalized to an interaction time scale
$t_{0}\propto 1/(\rho_{s}V_{\rm max}\sigma/m_{\text{DM}})$. The top axis shows
the self-interaction cross section per mass needed for the associated
dimensionless time on the bottom axis to correspond to 13 Gyr of evolution.
From Ref. [103]._Right panel:_ Evolution of the central density of a
$10^{10.5}M_{\odot}$ halo, relative to the density at the start of the
simulation, as a function of concentration $c$ and isolation criterion, for a
constant cross section of $\sigma_{\mathrm{T}}=6\,$cm2/g. Magenta lines denote
halos evolved in isolation. Black lines show halos on a radial orbit within a
group-scale halo, but with evaporation by the host turned off in order to
highlight the effects of gravitational tides on the halo evolution. Turning
off evaporation is also the correct thing to do for SIDM models with
significant velocity dependence where the cross section falls below a few $\rm
cm^{2}/g$ at velocities of $200\ \rm km/s$. From right to left, the
concentrations are 45, 60, 75 and 90. Figure by Z. C. Zeng, adapted from
[127].

Core collapse will contribute to the diversity of halo profile shapes, as
demonstrated in $N$-body simulations [128, 51, 129, 130, 127]. Moreover, the
process of core collapse is faster for halos with higher concentrations, and
the formation of compact, dense cores increases the probability of the halo
surviving the tidal disruption event [129]. While further studies and
simulations are needed to test these predictions of SIDM against data, there
is an intriguing anticorrelation between the central densities of bright MW
dwarf spheroidal galaxies and their inferred orbital pericenter distances from
Gaia [131]555https://sci.esa.int/web/gaia, which could be a window into the
core collapse process.

Tidal heating and SIDM-driven evaporation may also delay or halt core
collapse, potentially further enhancing the diversity of halo central
densities [127]. On the other hand, the presence of a deep baryonic potential
can accelerate core collapse as highlighted in Ref. [132, 121] and this can
have a range of astrophysical consequences [133, 134]. One of the relatively
less-explored consequences is the formation of central black holes [135]. It
is worth noting that black holes are expected as the end state in the core of
a SIDM halo at late times with or without the presence of baryons or a central
seed black hole [136, 103, 137].

The core collapse phase is the longest phase in the evolution of a SIDM halo
and subhalo [116, 103, 90] relative to the total time for collapse. In models
with large cross sections, it is natural to expect that many of the halos and
subhalos will be in this phase of evolution. A full investigation of the
different aspects of core collapse promises to be a fruitful avenue for
further research.

###  3.3 Lack of Dynamical Friction in Cored Halos

An object, like a satellite galaxy, passing through a field of discrete
particles, like stars or DM particles, experiences a drag force from the
gravitational interaction between the object and field particles—dynamical
friction [138]. In the context of galaxy formation theory, the main role
dynamical friction plays is to cause satellite galaxy orbits to shrink,
causing the satellite to merge with the host [139, 140]. Three different
interpretations of dynamical friction lead to the same quantitative
description of satellite orbital decay: 1. Slow field particles are
accelerated by the fast object, leading to a negative acceleration of the fast
object [138, 141] 2\. Slow field particles are deflected behind the fast
object, creating a wake that pulls back on the fast object [142]. 3. The fast
object loses energy and angular momentum through a series of resonances with
field particles [141, 143]. For cuspy DM halos, the decay time is short (of
order a Hubble time or less) for satellite-to-host mass ratios greater than
1:100 [144].

However, satellites orbiting in _cored_ DM halos do not experience dynamical
friction. Instead, satellites or central galaxies trapped in the core
experience (nearly) undamped oscillatory motion about the center of the halo.
This unusual effect was first revealed to explain the existence of globular
clusters in the MW dwarf spheroidal galaxy Fornax, which should have sunk and
dissolved if Fornax inhabited a cuspy DM halo [145, 146]. Ref. [145] showed
that this effect could be explained in the context of Interpretation #3. The
gravitational potential of a cored halo is identical to a harmonic oscillator
in the core; there is only one natural frequency for particle orbits.
Therefore, the satellite does not sweep through a series of particle
resonances, like the cusped halo case; there is only one resonance. The
resonance model of dynamical friction thus breaks down. This effect can also
be explained in terms of Interpretations 1 & 2 of dynamical friction. The
authors of Ref. [147, 148] showed that there are few particles in the cores
that are slower than a satellite on a circular orbit, owing to the isothermal
nature of the particle velocity distribution.

Ref. [149] found that mergers between halos can trigger the oscillation of
galaxies in the cores of SIDM halos. The effect is especially noticeable on
cluster scales because of the potentially large core sizes, and especially for
the brightest cluster galaxies (BCGs). Figure 5 from [149] demonstrates this
effect, the BCGs of the individual clusters in a merging system do not sink to
the center but instead oscillate about it. Thus, the sloshing of galaxies in
halo cores is potentially a smoking gun of cored DM halos. So far, this effect
has only been studied in the context of equal-mass mergers, and in
dissipationless simulations. Further work is required to determine if the BCG
sloshing is a generic prediction of SIDM for less extreme mergers.

Figure 5: Figure from [149] showing the evolution of DM and galaxies in a
merger of halos in SIDM. The top panel shows the distribution in both the
halos while the bottom panel shows the distribution in one of the merging
halos. The BCGs remain offset and do not merge due to reduced dynamical
friction in the core.

###  3.4 Enhanced Gravitational Tidal Stripping of Cored Subhalos

One of the original motivations for SIDM was the removal of subhalos through
the process of evaporation [25], in order to solve what at the time was a
pressing “missing satellite problem" [19, 20]. Although neither the original
problem—the missing satellite problem is no longer a problem [59, 68, 69,
150], and in fact, there may be a “too many satellites" problem in the MW [69,
151]—nor the physical mechanism—evaporation—are primary drivers for
contemporary subhalo evolution considerations outside core-collapse scenarios
[75, 127], it is useful to consider evaporation as a launch point to a more
general and observationally relevant discussion of subhalo evolution.

What happens when a smaller halo falls into a larger one? The original Spergel
& Steinhard (2000) picture [25] focused on evaporation. Because the typical
orbital speed of a subhalo is much larger than the escape velocity of a
particle from the subhalo, scattering between subhalo and host-halo DM
particles usually leads to neither remaining bound to the subhalo. This need
not be true if scattering by large angles is strongly suppressed relative to
small-angle scattering, but even in this case a large number of interactions
can push subhalo particle above the subhalo escape velocity, leading to
‘cumulative evaporation’ [91, 152]. Thus, the subhalo loses mass from these
evaporative scatters, a process that is accelerated as the subhalo expands and
becomes less bound.

While a theoretically attractive paradigm, simulations show that this effect
changes the subhalo mass function negligibly unless the momentum transfer
cross section is larger than
$\sigma/m_{\text{DM}}>10\,\mathrm{cm^{2}\,g^{-1}}$ for isotropic scattering
[74, 75] (although see [153] for an opposing view). Such cross sections are
excluded to high significance on cluster scales in particular [154, 55, 149].
To explain why such a large cross section is required, we note that a
subhalo’s orbital velocity will be similar to the velocity of individual host-
halo DM particles, such that the rate at which subhalo DM particles scatter
with host-halo DM particles will be similar to the rate at which host-halo DM
particles (at the radius of the subhalo) scatter with one another. Significant
evaporation is therefore only expected for subhalos at radii $\lesssim r_{1}$
(see Eq. 3.1), where we expect SIDM also to affect the DM density profile of
the host halo. For example for the MW the scale of an allowed DM core is much
smaller than the orbital radius of the classical dwarf spheroidal galaxies,
which are not therefore expected to be significantly affected by evaporation
due to SIDM interactions directly.

Although the subhalo mass function is essentially unaffected by self-
interactions, the evolution of the satellite contained within the subhalos may
be significantly altered on account of gravity alone. For dark-matter-
dominated systems like dwarf galaxies, self interactions lead to the formation
of cores. Cored halos are less bound than cusped halos with the same mass.
Well outside the cored region of the halo, the density profile of the halo is
nearly identical to a CDM halo of the same mass, as is the gravitational
potential. Mass loss from the halo is hence similar in the CDM and SIDM cases
until the tidal radius lies near the core radius, at which point the
divergence in tidal evolution accelerates [75]. This occurs in the absence of
scattering, and is a purely gravitational effect [155].

What does this mean for the baryonic satellite at the center of the halo?
There are actually two effects relevant for subhalo evolution, on account of
the baryons lying within the core region. First, even if the tidal radius is
well outside the core region, DM particles on eccentric orbits whose
pericenters lie within the core region may be tidally stripped at apocenter.
This reduces the mass of the core region, and leads to an expansion of the
satellite and dark core as the system equilibrates. The effect is more severe
for cored than cusped halos, leading to a greater expansion of the satellite
half-light radius [75]. The expansion of the half-light radius was also
invoked in the evaporation case (removal of core DM particles by scattering
rather than tides) to expand elliptical galaxies in cluster environments,
which would lead to a (so far unobserved) dependence of the fundamental plane
of galaxies on environment [156]. Second, once the tidal radius does approach
the core radius, tidal stripping of both DM and stars proceeds more quickly
than in the cusped case [155, 75]. These effects are illustrated in Fig. 6,
which shows the evolution of the DM halo (black) and stellar component (red)
of a dwarf galaxy simulated in a cluster environment, as a function of the
slope of the inner dark-matter halo profile. After 98% of the halo mass has
been stripped, the remaining stellar mass is much more extended if it is
embedded in the cored halo than if it were in a cusped halo.

In both cases, tidal effects will be more severe in the presence of a massive
galaxy at the center of the host galaxy beyond the dark-matter-only case
considered by Ref. [75]. The presence of a disk dramatically alters the host
gravitational potential within a few scale lengths of the disk, leading to a
substantial drop in the tidal radius of the satellite [155]. Because baryonic
feedback may also lead to core formation in and enhanced tidal evolution of
dwarf satellites [157], distinguishing between SIDM and baryonic core
formation requires a careful assessment of core size as a function of subhalo
mass for each effect.

An additional implication is that the subhalo radial distribution in the host
may be much less centrally concentrated than in dark-matter-only simulations,
on account of the extra tidal field from the host galaxy [158, 159, 160].
Satellite disruption will be enhanced preferentially for small-pericenter
orbits. Although this population is relatively small compared to the overall
satellite population, it is important for the interpretation of the missing
satellites problem [69] and prospects for substructure lensing probes of DM
physics [161].

Finally, we note that the global properties of stellar halos—the moderately
phase-mixed debris from disrupted satellites—are independent of SIDM model
[75]. Most of the mass of the stellar halo comes from the few most massive
satellites to fall into the host because of the steepness of the stellar
mass–halo mass relation and the halo mass function [162]. These satellites
quickly sink to the center of the host by dynamical friction. The mass loss is
driven by the rapidly shrinking tidal radius, accompanying the rapidly
decaying orbit. The orbital decay is independent of SIDM model, and hence, the
stellar halos have bulk properties identical to the CDM case. See Ref. [75]
for an extended discussion of the conditions under which SIDM might affect
stellar halo observables.

In short, the relative evolution of SIDM and CDM satellites is driven by
gravity rather than non-gravitational scattering. Self-interactions matter in
the sense that they set the core size of halos, and hence the tidal radius.
Once the core is set, though, the evolution can be modeled as if halos are
made up of collisionless particles, unless the evaporation rate is very high
[127].

Figure 6: Mass loss and stellar extent in cored (left) and cuspy (right)
subhalos. The halos originally have masses of $3.6\times 10^{10}$ M⊙ (solid
black line), with a stellar component of $7.5\times 10^{7}$ M⊙ embedded within
(red), and orbit a cluster-mass halo. After 98% of the halo mass has been
stripped (dotted lines), the stellar component is larger in half-light radius
if embedded in a cored halo (left) rather than a cusped halo (right).
Reproduced from Ref. [163].

###  3.5 Halo shapes

Our discussion above has focused on the spherically-averaged density profile
of SIDM halos. The shape of the SIDM halo is another key discriminant between
SIDM and CDM halos. This is because the self-interactions that thermalize the
inner part of the halo also isotropize the DM particle orbits and lead to axis
ratios closer to unity (meaning more spherical) in the inner regions. These
two predictions are in contrast to the predictions of CDM, where the velocity
dispersion anisotropy can be substantial and the shapes of halos are
distinctly triaxial in the center.

Figure 7: This figure from [164] shows the median minor-to-major axis ratios
of different components of simulated galaxy clusters in the BAHAMAS-SIDM
cosmological simulations. SIDM1 and SIDM0.1 correspond to $1$ and
$0.1\,\rm{cm^{2}\,g^{-1}}$, respectively, while vdSIDM is an SIDM model with a
velocity-dependent cross section. The panels show the shapes of the DM (left),
gas (middle) and stars (right), with dashed lines corresponding to the DM-only
results. While SIDM leads to significantly more spherical DM distributions
than CDM, this is not obviously reflected in the shapes of the stars or gas.

An example of the effects of SIDM on the shapes of halos for galaxy clusters
is shown in Fig. 7. SIDM makes halos rounder, especially towards the centre.
Baryons also affect halo shapes, but the difference between CDM and
$1\,\rm{cm^{2}\,g^{-1}}$ is significantly larger than the change when going
from DM-only to including baryons. While SIDM makes the DM distribution
rounder, it does not (at least for the clusters simulated in [164])
significantly affect the shapes of either the gaseous or stellar
distributions.

It is important to note that for a given cross section and halo age, the
degree of impact of self interactions on the structural properties of halos is
not the same in each case. Simulations have shown that of order one scattering
per particle is enough to significantly affect the density profiles (see
above). However, the central (typically cored) regions do not get fully
spherical (e.g. [165, 166]) and they retain a degree of orbital anisotropy
[166]. An interesting fact in this regard is that SIDM halos with a dominant
disk in the center will have a distinctly non-spherical DM halo [119] because
the halo shape will be set by the potential of the stellar disk.

###  3.6 Drag Force and Mass Loss due to Self-Interactions

In the presence of DM self-interactions, DM halos moving through a background
DM density experience both evaporation, as discussed in §3.4, and deceleration
[156, 167, 168]:

$\frac{\dot{M}}{M}=-R_{\mathrm{e}}\,,\qquad\frac{\dot{v}}{v}=-R_{\mathrm{d}}\;,$
(18)

where $R_{\mathrm{e}}$ and $R_{\mathrm{d}}$ denote the evaporation and
deceleration rates, respectively. In terms of the total scattering rate $R$
defined in eq. (1), these rates can be written as
$R_{\mathrm{e,d}}=\chi_{\mathrm{e,d}}\,R$, where $\chi_{\mathrm{e}}$ and
$\chi_{\mathrm{d}}$ denote the fraction of scattering events leading to
evaporation and the average fraction of momentum lost with each scattering
event, respectively. One then obtains the simple expressions

$\displaystyle\frac{\dot{M}}{M}=-\chi_{\mathrm{e}}\frac{\sigma}{m_{\text{DM}}}\rho\,v\,,\qquad
F_{\text{drag}}=-\chi_{\mathrm{d}}\,\sigma\,\rho\,v^{2}\;,$ (19)

where in the second expression we have defined $F_{\text{drag}}\equiv
m_{\text{DM}}\,\dot{v}$ in order to recover the well-known form of a drag
force.

The evaporation and deceleration fractions $\chi_{\mathrm{e}}$ and
$\chi_{\mathrm{d}}$ depend on the ratio of the escape velocity to the velocity
of the halo relative to the background density and on the differential cross
section $d\sigma/d\theta$ [152]. Roughly speaking, the evaporation fraction is
reduced for a tightly bound DM halo (high escape speed), while the
deceleration rate is enhanced (and vice versa) [91]. Conversely, the
evaporation fraction is enhanced for larger relative velocities, while the
deceleration fraction is suppressed [149]. The relative importance of drag and
evaporation depends sensitively on the angular and velocity dependence of the
self-interaction cross section. For example, frequent interactions with low
momentum transfer can lead to a a sizeable drag force on infalling objects
while preventing complete evaporation [152].

Furthermore, it was pointed out in Ref. [91] that the deceleration of a DM
halo due to DM self-interactions generically leads to a “heating” of the halo,
as the directed motion of the entire halo is converted into the random motion
of individual DM particles. This heating will cause some of the DM particles
to obtain a kinetic energy exceeding their binding energy and escape from the
DM halo, thus contributing to the evaporation rate. Writing the rate of this
_cumulative evaporation_ as

$\frac{\dot{E}}{E}\approx\frac{\mathrm{d}v^{2}/\mathrm{d}t}{v^{2}}=-R_{\mathrm{c}}\;,$
(20)

one typically finds $R_{\mathrm{c}}\approx R_{\mathrm{d}}$ [91]. In other
words, it is generally impossible to have deceleration without at least some
amount of evaporation.

While evaporation does not lead to a complete disruption of the subhalo, the
main consequence of evaporation due to self-interactions is a decrease of the
mass to light ($M/L$) ratios of halos in interacting systems [169]. This
phenomenon is relevant both in major and minor mergers. In minor mergers,
which are more common, the subhalo mass-to-light ratio is expected to be
smaller than the one of typical clusters after one pericenter passage. While
DM particles are scattered away into the ambient medium, the stars in the
galaxies do not experience scattering, reducing the $M/L$ ratio. Ref. [149]
points out however that evaporation of DM particles also can also lead to a
loss of stars due to the reduced gravitational binding energy. At least for
equal-mass mergers, this can potentially compensate the decrease in $M/L$,
making it more difficult to obtain robust predictions. One related consequence
of evaporation of subhalo particles is the thickening or puffing up of
galaxies that they harbor, due to the diminished gravitational binding energy
from the subhalo. However, it may be difficult to disentangle the thickenings
of disks and puffing up of ellipticals due to SIDM from that due to mass loss
from tidal forces in the cluster itself.

Finally, evaporation can also lead to a deformation of the shape of the DM
halo. The reason is that DM particles escaping from a DM halo typically do so
in the direction opposite to its velocity (relative to the background DM
density). Thus, SIDM may affect the formation of DM trails, which are much
more asymmetric than the ones induced by tidal forces [91]. This effect can be
searched for by measuring the skewness of DM distributions [170, 171].

The drag force experienced by SIDM halos has a number of observational
consequences. The main one is an offset between the centroids of stars and DM
in merging or interacting systems, most obvious in the case of merging
clusters where the DM would end up between the dissipational gas and
effectively interactionless galaxies [172]. Merging clusters are rare,
however, and a much more common situation is the infall of galaxies into more
massive halos, which is the process responsible for growing larger structures
in hierarchical structure formation.

The challenge with satellite galaxies is that the much smaller lensing signal
precludes statistically significant comparison of the stellar and lensing mass
centroids. A possible way forward is by statistically combining a large number
of minor and major mergers [112, 173], which can potentially lead to an
improved sensitivity (see section 4.4). Even relatively small star–halo
offsets may however be detectable indirectly by means of their effects on the
structure of stellar disks, and in particular the generation of a
characteristic U-shaped warping. We discuss this signal in more detail in Sec.
4.9.

In isotropic interactions the effect of drag, however, is subdominant.
Simulations of mergers of massive clusters [149] have shown that offsets
greater than $50\,\mathrm{kpc}$ are unlikely in merging systems where DM
interacts isotropically, these offsets are moreover short lived and disappear
in less than a dynamical time after pericenter passage. For observed offsets
as large as $100\,\mathrm{kpc}$, cross sections much higher than
$1\,\mathrm{cm^{2}\,g^{-1}}$ are needed which will cause the infalling subhalo
to evaporate in the first place. Interestingly they highlight an alternative
and more robust signature, whereby DM self-interactions alter the evolution of
cluster merger, triggering a core sloshing of galaxies. As clusters merge and
the system relaxes, the central BCG or the centroid of the galaxy distribution
is offset from the DM distribution and oscillates around the barycenter, a
phenomenon that does not arise in the collisionless DM scenario. The
oscillation of the BCG, that is mainly produced due to reduced dynamical
friction from the cored merger remnant (see section 3.3), can be as large as a
few $100\,\mathrm{kpc}$ for cross sections of
$\sigma_{\mathrm{T}}/m_{\mathrm{DM}}\leq 1\,\mathrm{cm^{2}\,g^{-1}}$.

While the main focus of observational tests in merging systems remains on the
offsets between stars and DM and on $M/L$ ratios, it is important to note that
DM self-interactions can affect the orbits of subhalos within their hosts,
thereby altering the distribution of galaxies in halos. For example subhalos
on radial orbits that have small pericenters, pass through denser regions of
the host, undergoing more interactions than subhalos that are on tangential
orbits. Depending on the nature of the cross section, this may cause radially-
orbiting objects to lose more energy due to drag forces than tangential ones.
This leads to objects on radial orbits being trapped near the centre changing
the velocity distribution of subhalos within their hosts in comparison to CDM
([95]). The splashback radius, which is the boundary of the multistreaming
region of a halo, is the apocenter of recently accreted objects [174, 175].
The slope of the density profile of a halo rapidly falls off in a narrow
localized region around this radius. Drag forces on subhalos due to self
interactions lead to a loss of momentum that may cause the splashback radius
to shift to smaller cluster-centric distances. The splashback radius is
observed as a minimum in the slope of the number density profile of galaxies
[176, 177, 178], the location of this feature can change depending on the
nature of the cross section of interactions and the relative importance of
drag forces.

###  3.7 Large-scale structure of SIDM models

For canonical SIDM models where the _only_ relevant dark sector interaction
following DM chemical decoupling is DM self-scattering, the success of the
standard $\Lambda$CDM paradigm on large cosmological scales is retained. For
these models, mass and momentum conservation ensures that the evolution of
large-scale density perturbations that are still in the linear regime are
unaffected by the presence of DM self-interactions [179]. The linear SIDM
power spectrum is thus identical to its CDM counterpart on large cosmological
scales for the canonical models.

In realistic particle physics implementations of SIDM models, however, new DM
interaction beyond self-scattering are often present, which could potentially
modify the evolution of DM density fluctuations. This occurs for instance in
models where DM self-interactions are mediated by a relativistic or massless
particle. For these, a large scattering rate between DM and the thermal bath
of relativistic mediator particles at early times can efficiently erase DM
density fluctuations in a process analogous to Silk damping [180]. A similar
phenomenon occurs in theories where the force mediator is massive but is
itself coupled to relativistic particles at early times.666From a model-
building perspective, such a coupling might be needed to avoid having the
mediator overclose the Universe (see e.g. Ref. [5]). In both cases, the key
ingredient responsible for the modified evolution of DM density fluctuations
is the presence of relativistic particles coupling (directly or indirectly) to
DM. Examples of SIDM models where such coupling occurs include:

  * DM interacting with a massless dark photon that mixes with the SM photon [181, 182, 183, 184].

  * DM interacting with a massless dark photon that does not mix with the SM photon [8, 185, 168, 186, 93, 187].

  * Atomic DM [188, 189, 190, 191, 192, 193, 194, 195, 84, 196, 197, 86, 198].

  * DM interacting with neutrinos [199, 200, 201, 202, 203, 204, 205, 206].

  * DM interacting with photons [199, 200, 201, 207, 208, 209, 210, 211, 205, 212].

  * DM interacting with the baryon-photon bath [213, 208, 214, 215, 216, 217, 218, 219, 220, 221, 57, 222, 223, 224, 225, 226]

  * DM interacting with a non-abelian gauge boson [85, 227, 228, 229, 230, 231, 232].

  * DM interacting with sterile neutrinos or other types of “dark” radiation [233, 234, 235, 236, 237, 5].

The evolution of DM fluctuations in the presence of significant interactions
with a relativistic species at early times mimics that of the SM baryons
before the epoch of recombination [238, 239, 240]. In the latter case, the
photons and baryons form a tightly-coupled fluid with a large sound speed,
allowing for the propagation of acoustic waves (known as the baryon acoustic
oscillation, BAO) to large cosmological distances. However, the finite photon
mean free path $\lambda_{\gamma}$ within the plasma weakens the pressure
support of the acoustic oscillations, leading to significant damping of their
amplitude on scales smaller than $\lambda_{\gamma}$. The net result of this is
the presence of damped baryonic oscillations in the matter power spectrum at
late times, which have been observed in galaxy redshift surveys [241, 242].

A similar situation occurs for DM interacting with a relativistic species at
early times, except that the amplitude of the “dark” acoustic oscillation
(DAO) in the matter power spectrum can be significantly larger than the BAO
due to the higher abundance of DM as compared to baryons. For the same reason,
the damping caused by the diffusion of the light species out of DM
overdensities can have a dramatic impact on structure formation, essentially
erasing all density perturbations on scales below the relativistic species’
mean free path. Of course, the exact scale at which this occurs as well as the
shape of the damping envelope of the matter power spectrum depend on the type
and strength of the coupling between DM and the relativistic species (e.g.,
see Refs. [8, 192]). More precisely, the shapes of the oscillations and
damping envelope strongly depend on the width of the DM drag visibility
function (defined in analogy with the baryon drag visibility function [243]),
which itself reflects the dependence of the scattering amplitude on the
momentum of the incoming relativistic species. This is illustrated in Fig. 8
where we observe that both the oscillation wavelength and the width of the
damping envelope decrease as the index $n$ describing the momentum dependence
of the squared matrix element is increased (here,
$|\mathcal{M}|^{2}\propto(p/m_{\rm DM})^{n-2}$, where $p$ is momentum of the
incoming relativistic species that is scattering with DM).

Figure 8: _Left panel_ : Interacting DM transfer function ($T(k)\equiv(P_{\rm
SIDM}(k)/P_{\rm CDM}(k))^{1/2}$) for two different types of interaction
between DM and a relativistic species. Each model is parametrized with an
index $n$, which corresponds to the momentum dependence of the squared matrix
element describing the scattering process, $|\mathcal{M}|^{2}\propto(p/m_{\rm
DM})^{n-2}$, where $p$ is the momentum of the incoming relativistic particle.
Here, the values of the coupling strength are chosen such that all models
shown have the same DM drag epoch. As shown in the legend, the $n=2$ model
corresponds to a Compton-like interaction between DM and the dark radiation,
while the $n=4$ case corresponds to DM interacting with dark radiation via a
massive (vector or scalar) mediator. For comparison, we also display the
transfer function for a warm DM model with thermal mass $m_{\rm WDM}=3.5$ keV.
_Right panel_ : DM drag visibility function plotted against conformal time
$\tau$ for the same models as the left panel. The DM drag visibility function
is essentially the probability distribution function for the time at which a
DM particle last scatters off the relativistic species. Figure adapted from
Ref. [179].

In SIDM models where DM remains coupled to the relativistic species until
close to the epoch of hydrogen recombination, the suppressed DM density
fluctuations affect the CMB temperature and polarization power spectra in
characteristic ways [244]. Generally speaking, DM must kinetically decouple
from the relativistic species before the epoch of matter-radiation equality to
avoid tensions with CMB measurements, unless only a small DM fraction couples
to the radiation [244, 245, 232]. More stringent constraints on the presence
of DAO and damping in the matter power spectrum can be put using Local Group
satellite galaxy counts and the Lyman-$\alpha$ forest [246, 211, 247, 205,
248, 5], and the reionization history of the Universe [249, 250]. These bounds
are primarily driven by the suppression in the DM power spectrum and the
correlated suppressed abundance of low-mass halos in the late Universe, which
lead to a deficit of satellite galaxies around the MW and could also delay the
formation of the first stars.

DM halos with masses below the scale at which the SIDM halo mass function
starts to deviate from its CDM counterpart, due to the early coupling between
DM and the relativistic species, experience a delay in their formation. Since
the typical inner density of DM halos is proportional to the mean cosmological
DM density at the epoch at which they form, these halos tend to be less dense
then their CDM counterparts. Therefore, SIDM models with significant damping
in their matter power spectrum on scales corresponding to dwarf galaxies yield
very diffuse DM halos at those mass scales that are difficult to reconcile
with observations [247]. An important corollary is that the constraints on the
SIDM cross section and on the allowed amount of suppression in the matter
power spectrum are not independent from each other [233, 86, 251, 247].
Ideally, for DM models where both self-interaction and damping are present,
these two quantities should be constrained simultaneously.

###  3.8 Are these physical effects unique?

To interpret the observations of any of the phenomena described in the
previous sections as evidence for self-interactions of DM, it is important to
consider whether they could arise by other means, for example baryonic
physics, modified gravity or alternative theories of DM. In this subsection,
we summarise the main alternative types of physics that could be
observationally degenerate with certain SIDM phenomenology, along with the
best methods to distinguish them from SIDM.

####  3.8.1 Degeneracies with baryonic physics

As we have seen, a key prediction of SIDM is the formation of DM cores in the
central regions of halos, in contrast to cuspy profiles seen in pure CDM
simulations. However, there are several reasons to expect the CDM prediction
of cusps to be altered by baryonic physics in the process of galaxy formation.
The simplest consequence of the dissipative formation of a cold baryonic disk
at halo centers is that DM is drawn inward gravitationally in a process known
as “adiabatic contraction” [252, 253, 254]. This process further steepens the
central density profile. However, stars, supernovae, and energetic processes
originating in the galactic center may subsequently inject energy and momentum
into the halo by means of feedback processes, causing it to expand (e.g. [255,
256]). Due to the highly non-linear nature of this feedback effect (and unlike
adiabatic contraction), its magnitude cannot be calculated analytically,
necessitating the use of hydrodynamical simulations which can track the
baryonic processes. However, even high-resolution simulations cannot cover the
full range of scales that are important for the processes leading to
feedback—from the collapse of giant molecular clouds at
$\mathcal{O}(\text{pc})$ to the baryon cycle of the circumgalactic medium at
$\mathcal{O}(\text{Mpc})$. The smallest relevant scales are therefore below
the grid size of all but the highest resolution simulations, introducing a
dependence on uncertain, empirically-calibrated “sub-grid models”. Different
simulations use different sub-grid models, such as the prescriptions for star
formation, the coupling of stellar radiation and winds to DM, supernova rates,
cooling, and gas re-accretion, and therefore end up with different predictions
for observables.

While some simulations predict adiabatic contraction to remain the dominant
effect on the halo density profile (e.g. [257, 258]), others produce cores. Of
notable interest are high-resolution galaxy simulations such as FIRE [259,
260], NIHAO [261] and MaGICC [262, 263], which have been designed to
investigate the effects of multiple interconnected feedback mechanisms in a
cosmological setting. FIRE is a suite of zoom simulations of galaxies ranging
from ultra-faint to MW mass, with sub-pc and $\sim 100M_{\odot}$ spatial and
mass resolution respectively. Significantly, it is found in these simulations
that both radiative and supernova feedback are necessary (and sufficient) for
reproducing the low galaxy formation efficiency expressed in the relations
between stellar and halo mass, and between stellar mass and specific star
formation rate. A similar conclusion is reached in the NIHAO suite of 100
galaxies using a comparable feedback scheme, as well as in the MaGICC
simulation which incorporates an early source of thermal stellar feedback.

Besides regulating star formation, these feedback schemes impart energy to the
DM, causing it to move out of the central regions inhabited by the galaxy. The
result is a flattening of the initial NFW cusp into a core. In particular, a
clear dependence is found between the central slope of the DM profile and
$M_{*}/M_{\text{halo}}$, which determines the amount of energy injected into
the halo relative to the depth of its potential well (e.g. [264]). Fig. 9
shows such an example from the MaGICC simulation. Core formation through
feedback is efficient at $10^{-3}\lesssim\log(M_{*}/M_{\text{halo}})\lesssim
10^{-2}$, corresponding to $10^{6}\lesssim M_{*}/M_{\odot}\lesssim 10^{10}$,
but ineffective on either side (see Fig. 9), and these results were used to
create parameterized halo profiles that incorporate baryonic effects (e.g.
[265]). A detailed Bayesian comparison of the fits with the parameterized halo
profile shows that these profiles do as well as the analytic SIDM profiles at
fitting the rotation curves in the SPARC database [266]. It would be
interesting to see how closely the density profile parameters fit to the
observations, track the same parameters for the galaxies in the MaGICC
simulations used to created the parameterized halo profile.

At low halo masses, the dominant source of feedback energy is from supernovae.
However, both star formation and supernova rates fall towards lower halo mass
as galaxy formation efficiency declines, so that below a certain mass
($M_{\text{vir}}\sim 10^{8}-10^{9}M_{\odot}$), it is unlikely that sufficient
supernovae energy is released over the entire formation history to flatten the
cusp [267]. Even if core formation were possible for the dwarf spheroidal
galaxies currently at the sensitivity limit of surveys, it would be very
unlikely for lower-$M_{*}$ ultra-faint galaxies that will be detected by next-
generation surveys such as the Legacy Survey of Space and Time (LSST; [268,
269])777http://www.lsst.org on the Vera C. Rubin Observatory. These would be
predicted unambiguously to have cusps even when accounting for the effect of
baryons (but see [270]), so that observations of cores would provide strong
evidence for new DM physics.

Figure 9: Slope of the halo density profile between 1 and 2% of
$R_{\text{vir}}$ vs $M_{*}/M_{\text{halo}}$ for 31 galaxies from the MAGICC
simulation. Reproduced from Ref. [264].

Hydrodynamic simulations in $\Lambda$CDM have traditionally had difficulty in
accounting for the diversity in the rotation curves of galaxies with a range
of surface brightness at fixed halo mass: the strong feedback necessary to
push DM out of the central regions also moves the stars, which are similarly
collisionless, out. This makes it difficult to form high surface brightness
galaxies with steeply rising rotation curves [271, 272, 273] (but see [274,
275] and Sec. 4.7). This is not a known issue for SIDM, however, because the
stars are not directly affected by the processes that set the DM profile. In
the SIDM model, thermalization leads to more compact DM cores in more luminous
galaxies, which naturally explains the flat outer rotation curves of high
surface brightness galaxies as well as the variation in halo profiles with
galaxy size [50, 48, 128]. The SIDM halo profile remains as steep as the
$\Lambda$CDM one in galaxies where the gravitational potential is dominated by
stars and similar to $V_{\rm max}^{2}$ [119, 132, 121]. In cases where the
stellar disk is not dynamically important, the SIDM predictions for the core
size are tied closely to the NFW scale radius and depend sensitively on the
concentration parameter [48].

Since high and low concentration halos behave differently for SIDM (the former
are subject to gravothermal collapse and the latter to core formation), a
moderately large interaction cross section increases the range of inner DM
densities in satellite galaxies, bringing them more in accord with the ultra-
faint dwarf spheroidals around the MW [276]. SIDM is also expected to affect
the orbital trajectories of the satellites within their host halos [277],
which may provide an additional means of differentiating the model from one
with CDM and baryons. In summary, mass profiles of ultra-faint dwarfs combined
with detailed kinematics of low and high surface brightness galaxies show
promise in being able to distinguish SIDM effects from those due to strong
feedback in CDM halos. In addition, age and metallicity gradients may also
prove useful as discriminators [278].

It is worth noting that it is possible for the effects of stellar feedback on
star formation and halo density profiles to be different in SIDM than CDM,
since energy exchange causes significant fluctuations in the inner potentials
of halos at early times. However, hydrodynamical SIDM simulations indicate
that the consequences of SIDM are more robust to the inclusion of baryonic
physics than those of CDM: in both $M_{*}=10^{8}M_{\odot}$ [279, 280] and
$10^{7}<M_{*}/M_{\odot}<10^{8}$ [281] galaxies, inner DM density slopes are
little affected by baryons.

####  3.8.2 Degeneracies with novel dark matter or gravitational physics

Besides baryonic effects, the redistribution of dynamical mass in SIDM could
be mimicked by other types of DM interaction, or by modified gravity. By
giving the DM particle a thermal energy, warm DM models reduce the clustering
on small scales as well as lower the amplitude of the mass function at the
dwarf scale by inhibiting the collapse of small halos [282]. Fuzzy and Bose-
Einstein condensate (e.g. axion) DM, and other scalar field models, can also
effectively lead to core formation, and would, like warm DM, reduce the amount
of substructure [283, 284]. However, the particle mass required to form cores
of the correct density is at odds with the mass required to have enough
substructure in many fuzzy DM models [285]. Interactions between DM and
baryons could also significantly alter halo density profiles [286, 287]. The
suppression of small-scale power these models imply is likely the best means
to break their degeneracies with SIDM, which in the canonical case would not
alter the halo mass function (though see §3.7). Small-scale power may be
effectively probed using the Lyman-$\alpha$ forest, which has been used to
constrain both the temperature and fuzziness of DM [288, 289, 290]. Fuzzy DM
also generates unique signatures such as interference patterns which heat
stellar orbits and alter stellar streams and lensing [291, 292, 293], and
presents more complex dynamical friction behaviour than CDM [294, 295].

Because the inference of the halo density profile is purely dynamical, an
alternative to altering the clustering properties of DM is to modify the law
of gravity. A paradigmatic example is Modified Newtonian Dynamics (MOND) [296,
297, 298], which alters either Newton’s second law or the gravitational
inverse square law to fit galaxy kinematics from the baryon mass alone. The
model postulates a critical acceleration scale, $a_{0}\simeq 1.2\times
10^{-10}\>\text{m/s}^{2}$: above this the dynamics is Newtonian, while much
below it $F\propto a^{2}$, leading to asymptotically flat rotation curves.
$a_{0}$ corresponds to a critical baryon surface density $\Sigma_{0}\equiv
a_{0}/G\simeq 1\>\text{kg/m}^{2}$ above which the “phantom” DM mass
(reconstructed by a Newtonian analysis) falls off. MOND has been used to fit
rotation curves and velocity dispersion profiles of galaxies over a wide range
of masses, and a number of relativistic parent theories have been proposed
(e.g. [299]; see [300, 301] for reviews). In terms of rotation curves, the
relation between the total radial acceleration and the radial acceleration due
to baryons is not a single curve in SIDM models but the scatter around the
average is small [50]. Thus, SIDM looks like MOND in an averaged sense and it
is not clear if any further insight could be teased out of rotation curves.
These theories could be distinguished if they predict deviations on large
scales, including the CMB, matter power spectrum and properties of clusters.

Another class of modified gravity theories that can impact mass profiles of
galaxies are the scalar–tensor theories (e.g. [302]), which supplements the
metric tensor with one or more scalar fields. These fields induce a new
interaction (“fifth force”), which boosts rotation and dispersion velocities.
In the well-studied $f(R)$ model [303, 304] (which is equivalent to a scalar-
tensor theory [305]), the scalar field effectively increases Newton’s constant
$G$—and hence the squared velocity for given enclosed mass—by a factor of
$1/3$ on scales below the Compton wavelength of the scalar field. Most viable
theories including $f(R)$ implement a “screening mechanism” to hide the scalar
field in high density regions such as the MW where observational constraints
are strong (e.g. [306, 307, 308, 309]. This still allows a fifth force in
dwarf galaxies of much lower density, however, and may make an underlying
cuspy matter distribution appear core-like [310, 311] or otherwise alter the
inferred dynamical density profile. Other modified gravity models such as
symmetrons may have important effects on galaxy dynamics too [312, 313].

Finally, we note that modified gravity may mimic or confuse other SIDM signals
too. In Sec. 4.9 we describe how SIDM can lead to warps, asymmetries and
thickening of galaxy disks: this may also occur in the presence of a
partially-screened fifth force. In particular, under chameleon, symmetron or
dilaton screening the DM (and gas) in a low-mass galaxy may be unscreened,
while the stars are sufficiently massive to screen themselves [314]. This
induces an offset between the stellar and halo centres of mass in the
direction of the external fifth-force field, generating a potential gradient
that warps the disk [315, 316, 317]. This can also cause asymmetries in
rotation curves [318]. However, the warp is in roughly the opposite direction
to that predicted by SIDM and also has a different dependence on gravitational
environment, enabling the two types of physics to be disentangled. These
effects are described in detail in the modified gravity arm of the Novel
Probes Project [319].

##  4 Observations of SIDM phenomenology

An impressive diversity of observations have been used to test for SIDM. These
include observations of satellite and dwarf galaxies, spirals and massive
ellipticals and clusters using techniques ranging from optical spectroscopy
and imaging (including strongly lensed images) to X-ray observations.

This section covers observations that have provided useful constraints and are
expected to deliver improvements in the coming years. §4.1 is focused on
strong lensing, starting with the ‘classic’ elliptical galaxy lenses and
galaxy clusters, and moving on to the subtler effects due to subhalos within
the lens galaxy and along the line of sight. §4.2 discusses stellar streams of
the MW. §4.3 discusses the shape of galaxies’ and clusters’ mass profiles as
probed by X-ray and weak lensing measurements. Cluster mergers, including the
Bullet cluster, are discussed in §4.4, while §4.5 describes the kinematics and
locations of the brightest galaxies in clusters (BCGs). Dwarf galaxies in the
Local Group are discussed in §4.6 and the rotation curves of spirals in §4.7.
Ultra-diffuse galaxies and warped disk galaxies are discussed in §4.8 and §4.9
respectively.

###  4.1 Strong Gravitational Lensing in Clusters, Groups and Large
Ellipticals

In a strong gravitational lens, a background source is multiply imaged and
magnified due to the distorting effects of a foreground deflector such as a
massive galaxy or cluster of galaxies. The positions and magnifications of the
images of the background source depend on the first and second derivatives of
the projected Newtonian gravitational potential, which in turn is related to
the mass distribution of the deflector projected along the line of sight. Here
we discuss two ways in which strong gravitational lensing can potentially test
SIDM. First we discuss strong gravitational lensing constraints on the _total_
mass distribution of massive elliptical galaxies as well as galaxy groups and
clusters. In the second subsection we discuss how strong gravitational lensing
is sensitive to the mass distribution of DM _substructure_.

####  4.1.1 The total deflector mass distribution

_Galaxy-scale lenses:_ Several dedicated surveys have focused particularly on
using strong gravitational lensing to study early-type galaxies, including the
Sloan Lensing Advanced Camera for Surveys (SLACS, [320]) survey, the Strong
Lensing Legacy Survey (SL2S, [321]), and the Boss Emission Line Survey (BELLS,
[322]). In the simplest approach, the slope of the mass distribution of a
gravitational lens can be measured using two quantities: the Einstein radius,
within which strong lensing precisely constrains the projected enclosed mass,
and the stellar velocity dispersion, which provides a three dimensional mass
measure within typically half of the effective radius of the deflector [323,
324]. Alternatively, with deep data and detailed modeling, the radial
distortion of lensed images provides a direct measurement of the local slope
of the mass distribution [325, 326, 327, 328]. A robust result arising from
these studies is that the _total_ density profile, including DM and baryons,
is well approximated by a power-law with $\rho_{\rm tot}\propto r^{\gamma}$
with $\gamma\sim-2$, and small intrinsic scatter of $\simeq 0.16$ [329, 330,
331].

Separating the baryonic and DM density profiles is challenging because of the
significant contribution of stars to the mass budget in the regime probed by
strong lensing. Uncertainties in the stellar mass-to-light ratio, arising
especially from the stellar initial mass function (IMF) in massive elliptical
galaxies, preclude a direct subtraction of the stellar mass. The combination
of the Einstein radius and velocity dispersion constrains the normalization of
the stellar and DM density profiles _if_ the latter’s slope is assumed [332,
333]. Measuring both the normalization and slope requires additional
information, which can come from several sources: (1) the relationship between
stellar and halo mass [334], (2) ensemble constraints derived from samples of
strong lenses assumed to have homologous structures [335, 336, 333], (3) the
local density slope as measured from the detailed structure of the multiple
images [331], (4) in rare cases, multiple Einstein rings [337], and (5)
spatially resolved stellar velocity dispersions and spectroscopic diagnostics
of the IMF [338, 339]. These studies have reported inner DM density slopes
ranging from $\beta_{\rm DM}\approx 0.5-2$, where $\rho_{\rm DM}\propto
r^{-\beta_{\rm DM}}$, which emphasizes the difficulty of the measurements and
the possibility of large galaxy-to-galaxy variations. Further progress in
understanding the DM profiles in these galaxies will require more precise
constraints on their stellar mass distributions, which may come from
improvements in detailed stellar population models [340] or microlensing of
quasars [341, 342].

_Group- and cluster-scale lenses:_ The aforementioned samples of strong lenses
are primarily galaxy-scale lenses residing in $\sim 10^{13}M_{\odot}$ halos
[343]. The DM profile is more readily constrained in more massive halos due to
their lower baryon fractions and richer observational constraints. Massive
galaxy clusters in $\sim 10^{15}M_{\odot}$ halos often present several
multiple image systems with Einstein radii $\sim 50$ kpc. Spatially resolved
stellar velocity dispersions within the BCG probe the mass distribution within
$\sim 20$ kpc of the cluster center, while weak lensing and X-ray emission
probe the mass distribution on $\sim$Mpc scales [344, 345, 346, 52, 53].
Although there is significant covariance in an individual cluster between the
inner DM profile and the stellar $M_{*}/L$ within the BCG, this can be
resolved by combining constraints from an ensemble of clusters under the
assumption that their BCGs share a common IMF that is not known _a priori_. An
analysis of 7 clusters in this way [53] showed that the DM density profile is
shallower than NFW within a typical radius of $\sim 30$ kpc and has a mean
inner slope of $\beta_{\rm DM}=0.50\pm 17$ (c.f. $\beta_{\rm DM}=1$ for NFW,
see Fig. 10). At smaller mass scales, new samples of group-scale lenses have
been identified in panoramic imaging surveys [347, 348], and these can be used
to investigate the DM distribution in intermediate-mass halos of $\sim
10^{14}M_{\odot}$. These group-scale systems appear to be consistent with NFW-
like profiles [54], suggesting that shallow DM profiles may be confined to
massive clusters. Ref. [331] also found tentative evidence that the inner DM
density slope may not be universal, but instead appears to vary with the halo
mass or environment.

A systematic variation in the DM density slope with halo mass might reflect
trends in the formation histories of the galaxies they host. BCGs are
primarily assembled through dissipationless mergers that might evacuate DM
from the inner halo via dynamical processes [349, 350], whereas in galaxy or
group-scale halos, dissipational processes are probably more important, and
these might instead compress the inner halo [351, 352]. The net effect of
these baryonic processes is not well understood, especially in high-mass
halos. Within the context of SIDM, the gravitational coupling of DM to the
baryon distribution could produce a trend of the DM density slope with halo
mass [119, 55, 48, 50]; this would arise from the underlying variation of the
central baryon fraction, not from baryon-DM interactions during galaxy
formation.

Figure 10: Marginalized posterior probability densities for the DM inner slope
from a sample of seven massive ($M_{200}=0.4-2\times 10^{15}$ M⊙) relaxed
galaxy clusters at $z=0.19-0.31$ (from Ref. [53]). The thick black line shows
the joint constraint from all clusters, while the thin lines show the
individual posteriors from each cluster. The vertical dotted line shows the
standard CDM expectation.

.

The measurements from cluster lenses and BCG stellar velocity dispersion
measurements make it possible to constrain the density profiles of galaxy
clusters down to scales of $\mathcal{O}(10\ \rm kpc)$. The resulting
constraints on the SIDM cross section are very strong because a cross section
of $\sim$1 $\rm cm^{2}/g$ can thermalize a cluster halo out to scales larger
than 100 kpc, which would be inconsistent with existing constraints [55, 132].
The fact that the cross section is so large means that the central dispersion
$v_{\rm 1D}\simeq V_{\rm max}/\sqrt{3}$, which is much larger than the stellar
potential well thereby reducing the effect of baryons on the SIDM density
profile (see Eq. 17).

These arguments were exploited to put stringent constraints on SIDM models
[55] at the level of about $0.2\ \rm cm^{2}/g$. This method has been refined
and extended by recent work to produce more robust constraints [353, 354].
Eight clusters with generally relaxed profiles and a large set of lensed
images were analyzed in Ref. [353]. The images in these systems were spread
over a large range of radii and allowed the strong lensing data by itself to
place constraints on the density profile. Seven out of the eight clusters had
median cross section values below $0.15\ \rm cm^{2}/g$ while one of the
clusters, MACS2129, was an outlier in many respects and preferred a larger
cross section of $0.28\ \rm cm^{2}/g$. Interestingly, MACS2129 was the only
system for which one of the images couldn’t be reproduced, and it preferred an
anomalously high mass-to-light ratio compared to the other clusters, which is
relevant because smaller mass-to-light ratios correlate with smaller core
(higher DM densities). It is likely that MACS2129 requires more lensed images
or stellar velocity dispersion data to provide a better model for the
cluster’s DM density profile. Combining all eight clusters, the upper limit on
the cross section is $0.13\ \rm cm^{2}/g$ at the 95% C.L. with the mean
$\langle v_{\rm rel}\rangle$ of about $1500\ \rm km/s$ [353]. It is important
to note that systematic errors at the level of $0.1\ \rm cm^{2}/g$ are
possible [353] and more work along these lines is warranted given the
stringent constraints that are now possible.

Ref. [354] analyzed the seven clusters from Ref. [54] the span $1000-2000\ \rm
km/s$ in relative velocities. They also include constraints from the stellar
velocity dispersion data. Their upper limit from the clusters is $0.35\ \rm
cm^{2}/g$ at 95% confidence. They take adiabatic contraction of the halo into
account (which Ref. [353] do not), which weakens the upper limit by $0.07\ \rm
cm^{2}/g$. The average $\langle v_{\rm rel}\rangle$ for the sample of seven
clusters is $1900\ \rm km/s$. Group scale lenses with average $\langle v_{\rm
rel}$ of $1150\ \rm km/s$ were also analyzed by the same method and a much
less stringent bound of $1.1\rm cm^{2}/g$ was obtained with the best-fit cross
section of $0.5\pm 0.2\ \rm cm^{2}/g$. We note that the constraints in Ref.
[354] are weaker than in Ref. [353]. This could be due to the fact that Ref.
[353] created new lens models that did as well as, or outperformed, existing
models for those clusters, and the full projected density profile was used to
constrain the SIDM halo profile. More work on this topic is clearly warranted
given the stringent constraints that are possible.

####  4.1.2 Constraints from resolved subhalos and line-of-sight halos

Because strong lensing observables depend on the mass distribution of the
deflector, perturbations to the smooth mass distribution from low-mass DM
halos in and out of the lens can alter observed image positions and fluxes
[355, 356, 357, 358, 359]. Of particular interest for constraining SIDM are
the perturbations caused by low mass ($M_{\rm vir}<10^{9}M_{\odot}$) subhalos
and line-of-sight halos which are believed to be largely dominated by DM and
thus less affected by baryons than higher mass structure (see Sec. 3.8). A
handful of low-mass perturbers have been identified along strongly lensed arcs
thus far [360, 361, 362, 363, 364] through a technique called gravitational
imaging [365, 327], and one system has been re-analyzed and shown to be
consistent with a line-of-sight halo [366].

The probability of observing significant perturbations to a smooth
gravitational lens depends in part on the mass distribution of DM subhalos
(i.e. their mass function) as well as their internal mass distributions.
References [161, 367] showed that an NFW perturber for example has a larger
region of influence relative to an SIS perturber ($\rho\propto r^{-2}$) for
fixed aperture mass. This can be understood by the fact that a shallower mass
profile has to have a larger total virial mass to achieve the same aperture
mass as a steeper mass profile. This higher total mass effectively extends the
region in which the perturber can influence lensed images, thus making it
possible to jointly infer typical perturber mass distributions and the total
subhalo mass function. Reference [368] demonstrated that if a perturber falls
on a resolved arc, high enough signal to noise data can distinguish between
the effects of an SIS and a generalized NFW mass profile.

The sensitivity of strong gravitational lensing to the perturber mass profile
is promising for the use of this method as a test of SIDM, which can predict
halos that have either steeper or shallower inner density profiles than NFW
halos at these low masses, depending if they are in the core expansion or core
collapse phase of their evolution (see Sec. 3.2). Since the core collapse of
SIDM subhalos is intrinsically linked to the amount of tidal disruption they
experience [103, 129, 127], it will primarily affect subhalos of the main lens
galaxy that have small pericenters. It is thus possible that in SIDM models
with a sufficiently large cross section, the relative lensing contribution
from subhalos compared to that from the line-of-sight is modified due to the
former potentially having denser density profiles (resulting in them being
better lenses) while the latter have (mostly) low-density cored profiles,
which are more difficult to detect at high significance through strong
lensing. Moreover, SIDM subhalos that are in their core expansion phase can
experience severe tidal disruption [75] once their tidal radii becomes of the
order of their core radii (see Sec. 3.4), resulting in a change to the subhalo
mass function that might be detectable through strong lensing.

Recent work has shown that for dense enough subhalos both the mass and
effective slope of the mass profile of the subhalos [369] or its concentration
[370] could be measured. The effective slope for the perturber would be
measured in the region where we expect them to have the largest observable
effect [371]. This is analogous to the situation for the perturber mass [372,
373]. The possibility of making measurements in the slope–mass plane opens up
a new discovery mode for dark sector dark matter since CDM and SIDM models
could look very different in the slope–mass plane; more work is needed to
flesh out the predictions. Recent work in Ref. [369] indicates that the dark
substructure detected in SDSSJ0946+1006 [374] has a density profile that is
steeper than expected for dark subhalos in CDM [369]. In this context, there
is a the exciting possibility of an excess of compact perturbers seen in
galaxy-galaxy strong lensing observed in galaxy clusters [375] that don’t seem
to be present in hydrodynamical simulations based on CDM [376] but which could
be consistent with SIDM models [133]. The connections of these lines of
evidence for compact perturbers to dark sector physics is still in its
infancy. In large surveys (e.g., Euclid [377]888http://sci.esa.int/euclid/ and
LSST), the information about individually detected perturbers in the slope-
mass plane could be combined with constraints on the subhalo mass function
[378] to increase the discriminatory power. Together, the prospects for
constraining dark matter physics with strong lensing is very promising.

####  4.1.3 Constraints from the collective effect of unresolved subhalos and
line-of-sight halos

Figure 11: Fisher forecast for the substructure convergence power spectrum in
three logarithmic wavenumber bins. We consider here observations with the
wide-field camera 3 (WFC3) aboard the Hubble Space Telescope (HST) using the
F555W filter, resulting in a point-spread function FWHM of $0.07$ arcsec. The
source is placed at $z_{\rm src}=0.6$ with an unlensed magnitude $m_{\rm
AB}=24$. The error bars show the $1$-$\sigma$ regions, while the green
rectangles display the sample variance contribution within each bin. The black
points display the forecasted error bars assuming an image exposure
corresponding to 7 HST orbits, while the red error bars show the forecast for
a 70-orbit total exposure (in both cases we conservatively assume that only
half of each orbit is available for observation). For clarity, the wavenumber
bin center for each observational scenario shown have been offset by $6\%$,
with the green rectangle showing the true wavenumber bin used in the analysis.
The blue solid line shows the fiducial substructure power spectrum model used
in the forecast, which corresponds to a CDM population of subhalos modeled
with truncated NFW profiles. The dotted magenta line shows the power spectrum
for SIDM, assuming a subhalo core size equals to $70\%$ of the scale radius.
For comparison, the orange dashed line shows the substructure power spectrum
for a thermal relic warm DM with mass of 3.5 keV. Adapted from Ref. [379].

Beyond the individual detection of massive DM halos through gravitational
imaging, it is also possible to constrain SIDM by looking at the collective
impact on lensed observables of the entire subhalo and line-of-sight halo
population. Indeed, most DM halos contained in the small cylinder stretching
from the source to the observer are either too light or lie too far in
projection from the lensed images to be directly detection with high
confidence. However, their collective effect might be observable at a
statistically significant level, especially when data from multiple
gravitational lenses are combined.

Such an approach was originally taken in Ref. [380] using flux-ratio anomalies
in quasar lenses, and recently refined using an updated sample of lenses in
Refs. [381, 382, 383]. Of particular relevance to SIDM, Ref. [384] showed that
a sample of eleven quasar lenses can provide constraints on the concentration
of subhalos and line-of-sight halos. With 50 lenses, it would be possible to
put constraints on the self-interaction cross section and its velocity
dependence [378]. Beyond quasars, extended lensed images could also be used to
constrain the properties of the low-mass DM halo population (see e.g. Refs.
[385, 386]). In this case, it is useful to describe the DM halo population in
terms of its convergence power spectrum, as originally proposed in Ref. [387]
(see also [388]).

To develop intuition about whether a convergence power spectrum measurements
could be used to distinguish between CDM and SIDM, Ref. [389] developed a
general formalism to compute from first principles the power spectrum of the
convergence field for different populations of subhalos. It was shown here
that the power spectrum can be mainly described by three quantities: a low-k
amplitude, that depends on the subhalo abundance and on specific statistical
moments of the subhalo mass function; on a turnover (intermediate) scale, that
probes the truncation radius of the largest subhalos in the system; and on the
small-scale slope, that probes the subhalo inner density profile. This latter
property appears key to distinguish SIDM from CDM since a population of cored
subhalos shows a small-scale ($\gtrsim 100$ kpc-1) power spectrum slope that
is much steeper than that of a population of truncated NFW subhalos (going as
$1/k^{8}$ as opposed to the $1/k^{4}$). These findings were confirmed
numerically in Refs. [390, 391], where it was also noted that one could in
principle detect the presence of both self-interaction and of a cutoff in the
subhalo mass function by measuring the amplitude and slope of the convergence
power spectrum on scales $k\sim$ 1-10 kpc-1. These conclusions did not,
however, take into account the possible impact of subhalo core collapse.

Ref. [379] performed a detailed forecast for the sensitivity of current and
future observations to the convergence power spectrum (see Fig. 11). There, it
was found that even deep space-based images at optical wavelengths are likely
not sufficient to distinguish between CDM and SIDM, although the possible
presence of a mass function cutoff could be detectable. Higher resolution
interferometric data are likely to be required to distinguish SIDM from CDM.
So far, only an upper bound [392] on the convergence power spectrum has been
measured.

###  4.2 Stellar streams in the Milky Way

Stellar streams, particularly those produced by tidally disrupted stars from
globular clusters have very small velocity dispersions. These cold streams can
be gravitationally perturbed by substructure in the Milky-Way halos making
them excellent probes of low mass subhalos [393, 394, 395, 396, 397, 398, 399,
400]. Subhalos are expected to alter orbits of the closest stream stars
creating characteristic gaps that can be observed to infer the time of
encounter and the mass of the perturber. Observable gaps can probe subhalos in
the mass range $10^{5}-10^{6}M_{\odot}$, pushing the limit significantly below
the galaxy formation threshold, allowing us to probe dark substructure. The
long thin stellar stream GD-1 shows evidence for such a gap and a correlated
spur of nearby stars, which could be consistent with an interaction with a
perturber [401]. This perturber has to be fairly compact with approximately
$10^{6}-10^{8}$ in mass for a Hernquist profile with a scale radius of $30\
\rm pc$ or smaller [398].

SIDM predictions can change the radial distribution, the numbers and the
internal densities of the subhalos, as discussed in the previous section on
strong lensing (Sec. 4.1). This will result in quantitatively different
predictions for the kinematic and density perturbations that are imprinted on
the streams thereby providing a venue to constrain SIDM. In particular, core
collapse can lead to dense inner cores that could survive the extreme tidal
effects near the disk. Densities of order $10\rm M_{\odot}/pc^{3}$ and higher,
similar to those required in Ref. [398], are possible within the inner part of
a core-collapsed SIDM subhalo [103] but a more detailed discussion of the
implications is lacking. Like with the strong lensing case, this preliminary
analysis highlights a real possibility for discovering SIDM if more thin
streams are discovered.

The most recent data from DES analysing 5000 deg2 of sky from their year 3