The redshift dependence of the structure of massive Λ cold dark matter haloes

Abstract

We use two very large cosmological simulations to study how the density profiles of relaxed Λ cold dark matter dark haloes depend on redshift and on halo mass. We confirm that these profiles deviate slightly but systematically from the NFW form and are better approximated by the empirical formula, d log ρ/d log r∝r^α, first used by Einasto to fit star counts in the Milky Way. The best-fitting value of the additional shape parameter, α, increases gradually with mass, from α∼ 0.16 for present-day galaxy haloes to α∼ 0.3 for the rarest and most massive clusters. Halo concentrations depend only weakly on mass at z= 0, and this dependence weakens further at earlier times. At z∼ 3 the average concentration of relaxed haloes does not vary appreciably over the mass range accessible to our simulations (M≳ 3 × 10¹¹h⁻¹M_⊙). Furthermore, in our biggest simulation, the average concentration of the most massive, relaxed haloes is constant at 〈c₂₀₀〉∼ 3.5–4 for 0 ≤z≤ 3. These results agree well with those of Zhao et al. and support the idea that halo densities reflect the density of the universe at the time they formed, as proposed by Navarro, Frenk & White. With their original parameters, the NFW prescription overpredicts halo concentrations at high redshift. This shortcoming can be reduced by modifying the definition of halo formation time, although the evolution of the concentrations of Milky Way mass haloes is still not reproduced well. In contrast, the much-used revisions of the NFW prescription by Bullock et al. and Eke, Navarro & Steinmetz predict a steeper drop in concentration at the highest masses and stronger evolution with redshift than are compatible with our numerical data. Modifying the parameters of these models can reduce the discrepancy at high masses, but the overly rapid redshift evolution remains. These results have important implications for currently planned surveys of distant clusters.

methods: N-body simulations, methods: numerical, galaxies: haloes, galaxies: structure, dark matter

1 INTRODUCTION

Over the past decade, cosmological N-body simulations have shown consistently that equilibrium dark matter haloes have spherically averaged mass density profiles which are approximately ‘universal’ in form; that is, their shape is independent of mass, of the values of the cosmological parameters, and of the linear power spectrum from which non-linear structures have grown. As a result, it is useful to parametrize halo profiles by simple empirical formulae, such as that proposed by Navarro, Frenk & White (1995, 1996, 1997, hereafter NFW):

1

where ρ_crit= 3H²/8πG is the critical density for closure,1δ_c is a characteristic density contrast, and r_s is a scale radius. Note that this formula contains two scale parameters but no adjustable shape parameter.

As discussed in some detail by NFW and confirmed by subsequent numerical work, the two parameters of the NFW profile do not take arbitrary values, but are instead correlated in a way that reflects the mass dependence of halo assembly times (e.g. Kravtsov, Klypin & Khokhlov 1997; Avila-Reese et al. 1999; Ghigna et al. 2000; Jing 2000; Bullock et al. 2001; Eke, Navarro & Steinmetz 2001; Klypin et al. 2001). The basic idea behind this interpretation is that the characteristic density of a halo tracks the mean density of the universe at the time of its formation. Thus, the later a halo is assembled, the lower its characteristic density, δ_c, or, equivalently, its ‘concentration’ (see Section 2.2 for a definition).

Although the general validity of these trends is well established, a definitive account of the redshift and mass dependence of halo concentration is still lacking, even for the current concordance cosmology. This is especially true at high masses, where enormous simulation volumes are required in order to collect statistically significant samples of these rare systems. Simulating large cosmological volumes with good mass resolution is a major computational challenge, and until recently our understanding of the mass profile of massive haloes has been rather limited, derived largely from small numbers of individual realizations or from extrapolation of models calibrated on different mass scales (NFW; Moore et al. 1998; Ghigna et al. 2000; Klypin et al. 2001; Diemand, Moore & Stadel 2004; Navarro et al. 2004; Tasitsiomi et al. 2004; Reed et al. 2005).

Individual halo simulations may result in biased concentration estimates, depending on the specific selection criteria used to set them up. In addition, they are unlikely to capture the full scatter resulting from the rich variety of possible halo formation histories. Extrapolation based on poorly tested models can also produce substantial errors, as recently demonstrated by Neto et al. (2007). These authors analysed the mass–concentration relation for haloes identified at z= 0 in the Millennium Simulation (MS) of Springel et al. (2005) and confirmed the earlier conclusion of Zhao et al. (2003b) that the models of Bullock et al. (2001, hereafter B01) and Eke et al. (2001, hereafter ENS) (which were calibrated to match galaxy-sized haloes) severely underestimate the average concentration of massive clusters, by up to a factor of ∼3.

Estimates of concentrations can also be biased by the inclusion of unrelaxed haloes. These often have irregular density profiles caused by major substructures. Smooth density profiles are often poor fits to such haloes, and the resulting concentration estimates are ill-defined because they depend on the radial range of the fit and choice of weighting. They can also lead to spurious correlations (see e.g. fig. 9 of Neto et al. 2007). Consequently, in this paper we follow Neto et al. and select only relaxed haloes for analysis. This is not without its own problems. Such selection biases against recently formed haloes, which may preferentially have lower concentrations. However, we believe that this is preferable to polluting the sample with meaningless concentration estimates of the kind that arise when smooth spherical models are force-fitted to lumpy, multimodal mass distributions. Hayashi & White (2008) stacked all haloes, regardless of dynamical state, in the MS and studied the resulting mean profiles as a function of halo mass. The relatively small differences between their results and those found below shows that the inclusion of unrelaxed haloes has rather little effect on the mean.

A further pre-occupation concerns indications that halo profiles deviate slightly but systematically from the NFW model (Moore et al. 1998; Fukushige & Makino 2001; Jing & Suto 2002; Navarro et al. 2004; Merritt et al. 2006; Prada et al. 2006), raising the possibility that estimating concentrations by force-fitting simple formulae to numerical data may result in subtle biases that could mask the real trends. This is especially important because of hints that such deviations depend systematically on halo mass (Navarro et al. 2004; Merritt et al. 2005). Evaluating and correcting for such deviations is important in order to establish conclusively the mass and redshift dependence of halo concentration.

These uncertainties are unfortunate since observations, especially at high redshift, often focus on exceptional systems. For example, massive galaxy clusters are readily identified in large-scale surveys of the distant universe, and understanding their internal structure will be critical for the correct interpretation of cluster surveys intending to constrain the nature of dark energy. These will make precise measurements of the evolution of cluster abundance in samples detected by gravitational lensing, by their optical or X-ray emission, or through the Sunyaev–Zel'dovich (SZ) effect (see e.g. Carlstrom, Holder & Reese 2002; Hu 2003; Majumdar & Mohr 2003; Holder 2006, and references therein).

There is at present no ab initio theory that can reliably predict the internal structure of cold dark matter (CDM) haloes. The models of Zhao et al. (2003a) and Wechsler et al. (2002) (see also Lu et al. 2006) are interesting, but as shown in Neto et al. (2007) they account for only a small fraction in the measured dispersion in concentration at a fixed mass. There is currently no substitute for direct numerical simulation when detailed predictions are needed for comparison with observation.

We address these issues here by combining results from the MS with results from an additional simulation which followed a substantially smaller volume but with better mass resolution. This allows us to extend the range of halo masses for which we can measure concentrations and to assess how these measures are affected by numerical resolution. Our analysis procedure follows closely that of Neto et al. (2007). In particular, we concentrate in this paper on the properties of haloes which are relaxed according to the criteria defined by these authors; mean density profiles for all MS haloes of given mass, regardless of dynamical state, are presented by Hayashi & White (2008). We begin in Section 2 by describing briefly the numerical simulations and the halo catalogue on which this study is based. In Section 3, we present our main results for the dependence of profile shape and concentration on halo mass and redshift. We conclude with a brief discussion and summary in Section 4.

2 THE NUMERICAL SIMULATIONS

The analysis presented here is based primarily on haloes identified in the MS of Springel et al. (2005). The halo identification and cataloguing procedure follows closely that described in detail by Neto et al. (2007). For completeness, we here recapitulate the main aspects of the procedure, referring the interested reader to the earlier papers for details.

2.1 Simulations

The MS is a large N-body simulation of the concordance ΛCDM cosmogony. It follows N= 2160³ particles in a periodic box of side L_box= 500 h⁻¹Mpc. The chosen cosmological parameters were Ω_m=Ω_dm+Ω_b= 0.25, Ω_b= 0.045, h= 0.73, Ω_Λ= 0.75, n= 1, and σ₈= 0.9. Here, Ω denotes the present-day contribution of each component to the matter-energy density of the Universe, expressed in units of the critical density for closure, ρ_crit; n is the spectral index of the primordial density fluctuations, and σ₈ is the rms linear mass fluctuation in a sphere of radius 8 h⁻¹ Mpc at z= 0. The particle mass in the MS is 8.6 × 10⁸h⁻¹M_⊙. Particle interactions are softened on scales smaller than the (Plummer-equivalent) softening length, ε= 5 h⁻¹kpc.

We also use a second simulation of a smaller volume (100³h⁻³ Mpc³) within the same cosmological model. This simulation followed N= 900³ particles of mass 9.5 × 10⁷h⁻¹M_⊙ and softened interactions on scales smaller than ε= 2.4 h⁻¹kpc. It thus has approximately nine times better mass resolution than the MS. Hereafter, we refer to it as the hMS.

2.2 Halo catalogues

Our halo cataloguing procedure starts from a standard b= 0.2 friends-of-friends list of particle groups (Davis et al. 1985) and refines it with the help of subfind, the subhalo finder algorithm described by Springel et al. (2001). Each halo is centred at the location of the minimum of the gravitational potential of its main subfind subhalo, and this centre is used to compute the virial radius and mass of each halo.

We define the virial radius, r_Δ, as that of a sphere of mean density Δ×ρ_crit. Note that this defines implicitly the ‘virial mass’ of the halo, M_Δ, as that enclosed within r_Δ. The choice of Δ varies in the literature. The most popular choices are: (i) a fixed value, as in NFW's original work, where Δ= 200 was adopted; or (ii) a value motivated by the spherical top-hat collapse model, where Δ∼ 178 Ω^0.45_m for a flat universe (Eke, Cole & Frenk 1996). The latter choice gives Δ= 95.4 at z= 0 for the cosmological parameters adopted for the MS. We keep track of both definitions in our halo catalogue, and will specify our choice by a subscript. Thus, M₂₀₀ and r₂₀₀ are the mass and radius of a halo adopting Δ= 200, whereas M_vir and r_vir are the values corresponding to Δ= 95.4 at z= 0. Quantities specified without subscript assume Δ= 200, so that, for example, M=M₂₀₀. The halo concentration is defined as the ratio between the virial radius and the scale radius: c=c₂₀₀=r₂₀₀/r_s. In this case, the characteristic density is related to the concentration by

2

Since haloes are dynamically evolving structures, we use a combination of diagnostics in order to flag non-equilibrium systems. Following Neto et al. (2007), we assess the equilibrium state of each halo by measuring: (i) the substructure mass fraction, that is, the total mass fraction in resolved substructures whose centres lie within ⁠; (ii) the centre of mass displacement, s= |r_c−r_cm |/r₂₀₀, defined as the normalized offset between the location of the minimum of the potential and the barycentre of the mass within r₂₀₀; and (iii) the virial ratio, 2T/|U|, where T is the total kinetic energy of the halo particles within r_vir and U is their gravitational self-potential energy.

Using these criteria, we will consider haloes to be relaxed if they satisfy all the following conditions: f_sub < 0.1, s < 0.07, and 2T/|U| < 1.35. (See Neto et al. 2007 for full details.) We will also impose a minimum number of particles in order to be able to say something meaningful about internal halo structure. We initially set this limit at 500 particles, but following the analysis of profile shapes presented in Section 3.2 we subsequently adopt a stricter criterion of 3000 particles. We consider only relaxed haloes in this study, since only for such systems can the mass profiles of individual objects be represented accurately by a simple fitting formula with a small number of parameters. Such formulae are also useful for characterizing the average profiles of large ensembles of haloes, since the fluctuations in the individual systems then average out. Hayashi & White (2008) present such mean profiles as a function of mass for all MS haloes regardless of their dynamical state. We compare their results with our own below.

With these restrictions, including the 3000 particle limit, our final MS catalogue contains 128 233, 77 190, 30 603 and 9392 relaxed haloes at z= 0, 1, 2 and 3, respectively. (The corresponding numbers for the hMS catalogue are 8131, 6652, 4112 and 2194.) The overall fractions of these haloes that are relaxed in the MS are 78, 60, 50 and 48 per cent at these redshifts, respectively. We note that these fractions also depend on halo mass: at z= 0 about ∼85 per cent of ∼10¹²h⁻¹M_⊙ haloes are relaxed by our criteria, but only ∼60 per cent of ∼10¹⁵h⁻¹M_⊙ haloes. In order to obtain usefully large statistical samples, we restrict our analysis to the redshift range 0 ≤z≤ 3 in the following.

3 HALO DENSITY PROFILES

For each halo in the samples described in Section 2, we have computed a spherically averaged density profile by measuring the halo mass in 32 equal intervals of log₁₀(r) over the range 0 ≥ log₁₀(r/r_vir) ≥−2.5.

Profiles may also be stacked in order to obtain an average profile for haloes of similar mass. This procedure has the advantage of erasing individual deviations from a smooth profile which are typically due to the presence of substructure. Such deviations increase the (already considerable) scatter in the parameters fitted to individual profiles and may mask underlying trends in the data. The left-hand panels in Fig. 1 show the profile that results from stacking 464 haloes of mass ∼2 × 10¹⁴h⁻¹M_⊙ identified at z= 0 in the MS. We choose to show r²ρ versus r rather than ρ versus r in order to remove the main radial trend and enhance the dynamic range of the plot. Similar stacked profiles for all haloes in the MS (regardless of dynamical state) are shown by Hayashi & White (2008).

Figure 1

Left-hand panels: the stacked density profile of 464 haloes of virial mass ∼2 × 10¹⁴h⁻¹M_⊙, identified at z= 0 in the MS. The curves in the upper left-hand panel show different NFW fits obtained by varying the radial range fitted: [0.02, 1]r_vir (solid red), [0.05, 1]r_vir (dashed green line) and [0.1, 1]r_vir (dot–dashed blue line). Note the small disagreement between the actual profile shape and the NFW model (upper panels). This leads to concentration estimates which depend slightly on the innermost radius of the fit. Fits using equation (4) are more robust to such variations in fitting range, as shown in the lower left-hand panel. Right-hand panels: the residuals corresponding to the fits shown in the left-hand panels.

3.1 Deviations from NFW and concentration estimates

The concentration of a halo is defined using the scale radius of the profile, r_s. This identifies the location of the maximum of the r²ρ profile. However, as shown in Fig. 1, the peak is rather broad, leading to some uncertainty in its exact location when noise is present. Typically, this problem is addressed by fitting the numerical data to some specified functional form over an extended radial range. The curves in the top left-hand panel of Fig. 1 show the result of fitting the stacked halo profile with the NFW formula (equation 1), but varying the radial range of the fit as shown by the labels. This results in slightly different estimates for r_s, and consequently for the concentration, c_vir=r_vir/r_s. Increasing the innermost radius of the fit from 0.02 to 0.1r_vir results in a concentration estimate that decreases from ∼7.5 to ∼6.7.

This variation is a result of the slight (but significant) mismatch between the shape of the NFW profile and that of the stacked halo, as shown by the residuals in the top right-hand panel of Fig. 1. The ‘S’ shape of the residuals implies that the stacked profile steepens more gradually with radius in a log–log plot than does the NFW profile, a result that has been discussed in detail by Navarro et al. (2004), Prada et al. (2006) and Merritt et al. (2006).

These results suggest that force-fitting NFW profiles may induce spurious correlations between mass and concentration. In particular, when haloes are identified in a single cosmological simulation, the numerical resolution varies systematically with halo mass (less-massive haloes are more poorly resolved) so that the radial range available for fitting is a strong function of halo mass.

One way to address this issue is to adopt a fitting formula that better matches the mean profile of simulated haloes. As discussed by Navarro et al. (2004), Prada et al. (2006) and Merritt et al. (2006), improved fits are obtained adopting a radial density law where the logarithmic slope is a power law of radius:

3

which implies a density law of the form

4

Here, ρ₋₂ is the density at r=r₋₂. This density law was first introduced by Einasto (1965) who used it to describe the distribution of old stars within the Milky Way. For convenience, we will thus refer to it as the Einasto profile. Note that according to our definition, r₋₂ corresponds to the radius where the logarithmic slope of the density profile has the ‘isothermal’ value, −2. In this sense, r₋₂ is equivalent to the NFW scalelength, r_s, and again marks the location of the maximum in the r²ρ profile shown in Fig. 1.

The bottom left-hand panel of Fig. 1 shows that, at the cost of introducing an extra shape parameter, the fits improve to the point that the sensitivity of concentration estimates to the fitted radial range is effectively eliminated. Thus, concentrations obtained by fitting equation (4) to the stacked halo profiles are robust against variations in fitting details. Hayashi & White (2008) show that the same is true for fits to stacks of all haloes of a given mass, rather than just the relaxed haloes used to make Fig. 1 and, indeed, the α and c values they find are very similar to the values we obtain here, showing that our restriction to relaxed haloes has relatively little effect on the mean.

We conclude that accounting for the subtle difference between halo profile shape and the NFW fitting formula is worthwhile in order to avoid possible fitting-induced biases in concentration estimates. In the remainder of this paper, we quote concentrations, c_Δ=r_Δ/r₋₂, which are estimated by fitting density profiles by equation (4). We discuss in the next section how α is chosen for these fits.

3.2 The mass and redshift dependence of profile shape

The above discussion suggests that the shape parameter, α, should be used to improve the description of the typical density profiles of simulated haloes and to eliminate possible biases in estimates of their concentration. To do this, it is necessary to understand whether (and how) α varies with redshift and/or halo mass.

We explore this in Fig. 2. The left-hand panel shows how the best-fitting value of α depends on halo mass for the average profiles of relaxed haloes stacked according to their mass. We consider only haloes with at least 500 particles within the virial radius, and stacks containing at least 10 haloes. The solid and dashed curves in this plot correspond to the MS and hMS simulations, respectively.

Figure 2

Left-hand panel: the best-fitting shape parameter, α (equation 4), as a function of halo mass and redshift, after binning and stacking haloes by mass. The solid and dashed lines correspond to the MS and hMS simulations, respectively. Only results corresponding to haloes with more than ∼500 particles and to stacks with more than 10 haloes are shown. Note the good general agreement between the two simulations. Differences are substantial only where the number of particles is less than 3000. Right-hand panel: values of α for haloes with more than 3000 particles plotted as a function of the dimensionless ‘peak height’ parameter ν(M, z), defined as the ratio between the critical overdensity δ_crit(z) for collapse at redshift z and the linear rms fluctuation at z within spheres containing mass M. The larger the value of ν, the rarer and more massive the corresponding halo is. Note that this scaling accounts satisfactorily for the redshift dependence of the mass–concentration relation shown in the left-hand panel.

This figure illustrates several interesting points. In the first place, we note that there is good agreement between the two simulations for haloes which are represented by at least 3000 particles, but that systematic differences are visible when the MS haloes contain fewer particles than this. In the rest of this paper, we will thus only present results for haloes containing at least 3000 particles within the virial radius.

A second interesting point is that there are well-defined trends for α to increase with mass at each redshift, and with redshift at each mass. The (weak) trend with mass was already visible in Navarro et al. (2004) and Merritt et al. (2005), although the small number of haloes in these studies made their results rather inconclusive in the face of the large halo-to-halo scatter. The use of stacked profiles in Fig. 2, together with the much larger number of haloes in the simulations used here, leads to a far more convincing demonstration of the trends than was previously possible.

The trend with redshift at a given mass may seem surprising, but its interpretation is made clear by the right-hand panel of Fig. 2. Here, we show α as a function of a dimensionless ‘peak-height’ parameter, defined as ν(M, z) ≡δ_crit(z)/σ(M, z); that is, as the ratio of the linear density threshold for collapse at z to the rms linear density fluctuation at z within spheres of mean enclosed mass M. The parameter ν(M, z) is related to the abundance of objects of mass M at redshift z. ν(M_*, z) = 1 defines the characteristic mass, M_*(z), of the halo mass distribution at redshift z. ν(M, z) ≫ 1 then corresponds to rare objects with M≫M_*, while ν(M, z) ≪ 1 corresponds to objects in the low-mass tail of the distribution. The parameter ν plays a key role in the Press–Schechter model for the growth of non-linear structure (see e.g. Lacey & Cole 1993).

The right-hand panel of Fig. 2 shows that all curves coincide when expressed in terms of ν(M, z) (and when considering haloes with N > 3000 particles). Thus, the redshift dependence in the left-hand panel merely reflects the fact that objects of given virial mass correspond to very different values of ν at different redshifts. It is interesting that the dependence of α on ν is very weak for ν < 1 (it is nearly constant at α∼ 0.16), but it increases rapidly for rarer, more massive objects, reaching α∼ 0.3 for ν∼ 3.5. A simple formula

5

describes the numerical results quite accurately.

Taylor & Navarro (2001), Austin et al. (2005) and Dehnen & McLaughlin (2005) investigated a halo model based on the assumption that the phase-space density, ρ/σ³, was a simple power law of radius. Interestingly, the density profile they found is almost identical to an Einasto profile with 0.14 < α < 0.18 and so the sharp upturn we find in α at ν > 1 could be taken as indication that such a model is not valid for rare and massive haloes.

The dependence of profile shape on ν is illustrated in Fig. 3, where we plot r²ρ profiles for halo stacks corresponding to three different values of ν: 1.0, 2.0 and 3.2. The larger ν, the larger α, and the more sharply the profile peaks. It is unclear what causes these systematic trends, but the fact that they depend on ν (rather than, say, on mass alone) is an important clue for models that attempt to explain the near-universality of dark halo density profiles. While our results here are based purely on relaxed haloes, very similar results were found by Hayashi & White (2008) for the average profiles of stacks of all MS haloes, regardless of dynamical state.

Figure 3

Three stacked halo density profiles with different values of ν, and at different redshifts, as labelled. The profiles are chosen to illustrate the variation in halo structure as a function of ν indicated in Fig. 2. The larger the halo mass, the larger is the value of α and the sharper is the curvature in the density profile as a function of radius.

3.3 Concentration estimates

As we discussed above, a fitting formula other than NFW is needed to obtain concentration estimates that are insensitive to the radial range fitted. Einasto's profile, equation (4), accomplishes this at the cost of introducing an additional shape parameter, α. Adjusting α to fit the detailed structure of individual haloes would negate the spirit of the NFW programme which attempts to predict the structure of dark haloes in any hierarchical cosmology from its initial power spectrum and the global cosmological parameters. Three-parameter fits to individual haloes are also susceptible to strongly correlated parameter errors. We have therefore explored whether fixing α as a function of halo mass and redshift according to the α(ν) relation of equation (5) results in significantly different concentration distributions than adjusting it freely to fit each halo.

We show the result in Fig. 4, which shows concentration distributions for the 464 individual haloes which were stacked to make Fig. 1. We compare the distribution obtained from full three-parameter Einasto fits to that obtained when α is set to the value predicted for each halo from its mass. We find that for most haloes the fits give essentially the same concentration, and the distributions of concentration are almost indistinguishable. The medians of the fixed-α and floating-α distributions are 6.85 and 6.90, respectively. The scatters also coincide, as noted in the labels of Fig. 4. Interestingly, the scatter around the best-fitting profile is typically only slightly smaller for the Einasto model than for the NFW model, showing that the difference between the two models is much smaller than the deviation of the profile of typical individual haloes from either model. In this sense, the three-parameter Einasto profile has little advantage over the two-parameter NFW profile when estimating dark matter halo concentrations.

Figure 4

Distributions of concentration parameters estimated using equation (4) for the 464 individual halo profiles stacked in Fig. 1. The solid (black) histogram corresponds to fits where α was adjusted separately for each individual halo, and the dashed (red) histogram corresponds to fits where α was set to the value implied by the α(ν) relation of Fig. 2 (equation 5). The numbers in the legend give the median and scatter of the distributions. The excellent agreement between the two distributions indicates that unbiased and accurate concentration estimates for relaxed haloes may be obtained by fixing α to the value determined by equation (5).

This is easily understood. Estimating the concentration of a halo is equivalent to locating the ‘peak’ of the r²ρ profile. Provided that the shape of the fitted profile is, on average, a good approximation to the simulated ones, no systematic offset is expected between concentrations estimated using the two procedures. We conclude that concentrations may be estimated robustly by fitting Einasto profiles to individual haloes with α fixed to the value obtained from equation (5). All values reported below were obtained using this procedure.

3.4 The mass and redshift dependence of halo concentration

The concentration–mass relation is shown in Fig. 5 for z= 0, 1, 2 and 3. The upper left-hand panel is equivalent to fig. 6 of Neto et al. (2007), except that our concentrations are estimated from fits covering the range [0.05, 1]r_vir using an Einasto rather than an NFW profile, and we only show results for haloes with more than 3000 particles. The differences are small, as may be judged from the power-law fit proposed by Neto et al., c₂₀₀= 5.26 (M₂₀₀/10¹⁴h⁻¹M_⊙)^−0.10, shown as a dotted line in the upper left-hand panel of Fig. 5. This power law fit is very similar (after correction for differing concentration and mass definitions) to that which Macciò et al. (2007) obtained from NFW fits to haloes in an independent set of (smaller) simulations, despite the fact that these authors included haloes with as few as 250 particles, which we would consider to be significantly underresolved on the basis of our own tests.

Figure 5

The mass dependence of concentration as a function of redshift. Concentration estimates are derived from Einasto fits to the density profile of ‘relaxed’ haloes over the radial range [0.05, 1]r_vir. The points, boxes and whiskers show the median and the 5, 25, 75 and 95 percentiles of the concentration distribution within each mass bin. The black symbols correspond to the MS, and red symbols to the smaller, but higher resolution hMS. Only results for haloes with more than 3000 particles are shown for each simulation. The solid dotted, dashed, dot–dashed and solid lines show the models of NFW, ENS, B01 and the modified NFW, respectively, for comparison with our results. The lower dotted line in the upper left-hand panel indicates the relation obtained by Neto et al. (2007) from NFW fits to relaxed haloes at z= 0. The remaining dotted lines show the power-law fits given in Table 1.

Both the concentration values and the trends with mass and redshift are very similar to those presented in fig. 2 of Zhao et al. (2003b) who analysed a set of ΛCDM simulations of varying resolution. The small offsets between their mean concentrations and our results are consistent with the slight differences we have noted when switching from NFW to Einasto models for determining concentrations. In addition to confirming these earlier results, the much larger volume of our simulations results in a better determination of the intrinsic scatter about the mean relation.

The red and black symbols in Fig. 5 correspond to the MS and hMS simulations, respectively, with the dots plotted at the median concentration for each mass bin. The boxes represent the lower and upper quartiles of the concentration distributions, while the whiskers show their 5 and 95 per cent tails. We also show the concentrations predicted by three previously proposed analytic models: NFW (solid dotted magenta lines), B01 (dot–dashed black lines) and ENS (dashed green lines). As discussed by Neto et al. (2007), none of these models reproduces the simulation results over the full mass range accessible at z= 0: the NFW and ENS models underestimate the concentration of galaxy-sized haloes, whereas the B01 model fails dramatically at the high-mass end, where it predicts a sharp decline which is not seen in the simulations.2

There is a hint in the z= 0 panel that the relation is flattening at the high-mass end, where the NFW predictions at z= 0 appear slightly better than those of ENS. This is because a constant concentration for very rare and massive objects is implicit in the NFW model, which assumes that the characteristic density of a halo reflects that of the universe at the time it collapsed. Very massive systems assembled very recently, and therefore share the same collapse time (i.e. they are being assembled today).

The near-constant concentration of the most massive haloes is considerably more obvious at higher redshift, presumably because well-resolved haloes (i.e. those with N > 3000 in the MS or hMS) become rarer and rarer with increasing lookback time. Indeed, whereas at z= 0 our MS halo catalogue spans the range 0.75 < ν < 3, at z= 3 all the haloes retained have ν≳ 3. As a result, the average concentration at z= 3 is almost independent of mass over the accessible mass range, that is, M≳ 3 × 10¹¹h⁻¹M_⊙. It is interesting that the average concentration of the most massive haloes (i.e. ν≥ 3) is similar at all redshifts, c∼ 3.5–4. This evokes the proposals of Zhao et al. (2003a,b), Tasitsiomi et al. (2004) and Lu et al. (2006), all of who argue that haloes undergoing rapid growth should all have similar concentration.

The evolution of the mass–concentration relation seen in our numerical simulations is not predicted by any of the published prescriptions. The original NFW model shows a flattening of the concentration–mass relation with increasing redshift, similar to that observed in the simulations, but it predicts insufficient evolution in concentration at a given mass. As a result, this model overestimates the concentrations by an increasing amount with increasing redshift, about 40 per cent at z= 3. The ENS and B01 models fail to reproduce the features of the mass–concentration–redshift relations at high mass, predicting a stronger mass dependence and much more evolution than is seen. In these two models, the concentration of haloes of given mass scales inversely as the expansion factor, so that shape of the mass–concentration relation remains fixed. While the high-mass discrepancy between the B01 model and our measurements can be reduced by changing the parameters to K= 2.8 and F= 0.001 (see Wechsler et al. 2006), neither for this model nor for the ENS model can parameter changes produce agreement with the weak redshift evolution seen at high mass both here and by Zhao et al. (2003b). On the other hand, the ENS and B01 models predict the concentration evolution of galaxy mass haloes substantially better than the NFW model. Note that our simulation data do not disagree significantly from those analysed by ENS and B01. The discrepancies result from extrapolation of their proposed relations beyond the range where they were reliably tested.

In the NFW model, the definition of formation time involves two physical parameters, F and f, and a proportionality constant, C, that relates the value of the characteristic halo density to the mean density of the universe at the time of collapse. The formation redshift of a halo is taken to be the redshift at which a fraction, F, of its mass was first contained in progenitors each individually containing at least a fraction f of its mass. In the original NFW prescription, F= 0.5, f= 0.01 and C= 3000. We find that the observed trend in the slope of the mass–concentration relation with redshift can be approximately reproduced by simply changing F to F= 0.1, keeping f= 0.01 as before; the normalization of the relation is then approximately reproduced by taking C= 600. The resulting curves are shown as the solid blue lines in Fig. 5. This modified NFW model succeeds well in matching the redshift evolution, but its mass dependence is still too shallow at z= 0 and too steep at z= 3, leaving room for improvement, especially at low masses where the B01 and ENS prescriptions give better predictions of the evolution rate. We have checked that the same model also gives an acceptable fit for other cosmologies, in particular for a simulation of a ΛCDM model similar to hMS but with the values of the cosmological parameters inferred from the three-year WMAP satellite data (Spergel et al. 2007) (Ω= 0.236, Ω_Λ= 0.764, n= 0.97 and σ₈= 0.74). We have also checked that this modified NFW model gives a good description of the concentrations found in the scale-free models presented in Navarro et al. (1997) which span a range of spectral indices and have either Ω= 1 or 0.1. If more accurate results for the particular cosmology assumed in this paper are desired, we provide the coefficients of power-law fits of the form

6

in Table 1 and plot these best-fitting curves in Fig. 5.

Table 1

Values of the constants A and B for the best straight-line fit (6) to the data shown in Fig. 5.

Redshift	A	B
0	−0.138	2.646
0.5	−0.125	2.372
1	−0.092	1.891
2	−0.031	0.985
3	−0.004	0.577

The disagreement between the simulation data and the original NFW prescription is up to a few tens of per cent. The BO1 and ENS prescriptions underestimate the concentration of the most massive haloes by factors of ∼2 to ∼3 at high redshift. Since the characteristic density of a halo scales roughly as the third power of the concentration (see equation 2), this implies that the characteristic densities of such haloes are underestimated by at least an order of magnitude by the latter two models, leading to dramatic changes in the expected lensing power, X-ray luminosity, and SZ detectability of massive clusters at high redshift.

It is also interesting to compare the z= 0 concentration–mass relation found here to the one which Hayashi & White (2008) obtained by fitting Einasto profiles to stacks of all haloes of a given mass, rather than to stacks of relaxed haloes. The results in their fig. 4 lie very close to the NFW relation plotted in the upper left-hand panel of our own Fig. 5. Thus, the restriction to relaxed haloes has little systematic effect on Einasto-based concentrations at masses above about 3 × 10¹³h⁻¹M_⊙, but at lower masses it results in somewhat larger concentrations. This differs slightly from the result of Neto et al. (2007) who found z= 0 concentrations obtained by fitting NFW (rather than Einasto) profiles to relaxed haloes to be greater by a larger amount at high mass. The median of their NFW-based concentrations is also closer to the NFW model prediction on galaxy scales (both for relaxed and for all haloes) than is the median of the Einasto-based concentrations which we plot in Fig. 5. This is at least in part a result of the bias illustrated in Fig. 1.

Neto et al. (2007), Hayashi & White (2008) and this paper are all based on the same simulations (although Hayashi & White do not use the hMS). The small but (statistically) significant differences in their derived concentration–mass relations show that at the 10–20 per cent level these relations depend on the details of the halo sample selection and profile-fitting procedures. The modified NFW model works reasonably well over the full range of mass and redshift studied here, certainly rather better than the revisions proposed by ENS and B01. Significant discrepancies remain, however, notably at galaxy masses where the NFW model underpredicts the rate of evolution, resulting in systematically low Einasto-based concentrations at z∼ 0 and systematically high values at z≥ 2. For such objects the evolution rate is better predicted by the ENS and B01 models. Given the good agreement between various simulations on this mass scale (see e.g. Zhao et al. 2003b; Macciò et al. 2007; Neto et al. 2007), we believe these concentration estimates for low-redshift ΛCDM galaxy haloes to be accurate. This implies that the NFW model needs further improvement for such objects. At higher masses and higher redshifts, however, the (modified) NFW formalism works quite well, and should be preferred to those of B01 or ENS.

With our modified parameters, the NFW prescription not only fits concentrations in the WMAP1 cosmology investigated here, but also our (less-extensive) set of data for the WMAP3 cosmology. Without detailed testing, however, it is unclear if the prescription can be extended to other regimes of interest. For example, what concentrations are expected for dwarf galaxy haloes or for ‘first object’ haloes at z∼ 30? The incorrect conclusions reached by applying the B01 and ENS formulae to massive clusters at high redshift are testimony to the dangers of using empirical formulae outside the range where they have been tested. Clearly, further theoretical effort aimed at understanding the factors which determine the internal structure of dark haloes would be a welcome complement to the numerical work presented here.

4 SUMMARY

We have used data from the MS and from a smaller but higher resolution simulation to examine how the density profiles of relaxed ΛCDM haloes vary with halo mass and redshift. We study profile shape and concentration for haloes with mass exceeding ∼3 × 10¹¹h⁻¹M_⊙ over the redshift range 0 ≤z≤ 3. Our main results may be summarized as follows.

• We confirm the conclusion of previous studies that, in the mean, the shape of spherically averaged ΛCDM halo density profiles deviates slightly but systematically from the two-parameter fitting formula proposed by Navarro et al. (1995, 1996, 1997). A more accurate description is provided by the three-parameter Einasto (1965) profile, for which the logarithmic slope is a power law of radius, d log ρ/d log r∝r^α. We show that this fitting formula gives concentration estimates which are insensitive to the radial range fitted, albeit at the price of an additional shape parameter. Although Einasto fits avoid small but significant biases that arise if all haloes are fitted to the NFW model, we emphasize that individual haloes typically deviate from either model by more than the difference between them.

• Using stacked profiles of haloes of similar mass, we show that the shape parameter, α, depends systematically on halo mass and on redshift. These dependencies can be collapsed into a dependence on the single ‘peak height’ parameter, ν(M, z) =δ_crit(z)/σ(M, z). Haloes with large ν are rare objects in the high-mass tail of the halo mass distribution and have large α values (see equation 5). This provides an important clue for models that attempt to interpret the dependence of halo density profiles on mass and redshift.

• The dependence of halo concentration on mass becomes progressively weaker with increasing redshift, as found earlier by Zhao et al. (2003b). At z= 3, concentrations are almost independent of mass over the mass range accessible to our simulations. Relaxed haloes with ν > 3 (the rarest and most massive systems) have similar concentrations, 〈c₂₀₀〉= 3.5–4, at all the redshifts we have studied.

• The models of B01 and ENS fail to reproduce our measured concentrations for high-mass and high-redshift objects, predicting a stronger mass dependence and more evolution than is seen in the simulations. Parameters in these models can be changed to reduce the strength of their mass dependence, but the predicted redshift evolution, while fitting galaxy mass haloes well up to redshift z= 1, remains substantially too strong at high mass and at higher redshifts. As a result, the predictions of these models for high-redshift galaxy clusters can be in error by large factors. The original model of Navarro et al. (1997) overpredicts the concentrations of such objects at redshifts beyond 1 (by up to ∼40 per cent at z= 3) but a modified NFW model with a different definition of formation redshift reproduces the simulation results substantially better over the redshift and mass ranges we have examined. Both the original and the modified NFW models underestimate the concentration evolution of relaxed 10¹²h⁻¹M_⊙ haloes, leading to 30 per cent discrepancies at z= 0 and z= 3.

We hope that our simulation results will stimulate theoretical work aimed at a deeper understanding of the factors which determine the internal structure of ΛCDM haloes. Such work may eventually result in simple recipes like those of NFW, ENS and B01. Substantial errors are found, however, when these published prescriptions are extrapolated beyond the regime where their authors tested them, in particular, to the regime relevant to high-redshift galaxy clusters. This demonstrates that careful numerical work is mandatory before any recipe can be applied in a new regime. When making forecasts for surveys of distant massive clusters, our simulations show our modified NFW recipe to give reliable results at least out to z= 3.

1

We express Hubble's constant as H(z) and its present-day value as H(z= 0) =H₀= 100 h km s⁻¹ Mpc⁻¹.

2

We note that NFW, ENS and B01 parametrize the initial power spectrum in slightly different ways, and that the predicted concentrations are sensitive to the exact choice of parameters. For example, NFW and ENS use the parameter Γ to characterize the shape of the linear power spectrum. We use Γ= 0.15 here since that provides the best fit to the actual power spectrum used in the MS. For B01 we have used the default values advocated in the latest version of their software (K= 2.8 and F= 0.001 in their notation), which is available from http://www.physics.uci.edu/~bullock/CVIR/. The differences between the predictions shown here and in fig. 4 of Hayashi & White (2008), or in Neto et al. (2007), for example, are due to slight differences in the values chosen for these parameters.

The MS was carried out as part of the programme of the Virgo Consortium on the Regatta supercomputer of the Computing Centre of the Max-Planck Society in Garching. The hMS simulation was carried out using the Cosmology Machine at Durham. We thank Adam Mantz for pointing out an error in the NFW97 curve plotted in the original preprint version of Fig. 5. JFN acknowledges support from the Alexander von Humboldt Foundation and from the Leverhulme Trust, as well as the hospitality of the Max-Planck Institute for Astrophysics in Garching, Germany, and the Institute for Computational Cosmology in Durham. CSF acknowledges a Royal Society Wolfson Research Merit Award. This work was supported in part by the PPARC Rolling Grant for Extragalactic Astronomy and Cosmology at Durham.