The impact of inelastic self-interacting dark matter on the dark matter structure of a Milky Way halo
Kun Ting Eddie Chua, Karia Dibert, Mark Vogelsberger, Jesús Zavala
MMNRAS , 1–16 (2020) Preprint 20 October 2020 Compiled using MNRAS L A TEX style file v3.0
The impact of inelastic self-interacting dark matter on the darkmatter structure of a Milky Way halo
Kun Ting Eddie Chua (cid:63) , Karia Dibert , Mark Vogelsberger , and Jesús Zavala Institute of High Performance Computing, 1 Fusionopolis Way Department of Physics, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA Center for Astrophysics and Cosmology, Science Institute, University of Iceland, Dunhagi 5, 107 Reykjavik, Iceland
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We study the effects of inelastic dark matter self-interactions on the internal structure of asimulated Milky Way (MW)-size halo. Self-interacting dark matter (SIDM) is an alternativeto collisionless cold dark matter (CDM) which offers a unique solution to the problemsencountered with CDM on sub-galactic scales. Although previous SIDM simulations havemainly considered elastic collisions, theoretical considerations motivate the existence of multi-state dark matter where transitions from the excited to the ground state are exothermic. In thiswork, we consider a self-interacting, two-state dark matter model with inelastic collisions,implemented in the Arepo code. We find that energy injection from inelastic self-interactionsreduces the central density of the MW halo in a shorter timescale relative to the elastic scale,resulting in a larger core size. Inelastic collisions also isotropize the orbits, resulting in anoverall lower velocity anisotropy for the inelastic MW halo. In the inner halo, the inelasticSIDM case (minor-to-major axis ratio s ≡ c/a ≈ . ) is more spherical than the CDM( s ≈ . ), but less spherical than the elastic SIDM case ( s ≈ . ). The speed distribution f ( v ) of dark matter particles at the location of the Sun in the inelastic SIDM model shows asignificant departure from the CDM model, with f ( v ) falling more steeply at high speeds. Inaddition, the velocity kicks imparted during inelastic collisions produce unbound high-speedparticles with velocities up to 500 km s − throughout the halo. This implies that inelasticSIDM can potentially leave distinct signatures in direct detection experiments, relative toelastic SIDM and CDM. Key words: cosmology: dark matter – galaxies: haloes – methods: numerical
The standard cold dark matter (CDM) model has been successfulin describing several properties of the universe, such as the Cos-mic Microwave Background, and the formation and evolution oflarge-scale structure in the distribution of galaxies throughout theuniverse (e.g. Spergel et al. 2003; Springel et al. 2005). On the otherhand, relevant discrepancies have been uncovered between CDMand observations at sub-galactic scales. Some of these include: (i)the cusp-core problem, arising from a disagreement between thepredicted dark matter (DM) density profiles of low-mass galaxiesand their observed density profiles (e.g. de Blok & McGaugh 1997;Walker & Peñarrubia 2011); (ii) the missing satellites problem, aris-ing from a discrepancy between the abundance of low-mass galaxiesin the Milky Way (MW) and that “naively" predicted by CDM sim-ulations (e.g. Moore et al. 1999; a similar problem appears in fieldgalaxies, e.g. Zavala et al. 2009; Klypin et al. 2015); (iii) the too-big-to-fail problem highlighted by Boylan-Kolchin et al. (2011), (cid:63)
Email: [email protected] which arises from the fact that low-mass CDM (sub)haloes, whichare expected to host low-mass galaxies, are seemingly too denseto explain the observed stellar kinematics of these systems; (iv) thedwarf rotation curve diversity problem, where the rotation curves ofsimulated dwarf galaxies at a fixed mass do not show the large vari-ety/diversity observed in real galaxies (Oman et al. 2015) ; and (v)conflicting expectations of the dark matter halo shape, where CDMsimulations predict average shapes far less spherical than those de-rived from tidal stream observations (e.g. Ibata et al. 2001; Law &Majewski 2010; Vera-Ciro & Helmi 2013; Bovy et al. 2016).It is important to remark that these small-scale disagreementsare only firmly established from the results of N -body simulations,which include only the effects of DM gravity. Recent efforts haveaccelerated the development of hydrodynamics simulations, e.g.Illustris (Vogelsberger et al. 2014), EAGLE (Crain et al. 2015),Illustris-TNG (Pillepich et al. 2018), where the incorporation of We note that with recent observations of ultra-faint galaxies, the too-big-to-fail problem also becomes a diversity problem for the broad distributionof stellar kinematics in dwarf spheroidals in the MW (Zavala et al. 2019) © a r X i v : . [ a s t r o - ph . GA ] O c t K.T.E. Chua et. al. galaxy formation processes allows the complex modelling of galax-ies to occur in a full cosmological context. Due to the coupling be-tween DM and baryons through gravity, it has been suggested thatbaryonic effects in hydrodynamic simulations can alleviate some,if not all, of the CDM challenges. For example, energy from su-pernova feedback can inject energy into the inner dark matter halo,alleviating the cusp-core problem (Navarro et al. 1996; Governatoet al. 2012; Oñorbe et al. 2015; Read et al. 2016), as well as thetoo-big-to-fail problem (Zolotov et al. 2012; Brook & Di Cintio2015; Sawala et al. 2016; Wetzel et al. 2016). It has also beenfound that while a combination of supernova and active galacticnuclei (AGN) feedback can reduce the central densities of haloes(e.g. Duffy et al. 2010), the condensation of baryons in halo cen-tres conversely increases the halo concentration through adiabaticcontraction (Blumenthal et al. 1986; Gnedin et al. 2004). The con-densation of baryons can also affect halo shapes, leading to morespherical DM haloes in hydrodynamic simulations compared to N -body simulations (e.g. Katz & Gunn 1991; Dubinski 1994; Tisseraet al. 2010; Abadi et al. 2010; Bryan et al. 2013; Butsky et al. 2016;Chisari et al. 2017; Chua et al. 2019). The impact of all of thesebaryonic processes certainly reduces the tension between the CDMmodel and observations of the dwarf galaxy population. However,it is important to note that a wide variety of physics implementa-tions exist, resulting in uncertain predictions. If these uncertaintiesin baryonic physics must be taken into account, we argue that is alsoprudent to consider another uncertainty in the physics of galaxies,which is the nature of the dark matter particle. Alternatives to darkmatter models can arise by relaxing the key assumptions of CDM- that dark matter is both cold and collisionless. For a recent re-view on different dark matter models and their impact on structureformation, see Zavala & Frenk (2019); for a review of the CDMchallenges and possible solutions, see Bullock & Boylan-Kolchin(2017).Relaxing the assumption that dark matter is cold for the pur-poses of galaxy formation implies that the primordial power spec-trum has deviations over CDM at galactic scales. This can happenfor example in models such as warm dark matter (WDM), in whichdark matter particles undergo free-streaming in the early Universe(Colín et al. 2000; Bode et al. 2001), as well as interacting darkmatter, in which dark matter particles interact with relativistic par-ticles in the early Universe (Bœhm et al. 2002; Buckley et al. 2014).The cutoff in the linear power spectrum reduces the severity of themissing satellite problem (Bœhm et al. 2014), the too-big-to-failproblem and the dwarf rotation curve diversity problem (Vogels-berger et al. 2016; Zavala et al. 2019). However, observations of theLy- α forest constrain the mass of thermal WDM to m WDM (cid:38) . keV (Viel et al. 2013; Iršič et al. 2017; although see Garzilli et al.2019). Although WDM also naturally predicts central density coresdue to the upper bound on the phase space density set by the pri-mordial thermal velocity dispersion (Dalcanton & Hogan 2001),allowed WDM models cannot create large dark matter cores with-out severely under-predicting the abundance of low-mass galaxies(Macciò et al. 2012). Similarly, the interacting dark matter modeldoes not alleviate the cusp-core problem since the interactions be-tween dark matter and the relativistic particles decouple long beforethe onset of dark matter haloes.Alternatively, models in which dark matter particles interactwith each other are known as Self-Interacting Dark Matter (SIDM,Spergel & Steinhardt 2000) models. SIDM has self-interaction crosssections with an amplitude near that of strong nucleon-nucleonelastic scattering, which are sufficient to reduce the central densitiesof DM haloes and alleviate the CDM challenges (for a review, see Tulin & Yu 2018). For example, by transferring energy from theouter regions of the halo to the inner regions, SIDM models cancreate cores on kpc scales (Colín et al. 2002). This alleviates boththe cusp-core and the too-big-to-fail problems (Vogelsberger et al.2012; Rocha et al. 2013; Zavala et al. 2013), and also the dwarfrotation curve diversity problem (Kamada et al. 2017; although seeSantos-Santos et al. 2020).Most SIDM simulations thus far have assumed that the scat-tering process between two SIDM particles is elastic, i.e. kineticenergy is conserved during the collision. In these purely elasticSIDM models, elastic scattering leads to a redistribution of darkmatter particles within the halo, which has been found to modifyDM haloes in terms of their phase-space structure (Vogelsberger &Zavala 2013) and halo shapes (Peter et al. 2013; Brinckmann et al.2018). However self-scattering interactions can also be inelastic andoccur in theoretical models which contain multi-state dark matter(Arkani-Hamed et al. 2009; Schutz & Slatyer 2015). For the caseof a two-state scenario, a transition from the ground ( χ ) to theexcited ( χ ) state (up-scattering) is an endothermic process, whilethe reverse process (down-scattering) is exothermic. Inelastic SIDMis especially interesting from a core-formation perspective, since adown-scattering event produces a kick in the velocities of the groundstate particles. For example, Todoroki & Medvedev (2019a,b, 2020)performed SIDM simulations using a multi-state dark matter modeland concluded that elastic scattering and energy injection are in-dependently sufficient to create isothermal cores in haloes of mass ≈ × M (cid:12) . The energy injected into the halo from exothermicdown-scattering can be on the order of 100 million Type II super-novae, which not only results in higher core formation efficienciesthan a purely elastic model, but also reduces the amount of substruc-ture present in the halo (Vogelsberger et al. 2019, hereafter V19).As such, inelastic SIDM can increase the allowed parameter spaceof self-interaction cross sections in SIDM models. So far, the im-pact of energy injection due to inelastic self-interactions on the darkmatter properties of galactic haloes remains largely unexplored.In this paper, we examine the effects of inelastic SIDM on halostructure and assembly using high-resolution DM-only simulationsof a MW-size halo. The study is based on a comparative analysisbetween a CDM model, an elastic SIDM model, and an inelasticSIDM model. By using DM-only simulations, we focus solely onthe impact of SIDM in the absence of gas physics and galaxy for-mation. The paper is structured as follows: Section 2 describes thenumerical simulations, including a description of the two-state in-elastic SIDM model and its implementation. We then present theresults of the simulations, focusing on the effects of inelastic SIDMon the structure of the simulated halo at z = 0 in Section 3, aswell as its assembly history in Section 4. Finally, we present ourconclusions in Section 5. Since the haloes analysed in this paper were previously introducedin V19, we refer the reader there for a full description of the modeland methods. We give a brief overview of the relevant points here.
Our SIDM simulations are based on the two-state model presentedby Schutz & Slatyer (2015), where the excited state is nearly degen-erate with the ground state, and analytic expressions for the elasticand inelastic s-wave cross sections have been derived. Although
MNRAS , 1–16 (2020) mpact of inelastic SIDM on a MW halo other inelastic SIDM models exist, this model is currently the onlyone in which an analytic description of scattering has proven to befeasible at dwarf galaxy velocity scales. The model is described bythe following parameters: the mass splitting δ between the ground( χ ) and excited ( χ ) states, the coupling constant α between thedark matter particle and the force mediator, the mass m χ of thedark matter particle, and the mass m φ of the force mediator. Inthis work, these parameters have the values δ = 10 keV, α = 0 . , m χ = 10 GeV, and m φ = 30 MeV. While these values representan arbitrary point within the entire parameter space, they have beenchosen to lie within the range for interesting cross sections for ve-locities of the order of 10 km s − (Schutz & Slatyer 2015). At thescale of the MW halo, the corresponding elastic cross section perunit mass is a few cm g − . Such cross sections have previouslybeen found to be capable of creating cores of size O (kpc) (e.g.Vogelsberger et al. 2012; Brinckmann et al. 2018)There are five possible interactions in this two-state SIDMmodel, namely:(i) elastic scattering of two ground state particles ( χ + χ → χ + χ ) ,(ii) elastic scattering of two excited state particles ( χ + χ → χ + χ ) ,(iii) elastic Yukawa scattering ( χ + χ → χ + χ ) ,(iv) inelastic endothermic up-scattering in which two groundstate particles transition to the excited state ( χ + χ → χ + χ ) ,and(v) inelastic exothermic down-scattering in which two excitedparticles transition to the ground state ( χ + χ → χ + χ ) .During down-scattering, our chosen model produces a velocity kick v kick = (cid:112) δ/m χ c (cid:39) km s − . Up-scattering can only occurfor relative velocities v rel > v kick (cid:39) km s − . For typicaldark matter velocities in a MW-size halo ( ≈ km s − ), inelasticendothermic scattering is essentially forbidden since the energysplitting δ is too large for this interaction to occur frequently.An additional parameter which must be specified in the simula-tions is the primordial fraction of dark matter in each of the groundand excited states. A primordial excitation fraction χ = 100% corresponds to having all particles being in the excited state, andleads to the maximum possible energy release during structureformation and hence the maximum possible effect on halo struc-ture. Conversely, an inelastic system with all particles initiallyin the ground state behaves essentially like purely elastic SIDM,since inelastic up-scattering is suppressed at galactic halo ve-locities. In this paper, we examine two initial configurations: (i) χ = 100% , where all SIDM particles begin in the excited state,and (ii) χ = 50% , where only half of the SIDM particles beginin the excited state. For possible theoretical justifications for thelarge primordial fraction of excited states, we refer the reader toV19. The two-state inelastic SIDM model is implemented within a gen-eral multi-state dark matter framework in the Arepo code (Springel2010). This framework is a generalisation of the probabilistic ap-proach presented in Vogelsberger et al. (2012) and Vogelsbergeret al. (2016), and is able to handle an arbitrary number of stateswith non-degenerate energy level splittings and all possible reac-tions between them, each with any given velocity-dependent crosssection.Each dark matter particle i is assumed to be in a specific state α . The simulation volume is populated by dark matter particles invarious states ( α, β, γ, δ ) with possible two-body scatterings: χ αi + χ βj → χ γi + χ δj . (1)Equation 1 represents the scattering of particles i and j from states α and β into states γ and δ . A dark matter particle in state (cid:15) ∈ ( α, β, γ, δ ) has mass m χ (cid:15) . The scattering rates for the possiblereactions are given by: R αβ → γδ = ρ β m χ β (cid:104) σ αβ → γδT ( v αβ ) v αβ (cid:105) , (2)where ρ β is the local mass density of particles in state β , σ αβ → γδT ( v αβ ) is the velocity-dependent transfer cross section forthe reaction, and v αβ is the magnitude of the relative velocity be-tween the interacting particles. Since the differential cross sectionhas no angular dependence in this model, the transfer cross sectionis also the same as the total cross section for the reaction.The scattering of particles i and j in states α and β withmasses m αi and m βj is performed in the centre of mass frame, andnew velocities are assigned to each particle: v i = m αi + m βj m γi + m δj v cm + m δj m γi + m δj A ij v ij ˆ e , (3) v j = m αi + m βj m γi + m δj v cm − m γi m γi + m δj A ij v ij ˆ e , (4)where v cm is the velocity of the centre of mass, ˆ e is a random vectoron the unit sphere, and A ij ( αβ → γδ ) is a dimensionless velocityscale factor that depends on the energy splitting of the reaction: ≤ A ij = 1 elastic < inelastic: endothermic > inelastic: exothermic (5)For inelastic scattering, each particle in the reaction gains or losesthe same amount of energy. The endothermic case has as a lowerlimit the completely inelastic collision, in which both particles moveat the centre of mass velocity. On the other hand, the exothermiccase has no upper limit - the amount of energy injected into thesystem in the exothermic case is determined by the mass splitting δ between the two states. We perform five high resolution DM-only simulations of a singleMW-size halo in a cosmological context. In addition to simulatingthe halo with the inelastic SIDM model, simulations with an elas-tic SIDM model as well as standard CDM are also performed forcomparison. For each SIDM model, we examine cases with the twoprimordial excited fractions χ = 100% , and χ = 50% .The elastic SIDM model we examine is obtained by neglectingenergy changes for the inelastic up and down-scattering reactionsi.e. keeping A ij = 1 in Equation 4. The resulting elastic SIDMmodel is otherwise identical to inelastic SIDM, which enables us toisolate the effects of energy injection into the halo.Our simulations used cosmological parameters consistent withPlanck ( Ω m = 0.302, Ω Λ = 0.698, Ω b = 0.046, h = 0.69, σ = 0.839, n s = 0.967, Planck Collaboration et al. 2014; Spergel et al. 2015).The gravitational softening length is fixed in comoving coordinatesuntil z = 9 and then fixed in physical units until z = 0 , resultingin a Plummer-equivalent softening length of 72.4 pc at z = 0 . The MNRAS , 1–16 (2020)
K.T.E. Chua et. al. dark matter particles have a mass resolution of . × M (cid:12) .Haloes are identified using a friends-of-friends (fof) algorithm witha linking length of 0.2 (Davis et al. 1985), and the subfind algorithmis subsequently used to identify gravitationally self-bound subhaloes(Springel et al. 2001; Dolag et al. 2009). In the CDM simulation,the virial mass of the halo at z = 0 is M = 1 . × M (cid:12) andthe virial radius is R = 243 kpc . In this section, we focus on the effect of inelastic SIDM on thestructure of the dark matter halo at z = 0 , particularly in terms ofthe density profile, velocity profile, phase-space structure, and haloshape. We begin by examining the spherically-averaged DM density profile ρ ( r ) of the MW-size halo. The particles are binned into 30 shellsspaced in logarithmic intervals in the range 1 kpc < r < 300 kpc. Theresults are shown in the upper panels of Fig. 1, with black, blue,and red lines representing the CDM, elastic SIDM, and inelasticSIDM models, respectively. The columns distinguish between theprimordial excitation fractions χ = 100% (left panels), and χ = 50% (right panels). In addition to the total mass density(solid lines), we also show the individual profiles for the ground(dotted lines) and excited states (dashed lines). To avoid shells withlow number counts, shells containing less than 20 particles areignored. This affects primarily the inelastic model with χ =100% (left panel, dotted red lines). We note that the density profileswere previously also presented in Fig. 8 of V19.The upper row of Fig. 1 shows that both elastic and inelasticSIDM models produce lower densities and extended cores near thehalo centre compared to CDM. At r = 1 kpc, elastic SIDM reducesthe total density by a factor of ∼
10 relative to CDM. Defining thecore size R core as the radius below which the SIDM density profilesdeviate from CDM, we estimate the elastic SIDM haloes to havedensity cores of size R core ∼ kpc.Inelastic self-interactions further increase the core size anddecrease the central density compared to elastic SIDM. For the χ = 100% excitation fraction (upper left panel), the total den-sity at r = 1 kpc is approximately times smaller than that ofthe CDM halo, with a core size of R core ∼ kpc. Furthermore,the ground state densities (dotted lines) exhibit a large differencebetween the inelastic and elastic SIDM models: the ground statedensity of the inelastic halo is significantly lower compared thanthat of the elastic counterpart, with a difference of approximatelytwo orders of magnitude near the halo centre. For the lower excita-tion fraction χ = 50% , the ground state density is only slightlyabove that of the excited state, since the de-excited particles onlyconstitute a small fraction of the total ground state population.For each SIDM case, we further plot the ratio of the SIDMprofiles compared to CDM ( ρ SIDM /ρ CDM ) in the narrow plot at thebottom of each panel. At intermediate radii (10 kpc (cid:38) r (cid:38)
100 kpc),we note that the density of the inelastic SIDM haloes actually exceedthat of the CDM counterpart. This is due to the redistribution of R is the radius within which the average density is ¯ ρ = 200 ρ crit , where ρ crit is the critical density of the universe. Other properties for the CDM andSIDM haloes can be found in Table 1 of V19. particles from the inner regions during the elastic scattering process.However, the particles remain bound to the halo in the absence ofan energy-injection mechanism. In contrast, the density profile ofthe inelastic χ = 100% case does not exceed that of CDM atthe same radii. Here, the down-scattering process imparts velocitykicks, which enable the de-excited ground state particles to travelfurther and even escape from the inelastic halo completely. As aresult, in addition to a redistribution of particles, there is a removalof particles from the halo ( halo evaporation ) which decreases itsoverall density and equivalently, mass. These velocity kicks injectenergy equivalent to hundreds of millions of SNII into the halo(V19), which further explains the more efficient core formationwith the inelastic SIDM model. The slight elevation of the inelasticdensity profile at r ≈ kpc compared to CDM, also reflectsdown-scattered particles on their way out towards the outskirts ofthe halo.The bottom panels of Fig. 1 show the slopes Γ( r ) ≡− d ln ρ/d ln r of the logarithmic density profiles. The slopes arecalculated from the density profiles ρ ( r ) using a central differencescheme, such that for the n -th bin, we have Γ( r n ) = − ln( ρ ( r n +1 )) − ln( ρ ( r n − ))ln( r n +1 ) − ln( r n − ) . (6)For the first and last radial bins, forward and backward schemesare used respectively to estimate the slopes. In the CDM case, theNavarro-Frenk-White (NFW) profile (Navarro et al. 1996), with acusp ( Γ ≈ ) at small radii and Γ ≈ near the virial radius, is agood fit to the halo density profile. On the other hand, in all SIDMmodels, we find that at small radii (1 kpc < r (cid:46) Γ ≈ , revealing the presence of density cores. The overallslopes of the inelastic models beyond the density cores remain belowthat of the elastic model, up to intermediate radii of r (cid:46) kpc.This effect is stronger and persists to a larger radius for the largerprimordial excitation fraction.Another key difference between the elastic and inelastic modelslies in the slopes of the ground state particles for χ = 100% (lower left panel). In the elastic model, the ground state densityprofile is steeper (larger Γ ) compared to the excited state beyondthe density core. The opposite trend occurs for the inelastic model,where the ground state density profile is shallower (smaller Γ ) thanthat of the excited state.In general, it is clear that the impact of inelastic SIDM dependsstrongly on the initial fraction of excited particles: a higher primor-dial excited fraction results in larger modifications to the densityprofile, a result of the larger amount of energy injected into thehalo. The velocity dispersion and anisotropy are useful in providing in-sights into the orbital structure of haloes since the collision of darkmatter particles can affect their velocity distributions due to the ther-malisation of the halo (Vogelsberger & Zavala 2013; Brinckmannet al. 2018). We calculate the velocity dispersion σ as σ ( r ) = (cid:10) | v − ¯ v | (cid:11) = 1 N p N p (cid:88) i =1 | v i − ¯ v | (7)where N p is the number of particles and ¯ v is the mean velocity ofparticles in a given spherical shell. The velocity anisotropy param-eter is defined as β ( r ) = 1 − σ t ( r )2 σ r ( r ) (8) MNRAS000
100 kpc),we note that the density of the inelastic SIDM haloes actually exceedthat of the CDM counterpart. This is due to the redistribution of R is the radius within which the average density is ¯ ρ = 200 ρ crit , where ρ crit is the critical density of the universe. Other properties for the CDM andSIDM haloes can be found in Table 1 of V19. particles from the inner regions during the elastic scattering process.However, the particles remain bound to the halo in the absence ofan energy-injection mechanism. In contrast, the density profile ofthe inelastic χ = 100% case does not exceed that of CDM atthe same radii. Here, the down-scattering process imparts velocitykicks, which enable the de-excited ground state particles to travelfurther and even escape from the inelastic halo completely. As aresult, in addition to a redistribution of particles, there is a removalof particles from the halo ( halo evaporation ) which decreases itsoverall density and equivalently, mass. These velocity kicks injectenergy equivalent to hundreds of millions of SNII into the halo(V19), which further explains the more efficient core formationwith the inelastic SIDM model. The slight elevation of the inelasticdensity profile at r ≈ kpc compared to CDM, also reflectsdown-scattered particles on their way out towards the outskirts ofthe halo.The bottom panels of Fig. 1 show the slopes Γ( r ) ≡− d ln ρ/d ln r of the logarithmic density profiles. The slopes arecalculated from the density profiles ρ ( r ) using a central differencescheme, such that for the n -th bin, we have Γ( r n ) = − ln( ρ ( r n +1 )) − ln( ρ ( r n − ))ln( r n +1 ) − ln( r n − ) . (6)For the first and last radial bins, forward and backward schemesare used respectively to estimate the slopes. In the CDM case, theNavarro-Frenk-White (NFW) profile (Navarro et al. 1996), with acusp ( Γ ≈ ) at small radii and Γ ≈ near the virial radius, is agood fit to the halo density profile. On the other hand, in all SIDMmodels, we find that at small radii (1 kpc < r (cid:46) Γ ≈ , revealing the presence of density cores. The overallslopes of the inelastic models beyond the density cores remain belowthat of the elastic model, up to intermediate radii of r (cid:46) kpc.This effect is stronger and persists to a larger radius for the largerprimordial excitation fraction.Another key difference between the elastic and inelastic modelslies in the slopes of the ground state particles for χ = 100% (lower left panel). In the elastic model, the ground state densityprofile is steeper (larger Γ ) compared to the excited state beyondthe density core. The opposite trend occurs for the inelastic model,where the ground state density profile is shallower (smaller Γ ) thanthat of the excited state.In general, it is clear that the impact of inelastic SIDM dependsstrongly on the initial fraction of excited particles: a higher primor-dial excited fraction results in larger modifications to the densityprofile, a result of the larger amount of energy injected into thehalo. The velocity dispersion and anisotropy are useful in providing in-sights into the orbital structure of haloes since the collision of darkmatter particles can affect their velocity distributions due to the ther-malisation of the halo (Vogelsberger & Zavala 2013; Brinckmannet al. 2018). We calculate the velocity dispersion σ as σ ( r ) = (cid:10) | v − ¯ v | (cid:11) = 1 N p N p (cid:88) i =1 | v i − ¯ v | (7)where N p is the number of particles and ¯ v is the mean velocity ofparticles in a given spherical shell. The velocity anisotropy param-eter is defined as β ( r ) = 1 − σ t ( r )2 σ r ( r ) (8) MNRAS000 , 1–16 (2020) mpact of inelastic SIDM on a MW halo D e n s i t y ρ [ M fl k p c − ] χ = 100% CDMelastic ( χ = 100% ) inelastic ( χ = 100% ) totalgroundexcited -1 r a t i o χ = 50% CDMelastic ( χ = 50% ) inelastic ( χ = 50% ) totalgroundexcited L o g . s l o p e Γ ( r ) χ = 100% r [kpc] r a t i o χ = 50% r [kpc] Figure 1. Spherically averaged radial density profile ( top panels ) and its logarithmic slope ( bottom panels ) for the Milky Way-size halo.
The columnsare for different primordial fractions of the excited states, χ = 100% (left) and χ = 50% (right). We show the results for CDM (black), elastic SIDM(blue) and inelastic SIDM (red), with the narrow plot at the bottom of each panel showing the ratio of each quantity relative to CDM. The dotted (dashed)lines show the profiles for only the ground (excited) state. Top panels : In the central regions, both SIDM models lead to the formation of density cores, withcentral densities depleted compared to CDM. Inelastic SIDM produces larger and lower density cores than elastic SIDM, since the exothermic down-scatteringreactions lead to an expulsion of ground state particles from the halo. At intermediate radii, the elastic density profile exceeds that of CDM while the inelasticprofile almost never does so. This is because particles are only redistributed from the central regions to the intermediate regions in the elastic case and arenot removed as in the inelastic case. For the excitation fraction χ = 100% , the removal of down-scattered ground state particles leads to a significantsuppression of the ground state inelastic density profile (blue dotted lines) relative to elastic SIDM. Bottom panels : Radial profiles of the logarithmic densityslope Γ( r ) = − d ln ρ/d ln r . In all SIDM cases, the slope reaches a value of zero, indicating clearly the presence of density cores. where σ t and σ r are the tangential and radial velocity dispersionsaveraged over spherical shells, calculated similarly to Equation 7.Haloes with β = 0 are considered isotropic, while β > and β < correspond to radially and tangentially biased haloes respectively.We ignore substructure in calculating the velocity anisotropy inorder to concentrate on the smooth component of the halo.Fig. 2 shows both the radial profiles of the velocity disper-sion σ ( r ) (top row) and the velocity anisotropy parameter β ( r ) (bottom row), with the columns distinguishing between the two pri-mordial excitation fractions. The velocity dispersion of the CDMhalo exhibits the well-known temperature inversion in the centralregion, where the velocity dispersion drops towards the centre, con-sistent with previous CDM simulations (e.g. Navarro et al. 2010;Tissera et al. 2010). By transporting energy into the centre, self-interactions substantially affect the dark matter velocities, erasingthe temperature inversion in the SIDM models. Thus, the SIDM MNRAS , 1–16 (2020)
K.T.E. Chua et. al. V e l d i s p l og σ [ ( k m / s ) ] χ = 100% CDMelastic ( χ = 100% ) inelastic ( χ = 100% ) totalgroundexcited r a t i o χ = 50% CDMelastic ( χ = 50% ) inelastic ( χ = 50% ) totalgroundexcited V e l a n i s o β ≡ − σ t / ( σ r ) χ = 100% r [kpc] r a t i o χ = 50% r [kpc] Figure 2. Velocity dispersion (top panels) and velocity anisotropy profiles (bottom panels) for the Milky Way-size halo . The narrow plot at the bottom ofeach panel shows the ratio of the SIDM quantity relative to CDM.
Top panels : The SIDM velocity dispersion profiles are flat (isothermal) in the inner regions,exceeding that of CDM which exhibits the well-known temperature inversion. Inelastic scattering results in a suppression of the velocity dispersion comparedto elastic SIDM, a difference which is most apparent with a primordial excited fraction of χ = 100% . Bottom panels : The velocity anisotropy β of theSIDM haloes is closer to β = 0 (more isotropic) compared to the CDM counterpart. Although the difference between SIDM and CDM decreases in generaltowards the virial radius, the inelastic SIDM halo with excitation fraction χ = 100% remains significantly more isotropic than the CDM halo even at thevirial radius. velocity dispersions flatten and become fairly constant (isothermal)towards the inner halo. Inelastic down-scattering further reducesthe velocity dispersions compared to the elastic counterparts, sincekinetic energy is carried away by escaping down-scattered parti-cles. At intermediate radii, the velocity dispersion of the inelastic χ = 100% halo (red curve, left panel) is ≈
10 per cent lower thanthe CDM counterpart, a difference which persists for a substantialfraction of the halo, up to r ≈ kpc.An additional impact of inelastic self-interaction can also beobserved from the ground state velocity dispersions. For the higherprimordial excitation fraction χ = 100% (Fig. 7, left panel), velocity kicks received by the ground state particles in the inelasticSIDM model increase the inelastic velocity dispersion of the groundstate (dotted line) over that of the excited state particles (dashed line).Conversely, for elastic SIDM, the velocity dispersion of ground stateparticles is lower than of the excited state. In this case, ground stateparticles are only ever the products of elastic interactions, and aretherefore more limited in their allowed velocities.The bottom panels of Fig. 2 present the radial profiles of thevelocity anisotropy parameter β ( r ) . In general, self-interactionsisotropize the orbits ( β closer to zero) compared to CDM, especiallyin the inner and intermediate regions of the halo. Near the halo cen- MNRAS000
10 per cent lower thanthe CDM counterpart, a difference which persists for a substantialfraction of the halo, up to r ≈ kpc.An additional impact of inelastic self-interaction can also beobserved from the ground state velocity dispersions. For the higherprimordial excitation fraction χ = 100% (Fig. 7, left panel), velocity kicks received by the ground state particles in the inelasticSIDM model increase the inelastic velocity dispersion of the groundstate (dotted line) over that of the excited state particles (dashed line).Conversely, for elastic SIDM, the velocity dispersion of ground stateparticles is lower than of the excited state. In this case, ground stateparticles are only ever the products of elastic interactions, and aretherefore more limited in their allowed velocities.The bottom panels of Fig. 2 present the radial profiles of thevelocity anisotropy parameter β ( r ) . In general, self-interactionsisotropize the orbits ( β closer to zero) compared to CDM, especiallyin the inner and intermediate regions of the halo. Near the halo cen- MNRAS000 , 1–16 (2020) mpact of inelastic SIDM on a MW halo Log. slope d ln ρ/ d ln r V e l o c i t y a n i s o t r o p y β χ = 100% totalgroundexcitedCDMelastic ( χ = 100% ) inelastic ( χ = 100% )Hansen & Moore (2006) Log. slope d ln ρ/ d ln r χ = 50% totalgroundexcitedCDMelastic ( χ = 50% ) inelastic ( χ = 50% )Hansen & Moore (2006) Figure 3. Velocity anisotropy parameter vs logarithmic slope of the density profile for the Milky Way-size halo.
The dashed black line represents thelinear relation proposed by Hansen & Moore (2006).The CDM halo agrees with the HM06 relation whereas the SIDM haloes diverge from this linear relation.The difference is particularly significant for d ln ρ/d ln r (cid:38) − , which corresponds to the core region in the SIDM haloes. In this region, a simple extrapolationof the CDM result leads to an under-prediction of the velocity anisotropy since the flattening of the density profile is not taken into account. tre, β approaches zero more rapidly in the SIDM models than inCDM, due to collisional relaxation driving the core region isother-mal (Colín et al. 2002). Whereas elastic SIDM differs from CDMonly in the inner and intermediate regions, the effects of inelasticSIDM persist up to and even beyond the virial radius, especiallyfor the χ = 100% excitation fraction. This result confirms onceagain that energy injection from inelastic down-scattering plays animportant role in the structure of the halo far from the halo centre.The effect of SIDM on halo structure is also evident whenplotting the anisotropy parameter β against the logarithmic slope ofthe density profile ( d ln ρ/d ln r ≡ − Γ ), shown in Fig. 3. A relationbetween these two quantities were first identified by Hansen &Moore (2006) (hereafter HM06) in CDM simulations. This relationis shown as black dashed lines in Fig. 3. Our CDM results (blackdots) is in good agreement with the HM06 relation, which predicts acentral velocity anisotropy and slope of ( β (0) , − Γ(0)) = (0 , − .In the following, we focus on results that include both DMstates (i.e. the filled circles). For steep slopes ( d ln ρ/d ln r (cid:46) − ),the SIDM results are close to, but lie slightly below the HM06relation. This is due to the increased velocity isotropy (smaller β ) at intermediate radii in the SIDM runs. In the core region,( d ln ρ/d ln r (cid:38) − ), a simple extrapolation of the HM06 rela-tion under-predicts the velocity anisotropy compared to the SIDMrun since the CDM results do not account for the flattening of thedensity profile. The Contrary to CDM, both SIDM models predicta central velocity anisotropy and slope of ( β (0) , − Γ(0)) ≈ (0 , ,regardless of the primordial fraction of excited states. In CDM N -body simulations, Taylor & Navarro (2001) discoveredthat the dark matter density and velocity dispersion of any halo can be combined to form a pseudo phase-space density Q ( r ) ≡ ρ ( r ) /σ ( r ) , a quantity which is inversely proportional to the localentropy. Q ( r ) is essentially universal for all haloes, and can beapproximated by a power law Q ( r ) ∝ r − α . The value of α obtained( α ≈ . ) was also predicted by the analytic spherical infallsolution of Bertschinger (1985).In Fig. 4, we present the pseudo phase-space density pro-files Q ( r ) of the main halo for the CDM and SIDM models. TheCDM case obeys and confirms the power law seen in previous N -body simulations (e.g. Taylor & Navarro 2001; Ludlow et al. 2011).With SIDM however, the power-law behaviour is broken since self-scattering leads to a flattening in Q ( r ) for r < kpc. This isdue to the central flattening of the individual density and veloc-ity profiles by self-interactions discussed in the previous sections.Inelastic interactions further modify the pseudo phase-space pro-file, resulting in a decrease in Q ( r ) in the inner halo compared tothe elastic case, most notably for the higher primordial excitationfraction χ = 100% .To further illustrate the differences in the SIDM models, weplot the radial phase-space distributions in Fig. 5. Here, the pixelcolour represents the particle count in each ( r, v r ) bin, with darkpoints representing low counts, enhancements representing highparticle counts and uncoloured portions representing unoccupiedregions of phase-space. For ground state particles, we distinguishbetween particles bound to the central subhalo (leftmost column,as identified by subfind), and all particles within the fof group(middle column), which includes both unbound particles as well assubstructure. In Fig. 5, CDM particles are considered as ground stateparticles. For the SIDM models, we include also the phase-spacedistribution of all excited state particles in the fof group (rightmostcolumn).We first note that for a primordial excitation fraction of MNRAS , 1–16 (2020)
K.T.E. Chua et. al. -6 -5 -4 -3 -2 -1 ρ / σ [ M fl s k p c − k m − ] ∝ r − . χ = 100% CDMelastic ( χ = 100% ) inelastic ( χ = 100% ) totalgroundexcited χ = 50% CDMelastic ( χ = 50% ) inelastic ( χ = 50% ) totalgroundexcited r [kpc] -2 -1 r a t i o r [kpc] Figure 4. Pseudo phase-space density as a function of radius for the Milky Way-size halo.
Whereas the CDM profile agrees with the power law profile( ρ/σ ∝ r − . ) found by Taylor & Navarro (2001), both elastic and inelastic SIDM lead to a flattening of the pseudo phase-space density profile towardsthe centre, due to a combination of core formation and flattening of the velocity dispersion. χ = 100% , the phase-space distributions show a clear deficit ofthe ground state abundances in the inelastic SIDM model (third row)compared to the elastic model (second row), which agrees with ourprevious results for the halo density profiles (Fig. 1). In the radialphase-space distributions, satellite subhaloes appear as enhance-ments in the number counts localised at particular radii and radialvelocities. The deficiency of ground state particles in the inelastic χ = 100% case carries over to the subhalo structure, whichis noticeably absent in the ground state phase-space histogram ofthe fof halo (third row, middle column). Due to the small poten-tial wells of low-mass satellite subhaloes, the ground state particlesformed from exothermic down-scattering can easily escape due tothe imparted velocity kicks. In the inelastic model, we find most ofthe resolved subhaloes do not contain any ground state particles:only 49 out of 3304 resolved subhaloes have a non-zero ground statefraction. In comparison, most of the subhaloes (4298 out of 4306)in the elastic model with the same primordial excitation fractionhost ground state particles. It is also important to note that sub-haloes might not correspond exactly between the various runs fortwo reasons: (i) inelastic self-interactions lead to the evaporation ofsubhaloes and thus a reduction in subhalo abundance (V19), and(ii) the temporal evolution of the counterparts deviate progressivelyas time goes by, since subhalo accretion histories and orbits areaffected by dynamical changes generated by self-interactions.For both inelastic SIDM cases, Fig. 5 shows a population ofground state particles moving outwards with high radial velocities of We consider a subhalo to be resolved if it contains at least 100 particlesin total. Model mean v median v [km s − ] [km s − ]CDM 187.3 179.4elastic SIDM ( χ = 100% ) 185.9 176.0elastic SIDM ( χ = 50% ) 183.3 172.7inelastic SIDM ( χ = 100% ) 168.5 159.4inelastic SIDM ( χ = 50% ) 177.1 166.7 Table 1.
Mean and median of the (total) DM velocity distributions shownin Fig. 6. up to ≈ km s − . These high-speed ground state particles, notice-ably absent in CDM and the elastic SIDM model, reflect the velocitykicks they received during inelastic down-scattering. Although themajority of these kicks result in outward-moving ( v r > ) parti-cles, they can occasionally result in inward-moving ( v r < ) par-ticles, as evidenced from the smaller population of particles withhigh-speed negative radial velocities in the inelastic ground statephase-space distributions (third and bottom rows). For the inelastic χ = 100% case, we find that 80 per cent of these inward mov-ing particles originate from accreted subhaloes, and the remainingone-fifth from the central halo. The local distribution of DM velocities is important for DM di-rect detection experiments, which depend on the shape of the localvelocity distribution f ( v ) . We measure f ( v ) within cubes of side-length 2 kpc, sampled at a distance of r = 8 kpc (the solar circle) MNRAS , 1–16 (2020) mpact of inelastic SIDM on a MW halo v r [ k m / s ] CDM
Ground (Central) Ground (FOF)5000500 v r [ k m / s ] Elastic ( χ = 100% ) Excited (FOF)5000500 v r [ k m / s ] Inelastic ( χ = 100% ) v r [ k m / s ] Elastic ( χ = 50% ) r [kpc]5000500 v r [ k m / s ] Inelastic ( χ = 50% ) r [kpc] 0 200 400 r [kpc] Figure 5. Radial phase-space histograms of all models for the Milky Way-size halo.
The pixel colour represents the particle count in each ( r, v r ) bin, withdark points representing low counts, enhancements representing high particle counts and uncoloured portions representing unoccupied regions of phase-space.We distinguish between bound ground state particles in the main subhalo (leftmost column), all ground state particles within the fof halo (middle column), andexcited state particles within the fof halo (rightmost column). Results for the fof halo consider all substructures as well as unbound particles. Inelastic SIDM(third and fifth rows) results in a population of ground state particles with radial speeds up to 500 km s − throughout the halo. For a primordial excitationfraction of χ = 100% , the inelastic SIDM model strongly suppresses substructure in the ground state (third row, middle column).MNRAS , 1–16 (2020) K.T.E. Chua et. al. f ( v ) [ − k m s − ] CDMelastic ( χ = 100% ) inelastic ( χ = 100% ) totalgroundexcited CDMelastic ( χ = 50% ) inelastic ( χ = 50% ) totalgroundexcited v [km s − ]01 R a t i o v [km s − ] Figure 6. Top panels: Local dark matter velocity distribution for the Milky Way-size halo at the solar circle.
We measure f ( v ) within cubes of side-length2 kpc, sampled at a distance of r = 8 kpc for 1000 randomly selected observers. For each cube, we calculate the histogram in v of all particles (ground+ excited states), as well as for the ground and excited states separately. The curves show the median distribution, while the shading denotes the 25th-75thpercentile across the 1000 random observers. The dotted-dashed line in each panel represents a Maxwellian distribution with mean velocity 187.3 km s − . Atthe bottom of each panel, we show the ratio of each SIDM case relative to CDM. Compared to CDM, the SIDM models decrease the overall peak velocity ofthe distribution and suppress the fraction of particles with velocities v > km s − . In the SIDM models, self-interactions cause the velocity distribution tobecome closer to Maxwellian. In the inelastic model with primordial excitation fraction χ = 100% (left panel), the down-scattering of particles into theground state results in population of high-speed ground state particles with speeds v (cid:38) km s − . Such a feature is absent in the velocity distribution of theexcited particles and for the elastic SIDM model. for 1000 randomly selected observers . For each cube, we calculatethe histogram in v for all particles (ground + excited states), aswell as the individual particle states, normalising each histogramsuch that (cid:82) f ( v ) dv = 1 . The resulting velocity distributions areshown in Fig. 6, with lines representing the median distributionsover the 1000 random observers, and shaded regions denoting the25th-75th percentile for the total (ground + excited state) distri-bution. For comparison, we plot also a Maxwellian distribution(dotted-dashed curve) with mean speed equal to that of the CDMhalo (187.3 km s − ). Table 1 reports the mean and median speedsof the total distributions for the five cases.For CDM, our results are similar to previous work, where thevelocity distribution is close to Maxwellian, with features that canbe traced to the halo assembly history, as well as a tail of high-speedparticles in excess of the best-fit Maxwellian distribution (e.g. Vo-gelsberger et al. 2009; Kuhlen et al. 2010; Pillepich et al. 2014;Butsky et al. 2016). In the SIDM models, the total velocity distri-butions are noticeably shifted to lower speeds compared to CDM.This is related to the mass redistribution within the inner regions ofthe halo, especially in the core region (see Fig. 2). The largest sup- The cubes are oriented aligned with the simulation box. Changing theorientation e.g. by rotating the cubes is not found to affect the overall result. pression occurs for the inelastic SIDM model with χ = 100% :compared to CDM, the mean speed decreases from 187.3 km s − to168.5 km s − , while the median speed decreases from 179.4 km s − to 159.4 km s − (see Table 1). Smaller changes are observed forelastic SIDM as well as for the cases with the smaller primordialexcited fraction.At higher velocities ( v (cid:38)
200 km s − ), the SIDM modelssuppress the high-speed tail and result in a more steeply fallingdistribution, similar to that observed in some hydrodynamic sim-ulations with baryonic physics (e.g. Pillepich et al. 2014; Butskyet al. 2016). While the elastic SIDM halo has a velocity distributionclose to Maxwellian at the solar radius, the inelastic halo distri-bution deviates more significantly from Maxwellian. This resultsfrom the lower number of scattering events in the inelastic SIDMmodel compared to the elastic model: within the inner 10 kpc, thereare twice (1.2 times) as many scattering events in the elastic modelrelative to the inelastic model with primordial excitation fraction χ = 100% (50%) (see the lower panels of Fig. 7 of V19).The distributions of individual particle states generally followthat of the overall distribution, with the exception of the inelasticmodel with primordial excited fraction χ = 100% (left panel).In this particular case, the model predicts a substantial fractionof ground state particles (red dotted line) with high speeds up to MNRAS000
200 km s − ), the SIDM modelssuppress the high-speed tail and result in a more steeply fallingdistribution, similar to that observed in some hydrodynamic sim-ulations with baryonic physics (e.g. Pillepich et al. 2014; Butskyet al. 2016). While the elastic SIDM halo has a velocity distributionclose to Maxwellian at the solar radius, the inelastic halo distri-bution deviates more significantly from Maxwellian. This resultsfrom the lower number of scattering events in the inelastic SIDMmodel compared to the elastic model: within the inner 10 kpc, thereare twice (1.2 times) as many scattering events in the elastic modelrelative to the inelastic model with primordial excitation fraction χ = 100% (50%) (see the lower panels of Fig. 7 of V19).The distributions of individual particle states generally followthat of the overall distribution, with the exception of the inelasticmodel with primordial excited fraction χ = 100% (left panel).In this particular case, the model predicts a substantial fractionof ground state particles (red dotted line) with high speeds up to MNRAS000 , 1–16 (2020) mpact of inelastic SIDM on a MW halo
50 0 50 r [kpc] χ = 100% Solid: InelasticDotted: Elastic χ = 100% Solid: InelasticDotted: Elastic Density [M fl kpc − ] Figure 7. Dark matter density of the MW-size halo simulated withinelastic SIDM for a primordial excitation fraction of χ = 100% , withina 20 kpc-thick slice. To better elucidate the shape of the DM density, weshow in solid lines three isodensity contours, which correspond to three halo-centric distances. For comparison, the corresponding isodensity contours forthe elastic SIDM halo with the same excitation fraction are also shown asdotted lines. The isodensity contours show that the inelastic halo is lessspherical than the elastic counterpart at these radii. v ∼ km s − at the solar radius. This result is coherent with theradial phase-space histogram (Fig. 5), where the presence of thesehigh-speed particles previously been noted. For the lower primordialexcited fraction χ = 50% , the phase-space histogram suggeststhe presence of similar high-speed ground state particles, but theseare suppressed in the velocity distribution since such particles onlyconstitute a tiny fraction of the ground state population (see Fig. 1). Hierarchical structure formation theory and CDM simulations pre-dict that dark matter haloes are triaxial, due to the anisotropic ac-cretion of matter during halo growth (Dubinski & Carlberg 1991;Warren et al. 1992; Bullock 2002; Jing & Suto 2002; Bailin &Steinmetz 2005; Allgood et al. 2006; Macciò et al. 2008). Sinceself-interactions isotropise the velocities of the dark matter par-ticles (see Section 3.2), we expect SIDM to have an impact onhalo shapes as well. The non-spherical nature of the dark matterhalo is illustrated in Fig. 7, which visualizes the dark matter den-sity of the inelastic SIDM halo with primordial excitation fraction χ = 100% , for r < kpc. We also include three isodensitycontours for the same model (solid lines), as well as for the elasticSIDM counterpart (dotted lines) to help elucidate the halo shape.To quantify and compare the departure from spherical symme-try, we calculate the halo shape profiles using an iterative algorithm(Allgood et al. 2006; Zemp et al. 2011), which we briefly describe as follows. We assume that the isodensity surfaces can be described byellipsoidal shells, and are interested in determining the halo shapeat a particular radius through the axis ratios b/a and c/a , where a > b > c are the lengths of the principal axes of the shell. For aset of particles with mass m k within a particular ellipsoidal shell,we define the components of the shape tensor as S ij = (cid:80) N p k =1 m k x ( i ) k x ( j ) k (cid:80) N p k =1 m k , (9)where N p is the number of particles within the shell and x ( i ) k refersto the i th component of the k th particle. Note that the shape tensor isdirectly proportional to the second moment of the mass distribution.For each iteration, we compute the eigenvalues and eigenvectors ofthe shape tensor S ij . The eigenvectors denote the orientation ofthe principal frame, while the eigenvalues ( λ > λ > λ ) givethe axis ratios q ≡ b/a = (cid:112) λ /λ and s ≡ c/a = (cid:112) λ /λ .In the next iteration, we deform the ellipsoidal shell using thenew values of q and s while keeping the length of the majoraxis constant. This means that only particles with elliptical radius r ell = (cid:112) x + y /q + z /s that fall within the bin width areused to calculate the new shape tensor. For a given elliptical radius r ell = R , we define the shell to be . R < r ell < . R , thusshells are allowed to overlap and particles can belong to multipleshells simultaneously. This procedure has been shown to producereliable estimates of the halo shape (e.g. Zemp et al. 2011; Brinck-mann et al. 2018; Chua et al. 2019), although it differs slightly fromother methods which keep the enclosed volume of the shell con-stant, or calculate shapes using the entire enclosed mass within theellipsoidal volume.We calculate the axis ratio profiles in 30 ellipsoidal shellsspaced logarithmically between 1 kpc < r <
300 kpc. To begin thealgorithm, we start with spherical shells i.e. q = s = 1 . The itera-tions are performed until convergence is obtained, when successivevalues of q and s differ by less than one per cent.Previous convergence studies have shown that convergence inthe DM halo shapes is more demanding that for the DM densityprofiles. Within a halo, κ ( r ) , the ratio of the two-body relaxationtime-scale to the circular orbit time-scale at the virial radius, can beexpressed as: κ ( r ) ≡ √ N ( r )ln N ( r ) (cid:20) ¯ ρ ( r ) ρ crit (cid:21) − / , (10)where N ( r ) is the number of particles enclosed within radius r and ¯ ρ ( r ) is the mean density within radius r (Power et al. 2003). Theanalysis of MW-size haloes from the Aquarius simulation (Springelet al. 2008; Vera-Ciro et al. 2011) as well as haloes from the Illustrisproject (Chua et al. 2019) have found that the convergence radius r conv , defined by κ ( r conv ) = 7 , gives a good indication of theminimum radius where the shape profiles remain reliable in N -body simulations. This choice of r conv is larger than the Powerradius r p , defined by κ ( r p ) = 1 , which is traditionally applied tohalo circular velocity profiles. In addition, we also consider onlyshells containing at least 1000 particles, which has been found tobe approximately the minimum number of particles required for theiterative shape algorithm to be reliable (e.g. Tenneti et al. 2014).We show the obtained shape profiles q ( r ) and s ( r ) from oursimulations in Fig. 8. For each profile, only the converged region( r > r conv ) is shown. From q and s (top and middle rows respec-tively), we find a clear separation between the SIDM cases andCDM for r (cid:46) kpc: all the SIDM cases result in dark matterhaloes which are more spherical (larger q and s ) than the CDM MNRAS , 1–16 (2020) K.T.E. Chua et. al. q ≡ b / a χ = 100% CDMelastic SIDMinelastic SIDM totalground ( χ )excited ( χ ) r a t i o χ = 50% s ≡ c / a χ = 100% r [kpc]12 r a t i o χ = 50% r [kpc] Figure 8. Halo shape profiles for the Milky Way-size halo , showing q ≡ b/a (top row) and s ≡ c/a (bottom) as a function of radius. The bottom plotsattached to each panel show the ratio of the SIDM profiles relative to CDM. All SIDM models lead to more spherical haloes (larger q and s ) than CDM for asignificant fraction of the halo. This increase in sphericity is strongest in the inner regions and decreases towards the virial radius. In the elastic models, the halowith χ = 100% (left panels) is more spherical than its counterpart with χ = 50% (right panels). The inelastic SIDM models lead to a less sphericalhalo compared to their elastic SIDM counterparts, especially for χ = 100% . model. For the primordial excited fraction χ = 50% , the elasticand inelastic SIDM models both have q ≈ . , s ≈ . at a radiusof r = 10 kpc, much higher than the values q = 0 . , s = 0 . forCDM. Self-interactions tend to isotropise particle orbits, makinghaloes more spherical. This effect is strongest near the centre anddecreases towards the viral radius where the SIDM shapes are verysimilar to that of CDM.Interestingly, the inelastic SIDM model results in haloes thatare less spherical than their elastic counterparts. For χ = 100% ,Fig. 8 shows that s = 0 . for the elastic halo, compared to s = 0 . for the inelastic halo at r = 10 kpc. Visually, this result can also beobserved by comparing the isodensity contours of the two modelsshown in Fig. 7. This separation in the halo shape profiles persistsup to 200 kpc, and provides an additional diagnostic which can potentially be used to distinguish between the elastic and inelasticSIDM models. The larger sphericity in the elastic SIDM case canbe explained by its larger number of scattering events compared tothe inelastic halo, (see Section 3.4 and Fig 7. of V19).For elastic SIDM, DM halo shapes have been previously stud-ied by Peter et al. (2013) and Brinckmann et al. (2018). Ourelastic results are qualitatively similar to the median results ofthe galaxy-size haloes studied by Peter et al. (2013), when com-pared to the largest self-interaction cross section they considered( σ/m = 1 cm g − ).For the MW, estimates of the shape of the inner halo ( r (cid:46) kpc) can be inferred observationally. For example, stellar kine-matics together with equilibrium modelling with the Jeans equa-tions has suggested that s = 0 . ± . (Loebman et al. 2012). MNRAS000
300 kpc. To begin thealgorithm, we start with spherical shells i.e. q = s = 1 . The itera-tions are performed until convergence is obtained, when successivevalues of q and s differ by less than one per cent.Previous convergence studies have shown that convergence inthe DM halo shapes is more demanding that for the DM densityprofiles. Within a halo, κ ( r ) , the ratio of the two-body relaxationtime-scale to the circular orbit time-scale at the virial radius, can beexpressed as: κ ( r ) ≡ √ N ( r )ln N ( r ) (cid:20) ¯ ρ ( r ) ρ crit (cid:21) − / , (10)where N ( r ) is the number of particles enclosed within radius r and ¯ ρ ( r ) is the mean density within radius r (Power et al. 2003). Theanalysis of MW-size haloes from the Aquarius simulation (Springelet al. 2008; Vera-Ciro et al. 2011) as well as haloes from the Illustrisproject (Chua et al. 2019) have found that the convergence radius r conv , defined by κ ( r conv ) = 7 , gives a good indication of theminimum radius where the shape profiles remain reliable in N -body simulations. This choice of r conv is larger than the Powerradius r p , defined by κ ( r p ) = 1 , which is traditionally applied tohalo circular velocity profiles. In addition, we also consider onlyshells containing at least 1000 particles, which has been found tobe approximately the minimum number of particles required for theiterative shape algorithm to be reliable (e.g. Tenneti et al. 2014).We show the obtained shape profiles q ( r ) and s ( r ) from oursimulations in Fig. 8. For each profile, only the converged region( r > r conv ) is shown. From q and s (top and middle rows respec-tively), we find a clear separation between the SIDM cases andCDM for r (cid:46) kpc: all the SIDM cases result in dark matterhaloes which are more spherical (larger q and s ) than the CDM MNRAS , 1–16 (2020) K.T.E. Chua et. al. q ≡ b / a χ = 100% CDMelastic SIDMinelastic SIDM totalground ( χ )excited ( χ ) r a t i o χ = 50% s ≡ c / a χ = 100% r [kpc]12 r a t i o χ = 50% r [kpc] Figure 8. Halo shape profiles for the Milky Way-size halo , showing q ≡ b/a (top row) and s ≡ c/a (bottom) as a function of radius. The bottom plotsattached to each panel show the ratio of the SIDM profiles relative to CDM. All SIDM models lead to more spherical haloes (larger q and s ) than CDM for asignificant fraction of the halo. This increase in sphericity is strongest in the inner regions and decreases towards the virial radius. In the elastic models, the halowith χ = 100% (left panels) is more spherical than its counterpart with χ = 50% (right panels). The inelastic SIDM models lead to a less sphericalhalo compared to their elastic SIDM counterparts, especially for χ = 100% . model. For the primordial excited fraction χ = 50% , the elasticand inelastic SIDM models both have q ≈ . , s ≈ . at a radiusof r = 10 kpc, much higher than the values q = 0 . , s = 0 . forCDM. Self-interactions tend to isotropise particle orbits, makinghaloes more spherical. This effect is strongest near the centre anddecreases towards the viral radius where the SIDM shapes are verysimilar to that of CDM.Interestingly, the inelastic SIDM model results in haloes thatare less spherical than their elastic counterparts. For χ = 100% ,Fig. 8 shows that s = 0 . for the elastic halo, compared to s = 0 . for the inelastic halo at r = 10 kpc. Visually, this result can also beobserved by comparing the isodensity contours of the two modelsshown in Fig. 7. This separation in the halo shape profiles persistsup to 200 kpc, and provides an additional diagnostic which can potentially be used to distinguish between the elastic and inelasticSIDM models. The larger sphericity in the elastic SIDM case canbe explained by its larger number of scattering events compared tothe inelastic halo, (see Section 3.4 and Fig 7. of V19).For elastic SIDM, DM halo shapes have been previously stud-ied by Peter et al. (2013) and Brinckmann et al. (2018). Ourelastic results are qualitatively similar to the median results ofthe galaxy-size haloes studied by Peter et al. (2013), when com-pared to the largest self-interaction cross section they considered( σ/m = 1 cm g − ).For the MW, estimates of the shape of the inner halo ( r (cid:46) kpc) can be inferred observationally. For example, stellar kine-matics together with equilibrium modelling with the Jeans equa-tions has suggested that s = 0 . ± . (Loebman et al. 2012). MNRAS000 , 1–16 (2020) mpact of inelastic SIDM on a MW halo Tidal stellar streams are also widely used to model the MW haloshape, providing estimates of: s > 0.7 (Ibata et al. 2001), s = 0.72(Law & Majewski 2010), s = 0.8 (Vera-Ciro & Helmi 2013), and s = ± So far, we have discussed the structural characteristics of inelasticSIDM haloes only at redshift z = 0 . It is also relevant to considerhow these characteristics change as the haloes evolve from highredshift to the present epoch.For a CDM halo, the density profile is well-described by anNFW profile (Navarro et al. 1996): ρ NFW ( r ) ρ crit = δ c ( r/r s ) (1 + r/r s ) , (11)where δ c is a characteristic density contrast and the scale ra-dius r s is the radius where the logarithmic density slope is d ln ρ/d ln r = − . Due to the presence of the constant densitycore in the SIDM haloes, the SIDM density profiles are not well-described by the NFW profile (see Section 3.1 and also Todoroki &Medvedev 2019b). Thus, we only calculate δ c and the scale radius r s, CDM for the CDM halo using a two-parameter least-squares fit.At redshifts z = 0 and z = 10 , the comoving scale radii for theCDM halo are r s, CDM = 24 . kpc and 14.8 kpc, respectively.We first examine the evolution of the average inner densitywithin . r s, CDM and . r s, CDM . This results in the two mea-sures: (i) ¯ ρ . ≡ ¯ ρ ( r < . r s, CDM ) , and (ii) ¯ ρ . ≡ ¯ ρ ( r < . r s, CDM ) . For each of the five models, we trace the main halofrom the final snapshot at z = 0 back to z = 10 , picking in eachsnapshot its progenitor based on mass and distance. Subsequently, ¯ ρ . and ¯ ρ . are calculated for each redshift using the scale radiusof the CDM counterpart. The standardisation of the radius (i.e. using r s, CDM ) at each redshift ensures that the inner densities calculatedcan be compared consistently across the models.Fig. 9 shows the evolution of ¯ ρ . and ¯ ρ . in physical unitsfor CDM and the SIDM models, with the bottom attached plotsshowing the ratio of the inner densities to that of the CDM model.For ¯ ρ . , the effect of dark matter self-interactions manifests withtime (decreasing redshift), and increasingly suppresses the innerdensity relative to CDM with time. Relative to the CDM halo, ¯ ρ . is suppressed by 15–45 per cent at z = 10 across the SIDMmodels. By z = 0 , this number has increased to 50–80 per cent. Atfixed primordial excitation fraction χ , the inelastic halo has thelowest inner density at each redshift. Thus, we conclude that energyinjection from inelastic down-scattering reduces the inner densityin a shorter timescale compared to the elastic scale.The effects of self-interactions are weaker on ¯ ρ . compared to ¯ ρ . , due to the larger enclosed mass within . r s, CDM . At z = 0 , ¯ ρ . is only suppressed by 15–55 per cent relative to CDM. At highredshifts, due to the smaller cumulative number of scattering events,only an inelastic model with high primordial excitation fractioninjects sufficient energy through inelastic interactions to evacuatedark matter from within . r s, CDM . As a result, only the inelastic SIDM model with χ = 100% (left panel, red curve) shows asubstantial decrease in ¯ ρ . relative to CDM (-14 per cent).Finally, we examine the inner slope of the logarithmic densityprofile, which provides more information about the growth of theconstant density core. Here, we measure the inner slope γ by fittingthe density profile in the region r < . r s, CDM to a power law ρ ( r ) ∝ r − γ . We emphasise again that at each redshift, the scaleradius for the CDM halo is used for consistency across the models.The resulting evolution of the inner slope γ with redshift is shown inFig. 10. In the CDM model, we find the inner slope is approximatelyconstant across redshift, with a value of γ ≈ . . This valueis consistent with the NFW profile , for which the local slope is Γ( r = 0) = 1 at the centre and Γ( r s ) = 2 at the scale radius.In the SIDM models, both elastic and inelastic self-interactionsalready flatten the density profiles at redshift z = 10 compared tothe CDM model. Between z = 10 and z = 0 , γ decreases from γ ≈ to γ ≈ . in the SIDM cases. As the halo grows, self-interactions cause the inner core to grow, which leads to a shallowerinner slope. At low redshifts, the inner slopes in the inelastic casesare slightly smaller (i.e. the density profiles are flatter) than theelastic cases, due to more efficient core formation associated withinelastic SIDM.Comparing these inner slopes with the central densities fromFig. 9, we note that at high redshifts, elastic self-interactions alonehave only a small effect on the average inner densities ¯ ρ . and ¯ ρ . relative to CDM. However, the effect of elastic and inelasticself-interactions on the inner slope is already considerable even at z = 10 . We have examined the results of high resolution simulations ofa Milky Way-size dark matter halo composed of inelastic self-interacting dark matter. Our inelastic SIDM model consists of nearlydegenerate two-state dark matter particles with an energy level split-ting of δ = 10 keV, which we implemented and simulated usingArepo. During the exothermic down-scattering process, the modelresults in velocity kicks of approximately 424 km s − . To under-stand the effects of such energy injection, we have examined andcompared inelastic SIDM simulations of a MW-size halo against anelastic SIDM model and the conventional CDM model. The elasticSIDM model simulated in this work suppresses the energy changeduring down-scattering and up-scattering but is otherwise identicalto inelastic SIDM. In addition, we have also distinguished betweenthe configurations where the simulation begins with all DM particlesin the excited state ( χ = 100% ) and where only half begin in theexcited state ( χ = 50% ). Using these five cases, we examinedthe effects of elastic and inelastic SIDM on the internal structureand assembly of dark matter haloes. In the following, we summariseour main findings, concentrating on the case with primordial excitedfraction χ = 100% :(i) Energy injection resulting from inelastic self-interactions re-duces the central density of the inelastic SIDM halo in the innerregions of the halo and results in a larger core compared to theelastic counterpart (Fig. 1). In the χ = 100% configuration,inelastic SIDM reduces the density by a factor of approximately 20relative to CDM, compared to a factor of 10 for elastic SIDM, at a For the NFW profile, the local slope of the logarithmic density profile is Γ( r ) ≡ − d ln ρ/d ln r = 1 + 2 r/ ( r + r s ) .MNRAS , 1–16 (2020) K.T.E. Chua et. al. a v e r a g e i nn e r d e n s i t y ¯ ρ [ M fl k p c − ] χ = 100% CDMelastic SIDMinelastic SIDM ¯ ρ . ≡ ¯ ρ ( r < . r s, CDM )¯ ρ . ≡ ¯ ρ ( r < . r s, CDM ) χ = 50%
10 8 6 4 3 2 1 0
Redshift z r a t i o
10 8 6 4 3 2 1 0
Redshift z Figure 9. Average inner density (in physical units) as a function of redshift for the MW-like halo , within . r s, CDM and . r s, CDM . At each redshift,the scale radius of the CDM halo ( r s, CDM ) is used across the models for consistency. The bottom plots attached to each panel show the ratio of the SIDMquantities relative to CDM. For a fixed primordial excitation fraction, the inner densities are decreased in a shorter timescale in the inelastic SIDM modelcompared to the elastic model, thus the inelastic halo has the lowest density at each redshift. At z = 10 , ¯ ρ . is suppressed relative to CDM only in the inelasticmodel with χ = 100% , due to the larger enclosed mass within . r s, CDM ,
10 8 6 4 3 2 1 0
Redshift z C e n t r a l s l o p e γ ( r < . r s , C D M ) CDMelastic ( χ = 100% ) inelastic ( χ = 100% )
10 8 6 4 3 2 1 0
Redshift z CDMelastic ( χ = 50% ) inelastic ( χ = 50% ) Figure 10. Inner slope as a function of redshift for the MW-like halo . At each snapshot, the inner slope γ is obtained from fitting the density profile ρ ( r ) within 1/4 of the NFW scale radius of the CDM halo counterpart to a power law r − γ . For the CDM case, γ is approximately constant across redshift, remainingat γ ≈ . . In all the SIDM cases, self-interactions already result in a flattened density profile at redshift z = 10 compared to the CDM counterpart. Between z = 10 and z = 0 , the inner slope becomes less steep due to the increasing impact of self-interactions, with γ decreasing roughly from 1 to 0.25 in all theSIDM cases. At low redshifts, the inelastic halo is slightly less steep than the elastic counterpart, due to the larger core formation efficiency in the former.MNRAS000
Redshift z CDMelastic ( χ = 50% ) inelastic ( χ = 50% ) Figure 10. Inner slope as a function of redshift for the MW-like halo . At each snapshot, the inner slope γ is obtained from fitting the density profile ρ ( r ) within 1/4 of the NFW scale radius of the CDM halo counterpart to a power law r − γ . For the CDM case, γ is approximately constant across redshift, remainingat γ ≈ . . In all the SIDM cases, self-interactions already result in a flattened density profile at redshift z = 10 compared to the CDM counterpart. Between z = 10 and z = 0 , the inner slope becomes less steep due to the increasing impact of self-interactions, with γ decreasing roughly from 1 to 0.25 in all theSIDM cases. At low redshifts, the inelastic halo is slightly less steep than the elastic counterpart, due to the larger core formation efficiency in the former.MNRAS000 , 1–16 (2020) mpact of inelastic SIDM on a MW halo radius r = 1 kpc. Although the density profiles in the inner 3 kpcare flat for both elastic and inelastic SIDM, the logarithmic slopein the elastic models approaches that of CDM more rapidly withincreasing radius.(ii) At intermediate radii, the inelastic SIDM radial density pro-files are suppressed relative to CDM, because the velocity kicksfrom inelastic down-scattering unbind and eject particles from theMW halo. In contrast, the density profiles of the elastic SIDM mod-els are slightly raised compared to CDM since particles are solelytransferred from the inner regions to the intermediate regions ofthe halo. The density of the ground state particles in the inelasticmodel is around two orders of magnitude lower in the core regioncompared to the elastic model, reflecting the characteristic particleejection unique to inelastic SIDM.(iii) We found that self-interactions flatten the velocity dispersionin the inner regions (Fig. 2, top panels), causing the halo to becomeisothermal, in contrast to the “temperature inversion" observed inthe CDM case. Inelastic self-interactions leads to lower velocitydispersions compared to elastic interactions, due to the expulsion ofhigh-speed particles from the inelastic halo.(iv) Inelastic scattering results in an inelastic halo which is moreisotropic ( β ≈ ) than the CDM and elastic SIDM counterparts(Fig. 2, lower panels). Although elastic SIDM results in a similarisotropisation in the inner and intermediate regions of the halo, itdoes not modify the velocity anisotropy substantially in the outerregions. This implies the energy released in the exothermic interac-tions could be important in modifying particle orbits near the halooutskirts.(v) Self-interactions flatten the pseudo phase-phase density(Fig. 4) in the inner regions, resulting in a profile substantially differ-ent from the universal power law followed by CDM haloes. Furtheranalysis of the radial phase-space distributions (Fig. 5) reveals asubstantial population of ground state particles with radial veloci-ties of up to 500 km s − which is uniquely present in the inelasticmodels. These particles are unbound from the central subhalo andreflect the velocity kicks associated with inelastic down-scattering.(vi) The local velocity distribution f ( v ) of DM particles at thesolar circle (Fig. 6) shows that the inelastic SIDM models predict thepresence of high-speed ground state particles which have receivedvelocity kicks. The decreased core density in the inelastic modelsresults in a total f ( v ) that differs from both CDM and elastic SIDM.In the SIDM haloes, f ( v ) is shifted to lower velocities, together asuppression of the high-speed tail ( v (cid:38) km s − ). For inelasticSIDM, we also found that the ground state particles show a distinctpopulation of high-speed particles ( v (cid:38) km s − ). These parti-cles are not present in the velocity distribution of the excited stateparticles or in the elastic SIDM counterpart. The unique presenceof high-speed particles in velocity distribution f ( v ) of the inelasticSIDM halo is a potential signature for direct detection experiments.(vii) Dark matter self-interactions result in haloes that are morespherical than the CDM counterpart (Fig. 8), and are thus in bet-ter agreement with observational estimates of the MW halo shape.Interestingly, we found that our simulated halo was more spheri-cal with elastic SIDM than with inelastic SIDM. For example, theminor-to-major axis s ratio at r = 10 kpc is s = 0 . for the inelastichalo, compared to s = 0 . for the elastic halo and s = 0 . for theCDM halo. The larger number of scattering events in the elastic halodrives the shape to more more spherical compared to the inelasticcounterpart.(viii) Tracing the halo assembly history, we found that inelas-tic self-interactions reduce the inner density of the MW halo in ashorter timescale relative to the elastic scale (Fig. 9), and also result in a shallower inner slope at low redshifts (Fig. 10). As such, in-elastic self-interactions affect significantly the innermost region ofthe halo ( r (cid:46) . r s ) at high redshifts ( z (cid:38) ). We found that bothelastic and inelastic self-interactions already result in a substantialflattening of the central density profile, even at a redshift as high as z = 10 .For simulations carried out with the lower primordial excitationfraction ( χ = 50% ), the results are in general similar to thatobserved for χ = 100% . With a lower χ , the differencesrelative to CDM were more subdued because: (i) less energy isinjected into the halo when less particles start out in the excitedstate, and (ii) the down-scattered ground state particles form only asmall fraction of the total ground state particles in the halo. Whileenergy injection from inelastic self-interactions drives many of thedifferences between the inelastic and elastic haloes, the numberof scattering events appears to be the primary factor for certainsituations. For example, the latter factor is responsible for the localvelocity distribution f ( v ) of the elastic SIDM haloes being closer toMaxwellian than their inelastic counterparts, and also for the elastichaloes being more spherical.In conclusion, we have found that inelastic self-interactionscan significantly impact dark matter haloes, causing the structure ofinelastic SIDM haloes to exhibit distinct differences with respect totheir elastic SIDM and CDM counterparts. By performing simula-tions with only dark matter particles, we have focused on the effectsof elastic and inelastic scattering through comparisons to the CDMmodel. Although we explored two primordial excitation fractions,we remark that our results presented in this work correspond to asingle choice for the parameters δ , α , m χ and m φ in the Schutz& Slatyer (2015) model. In future work, we plan to expand thestudy by exploring different choices in the parameter space. At thesame time, it is well known that hydrodynamics simulations withCDM and baryons have predicted that baryonic physics associatedwith galaxy assembly e.g. gas cooling, stellar feedback, black holefeedback etc. has a relevant impact on the structure of dark matterhaloes. Future simulations incorporating both inelastic SIDM andgalaxy formation physics could prove to be helpful in understandingtheir combined effects on both the dark matter halo as well as theluminous galaxy. ACKNOWLEDGEMENTS
We thank David Barnes for his constructive comments on the paper.The simulations were performed on the joint MIT-Harvard comput-ing cluster supported by MKI and FAS. JZ acknowledges supportby a Grant of Excellence from the Icelandic Research Fund (grantnumber 173929).
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010,MNRAS, 407, 435Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H.,Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781MNRAS , 1–16 (2020) K.T.E. Chua et. al.
Arkani-Hamed N., Finkbeiner D. P., Slatyer T. R., Weiner N., 2009, Phys.Rev. D, 79, 015014Bailin J., Steinmetz M., 2005, ApJ, 627, 647Bertschinger E., 1985, ApJS, 58, 39Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93Bœhm C., Riazuelo A., Hansen S. H., Schaeffer R., 2002, Phys. Rev. D, 66,083505Bœhm C., Schewtschenko J. A., Wilkinson R. J., Baugh C. M., Pascoli S.,2014, MNRAS, 445, L31Bovy J., Bahmanyar A., Fritz T. K., Kallivayalil N., 2016, ApJ, 833, 31Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40Brinckmann T., Zavala J., Rapetti D., Hansen S. H., Vogelsberger M., 2018,MNRAS, 474, 746Brook C. B., Di Cintio A., 2015, MNRAS, 450, 3920Bryan S. E., Kay S. T., Duffy A. R., Schaye J., Dalla Vecchia C., BoothC. M., 2013, MNRAS, 429, 3316Buckley M. R., Zavala J., Cyr-Racine F.-Y., Sigurdson K., Vogelsberger M.,2014, Phys. Rev. D, 90, 043524Bullock J. S., 2002, in Natarajan P., ed., The Shapes of Galaxiesand their Dark Halos. pp 109–113 ( arXiv:astro-ph/0106380 ),doi:10.1142/9789812778017_0018Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343Butsky I., et al., 2016, MNRAS, 462, 663Chisari N. E., et al., 2017, MNRAS, 472, 1163Chua K. T. E., Pillepich A., Vogelsberger M., Hernquist L., 2019, MNRAS,484, 476Colín P., Avila-Reese V., Valenzuela O., 2000, ApJ, 542, 622Colín P., Avila-Reese V., Valenzuela O., Firmani C., 2002, ApJ, 581, 777Crain R. A., et al., 2015, MNRAS, 450, 1937Dalcanton J. J., Hogan C. J., 2001, ApJ, 561, 35Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497Dubinski J., 1994, ApJ, 431, 617Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., BoothC. M., 2010, MNRAS, 405, 2161Garzilli A., Ruchayskiy O., Magalich A., Boyarsky A., 2019, arXiv e-prints,p. arXiv:1912.09397Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16Governato F., et al., 2012, MNRAS, 422, 1231Hansen S. H., Moore B., 2006, New Astronomy, 11, 333Ibata R., Lewis G. F., Irwin M., Totten E., Quinn T., 2001, ApJ, 551, 294Iršič V., et al., 2017, Phys. Rev. D, 96, 023522Jing Y. P., Suto Y., 2002, ApJ, 574, 538Kamada A., Kaplinghat M., Pace A. B., Yu H.-B., 2017, Physical ReviewLetters, 119, 111102Katz N., Gunn J. E., 1991, ApJ, 377, 365Klypin A., Karachentsev I., Makarov D., Nasonova O., 2015, MNRAS, 454,1798Kuhlen M., Weiner N., Diemand J., Madau P., Moore B., Potter D., StadelJ., Zemp M., 2010, J. Cosmology Astropart. Phys., 2, 030Law D. R., Majewski S. R., 2010, ApJ, 714, 229Loebman S. R., Ivezić Ž., Quinn T. R., Governato F., Brooks A. M., Chris-tensen C. R., Jurić M., 2012, ApJ, 758, L23Ludlow A. D., Navarro J. F., Boylan-Kolchin M., Springel V., Jenkins A.,Frenk C. S., White S. D. M., 2011, MNRAS, 415, 3895Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940Macciò A. V., Paduroiu S., Anderhalden D., Schneider A., Moore B., 2012,MNRAS, 424, 1105Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P.,1999, ApJ, 524, L19Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, L72Navarro J. F., et al., 2010, MNRAS, 402, 21Oñ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,2092Oman K. A., et al., 2015, MNRAS, 452, 3650 Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS,430, 105Pillepich A., Kuhlen M., Guedes J., Madau P., 2014, ApJ, 784, 161Pillepich A., et al., 2018, MNRAS, 473, 4077Planck Collaboration et al., 2014, A&A, 571, A16Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., SpringelV., Stadel J., Quinn T., 2003, MNRAS, 338, 14Read J. I., Agertz O., Collins M. L. M., 2016, MNRAS, 459, 2573Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-KimmelS., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81Santos-Santos I. M. E., et al., 2020, MNRAS, 495, 58Sawala T., et al., 2016, MNRAS, 457, 1931Schutz K., Slatyer T. R., 2015, J. Cosmology Astropart. Phys., 1, 021Spergel D. N., Steinhardt P. J., 2000, Physical Review Letters, 84, 3760Spergel D. N., et al., 2003, ApJS, 148, 175Spergel D. N., Flauger R., Hložek R., 2015, Phys. Rev. D, 91, 023518Springel V., 2010, MNRAS, 401, 791Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS,328, 726Springel V., et al., 2005, Nature, 435, 629Springel V., et al., 2008, MNRAS, 391, 1685Taylor J. E., Navarro J. F., 2001, ApJ, 563, 483Tenneti A., Mandelbaum R., Di Matteo T., Feng Y., Khandai N., 2014,MNRAS, 441, 470Tissera P. B., White S. D. M., Pedrosa S., Scannapieco C., 2010, MNRAS,406, 922Todoroki K., Medvedev M. V., 2019a, MNRAS, 483, 3983Todoroki K., Medvedev M. V., 2019b, MNRAS, 483, 4004Todoroki K., Medvedev M. V., 2020, arXiv e-prints, p. arXiv:2003.11096Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1Vera-Ciro C., Helmi A., 2013, ApJ, 773, L4Vera-Ciro C. A., Sales L. V., Helmi A., Frenk C. S., Navarro J. F., SpringelV., Vogelsberger M., White S. D. M., 2011, MNRAS, 416, 1377Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, Phys. Rev. D,88, 043502Vogelsberger M., Zavala J., 2013, MNRAS, 430, 1722Vogelsberger M., et al., 2009, MNRAS, 395, 797Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740Vogelsberger M., et al., 2014, MNRAS, 444, 1518Vogelsberger M., Zavala J., Cyr-Racine F.-Y., Pfrommer C., Bringmann T.,Sigurdson K., 2016, MNRAS, 460, 1399Vogelsberger M., Zavala J., Schutz K., Slatyer T. R., 2019, MNRAS, 484,5437Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20Warren M. S., Quinn P. J., Salmon J. K., Zurek W. H., 1992, ApJ, 399, 405Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D.,Quataert E., 2016, ApJ, 827, L23Zavala J., Frenk C. S., 2019, Galaxies, 7, 81Zavala J., Jing Y. P., Faltenbacher A., Yepes G., Hoffman Y., Gottlöber S.,Catinella B., 2009, ApJ, 700, 1779Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20Zavala J., Lovell M. R., Vogelsberger M., Burger J. D., 2019, Phys. Rev. D,100, 063007Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2011, APJS, 197, 30Zolotov A., et al., 2012, ApJ, 761, 71de Blok W. J. G., McGaugh S. S., 1997, MNRAS, 290, 533This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000