A Testable Conspiracy: Simulating Baryonic Effects on Self-interacting Dark Matter Halos

The dark matter (DM) paradigm has been tremendously successful in explaining the large-scale structure of our universe (see, e.g., Planck Collaboration et al. 2016; Rodríguez-Torres et al. 2016), though the precise nature of DM remains unknown. The simplest example of cold dark matter (CDM), consisting of a single, collisionless particle with negligible primordial thermal velocity dispersion, can match the large-scale data remarkably well. Alternatively, DM could be more complex, with nongravitational coupling to standard model particles (e.g., Boehm & Schaeffer 2005; Escudero et al. 2015) and/or new dark sector particles (e.g., Feng et al. 2010a; Khlopov et al. 2010; Lesgourgues et al. 2016); many models of this kind produce observable signatures in astronomical data sets (Mangano et al. 2006; Feng et al. 2009; Cyr-Racine et al. 2016). In this paper, we consider the possibility that DM has strong elastic self-scattering interactions and explore the implications of such interactions on the DM distributions within individual galaxies. We specifically focus on the back-reaction associated with galaxy formation.

Collisional or self-interacting dark matter (SIDM) was first explored in the context of galaxy formation by Spergel & Steinhardt (2000), who argued that SIDM models with the cross-section over mass $\sigma /m\sim 1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ should lead to observable constant-density cores in galaxies, which is in agreement with observations at that time. While early estimates suggested that SIDM models of this kind would significantly reduce substructure counts compared to CDM, more recent numerical investigations have shown that the substructure differences are minimal (Vogelsberger et al. 2012; Rocha et al. 2013). However, the original expectation that SIDM halos should have constant-density cores has been demonstrated robustly in cosmological simulations (Davé et al. 2001; Rocha et al. 2013; Zavala et al. 2013).

SIDM cores are generated by energy-exchange interactions, which heat the halo center until it becomes isothermal. The radial extent of this core is set by the requirement that a typical DM particle will experience at least one interaction per Hubble time (Rocha et al. 2013). This implies that larger SIDM cross-sections produce larger isothermal cores. If the cross-section is large enough, then the isothermal region can extend beyond the peak in the halo's velocity dispersion profile; in this case, energy-exchange interactions could extract heat from the core leading to core collapse, which increases the central density (Kochanek & White 2000; Balberg et al. 2002; Colín et al. 2002; Koda & Shapiro 2011; Vogelsberger et al. 2012). However, this effect is muted in cosmological simulations and Elbert et al. (2015) used DM (only) zoom cosmological simulations to show that core collapse behavior sets in only for very large cross-sections $\sigma /m\gtrsim 10\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$.

The tendency for SIDM models with $\sigma /m\lesssim 10\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ to produce constant-density cores with lower overall density is of special interest for comparisons to dwarf and low-surface brightness (LSB) galaxies. This is because many of these galaxies are observed to have cores on roughly the scales that are expected in SIDM (Flores & Primack 1994; Moore 1994; de Blok et al. 1996; Salucci & Burkert 2000; de Blok et al. 2001; Swaters et al. 2003; Gentile et al. 2004; Simon et al. 2005; Spekkens et al. 2005; de Blok et al. 2008; Kuzio de Naray et al. 2008; Donato et al. 2009; Oh et al. 2011; Adams et al. 2014), as opposed to the cusps predicted in dissipationless CDM simulations (Dubinski & Carlberg 1991; Navarro et al. 1997). SIDM cores also may provide a natural explanation for the unexpectedly low densities of local dwarf galaxies (Vogelsberger et al. 2012, 2014; Elbert et al. 2015), which is a problem known as "Too Big to Fail" (TBTF; Boylan-Kolchin et al. 2011, 2012; Ferrero et al. 2012; Garrison-Kimmel et al. 2014; Tollerud et al. 2014; Klypin et al. 2015; Papastergis et al. 2015). There are many in the galaxy formation community who believe these issues may be resolved by baryonic processes, such as supernova feedback (Navarro et al. 1996; Read & Gilmore 2005; Governato et al. 2012; Pontzen & Governato 2012; Di Cintio et al. 2014; Maxwell et al. 2015; Oñorbe et al. 2015; Dutton et al. 2016; Read et al. 2016; Katz et al. 2017), though not all authors necessarily agree (Peñarrubia et al. 2012; Garrison-Kimmel et al. 2013; Pace 2016). Tidal effects have been shown to solve TBTF in satellite galaxies (see e.g., Read et al. 2006; Zolotov et al. 2012; Arraki et al. 2014; Brooks & Zolotov 2014; Del Popolo et al. 2014), but the evidence for TBTF in the local field (Garrison-Kimmel et al. 2014; Kirby et al. 2014) necessitates another solution for these galaxies. This ongoing debate and the lack of DM detections in direct, indirect, and collider searches motivates a thorough exploration of the SIDM hypothesis.

The goal of this paper is to investigate the late-time effects of galaxy formation on SIDM halos, specifically the contraction of these halos due to the gravitational potential of the galaxy. To this end, we use a set of N-body simulations similar to those initially used to examine contraction in CDM halos. The work is organized as follows. In Section 2, we briefly describe the properties required of a viable SIDM model, and in Section 3, we sketch the physics of contraction of SIDM halos and motivate our work in this paper. We describe our simulations and analysis in Section 4. We present our results in Section 5, discussing our Milky Way analog halos in Section 5.1 and our elliptical and LSB simulations in Section 5.2, while in Section 5.3, we directly compare our simulations to the analytic model presented in Kaplinghat et al. (2016). Section 6 shows the results of our cluster simulations and compares these with the observations of Newman et al. (2013b). We summarize our results and conclude in Section 7.

Previous work has placed constraints on the SIDM cross-section over mass across a range of halo masses. Generally, $\sigma /m$ below $0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ have been found to be indistinguishable from CDM models (Rocha et al. 2013). In low-mass galaxies with maximum circular velocity of ${V}_{\max }\,\simeq 30\ \mathrm{km}\,{{\rm{s}}}^{-1},\sigma /m$, values ranging from 0.5 to 10 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ alleviate the core-cusp and TBTF (Vogelsberger et al. 2012; Zavala et al. 2013; Elbert et al. 2015; Fry et al. 2015). However, values significantly in excess of 1 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ may lead to the efficient tidal stripping of stars in the satellites of the Milky Way and Andromeda (Gnedin & Ostriker 2001; Peñarrubia et al. 2010; Dooley et al. 2016), providing a possible avenue for an upper limit on the cross-section in the future. Recent work by Kaplinghat et al. (2016) showed that SIDM models with cross-sections around 2 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ fit the rotation curves of the 12 analyzed dwarf and LSB galaxies well. This work used an analytic model built on arguments discussed previously (Rocha et al. 2013; Kaplinghat et al. 2014) and showed that the analytic model is a good match for the density profiles of halos in DM-only cosmological SIDM simulations. The summary of above constraints is that for collisional velocities of order 100 km s⁻¹ or smaller, a $\sigma /m$ values close to 1 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ is favored and is consistent with all existing constraints. A larger sample of rotation curves will reduce uncertainties in the determination of the cross-section on galactic velocity scales.

The DM velocities in galaxy clusters are an order of magnitude larger than in dwarf galaxies, and many techniques have been used to constrain the DM self-interaction cross-section at these velocities. Cluster mergers have been used by many studies to constrain SIDM (see e.g., Randall et al. 2008; Dawson et al. 2012; Kahlhoefer et al. 2014; Massey et al. 2015; Schaller et al. 2015; Kim et al. 2017; Robertson et al. 2017), with typical limits of $\sigma /m\leqslant 1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ on the self-interaction cross-section. However, recent work by Kim et al. (2017) showed that constraints based on the displacement of the stellar and DM centroids are overly stringent, weakening previous constraints. They find that the displacement of the brightest cluster galaxy (BCG) relative to the halo center may be a better observable, possibly allowing $\sigma /m$ values around $0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ to be tested. The Bullet cluster constraint based on the mass loss from the merging subcluster (Markevitch et al. 2004) also needs to be reevaluated using self-consistent SIDM merger simulations and taking into account cosmic variance in the initial conditions. In addition, we also need theoretical refinement to apply these constraints to velocity-dependent cross-sections. Cluster shapes have provided an orthogonal method of investigating self-interaction cross-section on these scales. Core formation in SIDM halos leads to more spherical inner density profiles, so measurements of the ellipticities of cluster halos have been used to constrain $\sigma /m$ values to below $0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ (Miralda-Escudé 2002). However, due to the large scatter in axis ratios, the ability of SIDM halos to retain some triaxiality and the observational methods used to constrain halo shapes, the ellipticity constraints are unlikely to be better than about $1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ for cluster velocities (see Peter et al. 2013 for a detailed discussion).

The most stringent constraint on cluster velocity scale arises from the fact that the measured DM density profiles are in substantial agreement with CDM outside the half-light radii of the BCGs. Within about 10–50 kpc (range of BCG half-light radii), however, the DM density profile is shallower than the CDM expectations (Newman et al. 2013b), as we discuss later. Kaplinghat et al. (2016) used these measurements to show that the preferred cross-section for relative velocities larger than 1000 km s⁻¹ is about 0.1 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, which is consistent with earlier results from Yoshida et al. (2000). If the inferred shallowness of the density profile is due to AGN feedback or some other baryonic process (e.g., Martizzi et al. 2013), then this value of 0.1 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ provides a stringent upper limit on the self-interaction cross-section for velocities in excess of 1000 km s⁻¹. This result and the large difference in DM velocities in dwarf galaxies and galaxy clusters demand that viable SIDM models must have a velocity-dependent self-interaction cross-section. The required velocity dependence—from 1 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ for velocities below about 100 km s⁻¹ to 0.1 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ for velocities above 1000 km s⁻¹—can be easily accommodated in a variety of particle physics models (Feng et al. 2010b; Loeb & Weiner 2011; Tulin et al. 2013; Boddy et al. 2014, 2016; Cline et al. 2014).

Our work on simulating SIDM halos including a baryonic component is important for two specific reasons. First, in galaxies dominated by baryons there is no systematic evidence for large cores or lowered DM density profiles (e.g., Cappellari et al. 2015). Second, in galaxy clusters the stars dominate the total mass budget within their half-light radii, yet the DM tends to be under-dense compared to predictions. There is a simple analytic model (which we discuss later) that can explain both these observations but it has not been tested against simulations including a stellar component. By testing the accuracy of the analytic model we are able to bolster the case for a velocity-dependent self-interaction cross-section. We also test the possibility of core collapse in systems that have extremely dense baryonic distributions.

In this work, our aim is to test and further elucidate the general physics of how SIDM halos become isothermal in the potential well of the baryons. Because of this, we have chosen to run idealized simulations with disks grown adiabatically. This is complementary to the approach of running full-fledged hydrodynamic simulations that include self-interactions between DM particles. Vogelsberger et al. (2014) and Fry et al. (2015) examine dwarf galaxies (${M}_{{\rm{v}}}\sim {10}^{10}\,{M}_{\odot }$) using fully self-consistent hydrodynamic simulations and find that observable cores are still formed in these dwarfs in SIDM. Fry et al. (2015) find that the cores in their simulated dwarfs are not substantially different from those formed purely via feedback in their CDM simulations. It is not clear how these results generalize to larger halo masses, where feedback is expected to be less important in driving core formation in CDM halos and halo contraction effects are expected to dominate (e.g., Di Cintio et al. 2014; Dutton et al. 2016; Fiacconi et al. 2016).

The fact that growing baryonic potentials can cause contraction of CDM halos was first investigated analytically by Blumenthal et al. (1986) and Ryden & Gunn (1987), who used an adiabatic invariant approach. They demonstrated that if a baryonic potential grew to dominate the central potential, then the entire dark halo would contract, increasing the central DM density by more than an order of magnitude in plausible cases. Other studies (e.g., Jesseit et al. 2002) ran numerical simulations of isolated DM halos with disk potentials and found the DM density in these halos reproduced analytical predictions. Gnedin et al. (2004) studied baryonic contraction in hydrodynamic cosmological simulations and found that both adiabatic model predictions and isolated simulations produced central densities that were roughly 50% too high in halos where baryons dominate and proposed an alternative model to encapsulate the adiabatic contraction effect.

Though previous work on halo contraction explored collisionless DM, the general results should apply to SIDM as well. The adiabatic contraction occurs at early times when the disk or bulge of a galaxy is forming. These timescales are much shorter than the typical timescale for self-interactions in galaxies, assuming a cross-section around 1 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$. At later times, the self-interaction process will allow the halo to become isothermal, and this process can make the DM density profile shallower, retain the steep density profiles created by adiabatic contraction, or steepen it further in the case of core collapse. As we show later, the outcome depends on the cross-section, the gravitational potential of the baryons and the outer (CDM) DM halo profile.

Feedback could change the density profile of DM for dwarf galaxies, but it is not expected to be important for the baryon-dominated systems we simulate, which are Milky Way size or larger (Di Cintio et al. 2014; Dutton et al. 2016). It is important to note that when the equilibration timescale is short, changes due to feedback need to happen very late in order to prevent the DM from moving toward the equilibrium solution.

Kaplinghat et al. (2014) discussed the response of SIDM halos to the formation of a stellar disk or bulge using analytic equilibrium models. They found that the resultant SIDM core radius and density should be linked closely to the underlying baryonic potential in systems where the baryons are important dynamically. In the limit where the stars dominate the gravitational potential, the DM density profile (in the region where it is isothermal) scales as $\exp (-{{\rm{\Phi }}}_{\star }(r)/{\sigma }_{v0}^{2})$ (Amorisco & Bertin 2010; Kaplinghat et al. 2014) where ${{\rm{\Phi }}}_{\star }$ is the gravitational potential of the stars and ${\sigma }_{v0}$ is the 1D central velocity dispersion of DM. Thus, stars and DM are tied together in terms of the density profile and the shape of the DM halo must follow the contours of the stellar gravitational potential. This result emphasizes the need to account for the baryons when exploring SIDM phenomenology in galaxies with significant gas or stellar components.

The goal of this paper is to numerically investigate the effects of baryonic contraction on SIDM halos using a set of isolated N-body DM simulations of Milky Way, elliptical, lower sideband (LSB), and cluster analog halos. These simulations are similar in spirit to those first used to test and confirm analytic contraction models in the context of CDM. All of our halos were simulated with fixed SIDM cross-sections of $\sigma /m=0.5\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, but we show below that our simulations reproduce the model of Kaplinghat et al. (2016), indicating we can use their model to extend our results to a wide range of cross-sections.

Because we are primarily interested in isolating the effects of baryonic contraction on SIDM halos, as opposed to modeling a fully realistic treatment of galaxy formation in an SIDM universe, we ignore the effects of halo growth and baryonic feedback in our simulations. To model halos after rapid early accretion has finished, we choose a spherical Navarro et al. (1997) density profile for our initial conditions because it describes the equilibrium state of DM halos well. We may then add the galaxy potential in order to simulate a system without violent starbursts or major mergers. These assumptions and approximations work well for galaxies with quiescent periods greater than SIDM thermalization times ($\lesssim 1$ Gyr).

Our code is a modified version of GADGET-2 (Springel 2005) that allows for the inclusion of hard-sphere scattering between DM particles (Rocha et al. 2013). The simulations consist of a series of 3–50 million particle DM halos, initialized as Navarro et al. (1997) profiles, and run in isolation with and without an analytic galaxy potential. The potentials are grown linearly in time from a mass of zero to a final mass ${M}_{\mathrm{gal}}$ in 1 Gyr at the start of our simulations. We also simulated our Milky Way halo forming the Fiducial Disk after the SIDM core stabilized and found no difference in the DM distribution. The particle initial conditions were generated by the public code SPHERIC,⁴ which was first introduced in Garrison-Kimmel et al. (2013). To increase our effective resolution we exponentially truncate the outer regions of our initial Navarro–Frank–White (NFW) halos. These truncation radii lie far outside the halo scale radius except in the case of the cluster, which we truncate exponentially at 500 kpc in order to resolve the central few kpc of the halo. We have compared our idealized, dark-matter-only halos to analogous cosmological halos presented in Rocha et al. (2013) and found their radial profiles to be in good agreement. This demonstrates that the rapid merging process that forms halos in cosmological simulations does not cause a significant deviation from NFW or isothermal profile.

Table 1 summarizes our simulations, which consist of four characteristic halo/galaxy mass⁵ combinations:

• 1.  
  Milky Way: ${M}_{{\rm{v}}}={10}^{12}\,{M}_{\odot }$ with disks ${M}_{\mathrm{gal}}=5\,\times {10}^{10}\,{M}_{\odot };$
• 2.  
  LSB: ${M}_{{\rm{v}}}=2\times {10}^{11}\,{M}_{\odot }$ with disk ${M}_{\mathrm{gal}}=6\times {10}^{8}\,{M}_{\odot };$
• 3.  
  Elliptical: ${M}_{{\rm{v}}}=1.8\times {10}^{12}\,{M}_{\odot }$ with ${M}_{\mathrm{gal}}=6.2\,\times {10}^{10}\,{M}_{\odot };$ and
• 4.  
  Cluster: ${M}_{{\rm{v}}}=1\times {10}^{15}\,{M}_{\odot }$ with ${M}_{\mathrm{gal}}=2.1\times {10}^{12}\,{M}_{\odot }$.

Table 1.  Summary of Simulated Halos

Name	${M}_{{\rm{v}}}$	${c}_{{\rm{v}}}$	Np	rcut		Convergence Radius	Potential Shape	${M}_{\mathrm{gal}}$	$a$	$b$	h
	(${10}^{12}\,{M}_{\odot }$)		(106)	(kpc)	(kpc)	(kpc)		(${10}^{10}\,{M}_{\odot }$)	(kpc)	(kpc)	(kpc)
MW	1.0	13	3	230	0.4	0.83	MN Disk	5.0	1.5, 3.0, 6.0	0.3	⋯
LSB	0.2	11.8	10	170	0.19	0.30	MN Disk	0.06	2.2	0.2	⋯
Elliptical	1.8	9.7	10	300	0.37	1.0	Hernquist Sphere	6.2	⋯	⋯	3.0
Cluster	103	3	50	500	2	3.4	Hernquist Sphere	210	⋯	⋯	28.5

Note. The first five columns list identifying names and general simulation properties: the halo mass, NFW concentration, particle number, exponential cut-off radius, and force softening. We define ${M}_{{\rm{v}}}$ following Bryan & Norman (1998), with a virial over-density of ${{\rm{\Delta }}}_{{\rm{v}}}=97$ with respect to the critical density. The sixth column lists the convergence radius for the SIDM runs, which we set to 0.6 times the Power et al. (2003) radius for CDM, as found in Elbert et al. (2015). The last four columns summarize the properties of the galaxy potentials grown in each case, where ${M}_{\mathrm{gal}}$ is the final galaxy mass and the other shape parameters are defined in Equations (1) and (2). Note that there are three separate disks of varying scale length for the Milky Way runs. We refer to these in the text and figures as "compact," "fiducial," and "extended."

Download table as:  ASCIITypeset image

Each of these simulations are run with CDM and SIDM, and also with and without the galaxy potentials for comparison. The Milky Way halo mass simulations include three separate galaxy disk potential runs, each of which has a fixed galaxy mass but a variable scale length (see below). We present 22 simulations in all. All SIDM halos were run with $\sigma /m=0.5\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ and used a self-interaction smoothing factor of 25% of the force softening length (Rocha et al. 2013). We also resimulated our cluster halo with $\sigma /m=0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$. The force softening and convergence radius (Power et al. 2003; Elbert et al. 2015) of each run are indicated in columns five and six.

For the disk potentials in the LSB and Milky Way runs, we adopt the form of Miyamoto & Nagai (1975, MN hereafter):

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{MN}}(R,z)=\displaystyle \frac{-{{GM}}_{\mathrm{gal}}}{\sqrt{{R}^{2}+{(a+\sqrt{{b}^{2}+{z}^{2}})}^{2}}},\end{eqnarray} \tag{ 1 }$$

where a defines a scale length, and b sets a scale height.

Each Milky Way simulation has a fixed scale height of $b=0.3\,\mathrm{kpc}$. We explore three scale lengths: $a=1.5,3.0$, and 6.0 kpc, which we refer to as "Compact Disk," "Fiducial Disk," and "Extended Disk," respectively. These values⁶ roughly span the lower two-sigma to upper one-sigma of disk sizes for ${M}_{\mathrm{gal}}=5\times {10}^{10}\ {M}_{\odot }$ galaxies (Reyes et al. 2011). In particular, the compact disk is extremely dense and was chosen to investigate whether or not core collapse occurs. Our LSB disk mimics a typical LSB from Kuzio de Naray et al. (2008), with ${M}_{\mathrm{gal}}=6.3\times {10}^{8}\ {M}_{\odot },a=2.2\,\mathrm{kpc}$, and $b=0.2\,\mathrm{kpc}$.

We use spherical Hernquist (1990) distributions for the elliptical and cluster galaxy runs:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{H}}}(r)=\displaystyle \frac{-{{GM}}_{\mathrm{gal}}}{(r+h)}.\end{eqnarray} \tag{ 2 }$$

For the elliptical, we adopt $h=3.0\,\mathrm{kpc}$ and ${M}_{\mathrm{gal}}=6.2\,\times {10}^{10}\,{M}_{\odot }$, which are motivated by matching the median of bin 28 in Graves et al. (2009). We relate the typical effective radius for galaxies of this size by demanding that the 3D half-light radii are equal: $h={R}_{{\rm{e}}}/(1+\sqrt{2})$. For the central cluster galaxy, we match the results for A2667, as quoted in Newman et al. (2013a), by fitting a Hernquist profile with the same half-light radius to the best-fit dual pseudo-isothermal elliptical (dPIE) profile, which yields ${M}_{\mathrm{gal}}=2.1\times {10}^{12}\,{M}_{\odot }$ and $h\,=28.5\,\mathrm{kpc}$.

Because we do not directly model rapid galaxy formation and feedback in our simulations, we restrict our analysis to the late-time results of contraction, after most of these effects have subsided in massive galaxies. Each simulation was analyzed after reaching dynamical equilibrium, such that the density profile was no longer evolving. This occurred within ∼5 Gyr for all cases except the "compact" SIDM Milky Way, which underwent core collapse and showed a slowly increasing core density until we stopped the simulation after 10 Gyr.

5.1. Milky Way Halos

The setup for the Milky Way analog simulations is shown in the left column of Figure 1. The black dotted line shows the circular velocity curve (${V}_{{\rm{c}}}(r)=\sqrt{{GM}(\lt r)/r}$) for our CDM-only run (equivalent to the initial conditions) and the solid black line shows ${V}_{{\rm{c}}}(r)$ for SIDM-only run, which is noticeably less dense in the center owing to SIDM core formation (e.g., Rocha et al. 2013). The colored lines show the implied in-plane circular velocities ($\sqrt{r\tfrac{d{\rm{\Phi }}}{{dr}}}{| }_{z=0}$) for the imposed Extended, Fiducial, and Compact Disk potentials. These show that the disk potential is dominant in our compact and fiducial runs. Our goal is to explore the halo back-reaction to the growth of each of these components.

Zoom In Zoom Out Reset image size

The middle column of Figure 1 shows the DM density profiles for all simulations. The solid lines correspond to the SIDM model and dotted lines correspond to the CDM model. The gray curves are the dark-matter-only runs, while the colored lines show what happens after the potentials are grown. For reference, the bottom panel plots the ratio of DM density in SIDM to CDM as a function of radius for each set of runs.

As expected, both CDM and SIDM halos contract in response to galaxy potential growth.⁷ The relative differences are enlightening. While the dark-matter-only simulations are quite distinct between the cuspy CDM run and the cored SIDM within ∼3 kpc, the DM profiles in the Fiducial Disk runs are almost identical down to the resolving limit. Specifically, the SIDM halo has responded more to the imposed potential than the CDM halo and this has driven the two profiles to a very similar end state. We only begin to see the formation of an SIDM core within ∼1 kpc, which is similar to the Milky Way core size measured in Portail et al. (2017). The Extended Disk runs, which impose a less severe potential, have maintained something closer to the original differences, with SIDM beginning to roll off toward a core within ∼3 kpc, but the differences between the SIDM and CDM are less severe than in the DM-only case (which disagree at 5 kpc). Finally, the Compact Disk has produced a dramatic change: the SIDM halo is now more dense than CDM at small radii, with a very cuspy distribution $\rho \sim {r}^{-2.5}$ at $r\sim 2\,\mathrm{kpc}$. This is a result of core collapse: the compact disk potential has heated the DM to such an extent that it is now hotter in the core than in the outer part. The SIDM particles are conducting heat outwards, resulting in a loss of core pressure and subsequent mass inflow.

The SIDM phenomenology is clarified in the right panel of Figure 1, which shows the velocity dispersion profile of the DM in each run. The two SIDM simulations with clear constant-density core behavior (DM-only and Extended Disk) are seen to have well-established isothermal velocity distributions at small radii. In these cases, the SIDM halos are hotter in their cores and colder in their outer regions than their CDM counterparts. This is exactly the situation that leads to heat transfer from the outside in. The same effect is seen, though much more mildly, in the Fiducial Disk case. In the Compact Disk runs, even the CDM halo is hotter in the core than in the outer part. Such a declining velocity dispersion profile is subject to outward heat flow in the SIDM simulation, and this drives core collapse.

We show in Section 5.3 that the resultant density profiles are well explained by the analytical predictions presented in Kaplinghat et al. (2014, 2016), with the exception of the Compact Disk. However this is to be expected, given the gravothermal core collapse occurring in the Compact Disk, which violates the assumption of isothermality the model is based on.

Figure 2 displays visualizations of the DM in three pairs of our simulations, with CDM runs shown along the top and SIDM along the bottom. The left row shows the dark-matter-only versions of each simulation. The middle and right rows show resultant DM distributions after the growth of our Fiducial and Compact disks, oriented such that the disks are seen edge-on. These results emphasize that the shapes of the SIDM halos have been altered substantially by the formation of the disk, mirroring the baryonic potential much more closely than the CDM cases within 10 kpc.

Zoom In Zoom Out Reset image size

5.2. LSB and Elliptical Halos

Figure 3 shows the density profiles and initial rotation curves of our elliptical and LSB halos. In the top left panel, we see that the elliptical halo is baryon-dominated within ∼5 kpc, and consequently the density profiles (bottom left panel) show significant contraction. Again we see minimal difference between the final SIDM and CDM halos, with no SIDM core resolved. This is not surprising, given our results from 5.1; denser galaxies will have greater impacts on their host halos and inhibit SIDM core formation. We also display the total (DM + baryons) density profiles in the bottom left panel of Figure 3, as well as an ${r}^{-2}$ power-law for comparison. As a result of contraction, both simulations have power-law slopes of $\alpha =-2$ around 3 kpc, and slightly steeper outside this region. This places our simulations in agreement with observations of elliptical galaxies (Gavazzi et al. 2007; Auger et al. 2010). In short, because they are centrally baryon-dominated, predictions for SIDM halos are largely the same as the CDM case for elliptical galaxies.

Zoom In Zoom Out Reset image size

Our LSB-analog simulations, however, do exhibit differences between CDM and SIDM runs. The right panel of Figure 3 demonstrates that the LSB disk has very little effect on our DM halos; the CDM halo undergoes barely any contraction, and the SIDM halo still displays a central core. This makes sense in light of our previous results that massive and centrally concentrated baryon densities generate the largest impacts on both SIDM and CDM halos. The diffuse nature of low-surface-brightness galaxies implies that their baryons have little effect on the host DM halos. Thus, these systems are the best laboratories to investigate DM self-interaction possibilities.

5.3. Analytic Model

In this section, we compare the results of our simulations to the analytic model presented in Kaplinghat et al. (2016). In this model, the DM is assumed to behave as an isothermal gas within a radius r₁, defined as the radius where particles interact at least once in the age of the system: ${\rm{\Gamma }}({r}_{1})\,{t}_{\mathrm{age}}=1$ or

$$\begin{eqnarray}&&\rho ({r}_{1})\ (\sigma /m)\ (4/\sqrt{\pi }){\sigma }_{v0}\,{t}_{\mathrm{age}}=1.\end{eqnarray} \tag{ 3 }$$

Here, $\sigma /m$ is the SIDM cross-section per unit mass, ${\sigma }_{v0}$ is the radial velocity dispersion in the core and the factor of $4/\sqrt{\pi }$ accounts for the average relative velocity for Maxwellian distribution. t_age is the time period over which the self-interactions have been operating. We take this to be the time after the stellar component is fully in place to be consistent with the boundary conditions imposed.

The analytic model we use is spherically symmetric. To apply it to our simulated galaxies, we follow the procedure of Kaplinghat et al. (2014) and construct a new spherical mass profile for the baryons ${M}_{\mathrm{baryon}}(r)$ by including all the mass in all the stars within a sphere of radius r. Inside r₁, the SIDM density is set by hydrostatic equilibrium, giving

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\mathrm{ln}{\rho }_{\mathrm{DM}}(r)=-\displaystyle \frac{4\pi }{{\sigma }_{v0}^{2}}G[{\rho }_{\mathrm{DM}}(r)+{\rho }_{\mathrm{baryon}}(r)],\end{eqnarray} \tag{ 4 }$$

where ${\rho }_{\mathrm{baryon}}$ is the density profile corresponding to the mass profile M_baryon. As boundary conditions at r₁, the isothermal mass and density profiles (from solving the hydrostatic equation above) are required to match the CDM halo profile after adiabatic contraction (essentially the z = 0 CDM profile in our simulations).

In all cases, we find two islands of solutions in the central density—dispersion plane, for the same boundary condition Kaplinghat et al. (2014). We choose the lower density island because it matches the simulated profiles well. Typically, the 1d central dispersion ${\sigma }_{v0}$ values picked out by the two solutions are very close. The value of ${\sigma }_{v0}$ value is set by the matching, and it recovers the correct inner dispersion seen in the simulations.

Figure 4 shows circular velocity profiles for our contracted SIDM halos along with the analytic predictions for rotation curves. Here the cyan bands indicate the range of solutions which match $\rho ({r}_{1})$ and $M({r}_{1})$ within 5%. We chose 5% for two reasons. One, it shows the sensitivity of the inner density profile of SIDM profiles to the matching (boundary) conditions. Two, the simulated SIDM-only and CDM-only halos are not identical and allowing for a small slop in the matching makes us less sensitive to their differences.

Zoom In Zoom Out Reset image size

The analytic solutions and the simulated halo profiles generally agree well, with the notable exception of the Compact Disk "Milky Way" case. The analytic mass profile for this case is systematically higher, even though the non-trivial shape of the mass profile (see Figure 1) is captured well by the analytic profile. The analytic profile is very sensitive to the value of ${\sigma }_{v0}$. For the Fiducial and Extended disks, the matching correctly infers the SIDM central velocity disperison $\sqrt{3}{\sigma }_{v0}$ to be 250 and 210 km s⁻¹ (see Figure 1) and the analytic solution does a good job of describing the simulated mass profiles. For the Compact disk case, the inferred velocity dispersion is close to 290 km s⁻¹. However, in the Compact disk simulation, the velocity dispersion rises this high only within about 2 kpc, while r₁ is about 14 kpc. As Figure 1 shows, the dispersion profile for the Compact disk case is increasing toward the center due to core collapse (see Section 5.1), so the mismatch between the simulated profile and the analytic model is not surprising.

Given the overall consistency of the analytic model (equilibrium) predictions with our simulated halo profiles when core collapse is not seen, we are confident that the Kaplinghat et al. (2016) model predicts the SIDM halos in realistic galaxies and clusters well.

We also simulate an analog of Abell cluster 2667. Our setup is described in the bottom row of Table 1. The virial mass and concentration of the initial NFW halo match the observations of Newman et al. (2013a), however they model the density of the BCG as a dPIE profile:

$$\begin{eqnarray*}&&\rho =\displaystyle \frac{{\rho }_{0}}{1+{r}^{2}/({r}_{\mathrm{core}}^{2})(1+{r}^{2}/{r}_{\mathrm{cut}}^{2})}\end{eqnarray*}$$

In order to model this in our simulations we use a Hernquist (1990) sphere with the same half-mass–radius as the measured profile, and use a least-squares fitting method to obtain a best-fit total mass. This is necessary because dPIE profiles exhibit central cores, while Hernquist profiles have central $1/r$ cusps, so a Hernquist distribution with the same mass and half-mass–radius as a dPIE distribution will be less dense at all radii outside of the core radius and over dense at inner regions, significantly impacting the contraction effect we wish to investigate. The fact that we are using a cuspy profile in our simulations may exaggerate any contraction that occurs in our halos, but both the observed dPIE and our Hernquist profiles are dominant over the DM halo in the central regions, so contraction will occur in either scenario. We grow the central potential over the course of 2 Gyr, and allow the halo to evolve for 3 Gyr after this, running for 5 Gyr total.

Figure 5 shows the density profiles of our simulations; uncontracted densities are plotted as black dashed lines, while the contracted halos are solid black. We also include a vertical arrow to indicate the smallest converged radius of our simulations, as defined in Elbert et al. (2015). We plot the inferred cored NFW DM density profile measured by Newman et al. (2013b) as a cyan band, though we note that the generalized NFW fits Newman et al. (2013b) reported are not significantly different. Finally, the effective density of the BCG Hernquist profile is shown in magenta. The leftmost figure displays our CDM simulations. Even before contraction, the initial NFW profile is too centrally dense to match observations by a factor of ∼1.5 at the limits of our resolution. Adiabatic contraction increases the density by another factor of 1.5, further worsening this discrepancy. If instead we assume the DM is self-interacting with a cross-section of $0.5\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, then we obtain the densities in the center plot. In this case, we see that such a cross-section results in a halo that is under-dense by more than a factor of 3 in its center, and still below the observed limits at 100 kpc in our SIDM-only simulation. While baryonic contraction does increase the DM density, it is not nearly strong enough to completely alleviate the problem. Indeed, over most of our resolved region the contracted halo is only 20%–30% higher than the SIDM-only simulation.

Zoom In Zoom Out Reset image size

This result may seem surprising at first given that the baryons are important in the center in terms of their dynamical mass. However, the key point is that the DM halo is so massive that its velocity dispersion dwarfs the stellar potential. In this respect, the contraction of the SIDM halo is different from the adiabatic contraction of CDM halos; applying the Blumenthal et al. (1986) adiabatic contraction formula would result in a halo profile much denser than the simulated result. The isothermal equilibrium solution, on the other hand, is a good match to the simulated halo profile.

We also simulate Abell 2667 with a cross-section of $0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1};$ our results are plotted in the right panel of Figure 5. In the SIDM-only case, the density is near the lower limits of the Newman et al. (2013b) data. After contraction has been accounted for, however, $0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ is consistent with observations, bordering the upper limits of the Newman et al. (2013b) measurement, and we estimate that a cross-section of $0.2\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ will border the lower limits, though a central dPIE galaxy would cause less contraction, making a $0.2\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ cross-section less viable. Our estimates are in agreement with the results of Kaplinghat et al. (2016), who find cross-sections of $0.1-0.2\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ (assuming ${t}_{\mathrm{age}}=5$ Gyr) by fitting to the Newman et al. (2013b) data using the analytic model (described previously). This implies that an SIDM model with a cross-section that falls from about $1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ on dwarf galaxy scales to about $0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ on cluster scales can resolve the small-scale puzzles (Zavala et al. 2013; Vogelsberger et al. 2014; Elbert et al. 2015; Kaplinghat et al. 2016), while also matching density profile constraints in clusters.

The Hernquist profile was chosen to have the same half-light radius; however, in detail there are differences especially in the inner parts. It is reasonable to expect that these differences will be reflected in the SIDM density profile for cross-sections around $0.1\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$. To test for the effect of the stellar profile, we used the analytic model to predict the SIDM profile using the dPIE stellar profile rather than Hernquist. As expected, we found no significant differences when the cross-section was set to $0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$. For $0.1\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, the dPIE profile showed more contraction by about 25% at 1 kpc. These differences are smaller than those when the cross-section is increased by a factor of 2, which provides a guideline to the effect this systematic has on the inferred cross-section. Regardless of the profile shape, a cross-section as large as $0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ cannot fit the inferred density profile for Abell 2667 by Newman et al. (2013b).

In this work, we have investigated the combined effects of baryonic gravitational potentials and DM self-interactions on DM halos using idealized simulations of DM halos with galactic potentials. By simulating halos of various sizes with many different potentials we have found the following:

• 1.  
  SIDM halo shapes are not inherently more resilient to effects from baryons than their CDM counterparts. For a Milky Way halo hosting a Milky Way analog disk, the SIDM halo is more compact along the disk axis than its CDM equivalent in agreement with the prediction of Kaplinghat et al. (2014). For an elliptical galaxy, whose stellar potential is markedly more spherical, we expect the SIDM halo to be correspondingly more spherical.
• 2.  
  Halos that host substantial baryonic populations display few differences in spherically averaged density profiles between CDM and SIDM models on observable scales. Even extended baryon populations can contract halos with respect to SIDM-only simulations, though these systems retain potentially observable constant-density cores and are less dense than CDM. In extreme cases, we find that potentials from dense baryonic structure can cause SIDM halos to core collapse and become denser than their CDM counterparts.
• 3.  
  Halos that host concentrated stellar populations may display few differences in spherically averaged density profiles between CDM and SIDM models on observable scales. Even extended stellar populations can contract halos with respect to SIDM-only simulations, though these systems retain potentially observable constant-density cores and are less dense than CDM. In extreme cases, we find that potentials from dense baryonic structure can cause SIDM halos to core collapse and become denser than their CDM counterparts.
• 4.  
  The densities of our contracted halos are in good agreement with the analytical predictions in Kaplinghat et al. (2014, 2016), with the exception of the core-collapsing Compact Disk because it no longer obeys the isothermal assumption of the model. In particular, we show that the spherically averaged density profiles are well approximated by the simple model in Kaplinghat et al. (2016), which has an isothermal core and an undisturbed CDM outer profile.
• 5.  
  We simulated a cluster halo for 3 Gyr after the BCG was fully in place to test against the mass measurements for A2667 Newman et al. (2013b). Our simulated CDM halo was denser than the observed central profile for A2667. On the other hand, SIDM with $\sigma /m\simeq 0.5\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ was too low in density compared to the measurements. The choice of $\sigma /m\simeq 0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ was in much better agreement with the measured normalization (and inner slope) of the A2667 density profile. Larger values like $\sigma /m\simeq 0.5\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ are ruled out, even allowing for a factor of 2 uncertainty in the age of the halo. These conclusions are in substantial agreement with the detailed analysis of 7 clusters (Newman et al. 2013b) by Kaplinghat et al. (2016), which found a average value of $\sigma /m\simeq 0.1\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ on cluster velocity scales for an assumed age of 5 Gyr.

Based on these results, an ideal scale to investigate possible DM self-interactions appears to be the dwarf galaxy scale with halo masses ${10}^{10-11}\ {M}_{\odot }$ scale, as they will have the largest interaction cross-sections and the least contracted halos. However, these are precisely the halos expected to be most vulnerable to stellar feedback (Governato et al. 2012; Pontzen & Governato 2012; Di Cintio et al. 2014; Oñorbe et al. 2015). Ongoing work (Vogelsberger et al. 2014; Fry et al. 2015; V. Robles et al. 2018, in preparation) is investigating the effects of feedback on the SIDM halos and their results suggest that dwarfs with stellar masses ${M}_{\star }\lesssim {10}^{6}\,{M}_{\odot }$ will have density profiles indistinguishable from the predictions of the DM-only simulations. This suggests that the faintest dwarf spheroidals provide excellent laboratories constraining SIDM models.

For halo masses much larger than ${10}^{11}\ {M}_{\odot }$ that host a large stellar disk or bulge, as the inner halo becomes isothermal the SIDM halo retains the high densities created by adiabatic contraction following the formation of the disk. Thus, in Milky Way-sized halos the CDM and SIDM halos densities are very similar beyond about a kpc, in marked contrast to the dark-matter-only predictions. As predicted by Kaplinghat et al. (2014), the self-interactions also force the SIDM halo to be more compact along the stellar disk axis. We find that the SIDM halo in the inner region of Milky Way analogs is more compact along the galactic disk axis than the CDM halo. Thus, it may be possible to use the shape of the DM halo in the inner regions of large spiral galaxies to provide a sharp test of the SIDM paradigm.

The predictive cross-talk between the DM and baryons in the SIDM models leads to a large diversity of halo profiles and halo shapes. This cross-talk is purely gravitational and the result of the DM becoming isothermal in the inner parts of the halos. This is fully explained by simple equilibrium models. The prospects for using these concrete predictions of the SIDM paradigm to rule in or rule out SIDM in the near future are excellent.

The authors thank Michael Boylan-Kolchin and Hai-Bo Yu for valuable discussions. O.D.E., A.S.G., and J.S.B. were supported by the National Science Foundation (grants PHY-1520921, AST-1518291, and AST-1009973) and by NASA through HST (GO-13343) awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under NASA contract NAS5-26555. M.K. is supported by NSF grant PHY-1620638. Support for S.G.K. was provided by NASA through Einstein Postdoctoral Fellowship grant number PF5-160136 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060.

This work also made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), matplotlib (Hunter 2007), numpy (van der Walt et al. 2011), scipy (Jones et al. 2001), ipython (Perez & Granger 2007), and NASA's Astrophysics Data System. This work also made use of the EXtreme Science and Engineering Discovery Environment (XSEDE; Towns et al. 2014).