- Split View
-
Views
-
Cite
Cite
Aaron D. Ludlow, Sownak Bose, Raúl E. Angulo, Lan Wang, Wojciech A. Hellwing, Julio F. Navarro, Shaun Cole, Carlos S. Frenk, The mass–concentration–redshift relation of cold and warm dark matter haloes, Monthly Notices of the Royal Astronomical Society, Volume 460, Issue 2, 01 August 2016, Pages 1214–1232, https://doi.org/10.1093/mnras/stw1046
- Share Icon Share
Abstract
We use a suite of cosmological simulations to study the mass–concentration–redshift relation, c(M, z), of dark matter haloes. Our simulations include standard Λ-cold dark matter (CDM) models, and additional runs with truncated power spectra, consistent with a thermal warm dark matter (WDM) scenario. We find that the mass profiles of CDM and WDM haloes are self-similar and well approximated by the Einasto profile. The c(M, z) relation of CDM haloes is monotonic: concentrations decrease with increasing virial mass at fixed redshift, and decrease with increasing redshift at fixed mass. The mass accretion histories (MAHs) of CDM haloes are also scale-free, and can be used to infer concentrations directly. These results do not apply to WDM haloes: their MAHs are not scale-free because of the characteristic scale imposed by the power spectrum suppression. Further, the WDM c(M, z) relation is non-monotonic: concentrations peak at a mass scale dictated by the truncation scale, and decrease at higher and lower masses. We show that the assembly history of a halo can still be used to infer its concentration, provided that the total mass of its progenitors is considered (the ‘collapsed mass history’; CMH), rather than just that of its main ancestor. This exploits the scale-free nature of CMHs to derive a simple scaling that reproduces the mass–concentration–redshift relation of both CDM and WDM haloes over a vast range of halo masses and redshifts. Our model therefore provides a robust account of the mass, redshift, cosmology and power spectrum dependence of dark matter halo concentrations.
1 INTRODUCTION
The scaling parameters of the NFW and Einasto profiles can be expressed in alternative forms, such as halo virial1 mass, M200, and concentration, c ≡ r200/r−2, defined as the ratio of the virial radius to that of the scale radius. At a given halo mass, the concentration provides an alternative measure of the characteristic density of a halo.
As discussed by NFW, M200 and c do not take on arbitrary values, but correlate in a way that reflects the mass dependence of halo formation times: those that assemble earlier have higher characteristic densities, reflecting the larger background density at that epoch. They used this finding to build a simple analytic model based on the extended Press-Schechter (EPS) theory (Bond et al. 1991) that reproduced the average mass and cosmology dependence of halo concentrations in their early simulations.
Subsequent work by Bullock et al. (2001) corroborated the general trends reported by NFW, but underscored a much stronger redshift dependence of the concentration–mass relation than expected from their model. These authors proposed a revised model that predicts concentrations which, at fixed mass, scale linearly with expansion factor (c ∝ (1 + z)−1) and, at fixed z, fall off rapidly with increasing mass. Later numerical work, however, found a much weaker mass and redshift dependence than predicted by this model. Most notably, the concentrations of very massive haloes are found to be approximately constant and to evolve little with redshift (Zhao et al. 2003, 2009; Gao et al. 2008).
Empirical models that link halo concentrations to the shape of their assembly histories fare better. The models of Wechsler et al. (2002) and Zhao et al. (2003), for example, assume that the concentration is set by the changing accretion rate of a halo, with the characteristic density tracing the time when haloes transition from an initial rapid accretion phase to a subsequent phase of slower growth. Very massive systems are still in their initial rapid-growth phase at present, thus explaining why they all have similar concentrations.
More recently, Ludlow et al. (2013, hereafter L13) used the Millennium simulations (hereafter MS) to investigate the connection between cold dark matter (CDM) halo assembly and structure. They pointed out that halo mass profiles and main-progenitor mass accretion histories (MAHs) are self-similar and have similar shapes. This becomes apparent when expressing mass profiles in terms of average enclosed density, M(〈ρ〉(r)), rather than radius, and MAHs in terms of main progenitor mass as a function of cosmic density rather than time or redshift, i.e. M(ρcrit(z)). Both follow the NFW profile. The two M(ρ) functions can thus be linked by a simple scaling that allows the characteristic density of a halo (or its concentration) to be inferred from the critical density of the Universe at a characteristic time along its MAH. Ludlow et al. (2014a, hereafter L14) showed how this result can be used to build a simple analytic model for the mass–concentration–redshift relation that accurately reproduced the trends obtained for CDM haloes in a large number of simulations, as well as the cosmology dependence of c(M, z) reported in previous work.
Although the model works well for CDM (see e.g. Correa et al. 2015c), its applicability to models with truncated power spectra, such as those expected for ‘warm’ dark matter (WDM), is unclear. Interest in such models has been revitalized by recent claims of detection of an ∼ 3.5 keV X-ray line, which is in principle consistent with the decay of a WDM particle (e.g. Bulbul et al. 2014; Boyarsky et al. 2014, 2015, but see, Malyshev, Neronov & Eckert 2014; Anderson, Churazov & Bregman 2015). These results have motivated observations of DM-dominated dwarf galaxies of the Local Group (Lovell et al. 2015) which, to date, have not provided compelling evidence for DM decay (Jeltema & Profumo 2016).
The structure of WDM haloes has been studied by Macciò et al. (2013) and, more recently, by Bose et al. (2016), who report that the NFW formula accounts well for their mass profile shape (see also e.g. Knebe et al. 2002; Villaescusa-Navarro & Dalal 2011; Polisensky & Ricotti 2015; González-Samaniego, Avila-Reese & Colín 2016). The resulting c(M) relation, however, is non-monotonic: concentrations reach a maximum at a halo mass about two decades above the truncation scale but decline gradually towards larger and smaller masses (see also e.g. Schneider et al. 2012; Macciò et al. 2013). This implies that fairly massive WDM haloes can be as concentrated as low-mass ones although the shape of their MAHs differ strongly, a feature that is difficult to reconcile with the MAH-based scenario discussed above for CDM.
Indeed, there have been relatively few attempts to model the c(M, z) relation for truncated power spectra. One exception is the work of Eke, Navarro & Steinmetz (2001), who assumed that both the normalization and shape of the power spectrum modulates c(M, z). By including a term proportional to d ln σ/d ln M in the definition of the collapse time they were able to reproduce the concentration–mass relation in both CDM models as well as several with truncated power spectra. This particular parametrization, however, does not lend itself to simple interpretation and predictions of their model were not borne out by more recent simulations (e.g. Gao et al. 2008; Diemer & Kravtsov 2015). Schneider (2015) provides empirical relations that can be used to map c(M, z) relations obtained for CDM haloes to those expected for warm or mixed DM models.
It is clear that a full picture of DM halo structure must address not only the mass and cosmological parameter dependence of halo concentrations, but also the effect of the initial density fluctuation spectrum. As more and more observations sensitive to the small-scale clustering of DM become available, theoretical tools such as these will be indispensable. This is the focus of the current paper. Using a large suite of CDM and WDM cosmological simulations we study the relationship between halo assembly and structure, paying particular attention to the signature of a potential WDM particle.
Our paper is structured as follows. We begin with a description of our simulations in Section 2; their analysis (including halo finding, merger tree construction and density profile estimates) are outlined in Section 3. In Section 4 we present our main results. The mass–concentration–redshift relation and its connection to halo assembly are presented in Sections 4.1 and 4.2. We then use these results to build a simple analytic model for c(M, z) in Section 4.3, which is compared to other available models in Sections 4.4 and 4.5. Finally, in Section 5, we provide a summary of our findings. We elaborate on various aspect of our c(M, z) model in Appendices A and B and provide an accurate fitting function for CDM haloes in the Planck cosmology in Appendix C.
2 NUMERICAL SIMULATIONS
Our analysis focuses on the growth histories and internal structure of collisionless DM haloes identified in a suite of cosmological numerical simulations. The majority of our results are based on the Copernicus Complexio (coco) simulations (Bose et al. 2016; Hellwing et al. 2016), supplemented by the Millennium (Springel et al. 2005; Boylan-Kolchin et al. 2009; Angulo et al. 2012) and Aquarius simulations (Springel et al. 2008; Lovell et al. 2014), and an additional suite of ΛCDM runs which vary the parameters of the background cosmological model.2 The main aspects of these models are detailed in Tables 1 and 2. We provide here a brief description of the runs, but refer the reader to the original papers for a more thorough discussion.
Simulation . | Model . | Vbox . | Np . | ϵ . | mp . | mWDM . | Mhm . |
---|---|---|---|---|---|---|---|
. | . | (h−3 Mpc3) . | . | (h−1 kpc) . | (h−1 M⊙) . | (keV) . | (1010 h−1 M⊙) . |
coco | CDM/WDM | 2.2 × 104 | 23443 | 0.23 | 1.135 × 105 | 3.3 | 0.025 |
color-1.5 | WDM | 3.5 × 105 | 16203 | 1 | 6.20 × 106 | 1.5 | 0.34 |
MS-II | CDM | 1 × 106 | 21603 | 1 | 6.89 × 106 | – | – |
MS-I | CDM | 1.3 × 108 | 21603 | 5 | 8.61 × 108 | – | – |
MS-XXL | CDM | 2.7 × 1010 | 67203 | 10 | 6.17 × 109 | – | – |
Aquarius | CDM/WDM | – | ∼8093 | 0.05 | 1.09 × 104 | 1.5, 1.6, 2.0, 2.3 | 0.34, 0.28, 0.13, 0.08 |
Cosmo-A | CDM | 5.2 × 107 | 10803 | 7.5 | 1.72 × 109 | – | – |
Cosmo-B | CDM | 8.6 × 106 | 10803 | 4.1 | 4.78 × 108 | – | – |
Cosmo-C | CDM | 1.1 × 107 | 10803 | 4.5 | 7.22 × 108 | – | – |
Cosmo-D | CDM | 5.5 × 106 | 10803 | 3.5 | 4.84 × 108 | – | – |
Simulation . | Model . | Vbox . | Np . | ϵ . | mp . | mWDM . | Mhm . |
---|---|---|---|---|---|---|---|
. | . | (h−3 Mpc3) . | . | (h−1 kpc) . | (h−1 M⊙) . | (keV) . | (1010 h−1 M⊙) . |
coco | CDM/WDM | 2.2 × 104 | 23443 | 0.23 | 1.135 × 105 | 3.3 | 0.025 |
color-1.5 | WDM | 3.5 × 105 | 16203 | 1 | 6.20 × 106 | 1.5 | 0.34 |
MS-II | CDM | 1 × 106 | 21603 | 1 | 6.89 × 106 | – | – |
MS-I | CDM | 1.3 × 108 | 21603 | 5 | 8.61 × 108 | – | – |
MS-XXL | CDM | 2.7 × 1010 | 67203 | 10 | 6.17 × 109 | – | – |
Aquarius | CDM/WDM | – | ∼8093 | 0.05 | 1.09 × 104 | 1.5, 1.6, 2.0, 2.3 | 0.34, 0.28, 0.13, 0.08 |
Cosmo-A | CDM | 5.2 × 107 | 10803 | 7.5 | 1.72 × 109 | – | – |
Cosmo-B | CDM | 8.6 × 106 | 10803 | 4.1 | 4.78 × 108 | – | – |
Cosmo-C | CDM | 1.1 × 107 | 10803 | 4.5 | 7.22 × 108 | – | – |
Cosmo-D | CDM | 5.5 × 106 | 10803 | 3.5 | 4.84 × 108 | – | – |
Simulation . | Model . | Vbox . | Np . | ϵ . | mp . | mWDM . | Mhm . |
---|---|---|---|---|---|---|---|
. | . | (h−3 Mpc3) . | . | (h−1 kpc) . | (h−1 M⊙) . | (keV) . | (1010 h−1 M⊙) . |
coco | CDM/WDM | 2.2 × 104 | 23443 | 0.23 | 1.135 × 105 | 3.3 | 0.025 |
color-1.5 | WDM | 3.5 × 105 | 16203 | 1 | 6.20 × 106 | 1.5 | 0.34 |
MS-II | CDM | 1 × 106 | 21603 | 1 | 6.89 × 106 | – | – |
MS-I | CDM | 1.3 × 108 | 21603 | 5 | 8.61 × 108 | – | – |
MS-XXL | CDM | 2.7 × 1010 | 67203 | 10 | 6.17 × 109 | – | – |
Aquarius | CDM/WDM | – | ∼8093 | 0.05 | 1.09 × 104 | 1.5, 1.6, 2.0, 2.3 | 0.34, 0.28, 0.13, 0.08 |
Cosmo-A | CDM | 5.2 × 107 | 10803 | 7.5 | 1.72 × 109 | – | – |
Cosmo-B | CDM | 8.6 × 106 | 10803 | 4.1 | 4.78 × 108 | – | – |
Cosmo-C | CDM | 1.1 × 107 | 10803 | 4.5 | 7.22 × 108 | – | – |
Cosmo-D | CDM | 5.5 × 106 | 10803 | 3.5 | 4.84 × 108 | – | – |
Simulation . | Model . | Vbox . | Np . | ϵ . | mp . | mWDM . | Mhm . |
---|---|---|---|---|---|---|---|
. | . | (h−3 Mpc3) . | . | (h−1 kpc) . | (h−1 M⊙) . | (keV) . | (1010 h−1 M⊙) . |
coco | CDM/WDM | 2.2 × 104 | 23443 | 0.23 | 1.135 × 105 | 3.3 | 0.025 |
color-1.5 | WDM | 3.5 × 105 | 16203 | 1 | 6.20 × 106 | 1.5 | 0.34 |
MS-II | CDM | 1 × 106 | 21603 | 1 | 6.89 × 106 | – | – |
MS-I | CDM | 1.3 × 108 | 21603 | 5 | 8.61 × 108 | – | – |
MS-XXL | CDM | 2.7 × 1010 | 67203 | 10 | 6.17 × 109 | – | – |
Aquarius | CDM/WDM | – | ∼8093 | 0.05 | 1.09 × 104 | 1.5, 1.6, 2.0, 2.3 | 0.34, 0.28, 0.13, 0.08 |
Cosmo-A | CDM | 5.2 × 107 | 10803 | 7.5 | 1.72 × 109 | – | – |
Cosmo-B | CDM | 8.6 × 106 | 10803 | 4.1 | 4.78 × 108 | – | – |
Cosmo-C | CDM | 1.1 × 107 | 10803 | 4.5 | 7.22 × 108 | – | – |
Cosmo-D | CDM | 5.5 × 106 | 10803 | 3.5 | 4.84 × 108 | – | – |
Model . | Ωbar . | ΩM . | |$\Omega _\Lambda$| . | h . | σ8 . | ns . |
---|---|---|---|---|---|---|
Planck | 0.0484 | 0.308 | 0.692 | 0.678 | 0.815 | 0.968 |
WMAP-7 | 0.0446 | 0.272 | 0.728 | 0.704 | 0.81 | 0.967 |
WMAP-1 | 0.045 | 0.25 | 0.75 | 0.73 | 0.9 | 1.0 |
Cosmo-A | 0.045 | 0.15 | 0.85 | 0.73 | 1.0 | 1.0 |
Cosmo-B | 0.045 | 0.25 | 0.75 | 0.73 | 0.6 | 1.0 |
Cosmo-C | 0.045 | 0.29 | 0.71 | 0.73 | 0.81 | 1.0 |
Cosmo-D | 0.045 | 0.40 | 0.60 | 0.73 | 0.7 | 1.0 |
Model . | Ωbar . | ΩM . | |$\Omega _\Lambda$| . | h . | σ8 . | ns . |
---|---|---|---|---|---|---|
Planck | 0.0484 | 0.308 | 0.692 | 0.678 | 0.815 | 0.968 |
WMAP-7 | 0.0446 | 0.272 | 0.728 | 0.704 | 0.81 | 0.967 |
WMAP-1 | 0.045 | 0.25 | 0.75 | 0.73 | 0.9 | 1.0 |
Cosmo-A | 0.045 | 0.15 | 0.85 | 0.73 | 1.0 | 1.0 |
Cosmo-B | 0.045 | 0.25 | 0.75 | 0.73 | 0.6 | 1.0 |
Cosmo-C | 0.045 | 0.29 | 0.71 | 0.73 | 0.81 | 1.0 |
Cosmo-D | 0.045 | 0.40 | 0.60 | 0.73 | 0.7 | 1.0 |
Model . | Ωbar . | ΩM . | |$\Omega _\Lambda$| . | h . | σ8 . | ns . |
---|---|---|---|---|---|---|
Planck | 0.0484 | 0.308 | 0.692 | 0.678 | 0.815 | 0.968 |
WMAP-7 | 0.0446 | 0.272 | 0.728 | 0.704 | 0.81 | 0.967 |
WMAP-1 | 0.045 | 0.25 | 0.75 | 0.73 | 0.9 | 1.0 |
Cosmo-A | 0.045 | 0.15 | 0.85 | 0.73 | 1.0 | 1.0 |
Cosmo-B | 0.045 | 0.25 | 0.75 | 0.73 | 0.6 | 1.0 |
Cosmo-C | 0.045 | 0.29 | 0.71 | 0.73 | 0.81 | 1.0 |
Cosmo-D | 0.045 | 0.40 | 0.60 | 0.73 | 0.7 | 1.0 |
Model . | Ωbar . | ΩM . | |$\Omega _\Lambda$| . | h . | σ8 . | ns . |
---|---|---|---|---|---|---|
Planck | 0.0484 | 0.308 | 0.692 | 0.678 | 0.815 | 0.968 |
WMAP-7 | 0.0446 | 0.272 | 0.728 | 0.704 | 0.81 | 0.967 |
WMAP-1 | 0.045 | 0.25 | 0.75 | 0.73 | 0.9 | 1.0 |
Cosmo-A | 0.045 | 0.15 | 0.85 | 0.73 | 1.0 | 1.0 |
Cosmo-B | 0.045 | 0.25 | 0.75 | 0.73 | 0.6 | 1.0 |
Cosmo-C | 0.045 | 0.29 | 0.71 | 0.73 | 0.81 | 1.0 |
Cosmo-D | 0.045 | 0.40 | 0.60 | 0.73 | 0.7 | 1.0 |
Note that in each of our WDM simulations, we can safely neglect the intrinsic thermal velocities of the particles which, at z = 0, are of the order of a few tens of m s−1. These particles will free-stream only a few kiloparsecs over a Hubble time, which is comparable to our interparticle separation (see Lovell et al. 2012).
2.1 The Copernicus Complexio simulations
The coco simulations track the evolution of DM in an approximately spherical high-resolution volume of radius ∼ 18 h−1Mpc embedded within a lower resolution periodic box of side-length 70.4 h−1Mpc. The high-resolution region contains approximately 13 billion particles and was chosen in order to provide a cosmologically representative sample of Milky Way-mass haloes whilst avoiding the unnecessary computational overhead of including substantially more massive systems. To this end, the high-resolution region was selected so that: (1) it includes no haloes more massive than 5 × 1013 h−1 M⊙; (2) has no haloes more massive than 5 × 1014 h−1 M⊙ within ∼5 h−1Mpc of its boundary, and (3) has a number density of ∼1012 h−1 M⊙ haloes that matches the universal halo mass function.
Linear perturbations were generated at z = 127 using second-order Lagrangian perturbation theory (Jenkins 2013) assuming a standard ΛCDM power spectrum, as well as with a truncated power spectrum compatible with a 3.3 keV thermal WDM particle. We will hereafter refer to these runs as coco-cold and coco-warm, respectively. Both simulations have identical phases and cosmological parameters, the latter adopting values consistent with the WMAP 7-year data release (Komatsu et al. 2011): ΩM = 0.272, ΩΛ = 0.728, σ8 = 0.81, h = 0.704 and ns = 0.967. Here Ωi is the present-day contribution to the energy density from component i; σ8 the linearly extrapolated rms mass-fluctuation in spheres of 8 h−1 Mpc; h is the current Hubble expansion rate in units of 100 km s− 1 Mpc− 1; and ns is the primordial spectral index of density perturbations. With these choices of cosmological and numerical parameters, the high-resolution particle-mass in the coco runs is mp = 1.135 × 105 h−1 M⊙.
We have also run a lower resolution version of the coco-warm simulation assuming a lighter WDM particle of mass 1.5 keV. We will refer to this run as color-1.5 (coco-low resolution, 1.5 keV). This run adopts the same set of WMAP-7 cosmological parameters, but samples the full 70.4 h−1Mpc box with 16203 particles of equal mass, mp = 6.2 × 106 h−1 M⊙. We will use this run to assess the effect of changing the thermal cut-off in the DM power spectrum on the internal structure of DM haloes.
2.2 The millennium and aquarius simulations
Because of the relatively small volume of the coco simulations, we will extend the dynamic range of our analysis using the MS suite. Each of these runs adopts cosmological parameters which were chosen to match the WMAP-1 ΛCDM values – ΩM = 0.25; ΩΛ = 0.75; σ8 = 0.9; h = 0.73; ns = 1 – but differ in both total particle number and in box size. The Millennium (Springel et al. 2005) and Millennium-II (Boylan-Kolchin et al. 2009) simulations evolve the DM density field using Np = 21603 particles in periodic boxes with side-lengths equal to 500 and 100 h−1 Mpc, respectively. The Millennium-XXL simulation (Angulo et al. 2012) adopts both a larger particle number, Np = 67203, and box size, Lbox = 3 h−1Gpc, than either MS or MS-II.
We also use the latest suite of CDM and WDM simulations from the Aquarius Project (Lovell et al. 2014). Like coco, the Aquarius simulations assumed a WMAP-7-normalized power spectrum, but focused computational resources on a single Milky Way-mass DM halo and its surroundings. Each run has the same high-resolution particle mass, mp = 1.09 × 104 h−1 M⊙ (equivalent to level-2 in the original nomenclature of the Aquarius Project), and identical phases of the initial Gaussian random field, but adopt transfer functions appropriate for CDM and thermal WDM models of mass mWDM = 2.3, 2.0, 1.6 and 1.5 keV.
2.3 Additional runs
We have also carried out four additional flat ΛCDM simulations which vary the parameters of the background cosmological model. Each run uses 10803 equal-mass particles, assumes h = 0.73 and ns = 1, but varies the matter density parameter, ΩM, and rms fluctuation amplitude, σ8. The cosmological parameters of these runs, which we have labelled Cosmo-A, B, C and D, are provided in Table 2.
3 HALO INVENTORY AND ANALYSIS TECHNIQUES
3.1 Halo identification
All simulation outputs were processed with a friends-of-friends (FoF) group finder (Davis et al. 1985) using a linking-length of b = 0.2 times the Lagrangian interparticle separation; at each snapshot, groups with fewer than 20 particles were discarded. The substructure finder subfind (Springel et al. 2001) was then run on the remaining FoF groups in order to identify their self-bound subhaloes. subfind decomposes each FoF group into a dominant (or central) subhalo and a contingent of less-massive subhaloes that trace the self-bound relics of past accretion events. For simplicity, we will refer to the assemblage of each central halo and its full subhalo population as a ‘main halo’.
For each main halo, subfind records a virial mass, M200, and associated radius, r200. For our analysis, we will retain only main haloes that exist as distinct objects at z0 = 0, 1, 2 or 3 and, additionally, contain at least N200 = 5000 particles within their virial radius.
3.2 Assembly histories
The halo catalogue is used to construct merger trees for each main halo following the procedure described in Jiang et al. (2014). This method tracks particles within each subhalo across simulation outputs in order to determine their descendants. Subhaloes and their descendants are then split into unique branches, with new branches growing when a subhalo first appears in the simulation and continuing until it has fully merged with a more massive system. The merger tree of a particular halo is then constructed by packaging the subfind merger trees of each of its surviving subhaloes.
Using these merger trees we construct MAHs for each halo, defined as the evolution of the virial mass, M200(z), of its main progenitor (hereafter MAH). The left-hand panel of Fig. 1 shows, for haloes identified at z0 = 0, the median MAHs computed in three separate mass bins. Solid (blue) curves correspond to coco-cold; dot–dashed (orange) curves to coco-warm. Note that the MAHs of CDM and WDM haloes differ strongly below the characteristic mass scale imposed by the free-streaming of the WDM particle, shown here as a horizontal grey line at the ‘half-mode’ mass.3
Although the MAH provides a useful proxy for the assembly history of a halo, it neglects the full spectrum of progenitors that contribute to its growth, motivating alternative measures. One possibility is to tally the mass of all progenitors that have collapsed by redshift z and that are above some fraction f of the halo's final mass, M0. This quantity, referred to as the ‘collapsed mass history’ (CMH) and denoted Mcoll(z), has a simple interpretation and is easily extracted from simulated or theoretical DM merger trees.
The thick dashed lines in Fig. 2 show the CMH of a 3.2 × 109 h−1 M⊙ halo for several values of the parameter f, and compares them to the median MAH (solid lines, repeated in each panel to aid the comparison). Thin lines show equation (3), adopting δsc = 1.26 for the collapse threshold.4 Note that this expression describes the shape of Mcoll(z) remarkably well for both coco-cold and coco-warm, independent of the value of f adopted. Unlike the main-progenitor MAH, the CMH provides a universal description of the halo assembly process, where the choice of f implicitly defines the halo collapse time: lower values of f imply earlier formation redshifts.
3.3 Mass profiles and concentration estimates
For each halo identified at z0 = 0, 1, 2 and 3 we have constructed spherically averaged density profiles, ρ(r), in 32 equally spaced steps in log r spanning −2.5 ≤ log r/r200 ≤ r200. Within each radial bin we also compute the total enclosed mass, M(r), and mean inner density profiles, 〈ρ〉(r) = M(r)/(4/3)πr3. To ensure that our halo mass profiles are well resolved, we restrict our analysis to those with N200 ≥ 5000 particles within their virial radius, r200.
We construct the c(M, z) relation by fitting median mass profiles after averaging over logarithmic mass bins of width Δlog M = 0.1. This smooths out any features unique to individual systems and dampens the influence of outliers allowing for a robust estimate of the average mass and redshift dependence of halo concentrations.
3.4 Relaxed versus unrelaxed haloes
DM haloes form hierarchically through a combination of smooth accretion, minor mergers and occasional major mergers with systems of comparable mass. These events can drive large but transient departures from quasi-equilibrium states during which the structural properties of DM haloes are rapidly evolving and ill defined. As a result, the majority of studies aimed at calibrating the c(M, z) relation have taken steps to identify and excise haloes believed to be far from equilibrium, thereby defining samples of ‘relaxed’ haloes with smooth mass profiles whose structural features can be meaningfully described with a few parameters. It is important to note, however, that relaxed haloes form a highly biased subsample of the full halo population, and their prevalence depends in non-trivial ways on both halo mass and on redshift (see e.g. Thomas et al. 2001; Macciò et al. 2007; Neto et al. 2007; Power, Knebe & Knollmann 2012; Ludlow et al. 2012; Angel et al. 2016; Klypin et al. 2016).
Identifying relaxed haloes is not without ambiguity, and a number of diagnostics have been proposed in the literature. Some are sensitive to geometric halo properties, such as the centre-of-mass offset parameter, defined as |$d_{\rm off}=|\boldsymbol {r}_{\rm p}-\boldsymbol {r}_{\rm CM}|/r_{200}$| (Thomas et al. 2001; Macciò et al. 2007; Neto et al. 2007), or the mass-fraction in substructure (e.g. Neto et al. 2007; Ludlow et al. 2012); others, such as the virial ratio, η = 2 K/|U| (e.g. Cole & Lacey 1996; Bett et al. 2007; Knebe & Power 2008), or spin parameter, λ (e.g. Klypin et al. 2016), gauge the internal dynamical state of the halo. Some authors reject haloes whose spherically average density profiles are poorly described by their chosen fitting formulae (e.g. Macciò et al. 2007; Macciò, Dutton & van den Bosch 2008; Dutton & Macciò 2014).
Neto et al. (2007) suggested a combination of three criteria that may be used to curtail haloes whose mass profiles are most likely to deviate from smooth spherical averages. These include: (i) the centre-of-mass offset, doff < 0.07, (ii) the substructure mass fraction, fsub = Msub( < r200)/M200 < 0.1, and (iii) the virial ratio, η < 1.35. During a merger each of these quantities fluctuates in predictable ways: doff, for example, traces the centre-of-mass of a merging system about its densest core, providing an estimate of the accuracy with which the halo centre can be defined; fsub monitors the mass contribution from undigested mergers, while η is sensitive to fluctuations in the gravitational potential as orbital energy is dissipated into binding energy. As discussed by Ludlow et al. (2012) and also Poole et al. (2016), merger-driven oscillations in these quantities are out of sync, making it unlikely that a halo will simultaneously fail all three at any point during a merger.
Are these criteria sufficient to ensure removal of all unrelaxed haloes? Arguably not. Because of its resolution dependence, only haloes with well-resolved substructure populations are sensitive to fsub. In simulations with uniform mass resolution (such as those used in this work), the least resolved haloes are also the most abundant; fsub is therefore a useful equilibrium statistic for only the most massive, best-resolved systems. Furthermore, since DM haloes are not truly isolated, the virial ratio should be corrected for external forces and surface pressure terms (see e.g. Poole et al. 2006; Knebe & Power 2008; Klypin et al. 2016, for a discussion) and only then can it be used to meaningfully assess departures from equilibrium.
Because of these uncertainties, we here adopt a simpler approach and use the dynamical age of a system as the primary diagnostic for equilibrium. We assume that any halo whose main progenitor has more than doubled in mass in under a crossing time cannot have had time to relax to an equilibrium configuration. More specifically, we require t50 ≳ 1.25 × tcross as a minimal but necessary condition for equilibrium. Here tcross ≡ 2 × r200/V200 is the characteristic crossing time of a halo and t50 is the lookback half-mass formation time of its main progenitor. This single criterion, however, will only flag haloes undergoing very rapid accretion or equal-mass mergers. For that reason, we additionally impose the familiar criterion doff < 0.1 to cull the remainder of the population. Our primary motivation for choosing these criteria is that they do away with uncertainties surrounding the importance of boundary terms in the virial ratio, and the resolution dependence of fsub. We will see in Appendix A that imposing these criteria on MS haloes results in a c(M, z) relation that decreases monotonically with mass over the redshift range probed by our simulations, removing the ‘upturn’ in the concentration of high-mass haloes reported by Klypin, Trujillo-Gomez & Primack (2011, see Ludlow et al. 2012, for further details).
It is worth noting, however, that different definitions of what constitutes a relaxed halo population lead to conflicting claims regarding the origin of the upturn, and whether or not it is a true property of the underlying structure of equilibrium DM haloes (see e.g. Ludlow et al. 2012; Correa et al. 2015c; Klypin et al. 2016, for discussions). A full assessment will likely require a detailed study of the perils and virtues of a variety of different equilibrium benchmarks, which we defer to future work.
The mass and redshift dependence of the relaxed halo fraction, frel (defined above) is shown in Fig. 3. The left-hand panel shows results for the coco simulations and the right for the MS, which extend to much higher masses. In CDM models halo collapse times decrease monotonically with increasing mass, which is reflected in the decreasing abundance of relaxed haloes amongst massive systems. For WDM models, however, collapse times are not monotonic: there is a maximum formation time for haloes at any given redshift, resulting in a non-monotonic relation between halo mass and the prevalence of relaxed systems.
4 RESULTS
4.1 The c(M, z) relation in CDM and WDM
The mass–concentration-redshift relations for equilibrium haloes in the coco and color-1.5 simulations are shown in Fig. 4. Dots show the best-fitting concentrations obtained for individual haloes, colour-coded to distinguish different runs. Symbols trace the median c(M, z) relations of the same sets of haloes, obtained by fitting the median mass profiles after averaging over logarithmic mass bins of width Δlog M = 0.1 (only bins containing at least 25 haloes are plotted in this figure).
These results confirm and extend previous work on the structure of WDM haloes. Unlike CDM, where concentrations increase monotonically with decreasing mass, in WDM models the c(M, z) relation has a characteristic shape: it first increases with decreasing mass, but reaches a well-defined maximum before decreasing again towards lower mass. Note that in the 3.3 and 1.5 keV models studied here, the mass scale at which the peak concentration is reached is roughly independent of redshift. Note also that differences between WDM and CDM are already evident at mass scales substantially larger than the free-streaming scale. For example, the ‘peak’ in the median concentration of WDM haloes occurs approximately two orders of magnitude above the half-mode mass, suggesting that differences in the very early stages of halo growth leave a lasting imprint on the final halo. We highlight this point using coloured arrows, which correspond to one hundred times the half-mode mass of each WDM run. This provides an important clue for models that aim to fully describe the c(M, z) relations from the power spectrum alone. The solid curves in Fig. 4 show the predictions of such a model, which we describe in more detail in Section 4.3.
4.2 Mass profiles and assembly histories
As discussed by L13, the MAHs and mass profiles of CDM haloes have, on average, the same NFW shape. This implies that halo concentrations can be obtained by simply rescaling their MAH by a fixed amount: the characteristic density of M(〈ρ〉) is simply proportional to that of M(ρcrit(z)). One consequence of this result is that two haloes of similar MAH must have the same characteristic density and vice versa, independent of their mass or identification redshift. This greatly simplifies the task of predicting concentrations from assembly histories when applied to CDM (see also Correa et al. 2015c). In WDM models, however, the suppression of gravitational collapse below the free-streaming scale breaks the scale invariance of the assembly process: it imprints a preferred scale on the MAHs, readily seen in the leftmost panel of Fig. 1. This implies that the mass profiles of WDM haloes cannot be obtained by simply rescaling the MAHs of their main progenitors.
We illustrate this in the upper panels of Fig. 5, where we show the median main-progenitor MAHs, M200(z) (solid lines), and the enclosed density profiles, M(〈ρ〉) (dashed lines), of coco-warm haloes for two different halo masses and at two different redshifts. The halo masses (log M200/[1010 h−1 M⊙] = 0.9 and log M200/[1010 h−1 M⊙] = −0.75) were selected so that their median concentration is roughly the same (c ≈ 9.5 at z0 = 0, and c ≈ 6.4 at z0 = 1), but fall on opposite sides on the ‘peak’ in the c(M) relations. To aid the comparison, masses have been normalized to the current mass, M0 = M(z0), critical densities to ρcrit(z0), and enclosed densities to 200 ρcrit(z0).
The dashed curves indicate that these haloes not only have the same c, but also similar mass profiles across the entire resolved radial range. The outsized symbols highlight the enclosed mass and mean density within r = 3 × r−2, r−2 and r−2/2, which are roughly equivalent for both masses. For comparison, the solid red line shows an NFW profile with the same concentration.
Despite the similarity of the halo mass profiles, it is clear from Fig. 5 that the shapes of the MAHs of the two haloes (solid lines) are substantially different. The dot–dashed curve shows the MAH obtained by rescaling the NFW mass profile, as described in L13 for CDM. This model describes quite well the MAH of massive WDM haloes, but fails dramatically at low mass, where the MAH shape differs substantially from NFW. The MAH of such haloes cannot be used then to infer the concentration of their mass profile in the same way as for CDM haloes.
Alternative descriptions of halo growth that preserve scale invariance may improve matters. One possibility, mentioned in Section 3.2, is to use the mass, Mcoll(f, z), in all progenitors (above a certain threshold f × M0) rather than just that of the main progenitor. The lower panels of Fig. 5 show, for the same two halo masses, the growth of the total mass in progenitors more massive than 2 per cent of the halo's final mass. The curves are now virtually indistinguishable, suggesting that the collapsed mass in progenitors other that the main one plays an important role in establishing a halo mass profile. Note also that the shape of Mcoll(z) is accurately described by equation (3), shown in the lower panels using a dot–dashed (red) line after rescaling to match each of the halo formation times. This suggests that it may be possible to use the CMH, Mcoll(f, z), to predict halo concentrations.
Fig. 6 shows that this is indeed the case, for both CDM and WDM haloes, regardless of mass or identification redshift. Here we compare the median Mcoll(z) (dot–dashed lines), constructed using f = 0.02 for haloes of mass M0 = 109, 1010 and 1011 h−1 M⊙ after rescaling each to the characteristic values of mass, M−2 = M( < r−2), and density, ρ−2.
Different panels show results for different identification redshifts, with blue and orange curves distinguishing haloes in the coco-cold and coco-warm runs, respectively. As anticipated, each curve has a similar shape, independent of M0, z0 or the shape of the DM power spectrum. Indeed, as alluded to above, the similarity of these curves is actually expected from EPS theory. The dashed grey line in each panel of Fig. 6, for example, shows equation (3), which describes the shapes of these curves remarkably well.
Fig. 6 also shows the median enclosed density profiles, M(〈ρ〉), for each set of haloes, again after normalizing to the characteristic values of M−2 and 〈ρ−2〉. The solid grey line shows an Einasto profile with α = 0.18, which provides an accurate approximation to the median mass profiles of these haloes over the mass and redshift range probed by our simulations.5
Because both Mcoll(ρcrit(z)) and M(〈ρ〉) have self-similar (albeit distinct) shapes, a single scaling relation between characteristic density and formation time is sufficient to anchor the two and to construct an analytic model for c(M, z), provided that, for some fixed value of f, the characteristic density of the CMH is simply proportional to that of the mass profile.
We show this in Fig. 7, where we plot the mean enclosed density within the halo scale radius, |$\langle \rho _{-2}\rangle =M_{-2}/(4/3)\pi r_{-2}^3$|, versus ρcrit(z−2), the critical density at the redshift z−2 when the enclosed mass was first assembled into progenitors more massive than 2 per cent of M0 (i.e. f = 0.02). Note that the linear scaling, 〈ρ−2〉 = C × ρcrit(z−2), between these two densities is independent of both mass and identification redshift. More importantly, however, the zero-point of this relation is independent of the DM particle model: both coco-cold and coco-warm have C ≈ 400. Note also that similar scalings can be found for characteristic densities measured within different fractions of r−2. This can be seen in the upper and lower panels of Fig. 7 which show the results of repeating this calculation within r−2/2 and 3 × r−2, respectively.
4.3 An analytic model for the c(M, z) relation
The results of the previous section suggest that the c(M, z) relation can be inferred from CMHs, Mcoll(z), provided those can be obtained from either simulations or theoretical models. Encouragingly, the most recent generation of algorithms accurately reproduce not only the evolution of the main progenitor halo (e.g. van den Bosch 2002; Jiang & van den Bosch 2014; Correa et al. 2015a,b), but the entire hierarchy of progenitors for both CDM (e.g. Parkinson, Cole & Helly 2008) and WDM (e.g. Benson et al. 2013) fluctuation power spectra.
In order to describe the average relation between concentration, mass and redshift, we use analytic arguments based on EPS theory to construct a simple model. Following the procedure laid out by NFW, the model assumes that a halo's characteristic density reflects the critical density of the Universe at a suitably defined collapse redshift. NFW adopted δc as the characteristic density (equation 1), and defined the collapse redshift as the time when half the virial mass of the halo was first contained in progenitors more massive than some fraction f of the final virial mass.
A simple modification to this procedure yields much better results which are applicable to both CDM and WDM initial power spectra. These modifications identify the characteristic density with the mean inner density within the scale radius, 〈ρ−2〉, rather than δc, and the characteristic halo mass with M−2 rather than the virial mass.
The solid curves in Fig. 4 show the resulting c(M, z) relations for f = 0.02 and using a normalization constant7 of C = 650, assuming that ρ(r) follows an Einasto8 profile with α = 0.18. Note that this model accurately reproduces the median trends for CDM, and both WDM cosmologies considered here without any additional tuning of the parameters.
Having calibrated the model parameters using these runs, we can now make predictions for c(M, z) that can be tested against simulations of different WDM and cosmological models. One such test is shown in Fig. 8, where we plot the c(M) relations for equilibrium haloes in four separate flat ΛCDM cosmological models. Each has H0 = 73 km s−1Mpc−1 and ns = 1, but varies the matter density parameter, ΩM, and rms fluctuation amplitude, σ8, as indicated in the legends. As in Fig. 4, dots show the best-fitting concentrations for individual haloes; filled symbols show the median relations. The solid lines show the predictions of equations (6) and (7) for C = 650 and f = 0.02, assuming an Einasto mass profile with α = 0.18. Note that these are the same parameters used to fit the c(M, z) relations in Fig. 4, which were obtained from simulations of a WMAP-7 cosmological model.
A further test of the model involves its extrapolation to much lower halo masses than can be resolved in coco-cold and coco-warm. We show this in Fig. 9 where we compare the Vmax–rmax relation (equivalent to the mass–concentration relation) for field haloes in the high-resolution region of the Aquarius simulations of Lovell et al. (2014). These runs follow the evolution of Milky Way-mass DM haloes and their immediate surroundings in CDM and a series of WDM models with thermal particle masses of ≈1.5, 1.6, 2.0 and 2.3 keV. We plot here isolated haloes, defined as main haloes that are farther than two virial radii from the largest halo in the simulation. This selection ensures that systems previously accreted and expelled from the most massive halo in the simulation are excluded from the analysis (see Ludlow et al. 2009). Grey dots show the Vmax–rmax relation obtained from the CDM run; coloured symbols show the results for the various WDM models, as indicated in the legend.
The predictions of the analytic model described above are shown in Fig. 9 for all DM models considered. As in previous plots, we compute each curve for C = 650 and f = 0.02; line and symbol colours were chosen to match the corresponding simulations. These curves accurately describe the systematic shift in concentration brought about by changes to the properties of the DM particle. The good agreement between model and simulation implies that our model provides a useful theoretical tool for studies aiming to constrain the nature of DM based on its imprint on the internal structure of low-mass haloes.
4.4 Concentrations of individual haloes
Although equations (6) and (7) provide a robust account of the median c(M, z) relation in each of our simulations, our methodology can also be applied to individual haloes. To do so, we select equilibrium systems with N200 > 104 particles from the coco-cold and coco-warm simulations and use their merger trees to construct Mcoll(z). Assuming that each halo's mass profile is reasonably well described by an (α = 0.18) Einasto profile, its concentration can be readily obtained from equation (6), with C = 400 in this case.
Fig. 10 plots the predicted versus the measured concentrations. The left-hand panel shows results for coco-cold, the right-hand panel for coco-warm. The coloured contours enclose 75 and 50 per cent of the data points, and are shown for z0 = 0, 1 and 2. This figure makes clear that our procedure faithfully reproduces the concentrations of relaxed haloes. The rms scatter about the one-to-one line, for example, is typically less than ∼0.09, indicating that the error on the predicted value of c is less than ∼23 per cent.
The fraction of scatter in c(M) that is due to different halo collapse times can be estimated by comparing the variance in the measured and predicted values of concentration for haloes in a given mass bin. For CDM haloes at z0 = 0, we find that variation in z−2 accounts for ∼40 per cent of the scatter for M200 ∼ 109h−1 M⊙, and ∼50 per cent at ∼109h−1 M⊙. For WDM haloes of the same mass and redshift, the reductions are ∼80 and ∼50 per cent, respectively. Thus, the scatter in collapse time does not fully account for the scatter in concentration. Future studies should assess the role of environment or initial conditions in establishing the mass profiles of DM haloes.
4.5 Comparison to previous work
Studies aimed at providing theoretical predictions for the c(M, z) relation have traditionally followed one of two routes. One class of models aims to connect the structural properties of haloes to some aspect of their assembly history (Bullock et al. 2001; Eke et al. 2001; Wechsler et al. 2002; Zhao et al. 2003; Lu et al. 2006; Macciò et al. 2008; L14; Correa et al. 2015c) most commonly by relating their characteristic densities to the mean or background density at some appropriately defined formation time. Other methods devise fitting formulae of varying complexity that are then calibrated to the results of numerical simulations over some range of halo mass and redshift (Prada et al. 2012; Dutton & Macciò 2014; Diemer & Kravtsov 2015; Klypin et al. 2016).
Each of these methods has their own virtues and weaknesses. Parametrized fits to simulation results, for example, generally yield simple and compact formulae that accurately describe the scaling relations between halo structural properties, but must be treated with caution when extrapolated outside the range of halo mass, redshift or cosmological parameters for which they were determined. Physically motivated models, however, attempt to link halo structure to some aspect of their assembly. As a result they are often more cumbersome, but arguably provide more reliable extrapolations of these relations.
In Fig. 11, we compare the predictions of our model (thin solid lines repeated in each panel) for the Planck cosmology with several empirical fits proposed in the literature (dashed lines). Different colours correspond to different redshifts – z0 = 0, 1, 2 and 3 – and all curves have been extrapolated down to ∼10−8 h−1 M⊙ to emphasize their differences.
Each model was calibrated, at z = 0, for masses 1010 ≲ M200/[h−1 M⊙] ≲ 1015 and over that range provide consistent predictions for the c(M) relation. Note, however, that the models show a larger variation at higher redshift, even over the mass scales at which they were calibrated (highlighted using thick lines). This reflects the fact that each group of authors used different equilibrium conditions (or none at all in the case of Diemer & Kravtsov 2015 and Prada et al. 2012) to define their halo samples, which can introduce subtle biases in the recovered c(M, z) relation, particularly at high redshift (Ludlow et al. 2012; Angel et al. 2016). In addition, different authors employ different techniques for estimating halo concentrations, which may also account for some of the differences (see e.g. Dutton & Macciò 2014, for a discussion).
More importantly, however, these models fail to account for the effects of CDM free-streaming on haloes of the smallest mass. This is shown in the lower-right panel, where the dashed lines show the effect of a 100 GeV neutralino CDM particle relative to pure CDM (solid lines) on the predicted c(M, z) relation (we model the neutralino transfer function according to Green, Hofmann & Schwarz 2004). Note that the concentrations of neutralino haloes peak at approximately 60 Earth masses, shown here using a downward pointing arrow. Regardless of the substructure mass function below this limit, the reduced concentration of haloes will substantially reduce the ‘clumpiness’ of matter on scales comparable to Earth's mass and smaller. This may have important implications for direct or indirect DM detection experiments.
5 SUMMARY
We have analysed the mass profiles and formation histories of a large sample of equilibrium DM haloes extracted from an ensemble of cosmological N-body simulations. These include the Copernicus Complexio simulations, the Millennium and Aquarius simulations and an additional set of runs with differing density fluctuation amplitudes and DM content. Our simulation suite includes both CDM and WDM runs, allowing us to assess the dependence of halo concentrations on cosmological parameters and on the shape of the linear matter power spectrum. Below we summarize our main results.
In agreement with previous work, we find that the spherically averaged density profiles of CDM and WDM haloes are approximately universal and are accurately described by the NFW or Einasto profile, regardless of the shape of the density fluctuation power spectrum. The main difference between CDM and truncated power spectra, such as WDM, is the mass–concentration–redshift relation of equilibrium haloes.
The concentrations of CDM haloes, at fixed z, decrease monotonically with increasing halo mass. At fixed mass, concentrations decrease monotonically with increasing redshift. We find no evidence for departures from this trend amongst our relaxed sample, providing further evidence that the recently reported ‘upturn’ in the concentrations of rare, massive haloes results from the inclusion of unrelaxed systems in the sample (see e.g. Ludlow et al. 2012, for details).
The WDM c(M) relation, on the other hand, is non-monotonic. At given z, concentrations peak at a multiple of the half-mode mass scale imposed by the power spectrum truncation. Concentrations decline above and below that peak mass scale (which seems to be independent of redshift) over the full mass range probed by the simulations.
The main-progenitor MAHs of CDM haloes are scale-free, and resemble the NFW shape when cast in terms of mass and critical density, rather than mass and time. The concentration of a CDM halo can therefore be inferred from a simple scaling of its MAH. A simple model based on this feature accounts successfully for the main properties of the CDM mass–concentration–redshift relations.
The ‘free-streaming’ truncation of the WDM linear power spectrum results in substantial delay and suppression of structure formation below a certain mass scale. This breaks the similarity of main-progenitor MAHs, whose shape now depends critically on halo mass. Low-mass WDM haloes near the truncation scale tend to form quickly and monolithically. Massive WDM haloes, on the other hand, form just like their CDM counterparts. The strong dependence on halo mass of the MAH shape precludes the application of the same model developed to infer CDM halo concentrations from their MAHs.
A simple extension of the model, however, provides an excellent description of the mass–concentration–redshift relation for both CDM and WDM haloes. The extension relies on using the full ‘CMH’ of a halo (i.e. Mcoll(z), the total mass in collapsed structures at given time) rather than just its main progenitor. Collapsed mass fractions are approximately scale-free, which implies that a simple model may be devised in order to infer halo concentrations directly from Mcoll(z). This model reproduces well the mass–concentration–redshift relation in all our simulations, CDM and WDM alike, over the whole range of halo masses and redshifts probed. Applied to CDM, this model may be used to extrapolate c(M, z) to very low masses, where the relation departs strikingly from a pure power law and differs also from the predictions of earlier models.
We provide upon request a simple code that computes the predicted c(M, z) relations for CDM and WDM power spectra with arbitrary free-streaming truncation masses. This should be useful for extrapolating the results we show here to halo masses and redshifts not well resolved by our simulations. Such extensions may be important, especially when evaluating predictions of these models for indirect and direct DM detection strategies.
We thank Mark Lovell for useful conversations and our referee for a constructive report that improved the paper. ADL is supported by a COFUND Junior Research Fellowship; SB by the STFC through grant ST/K501979/1. LW acknowledges support from the NSFC grants program (no. 11573031, no. 11133003), REA from AYA2015-66211-C2-2, and WAH from a Science and Technology Facilities Council grant ST/K00090/1 and the Polish National Science Centre under contract #UMO-2012/07/D/ST9/02785. This work was supported by the Science and Technology Facilities Council (grant number ST/F001166/1) and European Research Council (grant number GA 267291, ‘Cosmiway’). We also acknowledge financial support from the STFC consolidated grant ST/L00075X/1. This work used the COSMA Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.
We define the virial mass of a halo as that enclosed by a sphere (centred on the potential minimum) of mean density equal to 200 times the critical density, ρcrit = 3H2/8πG, where H(z) is the Hubble constant; the virial radius is therefore implicitly defined by |$200\rho _{\rm crit}={M}_{200}/(4/3)\pi r_{200}^3$|. Note that all particles in the simulation are used in calculating M200 and not only those deemed gravitationally bound to a particular halo.
Various aspects of the post-processed simulation data may be made available by the first author upon request.
The half-mode mass, Mhm, indicates the scale at which the WDM transfer function is reduced by half relative to a CDM model with the same cosmological parameters. In our coco-warm model this corresponds to Mhm = 2.46 × 108 h−1 M⊙, and in color-1.5 to Mhm = 3.4 × 109 h−1 M⊙. Half-mode masses for our remaining WDM runs are provided in Table 1.
Throughout the paper, we calculate σ(M) from the linear power spectrum using a real-space spherical top-hat filter. Although other possibilities exist, our tests show that the best results are obtained with this choice.
Note that this zero-point is different to C = 400 which is obtained from the simulations (see the middle panel of Fig. 7). This is due to inaccuracies of the spherical collapse model, which is implicitly assumed in the calculation of equation (7). Merger trees generated using Monte Carlo methods tailored to reproduce simulations results (e.g. Parkinson et al. 2008) should adopt C = 400 and f = 0.02.
We show in Appendix B that the predictions of the model depend only weakly on the assumed shape of the DM density profile, provided it resembles those of simulated haloes. For simplicity, we have adopted a single Einasto profile with shape parameter α = 0.18 throughout the paper.
REFERENCES
APPENDIX A: Comparison to other physical c(M, z) models
In Fig. A1, we compare the predictions of the model described in Section 4.3 to three common analytic recipes to compute c(M, z) (Bullock et al. 2001; Gao et al. 2008; Correa et al. 2015c) Each model adopts a WMAP-1-normalized cosmology, and its predictions are compared to the median c(M, z) relations obtained for relaxed haloes in the MS runs.
All three analytic prescriptions considered here agree reasonably well with our results at z = 0. The Bullock et al. (2001) model, including modifications suggested by Macciò et al. (2008, left-hand panel), deviates by less than ∼5 per cent at any mass scale ≲1015 h−1 M⊙. This model predicts a simple redshift dependence at fixed mass, c∝(1 + z)−1, which, at low masses (≲1010 h−1 M⊙), agrees well with our findings. At higher mass, however, it predicts a sharp decline that is inconsistent with the results of the MS. For WDM haloes this model predicts halo concentrations that increase monotonically with decreasing mass, approaching a constant for masses below the free-streaming scale (this is because Bullock et al. define formation times using D(zc)σ(F M0) = 1.686 which, for WDM, predicts a constant zc for M ≪ Mfs). Clearly this is not supported by the data in Fig. 4.
The modifications of the NFW model proposed by Gao et al. (2008) results in concentration estimates that are similar, but not identical to the predictions of the model we propose. The original NFW model has three free parameters: two physical parameters, F and f, and a scaling factor C that relates the characteristic density of a halo to the background density at its formation time. The formation time corresponds to the time which a fraction F of the halo's final mass, M0, was first assembled into progenitors each at least as massive as f × M0.
NFW chose F = 0.5, f = 0.01 and C = 3000, whereas Gao et al. found improved fits to the MS by choosing F = 0.1, f = 0.01 and C = 600. Our proposal instead relates F to the characteristic mass of the halo, i.e. F ≡ M−2/M0; it adopts f = 0.02 and C = 650. We compare the NFW predictions for these parameter choices in the middle panel of Fig. A1.
Note that defining F ≡ M−2/M0 reduces the number of model parameters while at the same time improves the agreement with the simulation results. Note also, that the Gao et al. model is expected to fare increasingly poorly at mass scales for which it predicts concentrations for which M−2/M0 ≪ or ≫0.1. As a result, we expect extrapolations of the Gao et al. model to very low or very high mass scales to become increasingly different. We provide a more detailed discussion of the importance of the parameters F and f in Appendix B.
The analytic predictions of Correa et al. (2015c), shown in the right-hand panel of Fig. A1, also agree well with our results. Like L14, this model relates the characteristic density of a halo to the critical density at the redshift when its main progenitor had first assembled the characteristic mass, M−2. The agreement is therefore not surprising, as their model is based explicitly on the results L14, which we deliberately reproduce using our new methodology. Indeed, the slight differences between the predictions of Correa et al. (2015c) and our own can be attributed to changes brought about by replacing Einasto fits with the NFW fits upon which their model was calibrated. Note, however, that the model of L14 (and Correa et al.) fails when applied to WDM power spectra because free-streaming modifies the proportionality constant between the characteristic density of a halo and the critical density at their particular choice of the collapse time.
APPENDIX B: PARAMETER DEPENDENCE OF THE CHARACTERISTIC DENSITY–FORMATION TIME RELATION
The model we propose here, like the NFW and Bullock et al. models, is afflicted by the introduction of two physical parameters, F and f, as well as a normalization constant, C, whose interpretation is unclear. The scaling relations in Fig. 7 indicate that the value of C will depend sensitively on the precise definition of ‘characteristic density’. If defined to be 〈ρ−2〉, the mean density within r−2, then C = 400. If, however, the characteristic density is defined as that within ∼3 × r−2 then C ≈ 200 describes our numerical results quite well. Note that this is precisely what is expected from the spherical top-hat collapse model. The simplest interpretation of the parameter C therefore is that it represents the contraction of the inner halo beyond what is expected from simple spherical collapse. Lacking a theoretical model to provide deeper insight, its numerical value must be calibrated using simulations, and may exhibit subtle dependences on the precise definition of halo mass used.
The success (of lack thereof) of previous attempts by NFW and Gao et al. to model the c(M, z) relation using methods similar to ours can be traced to each authors different definition of ‘formation time’. Both NFW and Gao et al. assumed that a halo had ‘formed’ when a fraction F of its final mass had first assembled into clumps more massive than a smaller fraction f ≪ F of the same mass. Each of these parameters affects the precise value of zf. NFW assumed F = 0.5, Gao et al. F = 0.1, and both found f = 0.01 as a suitable progenitor mass threshold.
In Fig. B1, we plot the mean enclosed density within the halo scale radius, |$\langle \rho _{-2}\rangle =M_{-2}/(4/3)\pi r_{-2}^3$|, versus ρcrit(zF), the critical density at the redshift zF when the total collapsed mass first exceeded F × M0. Different panels show results for different values of F, and all assume a small progenitor threshold, f = 0.02. The points corresponds to averages over single mass bins, and are shown for all four identification redshifts after normalizing densities by the critical value at z0, ρcrit(z0). Note that only mass bins corresponding to haloes with N200 ≥ 104 are plotted, which ensures that all progenitors are resolved with at least 100 particles.
Although all values of F (which span a factor of 50) result in a tight linear relation between the two densities, only the particular choice F = M−2/M0 provides a direct proportionality between 〈ρ−2〉 and ρcrit(z−2) (the thin grey line in each panel, for example, shows the direct scaling 〈ρ−2〉 = 400 × ρcrit(z−2), which provides a reasonable description of the data only for F = M−2/M0). All other values result in steeper slopes that gradually shallow as F is decreased but, even for F = 0.01, do not reach the natural linear scaling.
Formation times also depend on the progenitor mass threshold, f. Previous studies have hinted at puzzlingly small values: f = 0.01 in the case of NFW and Gao et al; f = 0.02 in our case. In Fig. B2 we show how the mean characteristic density, 〈ρ−2〉, of haloes varies as a function of the redshift at which their characteristic mass, M−2, had first assembled into clumps each larger than f × M0. Note that high values of f (e.g. f = 0.5, shown in the lower-right panel) result in steep power-law slopes. Decreasing f shifts all formation times to higher redshift, but the magnitude of the shift depends on 〈ρ−2〉. The net result, provided f becomes sufficiently small, provides a natural correspondence between characteristic density and the background density at the halo formation time. Note also that, at least in the CDM case, the precise values of f seems unimportant, provided it is ‘small enough’ (e.g. f = 0.005, 0.01 and 0.02 all result in scaling relations whose slopes do not deviate noticeably from one). As a rule of thumb, we suspect that values of f ≪ M−2/M0, of the order of a few per cent, will yield robust results.
Overall, we find that for F = M−2/M0 and f ≈ 0.02, the linear scaling, 〈ρ−2〉 = C × ρcrit(z−2), between these two densities is independent of both mass and identification redshift. More importantly, however, the zero-point of this relation is independent of the DM particle model: both coco-cold and coco-warm have C ≈ 400. We therefore advocate the use of these parameters in future modelling, but acknowledge that better data, alternative halo definitions, or a different variety of DM models may result in modifications.
Mapping between characteristic densities and concentrations requires an assumption regarding the halo mass distribution. Throughout the paper we have adopted an Einasto profile with α = 0.18. In Fig. B3, we plot the predicted concentration–mass relations at several redshifts and for a few other density profiles. The solid coloured lines correspond to our fiducial α = 0.18 Einasto profile, dashed lines to an NFW profile and the shaded regions highlight the range of c that is expected for Einasto models with a redshift-dependent scatter in α consistent with the findings of Dutton & Macciò (2014). Note that at high mass, above ∼1010 h−1 M⊙, the predicted concentrations depend weakly on the assumed mass profile. Towards lower mass, however, when concentrations become large, a modest halo-to-halo scatter in α can lead to a considerable scatter in the predicted value of c.
APPENDIX C: A fitting formula for the concentration–mass–redshift relation in the Planck cosmology
In Fig. C1, we compare the c(M, z) and c(ν) relations predicted by our model (coloured points) with the above fitting formula (solid lines). In agreement with previous studies (e.g. Dutton & Macciò 2014; Diemer & Kravtsov 2015; Hellwing et al. 2016), the c(ν) relation depends slightly but systematically on redshift. Note also that the residuals (shown in the lower panels) are small, typically less than ∼3 per cent at all mass scales relevant for the CDM cosmology, and show no systematic dependence on halo mass, M, or redshift.
In the Planck cosmology, the concentration of any halo of mass M at redshift z can therefore be estimated as follows.
Use the resulting values of these parameters to calculate c from equation (C1).