The impact of baryonic disks on the shapes and profiles of self-interacting dark matter halos
We employ isolated N-body simulations to study the response of self-interacting dark matter (SIDM) halos in the presence of the baryonic potentials. Dark matter self-interactions lead to kinematic thermalization in the inner halo, resulting in a tight correlation between the dark matter and baryon distributions. A deep baryonic potential shortens the phase of SIDM core expansion and triggers core contraction. This effect can be further enhanced by a large self-scattering cross section. We find the final SIDM density profile is sensitive to the baryonic concentration and the strength of dark matter self-interactions. Assuming a spherical initial halo, we also study evolution of the SIDM halo shape together with the density profile. The halo shape at later epochs deviates from spherical symmetry due to the influence of the non-spherical disk potential, and its significance depends on the baryonic contribution to the total gravitational potential, relative to the dark matter one. In addition, we construct a multi-component model for the Milky Way, including an SIDM halo, a stellar disk and a bulge, and show it is consistent with observations from stellar kinematics and streams.
keywords:methods: numerical-galaxies: evolution-galaxies: formation-galaxies: structure-cosmology: theory.
The cold dark matter (CDM) model, where DM is assumed to be collisionless, is the leading theory for the growth of structure and formation of galaxies in our Universe. It fits the spectrum of matter fluctuations in the early universe with extraordinary precision (Planck Collaboration et al., 2014) and explains many important aspects of galaxy formation and evolution (Springel et al., 2006; Trujillo-Gomez et al., 2011; Frenk & White, 2012; Vogelsberger et al., 2014a, c). However, the success of CDM does not preclude the possibility that DM may have strong self-interactions. When DM particles have a scattering cross section per unit mass, , DM collisions occur multiple times in the inner halo over the cosmological timescale and make distinct departures from the CDM predictions, while the outer halo remains collisionless, retaining the large-scale predictions of CDM (Spergel & Steinhardt, 2000; Yoshida et al., 2000; Davé et al., 2001; Colín et al., 2002; Vogelsberger et al., 2012; Rocha et al., 2013). In fact, SIDM has been motivated to address outstanding discrepancies between observations on galactic scales and CDM predictions, see Bullock & Boylan-Kolchin (2017) for a recent review on the CDM challenges and Tulin & Yu (2017) on their solutions within SIDM. In particular, it has been shown that the diverse shapes of galactic rotation curves (Oman et al., 2015; Kuzio de Naray et al., 2010) can be explained naturally in the SIDM model (Kamada et al., 2017; Creasey et al., 2017), while being consistent with observations in galaxy clusters (Kaplinghat et al., 2016).
DM self-interactions kinematically thermalize the inner halo and lead to distinct features in the halo properties. For dwarf galaxies, where DM dominates over all radii, SIDM thermalization leads to a large density core, and the stellar distribution is more extended in SIDM than CDM and there is tight correlation between the DM core size and the stellar one (Vogelsberger et al., 2014b). While for baryon-dominated systems, the thermalization can significantly increase the central SIDM density and the inner halo shape follows the baryon distribution due to the influence of the baryonic potential (Kaplinghat et al., 2014), a radical deviation from SIDM-only simulations. Moreover, Kaplinghat et al. (2014) argued that once the inner halo reaches equilibrium, the inner SIDM profile can be modeled as an isothermal distribution that is sensitive to the final baryonic potential, but not the formation history. This has motivated a number of isolated simulations to test the response of the SIDM halo to the baryonic potential (Elbert et al., 2016; Creasey et al., 2017), where they assumed a CDM halo and a stellar potential as the initial condition. Recently, Robertson et al. (2017) performed cosmological hydrodynamical SIDM simulations of galaxy clusters and explicitly confirmed this expectation. In addition, although in the presence of strong baryonic feedback both CDM and SIDM could lead to similar density profiles (Fry et al., 2015), the internal structure of the SIDM halo is more robust to the inclusion of baryonic feedback, compared to its CDM counterpart, due to the rapid energy redistribution caused by the DM collisions (Robles et al., 2017).
In this paper, we utilize isolated simulations to study the response of the SIDM halo in the presence of the baryonic potential for Milky Way (MW)-sized galaxies, where the baryonic contribution to the potential is important. Our goal is to understand the interplay between the DM self-scattering strength and the baryonic concentration in shaping the SIDM distribution, and the significance of the potential in altering evolution history of the SIDM halo. In the first two sets of simulations, we vary both the baryonic concentration and in the range of , and study the variation of the SIDM predictions in the density profile and shape as a function of the cross section. In the presence of baryons, the central density of an SIDM halo no longer decreases monotonically with increasing , as expected in the SIDM-only case for we take. Accordingly, the SIDM halo shape varies with even for the same baryonic potential. Our results indicate that inferring from stellar kinematics of luminous galaxies, where the baryons dominate the potential, could be challenging.
In our third set of simulations, we construct a realistic MW mass model, including an SIDM halo, a stellar bulge and disk. We fix and carefully adjust the model parameters to reproduce the mass model inferred from the stellar kinematics. We then make a detailed comparison between the halo shape predicted in our model and those inferred from observations.
The structure of this paper is as follows. In Sec. 2, we discuss the numerical details of our simulations and the methodology used to quantify the halo shapes. In Sec. 3, we use our code to explore the evolution of a MW-sized halo with a stellar disk and measure the effect of the radial length scale of the disk. In Sec. 4, we compare our predictions for an SIDM MW halo against those from CDM simulations and those inferred from observations of stellar streams. We conclude and summarize our results in Sec. 5.
2 Simulations and halo shape algorithms
2.1 Numerical Simulations
We carry out N-body simulations using the code Arepo (Springel, 2010). Gravity modules in Arepo are a modified version of GADGET-2 and Gadget-3 (Springel, 2005). We use the algorithm developed in Vogelsberger et al. (2012); Vogelsberger et al. (2016) to model DM self-interactions. This is a Monte Carlo-based method, where at each time step a particle may pairwise scatter with any of its nearest neighbors. We assume a velocity-independent constant cross section in our simulations. This is a good approximation, since the observationally self-scattering cross section varies mildly across galactic scales (Kaplinghat et al., 2016) and we mainly focus on isolated simulations for a given halo mass. We evolve our simulations for , slightly shorter than Hubble time scale () in order to account for the assembly of the primordial galactic halo.
Following Creasey et al. (2017), we model the baryonic component in our simulations as a static potential. This approach ignores the back-reaction of the halo evolution on the baryons, an effect expected to be sub-dominant, since we are interested in the systems that the final baryonic distribution is known. We consider two models for the baryonic potential. One is the Miyamoto-Nagai (MN) disk (Miyamoto & Nagai, 1975),
where is the disk mass, the disk scale length and the disk scale height. The implementation of the MN disk in Arepo is as described in Creasey et al. (2017). We also consider a Hernquist bulge potential (Hernquist, 1990)
where is the bulge mass and is the scale length.
We run three sets of simulations, varying the baryonic component and the strength of the cross section. In the first two sets, we only include the MN disk () with two disk scale lengths, (compact disk) and (extended disk). In both cases, we fix . The DM self-scattering cross section is chosen to be , and , i.e., simulations in total. For the initial halo component, we assume a spherical NFW profile Navarro et al. (1997), and take the halo parameters as and . The mass ratio of the disk to the halo is motivated by the baryonic Tully-Fisher relation (Lelli et al., 2016). We use the publicly available code SpherIC, introduced in Garrison-Kimmel et al. (2013), to generate the initial conditions. It truncates the outer halo profile exponentially at to avoid mass divergence. We take , close to the virial radius of our initial CDM halo. We fix the gravitational softening length to be and the mass resolution . We include million mass particles in our simulations, necessary for resolving the innermost regions, resulting in a halo mass of .
For the third set, we include both an MN disk and a Hernquist bulge to model the baryon distribution in the MW. The disk parameters are , and for the disk. The bulge ones are and . The initial halo parameters are and . We have chosen these parameters to reproduce the MW mass model presented in McMillan (2011) (hereafter McM11), see Sec. 4 for details. The baryon-model parameters used in our simulations are summarized in Table 1. We choose kpc, and . We simulate million mass particles and the total halo mass is .
Additionally, we have run cosmological zoom-in SIDM simulations for MW-mass Aquarius halos (Springel et al., 2008) with the initial conditions taken from (Vogelsberger et al., 2012; Zavala et al., 2013). We will present these simulation results in Sec. 4 for comparison.
|Component||Mass ()||Length scale (kpc)|
2.2 Halo shape algorithm
We use the method introduced in Dubinski & Carlberg (1991) (see also Allgood et al., 2006) to calculate the ellipticity of the simulated halos. It constructs the axial ratio for the best-fitting ellipsoid as a function of the major axis length. This method is an iterative one where at each iteration the reduced inertia tensor is determined for the set of particles within the previous ellipsoid, and then a new ellipsoid is determined from this tensor. Specifically, if we denote the major axis length , then at each iteration the reduced inertia tensor is given by
where denotes the coordinate of particle , and is the elliptical radius found from the previous inertia tensor. We have , where , and are the coordinates along the major-, intermediate- and minor-axes of the ellipsoid, and , are the axial ratios of the intermediate- and minor-to-major axes, respectively. After diagonalizing the inertia tensor with eigenvalues (ascending) , we have and . In the initial iteration, the ellipsoid is set to a sphere, i.e. . This process is continued until some convergence criteria, which we take it to be on the difference between successive iterations, is satisfied. We note that if the number of DM particles in an ellipsoid is too small, typically less than , the result from this method is not accurate (see Appendix A).
Fig. 1 shows a comparison between isodensity contours and ellipsoids for an example, where kpc and . We see the overall agreement between the two methods is excellent.
3 SIDM halo properties with a stellar disk
3.1 Density profiles
Fig. 2 shows the DM density (top) and velocity dispersion (bottom) profiles for (left) and (right). The solid curves are from our simulations for different values of in the presence of the stellar disk. For comparison, we also plot the SIDM density profiles (dashed) without the disk potential, calculated with the analytical method in Kaplinghat et al. (2014); Kaplinghat et al. (2016).
In both cases, the presence of a baryonic potential increases the SIDM density profile and reduces the core size, and the effect is more significant if the baryonic concentration is higher. For , a larger cross section leads to a higher DM density (), opposite to the case without baryons. For the extended disk with , the SIDM density profiles are almost identical, even though the value changes by a factor of . A deep baryonic potential also increases the DM velocity dispersion in the inner halo, as shown in the bottom panels. In the case of compact disk, it is evident that the SIDM halos are close to the threshold of mild core collapse, as the velocity dispersion profiles start to develop a negative gradient. The significance is continuously enhanced when changes from to . We have checked simulation results using the analytical method, where we assume a thin disk model and use the numerical templates developed in Kamada et al. (2017). They agree within in the central DM density and velocity dispersion.
The SIDM halo has a distinct evolution history. It first undergoes a core expansion phase, during which the DM collisions transport heat towards the inner region and a central density core forms. Since a self-gravitating system has a negative heat capacity, the core will eventually contract and collapse to a singular state (Balberg et al., 2002). In cosmological SIDM-only simulations, mild core collapse is observed within when (Vogelsberger et al., 2012; Elbert et al., 2015). We have also checked that for isolated SIDM-only ones with an NFW profile as the initial condition, the core contraction does not occur within for , consistent with the results in Koda & Shapiro (2011). However, the SIDM thermalization with a deep baryonic potential can speed up this process, as shown in our simulations (see also Elbert et al., 2016).
We see that the presence of the stellar potential breaks the monotonic relation between the value of and the central SIDM density. The effect depends on the baryonic concentration and the size of the self-scattering cross section. Our results indicate that it could be challenging to extract the information from stellar kinematics of galaxies dominated by baryons.
3.2 Halo shapes
In Fig. 3, we show the SIDM halo surface densities for (left) and (right). The density contrast for the compact case is higher for different cross sections, compared to the extended one, as expected from the density profiles shown in Fig. 2. It is also evident that the simulated halos are not spherically symmetric, although their initial conditions are exactly spherical.
Fig. 4 shows the ratio of minor-to-major axes vs. elliptical radius for and with different cross sections (solid). In all cases, the value deviates from and decreases towards the center ( remains close to ). However, the SIDM halos are more responsive to the presence of the baryonic disk than their collisionless counterpart, and their shapes are more aligned with the axisymmetric disk potential (dashed). Interestingly, for the compact case, increases when the cross section increases from to and the inner halo becomes rounder mildly. We can see a similar trend in the case of , although the errors in measuring , calculated using bootstrap method, for are large due to the lack of enough DM particles in the central region of the halos.
The behavior in Fig. 4 can be understood as follows. Since the DM self-interactions thermalize the inner halo, the DM density can modeled by the isothermal distribution (Kaplinghat et al., 2014), , where and are the DM and disk potentials, respectively. induces the deviation from spherical symmetry of the initial NFW halo, as indicated in Fig. 4 (dashed), and the significance depends on its magnitude relative to and . In the compact-disk case, the central DM density increases when increases from to , as well as the DM dispersion (very mildly), as shown in Fig. 2 (left). Accordingly, the baryonic potential becomes less dominant and the inner halo becomes more spherical. Note in the compact-disk case the simulated profile for agrees well with the isothermal profile due to the baryonic potential, because of the strong dominance of the disk in the inner regions. In addition, for (compact), both the inner DM density and velocity dispersion are higher, compared to the CDM case, but the SIDM halo is more aspherical and aligned with the disk than the CDM one. In the extended-disk case, the halo profiles also follow the disk one, but not as close as the compact case, since the disk does not dominate the potential at all radii, as shown in Fig. 2 (right).
3.3 Evolution history
In this section, we take a closer look at the evolution of the SIDM halo and explicitly show that the presence of the baryonic potential does speed up core contraction and shorten the expansion phase.
Fig. 5 shows the density and profiles at different epochs for three examples: the lowest (top) and highest (middle) cross sections in the simulations with the extended disk, and the highest cross section for the compact disk case (bottom). For and , the simulated halo is on the core expansion phase over the span of the simulation. In this case both central DM density and the ratio decrease continuously. When we increase the cross section to (middle), the duration of the core expansion phase becomes much shorter. After about , the halo enters the core contraction phase and the central DM density increases, as well as the ratio of minor-to-major axes. A more compact stellar disk can change the halo evolution even more dramatically, as shown in the bottom panel. The simulated halo almost never gets into the expansion phase and the central density and in the regions increases over time monotonically. In this case, the inner SIDM halo contains even more DM mass than its CDM counterpart.
We conclude that the evolution history of the SIDM halo is sensitive to the presence of the baryonic potential. The final halo properties, such as the density profile and the ellipticity, depend on the baryonic concentration and the strength of DM self-interactions.
4 Implications for the shape of the Milky Way Halo
The effect of baryons on the SIDM halo can potentially be tested with observations of the MW. Here, we construct a model for the MW potential consisting of an SIDM halo, a baryonic bulge and disk. The DM halo is chosen initially as a spherical NFW profile with and . We model the disk following an MN potential as in Eq. 1, with disk length scale and mass specified in Table 1 and bulge following a spherical Hernquist profile. We take the cross section as .
Top panel of Fig. 6 shows the DM density profiles for our SIDM halo (blue), the initial NFW model (gray), and the halo model in McM11 (black). Our MW model reproduces well the estimates for the local DM density near the solar neighborhood from Bovy & Tremaine (2012) (red). Note the initial halo concentration is about lower than the average for the MW mass object according to Dutton & Macciò (2014) and also lower than that in McM11. This is a necessary choice to be consistent with observations, since SIDM thermalization significantly increases the DM density in the inner regions due to the presence of the baryonic potential. Although the inner density profile of the SIDM halo deviates from the McM11 one, our MW mass model produces a circular velocity profile, consistent with the one in McM11 within uncertainties, as shown in Fig. 6 (bottom).
Fig. 7 shows the ratio of minor-to-major axes as a function of the elliptical radius, predicted in our SIDM MW halo model with (blue solid). Since we assume a spherical NFW initial halo, the deviation from is caused by the disk potential. We see the disk induces a mild asphericity in the inner regions of the halo, , in good agreement with the result based on the analytical model in Kaplinghat et al. (2014), where the spherical boundary condition is imposed at . Our simulations also show the effect is quite extended, with converging back to unity only at the distance . We also present cosmological SIDM-only simulations of Aquarius halos with . Our results are in agreement with the previous ones (Peter et al., 2013; Brinckmann et al., 2018). For the roundest halos (Aq-A and Aq-C), for , lending support to our assumption of an initially spherical NFW halo. Since the DM self-scattering rate increases by a factor of from to inner few , as indicated by the DM density profile shown in Fig. 6 (top), we expect the spherical assumption is well-justified in the inner halo.
The MW halo shape has been inferred from observations of stellar streams such as GD-1 and Pal 5 (e.g., Koposov et al., 2010; Bowden et al., 2015; Pearson et al., 2015; Küpper et al., 2015; Bovy et al., 2016). In Fig. 7, we show the results presented in Bovy et al. (2016) for GD-1, Pal 5 and the combined one, where an axisymmetric NFW density profile with was used to model the DM distribution in the MW. We see that the phase-space tracks of the streams are consistent with a spherical DM halo in the MW at intermediate radii, , which indicates any asphericity either intrinsic to the DM distribution or induced by the disk should be at most weak on that scale in order to accommodate the measurements. Our model predicts at , consistent with the combined constraint on within . Note that our scale height in the MN disk is , higher than the best-fit value in Bovy et al. (2016). We have estimated that taking their value would reduce by at most in the inner regions and the difference becomes negligible around . In addition, in our MW model, the total halo mass within is , consistent with measured in Bovy et al. (2016), although our initial NFW halo has lower concentration, compared to theirs (). It would be interesting to analyze the stream data with the SIDM halo model.
For comparison, we also plot cosmological CDM-only (gray) and hydrodynamical (red) simulations from the NIHAO project (Butsky et al., 2016). As is well-known, CDM halos from cosmological simulations are strongly triaxial (Frenk et al., 1988; Jing & Suto, 2002; Hayashi et al., 2007; Kuhlen et al., 2007; Vera-Ciro et al., 2011; Diemand & Moore, 2011). Taken at face value, CDM predictions (gray) could look at odds with the observations. However, baryons can make the CDM halo shapes more spherical (red) (see also, Dubinski, 1994; Debattista et al., 2008; Kazantzidis et al., 2010; Abadi et al., 2010; Tissera et al., 2010), an effect partially attributed to the change of orbits from boxy to tube or rounder loop as a result of the central concentration of baryons (Debattista et al., 2008). In the NIHAO simulations, the mean value of is for the CDM halos after including baryons, and it can reaches at the level of the scatter, consistent with the observations reasonably well.
Although the sphericity created by the baryons helps CDM to accommodate more easily the observational constraints, the trend with radius could be different in the two models. It seems that CDM halos plus baryons become more spherical at all radii, whereas the effects explored here in SIDM plus baryons would anticipate a flattening of the shapes towards the inner regions that follows that of the disk. Such premise, of course, ignores any effect of feedback or cosmological assembly, which may cause deviations of the system from equilibrium. Therefore, the exciting premise of using halo shape profiles to differentiate DM candidates awaits confirmation from cosmological hydrodynamical SIDM simulations. We hope such experiments will become available in the near future.
We use isolated N-body simulations of DM halos with static disk potentials to explore the gravitational effect of baryons on SIDM halos. We model the disk as a Miyamoto-Nagai potential embedded within an initially NFW halo with mass . We consider different self-scattering cross sections, , and besides the special case, . In addition, we vary the radial length scale of the disk, , to study in detail how the DM halo responds to the baryons as a function of how relevant their contribution is to the total potential at a given radius.
In the absence of baryons, SIDM halos develop a central flat core with its density and size that depend mostly on the self-scattering cross section. We confirm that the inclusion of a disk potential can change this behavior due to SIDM thermalization with the potential, resulting in a higher core density and a smaller core than expected without the disk, a crucial effect in solving the diversity problem in SIDM (Kamada et al., 2017; Creasey et al., 2017),
We highlight two phases of evolution during our numerical experiments: a first stage of core expansion, during which the density core gets established due to the turn-on of the self-interactions, and a second stage of core contraction due to the gravitational effects of the baryons. The timescale for these two phases of evolution is a function of both, the cross section and the relative importance of the baryons inside the core. Higher cross sections and more compact baryonic disks (encoded in a smaller length scale) speed up the transition between the two phases and make the timescale of core expansion shorter.
We have also studied the role of the disk potential in shaping the SIDM halo. To explore this subtle effect, we assumed an exact spherical initial NFW profile such that any departure from sphericity is due to the influence of the baryonic potential. Compared to the case of , the SIDM halos are more responsive to the potential due to the thermalization, and their final flattening is more aligned with the orientation of the disk, consistent with the expectation from the analytical method. Our simulations clearly demonstrate that the induced asphericity is mainly sensitive to the contribution of the disk to the total potential, relative to the DM one. We further confirmed this by checking the evolution history of the SIDM halos. The flattening effect is maximized during the epoch when the core has the lowest density, which coincides with the time when the disk contribution to the total potential is also maximized.
We have constructed a mass model for the MW and explored the shape prediction with observations. The model consists of a stellar disk and a bulge, embedded within a spherical SIDM halo. It reproduces observed stellar kinematics of the Galaxy within the uncertainties and the local DM density reported in the solar neighborhood. We find that the baryons are able to induce a mild flattening () in the inner regions but the effect weakens at larger radii. At kpc where observational constraints seem to suggest an almost spherical halo, the effects of the disk are not strong, in agreement with the observations. We propose that the quasi sphericity of the halo at large distance is easier to accommodate in SIDM models than within the strongly triaxial structures predicted by CDM, although considering the effects of baryons might help to reconcile CDM models with observed spherical halos. Furthermore, we argue that a study of halo shapes as a function of radius might be able to help distinguish the nature of DM, although a more stringent comparison to cosmological simulations are needed to confirm this last point. On the observational side, future surveys capable of inferring the shape of the Galactic halo within the inner regions are promising avenues to make progress on establishing the non-canonical nature of DM.
OS acknowledges support by NASA MUREP Institutional Research Opportunity (MIRO) grant number NNX15AP99A and HST grant HST-AR-14582. HBY acknowledges support from U. S. Department of Energy under Grant No. de-sc0008541 and the Hellman Fellows Fund. LVS is grateful for support from the Hellman Fellows Foundation and HST grant HST-AR-14582. MV acknowledges support through an MIT RSC award, the support of the Alfred P. Sloan Foundation, and support by NASA ATP grant NNX17AG29G. JZ acknowledges support by a Grant of Excellence from the Iceland Research Fund (grant number 173929-051).
Appendix A Convergence test for the halo shape algorithm
To evaluate accuracy of the shape measurement, we follow Vera-Ciro et al. (2011) to determine the convergence radius where shape measurements are robust in our simulations. We consider the convergence radius (Power et al., 2003; Navarro et al., 2010),
where is the two-body relaxation time scale due to gravitational encounters, is the circular orbit timescale at , is number of DM particles and is the average density inside the convergence radius. We take as in Vera-Ciro et al. (2011). In addition, we require of the particles to be inside the virial sphere and at least DM particles inside of convergence radius.
The first requirement is satisfied if we choose a large cutoff radius in the SpherIC code, at which the density profile transits from an NFW one to an exponentially decaying one to avoid the divergence of the mass. Then we run three simulations with different values of the DM particle number and the gravitational softening length as a convergence test with details summarized in Table 2. The convergence radius decreases with increasing the number of total DM particles. Thus, to probe the shape of the inner halo, down to few , we need at least million particles in simulations. Fig. 8 shows the (top) and (bottom) profiles for different resolutions, we take a static MN potential as stellar disk similar to Sec. 3, with a MW-sized halo with . We see that the convergence improves when increases. We take high-resolution run in the results presented in Sec. 3 and 4.
- Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010, MNRAS, 407, 435
- Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
- Balberg S., Shapiro S. L., Inagaki S., 2002, ApJ, 568, 475
- Bovy J., Tremaine S., 2012, ApJ, 756, 89
- Bovy J., Bahmanyar A., Fritz T. K., Kallivayalil N., 2016, ApJ, 833, 31
- Bowden A., Belokurov V., Evans N. W., 2015, MNRAS, 449, 1391
- Brinckmann T., Zavala J., Rapetti D., Hansen S. H., Vogelsberger M., 2018, MNRAS, 474, 746
- Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
- Butsky I., et al., 2016, MNRAS, 462, 663
- Colín P., Avila-Reese V., Valenzuela O., Firmani C., 2002, ApJ, 581, 777
- Creasey P., Sameie O., Sales L. V., Yu H.-B., Vogelsberger M., Zavala J., 2017, MNRAS, 468, 2283
- Davé R., Spergel D. N., Steinhardt P. J., Wandelt B. D., 2001, ApJ, 547, 574
- Debattista V. P., Moore B., Quinn T., Kazantzidis S., Maas R., Mayer L., Read J., Stadel J., 2008, ApJ, 681, 1076
- Di Cintio A., Tremmel M., Governato F., Pontzen A., Zavala J., Bastidas Fry A., Brooks A., Vogelsberger M., 2017, Mon. Not. Roy. Astron. Soc., 469, 2845
- Diemand J., Moore B., 2011, Advanced Science Letters, 4, 297
- Dubinski J., 1994, ApJ, 431, 617
- Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
- Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
- Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., Oñorbe J., Peter A. H. G., 2015, MNRAS, 453, 29
- Elbert O. D., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Graus A. S., Rocha M., 2016, preprint, (arXiv:1609.08626)
- Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507
- Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507
- Fry A. B., et al., 2015, Mon. Not. Roy. Astron. Soc., 452, 1468
- Garrison-Kimmel S., Rocha M., Boylan-Kolchin M., Bullock J. S., Lally J., 2013, MNRAS, 433, 3539
- Hayashi E., Navarro J. F., Springel V., 2007, MNRAS, 377, 50
- Hernquist L., 1990, ApJ, 356, 359
- Jing Y. P., Suto Y., 2002, ApJ, 574, 538
- Kamada A., Kaplinghat M., Pace A. B., Yu H.-B., 2017, Physical Review Letters, 119, 111102
- Kaplinghat M., Keeley R. E., Linden T., Yu H.-B., 2014, Physical Review Letters, 113, 021302
- Kaplinghat M., Tulin S., Yu H.-B., 2016, Physical Review Letters, 116, 041302
- Kazantzidis S., Abadi M. G., Navarro J. F., 2010, ApJ, 720, L62
- Koda J., Shapiro P. R., 2011, MNRAS, 415, 1125
- Koposov S. E., Rix H.-W., Hogg D. W., 2010, ApJ, 712, 260
- Kuhlen M., Diemand J., Madau P., 2007, ApJ, 671, 1135
- Küpper A. H. W., Balbinot E., Bonaca A., Johnston K. V., Hogg D. W., Kroupa P., Santiago B. X., 2015, ApJ, 803, 80
- Kuzio de Naray R., Martinez G. D., Bullock J. S., Kaplinghat M., 2010, Astrophys. J., 710, L161
- Lelli F., McGaugh S. S., Schombert J. M., 2016, ApJ, 816, L14
- McMillan P. J., 2011, MNRAS, 414, 2446
- Miyamoto M., Nagai R., 1975, PASJ, 27, 533
- Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
- Navarro J. F., et al., 2010, MNRAS, 402, 21
- Oman K. A., et al., 2015, MNRAS, 452, 3650
- Pearson S., Küpper A. H. W., Johnston K. V., Price-Whelan A. M., 2015, ApJ, 799, 28
- Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
- Planck Collaboration et al., 2014, A&A, 571, A16
- Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
- Robertson A., et al., 2017, preprint, (arXiv:1711.09096)
- Robles V. H., et al., 2017, MNRAS, 472, 2945
- Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81
- Spergel D. N., Steinhardt P. J., 2000, Physical Review Letters, 84, 3760
- Springel V., 2005, MNRAS, 364, 1105
- Springel V., 2010, MNRAS, 401, 791
- Springel V., Frenk C. S., White S. D. M., 2006, Nature, 440, 1137
- Springel V., et al., 2008, Mon. Not. Roy. Astron. Soc., 391, 1685
- Tissera P. B., White S. D. M., Pedrosa S., Scannapieco C., 2010, MNRAS, 406, 922
- Trujillo-Gomez S., Klypin A., Primack J., Romanowsky A. J., 2011, Astrophys. J., 742, 16
- Tulin S., Yu H.-B., 2017, preprint, (arXiv:1705.02358)
- Vera-Ciro C. A., Sales L. V., Helmi A., Frenk C. S., Navarro J. F., Springel V., Vogelsberger M., White S. D. M., 2011, MNRAS, 416, 1377
- Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740
- Vogelsberger M., et al., 2014a, MNRAS, 444, 1518
- Vogelsberger M., Zavala J., Simpson C., Jenkins A., 2014b, MNRAS, 444, 3684
- Vogelsberger M., et al., 2014c, Nature, 509, 177
- Vogelsberger M., Zavala J., Cyr-Racine F.-Y., Pfrommer C., Bringmann T., Sigurdson K., 2016, MNRAS, 460, 1399
- Yoshida N., Springel V., White S. D. M., Tormen G., 2000, ApJ, 544, L87
- Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20