The central densities of Milky Way-mass galaxies in cold and self-interacting dark matter models
Omid Sameie, Michael Boylan-Kolchin, Robyn Sanderson, Drona Vargya, Philip Hopkins, Andrew Wetzel, James Bullock, Andrew Graus
MMNRAS , 1–10 (2020) Preprint 4 March 2021 Compiled using MNRAS L A TEX style file v3.0
The central densities of Milky Way-mass galaxies in cold andself-interacting dark matter models
Omid Sameie, ★ Michael Boylan-Kolchin, Robyn Sanderson, , Drona Vargya Philip F. Hopkins, Andrew Wetzel, James Bullock, Andrew Graus Department of Astronomy, The University of Texas Austin, 2515 Speedway, Stop C1400, Austin, TX 78712 USA Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104, USA Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010 USA TAPIR, California Institute of Technology, Pasadena, CA 95616, USA Department of Physics and Astronomy, University of California, Davis, CA 95616, USA Department of Physics and Astronomy, University of California, Irvine, CA 92697, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present a suite of baryonic cosmological zoom-in simulations of self-interacting darkmatter (SIDM) haloes within the “Feedback In Realistic Environment” (FIRE) project. Thethree simulated haloes have virial masses of ∼ M (cid:12) at 𝑧 =
0, and we study velocity-independent self-interaction cross sections of 1 and 10 cm g − . We study star formation ratesand the shape of dark matter density profiles of the parent haloes in both cold dark matter(CDM) and SIDM models. Galaxies formed in the SIDM haloes have higher star formationrates at 𝑧 ≤
1, resulting in more massive galaxies compared to the CDM simulations. Whileboth CDM and SIDM simulations show diverse shape of the dark matter density profiles, theSIDM haloes can reach higher and more steep central densities within few kpcs compared to theCDM haloes. We identify a correlation between the build-up of the stars within the half-massradii of the galaxies and the growth in the central dark matter densities. The thermalizationprocess in the SIDM haloes is enhanced in the presence of a dense stellar component. Hence,SIDM haloes with highly concentrated baryonic profiles are predicted to have higher centraldark matter densities than the CDM haloes. Overall, the SIDM haloes are more responsive tothe presence of a massive baryonic distribution than their CDM counterparts.
Key words: methods: numerical-galaxies: evolution-galaxies: formation-galaxies: structure-dark matter
The leading theory of structure formation assumes dark matter (DM)is cold and collisionless (Peebles 1982; Bond et al. 1982; Blumen-thal et al. 1984; Davis et al. 1985). Predictions of the cold dark matter(CDM) model are well tested on large scales (Seljak et al. 2006;Percival et al. 2007; Vogelsberger et al. 2014a,c; Planck Collabo-ration et al. 2018), yet there are strong hints on the galactic scaleswhich suggest the CDM model might fail to reproduce a handfulof observations (see Bullock & Boylan-Kolchin 2017, for a com-prehensive review). These small-scale challenges include the CDMhaloes from numerical simulations that are more massive than theobserved population of galaxies around the Milky Way (MW) andin the Local Group (Boylan-Kolchin et al. 2011, 2012; Garrison-Kimmel et al. 2014), more flattened DM density profiles in theobserved galaxies (core) vs. steeper profiles (cusp) from CDM-only ★ E-mail: [email protected] simulations (Moore 1994; Kuzio de Naray et al. 2008; Walker &Peñarrubia 2011; Oh et al. 2015), and the diversity in the shape ofthe observed galactic rotation curves (Oman et al. 2015).There are two general approaches to address the shortcom-ings of CDM-only simulations to accommodate these observations.Gravitational potential fluctuations caused by stellar feedback, eitherin the form of single feedback events (e.g., Navarro et al. 1996) orfrom repeated cycles of gas inflow and outflow (e.g., Mashchenkoet al. 2008; Pontzen & Governato 2012; Governato et al. 2012),can account for the mass deficit in the observed dwarf galaxies.The critical requirements for the feedback to create extended DMcores are bursty and extended star formation histories (Benítez-Llambay et al. 2019), and high density threshold for star formation 𝑛 sf ∼ − − (Hopkins et al. 2014, 2018b; Bose et al. 2019)in galaxies that form their stars relatively late in their cosmologicalevolution, 𝑧 ≤ 𝑀 ∗ / 𝑀 h ∼ − ; galaxies with lower baryon con- © a r X i v : . [ a s t r o - ph . GA ] M a r Sameie et al. version efficiencies retain cuspier profiles, while more complicatedeffects occur at higher baryon conversion efficiencies. (Di Cintioet al. 2014; Chan et al. 2015; Tollet et al. 2016; Fitts et al. 2017;Lazar et al. 2020).In the second approach, basic assumptions about the nature ofDM are revisited. For example, two-body non-gravitational interac-tions of DM particles have been introduced as a potential solution tothe small-scale issues (Spergel & Steinhardt 2000; Ahn & Shapiro2005; Ackerman et al. 2009; Arkani-Hamed et al. 2009; Feng et al.2009; Loeb & Weiner 2011; Tulin et al. 2013). N-body simulationsof self-interacting dark matter (SIDM) models have shown that self-interactions would create constant density cores with spherical haloshapes in the core region (Vogelsberger et al. 2012; Zavala et al.2013; Rocha et al. 2013; Peter et al. 2013). Core formation in SIDMmodels is prompted by redistribution of energy and momentum ofDM particles: dark matter scattering leads to a cored central densityand an isothermal central temperature profile. As a result, coredDM densities are prevalent in DM-dominated SIDM systems (seeTulin & Yu 2018, for a complete review).Recent simulations of SIDM haloes in the presence of centralbaryonic components reveal that baryons modify this appealinglysimple picture: once baryons dominate the central gravitational po-tential, a DM core will transform into a more cuspy profile (Elbertet al. 2018; Sameie et al. 2018; Robles et al. 2019; Despali et al.2019). The stellar concentration, along with the DM cross section,therefore plays a critical role in establishing the density profile ofSIDM haloes (Elbert et al. 2018; Sameie et al. 2018). An importantresult from these simulations is that SIDM can lead to a large diver-sity in the shape of rotation curves (Creasey et al. 2017; Kamadaet al. 2017; Ren et al. 2019), similar to what is observed in galaxyrotation curves (Lelli et al. 2016). The observed spread in the cen-tral DM densities of Milky Way (MW) satellites (Read et al. 2019;Kaplinghat et al. 2019) could also be triggered by SIDM (Koda &Shapiro 2011; Kahlhoefer et al. 2019; Sameie et al. 2020b,a; Correa2020): core collapse in tidally-evolving SIDM subhaloes can occurin much shorter time scales than in isolated SIDM haloes owingto mass loss in their outskirts. Nonetheless, most SIDM simula-tions have relied on some idealized realization of isolated galaxiesor MW satellites in which either the effects of star formation andfeedback are ignored or the parent haloes are treated as fixed po-tentials without considering their mass assembly and cosmologicalevolution.Despite intriguing results, there have been comparably few at-tempts to run baryonic cosmological simulations in the context ofSIDM models. There are simulations in the mass scales relevant todwarf galaxies (Vogelsberger et al. 2014b; Fry et al. 2015; Fitts et al.2019), galaxies with halo masses few times 10 to 10 M (cid:12) , repre-sentative of MW-mass and elliptical galaxies (Di Cintio et al. 2017;Despali et al. 2019), galaxy clusters (Robertson et al. 2018, 2019),and large-box cosmological simulations with spatially-uniform res-olution (Lovell et al. 2018). Most of these works either havemedium-range mass resolutions or simulate their haloes with onlyone SIDM cross section. Thus far, there has not been any attempt ata suite of fully baryonic cosmological SIDM simulations with mul-tiple cross sections and high mass resolution devoted to MW-massgalaxies and their satellites.This paper studies the interplay of baryonic physics and DMself-interactions at the mass scales relevant to the MW. By employ-ing the state-of-the-art FIRE-2 implementation of star formationand feedback, our simulations benefit from an explicit treatment ofthe interstellar medium (ISM) that resolves its multi-phase struc-ture. Results from CDM simulations with this implementation agree with a broad range of observations including galaxy morphologies,the internal structure of the ISM, star formation histories, and theobserved mass-size relation (Wetzel et al. 2016; Hopkins et al.2018a,b). Our simulation suite also has the high mass and spatialresolution ( 𝑚 dm ∼ M (cid:12) , 𝜖 minp ∼
20 pc) required to resolve theinner structure of the parent haloes.In this work, we focus on the stellar mass assembly and DMdensity profiles in our SIDM simulations and comparisons with theirCDM counterparts; a companion paper (Vargya et al., in prepara-tion) employs the same simulation suite to study the shapes of darkmatter haloes in SIDM and CDM, both with and without baryonicphysics. We analyze in detail the cosmological evolution of bothDM and stellar density profiles, and highlight the importance ofconcentration of stellar distribution on the diversity in DM densityprofiles. Our paper is organized as follows: in Section 2 we givea brief description of our simulations. Section 3 discusses overallproperties of CDM and SIDM haloes (sub-section 3.1), the densityprofiles of DM and stars (3.2), the cosmological evolution of DMand stellar distribution (3.3), and the role of baryon concentrationto develop diverse DM central densities in SIDM (3.4). In section 4we summarize our results.
We perform a suite of cosmological simulations with virial masses ∼ M (cid:12) at 𝑧 =
0. Our simulations are part of the “FeedbackIn Realistic Environments" project (FIRE, Hopkins et al. 2014,2018b) , and are run using GIZMO (Hopkins 2015) . The gravity issolved with an improved version of the Tree-PM solver from GAD-GET3 (Springel 2005) and the hydrodynamical equations are treatedvia the mesh-free Lagrangian-Godunov (MFM) method which pro-vides adaptive spatial resolution while maintaining conservation ofmass, energy, and momentum.The baryonic physics implementation in FIRE-2 includes cool-ing, star formation, stellar feedback including SNe Ia & II, multi-wavelength photo-heating, stellar winds, radiation pressure, andUV-background radiation, all taken from stellar evolutionary mod-els. Star formation happens in molecular gas clouds that are locallyself-gravitating, self-shielding, Jeans unstable, and above the den-sity threshold 𝑛 H > 𝑛 crit = − . A detailed description ofthe baryonic physics implementations can be found in Hopkins et al.(2018b). Our simulations reach mass resolutions 𝑚 b = . × M (cid:12) and 𝑚 dm = . × M (cid:12) for the baryons and DM, and minimumphysical spatial resolutions 𝜖 gas = 𝜖 ★ = 𝜖 dm =
20 pc.Initial conditions are generated using the zoom technique(Katz & White 1993; Oñorbe et al. 2014) at 𝑧 =
99, embedded withinperiodic cosmological boxes of length 𝐿 =
60 Mpc / ℎ using the codeMUSIC (Hahn & Abel 2011). We adopt the following cosmologicalparameters: Ω Λ = . Ω m = − Ω Λ = . Ω b = . 𝑛 s = . 𝜎 = . ℎ = . https://fire.northwestern.edu MNRAS000
60 Mpc / ℎ using the codeMUSIC (Hahn & Abel 2011). We adopt the following cosmologicalparameters: Ω Λ = . Ω m = − Ω Λ = . Ω b = . 𝑛 s = . 𝜎 = . ℎ = . https://fire.northwestern.edu MNRAS000 , 1–10 (2020)
IDM haloes On FIRE Table 1.
Galaxy properties. Columns describe: (1) Name: name and DM type of each simulation, (2) 𝑀 : halo mass identified as the mass enclosed by theradius where cumulative density is 200 times the mean density of the Universe ( 𝜌 crit ), (3) 𝑟 : the radius in which the mean density is equal to 200 times 𝜌 mean (4) 𝑉 max : maximum circular velocity, (5) 𝑉 : total circular velocity at 2 kpc, (6) Φ pot : total gravitational potential at the center of the halo (see section 3.4),(7) 𝑀 ★ : galaxy stellar mass, computed as the sum of stellar masses inside of 0 . ∗ 𝑟 , (8) 𝑟 / : 3D stellar half-mass radius, (9) 𝜌 DM ( . ) : DM densityprofile at 500 pc, in M (cid:12) kpc − (10) 𝛾 ( . - . ) : logarithmic DM density slope measured between 0 . † : We stop this simulation at 𝑧 = . 𝑧 . We have checked that the density and velocity dispersion profiles are very similar between 𝑧 = . 𝑧 = 𝑀 ( M (cid:12) ) 𝑟 ( kpc ) 𝑉 max ( km / s ) 𝑉 Φ pot ( km s − ) 𝑀 ★ ( M (cid:12) ) 𝑟 / ( kpc ) 𝜌 DM ( . ) 𝛾 ( - . ) m12i-CDM 0 . ×
197 187 160 − . × . × . × -0.49m12i-CDM-only 1 . ×
204 162 63 − . × - - 1 . × -0.94m12i-SIDM1 0 . ×
197 214 207 − . × . × . × -0.97m12i-SIDM10(z=0.1) † . ×
190 211 207 − . × . × . × -1.3m12f-CDM 1 . ×
228 222 221 − . × . × . × -0.91m12f-CDM-only 1 . ×
228 176 71 − . × - - 2 . × -1.09m12f-SIDM1 1 . ×
228 225 202 − . × . × . × -0.71m12f-SIDM10 1 . ×
222 220 218 − . × . × . × -0.92m12f-SIDM1-only 1 . ×
228 177 27 − . × - - 1 . × -0.03m12f-SIDM10-only 1 . ×
228 185 19 − . × - - 4 . × . ×
217 199 118 − . × . × . × -0.30m12m-CDM-only 1 . ×
217 171 85 − . × - - 3 . × -1.03m12m-SIDM1 1 . ×
217 221 126 − . × . × . × -0.16m12m-SIDM10 1 . ×
217 226 204 − . × . × . × -0.81 re-simulate them with the SIDM model with two constant DMcross sections 𝜎 / 𝑚 = g − (labelled as SIDM1 andSIDM10). In our standard simulation suite any thermal-to-kineticenergy conversion from the stellar mass-loss processes are ignoredfor the sub-resolution regions (Hopkins et al. 2018a,b). As we dis-cuss it further in section 3.1, this implementation could have sig-nificant effect on stellar masses and DM central densities (the samesimulated galaxies with maximum stellar mass-loss energy conver-sion could have 2-3x higher stellar masses and denser DM densityprofiles). We also perform two DM-only (DMO) SIDM simulationswith the aforementioned DM cross sections for the m12f halo for thepurpose of comparison. The SIDM module in GIZMO (introducedin Rocha et al. 2013; Peter et al. 2013) uses a Monte-Carlo approachto assign scattering probability for nearest neighbours of each DMparticle via a spline kernel (Monaghan & Lattanzio 1985), and thenisotropically assigns velocities to the scattered particles such thatenergy and momentum are conserved. It has been widely discussed in the literature that SIDM preservesthe predictions of CDM on the scales larger than galaxies (e.g. seeRocha et al. 2013; Vogelsberger et al. 2016; Sameie et al. 2019). Inthis section, we discuss the impact of SIDM on the overall featuresof galaxies embedded in their DM haloes. Table 1 summarizes haloand galaxy properties for the simulation suite; CDM and SIDM sim-ulations have no notable difference in their virial masses, 𝑀 ,while SIDM simulations have higher stellar masses at 𝑧 = 𝑟 / , are smaller in the CDM sim-ulations compared to the SIDM1 runs, while the simulations with We define the virial mass to be the mass enclosed by the sphere with anaverage density of 200 times the mean density of the universe. The radiusof this sphere is called virial radius, 𝑟 . 𝜎 / 𝑚 =
10 cm g − show signs of contraction in their stellar distri-bution compared to the SIDM1 suite. The maximum total circularvelocity, 𝑉 max , is comparable or higher in the SIDM haloes thantheir CDM halo counterparts. A related quantity often used to char-acterize halo rotation curves, the total circular velocity computedat 2 kpc ( 𝑉 ), is also similar or higher in the SIDM haloes thanthe CDM ones. Introducing baryonic processes in both CDM andSIDM simulations causes a systematic increase in both 𝑉 max and 𝑉 compared to their DMO version. Interestingly, SIDM haloesrespond more significantly to the presence of baryons with a propor-tionally greater increase in their 𝑉 , DMO / 𝑉 , Hydro . of ∼ . . - .
35 for the CDM haloes.We further explore the growth of stellar mass in CDM versusSIDM by comparing the star formation rates (SFR) of the simulatedhaloes in Fig. 1. We compute the SFR by taking all star particleswithin 0 . ∗ 𝑟 of the host halo’s center in the last snapshot, andplotting their formation times in 100-Myr intervals. All the simu-lated galaxies show bursty star formation with two distinct phases:an early phase when the SFR increases with time, followed by a con-stant and non-zero SFR that extends to the present time. However,after early, rapid phase of star formation, the SIDM models formstars at similar or higher rates compared to the CDM runs, leadingto galaxies with higher stellar masses at 𝑧 =
0. The relative increaseof the SFR in the SIDM models is different for each halo. SIDMruns in m12m haloes have a larger difference in their SFR relative toCDM, while the SFR profiles of the m12i and m12f-SIDM haloesare closer to the corresponding CDM run.We also perform three more simulations with the CDM modelfor our m12i,f,m haloes with a slightly different implementationof the stellar mass-loss processes (e.g. O/B and AGB-star photo-spheric outflows) to examine their impact on the SFR and densityprofiles of the simulated galaxies. In our default simulations, weignore any conversion of thermal to kinetic energy (i.e., any “ 𝑃 d 𝑉 work” done) during the Sedov-Taylor phase of the expansion fromstellar mass-loss processes for the unresolved regions (see Hop-kins et al. 2018a,b, for more details). Alternatively, stellar mass-loss processes in unresolved regions can be treated as a prolongedenergy-conserving phase during which substantial 𝑃 d 𝑉 work wasdone, converting almost all of the thermalized/shocked ejecta en- MNRAS , 1–10 (2020)
Sameie et al. S t a r f o r m a t i o n r a t e ( M y r ) m12i Cosmic time (Gyr) m12f
Red : SIDM10Orange : SIDM1Black : CDMRed : SIDM10Orange : SIDM1Black : CDM m12m
Red : SIDM10Orange : SIDM1Black : CDM
Redshift
Figure 1.
Star formation rate profiles for the CDM (black), SIDM1 (orange), and SIDM10 (red) models in m12i (left), m12f (middle), and m12m (right)galaxies. In order to compute the SFR profiles, we take all the star particles inside of 0 . ∗ 𝑟 of the host in the last snapshot, and we plot their star formationtimes in 100 Myr intervals. ergy into kinetic energy (momentum) on large scales. This is thedefault approach taken in FIRE-2 simulations. We do this by treat-ing each stellar mass-loss event (which injects some Δ 𝑀 ≡ (cid:164) 𝑀 ∗ Δ 𝑡 ,with initial free-streaming kinetic luminosity/energy Δ 𝐸 ≡ (cid:164) 𝐸 Δ 𝑡 )as a “mini-supernova” and applying the exact same treatment as wedo for SNe following Hopkins et al. (2018a).The practical effect of adopting this (standard FIRE-2) treat-ment of stellar mass-loss, given the various scalings for, e.g., thecooling radii of SNe, is that most of total energy injection by stel-lar mass-loss is converted into momentum/kinetic energy on re-solved scales (i.e. ∼ 𝛾 = . − . 𝜌 DM ( . ) = ( . − . ) × M (cid:12) kpc − and 𝑉 ∼
250 kpcfor all three haloes (see table 1 for comparisons with our standardCDM runs) relative to the default model. Our results highlight theimportance of sub-grid implementation of feedback models on cre-ating diverse DM central densities, and we will present a morecomprehensive analysis in a forthcoming paper.
Fig. 2 shows density profiles for DM (solid) and stars (dashed) for thesimulated haloes in CDM (black), SIDM1 (orange), and SIDM10 (red). For comparison, DM density profiles of the available DMOruns are shown by dot-dashed lines. The DMO versions of theCDM runs, predictably, show cuspy inner regions with the DMdensity slope 𝛾 ( 𝑟 ) ≡ dlog 𝜌 ( 𝑟 )/ dlog r ∼ − . 𝑟 = . - − . ≤ 𝛾 ≤ − . 𝛾 ∼ − . ≤ 𝛾 ≤ − .
28 for the full physics simulations (Table 1).The build-up of the stars in the baryonic versions of our simu-lations results in stellar components that dominate the gravitationalpotential in the inner regions. This causes all the haloes with baryonsto develop central densities that are up to 100x higher than theirDMO counterparts (except m12m with CDM), with SIDM haloesresponding more dramatically compared to their CDM counterparts.The impact of stellar mass assembly is opposite for the CDM andSIDM models: in CDM haloes, stellar feedback flattens the densityprofile while in SIDM haloes central density becomes cuspier. Evi-dently, feedback does little to flatten the central density of the SIDMhaloes; the combination of SIDM thermalization and baryonic con-traction tends to make the DM density profiles less cored than inthe DMO case. We note that the shape of DM density profiles isnot universal among all SIDM haloes, and different realizations ofMW-mass galaxies could have different central DM densities andslopes (Table 1). As we discuss in the next few sections, a combina- DMO runs are identical to the hydro version except there are no baryons,so the baryonic mass is added to the DM particle mass. In order to havea direct comparison to DM density profiles from the hydro versions, wemultiply the density by 1 − Ω b / Ω m ∼ .
83 to account for the absence of thebaryons. MNRAS000
83 to account for the absence of thebaryons. MNRAS000 , 1–10 (2020)
IDM haloes On FIRE D e n s i t y p r o f il e ( M k p c ) m12i Black : CDMOrange : SIDM1Red : SIDM10 m12f
StarsDMonlyDM m12m V c i r , t o t a l ( k m s ) Distance (kpc)
Figure 2. Top : Spherically-averaged density profiles for DM (solid)and stellar (dashed) particles for m12i (left), m12f (middle), and m12m (right) systems.The color code is the same as in Fig. 1. We also plot available DMO runs for each halo (dot-dashed) for comparison.
Bottom : The corresponding total circularvelocity profiles for the simulated systems in the top panel. The dots specify the circular velocity at 2 kpc, i.e. 𝑉 . tion of SIDM cross section, and star formation history and baryonicconcentrations controls the evolution of the SIDM central densities.Haloes with different SIDM cross sections show quite differ-ent responses to the stellar mass assembly. In all the SIDM haloeswith 𝜎 / 𝑚 =
10 cm g − , the inner DM profile has higher cen-tral density and is cuspier than in the CDM and SIDM1 counter-parts. This reflects the fact that higher DM cross section leads tohigher interaction rate Γ ∝ 𝜎 / 𝑚 , and hence more efficient ther-malization process. On the other hand, the three SIDM haloes with 𝜎 / 𝑚 = g − show more comparable central densities to theCDM haloes: in m12i, the SIDM density profiles are both denserand cuspier than the CDM run, while the m12f-SIDM1 density pro-file is only slightly flatter compared to its CDM version, and them12m run has an extended core with similar density. In all cases,higher (lower) central DM densities are accompanied by higher(lower) amplitudes for the stellar density profiles. The co-evolutionof SIDM with 𝜎 / 𝑚 = g − and baryons in our baryonic sim-ulations generates slightly larger diversity in DM central density(measured by 𝜌 DM ( . ) ; see table 1) compared to the CDMruns (see also Elbert et al. 2018; Sameie et al. 2018; Robles et al.2019). This is very different from the DMO simulations, which result in constant central densities with amplitudes that decreasemonotonically with cross section.In the lower row of Fig. 2, total circular velocity profiles, 𝑉 tot = √︁ 𝐺 𝑀 tot ( < 𝑟 )/ 𝑟 = √︃ 𝑉 + 𝑉 ★ + 𝑉 , are plotted for theCDM and SIDM simulations. In order to reassess the diversity in thecentral DM distributions, we adopt the quantity 𝑉 , introducedby Oman et al. (2015), which is the total circular velocity at 2 kpc(filled circles). Again, SIDM models with 𝜎 / 𝑚 = g − show asimilar diversity in the distribution of 𝑉 compared to the CDMhaloes. The SIDM runs with cross section 10 cm g − have lessdiversity in their 𝑉 distribution, reflecting the very high centraldensities described above. The sample of our baryonic simulationssuggest that SIDM models with cross section 𝜎 / 𝑚 = g − , aswell as the CDM+feedback, can generate diverse DM distributionsin the MW-mass scale (see also Santos-Santos et al. 2018). A largersample of galaxies simulated in both CDM and SIDM models, withmore variety in their mass assembly and star formation history, andat different mass scales relevant to the diversity problem will becrucial to confirm our results. MNRAS , 1–10 (2020)
Sameie et al.
In order to better understand the interplay between baryons and DM,Fig. 3 shows the redshift evolution of the DM density (top), stellardensity (middle), and DM velocity dispersion (bottom) profiles (allin physical coordinates) for the CDM and SIDM models. Resultsfor 𝑧 > . .
75 dex for density anddispersion profiles respectively, for clarity. At 𝑧 = 𝑧 (cid:46) Γ ∝ 𝜎 / 𝑚 ) anda faster transition from core-expansion to core-contraction phase.A comparison of the three SIDM1 simulations reveals vary-ing amounts of halo contraction in the central region. We identifydifferent growth rates in the DM central densities to correlate withdifferent stellar formation histories of the inner regions. The mid-dle row of Fig. 3 shows the redshift-evolution of the stellar densityprofiles. At 𝑧 =
2, the stellar distributions in the CDM simulationshave vastly different central densities (nearly a factor 100) and innerslopes, reflecting the variety of star formation rates for each galaxy(Santistevan et al. 2020). The stellar profiles in the SIDM model areless centrally dense and more diffuse at 𝑧 = 𝑧 =
2, form more starsat later times; by 𝑧 =
0, their stellar profiles are denser near the cen-ter compared to the CDM run. The transition from cored to cuspyDM profiles happens faster by increasing SIDM cross section. Allthree SIDM10 models with lower DM central densities at 𝑧 = 𝑧 =
0, while two out of three SIDM1 models still have lower orsimilar central densities compared to the CDM models. The higherDM cross section also leads to higher stellar central densities in theSIDM10 galaxies (see section 3.4 for further discussion).The bottom row of Fig. 3 shows the evolution of the 1D DMvelocity dispersion profiles for the CDM and SIDM models. Forcomparison, we also show 𝑧 = 𝑧 =
2, both CDM and SIDM models havecomparable velocity dispersion in all radii with positive slopes, in-dicating that DM self-interactions are yet to fully thermalize thecentral regions. Subsequent assembly of stars and DM increase theamplitude of velocity dispersion for both CDM and SIDM simu-lations. However, redistribution of energy and momentum in theSIDM haloes creates isothermal cores (i.e. flat central velocity dis-persion). The relatively smaller halo-to-halo variation in the centralvelocity dispersion of the SIDM haloes signifies the role of ther-malization process: while the central regions in the CDM haloes isheated up in response to the assembly of the baryons, the SIDMthermalization process redistributes heat such that the DM densitycores remain isothermal. In addition, while the DM velocity disper- sion in the baryonic versions of the SIDM simulations is higher thanin their DMO versions (middle panel), the presence of baryons doesnot prevent the development of an isothermal core (i.e. flat velocitydispersion) induced by DM self-interactions.In summary, we find that while in CDM the inclusion ofbaryons flattens the central density profile compared to the DMOversion, in SIDM the dominance of baryons in the central regionscan make the SIDM haloes cuspier than the CDM haloes. Either adenser stellar components or higher DM cross section lead to morecuspy SIDM profile, making the SIDM models more responsive tothe presence of the baryons. Hence, the spread in the density ofSIDM haloes is tied to the baryonic concentration and local starformation in the central few kpc of MW-mass galaxies.
Comparisons within our sample of simulated galaxies suggest thatthe concentration of the stellar distribution is more important thanthe total disc mass in creating diverse SIDM density profiles. AmongSIDM1 simulations, the m12m galaxy is more massive than bothm12i and m12f, and yet it has the most extended stellar disc with thelargest stellar half-mass radius and lowest central DM density (table1). The lower DM and baryonic concentration in the m12m SIDM1halo gives rise to a central gravitational potential Φ pot ( 𝑟 = ) thatis shallower than the other two simulations (table 1). This picture isconsistent with the analytical model of (Kaplinghat et al. 2014), inwhich the SIDM density profile is related to the total gravitationalpotential by 𝜌 SIDM ∝ exp (− Φ pot / 𝜎 ) ( 𝜎 is the central DM ve-locity dispersion) and suggests that a larger (smaller) Φ tot leads toa higher (lower) central density in the SIDM haloes. In addition,controlled 𝑁 -body simulations of SIDM haloes have shown thatSIDM systems with higher baryonic concentration will transitionfaster from a “core-expansion" phase to “core-contraction" phase(Elbert et al. 2018; Sameie et al. 2018).As discussed in the previous section, SIDM10 models showsuch a transition from cored to cuspy profiles over the course ofcosmological evolution (see top row panels of Fig. 3). We find acorrelation in the assembly of baryons and DM around the centers ofour simulated galaxies. In the top panels of Fig. 4 we show the red-shift evolution of the mean stellar densities 𝜌 ★ = 𝑀 ( < 𝑟 )/( 𝜋𝑟 / ) of the CDM (black) and SIDM10 (red) haloes computed at 3D half-mass radius 𝑟 / . The middle and bottom rows show the redshiftevolution of DM central density computed at 𝑟 = . ± .
05 kpcand the DM density slope computed at 𝑟 = . - . At 𝑧 = 𝑧 (cid:46) 𝑧 ∼ . - .
5. The growth We compute the gravitational potential by direct summation of all particlesinside of virial radius, neglecting the contribution of boundary and distantparticles. We have checked that our results are not sensitive to the radius we choseto compute DM densities, and the cumulative DM densities computed, forexample, at 3D stellar half-mass radius also follow the same trend.MNRAS000
5. The growth We compute the gravitational potential by direct summation of all particlesinside of virial radius, neglecting the contribution of boundary and distantparticles. We have checked that our results are not sensitive to the radius we choseto compute DM densities, and the cumulative DM densities computed, forexample, at 3D stellar half-mass radius also follow the same trend.MNRAS000 , 1–10 (2020)
IDM haloes On FIRE D M ( M k m ) m12i SIDM10SIDM1CDM m12f m12m z = 0z = 1z = 21 1010 ( M k m ) v , d m ( k m s ) SIDM10
SIDM1CDM
Distance (kpc) z = 0z = 1z = 2
Figure 3.
Cosmological evolution of DM density (top), stellar density (middle) and 1D DM velocity dispersion (bottom) for m12i (left column), m12f (middlecolumn) and m12m (right column) haloes. Colored curves represent different redshifts. Density and dispersion profiles at 𝑧 > 𝑧 = of the stellar densities in SIDM10 haloes come from both increasein the central stellar mass and decrease in 𝑟 / relative to the CDMmodels; the latter effect is a response to SIDM halo contraction.A higher baryon concentration leads to higher SIDM interactionrates and more halo contraction. Our results are in agreement withthe previous works and affirm that in the presence of massive anddense baryonic distribution the SIDM haloes can develop higher and more cuspy DM density profiles. In fully isolated halos, thiscombination of galaxy and halo contraction generates a run-awayprocess that eventually leads to core collapse (Balberg et al. 2002;Koda & Shapiro 2011; Sameie et al. 2018; Essig et al. 2019).Thus far, we have discussed the importance of baryon massassembly in shaping the DM density profiles at late times ( 𝑧 (cid:46) 𝑧 >
2, the baryonic potential does not contribute
MNRAS , 1–10 (2020)
Sameie et al. ( r / ) m12i m12f m12m d m ( . k p c ) ( . k p c ) Redshift
CDMSIDM10
Figure 4. Top : The redshift evolution of the 3D stellar density, 𝜌 ★ , computed at 3D half-mass radius 𝑟 / for the CDM (black) and SIDM10 (red) galaxies since 𝑧 = Middle : DM density, computed at 𝑟 = . ± .
05 kpc for the CDM and SIDM10 haloes.
Bottom : DM density slope, computed between 𝑟 = . − . significantly to the total gravitational potential, and therefore it can-not affect the evolution of central densities. In Fig. 5 we compute theredshift evolution of DM central densities computed at 𝑟 = . ± . 𝑟 = . - 𝑧 =
6. We also compute these quantities for theDMO version of these haloes (dashed). At high redshifts, densityprofiles of the simulated haloes in full physics and DMO suitesshow good agreement in their amplitude and slope for both CDMand SIDM models. In the top panel, the central density in the CDMhalo remains quite similar between full physics and DMO suite forthe full redshift range plotted, while SIDM halo in the full physicssuite shows a significant boost in its central density after 𝑧 = 𝑧 = 𝑀 ★ / 𝑀 ∼ − atthat redshift) agrees well with its DMO version, while by 𝑧 = 𝑀 ★ / 𝑀 ∼ − ), the DM density profile becomes significantlyflattened owing to stellar feedback.Our results for the CDM halo are consistent with those inLazar et al. (2020), in which haloes with stellar-halo mass ratiosof 𝑀 ★ / 𝑀 ∼ − have DM density slopes roughly equal to theDMO simulations, while haloes with 𝑀 ★ / 𝑀 ∼ − have muchflatter DM density slopes. For the SIDM10 halo, the DM densityslope also agrees well with its SIDM-only version at 𝑧 >
3; as aresult of the self-interactions, both have much shallower inner den-sity profiles compared to the CDM halo (see also Fig. 3). At lowerredshifts, the CDM and SIDM10 density profiles evolve strikinglydifferently: around 𝑧 =
2, the time at which the stellar feedbackbegins to flatten the DM density profile, the SIDM10 halo starts to develop a steeper density profile (due to interplay between baryonsand SIDM thermalization), and it eventually becomes cuspier thanthe CDM version of the same galaxy. Evidently, the core forma-tion process is more controlled by the DM self-interactions than thestellar feedback in our SIDM simulation.In summary, while the CDM central densities at 𝑟 = . 𝑧 (cid:46) 𝜎 / 𝑚 =
10 cm / g)to the presence of the baryons results in enhanced overall centraldensities relative to CDM. We perform a suite of baryonic cosmological zoom-in simulationswith self-interacting dark matter (SIDM) models within the “Feed-back In Realistic Environment" (FIRE) project. The treatment of thebaryonic physics includes cooling, star formation, stellar feedbackfrom SNe, photo-heating, stellar winds, radiation pressure, and UV-background radiation. We take three Milky Way-mass dark matter(DM) haloes from the Latte suite and re-simulate them with SIDMmodels with constant DM cross sections 𝜎 / 𝑚 = ,
10 cm g − . Westudy the stellar mass assembly and the evolution of DM densityprofiles in our SIDM simulations compared to the CDM versions. MNRAS000
10 cm g − . Westudy the stellar mass assembly and the evolution of DM densityprofiles in our SIDM simulations compared to the CDM versions. MNRAS000 , 1–10 (2020)
IDM haloes On FIRE d m ( r = . k p c ) DM + baryonsDMOSIDM10CDM m12f
Redshift ( r = . k p c ) Figure 5.
Redshift evolution of DM central densities computed at 𝑟 = . ± .
05 kpc (top) and the DM density slope measured at 𝑟 = - . 𝑧 =
6, for the m12f simulated galaxies with full baryonic physics (solid)DMO (dashed) in CDM (black) and SIDM10 (red). While the SIDM haloshave cores even at the highest redshifts probed here, they end up with cuspyprofiles at 𝑧 = 𝑧 ≈ The CDM and SIDM haloes have almost identical halo assemblieswith similar virial masses at 𝑧 =
0. Star formation rates (SFR) inboth CDM and SIDM haloes are bursty and consist of two distinctphases: rapid increase in the SFR followed by a plateau from 𝑧 ≤ 𝑧 >
4, CDM and SIDM haloes in both FIREand DMO suite show good agreement in their central DM densitiesand slopes at high redshifts. However, at later times the CDM haloestend to become more flattened (due to stellar feedback) while theour SIDM haloes become denser and cuspier. This different evolu- tionary phases in the CDM and SIDM suite mark the importance ofthe thermalization process by the DM self-interactions.Overall, our results suggest that the SIDM haloes are creatingdiverse DM central densities from cored to cuspy density profilesand they could be even cuspier and denser than their CDM counter-parts. We observe a correlation between the growth rate of centralbaryonic distribution and the DM central densities, and we find thatSIDM haloes are more responsive to the presence of the baryonsthan the CDM haloes.
ACKNOWLEDGEMENTS
We thank Hai-Bo Yu for insightful discussions. MBK acknowl-edges support from NSF CAREER award AST-1752913, NSFgrant AST-1910346, NASA grant NNX17AG29G, and HST-AR-15006, HST-AR-15809, HST-GO-15658, HST-GO-15901, HST-GO-15902, HST-AR-16159, and HST-GO-16226 from the SpaceTelescope Science Institute, which is operated by AURA, Inc., un-der NASA contract NAS5-26555. RES acknowledges support fromNASA grant 19-ATP19-0068, the Research Corporation for Scien-tific Advancement through Scialog-TDA, and grant HST-AR-15809from the Space Telescope Science Institute, which is operated byAURA, Inc., under NASA contract NAS5-26555. Support for PFHwas provided by NSF Research Grants 1911233 & 20009234, NSFCAREER grant 1455342, NASA grants 80NSSC18K0562, HST-AR-15800.001-A. Numerical calculations were run on the Caltechcompute cluster “Wheeler,” allocations FTA-Hopkins/AST20016supported by the NSF and TACC, and NASA HEC SMD-16-7592. AW received support from NASA through ATP grants80NSSC18K1097 and 80NSSC20K0513; HST grants GO-14734,AR-15057, AR-15809, and GO-15902 from STScI; the Heising-Simons Foundation; and a Hellman Fellowship. ASG is supportedby the McDonald Observatory via the Harlan J. Smith postdoctoralfellowship. We ran simulations using: XSEDE, supported by NSFgrant ACI-1548562; Blue Waters, supported by the NSF; Pleiades,via the NASA HEC program through the NASA Division at AmesResearch Center. The analysis in this paper is carried out by pythonpackages Numpy (van der Walt et al. 2011), matplotlib (Hunter2007), scipy (Oliphant 2007), and h5py (Collette 2013).
DATA AVAILABILITY
The data supporting the plots within this article are available onreasonable request to the corresponding author.
REFERENCES
Ackerman L., Buckley M. R., Carroll S. M., Kamionkowski M., 2009, Phys.Rev. D, 79, 023519Ahn K., Shapiro P. R., 2005, MNRAS, 363, 1092Arkani-Hamed N., Finkbeiner D. P., Slatyer T. R., Weiner N., 2009, Phys.Rev. D, 79, 015014Balberg S., Shapiro S. L., Inagaki S., 2002, ApJ, 568, 475Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109Benítez-Llambay A., Frenk C. S., Ludlow A. D., Navarro J. F., 2019, MN-RAS, 488, 2387Blumenthal G. R., Faber S. M., Primack J. R., Rees M. J., 1984, Nature,311, 517Bond J. R., Szalay A. S., Turner M. S., 1982, Phys. Rev. Lett., 48, 1636Bose S., et al., 2019, MNRAS, 486, 4790Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40MNRAS , 1–10 (2020) Sameie et al.
Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2012, MNRAS, 422, 1203Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C. A., Quataert E., 2015, MNRAS, 454, 2981Collette A., 2013, Python and HDF5: unlocking scientific data. " O’ReillyMedia, Inc."Correa C. A., 2020, arXiv e-prints, p. arXiv:2007.02958Creasey P., Sameie O., Sales L. V., Yu H.-B., Vogelsberger M., Zavala J.,2017, MNRAS, 468, 2283Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371Despali G., Sparre M., Vegetti S., Vogelsberger M., Zavala J., Marinacci F.,2019, MNRAS, 484, 4563Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., DuttonA. A., Wadsley J., 2014, MNRAS, 437, 415Di Cintio A., Tremmel M., Governato F., Pontzen A., Zavala J., BastidasFry A. e., Brooks A., Vogelsberger M., 2017, MNRAS, 469, 2845El-Badry K., Wetzel A., Geha M., Hopkins P. F., Kereš D., Chan T. K.,Faucher-Giguère C.-A., 2016, ApJ, 820, 131Elbert O. D., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., GrausA. S., Rocha M., 2018, ApJ, 853, 109Essig R., McDermott S. D., Yu H.-B., Zhong Y.-M., 2019, Phys. Rev. Lett.,123, 121102Feng J. L., Kaplinghat M., Tu H., Yu H.-B., 2009, J. Cosmology Astropart.Phys., 2009, 004Fitts A., et al., 2017, MNRAS, 471, 3547Fitts A., et al., 2019, MNRAS, 490, 962Fry A. B., et al., 2015, MNRAS, 452, 1468Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Kirby E. N., 2014,MNRAS, 444, 222Garrison-Kimmel S., et al., 2019, MNRAS, 487, 1380Governato F., et al., 2012, MNRAS, 422, 1231Hahn O., Abel T., 2011, MNRAS, 415, 2101Hopkins P. F., 2015, MNRAS, 450, 53Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E.,Murray N., Bullock J. S., 2014, MNRAS, 445, 581Hopkins P. F., et al., 2018a, MNRAS, 477, 1578Hopkins P. F., et al., 2018b, MNRAS, 480, 800Hunter J. D., 2007, Computing in Science and Engineering, 9, 90Kahlhoefer F., Kaplinghat M., Slatyer T. R., Wu C.-L., 2019, J. CosmologyAstropart. Phys., 2019, 010Kamada A., Kaplinghat M., Pace A. B., Yu H.-B., 2017, Phys. Rev. Lett.,119, 111102Kaplinghat M., Keeley R. E., Linden T., Yu H.-B., 2014, Phys. Rev. Lett.,113, 021302Kaplinghat M., Valli M., Yu H.-B., 2019, MNRAS, 490, 231Katz N., White S. D. M., 1993, ApJ, 412, 455Koda J., Shapiro P. R., 2011, MNRAS, 415, 1125Kuzio de Naray R., McGaugh S. S., de Blok W. J. G., 2008, ApJ, 676, 920Lazar A., et al., 2020, MNRAS, 497, 2393Lelli F., McGaugh S. S., Schombert J. M., 2016, AJ, 152, 157Loeb A., Weiner N., 2011, Phys. Rev. Lett., 106, 171302Lovell M. R., et al., 2018, MNRAS, 477, 2886Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174Monaghan J. J., Lattanzio J. C., 1985, A&A, 149, 135Moore B., 1994, Nature, 370, 629Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, L72Oñorbe J., Garrison-Kimmel S., Maller A. H., Bullock J. S., Rocha M.,Hahn O., 2014, MNRAS, 437, 1894Oñorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., Kereš D.,Faucher-Giguère C.-A., Quataert E., Murray N., 2015, MNRAS, 454,2092Oh S.-H., et al., 2015, AJ, 149, 180Oliphant T. E., 2007, Computing in Science and Engineering, 9, 10Oman K. A., et al., 2015, MNRAS, 452, 3650Peebles P. J. E., 1982, ApJ, 263, L1Percival W. J., et al., 2007, ApJ, 657, 645Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS,430, 105 Planck Collaboration et al., 2018, arXiv e-prints, p. arXiv:1807.06209Pontzen A., Governato F., 2012, MNRAS, 421, 3464Read J. I., Walker M. G., Steger P., 2019, MNRAS, 484, 1401Ren T., Kwa A., Kaplinghat M., Yu H.-B., 2019, Physical Review X, 9,031020Robertson A., et al., 2018, MNRAS, 476, L20Robertson A., Harvey D., Massey R., Eke V., McCarthy I. G., Jauzac M., LiB., Schaye J., 2019, MNRAS, 488, 3646Robles V. H., Kelley T., Bullock J. S., Kaplinghat M., 2019, MNRAS, 490,2117Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-KimmelS., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81Sameie O., Creasey P., Yu H.-B., Sales L. V., Vogelsberger M., Zavala J.,2018, MNRAS, 479, 359Sameie O., Benson A. J., Sales L. V., Yu H.-b., Moustakas L. A., CreaseyP., 2019, ApJ, 874, 101Sameie O., Chakrabarti S., Yu H.-B., Boylan-Kolchin M., Vogelsberger M.,Zavala J., Hernquist L., 2020a, arXiv e-prints, p. arXiv:2006.06681Sameie O., Yu H.-B., Sales L. V., Vogelsberger M., Zavala J., 2020b, Phys.Rev. Lett., 124, 141102Samuel J., et al., 2020, MNRAS, 491, 1471Santistevan I. B., Wetzel A., El-Badry K., Bland-Hawthorn J., Boylan-Kolchin M., Bailin J., Faucher-Giguère C.-A., Benincasa S., 2020, MN-RAS, 497, 747Santos-Santos I. M., Di Cintio A., Brook C. B., Macciò A., Dutton A.,Domínguez-Tenreiro R., 2018, MNRAS, 473, 4392Seljak U., Makarov A., McDonald P., Trac H., 2006, Phys. Rev. Lett., 97,191303Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760Springel V., 2005, MNRAS, 364, 1105Tollet E., et al., 2016, MNRAS, 456, 3542Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1Tulin S., Yu H.-B., Zurek K. M., 2013, Phys. Rev. D, 87, 115007Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740Vogelsberger M., et al., 2014a, MNRAS, 444, 1518Vogelsberger M., Zavala J., Simpson C., Jenkins A., 2014b, MNRAS, 444,3684Vogelsberger M., et al., 2014c, Nature, 509, 177Vogelsberger M., Zavala J., Cyr-Racine F.-Y., Pfrommer C., Bringmann T.,Sigurdson K., 2016, MNRAS, 460, 1399Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D.,Quataert E., 2016, ApJ, 827, L23Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Scienceand Engineering, 13, 22This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000