NIHAO – IV: core creation and destruction in dark matter density profiles across cosmic time

Edouard Tollet, Andrea Macciò, Aaron Dutton, Greg Stinson, Liang Wang, Camilla Penzo, Thales Gutcke, Tobias Buck, Xi Kang, Chris Brook, et al.

► To cite this version:

Edouard Tollet, Andrea Macciò, Aaron Dutton, Greg Stinson, Liang Wang, et al. NIHAO – IV: core creation and destruction in dark matter density profiles across cosmic time. Monthly Notices of the Royal Astronomical Society, 2016, 456 (4), pp.3542-3552. 10.1093/mnras/stv2856. hal-02316499

HAL Id: hal-02316499

https://hal.science/hal-02316499v1

Submitted on 10 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# NIHAO – IV: core creation and destruction in dark matter density profiles across cosmic time

Edouard Tollet,<fnref n="1" />,<fnref n="2" />* Andrea V. Macciò,<fnref n="1" />,<fnref n="3" />* Aaron A. Dutton,<fnref n="1" />,<fnref n="3" /> Greg S. Stinson,<fnref n="1" />
Liang Wang,<fnref n="1" />,<fnref n="4" /> Camilla Penzo,<fnref n="1" /> Thales A. Gutcke,<fnref n="1" /> Tobias Buck,<fnref n="1" /> Xi Kang,<fnref n="4" />
Chris Brook,<fnref n="5" />† Arianna Di Cintio,<fnref n="6" />‡ Ben W. Keller<fnref n="7" /> and James Wadsley<fnref n="7" />

Accepted 2015 December 3. Received 2015 November 4; in original form 2015 July 13

## ABSTRACT

We use the NIHAO (Numerical Investigation of Hundred Astrophysical Objects) cosmological simulations to investigate the effects of baryonic physics on the time evolution of dark matter central density profiles. The sample is made of ≈70 independent high-resolution hydrodynamical simulations of galaxy formation and covers a wide mass range: $10^{10} \lesssim M_{\text{halo}}/M_{\odot} \lesssim 10^{12}$, i.e. from dwarfs to L*. We confirm previous results on the dependence of the inner dark matter density slope, $\alpha$, on the ratio between stellar-to-halo mass, $M_{\text{star}}/M_{\text{halo}}$. We show that this relation holds approximately at all redshifts (with an intrinsic scatter of ~0.18 in $\alpha$ measured between 1 and 2 per cent of the virial radius). This implies that in practically all haloes the shape of their inner density profile changes quite substantially over cosmic time, as they grow in stellar and total mass. Thus, depending on their final $M_{\text{star}}/M_{\text{halo}}$ ratio, haloes can either form and keep a substantial density core ($R_{\text{core}} \sim 1$ kpc), or form and then destroy the core and recontract the halo, going back to a cuspy profile, which is even steeper than cold-dark-matter predictions for massive galaxies ($10^{12} M_{\odot}$). We show that results from the NIHAO suite are in good agreement with recent observational measurements of $\alpha$ in dwarf galaxies. Overall our results suggest that the notion of a universal density profile for dark matter haloes is no longer valid in the presence of galaxy formation.

**Key words:** hydrodynamics – galaxies: evolution – dark matter.

## 1 INTRODUCTION

The cold-dark-matter (CDM) theory is very successful at describing the Universe's topology and evolution on large scales (e.g. Springel et al. 2005). This theory makes clear predictions of the distribution of DM on small scales, where all collapsed structures (haloes) are supposed to share an approximately self-similar DM density profile, as first pointed out by Navarro, Frenk & White (1997, NFW) and then confirmed via higher resolution numerical simulations by a large variety of works (e.g. Klypin et al. 2001; Power et al. 2003;

Diemand, Moore & Stadel 2004; Macciò et al. 2007; Neto et al. 2007; Macciò, Dutton & van den Bosch 2008; Navarro et al. 2010; Prada et al. 2012; Dutton & Macciò 2014).

This prediction of a universal cuspy profile for the DM seems to be in tension with observations of rotation curves of low-mass galaxies (e.g. Moore 1994; Salucci & Burkert 2000; Simon et al. 2005; de Blok et al. 2008; Kuzio de Naray, McGaugh & de Blok 2008; Kuzio de Naray, McGaugh & Mihos 2009; Oh et al. 2011a), which seem to prefer density profiles with a constant density core in the centre.

This cusp-core controversy suggests that either a modification of the whole CDM paradigm is required (e.g. self-interacting DM; Vogelsberger et al. 2014), or the inadequacy of pure N-body simulations to capture the DM dynamics on small scales, due to the absence of dissipative phenomena connected to baryonic physics.

<footnote n="1">¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="2">²GEPI, Observatoire de Paris, CNRS, Univ. Paris Diderot, 61 Avenue de l'Observatoire, F-75014 Paris, France</footnote>

<footnote n="3">³New York University Abu Dhabi, PO Box 129188, Abu Dhabi, UAE</footnote>

<footnote n="4">⁴Purple Mountain Observatory, the Partner Group of MPI für Astronomie, 2 West Beijing Road, Nanjing 210008, China</footnote>

<footnote n="5">⁵Departamento de Física Teórica, Universidad Autonoma de Madrid, E-28049 Cantoblanco, Madrid, Spain</footnote>

<footnote n="6">⁶Dark Cosmology Centre, NBI, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark</footnote>

<footnote n="7">⁷Department of Physics and Astronomy, McMaster University, Hamilton, ON L8S 4M1, Canada</footnote>

<footnote>* E-mail: edouard.tollet@obspm.fr (ET); maccio@nyu.edu (AVM)</footnote>

<footnote>† Ramon y Cajal Fellow.</footnote>

<footnote>‡ Dark Fellow.</footnote>
Recent and more accurate cosmological hydrodynamical simulations have indeed shown that baryons are able to alter the DM profile and to create substantial cores (up to ~1 kpc) in the DM distribution (e.g. Governato et al. 2010; Macciò et al. 2012; Martizzi, Teyssier & Moore 2013; Munshi et al. 2013; Di Cintio et al. 2014a,b; Oñorbe et al. 2015; Trujillo-Gomez et al. 2015).

As nicely described in Pontzen & Governato (2012) the formation of DM cores is linked to the rapid change of the total gravitational potential. In the central kiloparsecs the potential changes on subdynamical time-scales, as repeated, energetic feedback generate large underdense bubbles of expanding gas. The fluctuations in the central potential irreversibly (e.g. non-adiabatically) transfer energy into collisionless particles, thus generating a DM core (see also Ogiya & Mori 2014, for a detailed analysis of the effect of resonances between DM particles and the density wave excited by the oscillating potential). These strong gas outflows are triggered by stellar (and active galactic nuclei) feedback (Navarro, Eke & Frenk 1996; Read & Gilmore 2005; Mashchenko, Couchman & Wadsley 2006; Martizzi et al. 2013; Pontzen & Governato 2014).

Recently Di Cintio et al. (2014b, hereafter DC14) have shown that the modification of the initial DM profile (either leading to an expansion or a contraction) is clearly linked to the integrated star formation efficiency of the galaxy, which can be parametrized through the redshift zero stellar mass-halo mass ratio ($M_{\text{star}}/M_{\text{halo}}$), and that this holds for a large variety of stellar feedback implementations.

They clearly showed that the halo response to star formation is non-monotonic and that haloes with very low and very high star formation efficiency (or with very low and high $M_{\text{star}}/M_{\text{halo}}$) tend to preserve the initial cuspy profile, while haloes with $M_{\text{star}}/M_{\text{halo}} \approx 0.3$ per cent have a flat central density profile.

In this paper we want to expand the original work of DC14, by using a newer (and larger) sample of simulated galaxies. This new set of simulations has been performed using the new cosmological parameters from the Planck satellite (the Planck Collaboration XVI 2014), and an updated version of the GASOLINE code (Keller et al. 2014) which fixes some of the known problems of the SPH technique (Agertz et al. 2007). This newer and larger sample allows us to look at the evolution of the density profiles through cosmic time and to witness the creation (and destruction) of density cores as a function of stellar mass and halo mass.

This paper is organized as follows: in Section 2 we give an overview of the cosmological hydrodynamical simulations used in this work. Results are presented in Section 3, including core creation in Section 3.2 and core destruction in Section 3.3. Our conclusions are given in Section 4.

## 2 SIMULATIONS

In this study we use simulations from the Numerical Investigation of Hundred Astrophysical Objects (NIHAO) suite (Wang et al. 2015), based on an updated version of the MaGICC project (Stinson et al. 2012). All the simulations have adopted the latest compilation of cosmological parameters from the Planck satellite (the Planck Collaboration XVI 2014), namely $\Omega_m = 0.3175$, $H_0 = 67.1$, $\sigma_8 = 0.8344$, $n = 0.9624$ and $\Omega_b = 0.0490$. The haloes to be resimulated with baryons and higher resolution have been extracted from three different pure N-body simulations with a box size of 60, 20 and $15 h^{-1} \text{Mpc}$, respectively; more information on the collisionless simulations can be found in Dutton & Macciò (2014).

The aim of the NIHAO project is to study galaxy formation over a large mass range, from dwarf galaxies to massive spirals

such as the Milky Way. We have decided to keep the same relative mass resolution across the whole mass range, meaning that at high resolution we have (roughly) the same number of DM particles within the virial radius of our galaxies. This requirement sets the mass of the DM particle and the initial mass of the gas particles, since the latter is simply obtained by the DM mass multiplied by $(\Omega_m - \Omega_b)/\Omega_b$.

The zoomed initial conditions were created using a modified version of GRAFIC2 (Bertschinger 2001; Penzo et al. 2014). The starting redshift is $z_{\text{start}} = 99$, and each halo is initially simulated at high resolution with DM only using PKDGRAV (Stadel 2001). More details on the sample selection can be found in Wang et al. (2015). We refer to simulations with baryons as the hydro simulation, while we will use the term N-body for the DM-only simulation. Table 1 lists key properties for three galaxies across the mass range that are analysed in more detail; the complete list of NIHAO galaxies can be found in Wang et al. (2015).

The hydrodynamical simulations are evolved using an improved version of the SPH code GASOLINE (Wadsley, Stadel & Quinn 2004). The code includes a subgrid model for turbulent mixing of metals and energy (Wadsley et al. 2008), heating and cooling include photoelectric heating of dust grains, ultraviolet heating and ionization and cooling due to hydrogen, helium and metals (Shen, Wadsley & Stinson 2010).

For the NIHAO simulations we have used a revised treatment of hydrodynamics as described in Keller et al. (2014) that we refer to as ESF-GASOLINE2. Most important is the Ritchie & Thomas (2001) force expression that improves mixing and shortens the destruction time for cold blobs (see Agertz et al. 2007). ESF-GASOLINE2 also includes the time-step limiter suggested by Saitoh & Makino (2009), which is important in the presence of strong shocks and temperature jumps. We also increased the number of neighbour particles used in the calculation of the smoothed hydrodynamic properties from 32 to 50.

### 2.1 Star formation and feedback

The star formation and feedback modelling follows what was used in the MaGICC simulations (Stinson et al. 2013). Gas can form stars when it satisfies a temperature and a density threshold: $T < 15\,000\,\text{K}$ and $n_{\text{th}} > 10.3\,\text{cm}^{-3}$. Stars can feed energy back into the interstellar medium (ISM) via blast-wave supernova (SN) feedback (Stinson et al. 2006; see Agertz et al. 2013, for a discussion of the implementation of different energy and momentum stellar feedback in galaxy formation simulations) and via ionizing radiation from massive stars before they turn into SN. Metals are produced by Type II and Type Ia SN. These, along with stellar winds from asymptotic giant branch stars also return mass to the ISM. The metals affect the cooling function (Shen et al. 2010) and diffuse between gas particles (Wadsley, Veeravalli & Couchman 2008). The fraction of stellar mass that results in SN and winds is determined using the Chabrier (2003) stellar initial mass function.

There are two small changes from the MaGICC simulations. The change in number of neighbours and the new combination of softening length and particle mass means the threshold for star formation increased from $9.3$ to $10.3\,\text{cm}^{-3}$. The increased hydrodynamic mixing necessitated a small increase of pre-SN feedback efficiency, $\epsilon_{\text{ESF}}$, from $0.1$ to $0.13$. This energy is ejected as thermal energy into the surrounding gas, which does not have its cooling disabled. Most of this energy is instantaneously radiated away, and the effective coupling is of the order of $1$ per cent.
<table><caption>**Table 1.** Properties of our three selected galaxies at $z = 0$.</caption><thead><tr><th>Simulation ID</th><th>No. of particles<br>$n_{200}$</th><th>No. of DM particles<br>$n_{DM}$</th><th>No. of star particles<br>$n_{star}$</th><th>DM particle mass<br>$m_{DM}$</th><th>Virial mass<br>($M_{200}/M_{\odot}$)</th><th>Stellar mass<br>($M_{star}/M_{\odot}$)</th><th>$z=0$ SFR<br>($M_{\odot}$ yr$^{-1}$)</th><th>Conv. rad.<sup>a</sup><br>(kpc)</th></tr></thead><tbody><tr><td>g2.63e10</td><td>458 723</td><td>414 291</td><td>18 388</td><td>1.173 × 10<sup>4</sup></td><td>2.70 × 10<sup>10</sup></td><td>4.28 × 10<sup>7</sup></td><td>0.00</td><td>0.41</td></tr><tr><td>g2.19e11</td><td>920 447</td><td>557 247</td><td>113 958</td><td>2.169 × 10<sup>5</sup></td><td>1.31 × 10<sup>11</sup></td><td>9.27 × 10<sup>8</sup></td><td>8.2 × 10<sup>-2</sup></td><td>0.65</td></tr><tr><td>g8.06e11</td><td>1366 038</td><td>481 349</td><td>665 796</td><td>1.735 × 10<sup>6</sup></td><td>9.43 × 10<sup>11</sup></td><td>4.48 × 10<sup>10</sup></td><td>11.31</td><td>1.23</td></tr></tbody></table>

Note. $^a$Following Power et al. (2003); see text for details.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Stellar mass–halo mass relation for the NIHAO galaxies used in this work. All simulations have more than 400 000 particles in their virial radius (see Wang et al. 2015 for more details). The solid lines represent the most commonly used abundance matching results (see text).</figcaption></figure>

## 2.2 NIHAO galaxies’ stellar masses

Since our aim is to study the impact of galaxy formation on the DM distribution, it is very important to use realistic simulated galaxies, i.e. galaxies able to reproduce the observed scaling relations. Haloes in our simulations were identified using halo finder AHF<fnref n="1" /> (Gill, Knebe & Gibson 2004; Knollmann & Knebe 2009), adopting a spherical density contrast of 200 times the cosmic critical matter density. We dubbed the size of such a sphere $R_{vir}$; finally the stellar mass, $M_{star}$, is measured within a sphere of radius, $r_{gal} = 0.2R_{vir}$.

Fig. 1 shows the stellar mass–halo mass relation for the NIHAO galaxies used in this paper, compared with the most commonly used abundance matching relations from the literature: Brook et al. (2014), Kravtsov, Vikhlinin & Meshscheryakov (2014), Garrison-Kimmel et al. (2014), Behroozi, Wechsler & Conroy (2013) and Moster, Naab & White (2013). This plot is similar to the one showed in Wang et al. (2015) with the only difference that here only the most resolved halo in each simulation is shown. As detailed in Wang et al. (2015) NIHAO galaxies are able to reproduce the stellar mass–halo mass relation also at higher redshift, and have realistic star formation rates for their stellar masses. Overall the unprecedented combination of high resolution and large statistical sample of the NIHAO suite offers a unique tool to study the response of DM to galaxy formation.

## 2.3 Profile fitting

To construct and fit the DM density profiles we used the same methodology as introduced in DC14. Each DM density profile (in both hydro and N-body simulations) has been computed using regularly spaced shells on a logarithmic scale. The centre of the halo has been determined using the shrinking sphere method (Power et al. 2003) and we set the minimum radius for the first shell as twice the softening length of the DM particles. Each density profile is computed up to the virial radius, which is determined using the spherical overdensity criterion with a virial overdensity of $\Delta_{vir} = 200$ times the critical density of the universe. The error on the density in each bin is set according to the Poisson noise due to the finite number of particles each density bin.

Following DC14 the DM central density slope ($\alpha$) is subsequently fitted using a single power law, $\rho \propto r^{\alpha}$, over a limited radial range, namely $0.01 < r/R_{vir} < 0.02$, where $R_{vir}$ is the virial radius.

## 2.4 Profile convergence

The determination of the minimum radius above which the results of a simulation are not affected by the finite resolution is a non-trivial task (e.g. Gao et al. 2008). A quite commonly used criterion for convergence has been suggested by Power et al. (2003) for collisionless simulations, and it is based on the two body relaxation time-scale for particles in a gravitational potential. This empirical criterion radius ensures that the mean density inside the convergence radius is within 10 per cent of the value obtained in a simulation of much higher resolution (Schaller et al. 2015).

Due to choice of having similar particle resolution in all the NIHAO galaxies, this convergence radius is quite constant across the whole mass range and at $z = 0$ is $\approx 0.4–0.7$ per cent of $R_{vir}$. At $z = 2$, the inner region is marginally resolved if we use the original definition of Power et al. (2003), since we obtain a convergence radius of about 2 per cent of the virial radius. On the other hand as noticed by Schaller et al. (2015), this definition is quite restrictive and it was based on pure DM simulations. For hydrodynamical simulations it can be relaxed by requiring a 20 per cent convergence in the enclosed density (instead of 10 per cent). If we follow the prescription of Schaller et al. (2015) the convergence radius is of the order 1.0–1.5 per cent of $R_{vir}$; this means that while being ‘on the edge’ we can still use the same definition for $\alpha$ up to $z = 2$, as long as we do not try to be too quantitative at these redshifts.

# 3 RESULTS

Thanks to the large range of masses covered by our simulations we are able to see both expansion and contraction at $z = 0$ of the DM profile, depending on the halo mass and stellar content (see DC14 and references therein).

An example of such a dichotomy is presented in Fig. 2, where we show a comparison of the DM density profile at $z = 0$ in the hydro and in the N-body run for three galaxies with different total

<footnote n="1">¹ http://popia.ft.uam.es/AHF/Download.html</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Density profiles for three different simulations; on each plot we present the profile at $z = 0$ for the collisionless N-body (in black) and hydrodynamic (in red) simulation. The lines show power-law fits to the inner 1–2 per cent of the halo, and are shown to 0.1R$_{vir}$ for clarity. The upper panel shows an example where the DM halo contracts (halo mass of $9.43 \times 10^{11}$ M$_{\odot}$ and stellar mass of $4.48 \times 10^{10}$ M$_{\odot}$), the middle panel shows a galaxy with a DM halo that expands (host halo mass of $1.31 \times 10^{11}$ M$_{\odot}$ stellar mass of $9.27 \times 10^8$ M$_{\odot}$), while the bottom panel shows a galaxy with very little change (halo mass of $2.70 \times 10^{10}$ M$_{\odot}$ and stellar mass of $4.28 \times 10^7$ M$_{\odot}$).</figcaption></figure>

masses. The middle panel of Fig. 2 shows the DM density profile for a galaxy with a virial mass of $\approx 2.2 \times 10^{11}$ M$_{\odot}$. The profile in the hydro case (red points) shows a quite extended core. The situation is reversed for a more massive halo ($\approx 9.4 \times 10^{11}$ M$_{\odot}$, upper panel) where the DM profile is contracted in the hydro simulation with respect to the N-body run (black points). The lower panel shows a low-mass galaxy (halo mass $\approx 2.7 \times 10^{10}$ M$_{\odot}$), in this case both the profiles are cuspy, but it is clear that the hydro simulation has a shallower central slope than the N-body one.

This plot shows that simulations that include a realistic treatment of hydrodynamics do not share a universal DM density profile in contrast with DM-only simulations. We go into detail about how the universality is broken in Section 4.

As already pointed out by DC14 there is tight relation between the DM response and the star formation efficiency of a galaxy, defined as the ratio between stellar mass and halo mass. Fig. 3 shows the relation between $\alpha$ and the $M_{star}/M_{halo}$ ratio. The red points with error-bars are the NIHAO simulations results at $z = 0$.

The black open circles are the results from the N-body runs, and for these simulations we use the stellar mass of their hydrodynamical counterpart.

The NIHAO results show quite an interesting behaviour: the inner slope, $\alpha$, varies as a function of star formation efficiency. At low values of $M_{star}/M_{halo}$ the N-body and hydro results are quite similar and both predict a constant value of $\alpha \sim -1.5$. Then $\alpha$ increases steadily and reaches a maximum for $M_{star}/M_{halo} \sim 6 \times 10^{-3}$, in agreement with previous findings from DC14 (black dashed line). At the highest masses, after the maximum expansion, $\alpha$ decreases and for very large values of $M_{star}/M_{halo}$ of the order of few per cent, the hydro slopes become steeper than their N-body counterparts, suggesting halo contraction on those scales. This is in agreement with previous results from Di Cintio et al. (2014a), where it was also found that at the mass scales of the Milky Way, the DM profiles have a much higher concentration parameter in hydro simulations compared to collisionless results. It is also worth remembering that, as shown in Fig. 1, all our galaxies are on the observed stellar mass–halo mass relation (e.g. Behroozi et al. 2013; Kravtsov et al. 2014), and hence our halo contraction is not due to unphysical overcooling as it happened in older simulations (e.g. Gnedin et al. 2004) and it is consistent with recent results from other groups on the same mass scale (Schaller et al. 2015).

We tried to capture the behaviour of the slope of the density profile as a function of the star formation efficiency ($x = M_{star}/M_{halo}$) with a different fitting formula than DC14:

$$ \alpha(x) = n - \log_{10} \left[ n_1 \left( 1 + \frac{x}{x_1} \right)^{-\beta} + \left( \frac{x}{x_0} \right)^{\gamma} \right]. \quad (1) $$

This new fitting function has the advantage that it converges to a fixed value ($\approx -1.5$) for $x$ that goes to zero, in contrast to the power-law behaviour of the initial suggestion from DC14. Since, for continuity, we do expect to recover the slope of the N-body results for a negligible amount of stars ($x \lesssim 10^{-5}$), this new formula provides a better fit to the simulation data. The grey area is the 1$\sigma$ scatter around the mean value of $\alpha$ of $\approx 0.18$.

The difference between the two curves for $M_{star}/M_{halo} > 10^{-3}$ is mainly due to the different cosmological model adopted in the two studies (Planck versus WMAP3), while other differences, as for example the improved hydrodynamics, are just secondary effects (see Stinson et al. 2015 for a more detailed comparison between MaGICC and NIHAO galaxies). In the Planck cosmology the halo concentrations are higher by $\approx 45$ per cent (Dutton & Macciò 2014), which explains why we get lower (cuspier) $\alpha$-values with respect to DC14. Table 2 lists the values of our new relation for $\alpha$ versus $M_{star}/M_{halo}$, (updated to the Planck cosmology).

In Fig. 4 we show the same relation as in Fig. 3 but at redshift $z = 1$ (upper panel) and at redshift $z = 2$ (lower panel). The line in the plot is still the $z = 0$ fitting result. As the galaxy evolves, DM accretion and star formation generally move points from left to right. The points still cluster around the $z = 0$ relation, but scatter increases with redshift, so integrated star formation efficiency becomes a less useful metric when the galaxies are young.

The origin of the scatter is twofold. Part of it is due to resolution effects, at higher redshift the virial radius is smaller, this means that we are probing regions closer to our resolution limit which implies a lower number of particles in the inner region of the halo and hence a larger Poisson noise. Part of the scatter is also physical and related to the bursty nature of the star formation history at high redshift, that induces short time variations in the density profile (see next section).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Inner DM density slope, $\alpha$, as a function of $M_{\text{star}}/M_{\text{halo}}$ for galaxies at $z=0$. The filled red circles show results for hydro runs (with error-bars from the profile fit) while the black open circles show results for the N-body DM-only runs. The dotted line is the original function from DC14, the thick line is the refitted curve with our new simulations using the Planck cosmology with the shaded region showing the $1\sigma$ scatter of 0.18.</figcaption></figure>

<table><caption>**Table 2.** Best-fitting parameters for the value of $\alpha$ computed within 1 and 2 per cent of $R_{\text{vir}}$ as a function of $M_{\text{star}}/M_{\text{halo}}$, $M_{\text{halo}}$ and $M_{\text{star}}$.</caption><thead><tr><th></th><th>n</th><th>n<sub>1</sub></th><th>x<sub>0</sub></th><th>x<sub>1</sub></th><th>β</th><th>γ</th></tr></thead><tbody><tr><td>M<sub>star</sub>/M<sub>halo</sub></td><td>-0.158</td><td>26.49</td><td>8.77 × 10<sup>-3</sup></td><td>9.44 × 10<sup>-5</sup></td><td>0.85</td><td>1.66</td></tr><tr><td>M<sub>halo</sub></td><td>0.96</td><td>7161</td><td>9.96 × 10<sup>9</sup></td><td>1.85 × 10<sup>8</sup></td><td>1.01</td><td>1.13</td></tr><tr><td>M<sub>star</sub></td><td>0.53</td><td>61.44</td><td>2.48 × 10<sup>6</sup></td><td>6.83 × 10<sup>8</sup></td><td>29.5</td><td>0.42</td></tr></tbody></table>

To check whether star formation efficiency best correlates with $\alpha$ in our large sample of galaxies, Figs 5 and 6 show how $\alpha$ varies with the halo mass and the stellar mass of our galaxies, respectively. In both cases the behaviour of $\alpha$ is similar to the one seen in Fig. 3, with a maximum in ‘core’ creation around $10^{11}$ ($10^8$) M$_{\odot}$ in halo (stellar) mass. These plots show a similar scatter in $\alpha$ as in Fig. 3 for low values of the halo (stellar) mass: $M_{\text{halo}} < 5 \times 10^{11}$ M$_{\odot}$ ($M_{\text{star}} < 10^9$ M$_{\odot}$), but a larger scatter for $\alpha$ at high masses. This confirms the earlier results of DC14 where, using galaxies run with different stellar feedback prescription, they show that the integrated star formation efficiency is the best parameter to capture the effect of baryons on the DM distribution.

The role of the integrated star formation efficiency is also important to explain the similarities and (small, partial) differences between our results and some previous studies on the same topic from different groups.

Madau, Shen & Governato (2014) analysed the evolution of the density profile of a small group of dwarf galaxies, with only four of them with containing stars. The simulations were run with an earlier version of the GASOLINE code, with similar SN feedback. Two of those galaxies (dubbed Doc and Bashful) have masses around

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Inner DM density slope, $\alpha$, as a function of $M_{\text{star}}/M_{\text{halo}}$ for galaxies at $z=1$ (upper panel) and at $z=2$ (lower panel). The solid line and shaded region is $z=0$ relation from Fig. 3.</figcaption></figure>

$1-3 \times 10^{10}$ M$_{\odot}$ and do show extended cores, somehow in contrast with the results shown in Fig. 5. On the other hand the galaxies presented by Madau and collaborators have overmassive stellar bodies, that overpredict results from Local Group abundance matching
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Inner DM density slope, $\alpha$, as a function of the halo mass for galaxies at $z=0$. The solid line and shaded region are the fitted relation using equation (1) and its scatter. The parameters of the fitting are reported in Table 2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Inner DM density slope, $\alpha$, as a function of the stellar mass for galaxies at $z=0$. The solid line and shaded region are the fitted relation using equation (1) and its scatter. The parameters of the fitting are reported in Table 2.</figcaption></figure>

from Brook et al. (2014) and Garrison-Kimmel et al. (2014) by a factor >10 for ‘Doc’ and of about 3–4 for ‘Bashful’ (with respect to Garrison-Kimmel et al. (2014), even more for Brook et al.). This overcooling can be explained by the lack in their simulation of any other source of feedback besides SNe (see for example the results of Stinson et al. 2013). When the higher (may be too high) stellar mass-halo mass ratio of these galaxies, compared to our NIHAO sample, is taken into account the results are in very good agreement.

A similar agreement is also found when comparing our results with recent simulations presented by the FIRE collaboration in Chan et al. (2015). For example there is a quite remarkable similarity in the peak of the core formation and in the non-monotonic response of DM to galaxy formation. This similarity is even more striking when we take into account that the two simulation sets (NIHAO and FIRE) differ in resolution, numerical technique and stellar feedback implementation.

One place where our and their results do not agree is in the decline of the value of $\alpha$ at low masses and which brought them to suggest the possibility of a ‘threshold’ halo mass at around $10^{10} M_{\odot}$ needed to develop large cores. It is worth mentioning that while the FIRE simulations reach a much higher resolution than our current sample, they do not have large number statistics, with only nine galaxies across five orders of magnitude in stellar mass.

We do not see such a threshold in our sample at this halo mass (but note that we do see a substantial change in $\alpha$ for a stellar

<figure><img src="image_8.png" /><figcaption>**Figure 7.** DM density profiles for our three test galaxies at two redshifts. The profiles at $z=0$ are represented in red, and at $z=2$ in green. The dotted line is fitted for each profile between 1 and 2 per cent of the virial radius and is shown up to 10 per cent of the virial radius to help visualize.</figcaption></figure>

mass of $\approx 10^7 M_{\odot}$). This difference can have many explanations. The higher resolution of the FIRE simulations might model the interplay between star formation and the outflows they drive better. The partial evidence for a threshold could be due to the much lower sample size of Chan et al.: they only have three galaxies in the DM range between $5 \times 10^9$ and $5 \times 10^{10} M_{\odot}$, while we have 25, and hence we are less subject to artificial ‘jumps’ due to low number statistics and intrinsic scatter.

Finally as in the case of the Madau et al. results, the FIRE simulations predict a stronger star formation efficiency at low masses, which partially overproduces abundance matching results from Garrison-Kimmel et al. (2014) and Brook et al. (2014). This departure from abundance matching results, starts exactly above halo masses of $10^{10} M_{\odot}$ (Dutton et al. 2015, accepted fig. 1). Overall we would like to stress again the importance of having large samples and realistic stellar-to-halo-mass ratios in order to infer the behaviour of the DM distribution in hydrodynamic simulations.

## 3.1 Single profile evolution

To study the evolution of individual haloes, we make case studies of three simulations, g2.63e10, g2.19e11 and g8.06e11 (already shown in Fig. 2 and presented in Table 1), that show different behaviours as a function of redshift. In Fig. 7, we show the DM density profiles
at $z = 0$ (red points) and $z = 2$ (green points) for the three galaxies. At $z = 2$ the total halo masses are $1.48 \times 10^{10}$, $2.61 \times 10^{10}$ and $1.50 \times 10^{11}$, while the stellar masses are $3.29 \times 10^{7}$, $2.64 \times 10^{7}$ and $1.11 \times 10^{9}$, for g2.63e10, g2.19e11 and g8.06e11, respectively. The dashed line shows a power-law fit to the density between 1 and 2 per cent of the virial radius.

As expected, at $z = 2$ the middle mass galaxy already has an expanded profile (compared to the N-body run), with a slope of $\alpha = -0.65$, and then the inner density flattens by the end of the simulation. At $z = 2$, the DM profile of the higher mass halo looks nearly identical to the middle mass halo at $z = 0$. It evolves differently; its central density profile becomes steeper with time. The low-mass galaxy (g2.63e10) has a very mild evolution and the slope of the profiles does not change much from $z = 2$ to the present day.

This strong time variation of the size of the central density core are correlated to quick variations in the star formation rate of the galaxy (Pontzen & Governato 2012; see also Macciò et al. 2012 for the effect of the sloshing of the gas inside the potential well).

Fig. 8 provides example cases of how bursty star formation effects the inner DM profile slope for galaxies of three different masses spanning the NIHAO sample. The lowest mass case, g2.63e10, has early bursts of star formation before its star formation dies away entirely. The bursts of star formation are not reflected in the inner profile slope at early times, most likely because of the low resolution of the profiles that makes it hard to measure the slope. As time goes on and the resolution increases, the slope stays moderately flatter than NFW, at a value slightly higher than -1.

The middle mass case, g2.19e11, displays bursts of star formation throughout its evolution. The largest bursts start around 7 Gyr ($z \sim 1$). The bursts coincide with a flattening of the inner DM profile slope.

The highest mass galaxy, g8.06e11, shows a monotonically rising star formation history. Throughout the history, there are bursty periods. Its inner density profile evolves from a cusp to a core before contracting back into a cusp with $\alpha = -1.7$ at $z = 0$. It changes from a cusp to a core at the same time that its star formation rate increases by a factor of 4, around 7 Gyr.

## 3.2 Core creation

The variation of the central slope of the density profile is also linked to the creation (and destruction) of a constant density inner core. In order to have an estimate of the size of such a core, we have decided to follow the approach of Macciò et al. (2012) and fit the profiles with the following parametric description, originally presented in Stadel et al. (2009):

$$ \rho(r) = \rho_0 \exp(\lambda[\ln(1 + r/R_\lambda)]^2). \qquad (2) $$

In this parametrization, the density profile is linear down to a scale $R_\lambda$ beyond which it approaches the central maximum density $\rho_0$ as $r \to 0$. We also note that if one makes a plot of $d\ln\rho/d\ln(1 + r/R_\lambda)$ versus $\ln(1 + r/R_\lambda)$ then this profile forms an exact straight line with slope $2\lambda$. This fitting function is extremely flexible and makes it possible to reproduce at the same time both cuspy and cored density profiles. In the following we will use the fitted value of $R_\lambda$ as the core size.

In Fig. 9 we show the core size ($R_\lambda$) as function of time for the two more massive galaxies shown in Fig. 7. The grey band in the lower panel of the plot represents the softening of the simulations, any core value below this limit has to be interpreted as a non-cored profile. The lower mass galaxy (red line) shows a gradually increasing average core size from $z = 2$ to 0. The behaviour of the

<figure><img src="image_9.png" /><figcaption>Figure 8. Time evolution of the star formation rate (SFR, upper panels) and the central density slope ($\alpha$, lower panel), for our three test galaxies: g8.06e11 (blue), g2.19e11 (red) g2.63e10 (orange).</figcaption></figure>

more massive galaxy (blue line) is much more interesting. At high redshift ($z = 2$) the halo has a core of order $\sim 1.5$ kpc; in between $z = 1.5$ and 1 the core size fluctuates a lot and reaches quite large value up to 4 kpc. The core is then definitely destroyed after $z = 0.5$ and the galaxy shows no presence of a central constant density, as confirmed by the profile shown in Fig. 2.

It is interesting to look deeper into the origin of the spikes of the core value around $z \sim 0.5$ ($T \sim 8.5$ Gyr) for the g8.06e11 halo. In Fig. 10 we show the DM density profile of this galaxy at the time of the maximum size of the core ($z = 0.67$, black line), and right before the core creation (red line) and right after the core destruction (blue line). It is clear from the figure that the profile does change quite
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Core size $R_\lambda$ as a function of time $1/(1+z)$ for g2.19e11 (red) and g8.06e11 (blue). The grey shaded region corresponds to the resolution limit of the simulations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Evolution of the DM density profile for the g8.06e11 halo around $z \sim 0.7$. The profile changes significantly with clear evidence for the creation and subsequent destruction of a large density core. The time spanned by the three profiles is of the order of 1 Gyr. The two dotted lines mark the region within which $\alpha$ is computed.</figcaption></figure>

strongly in the time interval spanned by the spike. A large core (up to 4 kpc) is created and then destroyed in less than a Gyr.

This fluctuation is related to the rapid increase of an order of magnitude of star formation from $T = 6.5$ to 8.5 Gyr. The morphological analysis of the galaxy shows that this extended burst is due to disc instability, possibly triggered by a satellite accretion. During this phase a lot of gas is both accreted and ejected from the central region creating the perfect conditions for core creation and destruction (e.g. a la Pontzen & Governato 2012).

## 3.3 Core destruction

Our simulations clearly show that DM halo profiles are not a static entity. On the contrary they change significantly during the formation and evolution of the halo. While all haloes start as being cuspy, as predicted by pure CDM simulations, they can develop a core. Depending on the stellar-to-halo-mass ratio, the core can persist to the present day, or be destroyed, being replaced by a central cusp.

In massive haloes (e.g. with halo mass around $10^{12} M_\odot$) the cusp regeneration is mainly due to strong gas inflow in the central region, which builds up a large stellar body in the centre of the galaxy causing a deepening of the local potential well. We will now analyse these processes in more detail.

In Fig. 11 we show the relative abundance of DM, gas and stars in the inner 2 kpc for our test galaxies. In the case of the lowest

<figure><img src="image_12.png" /><figcaption>**Figure 11.** The contribution of DM (black), gas (blue) and stars (red) to the total mass in the inner 2 kpc for g8.06e11 (upper panel), g2.19e11 (middle panel) and g2.63e10 (lower panel), as a function of time. Note the log y-axis in the lower panel.</figcaption></figure>

mass galaxy (g2.63e10) the central potential is always dominated by DM which accounts for more than ~90 per cent of the total mass in the centre. The lack of evolution in $\alpha$ shown in Fig. 8 is then due to the very 'passive' evolution of the central region of this halo after star formation is over. Nevertheless the final effect of baryons on this low-mass halo is still a net expansion (as shown in Fig. 2) with $\alpha = -0.9$ in the Hydro run with respect to a value of ~-1.6 in the N-body one; despite stars and gas accounting for only ~10 per cent of the mass inside 2 kpc.

For the middle mass galaxy (which retains a cored profile all the way to $z=0$), the potential in the central region is dominated by DM and gas, while stars are subdominant. The amount of gas is rapidly varying with time and this causes quick potential variations and a rapid movement of the centre of mass of the gas compared to that of the DM; both these effects contribute in creating and maintaining a density core.

The situation is quite different for the most massive galaxy: after $T = 8$ Gyr the potential in central region is dominated by stars with the gas (and DM) being highly subdominant. At this point, the deeper potential well caused by the stars is anchoring the global potential and making any possible gas outflow a negligible effect. The DM reacts to the new enhanced central potential by contracting towards the centre. This contraction is causing the core size to drop
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Time evolution of the *total* mass enclosed within different radii, namely 5, 10, 50 and 100 kpc, normalized to the $z=0$ mass. From top to bottom the results for our three test galaxies: g8.06e11, g2.19e11 and g2.63e10. The strong evolution of the inner mass (2 and 5 kpc) for galaxy g8.06e11 is thought to be the cause of the halo recontraction (see text for more details).</figcaption></figure>

below our resolution (Fig. 9) and the profile slope $\alpha$ to drop below -1.5 (Fig. 8).

The deepening of the central potential (due to gas inflow and star formation) in our more massive halo can be seen in Fig. 12 where we plot the evolution of the *total* mass (DM+gas+stars) within different radii for our three galaxies. In the middle mass halo (g2.19e11) the mass is constant practically at all radii after 8 Gyr; this implies a constant gravitational potential. In the higher mass halo (g8.06e11) at large radii (100 kpc, still well within the virial radius) the total mass is also constant, but at small radii (5 and 10 kpc) the mass rapidly increases by a factor of 3 after $t=8$ Gyr. This higher mass caused the gravitational potential to deepen and reduces the ability of baryonic feedback to counteract the DM contraction.

It is worth noticing that the change in the mass distribution is mainly local (i.e. limited to the regions dominated by baryons) and not global, something that should be taken into account when computing analytic estimates of the total energy input from SN needed to alter the DM distribution. This is a still quite debated topic in the literature and our simulations can be used to determine the validity of different analytic approaches (e.g. Peñarrubia et al. 2012; Maxwell, Wadsley & Couchman 2015).

## 4 COMPARISON WITH OBSERVATIONS

There is a very large amount of literature on the DM density profile in observed galaxies (e.g Moore 1994; Salucci & Burkert 2000; Dutton et al. 2005; Simon et al. 2005; de Blok et al. 2008; Kuzio de Naray et al. 2008, 2009; Oh et al. 2011a, 2015). Despite the quite large scatter in the slope of the observed profiles there is a clear indication for $\alpha_{obs} > -1$ (e.g. Oh et al. 2011a), which is in contrast with naive expectations from collisionless simulations. This disagreement is strongly alleviated when the effects of baryons are taken into account, as recently shown in Karukes, Salucci & Gentile (2015), where the authors compared the measured DM density profile in NGC 3198 with the predictions from the DC14 model. Thanks to the large number of simulated galaxies in the NIHAO project we can for the first time extend previous observation-simulation comparisons that were based on a limited number of objects (e.g. Oh et al. 2011b).

In Fig. 13 we show the inner slope of the DM density profile $\alpha$ versus the radius $R_{in}$ of the innermost point within which $\alpha$ is measured for both observations and simulations. The observational data are split in three groups, the grey symbols are results from dwarf and low surface brightness (LSB) galaxies (de Blok et al. 2001; de Blok & Bosma 2002; Swaters et al. 2003), cyan filled squares are the results from the THINGS survey (Oh et al. 2011a), while open blue squares with error-bars are the results from the LITTLE THINGS survey (Oh et al. 2015). The three lines show theoretical predictions from a given profile: the dotted line is for an isothermal profile with a core size of 1 kpc, the solid line is for the Einasto profile (Einasto 1965; Merritt et al. 2005) with a concentration of 14 as expected for DM haloes around $10^{11} M_{\odot}$.

<figure><img src="image_14.png" /><figcaption>**Figure 13.** The inner slope of the DM density profile $\alpha$ versus the radius $R_{in}$ of the innermost point within which $\alpha$ is measured. Grey symbols are results from LSB galaxies (de Blok et al. 2001; de Blok & Bosma 2002; Swaters et al. 2003), cyan filled squares are the results from the THINGS survey (Oh et al. 2011a), open blue squares with error-bars are the results from the LITTLE THINGS survey (Oh et al. 2015). The red points are all NIHAO galaxies with stellar masses between $10^7$ and $10^{10} M_{\odot}$ to mimic the range of the observational points. The isothermal profile line is for a core of 1 kpc, while the NFW and Einasto lines are for a concentration of 14 (Dutton & Macciò 2014).</figcaption></figure>
(Dutton & Macciò 2014) and the dashed line is for the NFW profile with the same concentration as the Einasto one. It is interesting to note that on the radial scales of this plot, the Einasto profile (which is known to provide a better fit to simulations; Merritt et al. 2005; Navarro et al. 2010) predicts a steeper slope for DM haloes with respect to NFW.

The red circles with error-bars are the predictions for $\alpha$ from the NIHAO simulation suite at 1–2 per cent of $R_{vir}$, for galaxies with stellar masses between $10^7$ and $10^{10}$ M$_\odot$ to be consistent with the LITTLE THINGS sample range (Oh et al. 2015). All the NIHAO dark matter haloes are expanded with respect to collisionless simulation expectations. We do not have enough resolution to probe the inner part of the radial profile at the same extent as the LITTLE THING survey; nevertheless our results are in very good agreement with the observations where they overlap in $R_{in}$. Moreover thanks to the large number of simulated galaxies, we are also able to reproduce the observed scatter in $\alpha$ at a fixed radius.

This plot further demonstrates the importance of considering the effect of baryons on the DM distribution when comparing observations and simulations based on the CDM model.

## 5 CONCLUSIONS

We have used the NIHAO simulation suite (Wang et al. 2015) to study the effects of baryons on the radial distribution of DM in collapsed haloes. In comparison with a similar study performed earlier with the MaGICC sample (DC14) we have a larger and more coherent set of simulations, which have been performed with the latest compilation of cosmological parameters from the Planck Collaboration XVI (2014). Moreover we use a substantially improved version of the GASOLINE code, which thanks to a new formulation of the SPH force computation (Keller et al. 2014), aims to solve the known shortcomings of the classical SPH implementation (Agertz et al. 2007; Hopkins 2013).

Overall we confirm previous findings from DC14, namely that the effect of baryons on the slope ($\alpha$) of the central distribution of DM strongly depends on the ratio between the stellar mass and total mass of the halo. At the present time, haloes with an intermediate ratio of stellar-to-halo mass ($10^{-3} < M_{star}/M_{halo} < 10^{-2}$) tend to expand their profiles with respect to pure collisional simulations, haloes with lower values of $M_{star}/M_{halo}$ tend to retain the original DM cuspy profile, while haloes with large ($M_{star}/M_{halo} > 10^{-2}$) ratio exhibit contracted profiles. There is a smooth transition among these different behaviours and we have parametrized it with a new, simple but accurate fitting formula. The slopes predicted by the NIHAO simulation suite are in excellent agreement with high-resolution observations of dwarf galaxies from the THINGS and LITTLE THINGS surveys (Oh et al. 2015).

When the analysis is, for the first time, extended to higher redshifts, we find that haloes tend to cluster around the redshift zero $\alpha$ versus $M_{star}/M_{halo}$ relation, even though the scatter around the mean increases with redshift. Our findings imply that even haloes that have very cuspy profiles today (as for example for Milky Way mass haloes), likely had a cored DM distribution in the past.

We carefully looked at three specific haloes with masses around $10^{10}$, $10^{11}$ and $10^{12}$ M$_\odot$. The lightest shows a central slope of about $-0.9$ at $z=0$, which is consistent with a cuspy profile, but is also very much shallower than what a pure N-body simulation predicts at the same mass and radii ($\alpha \sim -1.6$), pointing to a clear halo expansion even at these low-mass scales. The density profile was even shallower at higher redshift ($z \sim 1$) when the halo went through several star formation episodes.

The other two haloes present a quite dynamic evolution in their central regions as a function of time, showing that not only cores are created and destroyed during the galaxy evolution but that this phenomenon happens on relatively (cosmologically) short time-scales of the order of a Gyr.

Overall our results show that one of the most firm predictions of the CDM theory, namely the existence of a ‘universal profile’, is no longer valid when galaxy formation is taken into account. The central part of any DM halo has had a much more tumultuous life than what dissipationless simulations predict.

## ACKNOWLEDGEMENTS

It is a pleasure to thank Se-Heon Oh, for sending us his data in electronic form. AVM, AAD, GSS, CP, TAG and TB acknowledge support from the Sonderforschungsbereich SFB 881 ‘The Milky Way System’ (subprojects A1 and A2) of the German Research Foundation (DFG). The simulations were performed on the THEO cluster of the Max-Planck-Institut für Astronomie and the HYDRA cluster at the Rechenzentrum in Garching; and the MILKY WAY supercomputer, funded by the Deutsche Forschungsgemeinschaft (DFG) through Collaborative Research Center (SFB 881) ‘The Milky Way System’ (subproject Z2), hosted and co-funded by the Jülich Supercomputing Center (JSC). We greatly appreciate the contributions of all these computing allocations. The authors acknowledge support from the MPG-CAS through the partnership programme between the MPIA group led by AVM and the PMO group led by XK. LW acknowledges support of the MPG-CAS student programme. XK is supported by the NSFC (no. 11333008) and the Strategic Priority Research Program the emergence of cosmological structure of the CAS (no. XDB09000000).

## REFERENCES

* Agertz O. et al., 2007, MNRAS, 380, 963

* Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

* Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

* Bertschinger E., 2001, ApJS, 137, 1

* Brook C. B., Di Cintio A., Knebe A., Gottlöber S., Hoffman Y., Yepes G., Garrison-Kimmel S., 2014, ApJ, 784, L14

* Chabrier G., 2003, PASP, 115, 763

* Chan T. K., Kereš D., Ónóorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 454, 2981

* de Blok W. J. G., Bosma A., 2002, A&A, 385, 816

* de Blok W. J. G., McGaugh S. S., Bosma A., Rubin V. C., 2001, ApJ, 552, L23

* de Blok W. J. G., Walter F., Brinks E., Trachernach C., Oh S.-H., Kennicutt R. C., Jr, 2008, AJ, 136, 2648

* Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014a, MNRAS, 437, 415

* Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Stinson G. S., Knebe A., 2014b, MNRAS, 441, 2986 (DC14)

* Diemand J., Moore B., Stadel J., 2004, MNRAS, 353, 624

* Dutton A. A., Macciò A. V., Frings J., Wang L., Stinson G. S., Penzo C., Kang X., 2015, accepted to MNRAS

* Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359

* Dutton A. A., Courteau S., de Jong R., Carignan C., 2005, ApJ, 619, 218

* Einasto J., 1965, Tr. Astrofizicheskogo Inst. Alma-Ata, 5, 87

* Gao L., Navarro J. F., Cole S., Frenk C. S., White S. D. M., Springel V., Jenkins A., Neto A. F., 2008, MNRAS, 387, 536

* Gill S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 399

* Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16

* Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Lee K., 2014, MNRAS, 438, 2578
Governato F. et al., 2010, Nature, 463, 203

Hopkins P. F., 2013, Astrophysics Source Code Library, record ascl:1305.006

Karukes E. V., Salucci P., Gentile G., 2015, A&A, 578, A13

Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013

Klypin A., Kravtsov A. V., Bullock J. S., Primack J. R., 2001, ApJ, 554, 903

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Kravtsov A., Vikhlinin A., Meshscheryakov A., 2014, ApJ, preprint (arXiv:1401.7329)

Kuzio de Naray R., McGaugh S. S., de Blok W. J. G., 2008, ApJ, 676, 920

Kuzio de Naray R., McGaugh S. S., Mihos J. C., 2009, ApJ, 692, 1321

Macciò A. V., Dutton A. A., van den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55

Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940

Macciò A. V., Stinson G., Brook C. B., Wadsley J., Couchman H. M. P., Shen S., Gibson B. K., Quinn T., 2012, ApJ, 744, L9

Madau P., Shen S., Governato F., 2014, ApJ, 789, L17

Martizzi D., Teyssier R., Moore B., 2013, MNRAS, 432, 1947

Mashchenko S., Couchman H. M. P., Wadsley J., 2006, Nature, 442, 539

Maxwell A. J., Wadsley J., Couchman H. M. P., 2015, ApJ, 806, 229

Merritt D., Navarro J. F., Ludlow A., Jenkins A., 2005, ApJ, 624, L85

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Moore B., 1994, Nature, 370, 629

Munshi F. et al., 2013, ApJ, 766, 56

Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, L72

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 (NFW)

Navarro J. F. et al., 2010, MNRAS, 402, 21

Neto A. F. et al., 2007, MNRAS, 381, 1450

Ogiya G., Mori M., 2014, ApJ, 793, 46

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt R. C., Jr, 2011a, AJ, 141, 193

Oh S.-H., Brook C., Governato F., Brinks E., Mayer L., de Blok W. J. G., Brooks A., Walter F., 2011b, AJ, 142, 24

Oh S.-H. et al., 2015, AJ, 149, 180

Oñ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, 2092

Peñarrubia J., Pontzen A., Walker M. G., Koposov S. E., 2012, ApJ, 759, L42

Penzo C., Macciò A. V., Casarini L., Stinson G. S., Wadsley J., 2014, MNRAS, 442, 176

Planck Collaboration XVI, 2014, A&A, 571, A16

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Pontzen A., Governato F., 2014, Nature, 506, 171

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Prada F., Klypin A. A., Cuesta A. J., Betancort-Rijo J. E., Primack J., 2012, MNRAS, 423, 3018

Read J. I., Gilmore G., 2005, MNRAS, 356, 107

Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743

Saitoh T. R., Makino J., 2009, ApJ, 697, L99

Salucci P., Burkert A., 2000, ApJ, 537, L9

Schaller M. et al., 2015, MNRAS, 451, 1247

Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581

Simon J. D., Bolatto A. D., Leroy A., Blitz L., Gates E. L., 2005, ApJ, 621, 757

Springel V. et al., 2005, Nature, 435, 629

Stadel J. G., 2001, PhD thesis, Univ. Washington

Stadel J., Potter D., Moore B., Diemand J., Madau P., Zemp M., Kuhlen M., Quilis V., 2009, MNRAS, 398, L21

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074

Stinson G. S. et al., 2012, MNRAS, 425, 1270

Stinson G. S. et al., 2013, MNRAS, 436, 625

Stinson G. S. et al., 2015, MNRAS, 454, 1105

Swaters R. A., Madore B. F., van den Bosch F. C., Balcells M., 2003, ApJ, 583, 732

Trujillo-Gomez S., Klypin A., Colín P., Ceverino D., Arraki K. S., Primack J., 2015, MNRAS, 446, 1140

Vogelsberger M., Zavala J., Simpson C., Jenkins A., 2014, MNRAS, 444, 3684

Wadsley J. W., Stadel J., Quinn T., 2004, New Astron., 9, 137

Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS, 387, 427

Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83

This paper has been typeset from a *TEX/LATex* file prepared by the author.