- Split View
-
Views
-
Cite
Cite
Miguel Rocha, Annika H. G. Peter, James S. Bullock, Manoj Kaplinghat, Shea Garrison-Kimmel, Jose Oñorbe, Leonidas A. Moustakas, Cosmological simulations with self-interacting dark matter – I. Constant-density cores and substructure, Monthly Notices of the Royal Astronomical Society, Volume 430, Issue 1, 21 March 2013, Pages 81–104, https://doi.org/10.1093/mnras/sts514
- Share Icon Share
Abstract
We use cosmological simulations to study the effects of self-interacting dark matter (SIDM) on the density profiles and substructure counts of dark-matter haloes from the scales of spiral galaxies to galaxy clusters, focusing explicitly on models with cross-sections over dark-matter particle mass σ/m = 1 and 0.1 cm2 g−1. Our simulations rely on a new SIDM N-body algorithm that is derived self-consistently from the Boltzmann equation and that reproduces analytic expectations in controlled numerical experiments. We find that well-resolved SIDM haloes have constant-density cores, with significantly lower central densities than their cold dark matter (CDM) counterparts. In contrast, the subhalo content of SIDM haloes is only modestly reduced compared to CDM, with the suppression greatest for large hosts and small halo-centric distances. Moreover, the large-scale clustering and halo circular velocity functions in SIDM are effectively identical to CDM, meaning that all of the large-scale successes of CDM are equally well matched by SIDM. From our largest cross-section runs, we are able to extract scaling relations for core sizes and central densities over a range of halo sizes and find a strong correlation between the core radius of an SIDM halo and the NFW scale radius of its CDM counterpart. We construct a simple analytic model, based on CDM scaling relations, that captures all aspects of the scaling relations for SIDM haloes. Our results show that halo core densities in σ/m = 1 cm2 g−1 models are too low to match observations of galaxy clusters, low surface brightness spirals (LSBs) and dwarf spheroidal galaxies. However, SIDM with σ/m ≃ 0.1 cm2 g−1 appears capable of reproducing reported core sizes and central densities of dwarfs, LSBs and galaxy clusters without the need for velocity dependence. Higher resolution simulations over a wider range of masses will be required to confirm this expectation. We discuss constraints arising from the Bullet cluster observations, measurements of dark-matter density on small scales and subhalo survival requirements, and show that SIDM models with σ/m ≃ 0.1 cm2 g−1 ≃ 0.2 barn GeV−1 are consistent with all observational constraints.
1 INTRODUCTION
There is significant evidence that some form of dark matter dominates the gravitating mass in the Universe and its abundance is known to great precision (Komatsu et al. 2011). The most popular candidate for dark matter is the class of weakly interacting massive particles (WIMPs), of which supersymmetric neutralinos are examples (Steigman & Turner 1985; Griest 1988; Jungman, Kamionkowski & Griest 1996). WIMPs are stable, with negligible self-interactions, and are non-relativistic at decoupling (‘cold’). It is important to recognize that of these characteristics, it is primarily their coldness that is well tested via its association with significant small-scale power. Indeed, WIMPs are the canonical cold dark matter (CDM) candidate. Cosmological models based on CDM reproduce the spatial clustering of galaxies on large scales quite well (Reid et al. 2010) and even the clustering of galaxies on ∼1 Mpc scales appears to match that expected for CDM subhaloes (Kravtsov et al. 2004; Conroy, Wechsler & Kravtsov 2006; Trujillo-Gomez et al. 2011; Reddick et al. 2012).
Beyond the fact that the Universe appears to behave as expected for CDM on large scales, we have few constraints on the microphysical parameters of the dark matter, especially those that would manifest themselves at the high densities associated with cores of galaxy haloes. It is worth asking what (if anything) about vanilla CDM can change without violating observational bounds. In this paper, we use cosmological simulations to explore the observational consequences of a CDM particle that is strongly self-interacting, focusing specifically on the limiting case of velocity-independent, elastic scattering.
Dark-matter particles with appreciable self-interactions have been discussed in the literature for more than two decades (Carlson, Machacek & Hall 1992; Machacek, Carlson & Hall 1993; de Laix, Scherrer & Schaefer 1995; Firmani et al. 2000; Spergel & Steinhardt 2000), and are now recognized as generic consequences of hidden-sector extensions to the Standard Model (Pospelov, Ritz & Voloshin 2008; Ackerman et al. 2009; Arkani-Hamed et al. 2009; Feng et al. 2009; Feng, Kaplinghat & Yu 2010; Loeb & Weiner 2011). Importantly, even if dark sector particles have no couplings to Standard Model particles they might experience strong interactions with themselves, mediated by dark gauge bosons (see Feng 2010; Peter 2012 for reviews). The implication is that astrophysical constraints associated with the small-scale clustering of dark matter may be the only way to test these scenarios.
Phenomenologically, self-interacting dark matter (SIDM) is attractive because it offers a means to lower the central densities of galaxies without destroying the successes of CDM on large scales. Cosmological simulations that contain only CDM indicate that dark-matter haloes should be cuspy and with (high) concentrations that correlate with the collapse time of the halo (Navarro, Frenk & White 1997; Bullock et al. 2001; Wechsler et al. 2002). This is inconsistent with observations of galaxy rotation curves, which show that galaxies are less concentrated and less cuspy than predicted in CDM simulations (e.g. Flores & Primack 1994; Salucci & Burkert 2000; Gentile et al. 2004; Simon et al. 2005; Kuzio de Naray, McGaugh & de Blok 2008; de Blok 2010; Dutton et al. 2011; Kuzio de Naray & Spekkens 2011; Oh et al. 2011a; Walker & Peñarrubia 2011; Castignani et al. 2012; Salucci et al. 2012). Even for clusters of galaxies, the density profiles of the host dark-matter haloes appear in a number of cases to be shallower than predicted by CDM-only structure simulations, with the total (dark matter + baryons) density profile in a closer match to the CDM prediction for the dark matter alone (e.g. Sand et al. 2004, 2008; Newman et al. 2009, 2011; Coe et al. 2012; Umetsu et al. 2012).
One possible answer is feedback. In principle, the expulsion of gas from galaxies can result in lower dark-matter densities compared to dissipationless simulations, and thus bring CDM models in line with observations (Governato et al. 2010; Oh et al. 2011b; Pontzen & Governato 2012; Brook et al. 2012; Governato et al. 2012). However, a new level of concern exists for dwarf spheroidal (dSph) galaxies (Boylan-Kolchin, Bullock & Kaplinghat 2011, 2013; Ferrero et al. 2012). Systems with M* ∼ 106 M⊙ appear to be missing ∼ 5 × 107 M⊙ of dark matter compared to standard CDM expectations (Boylan-Kolchin et al. 2012). It is difficult to understand how feedback from such a tiny amount of star formation could have possibly blown out enough gas to reduce the densities of dSph galaxies to the level required to match observations (Boylan-Kolchin et al. 2012; Peñarrubia et al. 2012; Teyssier et al. 2012; Zolotov et al. 2012; Garrison-Kimmel et al., in preparation).
Numerical simulations have confirmed the expected phenomenology of core formation (Burkert 2000) though Kochanek & White (2000) emphasized the possibility that SIDM haloes could eventually become more dense than their CDM counterparts as a result of eventual heat flux from the inside out (much like core-collapse globular clusters). However, this process is suppressed when merging from hierarchical formation is included (for a discussion see Ahn & Shapiro 2005). We do not see clear signatures of core collapse in the haloes we analysed for σ/m = 1 cm2 g−1.
The first cosmological simulations aimed at understanding dwarf densities were performed by Davé et al. (2001) who used a small volume (4 h−1 Mpc on a side) in order to focus computational power on dwarf haloes. They concluded that σ/m = 0.1-10 cm2 g− 1 came close to reproducing core densities of small galaxies, favouring the upper end of that range but not being able to rule out the lower end due to resolution. Almost concurrently, Yoshida et al. (2000) ran cosmological simulations focusing on the cluster-mass regime. Based on the estimated core size of cluster CL 0024+1654, they concluded that cross-sections no larger than ∼0.1 cm2 g−1 were allowed, raising doubts that constant-cross-section SIDM models could be consistent with observations of both dwarf galaxies and clusters.
These concerns were echoed by Miralda-Escudé (2002) who suggested that SIDM haloes would be significantly more spherical than observed for galaxy clusters. Similarly, Gnedin & Ostriker (2001) argued that SIDM would lead to excessive subhalo evaporation in galaxy clusters. More recently, the merging cluster system known as the Bullet cluster has been used to derive the limits (68 per cent C.L.) σ/m < 0.7 cm2 g−1 (Randall et al. 2008) based on evaporation of dark matter from the subcluster and σ/m < 1.25 cm2 g−1 (Randall et al. 2008) based on the observed lack of offset between the Bullet subcluster mass peak and the galaxy light centroid. In order to relax this apparent tension between what was required to match dwarf densities and the observed properties of galaxy clusters, velocity-dependent cross-sections that diminish the effects of self-interaction in cluster environments have been considered (Firmani et al. 2000; Colín et al. 2002; Feng et al. 2009; Loeb & Weiner 2011; Vogelsberger, Zavala & Loeb 2012).
There are a few new developments that motivate us to revisit constant SIDM cross-sections of the order of σ/m ∼ 1 cm2 g−1. For example, the cluster (CL 0024+1654) used by Yoshida et al. (2000) to place one of the tightest limits at σ/m = 0.1 is now recognized as an ongoing merger along the line of sight (Czoske et al. 2001, 2002; Zhang et al. 2005; Jee et al. 2007; Jee 2010; Umetsu et al. 2010). This calls into question its usefulness as a comparison case for non-merging cluster simulations. In a companion paper (Peter et al. 2012), we use the same simulations described here to show that published constraints on SIDM based on halo shape comparisons are significantly weaker than previously believed. Further, the results presented below clearly demonstrate that the tendency for subhaloes to evaporate in SIDM models (Gnedin & Ostriker 2001) is not significant for σ/m ∼ 1 cm2 g−1. Finally (and related to the previous point), the best numerical analysis of the Bullet cluster (Randall et al. 2008) used initial cluster density profiles that were unmotivated cosmologically with central densities about a factor of 2 too high for the SIDM cross-sections considered (producing a scattering rate that is inconsistently high). Based on this observation, the Bullet cluster constraint based on evaporation of dark matter from the subcluster should be relaxed since the amount of subcluster mass that becomes unbound is directly proportional to the density of dark matter encountered in its orbit. Moreover, their model galaxies were placed in the cluster halo potentials without subhaloes surrounding them, an assumption (based on analytic estimates for SIDM subhalo evaporation) that is not supported by our simulations. This could affect the constraints based on the (lack of) offset between dynamical mass and light. Thus, we believe that the Bullet cluster constraints as discussed above are likely only relevant for models with σ/m > 1 cm2 g−1. However, the constraints could be made significantly stronger by comparing SIDM predictions to the densities inferred from the convergence maps since the central halo densities for σ/m ≃ 1 cm2 g−1 are significantly lower than the CDM predictions, as we show later.
Given these motivations, we perform a set of cosmological simulations with both CDM and SIDM. For SIDM we ran σ/m = 1 and 0.1 cm2 g−1 models (hereafter SIDM1 and SIDM0.1, respectively), i.e. models that we have argued pass the Bullet cluster tests. Our simulations provide us with a sample of haloes that span a mass range much larger than either Davé et al. (2001) or Yoshida et al. (2000) both with and without self-interactions.
The major conclusion we reach based on the simulations and the analytic model presented here is that an SIDM model with a cross-section over dark-matter particle mass ∼0.1 cm2 g−1 would be capable of reproducing the core sizes and central densities observed in dark-matter haloes at all scales, from clusters to dSphs, without the need for velocity dependence in the cross-section.
In the next section, we discuss our new algorithm to compute the self-interaction probability for N-body particles, derived self-consistently from the Boltzmann equation. We discuss this new algorithm in detail in Appendix A. In Section 2, we show how this algorithm is implemented in the publicly available code gadget2 (Springel 2005). We run tests that show that our algorithm gets the correct interaction rate and post-scattering kinematics. The results of these tests are given in Section 3. The cosmological simulations with this new algorithm are described in detail in Section 4. In Section 5.1 we provide some preliminary illustrations of our simulation snapshots, and in Section 5.2 we demonstrate that the large-scale statistical properties of SIDM are identical to CDM. In Section 5.4 we present the properties of individual SIDM1 and SIDM0.1 haloes and compare them to their CDM counterparts. In Section 5.4 we discuss the subhalo mass functions in our SIDM and CDM simulations and show that SIDM1 subhalo mass functions are very close to that of CDM in the range of halo masses we can resolve. We provide scaling relations for the SIDM1 halo properties in Section 6, and in Section 7 we present an analytic model that reproduces these scaling relations as well as the absolute densities and core radii of SIDM1 haloes. We use these scaling relations and the analytic model to make a broad-brush comparison to observed data in Section 8. We present a summary together with our final conclusions in Section 9.
2 SIMULATING DARK MATTER SELF-INTERACTIONS
Our simulations rely on a new algorithm for modelling SIDM with N-body simulations. Here we introduce our approach and provide a brief summary. In Appendix A, we derive the algorithm explicitly starting with the Boltzmann equation and give details for general implementation.
In N-body simulations, the simulated (macro)particles represent an ensemble of many dark-matter particles. Each simulation particle of mass mp can be thought of as a patch of dark-matter phase-space density. In our treatment of dark matter self-scattering, the phase-space patch of each particle is represented by a delta function in velocity and a spatially extended kernel W(r, hsi), smoothing out the phase space in configuration space on a self-interaction smoothing length hsi. The value of hsi needs to be set by considering the physical conditions of the problem (see Section 3) as it specifies the range over which N-body particles can affect each other via self-interactions. In principle, hsi could be different for each particle and vary depending on the local density, but in the simulations presented here we fix hsi to be the same for all particles in a given simulation, setting the size of hsi according to the lowest densities at which self-interactions are effective for a given cross-section.
Our method for simulating scattering differs from previous approaches in a few key ways. It is most similar to that of Davé et al. (2001) in that we both directly consider interactions between pairs of phase-space patches and rely on a scattering rate similar in form to equation (3). The difference is that their geometric factor gji is not the same – our factor arises explicitly from the overlap in patches of phase space between neighbouring macroparticles, as derived from the collision term in the Boltzmann equation (see Appendix A for details). Other authors determine the scattering rate Γ of individual phase-space patches based on estimates of the local mass density (typically using some number of nearest neighbours or using a smoothed particle hydrodynamics kernel). The Monte Carlo is then based on an estimated scattering rate of an individual particle on the background, and a scattering partner is only chosen after a scattering event is determined to have occurred (Kochanek et al. 2000; Yoshida et al. 2000; Colín et al. 2002; Randall et al. 2008). The scattering probability in this latter approach is not symmetric. For macroparticles of identical mass, P(i|j) = P(j|i) explicitly in our approach, but not the other approach because the density estimated at the position of macroparticle i need not be the same as that estimated at the position of particle j. In the future, there should be a direct comparison among these scattering algorithms to determine if they yield consistent results.
We have implemented our algorithm in the publicly available version of the cosmological simulation code gadget2 (Springel 2005). gadget2 computes the short-range gravitational interactions by means of a hierarchical multipole expansion, also known as a tree algorithm. Particles are grouped hierarchically by a repeated subdivision of space, so their gravitational contribution can be accounted by means of a single multipole force computation. A cubical root node encompasses the full mass distribution. The node is repeatedly subdivided into eight daughter nodes of half the side length each (an octree) until one ends up with ‘leaf’ nodes containing single particles. Forces for a given particle are then obtained by ‘walking’ the tree, opening nodes that are too close for their multipole expansion to be a correct approximation to their gravitational contribution. In gadget2, spurious strong close encounters by particles are avoided by convolving the single-point particle density distribution with a normalized spline kernel (‘gravitational softening’).
The time-step criterion is also modified to assure that the scattering probability for any pair of particles is small, P = Γ δt ≪ 1. An individual particle time step is decreased by a factor of 2 if during the last tree-walk the maximum probability of interaction for any pair involving such a particle was Pmax > 0.2. Once a particle time step is modified due to the previous restriction, if Pmax < 0.1 for such a particle and its current time step is smaller than the one given by the standard criterion on gadget2, we increase it by a factor of 2.
3 TEST OF THE SIDM IMPLEMENTATION
Before performing cosmological simulations, we carried out a controlled test of the implementation in order to make sure the scattering rate and kinematics are correctly followed in the code, and to determine the optimum value of the SIDM softening kernel length hsi. The simplest and cleanest scenario for testing our implementation consists of a uniform sphere of particles moving through a uniform field of stationary background particles. The coordinate system is defined such as the sphere is moving along the positive z-direction with constant velocity vs. The particles forming the sphere and those forming the background field are tagged as different types within the code and here we will refer to them simply as sphere (s) and background (bg) particles, respectively. We only allow scatterings involving two different types of particles (i.e. sphere–background interactions only) and turn off gravitational forces among all of the particles. Furthermore, all particles have the same mass mp.
We also check the kinematics of the scatters in this test simulation and describe the results in Appendix B. The resulting kinematics and number of interactions from our test simulation agree well with the expectations from the theory as long as hsi(ρbg/mp)1/3 ≳ 0.2.
4 OVERVIEW OF COSMOLOGICAL SIMULATIONS
We initialize our cosmological simulations using the Multi-Scale Initial Conditions (MUSIC) code of Hahn & Abel (2011). We have a total of four initial condition sets, each run with both CDM and SIDM. The first two are cubic volumes of 25 and 50 h−1 Mpc on a side, each with 5123 particles. As discussed below, these simulations allow us to resolve the structure of a statistical sample of group (∼ 1013 M⊙) and cluster (∼ 1014 M⊙) haloes.
The second two initial conditions concentrate computational power on zoom regions (Katz & White 1993) drawn from the 50 h−1 Mpc box, specifically aimed at exploring the density structure of two smaller haloes, one with virial mass2Mvir = 7.1 × 1011 h− 1 M⊙ = 1 × 1012 M⊙ (Z12) and one with Mvir = 3.5 × 1011 h− 1 M⊙ = 5 × 1011 M⊙ (Z11). The Z12 run in particular is fairly high resolution, with more than five million particles in the virial radius. Table 1 summarizes the simulation parameters. The cosmology used is based on7-year Wilkinson Microwave Anisotropy Probe results for a ΛCDM universe: h = 0.71, Ωm = 0.266, ΩΛ = 0.734, Ωb = 0.0449, ns = 0.963 and σ8 = 0.801 (Komatsu et al. 2011).
Name . | Volume LBox ( h−1 Mpc) . | Number of particles Np . | Particle mass mp ( h−1 M⊙) . | Force softening ϵ ( h−1 kpc) . | Smoothing length hsi ( h−1 kpc) . | Cross-section σ/m ( cm2 g−1) . |
---|---|---|---|---|---|---|
CDM-50 | 50 | 5123 | 6.88 × 107 | 1.0 | – | 0 |
CDM-25 | 25 | 5123 | 8.59 × 106 | 0.4 | – | 0 |
CDM-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | – | 0 |
CDM-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | – | 0 |
SIDM0.1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 0.1 |
SIDM0.1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 0.1 |
SIDM0.1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 0.1 |
SIDM0.1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 0.1 |
SIDM1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 1 |
SIDM1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 1 |
SIDM1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 1 |
SIDM1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 1 |
Name . | Volume LBox ( h−1 Mpc) . | Number of particles Np . | Particle mass mp ( h−1 M⊙) . | Force softening ϵ ( h−1 kpc) . | Smoothing length hsi ( h−1 kpc) . | Cross-section σ/m ( cm2 g−1) . |
---|---|---|---|---|---|---|
CDM-50 | 50 | 5123 | 6.88 × 107 | 1.0 | – | 0 |
CDM-25 | 25 | 5123 | 8.59 × 106 | 0.4 | – | 0 |
CDM-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | – | 0 |
CDM-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | – | 0 |
SIDM0.1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 0.1 |
SIDM0.1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 0.1 |
SIDM0.1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 0.1 |
SIDM0.1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 0.1 |
SIDM1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 1 |
SIDM1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 1 |
SIDM1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 1 |
SIDM1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 1 |
aThe Z11 and Z12 runs are zoom simulations with multiple particle species concentrating on haloes of mass Mvir = 5 × 1011 and 1.0 × 1012 M⊙, respectively (no h). The volumes listed refer to the number of virial radii used to find the Lagrangian volumes associated with the zoom. The particle properties listed are for the highest resolution particles only.
Name . | Volume LBox ( h−1 Mpc) . | Number of particles Np . | Particle mass mp ( h−1 M⊙) . | Force softening ϵ ( h−1 kpc) . | Smoothing length hsi ( h−1 kpc) . | Cross-section σ/m ( cm2 g−1) . |
---|---|---|---|---|---|---|
CDM-50 | 50 | 5123 | 6.88 × 107 | 1.0 | – | 0 |
CDM-25 | 25 | 5123 | 8.59 × 106 | 0.4 | – | 0 |
CDM-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | – | 0 |
CDM-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | – | 0 |
SIDM0.1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 0.1 |
SIDM0.1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 0.1 |
SIDM0.1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 0.1 |
SIDM0.1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 0.1 |
SIDM1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 1 |
SIDM1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 1 |
SIDM1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 1 |
SIDM1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 1 |
Name . | Volume LBox ( h−1 Mpc) . | Number of particles Np . | Particle mass mp ( h−1 M⊙) . | Force softening ϵ ( h−1 kpc) . | Smoothing length hsi ( h−1 kpc) . | Cross-section σ/m ( cm2 g−1) . |
---|---|---|---|---|---|---|
CDM-50 | 50 | 5123 | 6.88 × 107 | 1.0 | – | 0 |
CDM-25 | 25 | 5123 | 8.59 × 106 | 0.4 | – | 0 |
CDM-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | – | 0 |
CDM-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | – | 0 |
SIDM0.1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 0.1 |
SIDM0.1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 0.1 |
SIDM0.1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 0.1 |
SIDM0.1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 0.1 |
SIDM1-50 | 50 | 5123 | 6.88 × 107 | 1.0 | 2.8ϵ | 1 |
SIDM1-25 | 25 | 5123 | 8.59 × 106 | 0.4 | 2.8ϵ | 1 |
SIDM1-Z11 | (3Rvir)a | 2.5 × 106a | 1.07 × 106a | 0.3 | 2.8ϵ | 1 |
SIDM1-Z12 | (3Rvir)a | 5.6 × 107a | 1.34 × 105a | 0.1 | 1.4ϵ | 1 |
aThe Z11 and Z12 runs are zoom simulations with multiple particle species concentrating on haloes of mass Mvir = 5 × 1011 and 1.0 × 1012 M⊙, respectively (no h). The volumes listed refer to the number of virial radii used to find the Lagrangian volumes associated with the zoom. The particle properties listed are for the highest resolution particles only.
Each of our four initial conditions has been evolved from redshift z = 250 to 0 with collisionless dark matter (labelled CDM) and with two types of SIDM: one with σ/m = 1 cm2 g−1 (labelled SIDM1) and another with σ/m = 0.1 cm2 g−1 (labelled SIDM0.1). We can use the same initial conditions for CDM and SIDM because at high redshift the low relative velocities of the dark matter make self-interactions insignificant. Table 1 lists all the simulations used for this study and details their force, mass and self-interaction resolution. In addition to the simulations listed in Table 1, we also ran the cosmological boxes with SIDM cross-sections σ/m = 0.03 cm2 g−1. We do not present results from these low cross-section runs here because no core density differences were resolved within the numerical convergence radii of our simulations.
As shown in Section 3, the self-interaction smoothing length hsi must be larger than 20 per cent of the interparticle separation in order to achieve convergence on the interaction rate. All the work for this paper was done with a fixed hsi for all particles carefully chosen for each simulation so that the self-interactions are well resolved at densities a few times to an order of magnitude lower than the lowest densities for which self-interactions are significant. We have run the cosmological boxes with different choices for hsi (changes by factors of 2–4) and have found that our results are unaffected. We have also run tests on isolated haloes with varying smoothing lengths and again find that the effects of self-interactions are robust to reasonable changes in hsi.
All of our halo catalogues and density profiles are derived using the publicly available code Amiga Halo Finder (ahf; Knollmann & Knebe 2009).
5 SIMULATION RESULTS
5.1 Preliminary illustrations
Before presenting any quantitative comparisons between our CDM and SIDM runs, we provide some simulation renderings in order to help communicate the qualitative differences.
The upper panels of Fig. 2 show a large-scale comparison: two (50 × 50 × 10) h−1 Mpc slices from the CDM-50 and SIDM1-50 boxes side-by-side at z = 0. The structures are colour-coded by local phase-space density (∝ ρ/v3rms). It is evident that there are no observable differences in the large-scale characteristics of CDM and SIDM1. We discuss this result in more quantitative terms in Section 5.2 but of course this is expected. The SIDM models we explore do not have appreciable rates of interaction for densities outside the cores of dark-matter haloes. The upper panels of Fig. 2 provide a visual reminder that the SIDM models we consider are effectively identical to CDM on large scales.
The differences between CDM and SIDM become apparent only when one considers the internal structure of individual haloes. The lower panels of Fig. 2 provide side-by-side images of a Milky Way mass halo (Z12) simulated with CDM (left) and SIDM1 (right). SIDM tends to make the cores of haloes less dense and kinetically hotter (see Section 5.3) and these two differences are enhanced multiplicatively in the phase-space density renderings. The central regions of the host halo are also slightly rounder in the SIDM case (Peter et al. 2012). Importantly, the difference in substructure characteristics is minimal, especially at larger radii. We return to a quantitative description of substructure differences in Section 5.4.
5.2 Large-scale structure and halo abundances
Fig. 3 provides a quantitative comparison of both the clustering properties (left) and halo abundance evolution (right) between our full-box CDM and SIDM1 simulations. The left-hand panel shows the two-point function of dark-matter particles in both cosmological runs for CDM and SIDM1. There are no discernible differences between SIDM and CDM over the scales plotted, though of course the different box sizes (and associated resolutions) mean that the boxes themselves only overlap for a limited range of scales. For a given set of initial conditions, however, SIDM and CDM give identical results.
The right-hand panel of Fig. 3 shows the cumulative number density of dark-matter haloes (including subhaloes) as a function of their peak circular velocity (Vmax) for the CDM-50 (solid) and SIDM1-50 (dashed) simulations at various redshifts. Remarkably, this comparison shows no significant difference either – indicating that SIDM with cross-sections as large as 1 cm2 g−1 does not strongly affect the maximum circular velocities of individual haloes. The two panels of Fig. 3 demonstrate that for large-scale comparisons, including analyses involving field halo mass functions, SIDM and CDM yield identical results. The implication is that observations of large-scale structure are just as much a ‘verification’ of SIDM as they are of CDM.
5.3 Halo structure
Before presenting statistics on halo structure, we focus on six well-resolved haloes that span our full mass range Mvir = 5 × 1011-2 × 1014 M⊙, selected from our full simulation suite, including our two zoom-simulation haloes (Z12 and Z11). Figs 4– 6 show radial profiles for the density, circular velocity and velocity dispersion for all three dark matter cases. In each figure, black circles correspond to CDM, green triangles to SIDM0.1 and blue stars to SIDM1. All profiles are shown down to the innermost resolved radius for which the average two-body relaxation time roughly matches the age of the Universe (Power et al. 2003).
The SIDM versions of each halo show remarkable similarity to their CDM counterparts at large radii. However, the SIDM1 cases clearly begin to roll towards constant-density cores at small radii. The best resolved haloes in the SIDM0.1 runs also demonstrate lower central densities compared to CDM, though the differences are at a factor of ∼2 level even in our best resolved systems. Clearly, higher resolution simulations will be required in order to fully quantify the expected differences between CDM and SIDM for σ/m ∼ 0.1 cm2 g−1.
As explained in Section 7, the fact that the SIDM1 profiles are reasonably well characterized by a single scale-radius Burkert profile may be a lucky accident, only valid for cross-sections near 1 cm2 g−1. It just so happens that for this cross-section the radius where dark-matter particles experience significant scattering sets in at r ∼ rs (see Fig. 7 and related discussion). For a smaller cross-section (with a correspondingly smaller core) a multiple parameter fit may be necessary. Given the beginnings of very small cores we are seeing in the SIDM0.1 runs, it appears that we would need one scale radius to define an rs bend and a second scale radius to define a distinct core.
Another qualitative fact worth noting is that the density profiles of the SIDM1 haloes overshoot the CDM density profiles near the Burkert core radius (not as much as the Burkert fits do, but the difference in the data points is noticeable). This is due to the fact that as particles scatter in the centre, those that gain energy are pushed to larger apocentre orbits. This observation invites us to consider a toy model for SIDM haloes where the effect of SIDM is confined to a region (smaller than a radius of about rb) wherein particles redistribute energy and move towards a constant-density isothermal core. We will develop this model further to explain the scaling relations between core size and halo mass in Section 6.
The circular velocity curves for the same set of haloes discussed above are shown in Fig. 5. The SIDM rotation curves rise more steeply and have a lower normalization than for CDM within the NFW scale radius rs. This brings to mind the rotation curves observed for low surface brightness galaxies and we will explore this connection later. Note though that the peak circular velocity Vmax actually is slightly higher for the SIDM1 case because of the mass rearrangement (evident in the density profiles in Fig. 4) briefly discussed in the last paragraph. At radii well outside the core radius, the rotation curves of the CDM and SIDM1 haloes converge, though this convergence occurs beyond the plot axes > rs for most of the haloes shown.
An appreciation of why the density profiles of SIDM haloes become cored can be gained from studying their velocity dispersion profiles compared to their CDM counterparts, as illustrated in Fig. 6. Here, vrms is defined as the root mean square speed of all particles within radius r. While the CDM haloes (black) are colder in the centre than in their outer parts (reflecting a cuspy density profile) the SIDM haloes have hotter cores, indicative of heat transport from the outside in. Moreover, the SIDM haloes are slightly colder at large radii, again reflecting a redistribution of energy. As discussed in the introduction, it is this heat transport that is the key to understanding why CDM haloes differ from SIDM haloes in their density structure (Balberg et al. 2002; Colín et al. 2002; Ahn & Shapiro 2005; Koda & Shapiro 2011). The added thermal pressure at small radii is what gives rise to the core. The SIDM1 simulations have sufficient interactions that they have been driven to isothermal profiles for r/rs ≲ 1, while for SIDM0.1 the vrms profiles typically begin to deviate from the CDM lines only at smaller radii, r/rs ∼ 0.2, reflecting the relatively lower scattering rate.
By comparing Fig. 7 to Fig. 6 (and to some extent to all Figs 4–6), we see that the effects of self-interactions do become evident at radii corresponding to ρ vrms ∼ 0.1 for SIDM1 (at r/rs ∼ 0.8) and ρ vrms ∼ 1 for SIDM0.1 (at r/rs ∼ 0.2). Interestingly, for the SIDM1 haloes this interaction radius is fairly close to the Burkert scale radius (shown by the blue arrows). It should be kept in mind, however, that the structure of haloes can be affected at larger radii because particles scattering in the inner regions can gain energy and move to larger orbits. A careful inspection of the density and rotation velocity profiles shows that this is indeed the case.
We will discuss these findings in more detail in Sections 6 and 7. In particular, in Section 7 we present an analytic model aimed at understanding how the central densities and scale radii of SIDM haloes are set in the context of energetics. Moreover, before moving on to those issues, we first explore halo substructure in SIDM.
5.4 Substructure
The question of halo substructure is an important one for SIDM. One of the original motivations for SIDM was to reduce the number of subhaloes in the Milky Way halo in order to match the relative dearth of observed satellite galaxies (Spergel & Steinhardt 2000). However, the overreduction of halo substructure is now recognized as a negative feature of SIDM compared to CDM, given the clear evidence for galaxy-size subhaloes throughout galaxy clusters (Natarajan et al. 2009) and the new discoveries of ultrafaint galaxies around the Milky Way (see Bullock 2010; Willman 2010 for reviews). In fact, one of the most stringent constraints on the self-interaction cross-section comes from analytic subhalo-evaporation arguments (Gnedin & Ostriker 2001).
Fig. 8 demonstrates that the effects of subhalo evaporation in SIDM are not as strong as previously suggested on analytic grounds. Here we show the cumulative number of subhaloes larger than a given Vmax for a sample of well-resolved haloes in our CDM (solid), SIDM0.1 (dotted) and SIDM1 (dashed) simulations. The associated virial masses for each host halo are shown in the artwork. The left-hand panel presents the Vmax function for all subhaloes within the virial radius of each host and the right-hand panel restricts the analysis to subhaloes within half of the virial radius. We see that generally the reduction in substructure counts at a fixed Vmax is small but non-zero and that the effects appear to be stronger at small radii than large. Similarly, there appears to be slightly more reduction of substructure in the SIDM cluster haloes compared to the galaxy size systems.
While the increase in destruction of subhaloes with host halo mass is not strong, it is clear from the above arguments that subhaloes in the inner parts of the halo (r/rs ≪ 1) should be destroyed but the bulk of the subhaloes around r/rs ∼ 1 and beyond should survive for σ/m = 1 cm2 g−1. This effect is strengthened by the fact that subhaloes in the innermost region of the halo were accreted much longer ago than subhaloes in the outskirts, so they have experienced many more orbits (Rocha, Peter & Bullock 2012). These arguments explain the comparisons between the subhalo mass functions plotted in Fig. 8. Our arguments demonstrate that a large fraction of the subhaloes found in CDM haloes (most of which are in the outer parts) would still survive in SIDM haloes for σ/m values around or below 1 cm2 g−1.
Overall in the previous two sections we have seen that the effects of self-interactions between dark-matter particles in cosmological simulations are primarily in the central regions of dark-matter haloes, leaving the large-scale structure identical to our non-interacting CDM simulations. Thus, we retain the desirable features of CDM on large scales while revealing different phenomenology near halo centres. In the following section, we will move to explore how the properties of SIDM haloes presented here scale with halo mass.
6 SCALING RELATIONS
In the previous section we saw that while SIDM preserves the CDM large-scale properties of dark-matter haloes, self-interactions in the central regions of haloes result in a decrease of central densities and the formation of cores in their density profiles. We found that the density profiles of haloes from our SIDM1 simulations can be relatively well fit by Burkert density profiles inside r ∼ 2-3rs (see Fig. 4). Here we define a sample of well-resolved haloes from all our SIDM1 simulations and use Burkert fits to their density profiles in order to quantify their central densities and core sizes. We then provide scaling relations of dark-matter halo properties with maximum circular velocity Vmax.
The sample of haloes used for the rest of this section consists of the two host haloes in our SIDM1-Z11 and SIDM1-Z12 simulations together with the 25 most massive haloes from our SIDM1-50 and the 25 most massive haloes from our SIDM1-25 simulations. That gives us a total of 52 haloes spanning a range Vmax = 30-860 km s− 1 or Mvir = 5 × 1011-2 × 1014 M⊙. For this set of haloes the innermost resolved radius, defined by equation (20) in Power et al. (2003), is always smaller than one-third of the Burkert scale radius from which we define the sizes of cores. It is vital that we do a conservative comparison to the Power et al. (2003) radius because both gravitational scattering and self-interactions lead to the same phenomenological result of constant-density cores. Most of the haloes (other than the 52 we select here) do not pass this test well enough for the core set by self-interactions to be resolved with confidence. This desire to be conservative in our presentation of scaling relations forces us to find these relations from only a small sample of haloes for SIDM1 and leaves us with basically no haloes to find scalings for SIDM0.1. Moreover, one has to keep in mind that our SIDM1 relations could be biased by selecting only the most massive haloes in our full-box simulations. Evidently higher resolution simulations are necessary to find definitive answers. It is reassuring however that the scaling relations derived from our analytical arguments in Section 7 agree so well with the ones presented here for σ/m = 0.1 cm2 g−1.
We have checked that for all of our haloes we resolve the scattering rate out to at least four times the Burkert scale radius. Outside of this point the scattering rate is underestimated because of our choice of the self-interaction smoothing length relative to the interparticle spacing (see Section 3). However, the expected scattering rate is negligible with respect to the Hubble rate outside that radius (Fig. 7). Moreover, we have re-run our 50 h−1 Mpc boxes for a range of SIDM smoothing values and found identical results. Thus, we consider our sample to be well resolved.
Eight haloes in our sample are undergoing significant interactions and have density profiles that are clearly perturbed even in the CDM runs. We include these eight systems in all of the following plots but indicate them with open symbols. We do not use them in the best fits for the scaling relations that we provide.
We start by examining the global structure of haloes as characterized by the maximum circular velocity Vmax and the radius where the rotation curve peaks, rmax. The relationship between Vmax and rmax provides a simple, intermediate-scale measure of halo concentration and we aim to investigate any differences between SIDM and CDM. Fig. 9 shows the Vmax-rmax relation for CDM (black) and SIDM1 (blue) haloes. We can see that small differences of about 10 per cent exist in both Vmax and rmax, with SIDM1 haloes having larger values for Vmax and smaller for rmax. This was already evident in Fig. 5, where the circular velocity curves of SIDM1 haloes seem to peak at slightly smaller radii and slightly larger velocities than their CDM analogues, even though SIDM1 curves decrease more steeply at the centre.
In this section, we have presented scaling relations for the properties of haloes in our SIDM1 simulations. Our limited resolution allows us to use only 52 haloes spanning a modest mass range, from which we throw out eight systems that are undergoing mergers. Admittedly, this sample is not large enough to be definitive, especially in regards to scatter. However, the strong correlation between the SIDM core radius rb and the counterpart CDM scale radius rs is clearly statistically significant and the general trends provide a useful guide for tentative observational comparisons – a subject we will return to in the final section below.
7 ANALYTIC MODEL TO EXPLAIN THE SCALING RELATIONS
In this section, we develop a simple model to understand the scaling relations shown in Section 6. This model is based on identifying an appropriate radius r1 within which self-interactions are effective and demanding that the mass as well as the average velocity dispersion within this radius are set by the mass and the average velocity dispersion (within the same radius) of the same halo in the absence of self-scatterings. The mass loss due to scatterings in the core should be insignificant because particles rarely get enough energy to escape and this implies that the mass within r1 should be close to what it would have been in the absence of self-interactions. This also implies that the potential outside r1 is unchanged from its CDM model prediction, but tends to a constant value faster inside r1. Within this set of approximations, the dominant effect due to scatterings is to redistribute kinetic energy in the core, while keeping the total kinetic energy within r1 the same as it would have had before self-interactions became important. We have looked at the kinetic energy profiles in the best-resolved haloes in our simulations and have confirmed that this is indeed a good approximation. Note that in this picture, there is a clear demarcation of time-scales such that the inner halo structure (say r ≲ rs) is set (the same way as in the CDM model) well before self-interactions become important. For cross-sections much larger than what we are interested in here, this need not hold.
We then set v2rms, 0 equal to the average velocity dispersion squared (i.e. two times kinetic energy divided by mass) within the region r1 in the absence of self-interactions. This basically demands that the kinetic energy within r1 is unchanged from the value it would have had in the absence of self-interactions. Note, however, that we are setting the average velocity dispersion squared equal to v2rms, 0 and not the corresponding average in the SIDM halo. This is an approximation, but one that is degenerate with choosing the ξ parameter.
To finish specifying this model, we need a density profile for the region inside r1. A Burkert profile has a velocity dispersion profile (assuming isotropy) that asymptotes very slowly to the central dispersion. For small radii, the radial dispersion profile is slowly increasing (with radius) because of the r/rb term in the Taylor expansion for the density profile. If we want a flatter central dispersion profile (as is observed for the SIDM1 haloes), we can fix this by either assuming an isothermal profile or something like 1/[1 + 1.52(r/r0)2]3/2. The final results turn out to be qualitatively similar for these profiles. Hence we adopt a Burkert profile for ease of comparison to the fits presented here and then check the results with more appropriate profiles later. Our two constraints (on the radial velocity dispersion and mass) fully specify the density and radial scales of the Burkert profile.
For the density profile in the absence of self-interactions, we assume an NFW profile and to fix the velocity dispersion we use the observed fact that the phase-space density is a power law in radius (Taylor & Navarro 2001). By noting that vrms, CDM(r) = [ρCDM(r)/Q(r)]1/3 and using a phase-space density profile Q(r) = Q(rs)(r/rs)−η (Taylor & Navarro 2001; Ascasibar et al. 2004; Rasia, Tormen & Moscardini 2004; Dehnen & McLaughlin 2005; Ascasibar & Gottlöber 2008), we may fully specify the dependence of r1 on the cross-section and halo parameters (say Vmax and rmax). For the phase-space density profile, we use a power-law index η = 2 and Q(rs) = 0.3/(GVmaxr2max) derived from jointly fitting our relaxed CDM haloes; these parameters are very similar to the fits provided in Ascasibar & Gottlöber (2008).
Let us first look at how r1 scales with rs in the NFW density profile. One notes that ρs = 1.72V2max/(Gr2max) and hence ρsVmax∝V3max/r2max which is a very mildly increasing function of Vmax as equation (15) shows. Thus, equation (22) implies that r1/rs should be roughly a constant. Numerically, we find that r1/rs ≃ 0.7-0.8 over the range of Vmax of interest for σ/m = 1 cm2 g−1.
In detail, the model predicts that rb/rs = 0.5-0.6 for dwarf to cluster haloes in good agreement with the fits to our SIDM1 haloes, but about 25 per cent smaller for Vmax ∼ 100 km s−1. It departs from the results of the simulation in predicting that rb/rs increases gently with Vmax, whereas Fig. 11 predicts that this ratio should decrease gently with Vmax. We find that this departure from simulations is likely related to the assumption of a constant age for all haloes. To generalize our model, we use the results of Wechsler et al. (2002) who show that the virial concentrations of haloes are correlated with their formation times, and in particular cvir = 4.1(1 + zform) for a particular definition of formation time. We invert this equation to derive an estimate of the halo age using zform. With the age thus specified in equation (22), we find that now rb/rs decreases gently with Vmax in substantial agreement with the fit to our simulations. Thus, the reason that larger haloes have a smaller rb/rs is because self-interactions have had less time to operate. We note that the values for the core radius in the analytic model with halo mass-dependent tage are uniformly about 25 per cent smaller, but this should not be a cause for concern given the approximation in demanding a sharp transition at r1.
Given the Burkert core radius rb and the central velocity dispersion vrms, 0, one can easily check that the central density ρb is about 0.01 M⊙ pc− 3 for Vmax = 300 km s−1 haloes and 0.005 M⊙ pc− 3 for Vmax = 1000 km s−1 in this analytic model. These numbers and the scaling with Vmax for ρb (when including the halo mass-dependent tage) are in good agreement with the densities in Fig. 12 and the fit in equation (19). As we have indicated before, the scaling relation for the central density should be interpreted with care given the large scatter. Given the tight correlation between core radius and rs, it is possible that the substantial scatter in the central density arises in large part due to the scatter introduced by the assembly history in the concentration–mass relation. This has important implications for fitting to the rotation velocity profiles of low surface brightness spirals (LSBs; Kuzio de Naray et al. 2010) and deserves more work.
The simple model constructed above also provides insight into the core-collapse time-scales. In particular, as long as the outer part (region outside r1) dominates the potential well and sets the average central temperature (or the total kinetic energy in the core), we do not expect core collapse. This is simply because core collapse requires uncontrolled decrease in temperature, which is prohibited here. Once r1 moves out well beyond rmax or to the virial radius, there is significant loss of particles and core collapse may occur if there are no further major mergers. The time-scale for this process is much longer than the age of the Universe for σ/m = 1 cm2 g−1 because the inner core is at r1 < rs after 10 Gyr for this self-interaction strength and we see no evidence for significant mass loss.
8 OBSERVATIONAL COMPARISONS
The goal of this section is to discuss our results in comparison to observationally inferred properties of dark-matter density profiles. In particular, we will focus on the core densities and core sizes. Section 8.1 presents our expectations for SIDM1 and SIDM0.1. Our predictions for σ/m = 1 cm2 g−1 are anchored robustly to our simulations, though they do require some extrapolation beyond the mass range directly probed by our simulations (Vmax = 130-860 km s− 1). For σ/m = 0.1 cm2 g−1 the predictions are much less secure because the associated core sizes are of the order of our resolution limit; thus, we rely on our analytic model more directly here. In Section 8.2, we discuss our predictions in light of observations of dark-matter haloes for a wide range of halo masses. In Section 8.3, we discuss our results on subhaloes in the context of past work and constraints on SIDM based on subhalo properties.
Before proceeding with this discussion we would like to clarify how we quantify core sizes. In this work, we have fit the σ/m = 1 cm2 g−1 haloes with Burkert density profiles. However, many observational constraints on cores on galaxy scales come from fitting pseudo-isothermal density profiles with core size rpi to data (e.g. Simon et al. 2005; Kuzio de Naray et al. 2008), although some constraints do come from Burkert modelling (Salucci et al. 2012). We found that pseudo-isothermal density profiles also give good fits to the inner regions of the SIDM1 haloes, but Burkert fits are better because of that profile’s ρ ∝ r−3 dependence at large radii. For a pseudo-isothermal density profile [∝1/(r2c + r2)], the density decreases to one-fourth the central density at 1.73 times its core radius rc. Thus, as a crude approximation, one may convert the Burkert radius to the equivalent pseudo-isothermal core radius by multiplying by a factor of 0.58 (rc ≃ rb/1.73).
8.1 Predicted core sizes and central densities in SIDM
8.1.1 SIDM with σ/m = 1 cm2 g−1
The central properties of dark-matter haloes have been inferred from observations from tiny Milky Way dSph galaxies (Vmax ≲ 50 km s−1) to galaxy clusters (Vmax ≳ 1000 km s−1). If we extrapolate the results from our set of SIDM1 simulations using equations (16)–(20), we predict that SIDM haloes with σ/m = 1 cm2 g−1 would have the following (Burkert) core sizes and central densities:
For galaxy clusters (Vmax ≃ 700–1000 km s−1),
For low-mass spirals (Vmax ≃ 50-130 km s− 1),
For dSph galaxies (Vmax ≃ 20-50 km s− 1),
Although we cannot completely determine the scatter in our scaling relations due to low number statistics, it is important to note from Figs 10 and 12 that a scatter of at least a factor of 2 in core sizes, and at least a factor of 3 in central densities, is expected for a given Vmax. We suspect that these differences are in large part a result of the diversity of merger histories of dark-matter haloes. Note that the Vmax-rmax and Q(rs) scalings assumed in the analytic model are the median values. The strong dependence of the SIDM halo profiles on these quantities makes it clear that the scatter in these relations will introduce significant scatter in the halo core sizes and core densities. Thus, the analytic model should also provide a simple way to understand (some of the) scatter seen for SIDM1 halo properties. In future work, we will characterize the relation between the core properties and merger history in the context of a detailed discussion of the scatter in the scaling relations, especially on scales that we do not resolve with our current simulations.
8.1.2 SIDM with σ/m = 0.1 cm2 g−1
As discussed in Section 5.4, our SIDM0.1 simulations are not well enough resolved to definitively measure a core radius for any of our haloes, much less define scatter in that quantity. Nevertheless, our best resolved systems do demonstrate some clear deviations from CDM and allow us to cautiously estimate individual core densities. Referring back to Fig. 4, we see that in our two best resolved cluster haloes (at Mvir ≃ 1014 M⊙) the SIDM0.1 core densities approach ∼ 0.01 M⊙ pc− 3 – each at least a factor of ∼3 denser than their SIDM1 counterparts. Similarly, in our Z12 Milky Way case, the SIDM0.1 core density appears to be approaching ∼ 0.1 M⊙ pc− 3 compared to ∼ 0.02 M⊙ pc− 3 in the SIDM1 case.
Given the lack of well-resolved halo profiles, it is worth appealing to the analytic model presented in Section 7 to estimate core radii for SIDM0.1. Using exactly the same arguments (including the halo mass-dependent age), we find that r1/rs ≃ 0.05-0.12 in the 100–1000 km s−1Vmax range and a corresponding Burkert core radius rb/rs ≃ 0.09-0.17. We note that the Burkert radius is close to but slightly larger than r1. It is important to keep in mind that in this analytic model we are only explicitly fitting the inner ‘self-interaction zone’ of r < r1. This does not imply that the entire halo has to be well fit by the Burkert profile. Recall that a single-scale Burkert profile only works as well as it does for σ/m = 1 cm2 g−1 because rb ≈ rs, such that to a good approximation there is only one relevant length scale. For the smaller cross-section that we are now considering we expect the core and NFW scale radii to be widely separated, suggesting that a generic functional form for SIDM haloes should have two scale radii. A wide separation between the SIDM0.1 core and rs does appear to be consistent with the highest resolved haloes presented in Fig. 4. However, we note that given the strong correlation between rb/rs, we still expect a one-parameter family of models for a given σ/m.
To see how dependent our results are on the shape of the inner halo profile, we modify the analytic model to include a density profile that decreases with radius as 1/[1 + (r/rc)2]α/2. For this density profile, the velocity dispersion profile has the right form to match our simulation results. The price we pay is the introduction of a new parameter α. We set this parameter α by additionally demanding that the slope of the mass profile (i.e. density) is continuous at r1, so that the mass profile joins smoothly with the NFW mass profile. This picks out a narrow range α = 5.5-7.0 as the solution over most of the Vmax range of interest (with smaller values corresponding to lower Vmax). Interestingly, this implies that at r1, the slope of the density profile is very close to −2 for the entire range of Vmax values of interest. Note that while the mass profile is continuous, the slope of the density profile is not matched smoothly at r1 (since the slope of the NFW profile would be closer to −1 at r1 ≪ rs). This probably signals that if the matching were not done sharply (at r1), the density profile of SIDM would overshoot that of CDM and catch up at some radius beyond r1 (as is seen in the comparison of SIDM1 and CDM density profiles).
As a check we apply this α model to the σ/m = 1 cm2 g−1 case and find that the results are qualitatively the same as the model with the Burkert profile. The quantitative differences are at 20 per cent level with the densities being smaller and inferred Burkert core radii (where density is one-fourth of the central density) larger compared to the Burkert profile model. The predicted slope of the density profile at r1 is close to −2.5 implying a smoother transition to the NFW profile (since r1 ∼ rs for σ/m = 1 cm2 g−1), as is seen Fig. 4.
For the σ/m = 0.1 cm2 g−1 case, we obtain rc/rs = 0.08-0.17 and an equivalent Burkert core radius (where the density is one-fourth of the central density) rb/rs = 0.06-0.14, in substantial agreement with the results we obtained using the Burkert profile. Thus, our analysis would suggest core sizes ∼0.1rs for σ/m = 0.1 cm2 g−1. The results from the analytic model for σ/m = 0.1 cm2 g−1 also seem consistent with our simulation results; see Fig. 6 where the vrms profiles for SIDM0.1 start to deviate from CDM at ∼ 0.2rs.
Based on the discussion above, we conclude that for σ/m = 0.1 cm2 g−1 we expect the following:
While these values are somewhat tentative compared to those presented above for SIDM1 (given our lack of direct simulation fits), two factors are reassuring. First, the analytic model is based on the simple assumption that scattering redistributes kinetic energy within the inner halo and the non-trivial aspect of the model is defining this ‘inner halo’ region. There is no reason to suspect that this assumption or the prescription breaks down for SIDM0.1 haloes when it works so well in describing the SIDM1 haloes. The predicted densities are in line with those inferred for the best resolved haloes in our SIDM0.1 simulations (shown in Fig. 4 and discussed above). For the core radii, we reiterate that the label ‘rb’ should be interpreted (according to its definition in the analytic model) as the radius where the density reaches one-fourth the asymptotic core density. The overall profile of a halo with such a small core compared to rs will not be fit by the Burkert form. Note that the strong correlations we predict between the core radius and the NFW scale radius raise the intriguing possibility that the SIDM haloes may be also well fit (modulo scatter) by a single parameter profile as is the case for CDM.
Next, we compare our predictions for SIDM core properties against data and show that the core radii and densities appear to be consistent with that seen in real data, motivating future simulations with high enough resolution to resolve cores in SIDM0.1 haloes.
8.2 Observed core sizes and central densities versus SIDM
In this section, we explore the predictions for the properties of density profiles with SIDM in the context of observational constraints on density profiles. We also revisit previous constraints on SIDM from observations in light of our simulation suite.
8.2.1 Clusters
One of the tightest SIDM constraints from the first generation of SIDM studies emerged from one cluster simulation and one observed galaxy cluster. Specifically, Yoshida et al. (2000) simulated an individual galaxy cluster with different SIDM cross-sections. When comparing the core size of this simulated cluster to the core size estimated by Tyson, Kochanski & dell’Antonio (1998) for CL 0024+1654, they found that the observed core in CL 0024+1654 would be consistent with SIDM only if σ/m ≲ 0.1 cm2 g−1. Since that time, evidence has emerged that this particular cluster is undergoing a merger along the line of sight (Czoske et al. 2001, 2002; Zhang et al. 2005; Jee et al. 2007; Jee 2010; Umetsu et al. 2010). Thus, this cluster is not the ideal candidate for SIDM constraints based on the properties of relaxed haloes, and the Yoshida et al. (2000) constraint is not valid in this context.
Using X-ray emission, weak lensing, strong lensing, stellar kinematics of the brightest cluster galaxy (BCG) or some combination thereof, the mass distributions within a number of galaxy clusters have been mapped in the past decade. Arabadjis, Bautz & Garmire (2002) placed a conservative upper limit of 75 kpc on the size of any constant-density core, and an average density within the inner 50 kpc of ∼ 0.025 M⊙ pc− 3 for a halo with an estimated mass M ∼ 4 × 1014 M⊙.
Sand et al. (2004), Sand et al. (2008), Newman et al. (2009) and Newman et al. (2011) all find central density profiles in clusters shallower than the NFW CDM prediction. The difference in the work between these authors and others is that they use stellar kinematics of the BCG to constrain the density profile of the cluster dark-matter halo on small scales. While this probe of the density profile is more sensitive on small scales than strong lensing is, proper inference of the dark halo properties depends on accurate modelling of the BCG density profile and equilibrium structure. They have typically assumed a ‘gNFW’ profile in order to constrain the central densities: ρ(r) ∝ 1/[xg(1 + x)3 − g] with r = xrs and the NFW form obtained when g = 1. The Newman et al. (2009, 2011) mass models of M ∼ 1015 M⊙ clusters show average dark-matter central densities within 10 kpc of ∼ 0.03-0.06 M⊙ pc− 3 and rs of the order of 100 kpc. Note that 10 kpc is typically the smallest radius our simulations can resolve.
Saha, Read & Williams (2006) and Saha & Read (2009) studied the mass structure of three cluster haloes from gravitational lensing and obtained density profiles that are consistent with ρ ∝ r−1 outside the inner 10-20 kpc regions. Similarly, Morandi, Pedersen & Limousin (2010) and Morandi & Limousin (2012) find that the radial mass distribution of cluster dark-matter haloes is consistent with NFW predictions outside 30 kpc in projection. The CLASH multi-cycle treasury programme on the Hubble Space Telescope is finding many new strongly lensed galaxies in about a set of 25 massive clusters (Postman et al. 2012). Initial results from this programme show that the total density profile of these clusters (or total density minus the BCG), if modelled as spherically symmetric, is consistent with NFW predictions for the halo alone if the gNFW functional form is used in the fit (Zitrin et al. 2011; Coe et al. 2012; Umetsu et al. 2012). However, Morandi et al. (2010) and Morandi, Pedersen & Limousin (2011) find that spherical mass modelling of galaxy clusters typically results in an overestimate of the cuspiness of the density profile, although axially symmetric modelling is found to lead to underestimates (Meneghetti et al. 2007). Thus, the present status of the density profiles of the CLASH clusters is unclear and clearly an interesting data set to look forward to.
We note here a complexity involved in using the lensing results to constrain SIDM models. Lensing provides mass in cylinders along the line of sight and this 2D mass profile is sensitive to mass from a large range of radii. As an example, let us consider mass within 30 kpc in projection. If we were to do something extreme and create a zero-density core inside 30-kpc sphere, the differences in the 2D mass profile would be less than a factor of 2 for clusters in the 1014-1015 M⊙ mass range. For SIDM0.1, the differences are comparatively benign. Our analytic model predicts that differences relative to CDM at about 0.1rs (which is 10-40 kpc for 1014-1015 M⊙ virial mass range) are 20–30 per cent, which implies that SIDM0.1 surface mass density profiles are very similar to CDM on these scales. Moreover, for SIDM1 the expected differences would be measurably large.
On a related technical note, we discourage the use of the gNFW functional form when thinking about models that deviate from the CDM paradigm. In the SIDM case, for σ/m < 1 cm2 g−1, there will generically be two scale radii: one is the NFW-like scale radius which is the result of hierarchical structure formation (Lithwick & Dalal 2011) and the second is the core radius from dark-matter self-scattering. For σ/m = 1 cm2 g−1, as we explained in detail in Section 7, the two scales are about the same. If most of the cluster data constrain the density profile beyond an SIDM core, as they may for weak lensing and X-ray studies, the gNFW or NFW fit is dominated by those data, and a core will not be ‘detected’ in the fit. In future work, we will simulate haloes with a broader range of σ/m and provide SIDM-inspired density profiles to the community.
The results discussed above seem to suggest that the density profile beyond about 25 kpc should be close to the predictions from the NFW profile. To test this we plot the average physical density within 25 kpc for well-resolved haloes in our CDM (black), SIDM0.1 (green) and SIDM1 (blue) simulations in Fig. 13. We see that for the most massive haloes, the σ/m = 1 cm2 g−1 run produces densities at 25 kpc that are ∼ 2-3 times lower than their CDM counterparts. Thus, it seems like the measured densities in clusters rule out the σ/m = 1 cm2 g−1 SIDM model. At the same time, the σ/m = 0.1 cm2 g−1 simulations are quite similar to CDM at these radii, though beginning to show some differences as we discussed earlier in this section. Analyses that combine information from X-rays, lensing and BCG stellar kinematics seem to suggest lowered densities (e.g. Newman et al. 2011) that would be compatible with SIDM0.1. Given this outlook, it is reasonable to conclude that estimates of the central dark-matter density in clusters will provide essential tests of interesting SIDM models.
8.2.2 Low-mass spirals
For low-mass spirals with maximum circular velocities in the range 50-130 km s− 1, constant-density cores with sizes of ∼ 0.5-8 kpc and central densities of ∼ 0.01-0.5 M⊙ pc− 3 have been observed (de Blok et al. 2001; Sánchez-Salcedo 2005; Simon et al. 2005; Kuzio de Naray et al. 2008, 2010; Oh et al. 2011a; Salucci et al. 2012). Similar to what we found for cluster scales, SIDM with σ/m = 1 cm2 g−1 would be able to reproduce the largest core sizes observed in low-mass galaxies but it predicts central densities that are too low. SIDM with σ/m = 0.1 cm2 g−1 would be much more consistent. Moreover, the predicted log-slope of the density profile at 500 pc for σ/m = 0.1 cm2 g−1 haloes in the 50–130 km s−1 range is −0.5 to 0, both facts consistent with results from THINGS (Oh et al. 2011a). Note that the slope at 500 pc for the σ/m = 1 cm2 g−1 model is 0 in the same Vmax range, which is not consistent with the scatter seen in the data.
We conclude, as before, that the observed densities and core radii are not consistent with SIDM1 but are fairly well reproduced in SIDM models with σ/m ≃ 0.1 cm2 g−1.
8.2.3 Dwarf spheroidals in the Milky Way halo
The least massive and most dark-matter-dominated galaxies provide an excellent setting to confront the predictions of different dark-matter models with observations. Recent work by Boylan-Kolchin et al. (2011, 2012) has found that the estimated central densities of the bright Milky Way dSph satellites are lower than the densities of the massive subhaloes in dark-matter-only simulations. SIDM offers a way to solve this problem because it reduces the central density of haloes. Thus in SIDM, the massive subhaloes do host the luminous dSph but have shallower density profiles than predicted in CDM simulations. This has recently been demonstrated by Vogelsberger et al. (2012). We do not directly compare to Vogelsberger et al. (2012) because their work is focused on the subhaloes of the Milky Way and the velocity-independent cross-section that they simulate (σ/m = 10 cm2 g−1) is larger than the cross-sections considered in our work.
Regardless of whether Milky Way dSphs have cuspy or cored dark-matter haloes, we may estimate the enclosed mass, and hence average density, around the half-light radius of the stellar distribution. Mass estimates within 300 pc and mass profile modellings using stellar kinematics together with chemodynamically distinct stellar subcomponents of Milky Way dSph galaxies suggest central densities of ∼ 0.1 M⊙ pc− 3 (Strigari et al. 2008; Wolf et al. 2010; Walker & Peñarrubia 2011; Amorisco & Evans 2012; Wolf & Bullock 2012). For the faintest dSph Segue 1, the density within the half-light radius (about 40 pc) is measured to be about 2.5+ 4.1− 1.9 M⊙ pc− 3 (Martinez et al. 2011; Simon et al. 2011). The errors on Segue 1 density are large but it is clear that if SIDM is to accommodate this result, it must allow for large scatter in the core sizes and densities for small Vmax haloes. With a factor of 2–3 scatter in the densities quoted earlier for SIDM0.1 haloes, Segue 1 would appear to be compatible with SIDM0.1 if its Vmax value is towards the lower end of the 20-50 km s− 1 range in Vmax.
For the two dSph galaxies that appear to have cored density profiles (Fornax and Sculptor), the core sizes must be of the order of ∼ 0.2-1 kpc (Walker & Peñarrubia 2011). For small haloes with circular velocities in the 20-50 km s− 1 range, which is close to the expected peak circular velocities of dSph haloes before infall into the Milky Way host halo, an SIDM with σ/m = 1 cm2 g−1 predicts core sizes of the order of ∼ 0.8-3.0 kpc, with central densities of ∼ 0.02-0.04 M⊙ pc− 3. Therefore, we find again that σ/m = 1 cm2 g−1 cannot reproduce the observed high central densities. On the other hand, our estimates suggest that an SIDM model with σ/m = 0.1 cm2 g−1 would produce central densities and core sizes consistent with the Milky Way dSph.
In this last section we have used the analytic results that explain the scaling relations for the core sizes and central densities of haloes in our SIDM1 and SIDM0.1 simulations, to extrapolate our results to scales ranging from galaxy clusters to dSph galaxies and to lower cross-sections. We have found that σ/m = 1 cm2 g−1 would be unable to reproduce the observed high central densities. Remarkably, we find that the observations should be consistent with the predictions of an SIDM with cross-section in the ballpark of σ/m = 0.1 cm2 g−1. These expectations are based on the scaling relations seen in SIDM1 simulations and our analytic model, which is consistent with the results from our direct σ/m = 0.1 cm2 g−1 simulations at the radii where we can trust our simulations. This deserves further study both in terms of simulations with SIDM cross-section values smaller than 1 cm2 g−1 and more detailed comparisons to observations. Our current look at the global data does not suggest a need for a velocity-dependent cross-section as has been previously suggested. In a companion paper (Peter et al. 2012), we show that these SIDM models are also consistent with observations of halo shapes.
8.3 Observed substructure versus SIDM
In Fig. 8, we show that the number of subhaloes for σ/m = 1 cm2 g−1 is not significantly different from CDM predictions, especially in galaxy-scale haloes. This is interesting because it means that SIDM fails to deliver on one of the original motivations for considering this model of dark matter. Recall that Spergel & Steinhardt (2000) originally promoted SIDM as a solution to the missing satellites problem (Klypin et al. 1999; Moore et al. 1999), stating that many subhaloes would be evaporated by interactions with the background halo. Given the new discoveries of ultrafaint galaxies around the Milky Way and the high likelihood of many more discoveries from surveys like LSST (Bullock et al. 2010; Willman 2010), a significant reduction in substructure counts may very well be a negative characteristic of any non-CDM model (Tollerud et al. 2008).
However, in Milky Way mass haloes, SIDM with σ/m = 1 cm2 g−1 will yield a significant probability for subhalo particle scattering only for systems that pass within ∼10 kpc of the host halo centre. Thus, for this cross-section, we can form interesting-sized cores but largely leave the subhalo mass function unaffected in Milky Way mass haloes. For smaller cross-sections, the differences between SIDM and CDM subhalo mass functions will be even smaller. We note that we are not the first to find that SIDM can form cores but not solve the missing satellites problem; it was first discussed in D’Onghia & Burkert (2003).
This finding is also interesting in the context of other alternatives to CDM. Warm dark matter (WDM) models, for which the outstanding difference from CDM is that dark-matter particles have high speeds at matter–radiation equality and a related free-streaming cut-off in the matter power spectrum, predict a suppression in the halo (and subhalo) mass function at small scales. Otherwise, the abundance and structure of haloes and subhaloes are nearly indistinguishable from CDM (Villaescusa-Navarro & Dalal 2011; Maccio’ et al. 2012). WDM haloes may be less concentrated than CDM haloes on scales not much larger than the free-streaming scale, but are still cusped. They are only significantly cored right at the free-streaming scale, at which the halo and subhalo abundance is highly suppressed. Thus, each of the two leading modifications to CDM can solve only one of the two historical motivations for looking beyond the CDM paradigm.
The lack of subhalo suppression for σ/m ≲ 1 cm2 g−1 has implications for another of the SIDM halo constraints from a decade ago. Gnedin & Ostriker (2001) set a constraint excluding the range of 0.3 < σ/m < 104 cm2 g−1 based on the Fundamental Plane of elliptical galaxies. The argument rests on the observation that there are no significant differences in the Fundamental Plane of field and cluster ellipticals (e.g. Kochanek et al. 2000; Bernardi et al. 2003; La Barbera et al. 2010). Elliptical galaxies have a significant amount of dark matter within their half-light radii, with more massive ellipticals having larger mass-to-light ratios, either caused by varying stellar mass-to-light ratios or varying dark-matter content (Padmanabhan et al. 2004; Tollerud et al. 2011; Conroy & van Dokkum 2012). Gnedin & Ostriker (2001) argue that elliptical galaxies falling into cluster-mass haloes should have dark matter evaporated from their centres if σ/m ≠ 0, which would cause the stars in the elliptical galaxy to adiabatically expand and hence move the galaxy off the Fundamental Plane.
However, in our simulations, we find that few subhaloes are fully evaporated, and that the subhalo Vmax function is not greatly different for σ/m = 1 cm2 g−1 from CDM. In addition, our analytic arguments show that the trend with (host) halo mass for the evaporation of subhaloes at fixed r/rs is mild. This suggests that the Gnedin & Ostriker (2001) constraints are overly conservative even at the σ/m ≃ 1 cm2 g−1 level. The main caveats are that the suppression of the subhalo Vmax function is higher in more massive clusters and that the suppression is highest at the centre of the cluster halo. It would also be interesting to see if there are any differences in the Fundamental Plane as a function of projected distance in the cluster, both observationally and in simulations. For all of these reasons, it would be worthwhile to perform simulations of elliptical galaxies in clusters with SIDM and explore the Fundamental Plane constraints in more depth.
To summarize, although we have not fully resolved the cores of σ/m = 0.1 cm2 g−1 SIDM haloes, the intuition gleaned from our analytic model (tested against the SIDM1 results) and our moderately resolved simulation results suggest that σ/m = 0.1 cm2 g−1 is an excellent fit to the data across the range of halo masses from dwarf satellites of the Milky Way to clusters of galaxies. Values of cross-section over dark-matter particle mass in this range are fully consistent with the published Bullet cluster constraints (cf. Section 1), measurements of dark-matter density on small scales and subhalo survival requirements. In a companion paper (Peter et al. 2012), we show that this model is also consistent with halo shape estimates. It is therefore important to simulate galaxy and cluster haloes with cross-sections in the 0.1 cm2 g−1 range.
9 SUMMARY AND CONCLUSIONS
We have presented a new algorithm to include elastic self-scattering of dark-matter particles in N-body codes and used it to study the structure of SIDM haloes simulated in a full cosmological context. Our suite of simulations (summarized in Table 1) relies on identical initial conditions to explore SIDM models with velocity-independent cross-sections, σ/m = 1 and 0.1 cm2 g−1 as well as a comparison set of standard CDM simulations (with σ/m = 0).
Our primary conclusion is that while SIDM looks identical to CDM on large scales, SIDM haloes have constant-density cores, with core radii that scale in proportion to the standard CDM scale radius (rcore ≃ ϵ rs). The relative size of the core increases with increasing cross-section (ϵ ≃ 0.7 for σ/m = 1 and ϵ ∼ 0.2 for σ/m = 0.1 cm2 g−1). Correspondingly, at fixed halo mass, core densities decrease with increasing SIDM cross-section. For both core radii and core densities, there is significant scatter about the scaling with Vmax of the halo. The scaling relationship is strong enough that measurements of dark-matter densities in the cores of dark-matter-dominated galaxies and large galaxy clusters likely provide the most robust constraints on the dark-matter cross-section at this time. In a companion paper (Peter et al. 2012) we demonstrate, in contrast to previous claims, that SIDM constraints from halo shape measurements may be less restrictive than (or at least similar to those from) measurements of absolute core densities alone.
Based on our simulation results, we conclude that the dark-matter self-scattering cross-section must be smaller than 1 cm2 g−1 in order to avoid underpredicting the observed core densities in galaxy clusters, LSBs and dSph galaxies. However, an SIDM model with a velocity-independent cross-section of about σ/m = 0.1 cm2 g−1 appears capable of reproducing reported core sizes and central densities of dwarfs, LSBs and galaxy clusters. Higher resolution simulations with better statistics will be needed to confirm this expectation.
An accounting of our results is as follows.
Outside of the central regions of dark-matter haloes (r ≳ 0.5Rvir), the large-scale properties of SIDM cosmological simulations are effectively identical to CDM simulations. This implies that all of the large-scale confirmations of the CDM theory apply to SIDM as well.
The subhalo Vmax function in SIDM with σ/m = 1 cm2 g−1 differs by less than ∼30 per cent compared to CDM across the mass range 5 × 1011-2 × 1014 M⊙ studied directly with our simulations. Differences in the Vmax function with respect to CDM are only apparent deep within the centres of large dark-matter haloes. Thus, although possible, it will be difficult to constrain SIDM models based on the effects of subhalo evaporation.
SIDM produces haloes with constant-density cores, with correspondingly lower central densities than CDM haloes of the same mass. For σ/m = 1 cm2 g−1, our simulated halo density structure is reasonably well characterized by a Burkert (1995) profile fit with a core size rb ≃ 0.7rs, where rs is the NFW scale radius of the same halo in the absence of self-interactions. Core densities tend to increase with decreasing halo mass (ρb∝M− 0.2vir) but demonstrate about a factor of ∼3 scatter at fixed mass (likely owing to the intrinsic scatter in dark-matter halo concentrations).
SIDM halo core sizes, central densities and associated scaling relations can be understood in the context of a simple analytic model. The model treats the SIDM halo as consisting of a core region, where self-interactions have redistributed kinetic energy to create an approximately isothermal cored density profile; and an outer region, where self-interactions are not effective. The transition between these regions is set by the strength of the self-interactions and this model allows us to make quantitative predictions for smaller cross-sections where the cores are not resolved by our simulations. Based on this model and a few of our best resolved simulated haloes, we find core sizes ∼0.1rs for σ/m = 0.1 cm2 g−1.
Halo core densities over the mass range from 1015 to 1010 M⊙ in SIDM with σ/m = 1 cm2 g−1 are too low (∼ 0.005–0.04 M⊙ pc− 3) to match the observed central densities in galaxy clusters (∼ 0.03 M⊙ pc− 3) and dSphs (∼ 0.1 M⊙ pc− 3).
Halo core central densities in SIDM with σ/m = 0.1 cm2 g−1 are in line with those observed from galaxy clusters to tiny dwarfs (0.02–0.5 M⊙ pc− 3) without the need for any velocity dependence. The densities are more consistent with observations than those predicted in dissipationless CDM simulations, which are generically too high. SIDM models with this cross-section over dark-matter particle mass value are consistent with Bullet cluster observations, subhalo survival requirements and, as we show in a companion paper (Peter et al. 2012), measurements of dark-matter halo shapes.
Future work is necessary to expand both the dynamic range of our simulations in halo mass and resolution as well as the dynamic range in cross-sections. These simulations are necessary in order to make detailed comparisons with observations given the exciting possibility that dark-matter self-interaction with σ/m in the ballpark of 0.1 cm2 g−1 could be an excellent fit to the central densities of haloes over 4–5 orders of magnitude in mass.
MR was supported by a CONACyT doctoral Fellowship and NASA grant NNX09AG01G. AHGP is supported by a Gary McCue Fellowship through the Center for Cosmology at UC Irvine, NASA Grant No. NNX09AD09G at UCI, National Science Foundation (NSF) grant 0855462 at UCI and the NSF under Grant No. NSF PHY11-25915 while visiting the Kavli Institute for Theoretical Physics. JSB was partially supported by the Miller Institute for Basic Research in Science during a Visiting Miller Professorship in the Department of Astronomy at the University of California Berkeley. JO was supported by a Fullbright-MICINN Post-doctoral Fellowship. MK is supported by NASA grant NNX09AD09G and NSF grant 0855462. This research was supported in part by the Perimeter Institute of Theoretical Physics during a visit by MK. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development and Innovation. The work of LAM was carried out at Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. LAM acknowledges NASA ATFP support. Simulations were performed in the Pleiades supercomputer of the NASA Advanced Supercomputing (NAS) Division, and the Kraken supercomputer of the National Institute for Computational Sciences (NICS) through an XSEDE allocation.
This equation applies only if hsi is the same for both particles. See Appendix A for the general form.
We define Mvir as |$M_{\mathrm{vir}}= \frac{4}{3} \pi \rho _{\rm b} \Delta _\mathrm{vir}(z) r_{\mathrm{vir}}^3$| and rvir as |$\tilde{\rho }(r_{\mathrm{vir}}) = \Delta _\mathrm{vir}(z) \rho _{\rm b}$|, where |$\tilde{\rho }(r_{\mathrm{vir}})$| denotes the overdensity within rvir, ρb is the background density and Δvir is the virial overdensity.
REFERENCES
APPENDIX A: DERIVATION OF THE HARD-SPHERE INTERACTION RATE IN N-BODY SIMULATIONS
APPENDIX B: TEST FOR THE SCATTERING KINEMATICS
We use the same set-up as described in Section 3 to test our implementation against the expected kinematics. For this we looked at the distributions of the post-scatter velocity magnitudes and directions for both the sphere and background particles. For the distributions of the velocity directions, we looked at the inclination and azimuthal angles of the post-scatter velocity vectors. The angles are defined such as the line of interaction is along the θ = 0 direction (i.e. the z-axis) and ϕ is the azimuthal angle about which the experiment is symmetric. The distributions resulting from our test simulation are compared to those obtained from the transformation of a uniform isotropic distribution in the centre-of-mass frame to the simulation/lab frame; the results are shown in Figs B1– B3. Fig. B1 shows that the distributions of the velocity magnitudes rise linearly from v = 0 to vs, followed by a sharp cut-off at v = vs, where vs is the relative speed between the sphere and the background. From conservation of energy, it is only possible to have particles with v > vs if they have interacted multiple times. Multiple interactions are not considered in our calculations of the theoretical distributions but they are allowed in our simulation; hence, in Fig. B1 one can observe a tail for velocities > vs on the distributions of both types of particle velocities, but not on the theoretical distribution. Looking at Fig. B2 one can see that most of the particles are scattered towards the θ = 45° directions, i.e. forming a 45° angle with respect to |$\boldsymbol {v}_{\rm s}$|. It is visible that the distributions resulting from the simulation in the left-hand panel of Fig. B2 are higher than expected for θ ≲ 20|${^{\circ}_{.}}$| This is again from the fact that multiple scatters are possible in the simulation and this is not included in the calculations of the theoretical histograms. We demonstrate that this is the case by showing in the right-hand panel of Fig. B2 the distributions obtained from the simulation when we exclude any particles with v > vs, excluding that way any particles that we know have interacted multiple times, as one can see doing this brings the distributions from the simulation to a better match with the theory. Fig. B3 shows the distributions of the velocities as a function of ϕ, these are flat as expected due to the symmetry of the experiment.