{'fraction_non_alphanumeric': 0.030976176275197238, 'fraction_numerical': 0.01962923631214753, 'mean_word_length': 4.819125821834702, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 10, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Astrophysical Tests of Dark Matter Self-Interactions

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.


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 velocitydependent (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 ∼ 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 ρ ∝ 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 (ρ ∼ 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 ∝ 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 problemhave 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 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 simulationtested 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 selfinteractions. 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 O(100) 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 , 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.

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 centerof-mass frame. This simple setup can be parameterized by one parameter, the total self-scattering cross section σ.
We typically, however, discuss the quantity σ/m DM , the DM-DM cross section per unit DM particle mass m 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 where ρ is the mass density (henceforth "density") of the background particles, such that n = ρ/m DM . The total probability for the moving particle to scatter is therefore given by where Σ ≡ v 0 ρ dt = ρ dx 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 Σ ∼ 1 g cm −2 on cluster scales (halo mass ∼ 10 15 M , virial velocity ∼ 10 3 km/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 σ/m DM ∼ 1 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 θ is where σ 0 denotes the total cross section for v → 0 and w = m med /m DM denotes the ratio of mediator mass to DM mass. Scattering is velocity-independent for v 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 α DM is small. If, however, we are in the regime w α DM , the computation of the differential cross section is no longer perturbative in α 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 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 velocityaveraged cross section where f rel (v 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 σ → σ 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: with the one-dimensional velocity dispersion σ v of the order of the virial velocity. The distribution of relative velocities is then given by a Maxwell-Boltzmann distribution with velocity dispersion √ 2σ v . For an even rougher estimate, it is often sufficient to approximate σ ≈ σ(v = v rel ) where v rel = 4/ √ π for a Maxwell-Boltzmann distribution.
For a velocity weighting of n = 1 in Eq. (4), the substitution σ → σ 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 σ 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 σ 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 w, the differential cross section is strongly peaked towards θ → 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 2 DM v 2 (1 − cos θ). 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 σ = (dσ/dθ)dθ 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 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 σ , one can then define the velocity-averaged momentum transfer cross section σ T [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 θ = 0 and θ = π. To address this issue, a modified definition is [91].
Indeed, explicit numerical simulations show this latter definition is more appropriate for a differential cross section of the form of Eq. (3) [92]. 1 Alternatively, Refs. [81,86] suggest considering the viscosity cross section which weighs the cross section by the fractional transverse energy transfer. Simulations show that σ V is also more appropriate than σ 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 σ/m DM ∼ 1 cm 2 g −1 in order for there to be a significant effect for halos, which typically have DM column densities of order Σ ∼ 1 g cm −2 .

Modeling the Effective Macroscopic Phenomena Arising from Microscopic Interactions
For most systems of interest in an SIDM cosmology, cross sections O(1) 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 which, evaluated for the DM-density in the solar-neighborhood (∼ 0.4 GeV cm −3 ≈ 10 7 M kpc −3 [98]), is ∼ 0.5 Mpc for σ/m DM = 1 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 λ (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 selfinteractions, which we now describe further. The generalization of equation (1) for the interaction rate to the case of a non-stationary background is where f v is the velocity distribution function, here normalized such that f v ( r, v) d 3 v = 1. In order to determine whether a particle in an N -body simulation should scatter, the quantities f v and ρ 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 ∼ 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 (σ T and σ 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]. 2 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, ∆t, is given by where ρ 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 ρ ij . For example: Vogelsberger [74] : Rocha [106] : Robertson [107] : where The other main difference between recent SIDM implementations has been the choice of h SI . Ref. [74] used a variable h SI that adapted to the local density to keep the number of neighbours ∼ 38, while both Refs. [106] and [107] used h SI that was fixed for all particles. Ref. [106] found that h SI 0.2 (ρ/m p ) −1/3 was required for scattering to be correctly implemented (i.e. h SI must be larger than 20% of the mean inter-particle separation), and set h 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 SI , similar to the Plummer-equivalent gravitational softening length) demonstrated that scattering can be correctly implemented when h SI 0.2 (ρ/m p ) −1/3 , and that the result in Ref. [106] was likely a case of pairwise scattering probabilities saturating above 1. 3 Ref. [107] also compared their fixed h 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 SI results in many neighbors being found, but the 2 A 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. 3 For more details see Appendix A in Ref. [107]. 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 θ 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.
probability of scattering from a particular one being small (in all cases ρ ij ∝ 1/h 3 SI ), while reducing h 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).
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. 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 cm 2 /g (red). The resolution levels correspond to particle masses of 4.9 × 10 4 M (solid), 3.9 × 10 5 M (dashed), and 3.1 × 10 6 M (dotted). The vertical lines mark the gravitational softening length (2.8 times the Plummer-equivalent softening length). The figure is adapted from Ref. [74].
A test with the Hernquist profile was proposed and used by Vogelsberger et al. [74] to validate their choice of scattering kernel (Eq. 12).
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. 4 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.

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 selfinteractions 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.

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 σ/m DM 10cm 2 /g, the predictions of ΛCDM for the large scale structure correlations and the growth of halos over time carries over unchanged to ΛSIDM models. For σ/m DM 10cm 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 age , is unity. This radius is known as r 1 , which is defined through: where v rel is the average relative velocity between DM particles andρ(x) = x −1 (1 + x) −2 is the functional form of the NFW profile ρ NFW (r) = ρ sρ (r/r s ). In the above equation, V max is the maximum circular velocity within the halo (i.e. the maximum value of G N M (< r)/r) and R 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 ΛCDM dictates that R max ∝ V 1.3−1.5 max and hence (16) implies that r 1 is an approximately fixed fraction of r s [106]. For a cross section of 1cm 2 /g and t age = 10 Gyr, r 1 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 1cm 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 ρ s [116,117,118,103]. As t age × σ/m 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 where v 1D is the one dimensional velocity dispersion of DM particles inside r 1 and can be related to v rel . The gravitational potential Φ 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 σ/m DM 1cm 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 half of the stars and r 1 [119]. The central density ρ 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 lowerpanels show the velocity dispersion profiles (σ v in the figure is v 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, ρ 0 , and temperature / velocity dispersion, v 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 1 but is a mild function of halo mass or V max [50].

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 velocitydependent 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 gc /t 0 ∼ 400, where t 0 = 1/(aρ 0 v 1D σ/m DM ) is the interaction time (i.e. 1/R, where R the scattering rate defined in equation 1) in the core with central density ρ 0 = 2.4ρ s , v 1D = 0.65V max and a = 4/ √ π [100,90]. For the median densities and velocities of galaxysized halos, the collapse time exceeds a Hubble time as long as the cross section per unit mass is O(50 cm 2 /g). Note that there are no constraints that prohibit such large cross sections for velocities of order 100 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-σ) higher would collapse within a Hubble time even for a cross section of about 10 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 σ/m DM ∼ 10 cm 2 /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 ρ(r) = ρ NFW (r t ) × (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.
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] 5 , 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.

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.

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 σ/m DM > 10 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 r 1 (see Eq. 16), 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 Figure 6: Mass loss and stellar extent in cored (left) and cuspy (right) subhalos. The halos originally have masses of 3.6 × 10 10 M (solid black line), with a stellar component of 7.5 × 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].
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 smallpericenter 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].  [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 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.

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.
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 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.

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]: where R e and R 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 e,d = χ e,d R, where χ e and χ 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 expressionṡ where in the second expression we have defined F drag ≡ m DMv in order to recover the well-known form of a drag force. The evaporation and deceleration fractions χ e and χ 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σ/dθ [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Ė one typically finds R c ≈ R 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 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 kpc, cross sections much higher than 1 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 kpc for cross sections of σ T /m DM ≤ 1 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.

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 Λ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. 6 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].
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 λ γ within the plasma weakens the pressure support of the acoustic oscillations, leading to significant damping of their amplitude on scales smaller than λ γ . 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, |M| 2 ∝ (p/m DM ) n−2 , where p is momentum of the incoming relativistic species that is scattering with DM).
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-α forest Interacting DM transfer function (T (k) ≡ (P SIDM (k)/P 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, 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 WDM = 3.5 keV.
Right panel: DM drag visibility function plotted against conformal time τ 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]. [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.

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.

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 O(pc) to the baryon cycle of the circumgalactic medium at O(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 ∼ 100M 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 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 log(M * /M halo ) 10 −2 , corresponding to 10 6 M * /M 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 vir ∼ 10 8 − 10 9 M ), 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 Figure 9: Slope of the halo density profile between 1 and 2% of R vir vs M * /M halo for 31 galaxies from the MAGICC simulation. Reproduced from Ref. [264].
next-generation surveys such as the Legacy Survey of Space and Time (LSST; [268,269]) 7 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.
Hydrodynamic simulations in Λ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 ΛCDM one in galaxies where the gravitational potential is dominated by stars and similar to V 2 max [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 ultrafaint 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 [279, 280] and 10 7 < M * /M < 10 8 [281] galaxies, inner DM density slopes are little affected by baryons.

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-α 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 1.2 × 10 −10 m/s 2 : above this the dynamics is Newtonian, while much below it F ∝ a 2 , leading to asymptotically flat rotation curves. a 0 corresponds to a critical baryon surface density Σ 0 ≡ a 0 /G 1 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].

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.

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.

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 ρ tot ∝ r γ with γ ∼ −2, and small intrinsic scatter of 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 β DM ≈ 0.5 − 2, where ρ DM ∝ r −β 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 ∼ 10 13 M 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 ∼ 10 15 M halos often present several multiple image systems with Einstein radii ∼ 50 kpc. Spatially resolved stellar velocity dispersions within the BCG probe the mass distribution within ∼ 20 kpc of the cluster center, while weak lensing and X-ray emission probe the mass distribution on ∼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 ∼ 30 kpc and has a mean inner slope of β DM = 0.50 ± 17 (c.f. β 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 ∼ 10 14 M . 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.
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 O(10 kpc). The resulting constraints on the SIDM cross section are very strong because a cross section of ∼1 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 1D V max / √ 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 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 cm 2 /g while one of the clusters, MACS2129, was an outlier in many respects and preferred a larger cross section of 0.28 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 cm 2 /g at the 95% C.L. with the mean v rel of about 1500 km/s [353]. It is important to note that systematic errors at the level of 0.1 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 km/s in relative velocities. They also include constraints from the stellar velocity dispersion data. Their upper limit from the clusters is 0.35 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 cm 2 /g. The average v rel for the sample of seven clusters is 1900 km/s. Group scale lenses with average v rel of 1150 km/s were also analyzed by the same method and a much less stringent bound of 1.1cm 2 /g was obtained with the best-fit cross section of 0.5 ± 0.2 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.

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 vir < 10 9 M ) 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 (ρ ∝ 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] 8 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.

Constraints from the collective effect of unresolved subhalos and line-of-sight halos
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 lowmass 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 ( 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 ∼ 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.

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 , 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 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 10M /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 deg 2 of sky from their year 3 data have found 11 new streams within 50 kpc of our location in the MW, these streams are observed down to surface brightness of 23 mag/arcsec [402]. While currently DES is the state-of-the-art, with the most sensitive wide-area view, with the advent of the Vera C. Rubin Observatory we are expected to see of the order of a hundred streams out to the virial radius of the MW halo [402]. Moreover, spectroscopic programs like the S 5 survey and instruments like the Maunakea Spectroscopic Explorer telescope will allow us to follow up on some of the streams with spectroscopic data to observe the disruptions of the streams in velocity space [403,404,405].

X-ray and Weak Lensing Observations of Clusters, Groups and Large Ellipticals
The overall distribution of DM in halos can be probed in a variety of ways. While rotation curves of galaxies can be measured where spectroscopic data is available, the gas distribution in X-ray and gravitational lensing have been used to probe the overall distribution and density profile of the DM out to the viral radius for larger systems. Current and future generation surveys like the Dark Energy Survey (DES; [406]) 9 , Hyper Suprime Camera (HSC; [407]) 10 survey and LSST will provide precision measurements of the DM halo density profiles through weak lensing of galaxies in the back ground of large, statistical samples of cluster and group mass halos. While the formation of a core is the most dramatic effect on the density profile of DM halos in the presence of self-interactions, mass conservation implies that the density at distances larger than the core becomes steeper than the density profile in the CDM scenario. While the latter is a relatively smaller effect, the scale at which this feature is expected, especially for clusters, is less affected by baryonic physics than the core. Further, the projected density profiles of clusters on these scales can be measured accurately through weak lensing measurements. Fig. 12 from [95] shows the stacked projected density profile, ∆Σ of cluster mass objects from N -body simulations of SIDM with different cross sections. The red curve in the figure corresponds to the current best constraints on the cross section from the Bullet cluster [107] of 2 cm 2 /g. The bottom panel shows the ratio of the expected weak lensing signal for various SIDM cross sections and CDM. The lightest grey bands show the measurement error bars for a similar sample of objects from the DES Year 1 data. The progressively darker bands show the projected error bars from DES Year 3 and LSST, taking into account their deeper coverage. This observable, therefore, can potentially provide much stronger constraints on the SIDM cross section at cluster scales, especially once data from LSST becomes available.
Apart from the overall density profile, as mentioned previously DM self-interactions can alter the shape of the DM halo, making the regions of higher density more spherically symmetric. Constraints on scattering cross sections from shapes of halos have been obtained using the X-Ray measurements of the shapes of the hot halo gas of NGC720. The shape and the twisting of X-Ray isophotes carry information about the three dimensional shape of the potential. Measurements from NGC720 constrain the halo ellipticity to be ∼ 0.35 ± 0.03, consistent with elliptical halos from CDM, this values is consistent with σ/m DM < 1 cm 2 /g [165].
While the shape of DM halos was previously thought to be a promising probe of the nature of DM, subsequent studies have found that introducing baryonic physics can significantly complicate the inferences of cross sections from halo shapes. For example, on galaxy scales [121] show that the asymmetric potentials from central disks can make the shape of the DM halos aspherical near the center on short timescales of a few Gyr depending on the size of the disk. Note also that on cluster mass scales, while the DM can remain spherical near the center for cross sections of σ/m DM < 1 cm 2 /g, as shown in Fig. 7 from [164], the stars and the gas within the halo remain aspherical for a longer time and their shapes are much slower to respond to the change in the DM potential.
Weak lensing measurements of shapes of DM halos, stacked on the galaxy light profiles or using three point functions can measure the shape of dark halo directly as a function of radius. [408,409,410] measure shapes of Luminous Red Galaxies and cluster mass halos using weak lensing and distribution of satellite galaxies and find that the ellipticities are consistent with CDM, [410] for example place strong constraints on the axis ratios of RedMaPPer clusters in SDSS at 0.573 ± 0.039, implying that the halo is overall quite aspherical.

Major and Minor Mergers in Groups and Clusters (and associated observables)
Just like terrestrial particle accelerators allow us to study the interactions of fundamental particles, merging galaxy clusters can be thought of as "cosmic accelerators" probing the properties of DM. Indeed, the Bullet Cluster is arguably among the most sensitive and certainly the most famous probe of DM self-interactions [169]. What makes this system so unique is that the gas in the smaller of the two merging clusters exhibits a shock wave, which enables us to conclude with certainty that the system is being observed shortly after core passage. Combining this information with a lensing reconstruction of the system's DM profile [411] enables us to infer the velocity and impact parameter of the collision [412,413,414]. Typical estimates yield a velocity of around v c ∼ 4000 km s −1 at core passsage and a projected (i.e. integrated) DM density of Σ ∼ 0.3 g cm −2 along the path of the subcluster. While the frequent claim that the Bullet Cluster demonstrates DM to be collisionless is clearly exaggerated, it is true that it places strong bounds on the DM self-interaction cross section. Given the projected DM density of the main cluster halo along the path of the subcluster, the main halo is expected to become optically thick for DM self-interactions with σ/m DM 3 cm 2 g −1 . This means that cross sections larger than this will lead to almost complete destruction of the subcluster's DM halo, as can be seen in Fig. 13 for the case of a simulated Bullet Cluster with σ/m DM = 10 cm 2 g −1 .
The most stringent SIDM constraint from the Bullet Cluster comes from the measurement of the mass-to-light ratio (M/L) within 150 kpc of the subcluster centre, which is found to be consistent both with that of the main cluster and with the typical range of M/L ratios for clusters [169]. Ref. [167] argues that this agreement places an upper bound on the integrated mass loss by the subcluster and that it cannot have lost more than 20-30 percent of its mass within the 150 kpc aperture. The resulting bound on isotropic, hard-sphere scattering is approximately σ/m DM 0.7 cm 2 g −1 . It should be noted that this result is subject to some caveats, most notably that there is significant intrinsic scatter in M/L of isolated clusters. Ref. [415] find a 55% scatter in M/L within the virial radius of clusters, which would presumably increase if only looking within a small aperture. The consistency in M/L between the main cluster and subcluster therefore does not exclude the subcluster having lost a significant fraction of its pre-collision DM mass, if it started with high M/L.
An independent constraint arises from the fact that the peak of the subcluster's projected density inferred from weak lensing measurements coincides with the galaxy centroid from visible light to within 1σ. This observation rules out the presence of a substantial drag force, which would lead to the deceleration of the DM halo and hence a separation between DM and galaxies [167]. The observational results is that the DM is lagging behind the galaxies by 25 ± 29 kpc [416,154]. The first numerical simulations of the system with SIDM translated this observation into the constraint σ/m DM < 1.25 cm 2 g −1 [154]. However, more recent work has found this constraint to be overstated, with a cross section of this magnitude producing offsets of only 10 − 15 kpc [91,107]. This difference is driven by the methods used to find the "positions" of the subcluster's DM and galaxy components, with the latter papers using methods more closely related to what was done observationally.
In spite of the discovery of numerous merging galaxy clusters in recent years, the Bullet Cluster remains the single most sensitive system, owing mostly to the superior projected density of the main cluster. Progress can potentially be made by the statistical combination of a large number of major and minor mergers, with the first attempt at this producing an upper limit on the SIDM cross section of σ/m DM < 0.47 cm 2 g −1 [112]. However, both the analytical model used to relate a DM-galaxy offset (strictly, the ratio of the DM-galaxy offset to the gas-galaxy offset) to a cross section, as well as the quality of the lensing maps used to infer DM positions, have been called into question [173]. It remains to be seen whether further discoveries of merging clusters can improve on the bounds obtained from the Bullet Cluster.
One system that has received particular interest in recent years is the Abell 3827 cluster. The advantage of this system is that strong gravitational lensing allows for a very precise reconstruction of the DM distribution. It has thus been possible to identify four galaxies within the core of the main cluster that all retain their own DM halos [172,418]. Intriguingly, early studies found a substantial separation between at least one of these DM halos and the associated stars [419], which has been interpreted as positive evidence for SIDM [420]. Unfortunately, more detailed observations have not confirmed this separation and seem to imply that Abell 3827 is fully consistent with the predictions for collisionless DM [421], in which the peaks in stellar and DM density should be coincident [422].
Another controversial system is the 'Train Wreck Cluster' Abell 520, which has been claimed to exhibit a dark core, i.e. a region of enhanced M/L ratio close to the collision point [423,424] (see however Ref. [425] for a contradicting claim). While the most plausible explanation for such a dark core may be a projection effect due to the presence of filaments [426], an interpretation in terms of SIDM is also conceivable. It is, however, unclear how elastic scattering of DM particles could lead to the formation of such a dark core, since the scattered DM particles would typically have too much kinetic energy to become gravitationally bound to each other. While it may be possible in principle to accommodate such an observation in more complex SIDM models (as discussed in Section 5), it remains challenging to bring the required self-interaction cross section into agreement with the bounds from the Bullet Cluster.
Finally, we mention the Musket Ball Cluster, which owes its name to the fact that it is both older and slower than the Bullet Cluster [427]. This system was initially considered to be particularly sensitive to DM self-interactions, due to the fact that it is much further progressed than other merging clusters. More detailed studies have however revealed, that it is typically not advantageous to observe a system long after core passage, because the DM halos quickly relax to their equilibrium position, so that the separation between DM and stars decreases rather than increases with time [91,149,92].
To summarize, systems of merging clusters with the following properties would make them ideal candidates for constraining DM self-interactions: 1. The integrated density of DM along the axis of merger for the system should be as large as possible, so that self-interactions have a larger effect.
2. The system should allow for a precise reconstruction of the DM profile via (strong) gravitational lensing so that any offset between the light profile and the center-of-mass of the DM can be accurately measured.
3. The system should be observed as shortly after core passage as possible, to prevent the DM Figure 14: Figure from [428] showing the evolution of the expected log-normal median DM-BCG offset in observations. The magenta symbol is the median BCG offset from the Local supercluster substructure survey measured in [429].
particles from relaxing back to their equilibirum positions.

The collision should be observed face on to minimize projection effects.
Even if no single system is discovered that surpasses the Bullet Cluster in terms of sensitivity to SIDM, progress can still be made through the statistical combination of many colliding galaxy clusters. Expanding the catalogue of such systems is therefore a highly promising avenue in the search for SIDM.

Incidence of Multiple BCGs, their Stellar Kinematics and their Spatial Separations
In Sec. 3.3, we discussed the counterintuitive result that galaxies do not virialize if bound to the core of a cored DM halo. In simulations of equal-mass mergers, Ref. [149] found that such galaxies oscillate about the center of the halo with an amplitude equivalent to the core radius. The easiest galaxies to definitively associate with cores are the biggest, the BCGs. If the models of Ref. [149] are correct and hold for more typical mergers where the masses of the two merging systems can be significantly different, then BCGs should be offset in both position and velocity from the center of a relaxed cluster. These studies from dark-matter-only simulations were followed up in [428] to show that the effect of the BCG "wobble" in a cored profile is present even in full hydrodynamic simulations. It is in this context that some key observational results about galaxy cluster centers become even more interesting. Because measurement of the density profile of clusters from weak lensing depends critically in accurate determination of cluster centers, the weak lensing community is highly invested in identifying the best observational tracer of halo centers. From numerous studies, the main finding is that the position of the central galaxy (usually, but not always, the BCG [430,431,432]) is by far the least noisy and most accurate tracer of the cluster center [433,434,431]. Ref. [431] finds that the median offset between the BCG and the X-ray center of massive clusters is about 10 kpc. If true, and if the equal-mass merger prediction of Ref. [149] holds more generally, the implied constraint on the SIDM hard-sphere cross section is of order 0.1 cm 2 /g. Ref. [429] claim a detection of BCG offsets, which can be verified with the large and well-studied cluster samples from current wide-field optical surveys. The measurement from [429] is presented in Fig. 14, along with the simulation results for the offset from [428]. The current constraints from the measurement of BCG wobbles in the cores of galaxy clusters is σ/m DM < 0.39 cm 2 /g at the 95% C.L. [428] and data from future surveys like Euclid are expected to further improve these tight bounds.
There are several important lines of investigation that need to be pursued to turn this promising probe of SIDM into a robust observable. First, more simulations are required to determine how generic the galaxy sloshing in cores is as a function of cluster assembly history (e.g., Ref. [49]). Second, it is important that the position and velocity offsets be evaluated in the context of both CDM and SIDM theory. This requires realistic mock observations to be created from simulations. Third, from the observational side, the position and velocity offsets for BCGs should be measured more precisely, and as a function of halo mass. In the case of position offsets, enormous datasets are already in hand from various imaging surveys, on account of the great interest in clusters in the dark energy/wide-field survey community. Fourth, these offset measurements should be performed not just for a single BCG, but for the handful of brightest cluster galaxies, to trace the merger history of the cluster. In the extreme case of an equal-mass merger with each component initially seeded with its own BCG, we expect two sloshing BCGs in the merger remnant sitting on opposite sides of the cluster center.

Dwarf Galaxies in the Local Group and beyond
Dwarf satellite galaxies are observed around the MW and Andromeda [435,436], and around galaxies in the local Universe [437]. In hierarchical structure formation scenarios, these dwarfs are hosted by subhalos bound to the parent galaxy [32,31]. There are about 60 known satellite galaxies of the MW [438], many of which have measurements of stellar kinematics from large ground-based telescopes and orbital properties (such as pericenter distances) constrained with Gaia data. Several surveys are now able to detect satellites in MW-like galaxies in the Local Volume and this number will grow in the coming years [439,437,440,441,56,442,443,444,445,446,447,448,449,450]. The implications of these measurements for general SIDM models remains to be fully explored.
The dispersion-dominated dwarf galaxies-known as dwarf Spheroidals (dSphs)-in the MW span a large range in stellar sizes and luminosities, from ultra-faint dwarfs such as Segue 1 with stellar mass of O(10 2 )M and half-light radius of 30 pc to Fornax with about 10 7 M of stellar mass within about a kpc. This is shown in the left panel of Fig. 15. The dynamical masses of these dSphs are significantly in excess of their stellar masses, making them some of the most dark-matter-dominated systems observed. A simple way to estimate the dynamical masses of the dSphs is to use the results of Refs. [451,452], who argue that the mass within the half-light radius can be related to the observed stellar dispersion for a wide range of halo profiles and anisotropies (see also [453]). In the right panel of Fig. 15, densities within the half-light radius estimated in this way are shown revealing a striking anti-correlation of density with half-light radius [131,454].
Comparisons between simulated subhalos and observed dSph populations reveal apparent discrepancies that have been the subject of intense scrutiny for the past two decades. The most obvious mismatch is between the number of simulated subhalos with virial mass larger than M vir 10 8 M (the mass corresponding to the smallest virial temperature, T vir 10 4 K, that enables collisional cooling of neutral hydrogen gas and hence star formation), and the much smaller number of observed Galactic satellites. This "missing satellites" problem [20,19] provided motivation to models like 'warm' dark matter (WDM), that invoke primordial velocity dispersions to suppress gravitational collapse on scales below the thermal free streaming length [282,455]. Other proposed solutions invoke the astrophysics of galaxy formation; in these scenarios, the vast majority of DM subhalos fail to host galaxy formation, due to a combination of factors including mass thresholds, orbital history and cosmic reionization [35,456]. Incompleteness and biases in x 2 x 2 Classicals, rP ≤ 40 kpc Classicals, rP > 40 kpc Ultra-faints, rP ≤ 40 kpc Ultra-faints, rP > 40 kpc Figure 15: Plot on the left from Ref. [438] shows the MW dSphs in the V-band magnitude vs. half-light radius plane showing how deeper surveys have led to more dwarfs being uncovered and how the dSphs properties compare against MW globular clusters. Plot on the right from Ref. [131] shows the inferred total mass density (almost all DM) within the half-light radius and the strong anti-correlation between the mass density and half-light radius and the large diversity in the properties of these dSphs. In the background are subhalo mass profiles from a MW halo (Kauket from ELVIS) [465].
observations have also been a significant factor in exacerbating the problem [457,65].
A satisfactory resolution of the luminosity function of the satellites in the context of the mass function of subhalos is still missing. Recent work has pointed out that there may be, in fact, no missing satellites problem [458,68] and that the radial distribution of the satellites can be explained in the context of CDM models. This requires, however, that satellites be hosted in subhalos with masses below the HI cooling temperature [69,459]-a puzzling aspect that deserves further study. This conclusion seems to critically depend on the radial distribution of satellites and the impact of the Large Magellanic Cloud, with Ref. [460] finding no significant evidence for populating satellites below the HI cooling limit. Another critical factor is the limited resolution of N-body simulations, which can lead to artificial destruction of what could be ultra-faint dSphs at the present time at radii below 50 kpc [461]. These are clearly important issues to sort out and the development of detailed models like those in Refs. [461,460] will facilitate this exciting science as more satellites are discovered by the Rubin Observatory.
The key take-away message is that there are recipes for matching the luminosity function of satellites (taking into account the observational incompleteness) to the subhalo mass function. However, a more detailed understanding of the physics responsible for this empirical matching is not yet available and without that, we cannot claim to have fully understood the "missing satellites" problem. Regardless of the eventual understanding, it is clear now that there is limited room for models that do not form subhalos below a certain mass scale (e.g., WDM, SIDM with additional hidden sector interactions or fuzzy DM [151,462,150,151]).
In the typical scenario where the luminosity function is abundance-matched to the subhalo mass functions in CDM simulations to assign dwarf galaxies to subhalos, the brighest dwarfs would reside in the most massive subhalos. However, the measurements of the central density of dwarf spheroidals (dSphs) are less than the central density of subhalos in CDM simulations (darkmatter-only). Thus, the stellar kinematics of dSphs imply that the most massive subhalos in the simulations are dark, which is a puzzle as they are "too big to fail" to form stars [463,464] (i.e., they would have retained their gas and cooled to form stars).
There are multiple resolutions to the too-big-to-fail problem for satellite galaxies. The most promising among them is that the discrepancy is just an artifact produced by comparing obser-vations to the dark-matter-only simulations. Simulations with baryons indicate that the subhalo density profiles are much flatter than in dark-matter-only CDM simulations. A flatter density profile goes a long way toward matching the central density, as measured observationally, to the total halo mass [466,467,468]. Moreover, including the disk of the MW is known to have an impact on the satellites that have pericenters close to the disk [469,467]. However, not all dSphs have come close to the disk of the MW, so the validity of this solution remains to be fully tested. All the key pieces of the too-big-to-fail problem in CDM and SIDM models, with the exception of core collapse is visible in Fig. 16 where the impact of the disk and pericenter distances in models with with σ/m DM = 0 and 1cm 2 /g are shown.
The presence of elastic self-interactions changes the evolution of a satellite within the host halo in many different ways as listed below and discussed in more detail in Refs. [152,134].
• Dark ram-pressure stripping: the high velocity interactions of the subhalo particles with the host halo particles can lead to mass being unbound from the subhalo.
• Deceleration of the subhalo due to self-interactions and enhanced dynamical friction because of orbital decay.
• Thermal evolution, including core collapse and expansion, of the core impacted by the mass removal by tides.
• Response of the DM and stellar density profiles to the removal of mass by tidal stripping and ram-pressure stripping.
The different aspects of the orbital and thermal evolution of a subhalo are sensitive to the velocity dependence of the self-interaction cross section [153,470,127] and the presence of the disk [471,472]. The thermalization of the subhalo is controlled by the cross section at relative velocities of about 4V max, subhalo / √ π, while the ram-pressure stripping is controlled by the cross section at relative velocities of order V max, main halo [75,155,153]. The gravothermal evolution of the SIDM subhalos is complicated by the fact that the subhalos' interactions drive the thermal evolution in the presence of mass loss and structural changes due to the tidal forces of the host halo [103,127]. It has been proposed that the too-big-to-fail problem can be resolved by cored halos in self interacting DM produced by cross sections of around 1 cm 2 /g [74,473,117], i.e, for cross sections lower than the threshold needed to trigger the gravothermal collapse phase. Refs. [474,475,51,131] have shown that this is unlikely given the large differences in the central densities of the dwarf galaxies. A major issue is that the least luminous satellites lie in dense halos, even as the more massive, classical satellite galaxies are low in central density [453,476,151]. Cross sections preferred by the more massive satellites dramatically underpredict the densities of the ultrafaint dwarfs [151]. Perhaps, the best match in terms of scatter to the MW satellite densities is achieved with a cross section of 10 cm 2 /g [74] or even larger once the gravothermal collapse mechanism is considered in velocity-dependent SIDM models [51,471,470]. In this context, it is important to keep in mind that the core collapse process only needs to increase the central densities by a factor of few typically to achieve consistency with measured densities. For a large enough cross sections, much of the evolution of the subhalo could be spent in the core collapse phase with central densities of between few to ten times times ρ s [103].
The discussion of the proposed solutions based on SIDM is incomplete without including the impact of the MW disk, which could have more dramatic effects for SIDM subhalos compared to CDM subhalos for two reasons. One, recent work shows that the core collapse of SIDM subhalos driven by the gravothermal catastrophe is accelerated by mass loss induced by tidal forces [103,276,127]. Two, the SIDM subhalos that do not undergo core collapse will have large cores that makes them more prone to tidal effects leading to further core expansion [478,479,153]. The combination of these two effects leads to an enhanced Figure 16: Circular velocity profiles for the ten subhalos with the largest peak velocity across time in the simulations with the disk and bulge potentials (right) and the DM-only simulations (left). Note that these are not necessarily the subhalos with the largest peak velocity today. The SIDM cross section over mass is set to 1cm 2 /g, which implies that there is no possibility of the subhalos entering the core collapse phase. The data for the classical MW satellites brighter than 2 × 10 5 L [452] are denoted with filled squares with point sizes proportional to the logarithm of their stellar masses. The points and curves are color coded by the inferred pericenter distances; for the MW satellites the pericenter distances calculated in [477] using Gaia DR2 have been used. Gray symbols are satellites whose inferred pericenter distances have errors of the same order as their mean value, and are hence uninformative. The thin dashed (faded) lines indicate unresolved regions in the simulations. The CDM halo has three subhalos as dense as Draco, while the σ/m DM = 1 cm 2 /g halo does not have any. Note that for all models, with and without the disk, there is a trend that subhalos with smaller pericenters (cyan lines) have lower densities. The MW satellites do not seem to adhere to this trend. Figure adapted from Ref. [277]. diversity of the inner DM distribution of SIDM subhalos for constant cross section SIDM models [276,127]. Notice that velocity-dependent SIDM models with large cross sections at the scale of dwarfs without the presence of the tidal effects of the MW disk already produce a diverse subhalo population [51]. It would then be expected that such environmental effects will enhance even further the diversity observed in such models. Recently, Ref. [471] modeled the combined effects of the gravothermal collapse, with the tidal interactions due to the MW disk using orbital information from the MW satellites to find that a velocity-dependent model similar to the one in Ref. [51] is able to fit the diversity of the MW satellites. The presence of evaporation can delay or halt core collapse [127], which provides an independent handle on the velocity dependence of cross section since evaporation is sensitive to the cross section. For sharp velocity dependence that imply that the cross section is 1 cm 2 /g at velocities of 200 km/s, evaporation is a small effect [472].
The too-big-to-fail problem relates to the mean DM density within galaxies' half-light radii. Complementary information can be obtained by measuring the density within other radii to get a better sense of the full DM profile. This is more challenging although a few different methods have been proposed: • The method of Ref. [480] to use multiple stellar populations with different half-light radii to infer the dynamical mass within different radii and hence infer slope of the DM density profile on those scales. Cored DM profiles are better fits than NFW to both Sculptor and Fornax according to some studies [481,482,480,475,476], but other studies argue the situation is not clear [483,484]. Careful mass modeling of classical dwarf galaxies indicate that the inner slopes are heterogeneous, with some dwarfs preferentially in cusped halos [475,476,485] • Using the projected separation of globular clusters (GCs) in Fornax and arguing that dynamical friction would have sunk these GCs to the center of the dSph if the DM subhalo density profile were cuspy [486,487]. This leads to the inference of a large (∼ kpc) core in Fornax [488]; although this issue remains controversial and recent studies have found consistency of a cuspy Fornax halo with its GC population, e.g. see Ref. [489,490].
• The existence of cold substructure in the dwarf spheroidals is another cited evidence for a cored halo due to the reduced tidal effects as compared to a cuspy central density [491].
Given the current situation, a key avenue of progress is to identify observables in dwarf galaxies that are truly distinct between CDM and SIDM (see also Sec. 3.8). A few promising examples can be mentioned: • A relevant observational trend reported recently is that of the correlation between the orbital pericenter distances and the central density measured for the MW dSphs [131,454]. The central densities are larger for dSphs that have closer pericenters, which makes it difficult to see how the MW disk is solely responsible for the diversity in central densities observed across the dSphs. On the other hand, this could fit well within the predictions of SIDM models if the high density dSphs (such as Draco and Ursa Minor) were in the core collapse phase accelerated by tidal effects, while low density dSphs (such as Fornax and Sextans) were in the usual core expansion phase [495,276].
• The efficiency of supernova feedback as a cusp-core transformation mechanism is strongly dependent of the energy released relative to the change of potential energy in the halo from being cuspy to cored [496,497,498,499]. The efficiency peaks at around M h ∼ 10 11 M , just around the upper mass end of dwarf galaxies, and drops sharply towards higher and lower masses [500,501]. For smaller dwarf galaxies, the stellar-to-halo mass ratio is expected to be too low to be efficient at forming large cores [492] (although see Ref. [502]). On the other hand, in the SIDM scenario with σ/m DM ∼ 1 cm 2 /g, cores of a size comparable to the scale radius of the halo are expected across all masses [106] in the core expansion Figure 17: Two potential ways to disentangle CDM + supernova feedback from SIDM as mechanisms of cusp-core transformations in dwarf galaxies. The upper panels from Ref. [492] show the slope α (within 0.5−1% of the virial radius) of the halo density profile as a function of stellar mass for a full baryonic physics implementation (solid symbols; FIRE implementation, see Ref. [493]) and the corresponding DM-only simulations (open symbols). CDM in the upper left and SIDM (σ/m DM = 1 cm 2 /g) in the upper right. Note that this assumes that all the subhalos are in the core expansion phase, which seems to be disfavored currently (see text). For SIDM models where some of the subhalos are in the core collapse phase, we would expect greater diversity in the slopes, possibly spanning the whole range from slopes from 0 to -2. For small stellar-to-halo mass ratios, the supernova feedback mechanism is (energetically) less efficient in producing sizeable cores, while the SIDM mechanism (constant cross section σ/m DM 1 cm 2 /g) remains efficient at all masses. The bottom panels are from the isolated controlled DM-only simulations of a dwarf-size halo from Ref. [494]. The halos in both panels result in a ∼ 0.5 kpc core after t ∼ 4 Gyr, with the simulation on the right having a toy supernova feedback model, while the one on the right has a standard SIDM implementation (σ/m DM = 2cm 2 /g). The plots show the radial 2D phase space density of 2000 tracers originally set as an orbital family (narrow energy and angular momentum distribution) at the end of the simulation. Supernova feedback is an impulsive mechanism that irreversible perturbs the orbital family creating distinct shells in phase space, while the SIDM case is an adiabatic mechanism, which leaves the orbital family of traces largely unperturbed.

SIDM CDM
phase. In the core collapse phase, we would expect much greater diversity in the values of the slopes as discussed in this section. For the lowest concentration halos or halos with larger V max where the cross section is small enough, we could have cored halos. For smaller mass halos or for halos with high concentrations or more tidally-stripped halos, one would expect slopes between 0 and -2 depending on the size of the core and the radius at which we are able to measure the slopes. Unambiguously identifying and characterizing sizeable cores for stellar masses 10 6 M is a promising avenue to constrain physical mechanisms of core formation (see top panel of Fig. 17).
• Cores formed due to DM self-interactions and due to supernova feedback have been shown to be widely degenerate (when it comes to the shape of the density profile, see Ref. [503]). However, the supernova-driven gas blowouts that are responsible for expanding the orbits of DM particles to create cores need to be impulsive (relative to the orbital timescales), while self-interactions are effectively an adiabatic process of energy redistribution leading to the formation of cores [256,494,499]. This intrinsically different nature leads to distinct signatures in the orbits of stars: shells in phase space are formed in the impulsive SNfeedback scenario, while they are not present in the SIDM case [494]; see bottom panel of Fig. 17. The presence (or absence) of such shells in the stellar distribution of dwarf galaxies with confirmed cores would strengthen (weaken) the supernova feedback scenario as the main mechanism of core formation.

Rotation curves of spiral galaxies
Problems associated with explaining the DM halo density profiles persist in higher mass galaxies, where they manifest themselves in several different but related forms. Unlike the situation with dSphs where it has proven difficult to measure unambiguously the presence of constant density cores, the situation in rotation-supported galaxies is clearer. The evidence for cores in many spiral galaxies is strong, particularly those of low surface brightness (e.g. [41,504,505,506,468]). This evidence is sufficiently ubiquitous to lead some to suggest a universal cored density profile with a well-defined relation between core density and size [507,508,509]. NFW fits to the rotation curves of individual galaxies of a wide range of type often perform poorly or require concentrations significantly outside the ΛCDM mass-concentration relation scatter (e.g. [510,511]). Fits to individual galaxies can suffer from systematic issues related to inclination angle estimation, non-circular motion and presence of bars [512,513,514,515]. The results discussed in the previous paragraph can be extended in a statistical sense across the galaxy population by employing dynamical scaling relations between the masses, sizes and characteristic velocities of galaxies. Empirical methods for assigning galaxies to the ΛCDM halo population from N -body simulations (e.g. subhalo abundance matching [516,517]) generate a larger ratio of dark to luminous matter than observed across the rotation curves of galaxies. This is manifest in the fact that the ΛCDM halo mass function with NFW (or similar) profile normalises both the Tully-Fisher and Faber-Jackson relations too high [518,519,520], and overestimates velocities in the lowacceleration outer regions of galaxies [521]. This may be seen as extending the Too Big To Fail problem [463] to higher masses. There is evidence that the overprediction of dynamical mass in ΛCDM is in fact more general than implied solely by these relations [522]. On the other hand, relaxing assumptions about the density profile of the dark-matter halos, allowing for cored halos, leads to better consistency between the halo masses preferred by rotation curves and those implied by abundance matching [523,524,468].
Several other aspects of galaxy phenomenology have proven challenging to ΛCDM and may hold clues to the nature of DM. Scaling relations like the baryonic Tully-Fisher [525] and mass discrepancy-acceleration [526,527] (or radial acceleration; [528]) relations and Fundamental Plane appear to have smaller scatter than the ΛCDM halo mass-concentration relation and galaxyhalo connection (e.g. M * − M vir ), from which they should receive contributions [529,530,521].  Figure 18: The radial acceleration relation from SIDM fits to low-redshift galaxies with σ/m DM = 3 cm 2 /g. The black solid and red dashed lines show an empirical fitting function (and its 1σ scatter) that describes the observational data well. The inset shows residuals with the fit function removed. From [50].

Number of Points
The diversity of rotation curve shapes has also been seen as problematic for CDM and may point towards DM interactions [45,48,46] or strong feedback [271]. However, recent work has argued that models with cuspy dark matter halos (i.e., ΛCDM halos without strong feedback) are likely to be more consistent with the observed rotation curves if the gas disks in late-type spirals are like those simulated in EAGLE simulations [274,275]. The origin of this conclusion lies in the discrepancy between recovered and actual rotation curves caused by systematic effects related to non-circular motion, thickness of gas disk and departures from dynamical equilibrium in a systematic way that makes cuspy halos more likely to look like cored halos. If this conclusion holds true, this may be an indication that core collapse due to self interactions is more common in dwarf halos than discussed in the SIDM investigations thus far.
Other observations that are potentially surprising in the context of ΛCDM include spatially and kinematically correlated planes of satellites around the MW and nearby galaxies [531], a dearth of DM in the central regions of high-mass ellipticals [532], and the strong point-by-point correlation of kinematics and baryonic surface density in spirals [300,528].
DM self-interactions may offer a solution to some of these problems. To explain the smallest core densities in the most extreme LSBs, the cross section needs to be 3cm 2 /g [50]. Such large cross sections are disfavored from cluster measurements as discussed in Sec. 4.1.1. These arguments therefore directly motivate the velocity-dependent cross section model, which we have argued is generic in the context of dark sector models. The velocity dependence allows for large effects on galaxy scales while maintaining consistency with constraints from cluster lensing studies. Small cross sections (well below 1cm 2 /g) may also work but will likely not be distinguishable from CDM.
The mass discrepancy-acceleration (or radial acceleration) relation has also been explicitly studied in the context of SIDM [50] (see Fig. 18); the good agreement with observations provides support for the redistribution of DM mass effected by self-interactions.
Solutions based on baryonic feedback, modified gravity and other DM physics are also possible however (see Sec. 3.8). High-redshift studies of galaxy kinematics (especially at low mass) may provide a way to disentangle SIDM from the effects of baryons by measuring the time scale of core formation. This is likely to differ between the different models, although generic predictions for SIDM are lacking.
Currently, the majority of detailed studies of kinematics in low mass galaxies are limited to the Local Volume, however there are a few cases where spectroscopic measurements of strongly lensed galaxies up to z 2 have been made (e.g. [533,534,535,536]). It is also possible to measure the Tully-Fisher relation as a function of redshift (z 2) to study the relative evolution of galaxies' stellar and DM components. Evidence for variation in the Tully-Fisher normalisation is however controversial, with some authors finding no significant trend (e.g. [537,538,539]) and other a trend towards relatively more DM at high z (e.g. [540,541,542,543]) The uncertainties in each of these methods are large due to the difficulty of reliably measuring the kinematic tracers and small sample sizes but they hold promise for testing the physics being implemented in hydrodynamic simulations.
Advancing this field will require a synergy of observational and modelling techniques. On the one hand, empirical knowledge will benefit both from the compilation of homogeneous datasets with high quality photometry and kinematics (e.g. [544]), and from the increase in statistics and precision afforded by current and next-generation IFU studies such as MaNGA [545]. Problems primarily based on photometry will benefit from DES, LSST, Euclid and Nancy Grace Roman Space Telescope [546] 11 data. Large future datasets will enable stringent cuts on quantities like inclination and the degree of pressure support, reducing their confounding effects. On the other hand, advances in modelling will sharpen the likelihood function for dynamical observables for given theoretical assumptions. This requires not only further modelling of a range of systems in SIDM, but also increased resolution in ΛCDM hydrodynamic simulations-to push more of the relevant physics above the grid scale and hence reduce uncertainty in the effect of feedback-and an increase in the sophistication with which models and data are compared.

Ultra-diffuse galaxies
Long suspected of dominating the low-luminosity galaxy population [547], LSBs with central surface brightness fainter than µ = 23 mag/ are experiencing a second wave of discovery after the initial discoveries in the 1990's enabled by CCDs [548,549,550,551,552]. With new telescope and survey designs, LSBs with surface brightnesses as faint as µ = 31 − 32 mag/ can be detected either as resolved stellar populations [553] or in diffuse light [554]. This population is heterogeneous. Many are "ordinary" dwarf galaxies, which are low surface brightness by virtue of the size-luminosity relation [555,556,557]. However, a number appear to be bigger outliers in either the size-luminosity relation or in the stellar mass-halo mass relation. Ultra-diffuse galaxies (UDGs) are an outlier population of LSBs with large sizes of r e > 1.5 kpc. Since 47 of these galaxies were found by the Dragonfly telescope in the Coma cluster [558], large populations have been observed in well-studied clusters, including Coma, Virgo, and Fornax [559,560,561,116,559,562]. In recent years, UDGs have been observed in a diversity of environments, including the field [563,550,564,551,552,565]. There are differences in the properties of UDGs in field and cluster environments-field galaxies tend to be gas-rich, blue, and star-forming, while those in clusters tend to be redder and quenched [550,561,562,566,567]. Moreover, a number of cluster UDGs have unusual properties, such as being surrounded by anomalously many globular clusters given their stellar mass, and inhabiting halos far larger than indicated by the mean stellar mass-halo mass relation [568,569,570,571,572,573,574]. Because of these differences in field and cluster environments, and because of the odd properties of cluster UDGs, an important question to answer is if these galaxies are true outliers in galaxy formation, or if they are part of a continuum of LSBs. It is important to understand if there are different formation mechanisms for the field and cluster populations [575], and how those Figure 19: Observed and simulated stellar masses and half light radii. From [163]. mechanisms arise in the context of CDM and SIDM galaxy formation models.
In the context of CDM simulations, separate formation pathways are identified for field and cluster UDGs, and SIDM may additionally play a role in both. In the field, although high spin was originally posited as the driver of UDG formation [576], the current conventional wisdom is that flybys and minor mergers drive the large size of M * ∼ 10 8 M UDGs [577,578]. Vigorous outflows can play a role too [579], and enhance the effect of flybys and mergers. However, the presence of a large reservoir of gas in the field UDGs makes the fly-by and concomitant tidal stripping unlikely explanations for field UDGs [567,580]. In addition, the rotation curves of these field UDGs present challenges for all models of galaxy formation (whether based on CDM or SIDM) making this one of the most interesting avenues to pursue in the coming years [581]. Selfinteractions of DM may play a role in UDG formation, by enhancing the distribution of galaxy sizes and rotation curve shapes at fixed stellar mass, as hinted at for the broader class of LSBs [48,50], but more work is clearly needed to flesh out the predictions of SIDM models for these outliers.
Cluster UDGs are shaped by their environment. CDM simulations show that although some cluster UDGs were previously UDGs in the field, many became UDGs after infall [163,582,583,584]. These simulations show that "typical" dwarf galaxies can become UDGs if they fall into clusters early and experience quenching (by ram-pressure stripping), fading of light as the stellar population ages [575], and tide-driven expansion of the stellar population. Intriguingly, Refs. [585,163] argued that galaxies in DM halos with large cores more easily become UDGs. Idealized simulation including DM self-interactions in Ref. [586] have shown that it is possible to reproduce the properties of NGC1052-DF2 and DF4 in the SIDM models.
The stellar component in dark-matter-dominated galaxies expands faster in cored halos rather than cusped halos in response to tides, which has been shown in the context of CDM [155,587] as well as SIDM [75]. In the semi-analytic model of Ref. [163], cored halos play the essential role in populating the tail encompassing large UDG sizes, as shown in Fig. 19. This model with cored halos at infall was able to reproduce the observed sizes and stellar masses of UDGs, as well as their abundance as a function of host halo mass. However, the assertion that cored halos are necessary has been questioned in Ref. [588] who argue that it is possible for a large number of UDGs to form in the field as LSBs and then be accreted into the cluster. In this scenario, there are still be some UDGs created close to the center of the cluster due to extreme tidal interactions in cuspy halos. More work on discriminating the two scenarios is warranted, as are cosmological CDM and SIDM simulations tailored to this problem.

Warping of stellar disks
We described in Sec. 3.6 how the drag force experienced in SIDM by halos falling through a medium of DM can cause a separation between the centroids of stars and DM in a galaxy. Both numerical [589] and simple analytic [590] models show that the potential gradient that this establishes across the stellar disk induces a U-shaped warp in the disk facing the direction of infall (Fig. 20, left), and a longer-lasting disk thickening. While S-shaped warps formed by tidal distortions of the stellar light profile are abundantly observed in cluster environments, indicating that they are readily produced by physics unrelated to DM self-interactions, such physics is unlikely to generate prominent U-shaped warps. Such warps are only formed by a differential force on the disk and its halo, e.g. due to SIDM drag.
The simulations of Ref. [589] indicate that the warping effect may be observable on a galaxyby-galaxy basis in current and next-generation photometric data when SIDM cross sections are 0.5 − 1 cm 2 g −1 . This has inspired work aimed at constraining SIDM by means of this signal [590]. Here, the r-band images of ∼3000 edge-on disk field galaxies with M * > 10 9 M and D < 250 Mpc from the Nasa Sloan Atlas are reduced to yield estimates of their degree of U-shaped warping, and a forward-model is developed to predict this in SIDM as a function of interaction cross section and range, halo properties, and DM environment. This is then combined with a simple model for astrophysical and measurement noise to create a galaxy-by-galaxy likelihood for the warp statistic, and the noise and other nuisance parameters are marginalised over to obtain constraints on the SIDM cross section. Due to the requirement that the system satisfies the fluid approximation and that drag rather than evaporation is the dominant result of the self-interactions, only long-range (light mediator) models can be constrained in this way, although a sample of galaxies in higher-density multi-streaming environments (e.g. clusters) may in the future be used to constrain shorter-range forces.
No correlation of the warp strength or direction with the SIDM expectation is found in the analysis, leading to the constraints shown in Fig. 20 (right) on the momentum transfer cross section at 300 km/s. These are plotted as a function of the median relative velocity of the galaxies and the background medium because this is considerably uncertain. The two models shown correspond to a) a constant velocity over the entire sample and b) velocities proportional to those from the CosmicFlows-3 catalogue, respectively. At face value these bounds are an order of magnitude tighter than those from other SIDM effects, e.g. dwarf galaxy evaporation. The model does however make a number of simplifying assumptions to render the inference tractable, for example that the halos are fully spherical, thermalised and equilibrated, and that all parts of the halo move together. Simulations of halos in environments similar to those of galaxies with measurable morphology will be needed to validate and refine this novel method for testing SIDM.

Large-scale structure constraints
As we have discussed in Sec. 3.7, there are several particle physics models that can modify the power spectrum of density fluctuations, thereby leaving signatures on the temperature and polarization of the CMB, the clustering of galaxies and the halo mass function. A necessary feature in these models is the presence of light mediators. As discussed in Section 3.7 dark matter models that tightly couple to dark radiation in the early universe can lead to dark acoustic oscillations and diffusion damping in dark matter (like the the BAO and Silk damping), modifying the linear power spectrum. In Ref. [244], cosmological data from the lensing of the CMB, the clustering of galaxies on linear scale and the BAO were analyzed to derive constraints on the fraction of dark matter that can couple to dark radiation and the strength of the coupling. Their results show that a fraction of partially interacting dark matter, f int > 5% is in severe tension with current observations of large scale structure.
The damping scale of the power spectrum due to interactions between dark matter and radiation is set by the horizon at the epoch of kinetic decoupling. For a kinetic decoupling temperature,  Figure 20: Left: Cartoon of the formation of a U-shaped stellar disk warp in SIDM. Right: Constraints on the momentum transfer cross section at 300 km/s for long-range (Rutherford) scattering from 3,213 measured warps in the nearby Universe, as a function of the median galaxy velocity relative to the background. Adapted from Ref. [590].
T kd , the suppression in the power spectrum occurs approximately for modes k > 10(T kd /keV) Mpc −1 [591,592,593,8,192]. SIDM coupling to a dark radiation component like dark photons or sterile neutrinos can significantly delay the decoupling of dark matter from the thermal bath leading to damping on scales of order 0.1 Mpc/h or smaller. An effective framework to describe these interactions and the associated phenomenology has been developed by the ETHOS collaboration [594,179]. One of the central predictions is a potentially observable cut-off in the halo mass function [233,234,235,251,245,5] arising from the damping of the power spectrum.
One of the ways to probe these scales is through the absorption spectrum of Lyman-α forest, which measures the density of neutral hydrogen on large scales. In Ref. [5] the Lyman-α forest power spectrum measurements was used to put strong constraints on the damping scale (by recasting the warm dark matter constraints) to set an upper limit on the decoupling temperature of T kd > 1 keV and a minimum halo mass of M min < 10 8 M /h. The origin of this constraint is rooted firmly in the physics of the dark sector, including the production mechanism for the dark matter. The massive gauge boson that mediates the self-interaction in simple SIDM models has to decay in order to avoid over-closing the universe and in Ref. [5] decay of the gauge boson to light DM particles were assumed. Decays to Standard Model particles have much stronger direct detection constraints [595] and indirect search constraints [4]. In the future, we expect the combined measurements of the Lyman-α power spectrum and the halo mass function from faint galaxies will strongly constrain the space of allowed parameter space for the models with dark radiation. Prospects for the future in this direction are bright, with CMB-S4 [596] 12 combined with LSST likely to improve constraints on ∆N eff in these models by an order of magnitude [231]. Another promising avenue to constrain simple SIDM models and those with additional light degrees of freedom and interactions is to probe the linear and quasi-linear power spectrum to higher k values using 21-cm data. Future measurements offer the possibility to constrain the power spectrum amplitude at the 10% level integrated over k = 40 − 80 Mpc −1 [597], which can probe decoupling temperatures well above a keV.

Summary of constraints on the self-interaction cross section
The various observational probes for SIDM can be thought of as a mapping between the velocity dependence of the cross section and observable processes at various energy scales within halos of 12 https://cmb-s4.org/ different masses during hierarchical structure formation. On the galaxy cluster scales, where the internal velocity dispersions are of the order of a few 1000s of km/s, the most stringent constraints currently come from measured core densities from strong lensing [55]. The small sizes of the cores and the high central densities in these clusters constrain cross section to σ/m DM < 0.13 cm 2 /g at 95% confidence in one study [353] and σ/m DM < 0.35 cm 2 /g at 95% confidence in another [354], as discussed in more detail in Sec. 4.1.1. Even though there is some evidence of a non-zero cross section in these analyses, it is not possible to conclude that we have a detection. In the future, detailed simulation work combined with joint analyses of strong lensing, weak lensing, X-ray and velocity dispersion data for a large number of clusters could provide a concrete path to measuring a non-zero cross section at relative velocities of about 1500 cm 2 /g. Similar analysis at the group scale gives an upper limit of 1.1 cm 2 /g, with a mild preference for a non-zero cross section (0.5 ± 0.2 cm 2 /g; [354]). A variety of other probes like the displacement between galaxies and dark matter in merging systems, the steepening of the profiles beyond the dark matter core and the existence of displaced BCGs will further help establish consistency between different observations once the measurement systematics and the interplay between baryonic and dark matter evolution are better constrained and understood. Among these, constraining BCG wobbles, which are correlated with the size of the constant density core, using offset BCGs seems highly promising with current constraints at σ/m DM < 0.39 cm 2 /g at the 95% C.L. [428].
At lower velocities, a wider range of cross sections is in principle still allowed as long as the model is consistent with the constraints at high velocities. As baryonic effects are more relevant for lower mass halos, the combination of baryons and SIDM leads to a wide range of halo mass profiles inside a region of about the NFW scale radius. In particular, for galaxies like the MW and large elliptical galaxies, the large concentration of baryons in the inner regions makes the density profiles of SIDM and CDM models indistinguishable [119,277,598,599]. We can expect constraints from evaporation of DM in satellites of the MW (and similar galaxies) due to scattering events between subhalo and host halo DM particles when the cross section at 200 km/s exceeds 5 − −10cm 2 /g (e.g., [153,276,470]) but the constraints have not been robustly inferred. At lower velocities, fitting the rotation curves in SPARC database has shown that σ/m DM > 3 cm 2 /g at about 100km/s is needed to adequately explain the low-surface brightness galaxies with the lowest central DM densities [50]. Note that the predictions for the cored profiles of SIDM halos do not change dramatically as the cross section is increased, so cross sections even an order of magnitude larger cannot be easily distinguished based on galaxy rotation curve data alone [117]. If the cross sections are 100 cm 2 /g, then core collapse will set in for median concentrations and it may become difficult to explain the low density DM cores of low-surface brightness galaxies.
Within the local universe, a self-consistent explanation for the distribution and inner structure of the satellites of the Milky way is not yet fully known. While the "missing satellites problem" is largely thought to be resolved given the current census of satellites, the exact nature of the galaxy-subhalo connection is still unclear and currently requires populating galaxies in subhalos with V peak ∼ 7 km/s that are below the collisional cooling limit. Recent work has shown that the minimum V peak will have to be smaller for moderate cross sections σ/m dm 5 cm 2 /g [472]. While the "Too Big To Fail" problem is resolved by moderate SIDM cross sections of ∼1 cm 2 /g [600,277], this underpredicts the central densities of ultra-faint satellites [601,51,276]. These conclusions are highly sensitive to the presence of the MW disk, but recent work including the disk potential has demonstrated the validity of these arguments [472]. The issue is not restricted to the ultra-faints dwarfs-using analytic methods it has been found that if all the bright Milky Way satellites are in the core expansion phase, it is difficult to fit their measured stellar dispersions with a single cross section [47]. A velocity-dependent cross section with σ/m DM > 10 cm 2 /g at velocities relevant for the internal dynamics of satellites (5 − 20 km/s) that drives the faint satellites to core-collapse is perhaps the best bet to explain these observations [103,276,51,471,470,472].
Putting together these constraints suggests that a significant velocity dependence is necessary to explain data constraining DM mass profiles from dwarf galaxy to clusters of galaxies. Compared to earlier work [55], the motivations for a strong velocity dependence have increased over time. An exciting possibility within models that feature large cross sections at low velocities is that many of the halos and subhalos could be in the core collapse phase of their evolution, with predictions for observables that are different from CDM as discussed in this review. Examples of such sharp velocity dependence are discussed in Refs. [471,130,470], where analytic methods and N-body simulations are used to predict the central densities of Milky Way satellites in the core expansion and core collapse phases. More comprehensive analytic and numerical treatments including the effect of the MW disk are required to make further progress. However, it seems clear that the satellites of the MW can provide a definite test of the SIDM models that aim to explain DM density profiles across the mass range from dwarf galaxies to galaxy clusters.

Beyond the simple SIDM framework
In the basic setup of SIDM, a single species of DM particles interact through the exchange of a light (or massless) scalar or vector mediator with a cross sections per unit mass of O(1 cm 2 /g), as discussed in Sec. 2. We have, however, largely neglected the implications of a hidden sector beyond the phenomenology of the DM particle itself. There are also more complex hidden sectors that predict interesting phenomenological signatures, not only in astrophysical systems, but in cosmology and potentially in laboratory experiments as well. While there are model-dependent constraints on any specific particle physics realization of SIDM, many models may produce similar phenomenological signatures and are thus subject to similar restrictions from data. In this section, we review the consequences of non-minimal SIDM models and of the hidden sectors themselves. We first discuss two broad generalizations of the simple SIDM scenario that we have assumed thus far. In Sec. 5.1, we relax the assumption that self interactions are elastic and allow a mechanism for SIDM to dissipate energy. In Sec. 5.2, we discuss the possibility that SIDM as a subcomponent of the overall DM content, allowing it to have very strong self interactions. The hidden sector contains light mediators that allow DM particles to self interact, and we discuss the consequences of the light mediator in Sec. 5.3. Finally, in Sec. 5.4, we mention various complementary constraints on DM sectors that interact with the SM.

SIDM with dissipation
If SIDM is not a simple single-state particle, it is possible for SIDM to undergo inelastic processes between its ground state and an excited state. The presence of exicited states can arise in models of composite SIDM, such as atomic-like DM [192,86] and nuclear-like DM [85], in which fundamental dark-sector particles are bound under a dark-sector force. The composite nature of DM permits a spectrum of internal energy levels, along with a naturally large self-interaction cross section. Excited states also arise if DM is comprised of different fundamental particles with a mass splitting and an off-diagonal interaction [602,603,604,605]. In either case, collisions may induce a transition to the excited state, which may then relax back to the ground state through collisional de-excitation or decay via the emission of some form of radiation. This radiation may be a massless force carrier, such as a dark photon, or a light mediating particle, with a mass smaller than the energy splitting between the ground and excited states.
Assuming there is sufficient kinetic energy available, collisions within DM halos can excite DM particles. If the decay width of the excited state is large and if the halo is optically thin to the emitted radiation, the halo loses energy and could eventually collapse [86]. As discussed in Sec. 3.2, the gravothermal collapse process with elastic SIDM interactions can be accelerated by environmental effects [103]; for an isolated halo, the collapse process can be accelerated through energy dissipation [118], assuming DM particles are not ejected from the halo [606]. The observation of uncollapsed dwarf galaxies and LSBs can place restrictions on the acceptable amount of halo cooling in those systems [118]. On the other hand, the gravothermal collapse of DM halos with dissipation may be able to form the seeds for the observed population of supermassive black holes [607].

Subcomponent SIDM
Some models of DM postulate that while most of the DM in the universe may be (approximately) collisionless, there may be a small component that has very large self-interactions, including dissipative interactions. This small component could have a large impact on the observational predictions. For example, even a very small amount of ultra-strong SIDM can cause a halo to experience gravothermal collapse, as discussed in Sec. 3.2, and form seed black holes [85,102]. Alternatively, some models also propose that that up to 15% of the DM of the Universe exhibit a behavior similar to baryons (i.e., experience energy dissipation), including disk collapse [195,194] and plasma instabilities [608,609,610]. These effects have for example been studied in the context of the formation of a dark core in major mergers such as Abell 520 [424,611]. Dissipation could also lead to qualitatively new predictions for the small-scale structure of DM, including asymmetric dark stars [612].

Phenomenology of light mediators
In many models of SIDM, the self interaction arises from the exchange of a light mediator, such as a 1-100 MeV scalar or vector particle with a Yukawa coupling to DM [81]. The mediating particles themselves can impact a wide range of observations, allowing for additional probes of SIDM hidden sectors.
The production mechanism of SIDM in the early Universe can play an important role throughout cosmic evolution (see, e.g., Ref. [613]). A common assumption is that DM particles were at some point in thermal equilibrium with the bath of CMB photons, through a common reheating mechanism. In a completely isolated hidden sector, however, DM may be produced in its own thermal bath through a separate reheating mechanism, resulting in isocurvature perturbations that are constrained by CMB observations [614], and thus place constraints on self interactions [615].
With a nontrivial dark sector, the SIDM relic abundance can be set through mechanisms analogous to the standard thermal freeze-out scenario. For instance, SIDM may undergo dark-sector freeze-out (in which SIDM particles annihilate into lighter dark-sector particles) or the SIMP mechanism (in which SIDM particles experience number-changing processes).
The mediator constitutes the thermal bath for DM and contributes an additional ∆N eff to the standard number of relativistic degrees of freedom in the early Universe. Measurements of the CMB and of nuclei abundances restrict the amount of ∆N eff permitted at the time of recombination and Big Bang Nucleosynthesis (BBN), respectively [1]. Moreover, for models of light DM (with masses in the MeV range), the DM itself can contribute to ∆N eff during BBN [616]. Since no interaction with the SM is required in a minimal model of SIDM, the hidden-sector bath can be slightly colder than the CMB and thus easily evade ∆N eff constraints [617,85]. However, if the hidden sector does communicate with the SM, ∆N eff will set constraints on such models.
The presence of a hidden thermal bath for SIDM also leads to the suppression of small-scale fluctuations in the CMB, as described in Sec. 3 for the case of dark photon mediators. Such signatures can be constrained by the CMB, as well as large-scale structure observables, such as the Lyman-α forest [251,5,231] as discussed previously in Sec. 4.10. If models predict a strong cutoff of the matter power spectrum, then additional powerful constraints at small scales arise from the existence of the smallest DM halo observed [221].

Connections to the Standard Model
While not strictly required, there are good reasons to consider non-gravitational interactions of SIDM with SM particles. These interactions allow SIDM to be produced from the SM thermal bath, which can constrain the model space stringently [4] and also provide simple viable models [613,618,619]. These interactions can impact the composition and cosmological evolution of the Universe and may give rise to observational signatures in laboratory experiments and astrophysical systems (see e.g., Refs. [620,621,622]).
While there are many different ways in which SIDM particles can couple to the SM, most attention has been paid to two benchmark scenarios, in which the couplings of SIDM particles are proportional to the electric charge and mass of the various SM particles, respectively [623]. The former case arises from the (kinetic) mixing between a new vector boson and the photon of the SM [182], while the latter case arises from the mixing between a new scalar boson and the Higgs boson of the SM [624]. The new exchange particle (often referred to as the mediator) is typically considered to also be responsible for the generation of strong self-interactions.
It may be possible to understand the cosmological abundance of SIDM particles in terms of their annihilation and scattering rates [613]. If these rates are large enough to bring SIDM particles into thermal equilibrium with the SM plasma, the relic abundance can be set through the standard WIMP freeze-out mechanism. Even if the interaction rates do not bring the two sectors into equilibrium, the relic abundance may be set, for example, via the freeze-in mechanism [3]. Moreover, many SIDM models require the presence of additional (lighter) states in the dark sector, as discussed in Sec. 5.3. In order to ensure that these additional particles do not constitute a large fraction of DM or dark radiation, they typically need to be unstable, such that they decay into SM particles. In order to evade constraints on ∆N eff from BBN and the CMB and in order not to modify the predictions of BBN through the photodisintegration of light elements, these decays should occur sufficiently early in cosmological history, placing a lower bound on the interaction strength between SIDM particles and SM states [625].
On the other hand, cosmological considerations also place an upper bound on the interaction strength for DM annihilation into and DM scattering with SM particles. DM annihilation can inject energy into the SM plasma, affecting the optical depth of CMB photons post-recombination [626,627,628,629,630,631,4,632,633,634]. CMB constraints on DM annihilation are strong for s-wave annihilation, but significantly weaken for p-wave etc. due to the velocity suppression of cold DM. Additionally, there are constraints on DM scattering with ordinary matter, which can induce spectral distortions [215,223] and a suppression of anisotropies in the CMB [213,208,214,216,217,218,219,220,224,635,226] and the matter fluctuations [221,57,222,225]. These constraints are complementary to the more traditional search methods of direct and indirect detection, described below.
Direct detection experiments place strong constraints on the interactions of DM particles with the SM and are known to exclude many DM models that exploit the freeze-out mechanism to explain the DM relic abundance. Another concrete way in which the direct detection of a SIDM particle can be motivated is the requirement in some models that the mediators decay before the beginning of BBN [636]. These experiments are particularly sensitive if the mediator of the interaction is light, such that scattering with small momentum transfer is enhanced [637,638,595]. In many models of SIDM, the DM particle mass is below 1 GeV, for which the most promising strategy is to search for DM-electron interactions [639].
A second source for observable signals of SIDM particles are annihilation processes occurring in regions of enhanced DM density, such as MW satellites or the Galactic Centre. The resulting fluxes of photons, neutrinos and/or changed anti-particles can be searched for with indirect detection experiments. These constraints are particularly strong if the annihilation cross section in enhanced for small relative velocities (so-called Sommerfeld enhancement [168,640,9,641,642,10]). We note, however, that annihilation signals may also be strongly suppressed or absent, for example if the annihilation cross section is suppressed in the non-relativistic limit, or if SIDM particles have a particle-antiparticle asymmetry similar to the one of baryons [634].
Finally we note that SIDM particles may also be produced at accelerators, although the resulting signatures are expected to look very similar to other DM models, such as WIMPs. A notable exception occurs in the case of the SIDM particles produced form a bound state (called WIMPonium or Darkonium), which subsequently decays into SM particles [643,644,645,646,647,7] 13 will allow us to observe the faint satellite galaxies that inhabit the DM halo of the MW, allowing a follow-up similar to the discovery of satellites and streams that DES, among other surveys, has already made possible. Apart from observing the satellites of the MW, there are several surveys that target the faint galaxies in the nearby universe, among them the SAGA survey [437,444] which searches for satellites around MW analogs, ELVES [649] which provides a volume-limited sample of classical satellites around MW analogs with confirmed distances, LBT-SONG [650] which searches for satellites around intermediate-mass hosts beyond the Local Group, and MADCASH [651] which observes Magellanic analogs in the Local Volume. Apart from focussing on satellite systems, surveys like DELVE [652]-a deep, multi-component survey that uses the Dark Energy Camera to study low-surface brightness, DM dominated galaxies in the local universe-will observe field dwarf galaxies. Similarly the Merian survey [653] 14 uses the HSC to measure weak lensing profiles of Small and Large Magellanic Cloud analogs in the local universe. Spectroscopic follow-up of these galaxies will be enabled by Keck, the Thirty Meter Telescope and proposed experiments like the Maunakea Spectroscopic Explorer (MSE; [654]) allowing us to characterize the detailed properties of these faint objects. Furthermore, surveys like the Widefield ASKAP Lband Legacy All-sky Blind surveY (WALLABY) [655]-an HI survey in the local universe-will illuminate the baryonic component of faint galaxies in the nearby universe.
Within the next decade, the Rubin Observatory and Euclid are expected to discover thousands of new galaxy-scale strong gravitational lenses [656]. In order to exploit the future data to the maximum possible extent, new methods have been developed to accelerate detections of such systems as well as perturbers in the images [657,658,659,660,661,662]. The vast increase in sample sizes will provide much stronger constraints on DM scenarios. Similar machine learningbased methods may also facilitate the analysis of other data from upcoming surveys.
A dedicated follow-up effort for a well-chosen subsample of these lens systems to obtain spectroscopic redshifts and (for Rubin-discovered systems) high-resolution imaging will enable improved constraints on SIDM. On the one hand, gravitational imaging analyses of extended arcs and rings are likely to lead to the detection of several DM perturbers (either substructure or lineof-sight halos) [360,361,362,364]. Detailed modeling of the density profiles of these perturbers [369,370] will provide constraints on the SIDM cross section at small velocities. On the other hands, a joint analysis of a large population of quasar lenses [384,378] could also be used to constrain the SIDM cross section. In parallel, follow-up ALMA observations of lensed dusty starforming galaxies discovered in South Pole Telescope (SPT; [663]) 15 maps [664,363] and Very Long Baseline Interferometry (VLBI) observations [665,666] of lensed radio sources will also be used to discover small DM halos and subhalos. In both cases, detailed studies of the density profile of these perturbers will lead to constraints on SIDM. As we have discussed, the possibility of the detection of a number of perturbers that are denser (or more compact) than expected in CDM can be a clean path to the discovery of dark sector DM. Finally, high-resolution imaging of galaxy-scale lenses could be used to measure the small-scale power spectrum [392,388,379,389,390,359] independent of how of the density profile of line-of-sight halos and subhalos are modeled, hence providing complementary constraints on SIDM.
The upcoming galaxy surveys across a wide range of halo masses, including both satellite and field galaxies, will help constrain the cross section on multiple velocity scales giving us a rich understanding of the varied phenomenology that derives from the interplay of hierarchical structure formation and the particle properties of DM. While observations discussed above will help chart out the nearby Universe, we note also that data from the James Webb Space Telescope [667] 16 and SphereX [668] 17 in combination with the Rubin Observatory data will give us novel constraints on the nature of DM from the first galaxies at high redshifts from the epoch of reionization.
In addition, constraints on SIDM cross sections from massive halos, cluster and group profiles inferred via weak lensing [95] and core sizes [669,354,353] determined via strong lensing and stellar velocity are expected to tighten significantly given the number of such massive clusters (M 1 × 10 14 M /h) that will be detected by a combination of optical surveys-such as Rubin and Euclid, as well as Y6 data from DES-and CMB surveys such as the SPT, Atacama Cosmology Telescope (ACT; [670]) 18 , Simons Observatory [671] 19 and CMB-S4. Given its sky coverage, Rubin alone is expected to detect more that 50, 000 clusters in this mass range, up from ∼ 5000 clusters detected in the DES Y1 data. At slightly lower mass scales (5 × 10 13 M /h M 10 14 M /h), the eROSITA mission [672] 20 is expected to map out ∼ 100, 000 low-mass clusters and galaxy groups. On those parts of the sky where eROSITA overlaps with optical surveys, weak lensing measurements around these objects will improve constraints on SIDM cross sections at the relevant velocity scale. Moreover, with an increasing sample of clusters and precise weak lensing maps from telescopes like Roman and Euclid cross section constraints from offsets between BCGs and potential centers of galaxy clusters are also expected to improve by a factor of two [428]. Experiments like MSE will allow complementary studies of positional offsets with velocity offsets.
A key advance in the future will be robustly combining weak lensing, strong lensing, X-ray and stellar velocity dispersion data with a dedicated analysis that includes a self-consistent halo model for SIDM. It is possible that we may be able to measure cross sections of order 0.1 cm 2 /g with the next generation of observations and analysis tools. A necessary requirement for this to become a reality is the availability of cosmological simulations of galaxy clusters using cross sections of 0.1 cm 2 /g and smaller. We discuss simulations in more detail in the next subsection.
Finally, the upcoming surveys, several of which are optimized to answering questions about the large-scale structure and evolution of the universe, will provide unprecedented data sets to test models that predict modifications to the power spectrum (Secs. 3.7, 4.10 and 5.3). Conservatively, CMB-S4 for example is expected to measure the energy density in weakly-coupled light particles, parametrized by ∆N eff , to a precision of σ(N eff ) = 0.02 − 0.03 [596]. As we have discussed, ∆N eff is a critical discriminant of SIDM models and dark sector physics in general. These experiments will also be contemporaneous with galaxy surveys from the Rubin Observatory, Euclid, DESI and the Roman Space Telescope, which will measure cosmological parameters from the clustering of galaxies, gravitational lensing and the abundance of clusters. Simons Observatory and CMB-S4 will also improve CMB lensing measurements by an order of magnitude compared to Planck [596]. Combined with measurements of Lyman-α forests from spectrographs like Prime Focus Spectrograph on the Subaru telescope [673], WEAVE on the William Herschel telescope [674], and in the future MOSAIC on the Extremely Large Telescope [675], WFOS on TMT [676] and MSE [654], we will have an opportunity to combine several different measurements to constrain these models of dark matter with light mediators. Moreover, 21cm surveys like HERA, SKA and LOFAR along with the James Webb Space Telescope will allow us to push to the epoch of reionization and measure the small-scale power spectrum [597]. Combined with the small-scale power spectrum, the mass function of dark matter halos and progress in understanding galaxy formation over a wide range of redshifts, we can be optimistic about making significant strides in pinning down the microphysical nature of dark matter in the near future.

Simulations and Modeling
With the wealth of data that will be available in the coming years, there is a need for simultaneous advancement in simulation and modeling methods of alternate DM models in order to use this data to constrain properties of DM. For SIDM models of the type considered in this work, much of the progress has been based on DM-only (i.e. without baryons) modified N -body simulations. Different groups involved in the simulation efforts broadly agree on the types of phenomena produced by the presence of elastic self-interactions, as well as the size of these effects as a function of interaction cross section. On the other hand, a number of powerful probes of SIDM presented in this work involve modeling of small-scale structure, especially the evolution of substructure inside a larger virialized object, where the different simulations techniques have not been rigorously tested against each other. For collisionless ΛCDM models, different implementations of N -body simulations have been stringently calibrated against each other over the years, ensuring that numerical artifacts from any specific implementation can be detected and corrected [677,678] (also see the AGORA code comparison project for hydrodynamic simulation algorithms [679,680,681]). A code comparison project of this type for different SIDM implementations will be hugely beneficial toward the goal of predicting the effects of SIDM on small-scale structure robustly (cf. [682]). Furthermore, since the small-scale structure of interest often lie close to the resolution limits of these simulations, a thorough understanding of whether the resolution affects the predictions of SIDM effects, even within individual implementations, is needed.
In recent years, there has been a growing focus on the mechanism of disruption of subhalos in collisionless CDM N -body simulations (see, e.g., [683,684,685,466,686]). It has been pointed out that much finer force resolution may be needed to faithfully follow the tidal stripping of subhalos inside a more massive host halo than what was estimated previously, and has been used in much of the existing literature. Lack of sufficient force resolution results in larger amounts of tidal stripping than is expected physically. A related issue is that of subhalos being stripped to below the mass resolution of the applied halo-finder and disappearing from the substructure catalog, even though the central cusp retains its identity, and can continue to host a luminous galaxy. While some of these issues can be overcome by additional modeling, e.g. including an "orphan model" for galaxies living in the substructure that fall below the halo finder detection threshold, these numerical effects can have major implications on predictions of strong lensing, satellite weak lensing signals and the phase space distributions of satellites in the nearby Universe. An initial exploration of this effect was done in [687], showing that accounting for disrupted orphans is important to model satellite weak lensing in clusters and to understand the overall distribution of subhalos. Additionally, it should be noted that tracking orphans in SIDM simulations is nontrivial as interactions can often scatter particles near the core of a subhalo, making methods like core-tracking inherently difficult. A detailed understanding of the effects of artificial disruption in SIDM N -body simulations is essential, especially as they can be degenerate with certain physical signatures of DM self-interactions. Ram pressure stripping of subhalos in SIDM simulations represent an additional mass-loss mechanism, and the presence of a central core in SIDM subhalos instead of a cusp can mean that the substructure can be completely destroyed.
Most algorithms that are currently employed to implement SIDM in N -body simulations are designed explicitly to work in the limit of rare self-interactions, i.e. when an individual simulation particle interacts via at most one other particle through self-interactions in a given time-step. This typically also corresponds to the regime where most particles do not scatter in any particular time-step. On the other hand, certain SIDM models have been proposed where it is necessary to model frequent self-interactions, i.e. where a single particle can, in principle, interact with multiple other particles in a single time step. The current simulation methods, therefore, need to be updated to correctly predict the phenomenology of these models, a direction in which there has already been some progress [114,688]. The development of these algorithms could also help robustly model the onset and subsequent evolution of the core-collapse phase of SIDM in N -body simulations, where interactions become progressively more frequent in the collapsing region and its outskirts. Since the core-collapse mechanism is potentially interesting for a whole range of observed phenomenology of dwarf satellites, advances in simulation techniques will have large science impact. In this context, it is also worth investigating whether N -body techniques are indeed the most appropriate tools to model this regime, or whether gravothermal fluid techniques in these regions can perform better.
DM-only simulations, of course, make predictions for the evolution of DM halos and subhalos. On the other hand, cosmological surveys generally measure the positions, and in some cases, velocities of luminous galaxies. The connection between the observed galaxies and the DM halos in which they live needs to be effectively modeled in order for the simulations to be used in the interpretation and analysis of observational data. In the context of CDM, the galaxy-halo connection has been studied extensively (see [33] for a recent review), and various effective parameterizations have been developed to marginalize over the unknown galaxy formation physics, and its potential impact back-reacting onto the dynamics and evolution of DM substructure. Such studies are also needed in the context of SIDM to promote results from DM-only SIDM simulations to analysis tools. An important question to answer in this context is whether the same parameterizations developed in the context of CDM can be applied to the SIDM context, or whether the galaxy-halo connection in SIDM is qualitatively different, and needs a completely different approach.
Along with a detailed scrutiny of cosmological simulations of DM it is also essential to understand the varied aspects of the interplay between baryonic physics and self-interacting DM over a wide range of halo mass scales. Ideally one would like to simultaneously explore the parameter space in SIDM along with varying the feedback mechanisms in hydrodynamic simulations.
At the very least, when using simulations to investigate how observable quantities are affected by the DM model one should compare this with the differences seen with a fixed DM model when different implementations of baryonic physics are used.
Another avenue for future code development is to include a larger number of possible DM interactions. To date, the bulk of work on SIDM has considered elastic scattering, typically with an isotropic differential cross section. There has been a small amount of work that considered anisotropic scattering cross sections [49,95], or inelastic scattering [606], but these scenarios have not been exhaustively explored and do not cover the full range of phenomenology that DM with self-interactions could possess.

Promising tests of simple SIDM models
We have described a variety of tests of SIDM models in this review, ranging from the least massive dark subhalos that could perturb lensed images to galaxy clusters and large-scale structure. It is fascinating that the simple modification of including elastic self-interactions to the CDM model, which is well-motivated from the particle physics perspective, could impact structure formation in so many different ways. It is also sobering to realize that despite these striking predictions, the remarkable progress in simulating galaxies, and the plethora of data that already exists, there are no definitive findings in favor of or against SIDM. The reasons include the fact that CDM models are an excellent fit to the majority of data, and where it may have trouble either star-formation physics (with large modeling uncertainties) is important or the data is not yet decisive. SIDM models predict deviations from the CDM model predictions predominantly in the central regions of galaxies, where baryonic physics is important. This is a feature of simple SIDM models and not a fine-tuning, as we have reviewed in Sec. 3.
We have identified multiple avenues to test SIDM through astrophysical observables. A key prediction of all DM models is the mass function of dark subhalos and their internal density pro-files. Strong lensing studies aimed at testing these predictions will mature in the coming years, as we have discussed. Observables associated with resolved and unresolved subhalos and their power spectrum should provide a definitive test of a wide swathe of SIDM models that predict halos and subhalos are mostly in the core-expansion phase of their evolution. We have also discussed the recent progress in fleshing out the predictions of SIDM models that predict some subhalos could be in the core-collapse phase. These subhalos with their steep density profiles in the inner regions (radii smaller than the NFW scale radius) provide a striking target for future observations of strongly-lensed systems and stellar streams. Much more work remains to be done to work out the predictions for these models and the critical correlations between subhalo's density profile, survival and mass loss due to stripping and evaporation.
We also highlighted the possibility of throwing the full power of lensing (strong and weak) reconstructions, X-ray observations, BCG offsets, mergers (major and minor) and stellar kinematics in groups and clusters of galaxies to put definitive constraints on the cross sections at relative velocities of 1000 − 2000 km/s. It is only possible to put upper limits currently but the massive influx of data in the coming years, along with the possibility of improvements in theory and modeling, makes this an exciting avenue to search for a positive detection.
The most concrete astrophysical motivations for SIDM have been the core-cusp and diversity issues of rotation curves, and the too-big-to-fail problem of the Milky Way satellites. We have reviewed these problems in detail along with the constraints imposed by the growing census of Milky Way satellites. In the near future, there will be an explosion of data on galaxies and satellites at low redshift, which will provide critical tests of SIDM models. In this context we have highlighted the unique opportunities provided by the ultra-faint dwarf spheroidals of the Milky Way, other dwarfs in the local volume and the ultra-diffuse galaxies in the field and clusters. It is safe to say that there is not currently a concrete understanding of the dark matter or stellar distribution of these objects in the context of either CDM or SIDM models. As with the dark subhalos, the prevalence of core expansion vs core collapse phases in the satellites will be critical. It seems clear that SIDM models where interactions are not strong enough to push satellites into the core collapse phase will be definitively tested in the near future.
SIDM models provide a viable, predictive and versatile framework for interpreting astrophysical observables on galactic and sub-galactic scales.
The predictions of these models need to be developed further if we are to maximize the science returns from the amazing data expected this decade and next. SIDM models provide a unique and compelling opportunity to detect the particle nature of DM through astrophysical observables, as we have highlighted throughout this review. They also provide a well-defined and predictive foil to the dominant CDM paradigm; with predictions for simple SIDM models in hand, the tests of the CDM model become clearer. However, progress will only be possible if there is a serious investment of resources in simulating SIDM models. We have discussed many aspects of the required simulations throughout this review. This investment in simulations is critical given the incomplete progress in mapping out the small-scale structure of even the simplest (single-component, elasticscattering-only) SIDM models. It is also vital because progress in understanding galaxy formation and the phenomenology of DM self-interaction are intertwined. We need a concerted effort on the theory side to take full advantage of the incredible data that we will have in the near future and develop the astrophysical observables we have reviewed into decisive probes of dark sector physics.