Dark Matter Self-interactions and Small Scale Structure

TL;DR

This paper reviews self-interacting dark matter theories, their astrophysical implications, and observational constraints, highlighting how SIDM can address small-scale structure issues while remaining consistent with large-scale cosmology.

Contribution

It provides a comprehensive overview of SIDM models, simulations, observational probes, and constraints, offering insights into their role in solving small-scale structure problems.

Findings

SIDM can explain core-cusp and diversity problems in galaxies.

Simulations show SIDM affects density profiles and halo shapes.

Current constraints limit the strength of dark matter self-interactions.

Abstract

We review theories of dark matter (DM) beyond the collisionless paradigm, known as self-interacting dark matter (SIDM), and their observable implications for astrophysical structure in the Universe. Self-interactions are motivated, in part, due to the potential to explain long-standing (and more recent) small scale structure observations that are in tension with collisionless cold DM (CDM) predictions. Simple particle physics models for SIDM can provide a universal explanation for these observations across a wide range of mass scales spanning dwarf galaxies, low and high surface brightness spiral galaxies, and clusters of galaxies. At the same time, SIDM leaves intact the success of $\Lambda$CDM cosmology on large scales. This report covers the following topics: (1) small scale structure issues, including the core-cusp problem, the diversity problem for rotation curves, the missing satellites problem, and the too-big-to-fail problem, as well as recent progress in hydrodynamical simulations of galaxy formation; (2) N-body simulations for SIDM, including implications for density profiles, halo shapes, substructure, and the interplay between baryons and self-interactions; (3) semi-analytic Jeans-based methods that provide a complementary approach for connecting particle models with observations; (4) merging systems, such as cluster mergers (e.g., the Bullet Cluster) and minor infalls, along with recent simulation results for mergers; (5) particle physics models, including light mediator models and composite DM models; and (6) complementary probes for SIDM, including indirect and direct detection experiments, particle collider searches, and cosmological observations. We provide a summary and critical look for all current constraints on DM self-interactions and an outline for future directions.

Tables

Table 1. Table 1: Summary of positive observations and constraints on self-interaction cross section per DM mass. Italicized observations are based on single individual systems , while the rest are derived from sets of multiple systems. Limits quoted, which assume constant σ / m 𝜎 𝑚 \sigma/m , may be interpreted as a function of collisional velocity v rel subscript 𝑣 rel v_{\rm rel} provided σ / m 𝜎 𝑚 \sigma/m is not steeply velocity-dependent. References noted here are limited to those containing quoted self-interaction cross section values. Further references, including original studies of observations, are cited in the corresponding sections below.

Positive observations	$\sigma /m$	$v_{\mathrm{rel}}$	Observation	Refs.
Cores in spiral galaxies	$\gtrsim 1{\mathrm{cm}}^2/g$	$30-200\mathrm{km}/s$	Rotation curves	Dave et al. (2001); Kaplinghat et al. (2016)
(dwarf/LSB galaxies)
Too-big-to-fail problem
Milky Way	$\gtrsim 0.6{\mathrm{cm}}^2/g$	$50\mathrm{km}/s$	Stellar dispersion	Zavala et al. (2013)
Local Group	$\gtrsim 0.5{\mathrm{cm}}^2/g$	$50\mathrm{km}/s$	Stellar dispersion	Elbert et al. (2015)
Cores in clusters	$\sim 0.1{\mathrm{cm}}^2/g$	$1500\mathrm{km}/s$	Stellar dispersion, lensing	Kaplinghat et al. (2016); Elbert et al. (2016)
Abell 3827 subhalo merger	$\sim 1.5{\mathrm{cm}}^2/g$	$1500\mathrm{km}/s$	DM-galaxy offset	Kahlhoefer et al. (2015)
Abell 520 cluster merger	$\sim 1{\mathrm{cm}}^2/g$	$2000-3000\mathrm{km}/s$	DM-galaxy offset	Jee et al. (2014a); Kahlhoefer et al. (2014); Sepp et al. (2016)
Constraints
Halo shapes/ellipticity	$\lesssim 1{\mathrm{cm}}^2/g$	$1300\mathrm{km}/s$	Cluster lensing surveys	Peter et al. (2013)
Substructure mergers	$\lesssim 2{\mathrm{cm}}^2/g$	$\sim 500-4000\mathrm{km}/s$	DM-galaxy offset	Harvey et al. (2015); Wittman et al. (2017)
Merging clusters	$\lesssim \mathrm{few}{\mathrm{cm}}^2/g$	$2000-4000\mathrm{km}/s$	Post-merger halo survival	Table 2
			(Scattering depth )
Bullet Cluster	$\lesssim 0.7{\mathrm{cm}}^2/g$	$4000\mathrm{km}/s$	Mass-to-light ratio	Randall et al. (2008)

Table 2. Table 2: Summary of merging cluster constraints on SIDM. All values for σ / m 𝜎 𝑚 \sigma/m are upper limits except for Abell 520 and 3827. For Abell 520, two quoted values of σ / m 𝜎 𝑚 \sigma/m are obtained under different assumptions in which τ = Σ dm ​ σ / m ≈ 0.25 𝜏 subscript Σ dm 𝜎 𝑚 0.25 \tau=\Sigma_{\rm dm}\sigma/m\approx 0.25   Mahdavi et al. ( 2007 ) or 0.13 0.13 0.13   Jee et al. ( 2014a ) .

Cluster	$\sigma /m$	Method used	Ref.
Bullet Cluster		Scattering depth $\left({\Sigma }_{\mathrm{dm}}\approx 0.3{\mathrm{cm}}^2/g\right)$	Markevitch et al. (2004)
(1E 0657-558)		Mass loss	Randall et al. (2008)
		DM-galaxy offset $25\pm 29$ kpc	Randall et al. (2008)
Abell 520	$3.8\pm 1.1{\mathrm{cm}}^2/g$	Scattering depth $\left({\Sigma }_{\mathrm{dm}}\approx 0.07{\mathrm{cm}}^2/g\right)$	Mahdavi et al. (2007)
	$0.94\pm 0.06{\mathrm{cm}}^2/g$	Scattering depth $\left({\Sigma }_{\mathrm{dm}}\approx 0.14{\mathrm{cm}}^2/g\right)$	Jee et al. (2014a)
Abell 2744		Offset	Jauzac et al. (2016)
		Scattering depth $\left({\Sigma }_{\mathrm{dm}}\approx 0.3{\mathrm{cm}}^2/g\right)$	Merten et al. (2011)
Musket Ball Cluster		Scattering depth $\left({\Sigma }_{\mathrm{dm}}\approx 0.15{\mathrm{cm}}^2/g\right)$	Dawson et al. (2012)
(DLSCL J0916.2+2951)
Baby Bullet		Scattering depth $\left({\Sigma }_{\mathrm{dm}}\approx 0.25{\mathrm{cm}}^2/g\right)$	Bradac et al. (2008)
(MACS J0025.4-1222)
Abell 3827	$\sim 1.5{\mathrm{cm}}^2/g$	Offset	Kahlhoefer et al. (2015)

Equations

$$R_{\rm scat}=\sigma v_{\rm rel}\rho_{\rm dm}/m\approx 0.1\;{\rm Gyr}^{-1}\times\Big{(}\frac{\rho_{\rm dm}}{0.1\;{\rm\,M}_{\odot}/{\rm pc}^{3}}\Big{)}\Big{(}\frac{v_{\rm rel}}{50\;{\rm km/s}}\Big{)}\Big{(}\frac{\sigma/m}{1\;{\rm cm^{2}/g}}\Big{)}\,,$$

$$\sigma /m\sim 1\mathrm{c}{\mathrm{m}}^2/g\approx 2\times 10^{-24}\mathrm{c}{\mathrm{m}}^2/GeV$$

$$\sigma/m=\frac{4\pi\alpha^{\prime 2}m}{m_{\phi}^{4}}\approx 1\;{\rm cm^{2}/g}\times\Big{(}\frac{\alpha^{\prime}}{0.01}\Big{)}^{2}\Big{(}\frac{m}{10\;{\rm GeV}}\Big{)}\Big{(}\frac{m_{\phi}}{40\;{\rm MeV}}\Big{)}^{-4}$$

$${\rho }_{\mathrm{NFW}}(r)=\frac{{\rho }_s}{(r/r_s)(1+r/r_s{)}^2},$$

$$V_{\mathrm{circ}}(r)=\sqrt{V_{\mathrm{halo}}(r{)}^2+{{\rm Y}}_{\ast }{V}_{\mathrm{star}}(r{)}^2+V_{\mathrm{gas}}(r{)}^2},$$

$${\beta }_{\mathrm{aniso}}=1-{\sigma }_t^2/{\sigma }_r^2,$$

$$P_i={\rho }_{\mathrm{dm}}^{(i)}{v}_{\mathrm{rel}}^{(i)}(\sigma /m)\Delta t,$$

$$P_{ij}={\rho }_{ij}|{\mathbf{v}}_i-{\mathbf{v}}_j|(\sigma /m)\Delta t,$$

$$w(x)=\frac{8}{\pi}\left\{\begin{array}[]{cc}1-6x^{2}+6x^{3},&0\leq x\leq\tfrac{1}{2}\\ 2(1-x)^{3},&\tfrac{1}{2}<x\leq 1\\ 0,&x>1\end{array}\right..$$

$${\rho }_{ij}=m_p\int d^{3}rW(|\mathbf{r}-{\mathbf{r}}_i|,h)W(|\mathbf{r}-{\mathbf{r}}_j|,h).$$

$$\frac{\partial f}{\partial t}+\mathbf{v}\cdot \nabla f-\nabla {\Phi }_{\mathrm{tot}}\cdot \frac{\partial f}{\partial \mathbf{v}}=\mathcal{C}[f]$$

$$\nabla ({\sigma }_0^2{\rho }_{\mathrm{dm}})=-{\rho }_{\mathrm{dm}}\nabla {\Phi }_{\mathrm{tot}}.$$

$${\sigma }_0^2{\nabla }^2\ln {\rho }_{\mathrm{dm}}=-4\pi G({\rho }_{\mathrm{dm}}+{\rho }_b),$$

$$\rho_{\rm dm}(r)=\left\{\begin{array}[]{cc}\rho_{\rm iso}(r)\,,&r<r_{1}\\ \rho_{\rm NFW}(r)\,,&r>r_{1}\end{array}\right.\,.$$

$$⟨\sigma v_{\mathrm{rel}}⟩/m=({\rho }_{\mathrm{dm}}(r_1)t_{\mathrm{age}}{)}^{-1}$$

$${\rho }_b(R,z)={\Sigma }_b{{\rm Y}}_{\ast }{e}^{-R/R_d}\delta (z),$$

$$\tau ={\Sigma }_{\mathrm{dm}}\sigma /m<1.$$

$${\beta }_{\parallel }=1-e^{-\frac{\sigma /m}{{\sigma }_{\ast }}}$$

$${\sigma }_T=2\pi {\int }_0^{\pi }\frac{d\sigma }{d\Omega }(1-\cos \theta )\sin \theta d\theta ,$$

$${\sigma }_V=2\pi {\int }_0^{\pi }\frac{d\sigma }{d\Omega }(1-{\cos }^2\theta )\sin \theta d\theta .$$

$${\mathcal{L}}_{\mathrm{int}}=-\frac{\lambda }{4!}{\phi }^4$$

$$m_{\phi }\approx 8\mathrm{MeV}\times {\lambda }^{2/3}{\left(\frac{1\mathrm{c}{\mathrm{m}}^2/g}{\sigma /m_{\phi }}\right)}^{1/3},$$

$${\cal L_{\rm int}}=\bigg{\{}\begin{array}[]{c l}g_{\chi}\bar{\chi}\gamma^{\mu}\chi\phi_{\mu}&\text{(vector mediator)}\\ g_{\chi}\bar{\chi}\chi\phi&\text{(scalar mediator)}\,,\end{array}$$

$$V(r)=\pm \frac{{\alpha }_{\chi }}{r}{e}^{-m_{\varphi }r}.$$

$$\frac{d\sigma }{d\Omega }=\frac{{\alpha }_{\chi }^2{m}_{\chi }^2}{{\left[m_{\chi }^2{v}_{\mathrm{rel}}^2(1-\cos \theta )/2+m_{\varphi }^2\right]}^2}.$$

$$\sigma /m\sim 3\mathrm{c}{\mathrm{m}}^2/g\times \left(\frac{{\Lambda }_{\mathrm{DM}}}{m}\right){\left(\frac{{\Lambda }_{\mathrm{DM}}}{a^{-1}}\right)}^2{\left(\frac{100\mathrm{MeV}}{{\Lambda }_{\mathrm{DM}}}\right)}^3.$$

$$ma\gg 1,$$

$$V(x)\approx {\alpha }^2{\mu }_H\left[e^{0.09678-1.10173x-0.003945x^2}-e^{-(10.04/x-1{)}^2}\left(\frac{6.5}{x^6}+\frac{124}{x^8}+\frac{3285}{x^{10}}\right)\right],$$

$$\sigma /m_H\approx 1\mathrm{c}{\mathrm{m}}^2/g\times {\left(\frac{0.1}{\alpha }\right)}^2{\left(\frac{8\mathrm{GeV}}{m_p}\right)}^3{\left(\frac{m_p/m_e}{15}\right)}^2.$$

$$V(r)=\left(\begin{array}{cc}0& \frac{-{\alpha }_{\chi }}{r}{e}^{-m_{\varphi }r}\\ \frac{-{\alpha }_{\chi }}{r}{e}^{-m_{\varphi }r}& 2\delta \end{array}\right)$$

Full text

1. I Introduction

2. I.1 The dark matter puzzle

3. I.2 Crisis on small scales

4. I.3 Self-Interacting dark matter

5. I.4 From astrophysics to particle theory

6. II Astrophysical observations

7. II.1 Core-Cusp Problem

8. II.2 Diversity Problem

9. II.3 Missing Satellites Problem

10. II.4 Too-Big-to-Fail Problem

11. II.5 Baryon feedback

12. III N-body simulations and SIDM halo properties

13. III.1 Implementing self-interactions in simulations

14. III.2 Halo density profiles

15. III.3 Halo shapes: ellipticity

16. III.4 Substructure

17. III.5 SIDM simulations with baryons

18. IV Jeans approach to relaxed SIDM halos

19. IV.1 Isothermal solutions to the Jeans equations

20. IV.2 Diverse rotation curves for SIDM

21. IV.3 Scaling relations

22. V Halo mergers

23. V.1 Observations

24. V.2 Self-interactions in merging clusters

25. V.3 Self-interactions in minor mergers

26. VI Particle physics models

27. VI.1 What cross section is relevant?

28. VI.2 Self-coupled scalar

29. VI.3 Light mediator models

30. VI.4 Strongly interacting dark matter

31. VI.5 Dark atoms

32. VI.6 SIDM with an excited state

33. VII Complementary searches

34. VII.1 Portals for light mediators

35. VII.2 Direct detection

36. VII.3 Indirect detection

37. VII.4 Collider searches

38. VII.5 Cosmological probes

39. VIII Conclusions

40. VIII.1 Summary

41. VIII.2 Outlook

I Introduction

I.1 The dark matter puzzle

It was long ago pointed out by Oort that the distributions of mass and light in galaxies have little resemblance to one another Oort (1940). At the time, observations of galaxy NGC 3115 found a rotation curve rising linearly with radius—indicating a constant mass density—despite its luminosity falling by over an order of magnitude. Oort derived a mass-to-light ratio that increased with distance up to $\sim 250$ in solar units111Modern distance estimates revise Oort’s value down to $\sim 40$. at the outermost radius. He concluded that the luminous disk of NGC 3115 is “imbedded in a more or less homogeneous mass of great density” Oort (1940). The same phenomenon was observed in M31 by Babcock a year earlier Babcock (1939). These observations stood in contrast to the Milky Way (MW), where, by studying vertical dynamics of nearby stars, the mass-to-light ratio within the disk was known to be only $\sim 2$ at the solar radius Öpik (1915); Kapteyn (1922); Jeans (1922); Oort (1932). Earlier pioneering observations of the Coma cluster by Zwicky Zwicky (1933) and Virgo cluster by Smith Smith (1936) inferred similarly striking discrepencies between mass and light on much larger scales.

These issues were not fully appreciated for several decades, when astronomers were able to study galactic rotational velocities at much larger distances. While many rotation curves exhibited a linearly increasing velocity at small distances, they were expected to turn over eventually and fall as $r^{-1/2}$ according to Kepler’s laws King (1963). Instead, optical observations by Rubin and Ford of M31 Rubin and Ford (1970)—soon after extended to larger radii using 21-cm radio observations Roberts and Whitehurst (1975); Roberts and Rots (1973)—revealed a circular velocity that did not fall off, but remained approximately constant. Many other spiral galaxies were found to exhibit the same behavior Freeman (1970); Roberts and Rots (1973); Rubin et al. (1980); Bosma (1981); Sofue and Rubin (2001), indicating that most of the mass in galaxies is found in massive, nonluminous halos extending far beyond the spatial extent of luminous stars and gas Ostriker et al. (1974); Einasto et al. (1974). Massive spherical dark halos could also explain the apparent stability of bulgeless spiral galaxies Ostriker and Peebles (1973), which by themselves are unstable to the formation of bars Toomre (1964); Hohl (1971).

The large amount of dark matter (DM) required by observations pointed toward its nonbaryonic nature.222In this context, the term “baryons” represents protons, neutrons, and electrons that constitute normal atomic matter. The total mass density found in halos was estimated to be around $20\%$ of the critical density Ostriker et al. (1974); Einasto et al. (1974), in remarkable agreement with present cosmological values Ade et al. (2016). If composed of baryons, the cosmological baryon density would be in tension with upper limits inferred from nucleosynthesis arguments Wagoner (1973), as well as being difficult to hide from astrophysical observations Hegyi and Olive (1983). Therefore, the “missing mass” puzzle in galaxies, as well as clusters, suggested the existence of a new dominant form of matter, such as an elementary particle leftover from the Big Bang.

At first, the neutrino appeared to be a promising DM candidate within the Standard Model (SM) of elementary particles Gershtein and Zeldovich (1966); Cowsik and McClelland (1972); Szalay and Marx (1976). A relic thermal bath of neutrinos, produced alongside the cosmic microwave background (CMB), could yield the required mass density if the neutrino mass was $\mathcal{O}(10\;{\rm eV})$. However, because neutrinos are “hot” DM—they decouple from photons and electrons around the nucleosynthesis epoch while still relativistic—free-streaming erases density fluctuations below supercluster scales (see, e.g., Ref. Primack and Gross (2000)). Numerical simulations have shown that top-down structure formation, where superclusters form first in the neutrino-dominated Universe and subsequently fragment to produce galaxies, is incompatible with galaxy clustering constraints White et al. (1983, 1984); Hut and White (1984).

Cosmological data has converged upon the $\Lambda$CDM paradigm as the standard model of cosmology (e.g., Ref. Bahcall et al. (1999)). Of the total mass-energy content of the Universe, approximately $26\%$ is cold dark matter (CDM) and $5\%$ is baryonic matter (while the remainder is consistent with a cosmological constant $\Lambda$), with a nearly scale-invariant spectrum of primordial fluctuations Ade et al. (2016). In this picture, structure in the Universe forms as primordial overdensities collapse under gravity. Since CDM, acting as a pressureless fluid, is more dominant and collapses more readily than baryonic matter, it provides the gravitational potential underlying the distribution of visible matter in the Universe. The observed matter power spectrum, as obtained from a variety of cosmological probes, is in remarkable agreement with $\Lambda$CDM cosmology, shown in Fig. 1. In addition, the $\Lambda$CDM model also explains many important aspects of galaxy formation Springel et al. (2006); Trujillo-Gomez et al. (2011).

Despite this success, all evidence to date for DM comes from its gravitational influence in the Universe. With no viable DM candidate within SM, the underlying theory for DM remains unknown. Many new particle physics theories proposed to address shortcomings of the SM simultaneously predict new particles that can be a DM candidate. Examples include weakly-interacting massive particles (WIMPs) motivated by the hierarchy problem, such as neutralinos in the supersymmetric models Goldberg (1983); Jungman et al. (1996) and Kaluza-Klein states in extra dimensional models Servant and Tait (2003); Cheng et al. (2002), as well as extremely light axion particles Preskill et al. (1983) associated with the solution to the strong CP problem in QCD Peccei and Quinn (1977). The comic abundance of these new particles can be naturally close to the DM abundance inferred from the cosmological observations (e.g., Ref. Lee and Weinberg (1977)). This coincidence has motivated decades of efforts to discover the particle physics realization of DM through experimental searches for new physics beyond the SM (e.g., see Bertone et al. (2005); Feng (2010); Arrenberg et al. (2013) and references therein).

On large scales, the structure of the Universe is consistent with DM particles that are cold, collisionless, and interact with each other and SM particles purely via gravity. WIMPs and axions have interactions with the SM that are potentially large enough to be detectable in the laboratory, while on astrophysical scales these interactions are negligible and these candidates behave as CDM. On the other hand, other particle physics candidates for DM may have interactions with the SM are too feeble to observe directly. Here, observational tests of structure and departures from the CDM paradigm play a complementary role in probing the particle physics of DM independently of its interactions with SM particles (as had been done for hot DM candidates).

I.2 Crisis on small scales

$\Lambda$CDM is an extremely successful model for the large scale structure of the Universe, corresponding to distances greater than $\mathcal{O}({\rm Mpc})$ today (see Fig. 1). On smaller scales, structure formation becomes strongly nonlinear and N-body simulations have become the standard tool to explore this regime. Cosmological DM-only simulations have provided several predictions for the structure and abundance of CDM halos and their substructure. However, it remains unclear whether these predictions are borne out in nature.

Since the 1990s, four main discrepancies between CDM predictions and observations have come to light.

• •

Core-cusp problem: High-resolution simulations show that the mass density profile for CDM halos increases toward the center, scaling approximately as $\rho_{\rm dm}\propto r^{-1}$ in the central region Dubinski and Carlberg (1991); Navarro et al. (1996a); Navarro et al. (1997). However, many observed rotation curves of disk galaxies prefer a constant “cored” density profile $\rho_{\rm dm}\propto r^{0}$ Flores and Primack (1994); Moore (1994); Moore et al. (1999a), indicated by linearly rising circular velocity in the inner regions. The issue is most prevalent for dwarf and low surface brightness (LSB) galaxies Burkert (1995); McGaugh and de Blok (1998); Côté et al. (2000); van den Bosch and Swaters (2001); Borriello and Salucci (2001); de Blok et al. (2001a, b); Marchesini et al. (2002); Gentile et al. (2005); Gentile et al. (2007a); Kuzio de Naray et al. (2006, 2008); Salucci et al. (2007), which, being highly DM-dominated, are appealing environments to test CDM predictions.

• •

Diversity problem: Cosmological structure formation is predicted to be a self-similar process with a remarkably little scatter in density profiles for halos of a given mass Navarro et al. (1997); Bullock et al. (2001). However, disk galaxies with the same maximal circular velocity exhibit a much larger scatter in their interiors Oman et al. (2015) and inferred core densities vary by a factor of ${\cal O}(10)$ Kuzio de Naray et al. (2010).

• •

Missing satellites problem: CDM halos are rich with substructure, since they grow via hierarchical mergers of smaller halos that survive the merger process Kauffmann et al. (1993). Observationally, however, the number of small galaxies in the Local Group are far fewer than the number of predicted subhalos. In the MW, simulations predict $\mathcal{O}(100-1000)$ subhalos large enough to host galaxies, while only 10 dwarf spheroidal galaxies had been discovered when this issue was first raised Moore et al. (1999b); Klypin et al. (1999). Nearby galaxies in the field exhibit a similar underabundance of small galaxies compared to the velocity function inferred through simulations Zavala et al. (2009); Zwaan et al. (2010); Trujillo-Gomez et al. (2011).

• •

Too-big-to-fail problem (TBTF): In recent years, much attention has been paid to the most luminous satellites in the MW, which are expected to inhabit the most massive suhalos in CDM simulations. However, it has been shown that these subhalos are too dense in the central regions to be consistent with stellar dynamics of the brightest dwarf spheroidals Boylan-Kolchin et al. (2011, 2012). The origin of the name stems from the expectation that such massive subhalos are too big to fail in forming stars and should host observable galaxies. Studies of dwarf galaxies in Andromeda Tollerud et al. (2014) and the Local Group field Garrison-Kimmel et al. (2014) have found similar discrepancies.

It must be emphasized, however, that these issues originally gained prominence by comparing observations to theoretical predictions from DM-only simulations. Hence, there has been extensive debate in the literature whether these small scale issues can be alleviated or solved in the $\Lambda$CDM framework once dissipative baryonic processes, such as gas cooling, star formation, and supernova feedback, are included in simulations Navarro et al. (1996b); Governato et al. (2010). We review these topics in detail in §II.

A more intriguing possibility is that the CDM paradigm may break down on galactic scales. One early attempt to solve these issues supposes that DM particles are warm, instead of cold, meaning that they were quasi-relativistic during kinetic decoupling from the thermal bath in the early Universe Colombi et al. (1996); Bode et al. (2001). Compared to CDM, warm DM predicts a damped linear power spectrum due to free-streaming, resulting in a suppression of the number of substructures. Warm DM halos are also typically less concentrated because they form later than CDM ones. Recent high-resolution simulations show that warm DM may provide a solution to the missing satellites and too-big-to-fail problems Lovell et al. (2012, 2014); Horiuchi et al. (2016). However, the favored mass range of thermal warm DM is in strong tension with Lyman-$\alpha$ forest observations Iršič et al. (2017); Viel et al. (2013)333The Lyman-$\alpha$ constraints may be weakened due to uncertainties in the evolution of intergalactic medium Kulkarni et al. (2015); Garzilli et al. (2017); Cherry and Horiuchi (2017). and the abundance of high redshift galaxies Menci et al. (2016). Also, while warm DM halos have constant density cores set by the phase space density limit, the core sizes are far too small to solve the core-cusp problem given Lyman-$\alpha$ constraints Maccio et al. (2012).

I.3 Self-Interacting dark matter

Another promising alternative to collisionless CDM is self-interacting dark matter (SIDM), proposed by Spergel & Steinhardt to solve the core-cusp and missing satellites problems Spergel and Steinhardt (2000). In this scenario, DM particles scatter elastically with each other through $2\to 2$ interactions. Self-interactions lead to radical deviations from CDM predictions for the inner halo structure, shown in Fig. 2. We summarize the expectations for SIDM halos (blue) compared to CDM halos (black) as follows:

• •

Isothermal velocity dispersion: Although a CDM halo is a virilized object, the DM velocity dispersion, indicating the “temperature” of DM particles, is not a constant and decreases towards the center in the inner halo. Self-interactions transport heat from the hotter outer to the cooler inner region of a DM halo, thermalizing the inner halo and driving the velocity dispersion to be uniform with radius (Fig. 2, left panel). The velocity distribution function for SIDM becomes more Maxwell-Boltzmann compared to CDM Vogelsberger and Zavala (2013).

• •

Reduced central density: Hierarchical structure formation leads to a universal density profile for CDM halos Navarro et al. (1996a); Navarro et al. (1997). In the presence of collisions, the central density is reduced as low-entropy particles are heated within the dense inner halo, turning a cusp into a core (Fig. 2, center panel).

• •

Spherical halo shape: While CDM halos are triaxial Dubinski and Carlberg (1991), collisions isotropize DM particle velocities and tend to erase ellipticity. The minor-to-major axis ratio $c/a$ is closer to unity toward the center of SIDM halos compared to CDM halos (Fig. 2, right panel).

Since the scattering rate is proportional to the DM density, SIDM halos have the same structure as CDM halos at sufficiently large radii where the collision rate is negligible.

The local collision rate is given by

[TABLE]

where $m$ is the DM particle mass, while $\sigma,v_{\rm rel}$ are the cross section and relative velocity, respectively, for scattering. Within the central region of a typical dwarf galaxy, we have $\rho_{\rm dm}\sim 0.1\;{\rm\,M}_{\odot}/{\rm pc}^{3}$ and $v_{\rm rel}\sim 50\;{\rm km/s}$ Oh et al. (2011a). Therefore, the cross section per unit mass must be at least

[TABLE]

to have an effect on the halo, corresponding to at least one scattering per particle over 10 Gyr galactic timescales. For $\sigma/m\sim 1\;{\rm cm^{2}/g}$, the mean free path of DM particles is larger than the core radii (Knudsen number larger than unity) and heat conduction is effective in the inner halo. Provided $\sigma/m$ is not dramatically larger than this value, $R_{\rm scat}$ is negligible during the early Universe when structure forms.444Although self-interactions may be active in the very early Universe, long before matter-radiation equality, they rapidly fall out of equilibrium due to the Hubble dilution and redshifting of DM particles. Therefore SIDM retains the success of large-scale structure formation from $\Lambda$CDM, affecting structure at late times and only on small scales in the dense inner regions of halos.

Cored density profiles lead to shallower rotation curves for dwarf and LSB spiral galaxies at small radii, in accord with observations, while cores in satellite galaxies ameliorate the too-big-to-fail problem by reducing the predicted stellar line-of-sight velocity dispersion for the largest subhalos. If $\sigma/m$ is fixed as in Eq. (2), the effect of self-interactions on the halo scales approximately in a self-similar fashion, and larger halos, such as those for massive elliptical galaxies and clusters, may be impacted by self-interactions at proportionally larger radii.

In principle, self-interactions can also affect substructure, reducing the subhalo mass function to solve the missing satellites problem. SIDM subhalos are prone to tidal disruption, being less concentrated than their CDM counterparts, as well as evaporation due to ram pressure stripping from the host halo Spergel and Steinhardt (2000). However, this mechanism requires a too-large cross section that is excluded for reasons we discuss in §III.

The effect of self-interactions on DM halos has been borne out through N-body simulations. Soon after Spergel & Steinhardt’s proposal, a first wave of SIDM simulations was performed Moore et al. (2000); Yoshida et al. (2000a); Burkert (2000); Kochanek and White (2000); Yoshida et al. (2000b); Dave et al. (2001); Colin et al. (2002). In conjunction with these simulations, a number of constraints on SIDM emerged, the most stringent of which were strong lensing measurements of the ellipticity and central density of cluster MS2137-23 Miralda-Escude (2002); Meneghetti et al. (2001). These studies limited the cross section to be below $0.1\;{\rm cm^{2}/g}$, which excludes self-interactions at a level to explain small scale issues in galaxies. The Bullet Cluster provides additional evidence for the collisionless nature of DM, requiring $\sigma/m\lesssim 1\;{\rm cm^{2}/g}$ Randall et al. (2008). Although ad hoc velocity dependencies were put forth to evade these constraints Yoshida et al. (2000b); Firmani et al. (2000, 2001); Colin et al. (2002), SIDM largely fell into disfavor in light of these difficulties.

More recently, new N-body simulations for SIDM, with dramatically higher resolution and halo statistics, have revived DM self-interactions Vogelsberger et al. (2012); Rocha et al. (2013); Peter et al. (2013); Zavala et al. (2013); Elbert et al. (2015); Vogelsberger et al. (2014); Fry et al. (2015); Dooley et al. (2016). Many of the constraints from larger scales are much weaker than previously thought Rocha et al. (2013); Peter et al. (2013). In particular, constraints based on ellipticity have been overestimated: self-interactions do not erase triaxiality as effectively as previously supposed, and moreover, an observed ellipticity has contributions along the line-of-sight from regions outside the core where the halo remains triaxial Peter et al. (2013). The conclusion is that $\sigma/m\sim 0.5-1\;{\rm cm^{2}/g}$ can solve the core-cusp and TBTF issues on small scales, while remaining approximately consistent with other astrophysical constraints on larger scales Rocha et al. (2013); Peter et al. (2013); Zavala et al. (2013); Elbert et al. (2015). However, more recent studies based on stacked merging clusters Harvey et al. (2015) and stellar kinematics within cluster cores Kaplinghat et al. (2016) suggest some tension with these values. Viable SIDM models are preferred to have a scattering cross section with a mild velocity-dependence from dwarf to cluster scales. We discuss these issues in further detail in §III.

On the theory side, a new semi-analytical SIDM halo model has been developed based on the Jeans equation Kaplinghat et al. (2014a); Kaplinghat et al. (2016); Kamada et al. (2016). It can reproduce the simulation results for SIDM profiles within 10–20% while being much cheaper computationally. Discussed in §IV, this approach provides insight for understanding the baryonic influence on SIDM halo properties Kaplinghat et al. (2014a), testing SIDM models from dwarf to cluster scales Kaplinghat et al. (2016), and addressing the diversity in rotation curves Kamada et al. (2016).

Furthermore, there has been important progress in particle physics models for SIDM. Both numerical and analytical methods have been developed to accurately calculate the cross section for SIDM models involving the Yukawa Feng et al. (2010a); Tulin et al. (2013a); Schutz and Slatyer (2015) or atomic interactions Cline et al. (2014a); Boddy et al. (2016). These studies make it possible to map astrophysical constraints on $\sigma/m$ to the particle model parameters, such as the DM and mediator masses and coupling constant. In addition, an effective theory approach has been proposed in parametrizing SIDM models with a set of variables that are directly correlated with astrophysical observations Cyr-Racine et al. (2016); Vogelsberger et al. (2016).

In Table 1, we summarize the present status of astrophysical observations related to SIDM . The positive observations indicate discrepancies with CDM-only simulations and the required cross section assuming self-interactions are responsible for solving each issue. Dwarf and LSB galaxies favor cross sections of at least $1\;{\rm cm^{2}/g}$ to produce large enough core radii in these systems, which is also consistent with alleviating the too-big-to-fail problem among MW satellites and field dwarfs of the Local Group. The cross section need not be particularly fine-tuned. Values as large as $50\;{\rm cm^{2}/g}$ provide consistent density profiles on dwarf scales Elbert et al. (2015). (The upper limit on $\sigma/m$ on dwarf scales due to the onset of gravothermal collapse, which we discuss in §III, is not well-known.)

Lastly, for completeness, let us mention that SIDM, in its original conception, was introduced much earlier by Carlson, Machacek & Hall Carlson et al. (1992) to modify large scale structure formation. In their scenario, DM particles (e.g., glueballs of a nonabelian dark sector) have $3\to 2$ (or $4\to 2$) number-changing interactions by which they annihilate and eventually freeze-out to set the relic density, as well as $2\to 2$ elastic self-interactions that maintain kinetic equilibrium. Unlike standard freeze-out Scherrer and Turner (1986), DM particles retain entropy as their number is depleted and therefore cool more slowly than the visible sector. When this theory was proposed, a mixture of both cold plus hot DM within a flat, matter-dominated Universe ($\Lambda=0$) was a viable and theoretically appealing alternative to $\Lambda$CDM Davis et al. (1992); Klypin et al. (1993); Primack (1997). Although this original SIDM was proposed as a hybrid between hot plus cold DM, it provided too much small scale power suppression in Lyman-$\alpha$ observations relative to larger scales and was found to be excluded de Laix et al. (1995).555We emphasize that this exclusion is for a matter-dominated cosmology. Recently, this SIDM scenario has been revived within an observationally consistent (i.e., $\Lambda$-dominated) cosmology Hochberg et al. (2014), discussed in §VI.

I.4 From astrophysics to particle theory

The figure of merit for self-interactions, $\sigma/m$, depends on the underlying DM particle physics model. WIMPs have self-interactions mediated through the weak force, Higgs boson, or other heavy states. Since WIMP interactions and masses are set by the weak scale, yielding $\sigma/m\sim 10^{-38}\;{\rm cm^{2}/GeV}$, they effectively behave as collisionless CDM. If self-interactions indeed explain the small scale issues, then DM cannot be a usual WIMP.

What sort of particle theory can have a large enough $\sigma/m$? An analogy is provided by nuclear interactions, mediated by pion exchange. In Fig. 3, we show $n$-$p$ elastic scattering data for $n$ kinetic energies 0.5 eV $-$ 10 MeV Obložinský et al. (2011), expressed in terms of $v_{\rm rel}$ and $\sigma/m$ (here, $m$ is the nucleon mass). The required cross section for SIDM is comparable in magnitude to nuclear cross sections for visible matter. The lesson here is that $1\;{\rm cm^{2}/g}$ or larger can be achieved if the interaction scale lies below $\sim 1$ GeV. However, unlike nuclear scattering, the theory of self-interactions need not be strongly-coupled, nor does the DM mass need to be below 1 GeV. For example, self-interactions can be a weakly-coupled dark force Ackerman et al. (2009); Feng et al. (2009, 2010a); Buckley and Fox (2010); Loeb and Weiner (2011); Tulin et al. (2013b, a), with the mediator particle denoted by $\phi$. A perturbative calculation gives (in the limit $v_{\rm rel}=0$)

[TABLE]

where $\alpha^{\prime}$ is the DM analog of the electromagnetic fine structure constant, $\alpha_{\rm em}\approx 1/137$. Self-interactions that are electomagnetic strength (or weaker) are sufficient, as are weak-scale DM masses, provided the mediator mass $m_{\phi}$ is light enough. Note that Eq. (3) is only valid in the weakly-coupled Born limit, $\alpha^{\prime}m/m_{\phi}\ll 1$, and there are important corrections outside of this regime Tulin et al. (2013a).

Another important point is that $\sigma$ need not be constant in velocity. DM particles in larger mass halos have larger characteristic velocities compared to those in smaller halos. Therefore, observations from dwarf to cluster scales provide complementary handles for constraining the velocity-dependence of $\sigma$. Along these lines, Refs. Yoshida et al. (2000b); Firmani et al. (2000, 2001); Colin et al. (2002) advocated a cross section of the form $\sigma\propto 1/v_{\rm rel}$ in order to evade cluster constraints and fit a constant central density across all halo scales. We caution that such a dependence is not motivated by any theoretic framework for SIDM.

On the other hand, many well-motivated particle physics scenarios predict velocity-dependent cross sections that are not described by simple power laws. Scattering through a dark force can be a contact-type interaction at small velocity, with constant $\sigma$; yet once the momentum transfer is larger than the mediator mass, scattering is described by a Rutherford-like cross section that falls with higher velocity. These models naturally predict a larger $\sigma/m$ on dwarf scales compared to clusters, depending on the model parameters. The velocity-dependence may be qualitatively similar to $n$-$p$ scattering, shown in Fig. 3, or it can have a more complicated form Tulin et al. (2013b, a). We discuss models for SIDM and their expected velocity-dependence in §VI.

The particle physics of SIDM is not just about self-interactions. While searching for their effect on structure is the only probe if DM is completely decoupled from the SM, this is unlikely from an effective field theory point of view. For example, a dark force can interact with the SM through hypercharge kinetic mixing Holdom (1986) if it is a vector or the Higgs portal operator Patt and Wilczek (2006) if it is a scalar. Such couplings allow mediator particles thermally produced in the early Universe to decay, avoiding overclosure Kaplinghat et al. (2014b). At the same time, these operators provide a window to probing SIDM in direct and indirect detection experiments, as well as collider searches. Similar to the WIMP paradigm, SIDM motivates a multifaceted program of study combining complementary data from both astrophysical and terrestrial measurements. Though the coupling between the dark and visible sectors must be much weaker than for WIMPs, this can be compensated by the fact that mediator states must be much lighter than the weak scale. We discuss in §VII these complementary searches for DM (and other dark sector states) within the context of the SIDM paradigm.

II Astrophysical observations

The kinematics of visible matter is a tracer of the gravitational potential in galaxies and clusters, allowing the underlying DM halo mass distribution to be inferred. Observations along these lines have pointed to the breakdown of the collisionless CDM paradigm on small scales. In this section, we review these astrophysical observations and discuss possible solutions from baryonic physics and other systematic effects. For other recent reviews of these issues, we direct the reader to Refs. de Blok (2010); Bullock (2010); Weinberg et al. (2014); Del Popolo and Le Delliou (2017).

The issues discussed here may share complementary solutions. For example, mechanisms that generate cored density profiles may help reconcile the too-big-to-fail problem by reducing central densities of MW subhalos, as well as accommodating the diversity of galactic rotation curves. On the other hand, the too-big-to-fail problem may share a common resolution with the missing satellites problem if the overall subhalo mass function is reduced.

II.1 Core-Cusp Problem

Collisionless CDM-only simulations predict “cuspy” DM density profiles for which the logarithmic slope, defined by $\alpha=d\ln\rho_{\rm dm}/d\ln r$, tends to $\alpha\sim-1$ at small radii Dubinski and Carlberg (1991); Navarro et al. (1996a); Navarro et al. (1997); Moore et al. (1999a). Such halos are well-described by the Navarro-Frenk-White (NFW) profile Navarro et al. (1996a); Navarro et al. (1997),

[TABLE]

where $r$ is the radial coordinate and $\rho_{s}$ and $r_{s}$ are characteristic density and scale radius of the halo, respectively.666High resolution simulations have found that collisionless CDM profiles become shallower than $\alpha\sim-1$ at small radii Navarro et al. (2010); Stadel et al. (2009), in better agreement with the Einasto profile Einasto (1965). However, the enclosed mass profile is slightly larger at small radii compared to NFW fits Navarro et al. (2010), further exacerbating the issues discussed here. On the other hand, as we discuss below, many observations do not find evidence for the steep inner density slope predicted for collisionless CDM, preferring “cored” profiles with inner slopes $\alpha\sim 0$ that are systematically shallower Flores and Primack (1994); Moore (1994). This discrepancy is known as the “core-cusp problem.”

A related issue, known as the “mass deficit problem,” emerges when observed halos are viewed within a cosmological context McGaugh and de Blok (1998); de Blok et al. (2001c); McGaugh et al. (2007); Oman et al. (2015). The mass-concentration relation predicted from cosmological CDM simulations implies a tight correlation between $\rho_{s}$ and $r_{s}$ (see, e.g., Refs. Ludlow et al. (2014); Dutton and Macci (2014); Rodriguez-Puebla et al. (2016)). Since cosmological NFW halos are essentially a single-parameter profile (up to scatter), determination of the halo at large radius fixes the halo at small radius. However, many observed systems have less DM mass at small radii compared to these expectations. Alternatively, if NFW halos are fit to data at both large and small radii, the preferred profiles tend to be less concentrated than expected cosmologically.

Rotation curves: Late-type dwarf and LSB galaxies are ideal laboratories for halo structure. These systems are DM dominated down to small radii (or over all radii) and environmental disturbances are minimized. Flores & Primack Flores and Primack (1994) and Moore Moore (1994) first recognized the core-cusp issue based on HI rotation curves for several dwarfs, which, according to observations, are well described by cored profiles Carignan and Freeman (1988); Lake et al. (1990); Jobin and Carignan (1990).777In these early studies, the Hernquist profile Hernquist (1990) was used to model collisionless CDM halos. It has the same behavior as the NFW profile in the inner regions, $\rho_{\rm dm}\propto r^{-1}$, but falls as $r^{-4}$ at large radii. On the other hand, a variety of cored profiles have been adopted in the literature, including the nonsingular isothermal profile Eddington (1926), Burkert profile Burkert (1995), and pseudo-isothermal profile $\rho_{\rm dm}(r)=\rho_{0}(1+r^{2}/r_{c}^{2})^{-1}$, where $\rho_{0}$ and $r_{c}$ are the core density and radius, respectively van Albada et al. (1985).

Indeed, rotation curve measurements for dwarfs and LSBs have been a long-standing challenge to the $\Lambda$CDM paradigm Burkert (1995); McGaugh and de Blok (1998); Côté et al. (2000); van den Bosch and Swaters (2001); Borriello and Salucci (2001); de Blok et al. (2001a, b); Marchesini et al. (2002); Gentile et al. (2005); Gentile et al. (2007a); Kuzio de Naray et al. (2006, 2008); Salucci et al. (2007). For axisymmetric disk galaxies, circular velocity can be decomposed into three terms

[TABLE]

representing the contributions to the rotation curve from the DM halo, stars, and gas, respectively. The baryonic contributions to the rotation curve are modeled from the respective surface luminosities of stars and gas. However, the overall normalization between stellar mass and light remains notoriously uncertain: stellar mass is dominated by smaller and dimmer stars, while luminosity is dominated by more massive and brighter stars. Estimates for the stellar mass-to-light ratio—denoted by $\Upsilon_{*}$ in Eq. (5)—rely on stellar population synthesis models and assumptions for the initial mass function, with uncertainties at the factor-of-two level Conroy (2013). Modulo this uncertainty, the DM profile can be fit to observations. For a spherical halo, the DM contribution to the rotation curve is $V_{\rm halo}(r)=\sqrt{GM_{\rm halo}(r)/r}$, where $G$ is Newton’s constant and $M_{\rm halo}(r)$ is the DM mass enclosed within $r$.

Fig. 4 illustrates these issues for dwarf galaxy DDO 154. The left panel shows the measured HI rotation curve Oh et al. (2015) compared to fits with cuspy (NFW) and cored profiles, which are shown in the right panel. The NFW halo has been chosen to fit the asymptotic velocity at large radii and match the median cosmological relation between $\rho_{s}$ and $r_{s}$ Kamada et al. (2016). However, this profile overpredicts $V_{\rm circ}$ in the inner region. This discrepancy is a symptom of too much mass for $r\lesssim 2~{}{\rm kpc}$, while the data favors a shallower cored profile with less enclosed mass. An NFW profile with alternative parameters can provide an equally good fit as the cored profile, but the required concentration is significantly smaller than preferred cosmologically Oh et al. (2015). We note not all dwarf galaxies require a density core, as we we will discuss later.

Recent high-resolution surveys of nearby dwarf galaxies have given further weight to this discrepancy. The HI Near Galaxy Survey (THINGS) presented rotation curves for seven nearby dwarfs, finding a mean inner slope $\alpha=-0.29\pm 0.07$ Oh et al. (2011a), while a similar analysis by LITTLE THINGS for 26 dwarfs found $\alpha=-0.32\pm 0.24$ Oh et al. (2015). These results stand in contrast to $\alpha\sim-1$ predicted for CDM.

However, this discrepancy may simply highlight the inadequacy of DM-only simulations to infer the properties of real galaxies containing both DM and baryons. One proposal along these lines is that supernova-driven outflows can potentially impact the DM halo gravitationally, softening cusps Navarro et al. (1996b); Oh et al. (2011b), which we discuss in further detail in §II E. Alternatively, the inner mass density in dwarf galaxies may be systematically underestimated if gas pressure—due to turbulence in the interstellar medium—provides radial support to the disk Dalcanton and Stilp (2010); Pineda et al. (2017). In this case, the observed circular velocity will be smaller than needed to balance the gravitational acceleration, as per Eq. (5), and purported cores may simply be an observational artifact.

In light of these uncertainties, LSB galaxies have become an attractive testing ground for DM halo structure. A variety of observables—low metallicities and star formation rates, high gas fractions and mass-to-light ratios, young stellar populations—all point to these galaxies being highly DM-dominated and having had a quiescent evolution de Blok and McGaugh (1997). Moreover, LSBs typically have larger circular velocities and therefore deeper potential wells compared to dwarfs. Hence, the effects of baryon feedback and pressure support are expected to be less pronounced.

Rotation curve studies find that cored DM profiles are a better fit for LSBs compared to cuspy profiles McGaugh and de Blok (1998); de Blok et al. (2001a, b); Kuzio de Naray et al. (2006, 2008). In some cases, NFW profiles can give reasonable fits, but the required halo concentrations are systematically lower than the mean value predicted cosmologically. Although early HI and long-slit H$\alpha$ observations carried concerns that systematic effects—limited resolution (beam-smearing), slit misalignment, halo triaxiality and noncircular motions—may create cores artificially, these issues have largely been put to rest with the advent of high-resolution HI and optical velocity fields (see Ref. de Blok (2010) and references therein). Whether or not baryonic feedback can provide the solution remains actively debated Kuzio de Naray and Spekkens (2011); Oman et al. (2015); Katz et al. (2016); Pace (2016). Cored DM profiles have been further inferred for more luminous spiral galaxies as well Gentile et al. (2004); Salucci et al. (2007); Donato et al. (2009).

Although observational challenges remain in interpreting the very inner regions of galaxies, other studies have shown that small scale issues persist out to larger halo radii as well. McGaugh et al. McGaugh et al. (2007) examined rotation curves restricted to $r\geq 1$ kpc for 60 disk galaxies spanning a large range of mass and types. Under plausible assumptions for $\Upsilon_{*}$, the inferred halo densities are systematically lower than predicted for CDM. Gentile et al. Gentile et al. (2007b) analyzed rotation curves of 37 spiral galaxies. For each galaxy, an NFW halo is identified by matching the total enclosed DM mass with the observed one at the last measured rotation curve point. The resulting NFW halos are overdense in the central regions, but less dense in the outer regions, compared to the DM distribution inferred observationally. These studies support the picture that cosmological NFW halos tend to be too concentrated to fit observed rotation curves.

Milky Way satellites: Stellar kinematics for dwarf spheroidal (dSph) galaxies around the MW support the existence of cored DM density profiles. These DM-dominated galaxies—unlike disk galaxies discussed above—are gas-poor, dispersion-supported systems in which the stellar kinematics are dominated by random motions. Since available observations only encode half of the full 6D phase space information (i.e., two spatial dimensions in the plane of the sky and one velocity dimension along the line of sight), there is a well-known degeneracy between halo mass and stellar velocity anisotropy, given by the parameter

[TABLE]

where $\sigma^{2}_{r,t}$ denotes the radial ($r$) and tangential ($t$) stellar velocity dispersions, which follows from the equilibrium Jeans equation. Fortunately, however, for a wide range of halo models and anisotropies, the halo mass can be robustly estimated at the half-light radius, the radius enclosing half of the luminosity of the galaxy, with little scatter Walker et al. (2009); Wolf et al. (2010).

Further studies have exploited the fact that Sculptor and Fornax dSphs have two chemo-dynamically distinct stellar subpopulations (metal-rich and metal-poor) Tolstoy et al. (2004). Both populations trace the same underlying gravitational potential—but at different radii—effectively constraining the halo mass at two distinct radii and allowing the DM halo profile to be determined. The analyses by Battaglia et al Battaglia et al. (2008), based on the Jeans equation, and by Amorisco & Evans Amorisco and Evans (2012), based on modeling stellar distribution functions, find that a cored DM profile provides a better fit to the velocity dispersion profiles of Sculptor compared to NFW. A third method by Walker & Peñarrubia Walker and Penarrubia (2011) reached a similar conclusion for both Sculptor and Fornax. Their method used a simple mass estimator to measure the slope of the mass profile, excluding a NFW profile at high significance.888Strigari et al. Strigari et al. (2014) performed a different analysis for the Sculptor galaxy and argued that an NFW profile is consistent with observations (see also Ref. Breddels and Helmi (2013)). The discrepancy could be due to different assumptions about the anisotropy of the stellar velocity dispersions. With the recent identification of three distinct stellar subpopulations in Fornax, Amorisco et al. Amorisco et al. (2013) further showed that the galaxy resides in a DM halo with a $\sim 1~{}{\rm kpc}$ constant density core, while an NFW halo is incompatible with observations unless its concentration is unrealistically low.

Observations of kinematically cold stellar substructures in dSphs lend indirect support to the presence of DM cores. Kleyna et al. Kleyna et al. (2003) detected a stellar clump with low velocity dispersion in the Ursa Minor dSph and argued that its survival against tidal disruption is more likely in a cored—rather than cuspy—host halo. Walker et al. Walker et al. (2006) found a similar clump in the Sextans dSph, also implying a cored halo Lora et al. (2013). Lastly, the wide spatial distribution of five globular clusters in the Fornax dSph again favors a cored host halo Sanchez-Salcedo et al. (2006); Goerdt et al. (2006); Read et al. (2006); Cole et al. (2012), as dynamical friction within a cuspy halo would cause these clusters to sink to the center in less than a Hubble time.

Galaxy clusters: There is evidence that the core-cusp problem may be present for massive galaxy clusters as well. Sand et al. Sand et al. (2002, 2004) have advocated that the inner mass profile in relaxed clusters can be determined using stellar kinematics of the brightest cluster galaxy (BCG) that resides there. In combination with strong lensing data, they find a mean inner slope $\alpha\approx-0.5$ for the three best-resolved clusters in their sample, shallower than expected for CDM, albeit with considerable scatter between clusters Sand et al. (2004). Further studies by Newman et al. Newman et al. (2009, 2011, 2013a, 2013b) have supported these findings. A joint analysis of stellar BCG kinematics and lensing data for seven clusters with masses $\sim(0.4-2)\times 10^{15}\,{\rm\,M}_{\odot}$ has found that, while the total mass density profiles are consistent with NFW, the density profiles for DM only become softer than NFW for $r\lesssim$ 30 kpc Newman et al. (2013a, b). The halos are described equally well by a shallower inner slope $\alpha\approx-0.5$ or a cored profile with core radius $\sim 10-20$ kpc.

Baryon dynamics from active galactic nuclei (AGN) can be important, similar to the effect of supernovae in smaller halos. Hydrodynamical simulations of clusters by Martizzi et al. Martizzi et al. (2012, 2013) have found that repeated outflows from gas accretion and ejection by AGN can produce a cored cluster halo. On the other hand, weaker feedback prescriptions in simulations by Schaller et al. Schaller et al. (2015a)—argued to produce a more realistic stellar mass function for cluster galaxies—do not produce cored DM profiles and yield total mass profiles that are somewhat steeper than NFW, seemingly in tension with the Newman et al. results Newman et al. (2013a, b).

Other systematic effects may also be relevant for reconciling this descrepency. BCGs are dispersion-supported galaxies dominated by the stellar density in their interiors and there is a systematic degeneracy between $\Upsilon_{*}$ and $\beta_{\rm aniso}$ that is important when disentangling the DM and baryonic densities. Schaller et al. Schaller et al. (2015a) argue that radially-biased $\beta_{\rm aniso}\approx 0.2-0.3$, as suggested by their simulations, can push their results into closer agreement with Refs. Newman et al. (2013a, b). On the other hand, Newman et al. Newman et al. (2009, 2013b) find that radially biased anisotropies only serve to make DM density profiles shallower in their fits. In any case, the comparison is not strictly apples-to-apples since the Newman clusters Newman et al. (2013a, b) are up to a factor of five larger than the largest simulated ones by Schaller et al. Schaller et al. (2015a).

II.2 Diversity Problem

In $\Lambda$CDM, hierarchical structure formation produces self-similar halos well-described by NFW profiles. Since the halo parameters (e.g., $\rho_{s}$ and $r_{s}$) are strongly correlated, there is only one parameter specifying a halo. For example, once the maximum circular velocity $V_{\rm max}$ (or any other halo parameter) is fixed, the halo density profile is completely determined at all radii including the inner density cusp (up to the scatter). On the other hand, the inner rotation curves of observed galaxies exhibit considerable diversity. Galaxies of the same $V_{\rm max}$ can have significant variation in their central densities. Any mechanism to explain the core-cusp issue must also accommodate this apparent diversity.

To illustrate this issue, Kuzio de Naray et al. Kuzio de Naray et al. (2010) fitted seven LSB galaxies with four different cored halo models, including, e.g., a pseudoisothermal profile $\rho_{\rm dm}(r)=\rho_{0}(1+r^{2}/r_{c}^{2})^{-1}$. Fig. 5 (left) shows the central DM density $\rho_{0}$—inferred by the inner slope of $V_{\rm circ}(r)$—versus $V_{\rm max}$ for these galaxies. Within the sample, there is no clear correlation between the inner ($\rho_{0}$) and outer ($V_{\rm max}$) parts of the halo. Moreover, the spread in $\rho_{0}$ can be large for galaxies with similar $V_{\rm max}$, up to a factor of ${\cal O}(10)$ when $V_{\rm max}\sim 80~{}{\rm km/s}$. The result is independent of the choice of halo model and mass-to-light ratio.

Instead of fitting to a specific halo profile, Oman et al. Oman et al. (2015) parametrized the diversity of rotation curves more directly by comparing $V_{\rm circ}(2~{}{\rm kpc})$ versus $V_{\rm max}$, which represent the inner and outer halos, respectively. Fig. 5 (right) shows the scatter in these velocities for observed galaxies (blue points) compared to the correlation expected from CDM-only halos (solid line) and CDM halos with baryons (red band). For $V_{\rm max}$ in the range of $50\textup{--}300~{}{\rm km/s}$, the spread in $V_{\rm circ}(2~{}{\rm kpc})$ is a factor of $\sim 3$ for a given $V_{\rm max}$. For example, when $V_{\rm max}\sim 70~{}{\rm km/s}$, CDM (only) predicts $V_{\rm circ}(2~{}{\rm kpc})\sim 50~{}{\rm km/s}$ (solid line), but observed galaxies span from $V_{\rm circ}(2~{}{\rm kpc})\sim 20~{}{\rm km/s}$ to $70~{}{\rm km/s}$. Galaxies at the low end of this range suffer from the mass deficit problem discussed above. For these outliers, including the baryonic contribution will make the comparison worse. On the other hand, galaxies at the upper end of this range could be consistent with CDM predictions once the baryonic contribution is included Oman et al. (2015). However, the spread in the baryon distribution plays a less significant role in generating the scatter in $V_{\rm circ}(2~{}{\rm kpc})$ for CDM halos since the enclosed DM mass in the cusp tends to dominate over the baryon mass.

II.3 Missing Satellites Problem

The hierarchical nature of structure formation predicts that CDM halos should contain large numbers of subhalos Kauffmann et al. (1993). Numerical simulations predict that a MW-sized halo has a subhalo mass function that diverges at low masses as $dN/dM_{\rm halo}\propto M^{-1.9}_{\rm halo}$ Springel et al. (2008). Consequently, the MW should have several hundred subhalos with $V_{\rm max}\sim 10-30~{}{\rm km/s}$ within its virial radius, each in principle hosting a galaxy Klypin et al. (1999); Moore et al. (1999b). However, only 11 dwarf satellite galaxies were known in the MW when the problem was originally raised in 1999, shown in Fig. 6 (left). The mismatch also exists between the abundance of observed satellites in the Local Group and that predicted in simulations. This conflict is referred as the “missing satellites problem.” We note that a similar descrepancy does not appear for galactic-scale substructure in galaxy clusters (shown in Fig. 6 (left) for the Virgo cluster).

One possibility is that these subhalos exist but are invisible because of the low baryon content. For low-mass subhalos, baryonic processes may play an important role for suppressing star formation. For instance, the ultraviolet photoionizing background can inhibit gas collapse into DM halos by heating the gas and reducing the gas cooling rate, which could suppress galaxy formation in halos with circular velocities less than $30~{}{\rm km/s}$ Thoul and Weinberg (1996); Bullock et al. (2000). In addition, after the initial star formation episode, supernova-driven winds could push the remaining gas out of the shallow potential wells of these low mass halos Dekel and Silk (1986).

The discovery of many faint new satellites in the Sloan Digital Sky Survey has suggested that as many as a factor of $\sim 5-20$ more dwarf galaxies could be still undiscovered due to faintness, luminosity bias, and limited sky coverage Tollerud et al. (2008); Walsh et al. (2009); Bullock et al. (2010). More recently, seventeen new candidate satellites have been found in the Dark Energy Survey Bechtol et al. (2015); Drlica-Wagner et al. (2015). Given these considerations, the dearth of MW subhalos may not be as severe as thought originally.

A similar abundance problem has arisen for dwarf galaxies in the field of the Local Volume. The velocity function—the number of galaxies as a function of their HI line widths—provides a useful metric for comparing to CDM predictions since HI gas typically extends out to large distances to probe $V_{\rm max}$ for the halo Cole and Kaiser (1989); Shimasaku (1993). While in accord with observations for larger galaxies, the velocity function for CDM overpredicts the number smaller galaxies with $V_{\rm max}\lesssim 80\;\,{\rm km/s}$ Zavala et al. (2009); Zwaan et al. (2010); Trujillo-Gomez et al. (2011). For example, Klypin et al. Klypin et al. (2015) find $\sim 200$ nearby galaxies within 10 Mpc with $V_{\rm max}\sim 30-50\;\,{\rm km/s}$, while CDM predicts $\mathcal{O}(1000)$. Unlike the satellites, which are considerably smaller and fainter, these galaxies are relatively bright dwarf irregulars where observations are essentially complete within this volume.

One explanation for this missing dwarf problem is that HI line widths may be biased tracers for $V_{\rm max}$. HI measurements for many dwarf galaxies are limited to the rising part of the rotation curve and therefore do not sample the full gravitational potential of the DM halo. Whether this bias can reconcile the observed velocity function with CDM Brook and Shankar (2016); Macciò et al. (2016); Brooks et al. (2017), or if the discrepancy still persists Trujillo-Gomez et al. (2016), remains to settled.

II.4 Too-Big-to-Fail Problem

Satellites within the Local Group: Boylan-Kolchin et al. Boylan-Kolchin et al. (2011, 2012) showed that the population of the MW’s brightest dSph galaxies, which are DM dominated at all radii, exhibit another type of discrepancy with CDM predictions. Since these satellites have the largest stellar velocities and luminosities, they are expected to live in the most massive MW subhalos. However, the most massive subhalos predicted by CDM-only simulations have central densities too large to host the observed satellites. As shown in Fig. 6 (right), simulations predict ${\cal O}(10)$ subhalos with $V_{\rm max}>30~{}{\rm km/s}$, whereas the bright MW dSphs have stellar dispersions corresponding to CDM subhalos with $12\lesssim V_{\rm max}\lesssim 25~{}{\rm km/s}$. It is puzzling that the most massive subhalos should be missing luminous counterparts since their deep potential wells make it unlikely that photoionizing feedback can inhibit gas accretion and suppress galaxy formation. Hence, these substructures should be too big to fail to form stars.

Several proposed mechanisms may address the TBTF problem without invoking DM physics. First, the MW halo mass may be underestimated. The TBTF discrepancy is based on simulated MW-like halos with masses in the range $1-2\times 10^{12}\,{\rm\,M}_{\odot}$ Boylan-Kolchin et al. (2011, 2012). However, since the number of subhalos scales with host halo mass, the apparent lack of massive subhalos might be accommodated if the MW halo mass is around $5\times 10^{11}~{}{\rm\,M}_{\odot}$, although Boylan-Kolchin et al. argue against this possibility (see also Ref. Wang et al. (2015) for a summary of different estimates.) Second, the MW may have less massive subhalos due to scatter from the stochastic nature of structure formation Purcell and Zentner (2012). However, a recent analysis finds that there is only $\sim 1\%$ chance that MW-sized host halos have a subhalo population in statistical agreement with that of the MW Jiang and Bosch (2015). In addition, similar discrepancies also exist for the brightest dwarf galaxies in Andromeda Tollerud et al. (2014) and the Local Group field Garrison-Kimmel et al. (2014), which further disfavor these explanations. Baryonic physics—including environmental effects from the MW disk—may also play a role, discussed in the next section.

Field galaxies: A similar TBTF problem also arises for dwarf galaxies in the field Ferrero et al. (2012); Papastergis et al. (2015); Garrison-Kimmel et al. (2014); Brook and Di Cintio (2015); Klypin et al. (2015). According to abundance matching, galaxies are expected to populate DM halos in a one-to-one relationship that is monotonic with mass (i.e., larger galaxies in larger halos). However, the galaxy stellar mass function, inferred by the Sloan Digital Sky Survey Baldry et al. (2008); Li and White (2009), is shallower at low mass than the halo mass function in $\Lambda$CDM, which suggests that galaxy formation becomes inefficient for halos below $\sim 10^{10}\,{\rm\,M}_{\odot}$ Guo et al. (2010). Hence, it is expected that most faint dwarf galaxies populate a narrow range of DM halos with $V_{\rm max}\sim 30\;{\rm km/s}$, while halos with $V_{\rm max}\lesssim 30\;{\rm km/s}$ would have no galactic counterpart Ferrero et al. (2012).

Ferrero et al. Ferrero et al. (2012) find that rotation curves for faint dwarf galaxies do not support these conclusions. A large fraction of faint galaxies in their sample appear to inhabit CDM halos with masses below $\sim 10^{10}\,{\rm\,M}_{\odot}$, which would imply that other more massive halos—which should be too big to fail—lack galaxies. On the other hand, optical galaxy counts may simply have missed a large fraction of low surface brightness objects.

More recently, Papastergis et al. Papastergis et al. (2015) have confirmed the findings of Ref. Ferrero et al. (2012) by utilizing data from the Arecibo Legacy Fast ALFA 21-cm survey. Since the majority of faint galaxies are gas-rich late-type dwarfs, these observations provide a more complete census without a bias against low stellar luminosity objects. Fig. 7 shows their main results, illustrating the TBTF problem in the field.

• •

The blue curve is the prediction of abundance matching between observed HI rotational velocities ($V_{\rm rot,H{\small{\sc I}}}$) and predicted maximum rotational velocities for $\Lambda$CDM halos ($V_{\rm halo,max}$).

• •

The various symbols indicate the observed correlation for a sample of 194 nearby galaxies. For each galaxy, $V_{\rm rot,H{\small{\sc I}}}$ is the rotational velocity at the outermost measured radius and $V_{\rm halo,max}$ is determined by matching a CDM halo rotation curve at that point.

• •

Since baryonic contributions to the rotation curve are neglected, each symbol is regarded as an upper limit on $V_{\rm halo,max}$. Hence, the symbols are indicated by leftward arrows.

Galaxies to the right of the blue line can potentially be reconciled with abundance matching predictions since baryons can make up the difference. On the other hand, the lowest mass galaxies to the left of the blue line seem inconsistent with these arguments. Abundance matching predicts they should inhabit halos with $V_{\rm halo,max}\sim 40\;{\rm km/s}$, yet this is incompatible with their measured HI circular velocities.

Common ground with the core-cusp problem: The essence of the TBTF problem is that low-mass galaxies have gas or stellar velocities that are too small to be consistent with the CDM halos they are predicted to inhabit. At face value, the issue is reminiscent of the core-cusp/mass deficit problems for rotation curves and other observations. Thus, one way to resolve the TBTF problem is if these galaxies have reduced central densities compared to CDM halos. By generating cored profiles in low-mass halos, self-interactions may resolve this issue for the dwarf galaxies in the MW Vogelsberger et al. (2012); Zavala et al. (2013), Local Group Elbert et al. (2015), and the field Schneider et al. (2016).

II.5 Baryon feedback

It has been more than 20 years, since the small scale “crisis” was first posed in the 1990s. Since then, there has been extensive discussion on whether all issues can be resolved within the CDM paradigm once baryonic feedback processes—gas cooling, star formation, supernovae, and active galactic nuclei—are accounted for. Here we give a brief review of feedback on galactic scales.

Cores in galaxies: As galaxies form, gas sinks into the inner halo to produce stars. The deepening gravitational potential of baryons causes the central density and velocity dispersion of DM to increase through adiabatic contraction Blumenthal et al. (1986). Although at first sight this makes the core-cusp problem worse, the situation can be very different due to non-adiabatic feedback from supernova-driven gas outflows Navarro et al. (1996b); Governato et al. (2010). If a sizable fraction of baryons is suddenly removed from the inner halo, DM particles migrate out to larger orbits. In this way, repeated bursts of star formation and outflows (followed by reaccretion) can reduce the central halo density through a purely gravitational interaction between DM and baryons. Non-adiabaticity is essential for the mechanism to work. This is achieved by assuming a “bursty” star formation history with $\mathcal{O}(10)$ variation in the star formation rate over time scales comparable to the dynamical time scale of the galaxy Teyssier et al. (2013). (See Refs. Pontzen and Governato (2014, 2012) for a pedagogical explanation.)

Governato et al. Governato et al. (2010) argued that bursty star formation in dwarf galaxies could be connected to another long-standing puzzle in galaxy formation: bulgeless disk galaxies. If strong outflows are necessary to remove low angular momentum gas to prevent bulge formation in certain galaxies, they may also induce DM cores. High-resolution hydrodynamical simulations show that supernova feedback can generate $\mathcal{O}({\rm kpc)}$ cores in CDM halos for dwarf galaxies Governato et al. (2010). Notably, the shallow inner slope of the DM distribution inferred from THINGS and LITTLE THINGS dwarf galaxies could be consistent with CDM halos once feedback is included (except DDO 210) Oh et al. (2011b); Governato et al. (2012); Oh et al. (2015) (see, e.g., Fig. 6 of Oh et al. (2015)).

Further studies by Di Cintio et al. Di Cintio et al. (2014a, b) investigated the effect of feedback as a function of halo mass. Their sample includes 31 galaxies from hydrodynamical simulations spanning halo mass $10^{10}\textup{--}10^{12}{\rm\,M}_{\odot}$. The corresponding stellar masses are fixed by abundance matching, with stellar-to-halo mass ratio in the range $M_{\rm star}/M_{\rm halo}\sim 10^{-4}\textup{--}10^{-1}$. Fig. 8 (left) shows how the DM inner density slope $\alpha$ (measured at $1\textup{--}2\%$ of the virial radius) depends strongly on $M_{\rm star}/M_{\rm halo}$. Di Cintio et al. Di Cintio et al. (2014a) conclude the following points:

• •

Maximal flattening of the inner DM halo occurs for $M_{\rm star}/M_{\rm halo}\sim 0.5\%$. According to abundance matching, this corresponds to stellar mass $\sim 3\times 10^{8}\,{\rm\,M}_{\odot}$, halo mass $\sim 6\times 10^{10}{\rm\,M}_{\odot}$, and asymptotic circular velocity $\sim 50~{}{\rm km/s}$. The THINGS dwarfs are right in the “sweet spot” for supernova feedback to be maximally effective.

• •

Supernova feedback is less effective in larger halos, which (by abundance matching) correspond to larger values of $M_{\rm star}/M_{\rm halo}$. Despite more star formation available to drive feedback, the deeper potential well for the halo suppresses the effect. For these galaxies, adiabatic contraction dominates and halos are cuspy.

• •

When the ratio is less than $10^{-4}$, there is too little star formation for feedback to affect the inner DM density. This implies that halos hosting galaxies with $M_{\rm star}\lesssim 3\times 10^{6}{\rm\,M}_{\odot}$ remain cuspy, consistent with earlier findings Governato et al. (2012)

These results have been further confirmed by the Numerical Investigation of Hundred Astrophysical Objects (NIHAO) simulations Tollet et al. (2016), where the same stellar feedback model was used for a larger halo sample. Other simulations with similar feedback prescriptions have also confirmed core formation in halos Teyssier et al. (2013); Madau et al. (2014).

Cosmological simulations with Feedback In Realistic Environments (FIRE) Hopkins et al. (2014) have emphasized that core formation in DM halos is tightly linked with star formation history O orbe et al. (2015); Chan et al. (2015). Oñorbe et al. O orbe et al. (2015) targeted low mass galaxies relevant for the TBTF problem and showed that kpc-sized cores may arise through feedback. However, cores form only when star formation remains active after most of the halo growth has occurred ($z\lesssim 2$). Otherwise, even if a core forms at early times, subsequent halo mergers may erase it if feedback has ceased. Lower mass galaxies, such as ultra-faint dwarfs, have too little star formation to affect the DM halo.

Chan et al. Chan et al. (2015) explored a wider mass range using the FIRE simulations, $M_{\rm star}\approx 10^{4}\textup{--}10^{11}{\rm\,M}_{\odot}$ and $M_{\rm halo}\approx 10^{9}\textup{--}10^{12}{\rm\,M}_{\odot}$. The DM density profiles become shallow for $M_{\rm halo}\approx 10^{10}\textup{--}10^{11}{\rm\,M}_{\odot}$ due to strong feedback, with the maximal effect at $M_{\rm halo}\sim 10^{11}{\rm\,M}_{\odot}$. The result is broadly consistent with Governato et al. Governato et al. (2012) and Di Cintio et al. Di Cintio et al. (2014a). Chan et al. also found that large cores can form only if bursty star formation occurs at a late epoch when cusp-building mergers have stopped, as pointed out in O orbe et al. (2015).

Although baryonic feedback may reconcile galactic rotation curves with the CDM paradigm (see, e.g., Refs. Brook (2015); Katz et al. (2016)), questions yet remain:

• •

Can feedback generate ultra-low density cores? Feedback-created cores are limited to the inner $\sim 1$ kpc where the star formation rate is high. However, some galaxies, such as dwarf IC 2574, are extreme outliers with huge core sizes and mass deficits that are far beyond what feedback can do. Fig. 8 (right) shows the rotation curve for IC 2574, together with simulated ones with similar values of $V_{\rm max}$. The observed core size is $\sim 6$ kpc Oh et al. (2011a)—set by the region over where $V_{\rm circ}(r)$ rises linearly with radius—which is too large to be created with bursty feedback (open symbols) Oh et al. (2011b) or smoother implementations (green band) Oman et al. (2015). Fig. 5 (right) shows that many galaxies have inner mass deficits out to 2 kpc and it is challenging for feedback to remove enough DM from their central regions to be consistent with observations.

• •

Does feedback explain the full spectrum of observed rotation curves? Galaxies exhibit a broad spread in rotation curves, encompassing both cored and cuspy profiles Oman et al. (2015). While some systems have large cores, there are galaxies with the same $V_{\rm max}$ which are consistent with a cuspy CDM halo. For fixed $V_{\rm max}$, the spread in the velocity at 2 kpc is a factor of $\sim 3$, as shown in Fig. 5 (right). It is unclear whether feedback prescriptions can account for such a scatter.

• •

Is bursty feedback required? Since star formation is far below the resolution of simulations, baryon dynamics depend on how feedback is modeled, especially the density threshold for star formation.999In various simulations, the gas density thresholds for star formation are taken to be $\gtrsim 10\textup{--}100$ (Brooks et al. Brooks et al. (2013)), $10\textup{--}100$ (FIRE, Chan et al. Chan et al. (2015)), $1000$ (FIRE-2, Wetzel et al. Wetzel et al. (2016)), and $0.1\textup{--}1~{}\rm{atoms/cm^{3}}$ (Oman et al. Oman et al. (2015)). Large thresholds lead to bursty star formation since energy injection to the dense medium causes strong outflows, which disturb the potential violently to generate DM cores Governato et al. (2010, 2012). However, Oman et al. Oman et al. (2015) have analyzed galaxies from the EAGLE and LOCAL GROUP simulation projects—which adopt a far smaller threshold—and find that feedback is negligible. These simulations have star formation occuring throughout the gaseous disk—without sudden fluctuations in the gravitational potential needed to produce DM cores—and yet produce galaxies that are consistent with other observational constraints.

• •

What is the epoch of core formation? Detailed comparisons between hydrodynamical simulations with bursty feedback would be useful to understand what systematic differences may be present (if any). One potential point of divergence is the epoch when star formation bursts need to occur to yield cores. Pontzen & Governato Pontzen and Governato (2012) argue that DM cores may form due to feedback at an early time $2<z<4$. However, Refs. O orbe et al. (2015); Chan et al. (2015) find that cores remain stable only if outflows remain active at later times, $z\lesssim 2$, once halo growth has slowed. This question is important for connecting observational tracers of starbursts in galaxies to feedback-driven core formation (e.g., Ref. Weisz et al. (2012); Kauffmann (2014)).

Substructure: A number of studies have found that baryon dynamics may solve the missing satellites and TBTF problem within $\Lambda$CDM, independently of whether star formation is bursty or not. Sawala et al. Sawala et al. (2016) performed a suite of cosmological hydrodynamical simulations of 12 volumes selected to match the Local Group. These simulations adopt a smooth star formation history, as in Ref. Oman et al. (2015), and do not yield cored profiles. Regardless, the number of satellite galaxies is reduced significantly in MW and M31-like halos (within 300 kpc), as well as in the broader Local Group (within 2 Mpc), in better accord with observations compared to the expectation from DM-only simulations. Supernova feedback and reionization deplete baryons in low mass halos and only a subset of them can form galaxies. In addition, the most massive subhalos in MW-like galaxies have an inner mass deficit due to ram pressure stripping and a suppressed halo growth rate due to the baryon loss. These massive subhalos could be consistent with kinematical observations of the MW dSphs, solving the TBTF problem in the MW, within measurement uncertainties Fattahi et al. (2016).

Earlier studies by Zolotov et al. Zolotov et al. (2012) and Brooks et al. Brooks and Zolotov (2014); Brooks et al. (2013) reached a similar conclusion in simulations with a bursty star formation scheme, as in Ref. Governato et al. (2012) (see Ref. Stinson et al. (2006) for details). Strong feedback produces cored profiles in subhalos, which in turn enhances tidal effects from the stellar disk of host galaxies (see also Penarrubia et al. (2010)). More recently, FIRE simulations in the Local Group environment by Wetzel et al. Wetzel et al. (2016) also showed that the population of satellite galaxies with $M_{\rm star}\gtrsim 10^{5}{\rm\,M}_{\odot}$ does not suffer from the missing satellites and TBTF problems.

Too-big-to-fail in the field: Papastergis & Shankar Papastergis and Shankar (2016) argued that the TBTF problem for field dwarf galaxies cannot be solved in $\Lambda$CDM even with baryonic feedback effects. They adapted the halo velocity function proposed in Sawala et al. Sawala et al. (2016), which includes the effects of baryon depletion and reionization feedback, and the modified CDM halo density profile due to supernova feedback from Di Cintio et al. Di Cintio et al. (2014a). While reionization effectively suppresses star formation for halos with $V_{\rm max}\lesssim 20~{}{\rm km/s}$, the TBTF problem here concerns halos with $V_{\rm max}\approx 25\textup{--}45~{}{\rm km/s}$, which are too massive to be affected significantly. In addition, cored profiles may help alleviate the tension for galaxies whose stellar kinematics are measured only in the very inner region, but not those with $V_{\rm circ}$ observed at large radii.

More recently, Verbeke et al. Verbeke et al. (2017) argued that the TBTF problem in the field may be resolved by noncircular motions in HI gas. According to their simulations, HI kinematics is an imperfect tracer for mass in small galaxies due to turbulence (as pointed out in Refs. Dalcanton and Stilp (2010); Pineda et al. (2017)). In particular, galaxies with $V_{\rm rot,H{\small{\sc I}}}\lesssim 30$ km/s (see Fig. 7) may be consistent with living in larger halos than would otherwise be inferred without properly accounting for this systematic effect.

Summary: The predictions of $\Lambda$CDM cosmology can be modified on small scales when dissipative baryon physics is included in simulations. Strong baryonic feedback from supernova explosions and stellar winds disturb the gravitational potential violently, resulting in a shallow halo density profile. Reionization can suppress galaxy formation in low mass halos and strong tides from host’s stellar disk can destroy satellite halos. These effects can also be combined in shaping galaxies. For example, if feedback induces cores in low mass halos, galaxy formation is more effectively suppressed by reionization and the halo is more vulnerable to tidal disruption. In simulations, the significance of these effects depends on the specific feedback models, in particular, the gas density threshold for star formation. It seems that strong feedback—leading to the core formation and significant suppression of star formation in low mass halos—is not required for reproducing general properties of observed galaxies in simulations (see, e.g., Refs. Oman et al. (2015); Sawala et al. (2016); Sales et al. (2017)). Thus, it remains unclear to what degree baryon dynamics affect halo properties in reality. Moreover, cored profiles generated by feedback, as proposed in Refs. Di Cintio et al. (2014a, b), may be inconsistent with the correlations predicted in $\Lambda$CDM cosmology, namely the mass-concentration and $M_{\rm star}\textup{--}M_{\rm halo}$ abundance matching relations Pace (2016) (however, see Ref. Katz et al. (2016)). Thus, it remains an intriguing possibility that small scale issues may imply the breakdown of the CDM paradigm on galactic scales, as we will focus on in the rest of this article.

III N-body simulations and SIDM halo properties

N-body simulations have been the primary tools for understanding the effect of self-interactions on structure Moore et al. (2000); Yoshida et al. (2000a); Burkert (2000); Kochanek and White (2000); Yoshida et al. (2000b); Dave et al. (2001); Colin et al. (2002); Vogelsberger et al. (2012); Rocha et al. (2013); Peter et al. (2013); Zavala et al. (2013); Elbert et al. (2015); Vogelsberger et al. (2014); Fry et al. (2015); Dooley et al. (2016). Early SIDM simulations used smoothed particle hydrodynamics to model DM collisions Moore et al. (2000); Yoshida et al. (2000a). This approach treats SIDM as an ideal gas described by fluid equations, which are valid in the optically-thick regime where the mean free path $\lambda_{\rm mfp}$ is much smaller than typical galactic length scales. Due to efficient thermalization, SIDM halos in this context form a singular isothermal profile, $\rho_{\rm dm}\propto r^{-2}$, which is steeper than collisionless CDM halos and exacerbates the core-cusp issue rather than solving it Moore et al. (2000); Yoshida et al. (2000a).101010The fluid equations approach has been used as a semi-analytical framework for modeling SIDM halos under simplifying assumptions (spherical symmetry and self-similar evolution) Balberg et al. (2002); Ahn and Shapiro (2005); Koda and Shapiro (2011). In fact, this model agrees well with N-body simulations for isolated halos, even for the optically-thin regime that is beyond the validity of the fluid approximation; however, the model is not yet able to reproduce N-body results for cosmological halos due to departure from self-similar evolution Koda and Shapiro (2011).

Consequently, most SIDM simulations have focused the more promising case of self-interactions in the optically-thin regime, with cross sections spanning $0.1-50\;{\rm cm^{2}/g}$  Burkert (2000); Kochanek and White (2000); Yoshida et al. (2000b); Dave et al. (2001); Colin et al. (2002); Vogelsberger et al. (2012); Rocha et al. (2013); Peter et al. (2013); Zavala et al. (2013); Elbert et al. (2015); Vogelsberger et al. (2014); Fry et al. (2015); Dooley et al. (2016). In this case, $\lambda_{\rm mfp}$ is larger than the typical $\mathcal{O}({\rm kpc})$ core radius over which self-interactions are active. Here we discuss these simulations and their implications for astrophysical observables.

The majority of simulations assume a contact-type interaction where scattering is isotropic and velocity-independent, described by a fixed $\sigma/m$. These studies have converged on $\sigma/m\approx 0.5-1\;{\rm cm^{2}/g}$ to solve the core-cusp and TBTF issues on small scales, while remaining approximately consistent with other astrophysical constraints on larger scales, such as ellipticity measurements.111111This value is consistent with Spergel & Steinhardt’s original estimate, $\sigma/m\sim 0.45-450\;{\rm cm^{2}/g}$, based on having $\lambda_{\rm mfp}=(\rho_{\rm dm}\sigma/m)^{-1}$ in the range $1\;{\rm kpc}-1\;{\rm Mpc}$ with DM density $\rho_{\rm dm}\approx 0.4\;{\rm GeV/cm^{3}}$ as in the local solar neighborhood Spergel and Steinhardt (2000). However, more recent studies based on massive clusters disfavor the constant cross section solution, prefering somewhat smaller values $\sigma/m\approx 0.1\;{\rm cm^{2}/g}$ on these scales Kaplinghat et al. (2016); Elbert et al. (2016).

Before we discuss these topics in greater detail, we emphasize a few key points:

• •

Self-interaction cross sections are generically velocity-dependent, as predicted in many particle models for SIDM (see §VI). Therefore, constraints on $\sigma/m$ must be interpreted as a function of halo mass, since DM particles in more massive halos will have larger typical velocities for scattering. (“Halo mass” refers to the virial mass, for which we note that different studies have adopted somewhat different conventions in defining.)

• •

The effect from self-interactions on a halo is not a monotonic function of $\sigma/m$, which can be understood as follows. If DM is collisionless, its velocity dispersion in the halo is peaked near the scale radius $r_{s}$, while particles in both the center and outskirts of the halo are colder. Once collisions begin to occur, only the inner region of the halo is in thermal contact. The center of the halo forms a core that grows larger as more heat flows inward due to the positive temperature gradient. Eventually, once thermal contact is reached with the outskirts of the halo where the temperature gradient is negative, heat is lost outward and gravothermal collapse of the core ensues Lynden-Bell and Wood (1968); Colin et al. (2002).

• •

Our discussion mostly focuses on DM-only simulations and is subject to the usual caveat that baryons have been ignored. While baryonic feedback may explain away the issues that SIDM aims to address, this depends strongly on the prescription for feedback adopted. Without a more definitive understanding, making quantitative statements about the evidence for SIDM remains limited by these systematics. At the end of this section, we discuss recent simulations for SIDM including baryonic dynamics Vogelsberger et al. (2014); Fry et al. (2015).

• •

The present section focuses on the quasi-equilibrium structure of SIDM halos. We discuss simulations of merging halos, such as Bullet Cluster-like systems, in a later section (§V).

III.1 Implementing self-interactions in simulations

In N-body simulations, DM particles are represented by “macroparticles” of mass $m_{p}\gg{\rm\,M}_{\odot}$, each representing a phase space patch covering a vast number of individual DM particles. Self-interactions, assumed to be a short-range force on galactic scales, are treated using a Monte Carlo approach Burkert (2000). Two macroparticles scatter if a random number between $[0,1]$ is less than the local scattering probability within a given simulation time step $\Delta t$, which must be small enough to avoid multiple scatterings. Outgoing trajectories preserve energy and linear momentum121212Angular momentum is not conserved for individual scatterings, due to finite separations between macroparticles, but any nonconservation is expected to average to zero over the halo since separation vectors are randomly oriented Kochanek and White (2000)., with a scattering angle chosen randomly, under the assumption that scattering is isotropic in the center-of-mass frame.

Different simulations have adopted different prescriptions for determining the scattering probability. Many earlier studies used a background density method Burkert (2000); Kochanek and White (2000); Yoshida et al. (2000b); Colin et al. (2002). In this approach, the probability for an individual macroparticle $i$ to scatter is

[TABLE]

where $\rho_{{\rm dm}}^{(i)}$ is the mean local background density, spatially averaged over its nearest neighbors, and $v_{\rm rel}^{(i)}$ is the relative velocity between $i$ and one or more of its neighbors. If a scattering occurs, particle $i$ is paired up with a nearby recoiling partner $j$.

An alternative approach treats scattering as a pair-wise process between macroparticles $i,j$ with velocities $\mathbf{v}_{i,j}$ Dave et al. (2001); Rocha et al. (2013). The probability is

[TABLE]

where $\rho_{ij}$ represents the target density from $j$ to be scattered by $i$. Rocha et al. Rocha et al. (2013) provide an insightful derivation of Eq. (8) starting from the Boltzmann collision term. Each macroparticle $i$, centered at $\mathbf{r}_{i}$, is coarse-grained over a finite spatial patch using a cubic spline kernel $W(|\mathbf{r}-\mathbf{r}_{i}|,h)$ with smoothing length $h$ Monaghan and Lattanzio (1985).131313The smoothing kernel is defined as $W(r,h)=h^{-3}w(r/h)$, where

w(x)=\frac{8}{\pi}\left\{\begin{array}[]{cc}1-6x^{2}+6x^{3},&0\leq x\leq\tfrac{1}{2}\\ 2(1-x)^{3},&\tfrac{1}{2}<x\leq 1\\ 0,&x>1\end{array}\right..

(9) $W$ is normalized as $\int d^{3}r\,W(|\mathbf{r}-\mathbf{r}_{i}|,h)=1$ and has dimensions of number density. The resulting collision rate between patches $i,j$ is determined by the overlap integral

[TABLE]

Numerical convergence studies show that $h$ must not be too small compared to the mean particle spacing $a=(m_{p}/\rho_{\rm dm})^{1/3}$, requiring $h/a\gtrsim 0.2$ Rocha et al. (2013). Simulations in Ref. Rocha et al. (2013) have fixed $h$ such that $0.2\lesssim h/a\lesssim 1$ throughout the inner halo where the self-interaction rate is large. However, self-interactions are artificially quenched in the outer halo where $h/a<0.1$, although scattering is not expected to be relevant here due to simple rate estimates.

Alternatively, Vogelsberger et al. Vogelsberger et al. (2012) follow a hybrid between these approaches. First, $P_{ij}$ is computed in Eq. (8) with $\rho_{ij}=m_{p}W(|\mathbf{r}_{i}-\mathbf{r}_{j}|,h)$. Here, $\rho_{ij}$ can be interpretted as the “background” density from $j$ at the position of $i$ (which follows from Eq. (10) by setting the smoothing kernel for $i$ to be a delta function). The smoothing length is also taken to be much larger than the Ref. Rocha et al. (2013) approach, with $h$ adjusted dynamically such that $\sim 40$ particles are in range of $i$. Since $h\sim\sqrt[3]{40}\times a$, scatterings are less localized. The total probability for $i$ to scatter is determined by $P_{i}=\tfrac{1}{2}\sum_{j}P_{ij}$, and if a scattering occurs, the scattering partner $j$ is chosen with probabilty weighted by $P_{ij}$. It is unknown what differences, if any, may arise between these various methods (although both Refs. Rocha et al. (2013); Vogelsberger et al. (2012) are in mutual agreement with the Jeans modeling approach Kaplinghat et al. (2016) within 10–20%).

Next, we turn to convergence issues. Since self-interactions affect mainly the innermost radii in the halo, simulations must have sufficient resolution to robustly model these dynamics. It is well-known that coarse-graining introduces an artificial relaxation process due to two-body gravitational scattering that can alter the inner halo profile, with a timescale that scales with particle number, $t_{\rm relax}\propto N/\log N$ Binney and Tremaine (1987). Power et al. Power et al. (2003) showed that this process is the main limit for resolving the innermost halo structure for collisionless CDM. The innermost radius of convergence, termed the Power radius, occurs where there are too few enclosed particles, $N(r)$, such that $t_{\rm relax}(r)$ becomes shorter than the Hubble time. This radius $r$ is determined by the condition $N(r)/\log N(r)\approx 0.3\sqrt{\bar{\rho}_{\rm dm}(r)/\rho_{\rm crit}}$, where $\bar{\rho}_{\rm dm}(r)/\rho_{\rm crit}$ is the enclosed mean DM density contrast relative to the critical density Power et al. (2003).

In fact, SIDM halos are more robust and better converged below the Power radius compared to their CDM counterparts. This is sensible since the rate for self-interactions is larger—by construction, the goal is to have more than one self-interaction per Hubble time—than the rate for gravitational scattering Vogelsberger et al. (2012); Elbert et al. (2015).

At the same time, the cosmological environment at the outermost radii is also important for the SIDM halo. Its thermal evolution is affected by heating from mergers and infall, which slows the gravothermal collapse of cores. Hence, simulations of isolated halos are not sufficient for SIDM, as we discuss below.

III.2 Halo density profiles

Cosmological simulations of SIDM halos show a mass deficit at small radii compared to collisionless CDM halos, provided self-interactions are in the optically-thin regime Dave et al. (2001); Yoshida et al. (2000b). Here we discuss these results and the implications for $\sigma/m$.

Galactic scales: Davé et al. Dave et al. (2001) performed the first cosmological SIDM simulations targeting dwarf galaxies ($M_{\rm halo}\sim 10^{10}\,{\rm\,M}_{\odot}$). Their simulation volume contained several dwarf galaxies, resolved down to $\mathcal{O}({\rm kpc})$ scales, simulated with $\sigma/m\approx 0.5$ and $5\;{\rm cm^{2}/g}$. Their results preferred $5\;{\rm cm^{2}/g}$ in order to produce galaxies with core densities $\rho_{0}\approx 0.02\;{\rm\,M}_{\odot}/{\rm pc}^{3}$, broadly consistent with observed galaxies Firmani et al. (2000), while $0.5\;{\rm cm^{2}/g}$ yielded densities a few times larger but still viable. Moreover, for their largest simulated halo ($M_{\rm halo}\sim 6\times 10^{11}\;{\rm\,M}_{\odot}$), no evidence of core collapse was seen with cross sections as large as $\sigma/m\approx 50\;{\rm cm^{2}/g}$.

More recently, Elbert et al. Elbert et al. (2015) simulated two (cosmological) dwarf galaxies, resolved down to $\mathcal{O}(100\;{\rm pc})$ scales and for many values of $\sigma/m$. Fig. 9 (left) shows their results for one dwarf (“Pippin”) and demonstrates that central density is not a monotonic function of $\sigma/m$. For $\sigma/m\lesssim 10\;{\rm cm^{2}/g}$, self-interactions are predominantly in the core-growth regime (heat flowing in), and the central density decreases with increasing $\sigma/m$. However, a larger self-interaction rate, $\sigma/m=50\;{\rm cm^{2}/g}$, leads to an increasing central density, indicating this halo has entered core collapse. Nevertheless, core collapse is mild. Density profiles with $\sigma/m=0.5-50\;{\rm cm^{2}/g}$, spanning two orders of magnitude, vary in their central densities by only a factor of $\sim 3$. Comparing with data for field dwarfs in the Local Group, Fig. 9 (right) shows that predicted SIDM rotation curves for $0.5-50\;{\rm cm^{2}/g}$ are consistent with the velocities and half-light radii inferred from several observed galaxies. This illustrates not only how SIDM affects both the core-cusp and TBTF problems simultaneously, but that $\sigma/m$ need not be fine-tuned to address these issues.

The conclusion from these studies is that $\sigma/m\gtrsim 0.5\;{\rm cm^{2}/g}$ can produce $\mathcal{O}({\rm kpc})$ cores needed to resolve dwarf-scale anomalies Elbert et al. (2015). However, the upper limit on $\sigma/m$ at these scales—due to core collapse producing a too-cuspy profile—remains unknown.

Cluster scales: Next, we turn to clusters ($M_{\rm halo}\sim 10^{14}-10^{15}\,{\rm\,M}_{\odot})$. The first cosmological simulations at these scales were performed by Yoshida et al. Yoshida et al. (2000b), which studied a single $10^{15}\,{\rm\,M}_{\odot}$ halo for $\sigma/m=0.1$, $1$, and $10\;{\rm cm^{2}/g}$. More recently, Rocha et al. Rocha et al. (2013) performed simulations targeting similar scales, but over much larger cosmological volume, for $\sigma/m=0.1$ and $1\;{\rm cm^{2}/g}$. The best-resolved halos in their volume span $10^{12}-10^{14}\,{\rm\,M}_{\odot}$. For $1\;{\rm cm^{2}/g}$, the central density profiles are clearly resolved for the Yoshida halo and for $\sim 50$ Rocha halos. On cluster scales, SIDM halos have $\mathcal{O}(100-200\;{\rm kpc})$ radius cores and central densities $\rho_{0}\sim{\rm few}\times 10^{-3}\;{\rm\,M}_{\odot}/{\rm pc}^{3}$. For $\sigma/m=0.1\;{\rm cm^{2}/g}$, the simulations lack sufficient resolution to fully resolve the cored inner halo, though $\mathcal{O}(30\;{\rm kpc})$ radius cores seem a reasonable estimate. For $\sigma/m=10\;{\rm cm^{2}/g}$, the Yoshida halo has a similar density profile compared to $1\;{\rm cm^{2}/g}$, although the former is considerably more spherical (ellipticity is discussed below).

It is important to note that SIDM halos exhibit variability in their structure. Within the Rocha et al. Rocha et al. (2013) halo sample, SIDM halos, with fixed $\sigma/m=1\;{\rm cm^{2}/g}$ and fixed $V_{\rm max}$, show an order-of-magnitude scatter in their central densities. The dwarf halo samples from Davé et al. Dave et al. (2001) show a similar scatter in central density, albeit with lower resolution. This variation reflects the different mass assembly histories for different halos.

Strong gravitational lensing data has been used to constrain the core size and density in clusters relevant for SIDM Firmani et al. (2000, 2001); Wyithe et al. (2001); Meneghetti et al. (2001). Mass modeling of cluster CL0024+1654 found a cored DM profile with radius $\sim 50$ kpc Tyson et al. (1998)—quite similar to the Yoshida $0.1\;{\rm cm^{2}/g}$ halo—which was interpretted as evidence for self-interactions at these scales Firmani et al. (2000, 2001). However, interpretation for this particular cluster is complicated by the fact that it has undergone a recent merger along the line of sight (see Ref. Umetsu et al. (2010) and references therein).

Meneghetti et al. Meneghetti et al. (2001) placed the strongest constraint on cluster cores by examining the ability of SIDM halos (specifically, the halo from Yoshida et al. Yoshida et al. (2000b, a)) to produce “extreme” strong lensing arcs, i.e., radial arcs or giant tangential arcs. Giant tangential arcs, with length-to-width ratio $l/w\geq 10$, are present in many lensing observations Wu and Hammer (1993), but simulated SIDM halos with $\sigma/m=1$ or $10\;{\rm cm^{2}/g}$ lack sufficient surface mass density, limiting arcs to $l/w\lesssim 3.5$. Observations of radial arcs pose a more severe constraint. Since they do not occur even for the SIDM halo with $0.1\;{\rm cm^{2}/g}$ (as well as for larger $\sigma/m$), Meneghetti et al. conclude $\sigma/m<0.1\;{\rm cm^{2}/g}$ on cluster scales. However, there are several caveats to keep in mind (as acknowledged in Meneghetti et al. (2001)). First, due to the aforementioned variation in SIDM density profiles, constraints based on a single simulated halo require caution. Second, and more importantly, the simulated SIDM halo does not include the baryonic density from a central galaxy. Recent studies of massive clusters by Newman et al. Newman et al. (2013b, a)—two of which exhibit radial arcs (MS2137-23 and A383)—demonstrate that an $\mathcal{O}(10\;{\rm kpc})$ radius DM core is consistent with lensing data provided the baryonic mass is included. The important point is that the total density profile is well-described by an NFW profile, which is well-known to permit radial arcs Bartelmann (1996), even though the DM density by itself may have a cored profile. Thus, we conclude that cluster cores are consistent with $0.1\;{\rm cm^{2}/g}$, but $1\;{\rm cm^{2}/g}$ is excluded. (We make these statements more precise in §IV).

Isolated vs cosmological simulations: Lastly, we note that simulations of isolated SIDM halos can evolve much differently than cosmological halos. Kochanek & White Kochanek and White (2000) found that cores in isolated halos can be short-lived, collapsing very soon after formation and evolving toward a steeper $\rho\propto r^{-2}$ profile. For example, an isolated halo similar to the Pippen halo with $50\;{\rm cm^{2}/g}$ is expected to form a core and recollapse all within $\sim 1\;{\rm Gyr}$, in contrast to results from Elbert et al. Elbert et al. (2015). On the other hand, isolated halos with smaller cross sections can have cores that persist over a Hubble time. This points toward the importance of cosmological infall for mitigating gravothermal collapse—especially for larger values of $\sigma/m$—since the influx of high entropy DM particles slows energy loss from the core. In fact, it is the occasional violent major merger, rather than the smoother continuous infall of smaller clumps, which is predominantly responsible for “resetting the clock” for an otherwise collapsing core Yoshida et al. (2000b); Colin et al. (2002). This further supports the fact that SIDM halos are expected to have scatter in their structure that reflects the stochastic halo assembly process Colin et al. (2002); Brinckmann et al. (2017).

In addition, core collapse for isolated halos is dependent on the initial density profile, which must be set by hand as an initial condition for non-cosmological simulations. For example, SIDM cores evolving from an initial Hernquist profile, as assumed in Ref. Kochanek and White (2000), collapse twice as quickly relative to an initial NFW profile Koda and Shapiro (2011).

III.3 Halo shapes: ellipticity

Self-interactions are expected to make DM halos more spherical compared to triaxial collisionless CDM halos, at least in the inner regions where the scattering rate is largest. In fact, halo shape observations of elliptical galaxies and clusters have provided some of the most stringent constraints on self-interactions that exist in the literature. These limits are based on the assumption that one scattering per particle is sufficient to substantially affect the ellipticity of a halo. Early simulations by Davé et al. Dave et al. (2001) supported this conclusion. For their largest and best-resolved halos ($\sim 10^{11-12}\,{\rm\,M}_{\odot}$), SIDM halos yield a minor-to-major axis ratio $c/a\gtrsim 0.85$ at radii where at least one scattering has occurred, while CDM halos have $c/a\approx 0.6-0.7$ at these radii. However, these conclusions have been revisited in light of recent simulations with larger halo masses and statistics, as we now discuss.

Cluster ellipticity: One of the strongest quoted limits on SIDM is due strong lensing measurements of cluster MS2137-23, based on its apparent ellipticity at distances down to $\sim 70$ kpc Mellier et al. (1993). Using these observations and assuming the halo should be spherical where $R_{\rm scat}^{-1}\lesssim 5$ Gyr, Miralda-Escudé obtained a limit $\sigma/m\lesssim 0.02\;{\rm cm^{2}/g}$ on cluster scales Miralda-Escude (2002), which utterly rules out a contact-type interaction for SIDM. However, it is crucial to clarify the extent to which these arguments are borne out in N-body simulations.

More recently, Peter et al. Peter et al. (2013) revisted these issues in detail. Their study (companion to Rocha et al. Rocha et al. (2013), discussed above) involves a sample of $\sim 550$ simulated SIDM halos in the range $10^{11}-10^{14}\,{\rm\,M}_{\odot}$ with $\sigma/m=0.03$, $0.1$, and $1\;{\rm cm^{2}/g}$. They conclude that the Miralda-Escudé limit at $0.02\;{\rm cm^{2}/g}$ is far overestimated, due to several reasons:

• •

For $0.03\;{\rm cm^{2}/g}$ (close to the Miralda-Escudé limit), ellipticities for SIDM and CDM halos were found to be virtually identical, at least down the inner $\mathcal{O}(10-20\;{\rm kpc})$ as limited by resolution. Although MS2137-23 is a factor of four larger than the largest Peter et al. halo, large deviations from CDM should be visible in this halo at $\sim 35$ kpc, which matches the same scattering rate in MS2137-23 at $\sim 70$ kpc, if Miralda-Escudé’s argument is correct.

• •

Even for larger cross sections, SIDM halos are not spherical at the radius where one or two scatterings have occured, but rather remain somewhat elliptical with a median value of $c/a\approx 0.7$ for halos in the range $10^{13}-10^{14}\,{\rm\,M}_{\odot}$. Brinckmann et al. Brinckmann et al. (2017) have largely corroborated these results for larger halo masses $\sim 10^{15}\,{\rm\,M}_{\odot}$ comparable to MS2137-23.

• •

Differing mass assembly histories lead to scatter in $c/a$ at the $\sim 10-20\%$ for individual halos, which complicates drawing a limit based on a single system.

• •

Lensing observables are sensitive to the projected mass density along the line of sight. Thus, the surface density at small projected radius includes contributions from the outer halo, which may remain elliptical since self-interactions are not efficient.

Peter et al. conclude that even $1\;{\rm cm^{2}/g}$ is not excluded by strong lensing for MS2137-23.141414This conclusion is further supported by the Jeans analysis in Sec. IV below, which combines lensing and stellar kinematics data to show MS2137-23 is consistent with $\sigma/m\approx 0.1\;{\rm cm^{2}/g}$. See Fig. 12 (right).

Statistical studies of ensembles of lensing observations, as opposed to single measurements, offer a more powerful probe of self-interactions. To illustrate the potential of this approach, Peter et al. Peter et al. (2013) consider a subset of five clusters from the Local Cluster Substructure Survey (LoCuSS) Richard et al. (2009), which are chosen based on having a parametrically similar surface density profile compared to their five most massive simulated SIDM halos. The distribution of ellipticities for SIDM halos with $0.1\;{\rm cm^{2}/g}$ is consistent with the LoCuSS observations, while those with $1\;{\rm cm^{2}/g}$ are not. Peter et al. conclude that these data tentatively suggest $\sigma/m\lesssim 1\;{\rm cm^{2}/g}$, albeit with caution due to the limited statistics, unknown selection bias, and lack of baryons included in the SIDM lens modeling. Indeed, for $0.1\;{\rm cm^{2}/g}$, noticible differences in ellipticity between SIDM and CDM halos become apparent at only small radii where the stellar density becomes important Peter et al. (2013); Brinckmann et al. (2017).

Elliptical galaxies: These massive systems offer the opportunity to constrain self-interactions in halos on $10^{12}-10^{13}\,{\rm\,M}_{\odot}$ scales. Constraints on SIDM have focused on the isolated elliptical galaxy NGC 720 Feng et al. (2010a, 2009); Lin et al. (2012), which has $M_{\rm halo}\sim 7\times 10^{12}\,{\rm\,M}_{\odot}$ Humphrey et al. (2006), based on X-ray shape measurements from the Chandra telescope Buote et al. (2002). Buote et al. Buote et al. (2002) have shown that the X-ray isophotes remain elliptical at least down to projected radii $\sim 5$ kpc. (At smaller radii, shape measurements suffer from systematic uncertainties due to point source subtraction.) Since the X-ray emissivity scales as $\rho_{\rm gas}^{2}$, where the gas density $\rho_{\rm gas}$ has a fairly steep radial dependence, emission along the line of sight is strongly weighted toward physical radii near the projected radius (in contrast with lensing measurements, as described above). To determine the shape of the DM density, Buote et al. treat the gas as a single isothermal component in hydrostatic equilibrium and model the total mass density with a spheroidal profile that is either NFW, Hernquist Hernquist (1990), or isothermal ($\rho_{\rm dm}\propto r^{-2}$). While the isothermal model provides the best $\chi^{2}$ fit, all models (whether prolate or oblate) produce a similar ellipticity $\epsilon\approx 0.4$.

Constraints from NGC 720 have been based on the assumption that one self-interaction per particle within the inner 5 kpc over 10 Gyr is enough sphericalize the inner halo Feng et al. (2010a, 2009); Lin et al. (2012). This yields a stringent constraint, $\sigma/m\lesssim 0.01\;{\rm cm^{2}/g}$, assuming a mean DM density $\rho_{\rm dm}\approx 0.1\;{\rm\,M}_{\odot}/{\rm pc}^{3}$ and relative velocity $\rm v_{\rm rel}\approx 540$ km/s Humphrey et al. (2006); Feng et al. (2010a).

However, the simulations by Peter et al. Peter et al. (2013) show such constraints to be substantially overestimated as well. As above, SIDM halos retain ellipticity even where $R_{\rm scat}\sim 0.1\;{\rm Gyr}^{-1}$. Based on their sample of simulated SIDM halos with mass $3-10\times 10^{12}\,{\rm\,M}_{\odot}$, the distribution of ellipticities for $0.1\;{\rm cm^{2}/g}$ is perfectly consistent with $\epsilon\approx 0.4$, and even $1\;{\rm cm^{2}/g}$ halos are marginally allowed. On the other hand, the mean central densities for SIDM halos with $1\;{\rm cm^{2}/g}$ are typically too small compared to NGC 720 (by a factor of a few), while $0.1\;{\rm cm^{2}/g}$ provides better agreement. But before any robust limit can be made, it is essential to consider the gravitational influence of baryons on the DM density Kaplinghat et al. (2014a). If baryons dominate the total mass density within the inner $\sim 5-10$ kpc (as in the case for several elliptical galaxies that have received detailed mass modeling Humphrey et al. (2006)), the SIDM density profile can be steeper and more elliptical than expected without baryons.

III.4 Substructure

In the collisionless CDM paradigm, hierarchical structure formation leads to an abundance of substructure within halos Kauffmann et al. (1993). Observationally, these “halos within halos” manifest as dwarf galaxies around larger galaxies, or galaxies within groups or clusters. Self-interactions tend to erase substructure, provided the scattering rate is sufficiently large, through two effects Spergel and Steinhardt (2000). First, self-interactions within subhalos lead to density profiles that are less concentrated and more prone to tidal disruption. Second, self-interactions lead to evaporation of subhalos via ram pressure stripping as they pass through their host halo, since the former have a much lower velocity dispersion compared to the latter. Suppressing the subhalo mass function on MW scales is relevant for addressing the missing satellites problem Moore et al. (1999b); Klypin et al. (1999), while reducing the central density profiles of subhalos is relevant for addressing the TBTF problem Boylan-Kolchin et al. (2011, 2012). However, substructure on cluster scales must be preserved Moore et al. (1999b); Gnedin and Ostriker (2001).

Milky Way substructure: Early simulations by Davé et al. Dave et al. (2001) reported modest reductions ($30-50\%$) in the subhalo mass function below $10^{9}\,{\rm\,M}_{\odot}$ for $\sigma/m\approx 0.5$ and $5\;{\rm cm^{2}/g}$, which could help, but not alleviate, the missing satellites problem. However, more recent, higher resolution SIDM simulations have reached a more pessimistic conclusion Colin et al. (2002); Vogelsberger et al. (2012); Rocha et al. (2013); Zavala et al. (2013); Dooley et al. (2016). Fig. 10 shows the projected densities for a MW-like halo from Ref. Vogelsberger et al. (2012) for different DM models. For constant cross sections of $\sigma/m=1\;{\rm cm^{2}/g}$ or less, the subhalo mass function remains unchanged compared to collisionless CDM Zavala et al. (2013), which is illustrated in Fig. 10 (left). For a larger cross section of $\sigma/m=10\;{\rm cm^{2}/g}$, the subhalo mass function could be reduced by $\mathcal{O}(30\%)$ for subhalos below $10^{8.5}\,{\rm\,M}_{\odot}$, while the number of larger subhalos remains unaffected. This scenario corresponds to Fig. 10 (center). While it is evident that substructure is suppressed—particularly for smaller halos located nearer to the center of the host halo—even this modest reduction in substructure comes at a price of making the host halo spherical out to distances $\sim{\rm 50\;kpc}$. This scenario is excluded according to the ellipticity constraints discussed above (albeit for halos a factor of a few more massive than the MW).

Alternatively, self-interactions may be velocity-dependent, which effectively allows the cross section to vary across halo mass scales. Fig. 10 (right) shows a MW-like SIDM halo with velocity-dependent self-interactions Vogelsberger et al. (2012). The cross section chosen, motivated by classical scattering from a Yukawa potential Feng et al. (2010a); Loeb and Weiner (2011), is sizable on dwarf scales but suppressed on MW and larger scales. It is clear from Fig. 10 (right) that the subhalo mass function is indistinguishable from collisionless CDM. Even though the scattering rate is large on dwarf scales, this does not translate into evaporation of substructure since the relative velocity between the subhalo and its host halo is set by the velocity dispersion of the latter, for which the cross section is suppressed.

Despite being one of its original motivations Spergel and Steinhardt (2000), the conclusion is that SIDM cannot solve the missing satellites problem.151515This conclusion only applies to the minimal SIDM scenario where the only interaction is elastic scattering between DM particles. Nonminimal variations of SIDM can have a substantial impact on substructure (see §VII). Having a constant cross section around $0.5\;{\rm cm^{2}/g}$ is insufficient to erase MW substructure Colin et al. (2002); D’Onghia and Burkert (2003). Reducing substructure requires a much larger cross section on MW scales, which is excluded based on halo shape constraints, while merely having a large cross section on dwarf scales is insufficient to erase substructure Vogelsberger et al. (2012); Zavala et al. (2013). However, allowed SIDM models can still impact the stellar mass function for satellite galaxies Dooley et al. (2016). SIDM subhalos, while negligibly impacted by DM evaporation, have cored density profiles, resulting in a stellar density that is less tightly bound and more prone to tidal stripping.

Next, we turn to the TBTF problem. For collisionless CDM, the most massive subhalos within a simulated MW-like halo are too dense, with too-large predicted stellar velocity profiles, to match observed velocity dispersions for the MW dSphs Boylan-Kolchin et al. (2012). Despite having little effect on the abundance or total mass of subhalos, self-interactions can reduce their central densities and thereby reduce their velocity profiles in accord with observations, as shown in Fig. 11. Self-interactions at the level of $1\;{\rm cm^{2}/g}$ (or the velocity-dependent models considered therein) yield a reduced central density of $\rho_{0}\sim 0.1\;{\rm\,M}_{\odot}/{\rm pc}^{3}$ Vogelsberger et al. (2012); Zavala et al. (2013). These SIDM scenarios are consistent with the observed stellar kinematics for the MW dSphs Strigari et al. (2008); Wolf et al. (2010), while for $0.1\;{\rm cm^{2}/g}$ self-interactions are insufficient Zavala et al. (2013).

Cluster-scale substructure: In contrast to MW scales, there is no “missing galaxies” problem for clusters. Moore et al. Moore et al. (1999b) showed that the abundance of substructure within the Virgo cluster is well-described by simulations of a collisionless CDM halo with mass $5\times 10^{14}\,{\rm\,M}_{\odot}$. However, this effect does not provide a strong constraint on self-interactions. SIDM simulations with $1\;{\rm cm^{2}/g}$ find only a modest effect from subhalo evaporation on cluster scales, while observables related to the central cluster densities and core sizes are in principle much more stringent Rocha et al. (2013).

Using analytic estimates, Gnedin & Ostriker Gnedin and Ostriker (2001) excluded $\sigma/m$ between $0.3-10^{4}\;{\rm cm^{2}/g}$ based on the effect of DM evaporation on elliptical galaxies in cluster halos. As an elliptical galaxy passes through its host cluster halo, stars in the elliptical galaxy expand adiabatically as DM mass is lost through self-interactions, while the luminosity remains unchanged. This causes a shift in the fundamental plane for elliptical galaxies found in clusters and in the field, since the latter do not experience evaporation, and no significant environmental dependence is observed Kochanek et al. (2000).161616More recent studies have found a weak dependence to the fundamental plane intercept dependent on the local galaxy density. However, this dependence correlates with the stellar mass-to-light ratio, while the total dynamical mass (which would be affected by evaporation in denser environments) shows no such correlation La Barbera et al. (2010). However, based on simulation results for $1\;{\rm cm^{2}/g}$, Gnedin & Ostriker’s constraint may be overestimated, although detailed comparisons between N-body simulations and the fundamental plane relation are lacking Rocha et al. (2013).

III.5 SIDM simulations with baryons

While DM-only simulations have played a crucial role toward understanding the impact of self-interactions on DM halos, the next step is to incorporate self-interactions within hydrodynamic simulations that include the dynamics of baryons. Feedback processes remain the leading “vanilla” explanation for small scale structure issues apparent in CDM-only simulations, despite considerable debate over precisely how they are implemented. Hydrodynamic simulations with SIDM may provide guidance for how to disentangle and distinguish the effects of self-interactions and feedback. Additionally, it is important to understand how the predictions of SIDM are influenced by the presence of baryons, as well as conversely how baryonic tracers for DM are influenced by self-interactions.

Recently, Vogelsberger et al. Vogelsberger et al. (2014) and Fry et al. Fry et al. (2015) have performed the first N-body simulations including both self-interactions and baryons, both targeting dwarf scales. However, both sets of simulations have implemented star formation and feedback to opposite effect. In Ref. Fry et al. (2015), galaxies have bursty star formation histories (provided enough baryons are present) marked by episodic supernovae-driven gas outflows over times much shorter than the dynamical time scale of the galaxy. This process—motivated in part by explaining the formation of bulgeless dwarf galaxies by expelling low angular momentum gas that would otherwise form a stellar bulge—leads to nonadiabatic changes in the gravitational potential that can significantly impact the central DM distribution Governato et al. (2010); Pontzen and Governato (2012); Governato et al. (2012). On the other hand, Ref. Vogelsberger et al. (2014) has implemented a smoother prescription for star formation and feedback, that, while successfully suppressing bulge formation, provides a much smaller influence on the DM halo.

Vogelsberger et al. Vogelsberger et al. (2014) present results for two dwarf galaxies, simulated for a range of self-interaction cross sections. These galaxies are DM-dominated at all radii and are comparable to observed THINGS dwarfs Oh et al. (2011a) in terms of their stellar mass and maximum circular velocity, $V_{\rm max}\sim 50-80\;{\rm km/s}$. Whether DM is self-interacting or collisionless, baryonic feedback provides a negligible impact on the DM halo, and the DM central density is primarily influenced by the effect of self-interactions as in SIDM-only simulations. However, the presence of an $\mathcal{O}({\rm kpc})$ core in SIDM halos can affect the baryon component compared to a more concentrated DM profile, resulting in an $\mathcal{O}(30\%)$ reduction in star formation and reducing the central densities of stars and gas. In particular, the stellar distribution inherits a core radius that correlates with the core radius for DM, which may provide a useful observational handle for SIDM.

Fry et al. Fry et al. (2015) have simulated somewhat smaller systems, including three dwarf galaxies with $V_{\rm max}\sim 30-60\;{\rm km/s}$—comparable to field dwarfs of the Local Group Wolf et al. (2010)—as well as making predictions for even smaller dwarfs that are below current observational sensitivities. For the larger dwarfs, baryonic feedback by itself produces reduced central densities and $\mathcal{O}({\rm kpc})$ cores for collisionless CDM, upon which the further inclusion of self-interactions with $2\;{\rm cm^{2}/g}$ has only marginal effect. For smaller dwarfs with $V_{\rm max}\lesssim 20\;{\rm km/s}$, baryonic feedback is insufficient to form cores due to a lack of star formation. At the same time, self-interactions with $2\;{\rm cm^{2}/g}$ are unable to form cores larger than 500 pc in these smaller halos (consistent with simple rate arguments).

It is evident that the conclusions of Refs. Vogelsberger et al. (2014); Fry et al. (2015) are quite different due to the feedback prescriptions adopted therein. This underscores the role of baryonic physics as a systematic uncertainty in the small scale structure puzzle, as well the importance of taking new approaches to test the consistency of feedback models with collisionless CDM (e.g., Refs. Katz et al. (2016); Pace (2016); El-Badry et al. (2016, 2017)).

IV Jeans approach to relaxed SIDM halos

IV.1 Isothermal solutions to the Jeans equations

Despite the importance of N-body simulations, these methods are limited by their intensive computational nature. Even minimal particle models for SIDM exhibit rich dynamics for elastic scattering Tulin et al. (2013b, a), and it is not feasible to explore the full range of possibilities with simulations. To complement these studies, there is a useful semi-analytic method based on the Jeans equation for understanding SIDM halo profiles in relaxed systems Rocha et al. (2013); Kaplinghat et al. (2014a); Kaplinghat et al. (2016). This approach is well-suited to the intermediate cross section regime, $\sigma/m\sim 1\;{\rm cm^{2}/g}$, where DM is neither fully collisionless nor collisional throughout an entire halo. The Jeans method can be fit directly to observations for individual systems, from dwarf galaxies to massive clusters, including the gravitational effect on the halo from the observed baryonic mass distribution. Therefore, this method provides a bridge between simulations, astrophysical observations, and SIDM particle models.

The Jeans equation may be derived from the collisional Boltzmann equation (see, e.g., Binney and Tremaine (1987))

[TABLE]

where $f(\mathbf{r},\mathbf{v},t)$ is the DM distribution function, $\Phi_{\rm tot}=\Phi_{\rm dm}+\Phi_{b}$ is the total gravitational potential for both DM and baryons, and $\mathscr{C}$ is the collision term from self-interactions. In the inner halo, the collision term drives the distribution function toward kinetic equilibrium, $f\propto\exp(-\tfrac{1}{2}|{\mathbf{v}}|^{2}/\sigma_{0}^{2})$, where $\sigma_{0}$ is the isotropic one-dimensional velocity dispersion. Taking the first moment of Eq. (11) and searching for quasi-equilibrium solutions in which time derivatives can be neglected, we have the (time-independent) Jeans equation

[TABLE]

This is simply the condition for hydrostatic equilibrium for an ideal gas of pressure $p_{\rm dm}=\sigma_{0}^{2}\rho_{\rm dm}$. The velocity anisotropy is assumed to vanish since collisions isotropize trajectories. Moreover, N-body simulations have shown that DM particles are approximately isothermal within the inner halo, such that $\sigma_{0}$ can be taken as a constant.171717The velocity dispersion is expected to increase (decrease) with halo radius for heat transfer flowing in (out), corresponding the core growth (collapse) phase of the SIDM halo, depending on the value of $\sigma/m$. However, for a range of $\sigma/m$, N-body simulations have shown that the spatial variation in $\sigma_{0}$ is only $\mathcal{O}(10\%)$ within the inner halo Elbert et al. (2015). Combining with Poisson’s equation, we have

[TABLE]

where $\rho_{b}$ is the baryon mass density and $G$ is Newton’s constant.

On the other hand, DM is effectively collisionless in the outer halo due to the reduced particle number density. The delineation between the inner and outer halo occurs at radius $r_{1}$ where, on average, one collision per particle has occurred over the age of the halo, $t_{\rm age}\sim 5-10\;{\rm Gyr}$. This condition is $R(r_{1})t_{\rm age}\sim 1$, where the scattering rate is given in Eq. (1). The full density profile $\rho_{\rm dm}$ is taken to be a hybrid of collisional and collisionless (NFW) profiles, matched together at $r=r_{1}$:

[TABLE]

Here, $\rho_{\rm iso}$ is the isothermal density profile defined as the solution to Eq. (13). To match the two regions, the density profiles and enclosed mass are assumed to be equal at $r_{1}$. The physical picture is that self-interactions simply rearrange and thermalize the DM particle distribution in the inner halo, while leaving the outer halo as it would be in the absence of collisions.

With this simple halo model in mind, we turn to several important questions: Do self-interactions provide a consistent solution to the core-cusp problem in all astrophysical systems? What are the particle physics implications for observations spanning widely different halo mass scales, from dwarf galaxies to clusters? And lastly, does the Jeans approach for SIDM halos agree with N-body simulations?

To address these questions, Kaplinghat, Tulin & Yu Kaplinghat et al. (2016) applied the Jeans method to an astrophysical data set spanning DM halo masses in the range $\sim 10^{9}-10^{15}\;{\rm\,M}_{\odot}$. These data include rotation curves from twelve DM-dominated galaxies, including those from THINGS Oh et al. (2011) and LSB galaxies from Kuzio de Naray et al. Kuzio de Naray et al. (2008), combined with kinematic and lensing studies from six massive clusters by Newman et al. Newman et al. (2013b, a), all of which exhibit cored profiles. For each system, the cross section is obtained by

[TABLE]

following Eq. (1). Here, the quantity on the left-hand side represents the velocity-weighted cross section per unit mass, statistically averaged over velocity, while the $\rho_{\rm dm}(r_{1})$ on right-hand side is obtained by fitting Eq. (14) to astrophysical data for each system. Since more massive halos correlate with higher average relative velocities $\langle v_{\rm rel}\rangle$ for DM particles, the range of halos provides an important probe of the velocity-dependence of self-interactions. Analogous to tuning the beam energy in a particle collider, the energy-dependence of scattering is crucially important for probing the underlying particle physics of SIDM.

To illustrate the Jeans method, Fig. 12 shows results for two clusters from Ref. Kaplinghat et al. (2016). The full SIDM profile has been fit to the stellar velocity dispersions for the brightest central galaxy at small radii ($\lesssim 10$ kpc) and strong and weak lensing data at larger radii ($\gtrsim 10$ kpc). These data prefer cluster profiles with cores. A cuspy (NFW) profile fit only from lensing data does not agree with stellar data in the central halo (see inset). By matching the collisional and collisionless regions of the halo together to determine $r_{1}$, the preferred cross section for these clusters is $\sigma/m\approx 0.1\;{\rm cm^{2}/g}$. These conclusions may be weakened if stellar anisotropies are far more significant that assumed Schaller et al. (2015a) or if AGN feedback is relevant Martizzi et al. (2012), in which case $\sigma/m\lesssim 0.1\;{\rm cm^{2}/g}$.

Jointly analyzing both galaxies and clusters, Fig. 13 (left) illustrates how $\langle\sigma v_{\rm rel}\rangle/m$ depends on the average collision velocity $\langle v_{\rm rel}\rangle$, assuming self-interactions are responsible for the observed cores in these systems. While galaxy-scale observations favor $\sigma/m\approx 2\;{\rm cm^{2}/g}$, data from clusters prefers a much smaller cross section, $\sigma/m\approx 0.1\;{\rm cm^{2}/g}$ Kaplinghat et al. (2016). Taken at face value, these data imply that SIDM can provide a consistent solution to the core-cusp problem, provided self-interactions are relatively suppressed in clusters compared to dwarf galaxies. Such a behavior is well-motivated from a particle physics perspective, as discussed below. The data given here may be fit by a massive dark photon model (dashed orange curve in Fig. 13). Lastly, to verify the validity of the Jeans approach, Ref. Kaplinghat et al. (2016) analyzed mock rotation curves produced from eight SIDM halos in a similar mass range from N-body simulations Rocha et al. (2013); Elbert et al. (2015), reproducing the input cross section value $1\;{\rm cm^{2}/g}$ in those simulations.

Next, we turn to another important question: what is the interplay between self-interactions and baryons? The first simulations including both baryons with feedback processes and self-interactions for DM have only recently been performed (see §III.5), mainly targeting dwarf halos Vogelsberger et al. (2014); Fry et al. (2015). The Jeans method provides complementary insights into the effect of baryons on SIDM, especially in systems like the MW or larger that have a significant baryon fraction Kaplinghat et al. (2014a).

Although the Jeans approach is limited to quasi-equilibrium solutions (hence the dynamics of feedback is ignored), the static gravitational potential from baryons can dramatically change the predictions for observations compared to SIDM-only simulations. The baryon density enters through $\rho_{b}$ in Eq. (13), modifying the solution for $\rho_{\rm dm}$ from the usual cored isothermal profile. If $\rho_{b}$ dominates over $\rho_{\rm dm}$ in the inner halo, the core radius shrinks substantially compared to the SIDM-only halo without baryons. This effect is shown in Fig. 13 (right) for a MW-like halo and a self-interaction cross section of $\sigma/m\approx 0.5\;{\rm cm^{2}/g}$. While the SIDM-only halo has a core of size $\sim 5$ kpc, the core size is reduced by an order of magnitude due to the baryonic potential from stars. Except for the innermost $\sim 0.5$ kpc, the density profiles for both collisionless and collisional DM are virtually identical. Since $r_{1}\approx 8.5\;{\rm kpc}$ in this case, where $\rho_{b}$ and $\rho_{\rm dm}$ are comparable, adiabatic contraction Blumenthal et al. (1986) may modify the initial NFW profile matched to the inner isothermal profile, as indicated in Fig. 13 (right). However, it is negligible for the SIDM fits of the clusters shown in Fig. 12, because the contraction effect is very mild and it occurs $r\lesssim r_{1}$ Newman et al. (2013b).

The baryon density may affect the shape of the DM halo as well. While SIDM-only simulations predict halos that are spherical within the core radius, this conclusion changes once baryons are present. For a halo like the MW, the gravitational potential of the baryonic disk causes the SIDM halo to become oblate, with aspect ratio $\sim 0.6$, out to $\sim 5$ kpc, the would-be core radius for SIDM without baryons Kaplinghat et al. (2014a).

To summarize, recent studies using the Jeans method have challenged the conventional SIDM paradigm in several important ways.

• •

Astrophysical data from galaxies to clusters disfavor a constant self-interaction cross section. The velocity-dependence of scattering is important for constraining the particle physics underlying SIDM.

• •

Expectations for SIDM halo profiles can change dramatically from SIDM-only halos if there is a sizable baryonic component. The central regions of halos need not have large spherical cores with slope $d\log\rho_{\rm dm}/d\log r\approx 0$. Constraints based on halo shapes Peter et al. (2013); Feng et al. (2010a, 2009) need to be re-evaluated in light of the role of baryons in those systems.

In modeling the SIDM halo properties using the Jeans method, we have neglected dynamical effects, such as baryon feedback or environmental interactions, but we expect the result is robust to the galaxy formation history. Provided galaxy formation dynamics occurs interior to $r_{1}$, self-interactions will re-equilibrate the halo in response to any variations in the baryon distribution. For $\sigma/m\sim 1~{}{\rm cm^{2}/g}$, $r_{1}$ is close to the scale radius $r_{s}$, which is well outside typical radius for the stellar disk or bulge. It is hard to imagine that feedback processes would change the halo beyond $r_{1}$. Thus, the final SIDM density profile should be close to its equilibrium prediction, which depends on the baryon distribution and not the formation history.

IV.2 Diverse rotation curves for SIDM

The observed diversity of rotation curves poses a challenge for feedback solutions to the core-cusp problem. The issue is made clear in recent simulations by Oman et al. Oman et al. (2015), shown in in Fig. 14. Each panel shows the rotation curve (data points) for one of four galaxies. All four galaxies have comparable $V_{\rm max}\sim 80\;{\rm km/s}$, and therefore, within $\Lambda$CDM cosmology, are expected to have similar rotation curves. The colored bands show the scatter in rotation curves from hydrodynamical CDM simulations for similarly sized halos. Not only do the observations lie outside the expected range, they do so in different ways: IC 2574 has a slowly rising velocity indicating a large shallow core, while UGC 5721 has a steeply rising velocity due to a cuspy profile. Although this tension is somewhat dependent on the feedback prescription—Ref. Oman et al. (2015) adopt a smooth prescription—even bursty star formation cannot explain a core as large as IC 2574 (see Fig. 8).

On the other hand, DM self-interactions can explain this diversity. There are two effects at play. First, halos have intrinsic scatter due to variation in their mass assembly history. For CDM halos, this is reflected in the scatter in the mass-concentration relation. Since self-interactions rearrange the distribution of the halo interior to $r_{1}$, this scatter is reflected in the core size and central density for SIDM halos Kaplinghat et al. (2016); Kamada et al. (2016). However, this cannot be the whole story since the observed scatter is larger than can be explained by assembly history alone.

Second, the baryon distribution plays an important role in SIDM halos and thermalization can lead to very different halo profiles depending on whether the inner halo is dominated by DM or baryons. If DM-dominated, the density profile is $\rho_{\rm dm}\propto\exp(-\Phi_{\rm dm}/\sigma_{0}^{2})$, the solution to which is the nonsingular isothermal profile Eddington (1926), while if baryon-dominated, the halo profile is $\rho_{\rm dm}\propto\exp(-\Phi_{b}/\sigma_{0}^{2})$ and is largely set by the baryon density. In the latter case, the DM core size is controlled by the baryonic scale radius Kaplinghat et al. (2014a). Late-type galaxies span both cases and therefore scatter in the baryon profile translates into diverse DM profiles Kamada et al. (2016).

To investigate these effects, Kamada et al. Kamada et al. (2016) used the Jeans approach to model rotation curves for 30 spiral galaxies with asymptotic velocities $25-300\;{\rm km/s}$ and fixed $\sigma/m=3\;{\rm cm^{2}/g}$. The solution to Eq. (13) provides the density profile in the inner isothermal region. The baryon density is assumed to be a flat exponential disk

[TABLE]

where $\Sigma_{b}$ is the central surface density, $\Upsilon_{*}$ is the stellar mass-to-light ratio, $R_{d}$ is the disk scale radius, and $(R,z)$ are cylindrical coordinates. In the outer halo, the density profile is matched onto an NFW profile at $r_{1}$, according to Eq. (14), that is assumed to lie on the CDM mass-concentration relation within scatter Dutton and Macci (2014). In practice, the virial mass $M_{200}$ is fixed to match the asymptotic circular velocity, while $\Upsilon_{*}$ and the concentration parameter $c_{200}$ are determined to fit the inner rotation curve Kamada et al. (2016). According to the Jeans model, the isothermal halo parameters $\rho_{0},\sigma_{0}$ are determined once the outer NFW halo is fixed.

Fig. 15 shows how SIDM can accommodate rotation curve diversity. The four galaxies shown here are the same as those in Fig. 14. The solid red curves show total rotation curves that have been fit to these galaxies. It is quite remarkable that cosmological halos with fixed $\sigma/m$ can accommodate such disparate galaxies. For UGC 5721, the combination of a high concentration and dense baryon distribution yields a steep DM profile (solid blue) that is closer to NFW (dashed blue) than to a typical cored SIDM profile without baryons (stars). On the other hand, for IC 2574, a shallower core for SIDM can result from a low concentration $\sim 2.5\sigma$ below the median plus a diffuse baryon density. Despite having the same $V_{\rm max}$, SIDM can account for their wildly different inner velocities, $V_{\rm circ}({\rm 2\;kpc})\approx 70\;{\rm km/s}$ and $25\;{\rm km/s}$, respectively. These galaxies represent extremes of the scatter shown in Fig. 5 (right). Other galaxies in Ref. Kamada et al. (2016), such as those shown in Fig. 14, are intermediate to these extremes and can also be accommodated.

To further investigate these effects, Creasey et al. Creasey et al. (2017) performed a set of controlled simulations for isolated disk galaxies. In these simulations, the density profile of the initial halo is set up with a Hernquist profile, while the baryonic component is modeled as a fixed disk potential, and the disk mass is the sum of a stellar component and a gaseous one. The stellar mass is set by the abundance matching relation ($M_{\rm star}\textup{--}M_{200}$) in Behroozi et al. (2013) and the gas mass computed from the $M_{\rm gas}\textup{--}M_{200}$ relation proposed in Huang et al. (2012). Each halo—initialized to be in equilibrium with respect to the collisionless Boltzmann equation—is simulated with 18 different permutations: collisionless or $\sigma/m=2\;{\rm cm^{2}/g}$; concentration parameter at its median value or $\pm 2\sigma$ extreme values; and baryonic disk scale lengths $R_{d}=0.5$, $1.5$, and $6.0~{}{\rm kpc}$, representing ultra-compact, compact, and extended baryon distributions for observed galaxies, respectively. As shown in Fig. 16 (left), variations in the baryon distribution induce greater diversity in SIDM halos (solid magenta symbols) compared to CDM halos (open blue symbols). The latter has too little scatter to accommodate the diversity of observed rotation curves (green points).

The reason SIDM halos have greater diversity is as follows Creasey et al. (2017). Cored SIDM profiles have lower DM densities in their inner regions and the potential contribution of baryons is more important than for CDM, creating a larger range for $V_{\rm circ}({\rm 2\;kpc})$. For the compact disks ($R_{d}=0.5$ and $1.5~{}{\rm kpc}$), the predictions from SIDM and CDM are not very different, since the disk is the dominant contribution to $V_{\rm circ}({\rm 2\;kpc})$. For an extended disk, $R_{d}=6.0~{}{\rm kpc}$, $V_{\rm circ}({\rm 2\;kpc})$ in SIDM is lower than predicted in CDM and this difference becomes more prominent when $V_{\rm max}$ increases. Lastly, we note that there is a gap between compact ($R_{d}=1.5~{}{\rm kpc}$) and extended ($6~{}{\rm kpc}$) bayon disks in Fig. 16 (left). Intermediate baryon disks, which have not been included in these simulations, would populate this region.

IV.3 Scaling relations

In hierarchical structure formation, CDM halo parameters have a tight correlation known as the mass-concentration relation, which emerges from cosmological N-body simulations (see, e.g., Ref. Dutton and Macci (2014)). Since DM self-interactions only change the inner halo, while preserving the structure of the collisionless outer halo, we expect that the parameters of the inner halo—namely, the core density $\rho_{0}$ and core radius $r_{c}$—are also correlated for a given value of $\sigma/m$.

The right panel of Fig. 16 shows the surface density of the DM halo, given by the product $\rho_{0}r_{c}$, as a function of $V_{\rm max}$. Here, $\rho_{0}$ is the central density obtained from SIDM fits to 30 galaxies in Ref. Kamada et al. (2016) and $r_{c}$ is the radius at which the density is $\rho_{0}/2$. In the absence of baryons, the Jeans model predicts $\rho_{0}r_{c}\propto V^{0.7}_{\rm max}$ (dashed line), which can be traced to the mass-concentration relation Boyarsky et al. (2009, 2010); Lin and Loeb (2016). The “minimal” sample in Fig. 16 corresponds to galaxies where the stellar disk has negligible effect and is consistent with this relation. For the “moderate” and “maximal” samples, there is a large effect on the SIDM halo from baryons and these galaxies show a deviation from the simple scaling relation.

A number of other studies have shown that the DM halo surface density is nearly a constant Spano et al. (2008); Donato et al. (2009); Salucci et al. (2012); Kormendy and Freeman (2016); Burkert (2015); Hayashi and Chiba (2015); Karukes and Salucci (2016). These works adopt pseudo-isothermal or Burkert profiles to fit rotation curves (and stellar dispersions of dwarf spheroidals).181818In the case of a Burkert profile, the DM density at the core radius is one quarter of its central value. Although these studies do not impose the halo mass-concentration relation in their fits, the inferred DM surface density $\rho_{c}r_{c}\sim 100{\rm\,M}_{\odot}{\rm pc^{-2}}$ is consistent with the result shown here.

V Halo mergers

Merging systems are another avenue for exploring the collisional nature of DM. These systems include merging clusters, with the most famous being the Bullet Cluster Tucker et al. (1998); Markevitch et al. (2002, 2004); Clowe et al. (2006), as well as minor mergers from substructure (galactic or group halos) infalling into clusters Furlanetto and Loeb (2002); Natarajan et al. (2002); Massey et al. (2011). These systems probe the non-equilibrium dynamics of SIDM acting over timescales $\lesssim 1$ Gyr, which is complementary to the observables discussed above that test the quasi-equilibrium distribution of SIDM reached over $\sim$ 10 Gyr timescales.

Major mergers between clusters are characterized by three components:

• •

Galaxies: Since galaxies have negligible cross section for scattering, they act as collisionless test particles. (For minor mergers involving galactic-scale substructure, the stellar density plays the same role.)

• •

Intracluster medium (ICM) gas: Ram pressure during the merger causes ICM gas to get shocked and dissociated from any collisionless components, namely galaxies and (potentially) DM.

• •

DM: The total mass density, which is predominantly DM, is inferred through strong and weak gravitational lensing.

Although the initial state of the merger is not observed, X-ray observations of the shocked ICM, along with dynamical timing arguments Dawson (2013) and synchotron emission associated with turbulence (radio halos) or shocks (radio relics) Ensslin et al. (1998); Skillman et al. (2013), allow reconstruction of the three-dimensional geometry and velocity of the merger. If DM is collisionless, it should remain coincident with the galaxies after the merger, while if it is strongly self-interacting, it behaves qualitatively more like collisional gas. In the intermediate case in which the optical depth for self-scattering is nonzero but less than unity, self-interactions are expected to induce offsets between DM and galactic components due to collisional drag incurred between SIDM halos during pericenter passage Markevitch et al. (2004). The magnitudes of such offsets have been addressed in recent numerical studies Kahlhoefer et al. (2014); Robertson et al. (2017a); Kim et al. (2016); Robertson et al. (2017b), which we discuss below.

The same argument can be applied to minor mergers as well Massey et al. (2011); Harvey et al. (2014). Since substructure infall is the dominant mechanism for structure growth, these events are far more ubiquitous than cluster mergers. However, the expected offsets are much smaller: for collisionless DM, the expected gas-DM offsets are $\sim 20$ kpc (versus $\gtrsim 100$ kpc for major mergers) and clearly sensitivity below this scale is required to test SIDM Harvey et al. (2013, 2014). With the exception of rare systems with strong lensing data (such as Abell 3827), substructure merger studies must rely on weak lensing maps and require a stacked analysis of many clusters. Harvey et al. Harvey et al. (2014) proposed the “bulleticity ratio” $\beta=d_{\rm SD}/d_{\rm SG}$ as a suitable quantity to be averaged over many systems. Here, $d_{\rm SG}$ ($d_{\rm SD}$) is the projected distance between the stellar and gas (DM) peaks. The virtue of $\beta$ is that the orientation direction of the merger cancels in the ratio and that it has a relatively simple dependence on $\sigma/m$, under certain assumptions, as we discuss below. In subsequent work, Harvey et al. Harvey et al. (2015) determined the first constraints on self-interactions using a statistical analysis of $\beta$ from 30 merger systems. Their study is consistent with $\beta=0$, as expected for collisionless DM, resulting in the strongest merger-derived limit to date (however, see Ref. Wittman et al. (2017)).

V.1 Observations

The Bullet Cluster remains the archetypal dissociative merger. In this system, a smaller subcluster (the “Bullet”) has passed through a much larger main subcluster, evidenced by the dramatic bow shock left in the wake of the merger, shown in Fig. 17. From the orientation and sharp brightness edges of the shock, it is clear that the merger was a nearly head-on collision in the plane of the sky Markevitch et al. (2002). The centroids for the Bullet’s total mass and galactic distributions are offset by $25\pm 29$ kpc, implying that DM behaves as collisionless CDM within uncertainties Randall et al. (2008). Several hydrodynamical simulations have been performed with CDM plus gas in order to reconstruct the merger Springel and Farrar (2007); Mastropietro and Burkert (2008); Lage and Farrar (2014). Most recently, Lage & Farrar Lage and Farrar (2014) have simulated over a thousand Bullet Cluster realizations, including gas and magnetic fields, which were confronted against lensing, X-ray, and radio observations, as well as CMB measurements of the Sunyaev-Zel’dovich effect. In their best-fit model, the main and Bullet subclusters are consistent with triaxial NFW profiles of $2\times 10^{15}$ and $2\times 10^{14}\,{\rm\,M}_{\odot}$, respectively, and infall velocity consistent with $\Lambda$CDM cosmology Lage and Farrar (2014, 2015); Kraljic and Sarkar (2015). Constraints on self-interactions must be interpretted within the context of the large relative velocity between the subclusters, corresponding to $v_{\rm rel}\sim 4000$ km/s at centroid crossing Springel and Farrar (2007); Lage and Farrar (2015).

After the Bullet Cluster’s discovery Tucker et al. (1998), several other dissociative cluster mergers have been found. Those for which direct constraints on self-interactions have been quoted include DLSCL J0916.2+2951 (Musket Ball Cluster) Dawson et al. (2012); Dawson (2013), MACS J0025.4-1222 (Baby Bullet) Bradac et al. (2008), CIZA J2242.8+5301 (Sausage Cluster) van Weeren et al. (2011); Dawson et al. (2015); Jee et al. (2015), and ACT-CL J0102-4915 (El Gordo) Menanteau et al. (2012); Jee et al. (2014b); Ng et al. (2015), Abell 520 (Train Wreck Cluster) Markevitch et al. (2005); Mahdavi et al. (2007), and Abell 2744 (Pandora’s Cluster) Kempner and David (2004); Owers et al. (2011); Merten et al. (2011). Among substructure mergers, Abell 3827 has received considerable attention for its implications for SIDM Carrasco et al. (2010); Williams and Saha (2011); Massey et al. (2015); Kahlhoefer et al. (2015).

Musket Ball Cluster: This system is a binary merger between roughly equal mass subclusters ($2-3\times 10^{14}\,{\rm\,M}_{\odot}$) along a collision axis inclined $\sim 45^{\circ}$ to the plane of the sky Dawson et al. (2012); Dawson (2013). A feature of this merger is the large physical separation $\sim$ 1.3 Mpc between the subclusters, which implies that this system has evolved post-collision $2-5$ times longer than the Bullet Cluster Dawson et al. (2012). Unfortunately this does not directly translate into an enhanced sensitivity to self-interactions since, according to SIDM simulations, galaxy-DM offsets do not continue to grow with time after halo passage Kahlhoefer et al. (2014); Robertson et al. (2017a); Kim et al. (2016). One subcluster in this system appears to have an offset $\sim 130$ kpc between the weak lensing mass peak and its corresponding galactic centroid, with the former trailing the latter along the merger axis Dawson (2013). However, the statistical significance of this offset is not high (less than $2\sigma$) and does not exclude collisionless CDM Dawson et al. (2012); Dawson (2013).

Baby Bullet Cluster: Two subclusters with equal mass ($2.5\times 10^{14}\,{\rm\,M}_{\odot}$) have undergone a merger oriented in the plane of the sky. The mass lensing peaks for both subclusters are concident with their respective galactic luminosity centroids, implying that DM appears collisionless, while the peak of the gas density is offset from both subclusters and is located between them Bradac et al. (2008).

Sausage Cluster: This system is dominated by two massive subclusters ($10^{15}\,{\rm\,M}_{\odot}$) that have undergone a merger in the plane of the sky, as evidenced by the high polarization fraction of the prominent radio relics that trace the shock fronts van Weeren et al. (2011); Dawson et al. (2015); Jee et al. (2015). Weak lensing observations reveal two distinct mass peaks separated by $\sim 1$ Mpc, while the X-ray emission region is mainly located between the mass peaks and is highly elongated along the merger axis Dawson et al. (2015); Jee et al. (2015). Interestingly, the mass peaks for both subclusters lag behind their respective galactic luminosity centroids by $\sim 200$ kpc Jee et al. (2015). However, the statistical significance of these offsets is not high (below $2\sigma$) and becomes lower if number density (instead of luminosity) is used to weight the galactic centroids. Moreover, a number of systematic effects need to be clarified before this may be claimed as an effect due to SIDM, e.g., the total mass peaks getting skewed closer to the cluster center due to the gas density Jee et al. (2015).

El Gordo: This system is another binary merger between massive subclusters ($10^{15}\,{\rm\,M}_{\odot}$) Menanteau et al. (2012); Jee et al. (2014b), with a merger axis inclined $\sim 20^{\circ}$ with respect to the plane of the sky, inferred from the polarization fraction of the radio relics Ng et al. (2015). El Gordo, one of the most massive clusters found to date at redshifts $z>0.6$, has a prominent “bullet” of cool gas (similar to the Bullet Cluster) followed by twin-tailed wake Menanteau et al. (2012). The mass and luminosity peaks for each subcluster are coincident with one another, as expected for collisionless CDM. However, the gas bullet does not lag its corresponding mass peak as expected if the two subclusters are receding from one another after first pericenter passage Jee et al. (2014b). This may indicate that El Gordo lies on a return trajectory post apocenter Ng et al. (2015).

Abell 520: For this system, the interpretation of collisionless CDM is less clear. Like the Bullet Cluster, this system has a prominant bow shock indicating a recent, high-velocity merger Markevitch et al. (2005). However, the ICM gas has a much more complicated morphology Wang et al. (2016a), while weak lensing data reveals at least five distinct mass peaks Mahdavi et al. (2007). Hence, Abell 520 is no mere binary collision, but rather a “cosmic train wreck” of multiple simultaneous mergers Mahdavi et al. (2007); Okabe and Umetsu (2008). Of particular interest for SIDM is the observation of a dark core in Abell 520, i.e., a peak in the lensing mass map without a significant galactic counterpart Mahdavi et al. (2007); Okabe and Umetsu (2008); Girardi et al. (2008); Jee et al. (2012). However, the dark core was not confirmed in other lensing studies Clowe et al. (2012). More recently, Jee et al. Jee et al. (2014a) revisited this issue with the combined datasets used in Refs. Jee et al. (2012); Clowe et al. (2012), re-confirming the presence of the dark core albeit at a shifted location.191919The discrepancy with Ref. Clowe et al. (2012) was reconciled by systematic differences between lensing datasets, including source image coverage and differences in charge transfer inefficiencies in the CCDs used for these observations, which can skew galaxy ellipticities used for weak lensing studies (see, e.g., Ref. Rhodes et al. (2010)). Since the dark core lies coincident with the gas density, this suggests that DM (or a fraction thereof) may have self-interacted during the merger. Alternatively, the dark core may be a narrow DM filament oriented along the line of sight or a subcluster with an extreme mass-to-light ratio Mahdavi et al. (2007).

Abell 2744: This complex system appears to be a simultaneous merger between at least four clusters Owers et al. (2011); Merten et al. (2011). Previous analyses revealed a Pandora’s box of perplexing features in this system, such as a clump of X-ray-emitting gas leading a dark core Owers et al. (2011); Merten et al. (2011); Medezinski et al. (2016). With more recent observations, Jauzac et al. Jauzac et al. (2016) found eight mass peaks in Abell 2744 with masses $\sim 10^{14}\;{\rm\,M}_{\odot}$. With higher resolution mass maps, the purported dark core was fully resolved into two separate peaks coincident with galaxies, while the gas clump was seen to be associated with a previously unresolved DM halo in the process of merging Jauzac et al. (2016). Aside from being a remarkable coincidence of mergers Schwinn et al. (2017), Abell 2744 is consistent with CDM. Constraints on self-interactions have been obtained for two of the larger subclusters, which have undergone a Bullet Cluster-like merger, albeit significantly inclined along the line of sight Merten et al. (2011); Jauzac et al. (2016). For both subclusters, the total mass and galactic distributions remain coincident within uncertainties, while the gas is stripped, implying that DM appears collisionless.

Abell 3827: This remarkable cluster is a nearby late-stage merger of four elliptical galaxies, all located within $\sim 10$ kpc of the cluster core Carrasco et al. (2010). The fortutious proximity of strong lensing images near these substructures allows for a detailed reconstruction of the mass sub-peaks. Interestingly, the best resolved sub-peak shows a significant offset from its stellar counterpart, perhaps explained by drag due to self-interactions Williams and Saha (2011); Massey et al. (2015); Taylor et al. (2017). The size of the offset, $1.6\pm 0.5$ kpc Massey et al. (2015), is inconsistent with offsets seen in hydrodynamical simulations for collisionless CDM at $>99\%$ C.L (statistics only) Schaller et al. (2015b).202020After this report was completed, a new study was released for Abell 3827 that erased this discrepancy with CDM Harvey et al. (2017). Improvements in foreground subtraction and image identification allow for an improved lensing mass map. The position of the sub-peak is now only $0.54^{+0.22}_{-0.23}\;{\rm kpc}$ offset from its stellar density.

V.2 Self-interactions in merging clusters

To place constraints on self-interactions, merger studies have relied on three approaches: scattering depth, mass loss, and offsets. For clarity, we focus our discussion on the Bullet Cluster, to which all three methods have been applied. Table 2 summarizes constraints on self-interactions for all mergers quoted in the literature. However, as we now discuss, not all limits are equally robust.

The first (and most conservative) approach is based on the optical depth argument. For the Bullet Cluster, self-interactions must not be optically thick since the Bullet halo has survived the merger. The scattering depth $\tau$ must satisfy

[TABLE]

Here, $\Sigma_{\rm dm}\approx 0.3\;{\rm g/cm^{2}}$ is the peak projected mass density along the line-of-sight, which is assumed to be the same as the column density of the larger main halo along the trajectory of the Bullet Markevitch et al. (2004). For $\tau>1$, the Bullet halo would have interacted more like a fluid, experiencing similar drag and stripping as the gas, which implies $\sigma/m\lesssim 3\;{\rm cm^{2}/g}$.

The second method for constraining SIDM is based on mass loss Markevitch et al. (2004). For the Bullet Cluster, both subclusters have similar mass-to-light ratios, which is consistent with general expectations for clusters (e.g., see Ref. Carlberg et al. (1997)). Assuming both subclusters began with equal mass-to-light ratios prior to merger, present observations require that the Bullet halo could not have lost more than $\sim 23\%$ of its initial DM mass within its innermost 150 kpc at 68% CL Randall et al. (2008). SIDM simulations by Randall et al. Randall et al. (2008) show that $\sigma/m<0.7\;{\rm cm^{2}/g}$ is required to satisfy this constraint. While this is the strongest quoted bound from the Bullet Cluster, it relies on theoretical priors for the unobserved initial condition of the subclusters, as well as the merger itself not substantially affecting the star formation rate, which may not be the case (see, e.g., Ref. Sobral et al. (2015)).

The third method for probing self-interactions is based on offsets. Markevitch et al. Markevitch et al. (2004) proposed that self-interactions can lead to an effective drag force for the Bullet halo as it traverses the main halo, with DM particles losing momentum through scattering. This causes an apparent offset between the galactic and DM centroids, provided the drag force is sufficient to exceed the gravitational restoring force attracting the components together. From the measured null offset in the Bullet subcluster, Randall et al. Randall et al. (2008) obtained a limit $\sigma/m<1.25\;{\rm cm^{2}/g}$ at 68% CL according to their simulations.

In principle, offsets provide more robust constraints on SIDM free from theoretical priors for mass-to-light ratios and with sensitivity to smaller $\sigma/m$ outside the optically thick regime. Consequently, this method has received much theoretical attention, including N-body simulations Randall et al. (2008); Robertson et al. (2017a); Kim et al. (2016); Robertson et al. (2017b), analytic methods based on an effective drag force description for self-interactions Harvey et al. (2014), and hybrid approaches Kahlhoefer et al. (2014).

The offset effect depends on the type of self-interaction assumed Kahlhoefer et al. (2014); Robertson et al. (2017b). Due to the strong directionality inherent for mergers, it is important to distinguish between isotropic hard-sphere scattering and long-range interactions where scattering is forward-peaked (a la Rutherford scattering).212121Long-range does not necessarily mean macroscopic over galactic distances, but rather long-range with respect to the de Broglie wavelength of DM particles. As discussed in §VI, SIDM models of this type can have self-interactions through light force mediators with MeV-scale masses, which are clearly microscopic compared to galactic scales. The drag force description generally applies only for long-range interactions, in which the accumulation of many small-angle collisions retards the Bullet halo as it passes through the main halo.

On the other hand, the drag force description can break down for hard-sphere scattering, especially for large mass ratio systems like the Bullet Cluster where the merger velocity exceeds the escape velocity of the smaller halo. Large-angle collisions tend to eject particles from the Bullet halo, causing a centroid shift due to the expulsive tail of backward-going particles. (The tail also exerts a gravitational pull that slows down the remaining Bullet halo. However, since this effect is gravitational, it is the same for the galactic component and does not lead to a net offset Kahlhoefer et al. (2014).) For mergers with $\mathcal{O}(1)$ mass ratios, offsets can arise through a combination of both drag force and expulsive tail effects Kim et al. (2016).

For both types of interactions, the Bullet Cluster offset constraint from Randall et al. Randall et al. (2008) appears to be overestimated Kahlhoefer et al. (2014); Robertson et al. (2017a). For hard sphere scattering, Randall et al. Randall et al. (2008) find a 54 kpc offset for $\sigma/m=1.25\;{\rm cm^{2}/g}$, which is excluded at 68% CL. However, recent simulations by Robertson et al. Robertson et al. (2017a) find a much smaller offset $\sim 20$ kpc for $2\;{\rm cm^{2}/g}$. The discrepancy arises because DM halo positions within simulations are highly sensitive to the method for how they are measured Kahlhoefer et al. (2014); Robertson et al. (2017a). Ref. Randall et al. (2008) used the “shrinking circles” method: the position of the Bullet halo is given as the mean position of all particles within a circle of a given radius, which is then repeated with iteratively smaller circles (down to 200 kpc radius) each centered at the previous mean position. While this type of method is often used for locating density peaks in N-body simulations Power et al. (2003), the centroid shift in the Bullet halo is largely due to the expulsive tail, not due to a shift in the density peak Kahlhoefer et al. (2014). The shrinking circles method appears less robust for this situation, generating far larger offsets and with sensitivity to the initial circle chosen, compared to more observationally-motivated methods for determining positions based on parametric fits to the projected density or lensing shear map Robertson et al. (2017a).

Somewhat larger offsets can arise for anisotropic scattering, e.g., forward-peaked long-range interactions. Recently, Robertson et al. Robertson et al. (2017b) have performed the first N-body simulations with anisotropic self-interactions. Comparing models with forward-peaked and isotropic angular dependencies that both give rise to the same core radius in isolated SIDM halos, the former case produces a $\sim 50\%$ larger galaxy-DM offset in a Bullet Cluster-like merger compared to the latter case. However, this enhancement is reduced if the anisotropic cross section is associated with a Rutherford-like velocity dependence that falls with $v_{\rm rel}$ Robertson et al. (2017b), as is the case in light mediator models Ackerman et al. (2009); Feng et al. (2009); Buckley and Fox (2010); Feng et al. (2010a); Loeb and Weiner (2011); Tulin et al. (2013b, a). Since smaller $v_{\rm rel}$ corresponds to pairs of DM particles that move opposite to their subhalos’ bulk velocity during core passage, the remaining unscattered particles will be moving faster than the bulk velocity, hence reducing any offset Robertson et al. (2017b).

In any case, it is unclear whether any SIDM model can generate offsets larger than $\sim 50-100$ kpc needed to reach present observational sensitivities. Such values require $\sigma/m$ so large that the scattering depth approaches unity, in which case self-interactions enter the optically thick regime and merging subhalos would coalesce on impact Kim et al. (2016). Alternatively, Kim et al. Kim et al. (2016) have proposed searching for misalignments between the DM centroids and the brightest cluster galaxy (BCG) in merger remnants. Due to the cores in SIDM halos, BCG misalignment and oscillation induced during the merger can persist for Gyrs, even after the cluster has otherwise relaxed. Sensitivities down to $0.1\;{\rm cm^{2}/g}$ may be achievable, although further simulations (particularly, including gas physics) are required for more definitive statements.

V.3 Self-interactions in minor mergers

The most conservative constraint is based on the survival of substructure in clusters Gnedin and Ostriker (2001). As discussed in §III, the subhalo mass function in clusters is not substantially affected by self-interactions at $1\;{\rm cm^{2}/g}$, at least outside the central part of the cluster Rocha et al. (2013).222222For Abell 3827, we estimate a conservative constraint based on scattering depth. Taking a Hernquist profile for the main halo Kahlhoefer et al. (2015), we find $\Sigma_{\rm dm}\approx 0.3\;{\rm g/cm^{2}}$ for radial infall of a subhalo to 15 kpc, implying $\sigma/m\lesssim 3\;{\rm cm^{2}/g}$. Consequently, recent substructure merger studies have focused on offsets.

Notably, Harvey et al. Harvey et al. (2015) performed a stacked analysis for the bulleticity ratio $\beta$, using observations for 72 substructures in 30 systems, including both major and minor mergers. If self-interactions cause DM subhalos to lag their stars/galaxies, the parallel component $\beta_{\parallel}$ along the merger axis will be nonzero and positive. On the other hand, the orthogonal component $\beta_{\perp}$ is expected to be zero and can be used as a measure of systematics. Harvey et al. Harvey et al. (2015) find $\beta_{\parallel}=-0.04\pm 0.07$ (and $\beta_{\perp}=-0.06\pm 0.07$), consistent with collisionless CDM. This result is interpreted as a constraint on SIDM using an analytic relation

[TABLE]

where $\sigma_{*}$ is the characteristic cross section per unit mass when the halo becomes optically thick Harvey et al. (2014). Eq. (18) is derived assuming (i) a drag force description of self-interactions, which is valid for forward-peaked cross sections, and (ii) that the offsets between stars, gas, and DM are all small ($\lesssim 30$ kpc), which is in line with collisionless CDM simulations for minor mergers Harvey et al. (2014). The resulting limit is $\sigma/m<0.47\;{\rm cm^{2}/g}$ at $95\%$ C.L. Harvey et al. (2015).

Recently, Wittman et al. Wittman et al. (2017) raised serious concerns about this result. It is important to note that the highest weighted systems in the Harvey et al. Harvey et al. (2015) dataset are those with the largest stellar-gas offsets, i.e., the major mergers between clusters.232323To put the stringent limit of Ref. Harvey et al. (2015) in context, we note that their sample is heavily weighted ($67\%$) by five dissociative mergers whose constraints are listed in Table 2 Wittman et al. (2017). Moreover, the two highest weighted systems, Abell 520 ($21\%$) and Abell 2744 ($16\%$), have the most complicated morphologies. Since these systems have far larger offsets than 30 kpc, Eq. (18) is simply not valid and limits must be interpreted within the context of simulations. Setting aside one’s qualms with Eq. (18), Wittman et al. Wittman et al. (2017) also found that many of the stellar-DM offsets in Ref. Harvey et al. (2015), which utilize only single-band imaging for their weak lensing maps and galaxy identification, are inconsistent with values from other more comprehensive studies, which have obtained more accurate positions using strong lensing and multiband imaging data. Accounting for these new offsets, as well as other issues, the constraint is relaxed to $\sigma/m\lesssim 2\;{\rm cm^{2}/g}$ Wittman et al. (2017).

Lastly, we turn to Abell 3827, a unique cluster in which strong lensing observations reveal a significant DM-stellar offset for one of its infalling elliptical galaxies Williams and Saha (2011); Massey et al. (2015); Taylor et al. (2017). According to simulations, this offset is consistent with a self-interaction cross section $\sim 1.5\;{\rm cm^{2}/g}$ for a contact interaction or $\sim 3\;{\rm cm^{2}/g}$ for a drag force interaction Kahlhoefer et al. (2015). Taken at face value, however, such large cross sections appear to be inconsistent with constraints from cluster core sizes Kaplinghat et al. (2016).

VI Particle physics models

SIDM provides a compelling solution to the long-standing issues in galactic systems, while keeping all the success of CDM on larger scales. The preferred value of the DM self-scattering cross section is $\sigma/m\sim 1~{}{\rm cm^{2}/g}$ in galaxies, which is much larger than the weak-scale cross section expected for a usual WIMP DM candidate. In this section, we discuss particle physics models for SIDM and show that astrophysical observations over different scales can provide complementary information on the particle nature of DM self-interactions.

VI.1 What cross section is relevant?

Before we delve into particle physics models, we note that there are subtleties in mapping the actual observational constraints to the particle physics parameters for a given model. Since numerical simulations are the primary tools in studying the effect of self-interactions on DM structure formation, we would like to map the particle parameters to the simulation results, i.e., the preferred values of $\sigma/m$ in the scattering probability defined in Eqs. (7) and (8). If the DM scattering process is isotropic, such as in the $s$-wave limit, the cross section is independent of the scattering angle, the mapping is straightforward. However, if DM self-scattering is mediated by a long-range Coulomb-like interaction, it becomes difficult. In this long-range interaction limit, the total scattering cross section can be enhanced dramatically due to the singularity in the forward scattering direction, but it does not capture all relevant physics in this case. Since the small-angle forward scatterings actually play little role in conducting heat and changing the inner halo properties accordingly, the enhancement is spurious.

On the simulation side, it is time-consuming to simulate a large $\sigma/m$, enhanced by the forward scatterings, because a small time step is required. In addition, simulations need to take into account the angular distribution of outgoing DM particles. Hence, it is useful to consider a proxy that properly captures important physics and is also computationally cheap. A commonly-used prescription is to assume scattering is isotropic and that the scattering probability is proportional to the so-called transfer cross section, defined as Mohapatra et al. (2002); Feng et al. (2010a); Buckley and Fox (2010); Ibe and Yu (2010); Loeb and Weiner (2011); Tulin et al. (2013b, a); Schutz and Slatyer (2015)

[TABLE]

where $\theta$ is the scattering angle in the center of mass frame. Since the longitudinal momentum transfer is $\Delta p=m_{\chi}v_{\rm rel}(1-\cos\theta)/2$, $\sigma_{T}$ estimates the average forward momentum lost in collisions. The transfer cross section has been used to model the long-range interactions in simulations Vogelsberger et al. (2012); Zavala et al. (2013); Vogelsberger et al. (2016).242424Tulin, Yu & Zurek Tulin et al. (2013a) showed that the DM scattering is almost isotropic in dwarf galaxies with the model parameters adopted in Refs. Vogelsberger et al. (2012); Zavala et al. (2013).

Tulin, Yu & Zurek Tulin et al. (2013a) (see also Cline et al. (2014a); Boddy et al. (2016)) suggested that the viscosity (conductivity) cross section is a better proxy in capturing the DM self-scattering effects,

[TABLE]

It emphasizes scattering in the perpendicular direction ($\theta=\pi/2$) and measures the rate of energy equalization Schultz et al. (2008). Since both forward and backward scatterings lead to little energy exchange, they do not contribute to thermalization of the inner halo. In addition, $\sigma_{V}$ regulates scatterings in both forward and backward directions at the same time and it can be unambiguously calculated even when the quantum interference between $t$- and $u$ channels becomes relevant for identical particles.

The difference between $\sigma_{T}$ and $\sigma_{V}$ in interpreting the observational constraints is not significant. In the $s$-wave limit, scattering is isotropic and $\sigma_{T}=\sigma_{\rm tot}=(3/2)\sigma_{V}$, where $\sigma_{\rm tot}$ is the total scattering cross section. In the case with higher partial waves, $\sigma_{T}$ differs from $\sigma_{V}$ by an ${\cal O}(1)$ factor for Yukawa interactions between distinguishable particles Tulin et al. (2013a). For identical particles, $\sigma_{T}$ may have a backward-scattering enhancement and $\sigma_{V}$ must be used, though it differs from the distinguishable case only by an ${\cal O}(1)$ factor. Given systematic uncertainties in astrophysical observations, we expect both $\sigma_{T}$ and $\sigma_{V}$ provide a good measure for DM self-interactions, and both of them are widely used in the literature.

Recently, Robertson et al. Robertson et al. (2017b) performed the first SIDM simulations with a full treatment of anisotropic scattering. Their results for core sizes in isolated cluster halos show that the simplified prescription of taking $\sigma_{T}$ with isotropic scattering agrees with a full angular-dependent treatment within $20\%$.252525Kahlhoefer et al. Kahlhoefer et al. (2014) advocated an alternative definition of the transfer cross section, $\sigma^{\prime}_{T}=\int(1-|\cos\theta|)(d\sigma/d\Omega)d\Omega$. In the $s$-wave case, $\sigma^{\prime}_{T}=\sigma_{T}/2$, while $\sigma^{\prime}_{T}\approx\sigma_{T}$ if large-angle scattering can be neglected. Robertson et al. Robertson et al. (2017b) found that $\sigma^{\prime}_{T}$ works better than $\sigma_{T}$ in reproducing the full simulation results. The comparison is made in a cluster-like halo with the Hernquist DM density profile as the initial condition. Robertson et al. Robertson et al. (2017b) also showed that it is important to take into account full angular dependence in simulating merging clusters because there is a preferred direction when clusters collide Kahlhoefer et al. (2014). It will be interesting to further compare how well $\sigma_{T}$ and $\sigma_{V}$ serve as effective parametrizations for anisotropic scattering models across different halo mass scales within a cosmological environment.

VI.2 Self-coupled scalar

The simplest SIDM model consists of a real scalar field $\varphi$, which is a suitable DM candidate provided it is cosmologically stable Bento et al. (2000); McDonald (2002). Self-interactions arise through the quartic coupling

[TABLE]

and the self-scattering cross section is $\sigma(\varphi\varphi\to\varphi\varphi)=\lambda^{2}/(128\pi m_{\varphi}^{2})$, where $\lambda$ is the self-coupling constant and $m_{\varphi}$ is the DM mass in this model. Writing

[TABLE]

it is clear that $m_{\varphi}$ must lie near the MeV scale for $\lambda$ of order unity and $\sigma/m\sim 1~{}{\rm cm^{2}/g}$, as required to create cores in dwarf galaxies.

In the early Universe, thermal production and pair annihilation to SM fermions can achieve the correct DM relic abundance for $\varphi$ if, e.g., it couples to the Higgs boson Bento et al. (2000); Burgess et al. (2001). However, this scenario is subject to strong constraints from precision Higgs decay measurements Khachatryan et al. (2016) since it contributes to invisible Higgs decay. Alternatively, the relic abundance can be set by “cannibalization” Carlson et al. (1992); Hochberg et al. (2015). If the Lagrangian has a $\varphi^{3}$ term, then the self-annihilation process $\varphi\varphi\varphi\to\varphi\varphi$ can occur in the early Universe and effectively deplete the $\varphi$ abundance, resulting in the correct relic density Hochberg et al. (2015). The DM candidate $\varphi$ could be a composite state, such as glueballs of non-Abelian gauge fields in the hidden sector Carlson et al. (1992). In this case, both $\varphi^{3}$ and $\varphi^{4}$ terms are present naturally.

The major shortcoming of this simple model is that it predicts a constant cross section over all scales. As we discussed in §IV, $\sigma/m\sim 1~{}{\rm cm^{2}/g}$ is favored in dwarf galaxies, but it must be $\sim 0.1~{}{\rm cm^{2}/g}$ (or less) on cluster scales Kaplinghat et al. (2016), otherwise DM self-interactions create a ${\cal O}(100)~{}{\rm kpc}$ density core in galaxy clusters, too large to be consistent with observations Newman et al. (2013b). Nevertheless, this model does illustrate a key ingredient in building a successful SIDM model: there must exist a low mass scale much below the weak scale.

VI.3 Light mediator models

Self-interactions can arise through the exchange of additional weakly-coupled states that interact with DM particles. Perhaps the best-motivated setup along these lines is if DM is charged under a spontaneously broken $U(1)$ gauge symmetry. In this case, DM stability is naturally explained by charge conservation. As a by-product, gauge boson exchange mediates self-interactions, analogous to Rutherford scattering Feng et al. (2010a); Buckley and Fox (2010); Ibe and Yu (2010); Loeb and Weiner (2011); Tulin et al. (2013b, a); Schutz and Slatyer (2015). Historically, this type of model had been motivated earlier by purported anomalies in indirect detection observations, such as the 511 keV $\gamma$-rays from the galactic bulge Boehm et al. (2004); Boehm and Fayet (2004), the PAMELA positron excess Arkani-Hamed et al. (2009); Pospelov and Ritz (2009); Pospelov et al. (2008), and hidden sector DM model building (e.g., Refs. Chen and Tye (2006); Feldman et al. (2007a, b); Feng and Kumar (2008); Feng et al. (2008); Chen (2009)). We also consider self-scattering mediated by scalar bosons as well.

The model is described by the following interaction

[TABLE]

where $\chi$ is the DM particle (assumed to be a fermion), $\phi$ is the mediator, and $g_{\chi}$ is the coupling constant. In the non-relativistic limit, self-interactions are described by the Yukawa potential

[TABLE]

The model parameters are the dark fine structure constant $\alpha_{\chi}\equiv g^{2}_{\chi}/4\pi$, the mediator mass $m_{\phi}$, and the DM mass $m_{\chi}$. For a vector mediator, the potential is attractive ($-$) for $\chi\bar{\chi}$ scattering and repulsive ($+$) for $\chi\chi$ and $\bar{\chi}\bar{\chi}$ scattering, while for a scalar mediator, the potential is purely attractive. More detailed models Boddy et al. (2014a); Ko and Tang (2014a); Kang (2015); Kainulainen et al. (2016); Wang et al. (2016b); Ma (2017a); Kitahara and Yamamoto (2017); Ma (2017b), as well as potentials arising from other types of mediators Bellazzini et al. (2013), have been considered as well.

In general, this framework admits a rich cosmological phenomenology, similar to electromagnetism in the visible sector. If the mediator is sufficiently light, bound state formation and dark recombination Wise and Zhang (2014); Petraki et al. (2014, 2015, 2017); Cirelli et al. (2016), delayed kinetic decoupling Feng et al. (2009); van den Aarssen et al. (2012a); Chu and Dasgupta (2014); Buckley et al. (2014); Tang (2016); Agrawal et al. (2017), dark acoustic oscillations Feng et al. (2009); Cyr-Racine and Sigurdson (2013); Pearce and Kusenko (2013); Buckley et al. (2014), and dissipation could be relevant Fan et al. (2013a, b).

Here, we follow Refs. Tulin et al. (2013b, a) and restrict our attention to the effect of elastic DM self-scattering. In the perturbative limit ($\alpha_{\chi}m_{\chi}/m_{\phi}\ll 1$), the Born differential cross section is

[TABLE]

In the limit of $m_{\phi}\gg m_{\chi}v_{\rm rel}$, scattering is a contact interaction and the cross section is independent of $v_{\rm rel}$, as given in Eq. (3). In the opposite limit when $m_{\phi}\ll m_{\chi}v_{\rm rel}$, it scales as $1/v_{\rm rel}^{4}$ à la Rutherford scattering Feng et al. (2009); Ackerman et al. (2009); Agrawal et al. (2016). Neither limit provides the mildly velocity-dependent cross section favored by observations if self-interactions are to solve small scale structure issues on galactic and cluster scales Kaplinghat et al. (2016). However, a small but finite mediator mass can provide the right velocity dependence. This requires the transition between contact and Rutherford limits to occur around $v_{\rm rel}\sim 300$ km/s (between dwarf and cluster scales), implying $m_{\phi}/m_{\chi}\sim v_{\rm rel}/c\sim 10^{-3}$. For instance, for DM with 10 GeV mass, the required mediator mass is $m_{\phi}\sim 10~{}{\rm MeV}$ and the corresponding range of the force is $\sim 20~{}{\rm fm}$. It is remarkable that a dark force with femtoscale range can affect the dynamics of galaxies.

Exploring the full parameter space of this model requires calculating $d\sigma/d\Omega$ beyond the pertubative limit, where multiple scattering—ladder diagrams in Feynman diagramatic language—become important. Buckley & Fox Buckley and Fox (2010) and Tulin, Yu & Zurek Tulin et al. (2013b, a) developed a numerical procedure using a nonrelativistic partial wave analysis to study the nonperturbative regime. In some regimes, analytical expressions are also available. Feng, Kaplinghat & Yu Feng et al. (2010a) introduced formulae for the transfer cross section in the Born regime and also in the classical regime, the latter based on an empirical formula originally developed for plasma collisions Khrapak et al. (2003); Khrapak (2014).262626In the plasma medium, the photon obtains an effective mass due to the Debye screening effect and ion-ion (electron-ion) interactions can be modeled with a repulsive (attractive) Yukawa potential. Feng, Kaplinghat & Yu Feng et al. (2010a) first applied the semi-analytical fitting formula from Khrapak et al. (2003) to study DM self-interactions. Vogelsberger et al. Vogelsberger et al. (2012) performed the first Yukawa SIDM simulations, based on the plasma formula introduced in Feng et al. (2010a) and its reparametrized form suggested by Loeb & Weiner Loeb and Weiner (2011). Ref. Tulin et al. (2013a) also found an analytical expression for $s$-wave scattering, valid in the nonperturbative regime, using the Hulthén potential as an approximation, while Braaten & Hammer Braaten and Hammer (2013) studied $s$-wave resonant scattering on more general grounds. See Ref. Tulin et al. (2013a) for a map of the full parameter space and summary of numerical methods and analytical formulae.

Even this simple SIDM model exhibits many possible velocity dependencies for the scattering cross section. Fig. 18 (left) shows the transfer cross section $\sigma_{T}/m$, as a function of relative velocity. Each curve represents a different choice of parameters $(\alpha_{\chi},m_{\chi},m_{\phi})$ for an attractive Yukawa potential. In these cases, $\sigma_{T}$ is suppressed on clusters scales and larger ($v_{\rm rel}\gtrsim 1500\,{\rm km/s}$), corresponding to the Rutherford limit. However, the dependence becomes much more complicated on galactic scales with the potential for quantum mechanical resonances (and anti-resonances). For example, for $\alpha_{\chi}m_{\chi}/m_{\phi}\approx 1.6\,n^{2}$ where $n$ is a positive integer, scattering has an $s$-wave resonance and $\sigma_{T}/m_{\chi}$ becomes enhanced at low velocity.

Since the velocity dependence of $\sigma_{T}/m$ is quite sensitive to the choice of the particle physics parameters, astrophysical observations on different scales can constrain or even discover these parameters. Kaplinghat, Tulin & Yu Kaplinghat et al. (2016) showed that SIDM with a repulsive Yukawa interaction could yield a suitable cross section inferred from both galaxy and cluster data; see Fig. 13 (left). Fig. 18 (right) illustrates the favored range of ($m_{\chi},~{}m_{\phi}$) for fixed coupling, chosen to be $\alpha_{\chi}\approx 1/137$. Colored bands denote the regions preferred by dwarf galaxies (red), LSB galaxies (blue), and galaxy clusters (green). Remarkably, there exists a closed common region (black contours) that points to DM mass of $\sim 15~{}{\rm GeV}$ and mediator mass of $\sim 17~{}{\rm MeV}$. The result is independent of whether the dark and visible sectors are coupled via interactions beyond gravity.

An interesting question is whether the mediator $\phi$ can be massless Feng et al. (2009); Ackerman et al. (2009); Agrawal et al. (2016). In the absence of recombination into bound states, scattering is described by the Coulomb potential. Since the scattering cross section scales as $1/v_{\rm rel}^{4}$, fixing $0.1~{}{\rm cm^{2}/g}$ in galaxy clusters Kaplinghat et al. (2016) leads to an enormous cross section $\sim 10^{6}~{}{\rm cm^{2}/g}$ on dwarf scales. Such large cross sections are expected to lead to gravothermal collapse of a dwarf halo. However, Agrawal et al. Agrawal et al. (2016) argue that core formation and collapse are inhibited for large cross sections since heat conduction is suppressed by the small mean free path, suggesting an approximate duality between strongly and weakly self-interacting regimes. Ref. Ahn and Shapiro (2003) predicted that a large cross section value $\sigma/m\sim 10^{4}\;{\rm cm^{2}/g}$ leads to similar cored profiles as a smaller cross section $\sigma/m\sim 1\;{\rm cm^{2}/g}$. On the other hand, cosmological simulations that had been performed in the fluid limit for SIDM yield singular isothermal halos that are steeper than CDM halos Moore et al. (2000); Yoshida et al. (2000a), disfavoring this picture. Nevertheless, further simulations are required to make these statements more quantitative. Lastly, another complication is that collective plasma effects due to dark electromagnetic fields may be important in affecting the halo Ackerman et al. (2009); Heikinheimo et al. (2015); Sepp et al. (2016).

In the early Universe, light mediator states can play an important role in setting the DM relic abundance through $\chi\bar{\chi}\rightarrow\phi\phi$ annihilation. For symmetric DM, the required annihilation cross section is $\langle\sigma_{\rm ann}v_{\rm rel}\rangle\approx 5\times 10^{-26}{\rm~{}cm^{3}/s}$ assuming the dark and SM sectors are thermalized. For asymmetric DM, the relic density is determined by a primordial DM asymmetry and the annihilation cross section has to be larger than this value to deplete the symmetric density (see Refs. Davoudiasl and Mohapatra (2012); Petraki and Volkas (2013); Zurek (2014) and references therein). These considerations imply $\alpha_{\chi}\gtrsim 4\times 10^{-5}\,({m_{\chi}/{\rm GeV}})$ for the vector mediator case. On the other hand, smaller couplings are viable if the two sectors are thermally decoupled and the SM has a higher temperature Feng et al. (2008). In addition, the SIDM abundance could also be set by non-thermal production mechanisms Bernal et al. (2016); Reece and Roxlo (2016).

VI.4 Strongly interacting dark matter

SIDM may be a composite state of a confining non-Abelian gauge theory in the dark sector Carlson et al. (1992); Faraggi and Pospelov (2002); Frandsen et al. (2011); Cline et al. (2014b); Boddy et al. (2014a, b); Soni and Zhang (2016); Forestell et al. (2016); Ko and Tang (2017); Prilepina and Tsai (2016). Since such theories are already known to exist in nature—namely, quantum chromodynamics (QCD)—it is appealing that similar physics may exist for DM. QCD enforces that the proton is long-lived (due to an accidental symmetry) and accounts for most of the mass of visible matter in the Universe. Similarly, DM may be a dark hadron whose mass arises through nonperturbative physics and whose stability is imposed automatically by symmetry. Proposed candidates typically fall into the categories of dark baryons, dark mesons, or dark glueballs. However, the gauge group and properties of the dark quarks (number of flavors, masses, representations) are unknown a priori and many possible theories exist (see Ref. Kribs and Neil (2016) for a review).

Due to the nonperturbative nature of the theory, it is nontrivial to compute the self-interaction cross section and mass spectrum. In the low-energy limit, the self-interaction cross section can be expressed as $\sigma=4\pi a^{2}$ where $a$ is the scattering length. On dimensional grounds, one estimates $a\sim\Lambda_{\rm DM}^{-1}$, where $\Lambda_{\rm DM}$ is the dark confinement scale (in analogy with $\Lambda_{\rm QCD}\approx 300$ MeV). The DM mass $m$ depends in some detail on the nature of the DM state. For instance, for $SU(N)$ gauge theory, dark baryons have mass $m\sim N\Lambda_{\rm DM}$ if composed of effectively massless constituents (nucleon-like), while $m$ can be larger if one or more constituents are massive (like a heavy flavor baryon). If DM is a meson-like state, $m$ is somewhat arbitrary and DM can be lighter than $\Lambda_{\rm DM}$ (pion-like) or heavier (like a heavy flavor meson). For dark glueballs, the mass of the lightest state typically scales as $m\sim{\rm few}\times\Lambda_{\rm DM}$ Ochs (2013).

The cross section per unit mass for self-interactions can be expressed as follows:

[TABLE]

If all dimensionful scales are set by $\Lambda_{\rm DM}$, then taking $\Lambda_{\rm DM}\approx\mathcal{O}(0.2-0.6)\times\Lambda_{\rm QCD}$ gives a cross section in the right ballpark for addressing small scale issues on dwarf scales. However, in this simple picture, the nonrelativistic cross section is constant with velocity, which is disfavored by cluster bounds.

If self-interactions are the solution to small scale issues, then a velocity-dependent cross section is required. The cross section is expected to fall with velocity once the de Broglie wavelength becomes smaller than the scattering length, i.e., when $mv_{\rm rel}\gtrsim a^{-1}$. Clearly if both $m,a^{-1}\sim\Lambda_{\rm DM}$, this condition can never be satisfied for nonrelativistic scattering. On the other hand, if

[TABLE]

it may be possible to achieve the desired cross section on both galactic and cluster scales.

There are two simple ways for a composite DM model to satisfy Eq. (27). First, the constituents (e.g., dark quarks) may be very heavy such that $m\gg\Lambda_{\rm DM}$. As an example along these lines, Boddy et al. Boddy et al. (2014a) proposed a model consisting of heavy fermions in the adjoint representation of a confining $SU(N)$ gauge theory. DM particles consist of single fermions—with mass $\sim 10$ TeV to yield the correct relic density—plus dark gluons (“glueballinos”). The range of the interaction is set by $a^{-1}\sim\Lambda_{\rm DM}\sim 100$ MeV to give an acceptable self-scattering cross section.

Second, self-interactions may have a larger-than-expected scattering length $a\gg\Lambda_{\rm DM}^{-1}$. Precisely this effect occurs in QCD. Neutron-proton scattering has a resonant enhancement at low energy due to the weakly bound deuteron, such that $a^{-1}\approx 15\;{\rm MeV}\ll\Lambda_{\rm QCD}$. As shown in Fig. 3, the self-scattering cross section falls with velocity for $v_{\rm rel}\gtrsim(m_{p}a)^{-1}c\approx 5000\;{\rm km/s}$, while it is constant at low velocity (neglecting electromagnetism). The same physics may provide a viable SIDM candidate Cline et al. (2014b). However, to satisfy constraints on cluster scales ($v_{\rm rel}\sim 1500\;{\rm km/s}$), larger values of $ma$ are required so that the turn-over in the cross section is at lower velocities, ${\cal O}(100\;{\rm km/s})$. Other models with long range forces mediated by pseudo-Goldstone bosons may also be viable.

The observed abundance of composite SIDM particles can be realized in different ways. Since the gauge coupling is large, the usual freeze-out mechanism only works if $m\sim{\rm TeV}$. For example, in Ref. Boddy et al. (2014a), massive dark quarks annihilate into dark gluons before the dark confinement phase transition. On the other hand, if $m\ll{\rm TeV}$, other mechanisms must be utilized since the DM abundance would be too small in the usual freeze-out scenario. For DM carrying a conserved charge—analogous to baryon number—it is natural to assume that its abundance is set by a primordial asymmetry similar to the baryon asymmetry of the Universe (see, e.g., Davoudiasl and Mohapatra (2012); Petraki and Volkas (2013); Zurek (2014) and references therein). Here, a large annihilation cross section is important to deplete the symmetric density efficiently. One needs make sure that there are additional light degrees of freedom into which the symmetric component can annihilate, such as dark pions, which must subsequently get depleted. Dark pions may decay into SM fermions through massive mediators Cline et al. (2014b). Alternatively, if dark quarks are massless, dark pions are Goldstone bosons and get redshifted away as dark radiation, which remains allowed if the dark sector is slightly colder than the visible sector Feng et al. (2008); Cline et al. (2014b). In the case of glueball DM, $3\rightarrow 2$ cannibalization can achieve the correct DM abundance Carlson et al. (1992); Hochberg et al. (2015); Soni and Zhang (2016) and excited states may be important Forestell et al. (2016). A strongly-coupled SIDM sector could also be realized in the dynamical DM framework Dienes and Thomas (2012); Dienes et al. (2016). In this scenario, the DM density comprises an ensemble of different hadron-like states and the DM relic abundance evolves with time as some states decay away.

VI.5 Dark atoms

Composite SIDM candidates may arise in weakly coupled theories in the form of atomic (or molecular) bound states Mohapatra and Teplitz (2000); Mohapatra et al. (2002); Foot (2001); Khlopov (2006); Kaplan et al. (2010); Behbahani et al. (2011); Cyr-Racine and Sigurdson (2013); Cline et al. (2012); Fan et al. (2013a, b); Cline et al. (2014a); Foot (2014); Boddy et al. (2016). The simplest atomic DM model consists of a hydrogen-like bound state composed of two distinct species with opposite charges under a massless $U(1)$ gauge symmetry. Following the same nomenclature as for normal hydrogen, there is a dark electron $e$ and dark proton $p$, with masses $m_{e,p}$, and a dark fine structure constant $\alpha$. The DM particle is denoted $H$, with mass $m_{H}\approx m_{e}+m_{p}$. The effective potential between hydrogen atoms is known in the atomic physics literature and depends on the total electronic spin configuration (singlet or triplet) Silvera (1980). Refs. Cline et al. (2014a); Boddy et al. (2016) adopted these results to compute self-interactions between dark atoms. As an example, the potential for the triplet state is

[TABLE]

where $\mu_{H}=m_{e}m_{p}/m_{H}$ is the atomic reduced mass and the numerical factors are fitted from experimental atomic data Silvera (1980); Cline et al. (2014a). The complicated expression in brackets—the first term models repulsive exchange and the second term represents the Van der Waals interaction—is easily adapted to the dark sector by expressing Eq. (28) in terms of the radius $x=r/a_{0}$ normalized in units of the Bohr radius $a_{0}=(\alpha\mu_{H})^{-1}$.

To calculate the atomic self-scattering cross section, Cline et al. Cline et al. (2014a) used a partial wave analysis. Boddy et al. Boddy et al. (2016) further included inelastic scattering due to hyperfine transitions. At low energy, $s$-wave scattering dominates and the cross section scales with the geometric size of the dark atom, $\sigma\approx 100\,a^{2}_{0}$ Cline et al. (2014a). In the limit $m_{p}\gg m_{e}$, the cross section is numerically given by

[TABLE]

Dark atom scattering can produce a velocity-dependent scattering cross section, which is preferred for viable SIDM models. Fig. 19 (left) shows that the cross section decreases once the collision energy exceeds about one-tenth of the atomic binding energy, corresponding to the de Broglie wavelength becoming smaller than $a_{0}$ Cline et al. (2014a); Boddy et al. (2016). In this region, the total (“elastic”) cross section is larger than the viscosity (“transport”) cross section due to the spurious enhancement from forward and backward scattering, even although they are normalized to be identical in the low-energy limit. Boddy et al. Boddy et al. (2016) showed that atomic scattering can have the right velocity dependence to explain density cores inferred in dwarf galaxies while being consistent with observations on cluster scales. Fig. 19 (right) maps out the parameter space favored by these observations, which imposes $\alpha>0.01$ and the DM mass to be in the $10\textup{--}100~{}{\rm GeV}$ range.

There are several realizations of atomic DM. In the mirror DM model, the dark sector is a mirror copy of the SM sector Mohapatra and Teplitz (2000); Mohapatra et al. (2002); Foot and Vagnozzi (2015). In the double disk DM model Fan et al. (2013a), dark atoms only compose a very small fraction of DM, and they may form a dark disk in DM halos through dissipative processes. For atomic DM, the relic density should be set by a primordial asymmetry. Since $\alpha$ is typically large for atomic bound states to form in the early Universe, the annihilation cross section of baryons and anti-baryons to massless dark photons is significant enough to deplete the symmetric component. One important constraint on atomic DM is that the energy dissipation rate should be small enough such that DM halos can form Mohapatra et al. (2002); Cyr-Racine and Sigurdson (2013); Boddy et al. (2016), although Foot & Vagnozzi Foot and Vagnozzi (2015, 2016) argued that the energy injection from supernova explosions in the visible sector could compensate the energy loss due to the dissipative cooling process if the two sectors are connected with a kinetic mixing term, and both the core-cusp and the missing satellites problems could be resolved.

VI.6 SIDM with an excited state

Self-interactions need not be a purely elastic process. The dark sector may include a spectrum of closely lying states and self-interactions can mediate transitions between them if allowed kinematically. In the literature, DM models with inelastic transitions have been motived for other phenomenological studies, such as direct detection Han and Hempfling (1997); Tucker-Smith and Weiner (2001); Batell et al. (2009a); An et al. (2012); Bramante et al. (2016), indirect detection Finkbeiner and Weiner (2007); Arkani-Hamed et al. (2009); Tulin et al. (2013c); Finkbeiner and Weiner (2016), and collider signals Baumgart et al. (2009); Bai and Tait (2012). For SIDM, such a spectrum of states may arise as energy levels of atomic DM Cline et al. (2014c); Boddy et al. (2016) or strongly interacting composite DM Boddy et al. (2014b).

The minimal model along these lines consists of two states, the DM ground state $\chi$ and an excited state $\chi^{*}$, with an interaction that mediates transitions between them. If DM particles are charged under a $U(1)$ gauge symmetry, symmetry breaking can generate a mass term that splits the charged state into two real states $\chi$ and $\chi^{*}$ with a coupling to the gauge boson $\phi$ that is purely off-diagonal. In the non-relativistic limit, the potential is

[TABLE]

in the basis of ($\chi\chi$, $\chi^{*}\chi^{*}$) two-body scattering states, where $\delta$ is the mass splitting Slatyer (2010); Hisano et al. (2005). In general, all possible combinations of initial and final states need to be considered, i.e.,

[TABLE]

since multiple interactions are important outside the Born regime. Schutz & Slatyer Schutz and Slatyer (2015) derived analytic expressions for these cross sections including $s$-wave contributions only. Zhang Zhang (2017) used an adiabatic approximation to investigate parameters with a large mass splitting and small coupling, while Blennow, Clementz & Herrero-Garcia Blennow et al. (2017) obtained numerical solutions covering the full parameter range.

Fig. 20 (left) shows the elastic (blue) and inelastic (red) scattering cross sections as a function of the DM velocity, for the benchmark models presented in Ref. Blennow et al. (2017), where the coupling constant is fixed as $\alpha_{\chi}=0.01\,(m_{\chi}/270~{}{\rm GeV})$ from relic density. In all cases, the viscosity cross sections for $\chi\chi\rightarrow\chi\chi$ are larger than $\sim 1~{}{\rm cm^{2}/g}$ in dwarf halos, but decrease below $\sim 0.1~{}{\rm cm^{2}/g}$ in clusters. Inelastic scattering $\chi\chi\rightarrow\chi^{*}\chi^{*}$ only occurs in cluster halos, where the DM velocity is large enough to overcome the mass gap $2\delta$. Fig. 20 (right) illustrates the parameter regions, where the DM self-interactions create cores in dwarf galaxies. Compared to the elastic case, the inelastic model requires a smaller mediator mass ($m_{A^{\prime}}$) to compensate for the mass splitting.

Down-scattering $\chi^{*}\chi^{*}\rightarrow\chi\chi$ may also play an interesting role in heating a halo by converting excitation energy into kinetic energy. Typically, however, the abundance of $\chi^{*}$ is negligible in DM halos even if it is stable because down-scattering is in chemical equilibrium in the early Universe, which depletes its abundance Finkbeiner et al. (2009); Batell et al. (2009a); Blennow et al. (2017).

VII Complementary searches

The unifying feature of all SIDM models is the requirement of a new mass scale below the GeV-scale. For example, in weakly coupled theories, self-interactions arise through dark sector mediators whose mass must be sufficiently light to give a large enough cross section. In this section, we consider how complementary searches for light dark sectors can shed light on the physics of SIDM. Our discussion focuses on two general themes: (1) interactions between DM and SM particles through sub-GeV mediators and (2) cosmological probes of decoupled dark sectors that do not interact with the SM.

VII.1 Portals for light mediators

Self-interactions through the Yukawa potential provide one of the simplest viable frameworks for SIDM. Typically, the mediator $\phi$ has mass in the range of $\sim 1-100$ MeV to yield the correct self-scattering cross section. If $\phi$ was thermally produced in the early Universe, it must decay, otherwise its abundance would dominate the total energy density, leading to a cosmic overclosure problem Kaplinghat et al. (2014b); Boddy et al. (2014a). Kaplinghat, Tulin & Yu Kaplinghat et al. (2014b) suggested that the mediator may decay to SM particles if it couples to the SM sector. The same coupling, therefore, may lead to signals in DM direct and indirect detection searches, as sketched in Fig. 21, as well as collider searches.

There are several ways that $\phi$ may couple to the SM. Depending on its spin, $\phi$ can mix with photon or $Z$ boson if it is a vector, or it can mix with the Higgs boson if it is a scalar. This mixing allows relic $\phi$ particles to decay, typically into $e^{+}e^{-}$ and $\nu\bar{\nu}$ pairs. Here, we follow Ref. Kaplinghat et al. (2014b) and discuss these examples.

For the vector mediator case, mixing with photon or $Z$ boson is described by the Lagrangian

[TABLE]

where $\phi_{\mu}$ is the vector mediator field and $\phi_{\mu\nu}\equiv\partial_{\mu}\phi_{\nu}-\partial_{\nu}\phi_{\mu}$ is its field strength, and $F_{\mu\nu}$ is the photon field strength. The first term corresponds to photon kinetic mixing Holdom (1986), parametrized by $\varepsilon_{\gamma}$, which has been widely studied as a DM portal, e.g., Foot (2004); Feldman et al. (2007a); Pospelov and Ritz (2009); Arkani-Hamed et al. (2009). The second term corresponds to mass mixing with the $Z$ Babu et al. (1998); Davoudiasl et al. (2012, 2013), and it can be parametrized by $\varepsilon_{Z}\equiv\delta m^{2}/m_{Z}^{2}$, where $m_{Z}$ is the $Z$ boson mass. We assume the mixing parameters are small: $\varepsilon_{\gamma,Z}\ll 1$.

The mixing terms in Eq. (32) induce a coupling of $\phi$ to SM fermions

[TABLE]

where $J_{\rm em}^{\mu}=\sum_{f}Q_{f}\bar{f}\gamma^{\mu}f$ and $J_{\rm NC}^{\mu}=\sum_{f}\bar{f}\gamma^{\mu}(T_{3f}P_{L}-Q_{f}\sin^{2}\theta_{W})f$ are the usual electromagnetic and weak neutral currents, respectively. Additionally, $g_{2}$ is the SU($2$)L coupling, $\theta_{W}$ is the weak mixing angle, $P_{L}$ is the left-handed projection operator, and $Q_{f}$ and $T_{3f}$ denote the charge and weak isospin for SM fermion $f$, respectively.

If the dark mediator $\phi$ is a scalar, the leading renormalizable couplings to SM particles arise through the Higgs portal Patt and Wilczek (2006); March-Russell et al. (2008); Ahlers et al. (2008); Andreas et al. (2008); Arina et al. (2010); Chu et al. (2012); Djouadi et al. (2012); Bhattacherjee et al. (2013); Greljo et al. (2013); Bian et al. (2014); Choi et al. (2013); Kouvaris et al. (2015); Kahlhoefer et al. (2017). Assuming $\phi$ is a real scalar singlet, the relevant terms in the scalar potential are

[TABLE]

where $H$ is the Higgs doublet and $a,b$ are coupling constants. After electroweak symmetry breaking, mixing arises between $\phi$ and the physical Higgs boson $h$ due to the Higgs vacuum expectation value (vev) $v\approx 246$ GeV. In the limit $a,m_{\phi}\ll v,m_{h}$, this mixing angle is $\varepsilon_{h}\approx av/m_{h}^{2}$. This generates an effective $\phi$ coupling to SM fermions

[TABLE]

In the early Universe, $\phi$ decays must not affect Big Bang Nucleosynthesis (BBN) and dilute the baryon density, which puts a constraint on the mixing parameter Kaplinghat et al. (2014b). For example, in the kinetic mixing case ($\phi$ is known as a “dark photon”), the final states are $e^{+}e^{-}$ and the decay rate is $\Gamma_{\phi}=\alpha_{\rm em}m_{\phi}\varepsilon_{\gamma}^{2}/3$, resulting in a lifetime

[TABLE]

Kaplinghat, Tulin & Yu Kaplinghat et al. (2014b) estimated a lower bound on the mixing parameter, $\epsilon_{\gamma}\gtrsim 10^{-10}$, by assuming that the mediator decays before weak freeze-out, $\sim 1~{}{\rm s}$. Fradette et al. Fradette et al. (2014) studied BBN and CMB constraints on the $(m_{\phi},\varepsilon_{\gamma})$ plane with the assumption that the cosmic abundance of the mediator particle is via the freeze-in process, i.e., the inverse decay process $f\bar{f}\rightarrow\phi$. Since the mixing parameter controls both production and decay rates in this case, the constrained parameter regions are closed. The injected photon energy released from $\phi$ decays could disassociate the light elements ($\sim 10^{4}\textup{--}10^{6}~{}{\rm s}$) and reduce their abundances. The BBN constraint is sensitive to the mixing parameter in the range of $10^{-10}-10^{-14}$ for $m_{\phi}\gtrsim 10~{}{\rm MeV}$ Fradette et al. (2014). If the mediator has a longer lifetime and decays after cosmic recombination ($\sim 10^{13}~{}{\rm s}$), it could leave an imprint on the CMB power spectrum Adams et al. (1998); Chen and Kamionkowski (2004); Zhang et al. (2007); Galli et al. (2009); Slatyer et al. (2009). The Planck experiment has excluded $\epsilon_{\gamma}\sim 10^{-16}\textup{--}10^{-17}$ for the mediator mass in the range of $1\textup{--}300~{}{\rm MeV}$ Fradette et al. (2014). Proposed CMB missions, such as PIXIE Kogut et al. (2011), could be sensitive to $\epsilon_{\gamma}\sim 10^{-13}\textup{--}10^{-17}$ for a similar mediator mass range Berger et al. (2016). We expect these constraints will be stronger if the mediator abundance is also produced thermally, in addition to the freeze-in process.

VII.2 Direct detection

Direct detection experiments search for nuclear recoils due to particle collisions with DM in the local halo. Such interactions necessarily arise if mediator states interact with SM fermions as well as DM, as illustrated in the right panel of Fig. 21. For typical weakly interacting massive particles (WIMPs), direct detection cross sections are suppressed by the mass scale of heavy mediators (such as the Higgs boson). However, since the SIDM framework requires mediators to be light, the cross section could be strongly enhanced. Therefore, current null results place strong constraints on the mixing parameter connecting the two sectors Fornengo et al. (2011); Kaplinghat et al. (2014b). In addition, since DM self-scattering erases the non-thermalized features seen in CDM simulations, a Maxwell-Boltzmann DM velocity distribution, often used in calculating the direct detection rate, is more justified in SIDM than CDM Vogelsberger and Zavala (2013).

The low mediator mass scale makes an interesting prediction for direct detection. In contrast to WIMPs, for which scattering is a contact interaction, the mediator mass for SIDM is comparable to the typical momentum transfer $q$ in DM-nucleus scattering. Thus, SIDM may interact with nuclei through a long-range force, which yields an event spectrum that is more peaked toward lower recoil energies Fornengo et al. (2011); Li et al. (2015); Del Nobile et al. (2015). This feature may help distinguish between SIDM and WIMP recoil spectra in experiments Del Nobile et al. (2015). If the SIDM candidate is a bound state, the form factor of the bound state and the possibility of breakup of the bound state could produce new signatures in the recoil spectrum Laha and Braaten (2014); Laha (2015).

Del Nobile, Kaplinghat & Yu Del Nobile et al. (2015) recast the LUX Akerib et al. (2014) and SuperCDMS Agnese et al. (2014) constraints on the WIMP-nucleon cross section to set a limit on SIDM parameter space, as shown in Fig. 22. The model considered here consists of DM particle $\chi$ with a vector mediator $\phi$ that couples to the SM by kinetic mixing (dark photon). The blue band denotes the parameter space where $\sigma/m=0.1-10\;{\rm cm^{2}/g}$, as relevant for solving small scale anomalies, for $\alpha_{\chi}=0.01$. At present, LUX excludes this entire parameter space for $\varepsilon_{\gamma}\gtrsim 10^{-9}$ and $m_{\chi}\gtrsim 7~{}{\rm GeV}$. The SuperCDMS limit is about an order of magnitude weaker, but it is more sensitive to low mass SIDM models because of its lower energy threshold and lighter target nucleus (germanium compared to xenon).

In the event of a positive detection, detailed studies of the signal spectrum and annual modulation can distinguish between SIDM and WIMP models. Fig. 23 (left) shows the signal spectrum for SIDM in LUX, for fixed DM mass, taking into account realistic efficiency, acceptance, and energy resolution of the detector. The three spectra indicate models with different mediator masses and have been normalized to give the same rate within the signal range (area enclosed by the vertical dashed lines). SIDM models with lower mediator masses have spectra peaked toward lower nuclear recoil energy, due to the long-range nature of the interaction, compared to a contact interaction, e.g., in WIMP models.

However, in some cases, there is a degeneracy in the signal spectrum between mass and the nature of the interaction. For example, a $20~{}{\rm GeV}$ WIMP with contact interactions can yield the same spectrum at LUX compared to SIDM with $m_{\chi}=100~{}{\rm GeV}$ and $m_{\phi}=15~{}{\rm MeV}$. One way to break the degeneracy is to compare signals between experiments employing different target nuclei, such as SuperCDMS. In addition, these two models have very different annual modulation spectra, as shown in Fig. 23 (right). Remarkably, they have opposite phase: light WIMPs are peaked when the Earth moves into the galactic DM wind, since the number collisions above threshold is larger, while the SIDM model is peaked moving with the wind, since the cross section is enhanced for smaller velocities.

Direct detection experiments have already placed stringent constraints on possible couplings between SIDM and SM particles. For the minimal models considered here, upcoming searches have the opportunity to cover SIDM parameter space down to the BBN floor around $\varepsilon_{\gamma}\sim 10^{-10}$ Kaplinghat et al. (2014b). However, these limits may be avoided if SIDM is light ($m_{\chi}\lesssim{\rm GeV}$) or if the couplings between SM particles and mediators are different than those considered here.

Inelastic self-interactions provide another possibility for avoiding stringent direct detection constraints on the couplings between DM and the visible sector. Zhang (2017); Blennow et al. (2017). These limits may be evaded completely if scattering from nuclei requires an inelastic transition that is kinematically forbidden. For example, for $m_{\chi}=10~{}{\rm GeV}$, endothermic DM-nucleon scattering is kinematically forbidden in both xenon- and germanium-based detectors when $\delta\gtrsim 30~{}{\rm keV}$ Blennow et al. (2017). The required minimal DM velocity is $v_{\rm min}\approx\sqrt{2\delta/\mu_{\chi A}}$, where $\mu_{\chi A}$ is the DM-nucleus reduced mass, which is larger than the escape velocity of the MW halo, $\sim 750~{}{\rm km/s}$. For the same reason, up-scattering $\chi\chi\rightarrow\chi^{*}\chi^{*}$ is suppressed or forbidden in galactic halos. However, elastic self-scattering $\chi\chi\rightarrow\chi\chi$, with $\chi^{*}$ as an intermediate state, is not suppressed due to nonperturbative effects Schutz and Slatyer (2015).

VII.3 Indirect detection

A compelling feature of WIMP models is that the same annihilation processes that set the relic abundance can also give rise to indirect detection signals in DM halos at present. The signal strength is determined by the requirement that the WIMP relic density should be consistent with the observed DM abundance. For SIDM models with a light mediator $\phi$, the expectation is similar. First, DM will annihilate into $\phi\phi$ pairs, shown in the middle panel of Fig. 21, each of which in turn will decay to produce a detectable signal, e.g., $\phi\to e^{+}e^{-}$. The signal strength is largely determined by dark sector couplings only and do not turn off in the limit that the portal couplings to the SM become tiny (but nonzero). In addition, there are a number of characteristic features of SIDM indirect detection, which we now discuss.

Light mediators can boost the annihilation cross section in low-velocity environments due to the Sommerfeld enhancement Sommerfeld (1931). Attractive self-interactions distort the plane wave of the incoming DM particles and increase the probability of annihilation, illustrated diagramatically in Fig. 24 (left). The DM annihilation cross section can be factorized as

[TABLE]

where $(\sigma_{\rm ann}v_{\rm rel})_{0}$ is the velocity-weighted Born annihilation cross section and $S(v_{\rm rel})$ represents the nonperturbative Sommerfeld enhancement factor. It can be computed by numerically solving the non-relativistic Schrödinger equation (see, e.g., Refs. Hisano et al. (2005); Arkani-Hamed et al. (2009)) or analytically by using the Hulthèn potential as an approximation for the true Yukawa potential Cassel (2010); Slatyer (2010); Feng et al. (2010b); Blum et al. (2016).

This effect has important implications for indirect detection Baer et al. (1999); Hisano et al. (2003, 2005, 2006); Cirelli et al. (2007); March-Russell et al. (2008); Cirelli et al. (2009); Arkani-Hamed et al. (2009); Lattanzi and Silk (2009); Feng et al. (2010b). During freeze-out, the DM velocity is $v_{\rm rel}\sim 0.3\,c$, while at later times—in DM halos today and during the recombination epoch—the velocity is far smaller. For light mediator models, this implies that the annihilation signal can be enhanced compared to the value needed to fix the relic density.272727To accurately calculate the relic density in DM models with a light mediator, one needs to consider several subtle effects Feng et al. (2009); Dent et al. (2010); Zavala et al. (2010); Feng et al. (2010b); van den Aarssen et al. (2012b); von Harling and Petraki (2014), such as enhanced annihilation after freeze-out, chemical recoupling, and bound state formation. If the two sectors were not in thermal equilibrium during DM freeze-out, the relic abundance also depends on the temperature ratio of the two sectors, see, e.g., Feng et al. (2008). Fig. 24 (right) shows that the thermally-averaged Sommerfeld enhancement can provide a large boost to annihilation on dwarf (10 km/s), MW (200 km/s), and cluster (1000 km/s) scales compared to the Born cross section. The effect becomes important for $\alpha_{\chi}m_{\chi}/m_{\phi}\gtrsim 1$. The peaks are resonant enhancements due to quasi-bound state formation, which is the same mechanism that leads to the resonant DM self-scattering discussed in §VI.3. When the mass ratio $m_{\chi}/m_{\phi}$ is larger than $\sim 10^{3}\textup{--}10^{5}$, the mediator is effectively massless and the enhancement factor is simply Coulomb-like, given by $S\approx\pi\alpha_{\chi}/v_{\rm rel}$. On the other hand, the case with $m_{\chi}/m_{\phi}\sim 1$ corresponds to the usual WIMP and the enhancement is typically negligible.

Models along these lines were originally proposed by Arkani-Hamed et al. Arkani-Hamed et al. (2009) and Pospelov & Ritz Pospelov and Ritz (2009) to explain the increasing positron fraction in high-energy cosmic rays reported by the PAMELA collaboration Adriani et al. (2009), later confirmed by the Fermi-LAT Ackermann et al. (2012) and AMS-02 Aguilar et al. (2013) collaborations. The signal can be fitted if the DM mass is $\sim{\rm TeV}$ and the annihilation cross section is a factor of ${\cal O}(100\textup{--}1000)$ larger than the standard WIMP cross section due to the Sommerfeld enhancement Cirelli and Strumia (2008). However, this explanation remains in tension with other observations. For instance, energy injection from DM annihilation may increase the ionization fraction during recombination, resulting in a distorted CMB spectrum Chen and Kamionkowski (2004); Kamionkowski and Profumo (2008); Galli et al. (2009); Slatyer et al. (2009); Hisano et al. (2011); Hutsi et al. (2011); Madhavacheril et al. (2014); Lopez-Honorez et al. (2013). Since the DM velocity is very low during the recombination epoch, $v_{\rm rel}\sim 10^{-8}c$, the predicted Sommerfeld effect is very large at this epoch Feng et al. (2010b); Ade et al. (2016).

More recently, Fermi-LAT observations have provided a purported excess in GeV $\gamma$-rays coming from the Galactic Center (GC) with an energy spectrum and spatial morphology consistent with WIMP annihilation Goodenough and Hooper (2009); Hooper and Goodenough (2011); Abazajian and Kaplinghat (2012); Ajello et al. (2016). Kaplinghat, Linden & Yu Kaplinghat et al. (2015) proposed that SIDM annihilation, shown in Fig. 21 (center), may explain this excess as well. In this model, DM annihilates into mediator particles, which then decay into energetic $e^{+}e^{-}$ pairs, which in turn yield GeV $\gamma$-rays by Compton upscattering starlight in the GC. In contrast to WIMPs, which produce $\gamma$-rays more directly in decay cascades, $\gamma$-ray signals from SIDM are negligible in dwarf galaxies due to their small stellar density. On the other hand, the clustering of GC photons suggests that unresolved point sources are a more likely explanation of this excess Bartels et al. (2016); Lee et al. (2016).

Recently, Bringmann et al. Bringmann et al. (2017) and Cirelli et al. Cirelli et al. (2016) have found that SIDM annihilation is excluded over a wide range of parameters, $100\;{\rm MeV}<m_{\chi}<100\;{\rm TeV}$ and $1\;{\rm MeV}<m_{\phi}<1\;{\rm GeV}$, based on indirect detection bounds from Fermi, AMS-02, and Planck. However, there are number of caveats to keep in mind:

• •

The bound can be relaxed if the dark and visible sectors have different temperatures during the freeze-out. In fact, given present constraints from direct detection, the couplings between the two sectors are likely too small to achieve thermalization Kaplinghat et al. (2014b); Del Nobile et al. (2015). To fit the GC excess, the dark-to-visible temperature ratio is $\sim 0.1$ Kaplinghat et al. (2015).

• •

It is assumed that $\phi$ is a vector mediator that decays via kinetic mixing. If $\phi$ decays via $Z$ mixing, then a sizable fraction of decays produce neutrinos, which are far less constrained. Alternatively, if $\phi$ is a scalar, annihilation is a $p$-wave process and thus the cross section is velocity-suppressed. However, CMB bounds may still be relevant for heavy SIDM due to bound state formation An et al. (2016a).

• •

Indirect detection constraints are avoided if SIDM is asymmetric. (Late-time particle-antiparticle oscillations can in principle repopulate the symmetric density Cohen and Zurek (2010); Buckley and Profumo (2012); Cirelli et al. (2012); Tulin et al. (2012).)

Inelastic self-interactions can give rise to a different class of indirect detection signals. If DM $\chi$ up-scatters to an excited state $\chi^{*}$, it may subsequently de-excite through a decay, e.g., $\chi^{*}\to\chi\gamma$, producing $\gamma$-ray signatures. Finkbeiner & Weiner Finkbeiner and Weiner (2016) suggested that the recently claimed 3.5 keV line signal in M31 and several clusters Bulbul et al. (2014); Boyarsky et al. (2014) could arise from this effect. Since up-scattering is restricted by the kinematics of the mass gap, this scenario explains why the signal is seen in larger systems but not in dwarf galaxies. Boddy et al. Boddy et al. (2014b) showed that composite SIDM can produce the same effect, where the mass gap arises as a hyperfine splitting that is naturally in the right energy range.

We may also expect indirect detection signals from the centers of large astrophysical bodies such as the Sun and the Earth, which may capture ambient DM particles Gould (1987). Annihilation of those captured DM particles may lead to some striking signals. For example, WIMP annihilation in the Sun could produce high-energy neutrinos detectable at the IceCube Aartsen et al. (2016) and SuperK Choi et al. (2015) experiments, and annihilation in the Earth may produce anomalous heat Mack et al. (2007). Compared to WIMPs, DM models with a light mediator have several interesting features in generating such signals. First, mediators may be so long-lived that their decays occur outside of the astrophysical object, allowing for the detection of charged final states Batell et al. (2010); Schuster et al. (2010). Second, the capture rate can be enhanced in the Sun due to the self-scattering Zentner (2009); Albuquerque et al. (2014); Chen et al. (2014, 2016). Third, the Sommerfeld effect can increase the annihilation cross section and shorten the equilibrium timescale of DM particles captured in the Earth Delaunay et al. (2009); Feng et al. (2016a), resulting in a maximal signal flux. Fourth, the capture cross section is momentum dependent Feng et al. (2016a, b); Kouvaris et al. (2016); Smolinsky and Tanedo (2017), in contrast a contact interaction for WIMPs.

Lastly, even if annihilation signals are absent (e.g., in the asymmetric case), DM accumulation in astrophysical objects may still cause observable effects. For example, captured asymmetric SIDM in stars can affect their evolution Frandsen and Sarkar (2010); Zentner and Hearin (2011); Bramante et al. (2014); Kouvaris (2012); Eby et al. (2016). Interestingly, a Rutherford-like DM-baryon interaction can significantly enhance the energy transport and produce an impact on the helioseisomology data, improving the agreement between the best solar model and the helioseismic data Lopes et al. (2014).

VII.4 Collider searches

DM particles produced in colliders, although not detected directly, are usually inferred through a net momentum imbalance when they escape the detector. The typical WIMP signal consists of initial-state radiation recoiling against missing energy Birkedal et al. (2004); Feng et al. (2006); Beltran et al. (2010); Goodman et al. (2011); Bai et al. (2010); Goodman et al. (2010); Bai and Tait (2013); Zhou et al. (2013); Carpenter et al. (2014); Hochberg et al. (2016a), which is being actively searched for at the Large Hadron Collider (LHC) Chatrchyan et al. (2012); Aad et al. (2013). However, if the dark and visible sectors are coupled through sub-GeV mediators only, the LHC is not the best environment to test these models (compared to direct detection and other types of searches discussed) since the couplings must be small and the light mediator cannot produce DM through an on-shell decay.

On the other hand, production can be enhanced at high energy colliders if there exist other types of interactions distinct from what mediates self-interactions. For example, if the light mediator $\phi$ couples to the SM through Higgs or $Z$ mixing, then on-shell decays of the Higgs or $Z$ boson can produce DM particles if kinematically allowed, contributing to their invisible widths Kouvaris et al. (2015). Alternatively, Tsai et al. Tsai et al. (2016) proposed that SIDM may couple to quarks through non-renormalizable contact operators. In particular, for operators that are poorly constrained by direct detection, these interactions may be large enough to allow sizable production of SIDM at the LHC, yielding mono-jet signatures similar to WIMP searches.

Additionally, SIDM may have unique signals at the LHC due to the light mediator $\phi$. First, since $m_{\phi}$ is much smaller than the typical parton energy, highly boosted DM particles will emit collinear $\phi$ radiation, analogous to collinear photon emission in QED Buschmann et al. (2015). Since these $\phi$ particles will decay to produce leptons, the signature of this effect is pairs of collimated lepton-jets. If $\phi$ is a dark photon, this type of search is complementary to other dark photon probes (discussed below), albeit with an additional model dependence for production. Second, if two SIDM particles are produced near the threshold, they may form a bound state due to strong self-interactions mediated by $\phi$, provided the interaction range is larger than the Bohr radius Shepherd et al. (2009); An et al. (2016b); Tsai et al. (2016); Bi et al. (2017). The lifetime of the bound state is short and it subsequently decays back to mediators or SM particles inside the detector. This process does not have a large SM background and provides a direct way for testing the self-interacting nature of DM at particle colliders An et al. (2016b); Tsai et al. (2016).

Recently, there has been considerable progress in searches for light weakly-coupled dark sectors at high-luminosity experiments, such as $e^{+}e^{-}$ colliders and fixed target experiments Batell et al. (2009b); Essig et al. (2009); Bjorken et al. (2009); Batell et al. (2009c); Essig et al. (2011) (see, e.g., Ref. Alexander et al. (2016) and references therein). These studies may be able to discover the mediator $\phi$ for self-interactions purely through its coupling to the SM. For the case that $\phi$ is a dark photon (denoted $A^{\prime}$), Fig. 25 shows current and prospective constraints on the parameter space of kinetic mixing parameter $\varepsilon$ and mass $m_{A^{\prime}}$. In these studies, it is assumed that $A^{\prime}$ will decay into SM fermions, resulting in distinctive dilepton resonances. Alternatively, if DM is sufficiently light, low-mass mediators may decay predominantly into DM particles. This scenario may be detected as missing energy signatures at fixed target experiments Izaguirre et al. (2015a, b).

An et al. An et al. (2016b) showed that SIDM bound states may form at low-energy $e^{+}e^{-}$ colliders if the DM mass is sufficiently light and explored constraints from B-factories. Fig. 26 illustrates the production and decay processes. When bound states form, they subsequently decay to charged leptons, resulting in six-lepton final states with low SM backgrounds. Compared to direct $A^{\prime}$ searches at $e^{+}e^{-}$ colliders, these bound state signatures improve the sensitivity to SIDM by $\mathcal{O}(10)$ for very large DM self-couplings ($\alpha_{\chi}=0.5$) An et al. (2016b).

VII.5 Cosmological probes

If the dark sector is completely decoupled from the SM, all of the complementary signals of SIDM discussed above do not apply. However, in this case, many SIDM models predict a suppression of cosmological structure on small scales, which can be tested observationally.

Most SIDM studies follow the assumption that the linear matter power spectrum remains unchanged compared to CDM. However, several SIDM models—especially those that are decoupled from the SM—require additional model ingredients that can change this picture. For example, in light mediator models, the mediator $\phi$ must decay to avoid overclosure. Hence, there must exist new light states for $\phi$ to decay into if it is decoupled from the SM. Similarly, for atomic SIDM models, massless mediators responsible for atomic binding have a thermal density in the early Universe. Tight coupling between DM and a dark radiation component leads to acoustic oscillations and diffusion damping for DM, analogous to the baryon-photon fluid, modifying the linear matter power spectrum Feng et al. (2009); Cyr-Racine and Sigurdson (2013); van den Aarssen et al. (2012a); Cyr-Racine et al. (2014); Buckley et al. (2014); Boddy et al. (2016). Massless particles may also contribute to the effective degrees of freedom in the early Universe, potentially detectable in CMB observations (see, e.g., Refs. Baumann et al. (2016); Chacko et al. (2015, 2016); Abazajian et al. (2016); Ko and Tang (2017, 2016); Banerjee et al. (2016); Tang (2017); Brust et al. (2017)). Alternatively, in cannibalization models, DM particles cool more slowly during freeze-out if entropy is conserved in the dark sector, which can also suppress small scale power Carlson et al. (1992); de Laix et al. (1995).

Damping caused by interactions between DM and radiation can lead to observable cosmological signals, as first noticed by Boehm et al. Boehm et al. (2001, 2002); Boehm and Schaeffer (2005).282828These works considered interactions between DM and SM radiation, i.e., photons or neutrinos. Subsequent works have shown that both the missing satellites and TBTF problems can be alleviated due to damping if the DM-radiation scattering cross section is $2\times 10^{-9}m_{\chi}/{\rm GeV}$ times the Thomson cross section Boehm et al. (2014); Schewtschenko et al. (2016). Furthermore, van der Aarssen, Bringmann & Pfrommer van den Aarssen et al. (2012a) suggested that a light mediator coupled to both DM and neutrinos will generate DM self-interactions and DM-neutrino interactions to solve these issues plus the core-cusp problem. However, couplings between SM neutrinos and the dark sector are subject to various experimental constraints Ahlgren et al. (2013); Laha et al. (2014); Bertoni et al. (2015). The damping scale is set by the horizon size when kinetic decoupling occurs, which is of order $10\;{\rm pc}\times(T_{\rm kd}/{\rm 10\;MeV})$ for kinetic decoupling temperature $T_{\rm kd}$ Chen et al. (2001); Hofmann et al. (2001); Loeb and Zaldarriaga (2005); Profumo et al. (2006); Bertschinger (2006); Bringmann (2009). To estimate $T_{\rm kd}$, one equates the Hubble rate to the momentum transfer rate for scattering

[TABLE]

where $n_{r}\sim T^{3}$ is the radiation number density, $\left<\sigma_{\chi r}v\right>$ is the velocity-weighted scattering cross section between DM and radiation, and $T$ is the temperature of the thermal bath. Parametrizing the cross section to be $\left<\sigma_{\chi r}v\right>=T^{2}/m_{\phi}^{4}$, where $m_{\phi}$ represents the mediator mass scale for DM-radiation interactions, yields292929If the two sectors have different temperatures, $T_{\rm kd}$ also depends on the temperature ratio, see, e.g. Feng et al. (2009); Cyr-Racine et al. (2016).

[TABLE]

For WIMPs, the mediator $\phi$ represents weak-scale degrees of freedom that couple DM $\chi$ to SM radiation, with all masses $m_{\chi},m_{\phi}$ set by the weak-scale. Thus, the damping scale is $\mathcal{O}(10\;{\rm pc})$, much smaller than current observational limits.

On the other hand, SIDM may couple to a dark radiation component, such as light dark photons or sterile neutrinos. In this case, the relevant mediator mass scale is likely $m_{\phi}\sim 10\;{\rm MeV}$, which significantly delays DM kinetic decoupling from the thermal bath Hooper et al. (2007); van den Aarssen et al. (2012a); Dasgupta and Kopp (2014); Bringmann et al. (2014); Ko and Tang (2014b); Chu and Dasgupta (2014); Cherry et al. (2014); Bertoni et al. (2015); Binder et al. (2016); Bringmann et al. (2016). In this case, even for a weak-scale DM mass, one has $T_{\rm kd}\sim 1~{}{\rm keV}$ and damping scales of $\mathcal{O}(\sim 0.1\;{\rm Mpc})$. Through this enhanced damping effect, the missing satellites problem may be addressed within the SIDM framework van den Aarssen et al. (2012a).

The ETHOS collaboration performed cosmological N-body simulations of SIDM with damping due to dark radiation Vogelsberger et al. (2016). Their study considers four SIDM models (ETHOS-1–4) with different self-scattering cross sections and damping scales. Their simulations assume DM couples to a massless hidden fermion via the same mediator $\phi$ that generates self-interactions, however they also explore a more general parametrization of damping within a broader class of models Cyr-Racine et al. (2016). Fig. 27 (left) shows the linear DM power spectra for the ETHOS models, compared to those for warm DM and CDM. For each ETHOS model, the power spectrum and the self-scattering cross section are calculated in a self-consistent way for a given set of model parameters. As expected, the power spectra of the ETHOS models exhibit acoustic oscillations, which are different from free-streaming damping for warm DM. Fig. 27 (right) shows that the number of subhalos in a MW-sized halo is significantly reduced in ETHOS models. These models agree with observational results better than CDM, thus potentially alleviating the missing satellites problem.

It is of great interest to reassess all small scale issues in SIDM scenarios that include dark acoustic damping. From the model building perspective, dark radiation (or other modifications to the power spectrum) is a generic feature of hidden sectors that are completely decoupled from the visible Universe. Even in this “nightmare” scenario, SIDM may still be testable. In fact, the ETHOS studies have made clear that all small structure issues are highly interconnected: the combined effect of DM self-scattering and acoustic damping may provide too large of an effect on dwarf scales to solve all issues simultaneously. Future work remains to find the optimal model to solve all issues.

VIII Conclusions

VIII.1 Summary

Studies of astrophysical structure, combining both observations and N-body simulations, are confronting the paradigm of cold, collisionless DM. While this paradigm works exceedingly well to explain data on large scales, the situation is less clear on smaller scales. The inner densities of dwarf and LSB galaxies, as well as galaxy clusters, exhibit a mass deficit in their inner halos compared to expectations from CDM-only simulations. While this may signal the importance of baryon dynamics, or other limitations in interpreting these observations, the particle dynamics of self-interacting DM can provide a viable explanation. The figure of merit for self-interactions is $\sigma/m$, the cross section per unit DM mass. Provided that $\sigma/m$ is nonzero, but not too large, heat is transferred inward to the center of the halo, inducing a mass deficit as DM particles are energized and move outward to larger orbits.

The past few years have seen a renaissance for SIDM. On the simulation side, there have been several important developments:

• •

Halo shape observations of clusters Miralda-Escude (2002) and elliptical galaxies Feng et al. (2010a, 2009); Lin et al. (2012) had previously provided stringent upper limits on $\sigma/m$, at the level of $\sim 0.01\;{\rm cm^{2}/g}$. In light of higher resolution simulations with greater halo statistics, these limits have been reassessed and weakened to $\sim 1\;{\rm cm^{2}/g}$ Peter et al. (2013).

• •

Simulations have begun to explore the importance of baryons in SIDM halos. However, when it comes to stellar feedback on dwarf scales, there is little consensus. For a bursty star formation history, CDM and SIDM halos both have large density cores and are virtually indistinguishable Fry et al. (2015), while for a smoother history, CDM and SIDM halos are very different since baryonic effects are insufficient to generate cores Vogelsberger et al. (2014).

• •

Aside from feedback, the static gravitational potential of baryons can greatly impact halo profiles compared to SIDM-only expectations, especially in regions that are baryon-dominated Kaplinghat et al. (2014a). This effect, recently verified through simulations, can result in SIDM profiles as cuspy as their CDM counterparts and leads to a predictive correlation between baryonic and SIDM density profiles Elbert et al. (2016); Creasey et al. (2017). The scatter in rotation curves for a fixed halo mass is increased compared to CDM halos, lessening the diversity problem Kamada et al. (2016); Creasey et al. (2017).

• •

Recent cluster merger simulations have downgraded the magnitude of DM-stellar offsets induced by self-interactions during core passage Robertson et al. (2017a); Kim et al. (2016). Without $\mathcal{O}(10\;{\rm kpc})$ resolution on the DM and galactic centroids of the subclusters, it may not be feasible to improve sensitivities much beyond simple scattering depth arguments. On the other hand, post-merger misalignment of brightest central galaxies may provide greater sensitivities Kim et al. (2016).

• •

Simulations have moved beyond the minimal SIDM model with a constant cross section, albeit for only a few benchmark models. The scenarios considered include both velocity-dependent Vogelsberger et al. (2012) and angular-dependent Robertson et al. (2017b) cross sections, as well as models that include a dark radiation component as an additional source of damping on small scales Buckley et al. (2014); Cyr-Racine et al. (2016); Vogelsberger et al. (2016). However, the model-space of potential dependencies remains large Tulin et al. (2013a).

In parallel, semi-analytic methods based on the Jeans equation have been developed for relaxed halos, providing a complementary path for confronting SIDM against observations that does not rely on computationally expensive simulations Kaplinghat et al. (2014a). These methods have been tested against and agree with results from SIDM simulations, both without and with baryons Kaplinghat et al. (2016); Elbert et al. (2016); Creasey et al. (2017).

As summarized in Table 1, dwarf and LSB galaxies prefer a cross section of at least $\sim 1\;{\rm cm^{2}/g}$ to generate a density core consistent with observations. The strongest constraint on SIDM comes from cluster scales. For fixed $\sigma/m$, the self-scattering rate in clusters is boosted compared to dwarf galaxies by a larger collisional velocity, providing an effect out to larger physical radii. The absence of $\mathcal{O}(100\;{\rm kpc})$ DM cores in massive clusters disfavors a constant cross section of $1\;{\rm cm^{2}/g}$. Interestingly, a recent study of seven massive relaxed clusters, including data from lensing and stellar kinematics, supports the presence of $\mathcal{O}(10\;{\rm kpc})$ cores Newman et al. (2013b, a). These observations are consistent with $0.1\;{\rm cm^{2}/g}$, favoring velocity-dependent self-interactions as a unified explanation across all scales Kaplinghat et al. (2016); Elbert et al. (2016). Earlier cluster studies Miralda-Escude (2002); Meneghetti et al. (2001), which constrained $\sigma/m<0.1\;{\rm cm/g}$ based on the inner halo ellipticity and mass density of cluster MS2137-23, were overestimated for reasons discussed in Sec. III and additionally were based on lensing only, which lacks sensitivity to small radii compared to stellar kinematics data Sand et al. (2002). Although the Bullet Cluster famously confirms the collisionless nature of DM, null results in this and other merging systems are fully consistent with $0.1\;{\rm cm^{2}/g}$ on cluster scales. However, Abell 3827 and Abell 520 are exceptions, exhibiting significant DM-stellar offsets. If these offsets are explained by SIDM, larger cross sections are required, in tension with cluster core data.

Self-interactions are a general feature of hidden sector models for DM. The particle physics of self-interactions is not particularly exotic. As proof of principle, the self-scattering cross section for nucleons has a similar order of magnitude to what is required for SIDM. In fact, many SIDM models are quite reminiscent of the SM, such as QED- and QCD-like theories for DM involving dark photons or other hidden gauge forces. Although the mass of DM itself can take a range of values depending on the model, the hidden sector is typically characterized by additional mass scales that govern the strength of self-interactions. These scales, which represent the mass of mediator particles or binding/confinement energies in composite DM models, must lie below the GeV-scale in order to achieve $\sigma/m\sim 1\;{\rm cm^{2}/g}$. In contrast, WIMP DM cannot generate a large enough self-interaction cross section since all mass scales are assumed to lie at or above the weak scale.

Where it comes to experimental signatures, the SIDM paradigm is no less rich than the WIMP paradigm. Mediators for self-interactions may couple to SM particles through renormalizable interactions, such as kinetic mixing Holdom (1986) or Higgs portal Patt and Wilczek (2006) operators. This has motivated, in part, experimental efforts to discover new sub-GeV states in high-luminosity colliders Alekhin et al. (2016); Alexander et al. (2016).

SIDM signals may also arise in traditional WIMP searches. For DM mass $m\gtrsim 10\;{\rm GeV}$, DM-nucleon scattering via light mediators is strongly constrained by direct detection searches, requiring very feeble couplings to the visible sector Fornengo et al. (2011); Kaplinghat et al. (2014b). Interestingly, the phase of the annual modulation signal can be reversed compared to WIMP DM, due to the momentum-dependence of the interaction, which provides an avenue for model discrimination in the event of a positive signal Del Nobile et al. (2015). Alternatively, lower mass SIDM can be detected via electron recoils Essig et al. (2012a, b, 2017).

For indirect detection, SIDM annihilation can be impacted by nonperturbative effects from the light mediator, including Sommerfeld enhancement and bound state formation An et al. (2016c). $s$-wave annihilation of SIDM into mediator states, which subsequently decay into SM states, is constrained by null results from AMS-02, Fermi-LAT, and Planck Bringmann et al. (2017). However, annihilation limits may be evaded in a variety of ways: e.g., if the dark sector is thermally decoupled from the visible sector (and colder) during freeze-out, a smaller annihilation cross section can produce the observed relic density. Self-interactions can also enhance the capture rate of DM in the sun, relevant for neutrino signals Zentner (2009), as well as induce collapse of neutron stars Kouvaris (2012); Güver et al. (2014).

On the other hand, self-interactions may be completely sequestered from the visible sector, in which case the observational probes discussed above do not apply. However, we argue that this is not a “nightmare” scenario for SIDM. Getting the correct relic density in SIDM requires additional model-building ingredients that may be tested cosmologically. Well-known scenarios involve either a dark radiation component (e.g., massless dark photons in atomic DM models Cyr-Racine and Sigurdson (2013)) or additional heating due to $3\to 2$ cannibalization (e.g., glueball DM models) Carlson et al. (1992); de Laix et al. (1995), both of which affect the matter power spectrum. For atomic DM, the effect of damping and dark acoustic oscillations at small scales may also address the dearth and diversity of MW satellites Buckley et al. (2014); Vogelsberger et al. (2016).

If SIDM makes up only a subleading fraction of the total DM density, additional possibilities may arise. For example, strongly self-interacting DM ($\sigma/m\gg 1$) can seed the initial growth of supermassive black holes Pollack et al. (2015) or undergo dissipative collapse to form compact dark disks Fan et al. (2013a).

VIII.2 Outlook

We conclude this report with a glimpse toward the horizon for astroparticle physics, small scale structure, and DM self-interactions. On the astrophysics side, this is already a data-rich field and will become even more so once planned instruments come online in the near future. On the particle physics side, a variety of experimental and observational techniques have the potential to access SIDM, including both traditional WIMP searches and more novel directions. We highlight several of the latter approaches below.

Future astronomical surveys: Over the next decade, upcoming instruments will boost the statistics and resolution for all observations relevant for SIDM. The Square Kilometer Array (SKA) radio telescope has the potential to observe hundreds of nearby spiral galaxies at resolutions below 100 pc, providing a deep look into their gas kinematics de Blok et al. (2015). This will provide a large and detailed sample of rotation curves, clarifying statistical arguments for the consistency of feedback-generated cores with $\Lambda$CDM cosmology Katz et al. (2016); Pace (2016), with an eye on potential outliers with large core sizes, $\gtrsim$ few kpc, that appear inconsistent with feedback (see Fig. 8). SKA will also provide new insights into the cooling and collapse of H gas—an important consideration for the conversion of gas into stars in numerical simulations de Blok et al. (2015). In addition, the Giant Magellan Telescope (GMT) 200 (2001) and the European Extremely Large Telescope (E-ELT) Hook (2005) can provide improved measurements of stellar kinematics of nearby dwarf galaxies.

Optical/near-infrared surveys with the Large Synoptic Survey Telescope (LSST) Abell et al. (2009), Euclid Amendola et al. (2013) and Wide Field Infrared Survey Telescope (WFIRST) Spergel et al. (2013) will observe thousands of strong lensing systems on group and cluster scales. Both space-based (Euclid, WFIRST) and ground-based (LSST) instruments play a complementary role Jain et al. (2015). For probing self-interactions, statistical tests of halo ellipticity, as advocated in Ref. Peter et al. (2013), will improve. In addition, the prospects for testing SIDM using minor mergers (substructure infalling into cluster halos) will be improved as well, with increases both in cluster statistics and in background source densities used for weak lensing studies Massey et al. (2011); Harvey et al. (2014). For major cluster mergers, mass maps for strong and weak lensing will become more accurate than before and moreover will benefit from a wider field of view compared to the current state-of-the-art observations from the Hubble Space Telescope Spergel et al. (2013).

Substructure: Certain classes of SIDM models (e.g., coupling to a radiation component van den Aarssen et al. (2012a); Bringmann et al. (2014); Buckley et al. (2014); Cyr-Racine et al. (2016); Boddy et al. (2016)) may damp the matter power spectrum on small scales. Lyman-$\alpha$ forest data exclude significant damping for halos above $3\times 10^{8}{\rm\,M}_{\odot}$ Viel et al. (2013). Smaller halos predicted by $\Lambda$CDM would be devoid of baryons and must be discovered gravitationally, either through lensing Mao and Schneider (1998); Chiba (2002); Dalal and Kochanek (2002) or perturbations to dynamically cold structures in the MW, such as tidal streams Johnston (1998); Ibata et al. (2002); Johnston et al. (2002); Siegal-Gaskins and Valluri (2008). Recently, studies of lensed dusty star-forming galaxies have shown great promise to probe substructure Hezaveh et al. (2013). Observations from the Atacama Large Millimeter/submillimeter Array (ALMA) can exploit individual star-forming clumps as sources to probe the lensing potential on small angular scales. The first studies based on ALMA science verification data have already found evidence for a $10^{9}{{\rm\,M}_{\odot}}$ subhalo in one elliptical galaxy, with a potential to reach down to $10^{7}{\rm\,M}_{\odot}$ in the future Hezaveh et al. (2016) (see also Inoue et al. (2016)). Another method for finding substructure uses tidal streams from globular clusters. These low-dispersion stellar systems are prone to disruption from subhalo fly-bys, leaving gaps in the stream. While observed tidal stream gaps are qualitatively consistent with the rate for subhalo encounters expected from $\Lambda$CDM Yoon et al. (2011); Carlberg et al. (2012), improved stream modeling, coupled with future observations Erkal and Belokurov (2015), offer the potential to probe down to $10^{5}{\rm\,M}_{\odot}$ subhalos Bovy et al. (2016).

Testing bursty feedback models: In principle, observations can test whether star formation histories in dwarf and LSB galaxies are bursty enough to transform cusps into cores. One method is to compare tracers of the star formation rate that have sensitivity to different time scales: e.g., H$\alpha$/H$\beta$ and far-ultraviolet (FUV) emission, which probe the mean rate over the past 10 Myr and 200 Myr, respectively Weisz et al. (2012) (see also Ref. Kauffmann (2014)). Nearby dwarf galaxies exhibit considerable scatter in their H$\alpha$/FUV ratios, which is interpretted as evidence for burstiness Weisz et al. (2012). Recent FIRE simulations Sparre et al. (2017) have been able to reproduce these observations for the most part (except for a non-negligible fraction of low-mass galaxies that have very small H$\alpha$/FUV ratios, which may indicate that the simulations are overestimating somewhat the peak-to-trough variation in the star formation rate). Additionally, if outflows from starbursts drive rapid fluctuations in the gravitational potential, the stellar distribution is heated gravitationally alongside the DM El-Badry et al. (2016). This yields a testable correlation between star formation rate and stellar velocity dispersion El-Badry et al. (2017).

LSBs have quiescent star formation today de Blok and McGaugh (1997), although bursty feedback may have occurred in the past to form a core. It is important, therefore, to understand the redshift dependence of these processes. Studies of H$\beta$/FUV ratios in low-mass galaxies at $0.4<z<1$ find evidence that star formation is burstier at higher redshifts Guo et al. (2016). However, applying these probes to faint DM-dominated galaxies at distant redshifts is challenging. An alternative method is to study the radial distributions of stellar ages and chemical compositions in nearby galaxies. Here, feedback can imprint itself dramatically: since feedback causes dynamical heating of the stellar (as well as DM) density, stars undergo outward radial migration El-Badry et al. (2016). Since this effect builds up the most for the oldest (metal-poor) stars, it can potentially invert the age and metallicity gradients in galaxies compared to their intrinsic profiles. This process is expected to be most pronounced for DM halos with the largest cores if both arise through feedback. On the observational side, Young et al. Young et al. (2015) have initiated a program of MUltiwavelength observations of the Structure, Chemistry, and Evolution of LSB galaxies (MUSCEL) with promise to shed light on these issues.

Experimental detection of light dark matter: Searches for DM at colliders and in direct detection experiments have focused primarily on the weak scale, motivated by theoretical expectations. With the absence of any definitive signals thus far, theorists have turned toward new directions and the idea of light dark matter (LDM), in the keV–GeV mass range, has garnered much attention recently Alekhin et al. (2016); Alexander et al. (2016).303030If DM is a thermal relic, it cannot be below the keV scale based on Lyman-$\alpha$ constraints on warm dark matter (see, e.g., Ref. Viel et al. (2013)). Even if DM is a thermal relic of a hidden sector that is much colder than the visible sector, its mass cannot be below 1.5 keV due to free streaming and phase space constraints Das and Sigurdson (2012). Models along these lines typically require light mediators to be viable, and hence these scenarios are a natural framework for SIDM Lin et al. (2012). There are many proposals to detect LDM experimentally, offering the potential to study the physics of DM self-interactions in the laboratory. Accelerator-based proposals include missing-energy signatures at GeV-scale colliders Essig et al. (2013); Izaguirre et al. (2015a), as well as the production and detection of dark matter beams produced in fixed target and beam dump experiments Batell et al. (2009c); Izaguirre et al. (2013); Kahn et al. (2015); Alekhin et al. (2016); Battaglieri et al. (2016). The MiniBoone collaboration recently performed a dedicated search along these lines Aguilar-Arevalo et al. (2017). Alternatively, there are many ideas for the direct detection of LDM as well, such as WIMP search constraints repurposed in terms of electron recoils Essig et al. (2012a, b, 2017), as well as new detector technologies with ultra-low thresholds Graham et al. (2012); Hochberg et al. (2016b); Essig et al. (2016); Schutz and Zurek (2016); Derenzo et al. (2016).

Theoretical modeling and numerical simulations: On the theory side, it is important to understand the role of DM self-interactions in galaxy formation together with baryon physics. Some issues can be studied using a quasi-equilibrium analysis, as we discussed in §IV. One example is the tight connection predicted between DM and baryon distributions Kaplinghat et al. (2014a); Elbert et al. (2016); Creasey et al. (2017), which may explain the long-mysterious halo-disk conspiracy van Albada and Sancisi (1986); Kamada et al. (2016). On the other hand, there are dynamical questions of great interest as well that must be addressed using hydrodynamical SIDM simulations (see Ref. Di Cintio et al. (2017)). For example, SIDM inner halos are expected to be more stable than for CDM in response to violent baryonic outflows since self-interactions rapidly redistribute energy deposited from feedback. Additionally, the star formation rate may be more suppressed in dwarf galaxies with SIDM halos because they have shallower gravitational potential wells from which UV heating can remove gas more effectively. This suggests that the missing satellites problem and the TBTF problem in the field may be more easily solved in SIDM than CDM with baryon feedback.

Acknowledgments: We are indebted to James Bullock, Peter Creasey, Jonathan Feng, Ran Huo, Manoj Kaplinghat, Rachel Kuzio de Naray, David Morrissey, Laura Sales, Kristine Spekkens, Philip Tanedo, and Kathryn Zurek for helpful discussions and collaborations on this subject. This work is supported by the Natural Science and Engineering Research Council of Canada (ST), the U.S. Department of Energy under Grant No. de-sc0008541 (HBY), and the Hellman Fellows Fund (HBY).