The shape of dark matter haloes: dependence on mass, redshift, radius and formation

Abstract

Using six high-resolution dissipationless simulations with a varying box size in a flat Lambda cold dark matter (ΛCDM) universe, we study the mass and redshift dependence of dark matter halo shapes for M_vir= 9.0 × 10¹¹− 2.0 × 10¹⁴h⁻¹ M_⊙, over the redshift range z= 0–3, and for two values of σ₈= 0.75 and 0.9. Remarkably, we find that the redshift, mass and σ₈ dependence of the mean smallest-to-largest axis ratio of haloes is well described by the simple power-law relation 〈s〉= (0.54 ± 0.02)(M_vir/M_*)^{−0.050±0.003}, where s is measured at 0.3R_vir, and the z and σ₈ dependences are governed by the characteristic non-linear mass, M_*=M_*(z, σ₈). We find that the scatter about the mean s is well described by a Gaussian with σ∼ 0.1, for all masses and redshifts. We compare our results to a variety of previous works on halo shapes and find that reported differences between studies are primarily explained by differences in their methodologies. We address the evolutionary aspects of individual halo shapes by following the shapes of the haloes through ∼100 snapshots in time. We determine the formation scalefactor a_c as defined by Wechsler et al. and find that it can be related to the halo shape at z= 0 and its evolution over time.

galaxies: formation, galaxies: haloes, cosmology: theory, large-scale structure of Universe

1 INTRODUCTION

A generic prediction of cold dark matter (CDM) theory is the process of bottom up dark matter halo formation, where large haloes form from the mergers of smaller haloes, which are in turn formed from even smaller haloes. This is often a violent process violating most of the assumptions that go into the spherical top-hat collapse model of halo formation often used to describe haloes. Since mass accretion on to haloes is often directional and tends to be clumpy, one should not expect haloes to be spherical, if the relaxation time of the haloes were longer than the time between mergers and/or the infalling haloes came along a preferential direction (such as along a filament). In both theoretical modelling of CDM and observations, haloes are found to be very non-spherical. In fact, spherical haloes are rare. Therefore, the analysis of halo shapes may shed light on the nature of the dark matter and the process of halo and galaxy formation.

In order to more precisely approximate and describe the shape of haloes, going one step beyond the spherical approximation, haloes can be described by associating ellipsoids with their mass distribution. Ellipsoids are characterized by three axes, a, b, c, with a≥b≥c, which are normally represented in terms of scale-independent ratios, s≡c/a, q≡b/a and p≡c/b. There are three classes of ellipsoids which have corresponding ratio ranges: prolate (sausage shaped) ellipsoids have a >b≈c leading to axial ratios of s≈q < p, oblate (pancake shaped) ellipsoids have a≈b > c leading to axial ratios of s≈p < q and triaxial ellipsoids are in between prolate and oblate with a > b > c. Additionally, when talking about purely oblate ellipsoids, a=b, it is common to just use q, since s and q are degenerate.

There have been many theoretical papers published over the years examining the subject of halo shapes. The early works on the subject include Barnes & Efstathiou (1987), Dubinski & Carlberg (1991), Katz (1991), Warren et al. (1992), Dubinski (1994), Jing et al. (1995), Tormen (1997) and Thomas et al. (1998). All of these works agreed that haloes are roughly ellipsoidal, but otherwise differed on the details. Dubinski & Carlberg (1991) found that haloes have axial ratios of s∼ 0.5 in the interior and become more spherical at larger radii, while Frenk et al. (1988) found that haloes are slightly more spherical in the centres. Thomas et al. (1998) claimed that larger mass haloes have a slight tendency to be more spherical, whereas more recent simulation results found the opposite. Despite these disagreements, many of the early authors give us clues about the nature of haloes' shapes. Warren et al. (1992) and Tormen (1997) showed that the angular momentum axis of a halo is well correlated with the smallest axis, c, although, as most of these authors pointed out, haloes are not rotationally supported. This therefore has led many to conclude that the shapes are supported by anisotropic velocity dispersion. Tormen (1997) took it one step further and found that the velocity anisotropy was in turn correlated with the infall anisotropy of merging satellites. Most authors found that the axial ratios of haloes are ∼0.5 ± 0.1 and that haloes tend to be prolate as opposed to oblate in shape. The most likely source of disagreement in these works and in the more recent works described below is the different methods used often coupled with inadequate resolution.

More recent studies of haloes' shape were performed by Bullock (2002), Jing & Suto (2002) (JS), Springel, White & Hernquist (2004), Bailin & Steinmetz (2005), Kasun & Evrard (2005) and Hopkins, Bahcall & Bode (2005). The results of these authors differ, in some cases, considerably. One of the goals of this paper is to carefully examine the differences in the findings presented by the above mentioned authors.

All of the aforementioned publications [except Springel et al. (2004)] analyse simulations with either no baryons or adiabatic hydrodynamics. This is both due to the cost associated with performing self-consistent hydrodynamical simulations of large volumes with high mass resolution and due to the fact that very few cosmological simulations yet produce realistic galaxies. None the less, we know that the presence of baryons should have an effect on the shapes of haloes due to their collisional behaviour. Three recent papers have attempted to examine the effect of baryons (Kazantzidis et al. 2004; Springel et al. 2004; Bailin et al. 2005). In Springel et al. (2004), the same simulations were performed using no baryons, adiabatic baryons and baryons with cooling and star formation. In the first two cases, there was very little difference, but with the presence of cooling and star formation the haloes became more spherical. The radial dependence of shape also changes such that the halo is more spherical in the centre. At R > 0.3R_vir, the axial ratio s increases by ≤0.09, but in the interior the increase Δs is as much as 0.2. In an independent study, Kazantzidis et al. (2004) found an even larger effect due to baryonic cooling in a set of 11 high-resolution clusters. At R= 0.3R_vir, the authors found that s can increase by 0.2 − 0.3 in the presence of gas cooling. The extent of the over-cooling problem plaguing these simulations is still uncertain. For this reason, the amount of change in the shape should be viewed as an upper limit. The most recent work on the subject is by Bailin et al. (2005), who concentrate more on the relative orientation of the galaxy formed at the centre of eight high-resolution haloes than on the relative sphericity of the haloes. Despite this, from fig. 1 of Bailin et al. (2005), it seems that they would also predict an increase of ∼0.2 for s. In light of these studies, it is still useful to study shapes of haloes without baryonic cooling. Cooling and star formation in simulations is still a very open question, making the effect of the cooling uncertain, and we show in Flores et al. (2005), hereafter referred to as Paper II, that our simulations without cooling match shapes of X-ray clusters.

Measurements of the shapes of both cluster and galaxy mass haloes through varied observational techniques are increasingly becoming available. There have been many studies of the X-ray morphologies of clusters (see McMillan, Kowalski & Ulmer 1989; Mohr et al. 1995; Kolokotronis et al. 2001) which can be directly related to the shape of the inner part of the cluster halo (Buote & Xu 1997; Lee & Suto 2003). For a review of X-ray cluster shapes and the latest results, see Paper II.

In addition, there has been an addition of important new information on the shape of galaxy mass haloes, in particular our own Milky Way (MW) halo. Olling & Merrifield (2000) concluded through the examination of disc warping that the host halo around the Galactic disc is oblate with a short-to-long axial ratio of 0.7 < q < 0.9. Investigations of Sagittarius' tidal streams have led to the conclusion that the MW halo is oblate and nearly spherical with q≳ 0.8 (Ibata et al. 2001; Majewski et al. 2003; Martínez-Delgado et al. 2004). However by inspecting M giants within the leading stream, Helmi (2004) and Law, Johnston & Majewski (2005) found a best-fitting prolate halo with s= 0.6. Merrifield (2004) summarizes the currently reliable observations for galaxy host halo shapes using multiple techniques and find that the observations vary a lot.

Another method for studying the shapes of haloes at higher redshift is galaxy–galaxy weak lensing studies. Analysing data taken with the Canada–France–Hawaii telescope, Hoekstra, Yee & Gladders (2004) find a signal at a 99.5 per cent confidence level for halo asphericity. They detect an average projected ellipticity of 〈ε〉≡〈1 −q_2D〉= 0.20^+0.04_−0.05, corresponding to s= 0.66^+0.07_−0.06, for haloes with an average mass of 8 × 10¹¹h⁻¹ M_⊙. Ongoing studies of galaxy–galaxy weak lensing promise rapidly improving statistics from large-scale surveys such as the Canada–France Legacy Survey.

This paper is organized as follows. In Section 2, we describe the simulations, halo finding method and halo property determination methods used in this study. In Section 3, we discuss the method used to determine the shapes of haloes. In Section 4, we examine the mean halo axial ratios from our simulations and the ratios' dependence on mass, redshift and σ₈. We then examine the dispersion of the axial ratio versus mass relation. We briefly discuss the shape of haloes as a function of radius, and then examine the relationship of the angular momentum and velocity anisotropy to the halo shape. In Section 5, we examine the relationship between the formation history of haloes and their present-day shapes. In Section 6, we compare our results to those of previous authors and explain the sources of the differences. In Section 7, we examine the observational tests and implications of our findings. Finally, Section 8 is devoted to summary and conclusions.

2 SIMULATIONS

2.1 The numerical simulations

All of the simulations (see Table 1) were performed with the Adaptive Refinement Tree (ART) N-body code of Kravtsov, Klypin & Khokhlov (1997) which implements successive refinements in both the spacial grid and time-step in high-density environments. We analyse the shapes of haloes and their merger histories in the concordance flat ΛCDM cosmological model: Ω_m= 0.3 = 1 −Ω_Λ, h= 0.7, where Ω_m and Ω_Λ are the present-day matter and vacuum energy densities in units of critical density and h is the Hubble parameter in units of 100 km s⁻¹ Mpc⁻¹. The power spectra used to generate the initial conditions for the simulations were determined from a direct Boltzmann code calculation (courtesy of Wayne Hu).

Table 1

Simulation parameters.

Name	σ8	Ωb	Lbox (h−1 Mpc)	Np	mp (h−1 M⊙)	hpeak (h−1 kpc)	M*(1012h−1 M⊙)
							z= 0	z= 1	z= 2	z= 3
L800.75	0.75	0.030	80	5123	3.16 × 108	1.2	3.0	0.11	0.0046	0.00027
L800.9a	0.9	0.045	80	5123	3.16 × 108	1.2	8.0	0.35	0.019	0.0013
L800.9b	0.9	0.045	80	5123	3.16 × 108	1.2	8.0	0.35	0.019	0.0013
L2000.9	0.9	0.030	200	2563	3.98 × 1010	5.0	8.6	0.41	0.023	0.0018
L1200.9	0.9	0.045	120	5123	1.07 × 109	1.8	8.0	0.35	0.019	0.0013
L1200.9r	0.9	0.045	40 sphere	∼2563	1.33 × 108	0.9	8.0	0.35	0.019	0.0013

To study the effects of the power spectrum normalization and resolution, we consider five simulations of the ΛCDM cosmology. The first simulation (L80_0.75) followed the evolution of 512³= 1.34 × 10⁸ particles in a 80 h⁻¹ Mpc = 114.29 Mpc box and was normalized to σ₈= 0.75, where σ₈ is the rms fluctuation in spheres of 8 h⁻¹ Mpc comoving radius. The second simulation (L80_0.9) is an exact replica of the L80_0.75 simulation with the same random number seed, but the power spectrum was normalized to σ₈= 0.9. The first simulation was also used to study the halo occupation distribution and the physics of galaxy clustering by Kravtsov et al. (2004) and Zentner et al. (2005). Unfortunately, both of these simulations were generated with a power spectrum which had a little more than average power on large scales. This may happen when generating power spectra due to cosmic variance. The simulation is still a good representation of a volume in the Universe, but to avoid becoming non-linear on large scales, the second simulation was stopped at z= 0.1. Due to the lower normalization of the L80_0.75 box, it was allowed to run until z= 0. We use these two simulations to study the effects of the spectrum normalization, but to achieve better statistics and make predictions for σ₈= 0.9 at z= 0 we also include another simulation of the same size and resolution (L80_0.9b). The fourth simulation (L200_0.9) followed the evolution of 256³= 1.68 × 10⁷ particles in a 200 h⁻¹ Mpc = 285.7 Mpc box. The fifth simulation L120_0.9 followed the evolution of 512³ particles in a 120 h⁻¹ Mpc = 171.43 Mpc box and was normalized to σ₈= 0.9. This simulation is used for several purposes: first to achieve better statistics for rare high-mass objects, and secondly as the basis for the sixth simulation. The sixth simulation is a resimulated Lagrangian subregion of the L120_0.9 box corresponding at z= 0 to a sphere in position space of diameter D= 40 h⁻¹ Mpc. The initial conditions of the L120_0.9 box contain 1024³ particles which were combined into the 512³ particles used for the fifth simulation. The Lagrangian subregion was then chosen and the original higher resolution particles of mass m_p= 1.33 × 10⁸h⁻¹ M_⊙ within this region, corresponding to 1024³ particles in the box, were followed from the initial time-step, z_i= 40. The high mass resolution subregion was surrounded by layers of particles of increasing mass, with a total of five particle species in order to preserve the large-scale gravitational field. Only the regions containing the highest resolution particles were adaptively refined. The maximum level of refinement in the simulation corresponded to a peak formal spatial resolution of 0.9 h⁻¹ kpc. For more details about the multi-mass technique, consult Klypin et al. (2001). The subregion was chosen to not contain haloes above M_vir > 10¹³h⁻¹ M_⊙ in order to increase the statistics of isolated galaxy mass haloes.

2.2 Halo identification and classification

A variant of the Bound Density Maximum (BDM) algorithm is used to identify haloes and subhaloes in our simulations (Klypin et al. 1999). The details of the algorithm and parameters being used in the halo finder can be found in Kravtsov et al. (2004). We briefly describe the main steps in the halo finder here. First, all particles are assigned a density using the smooth algorithm,1 which uses a symmetric smoothed particle hydrodynamics (SPH) smoothing kernel on the 32 nearest neighbours. Density maxima which are separated by a minimum distance of r_min= 50 h⁻¹ kpc are then identified. r_min defines the minimum distinguishable separation of haloes and subhaloes in the catalogue. Using the maxima as centres, profiles in circular velocity and density are calculated in spherical bins. Unbound particles are removed iteratively as described in Klypin et al. (1999). The halo catalogue used is complete for haloes with ≳50 particles. This corresponds to a mass below which the cumulative mass and velocity functions begin to flatten (see Kravtsov et al. 2004 for details).

The halo density profiles are constructed using only bound particles, and they are fit by an Navarro, Frenk and White (NFW) profile (Navarro, Frenk & White 1996):

1

where r_s is the radius at which the log density profile has a slope of −2 and the density is ρ_s/4. One of the parameters, r_s or ρ_s, can be replaced by a virial parameter (R_vir, M_vir, or V_vir) defined such that the mean density inside the virial radius is Δ_vir times the mean universal density ρ_o(z) =Ω_m(z)ρ_c(z) at that redshift:

2

The parameter ρ_c(z) is the critical density, and

3

is the virial overdensity from Bryan & Norman (1998) with Δ_vir(0) ≈ 337 for the ΛCDM cosmology assumed here. The NFW density profile fitting is performed using a χ² minimization algorithm. The profiles are binned logarithmically from twice the resolution length (see Table 1) out to R₅₀₀, the radius within which the average density is equal to 500 times the critical density of the universe. The choice of this outer radius is motivated by Tasitsiomi et al. (2004) who showed that haloes are well relaxed within this radius. The binning begins with 10 radial bins, with the number of bins being reduced if any bin contains fewer than 10 particles or is radially smaller than the resolution length. This reduction of bins is continued until both criteria are met. Fits using this method have been compared to fits determined using different merit functions, such as the maximum deviation from the fit as described in Tasitsiomi et al. (2004), and it was found that they give very similar results for individual haloes. After fitting the haloes, the host halo and subhalo relationship are determined very simply. If a halo's centre is contained within the virial radius of a more massive halo, that halo is considered a subhalo of the larger halo. All halo properties reported here are for haloes which are determined to be isolated or host haloes (i.e. not subhaloes).

3 METHODS OF DETERMINING SHAPES

There are a few different ways found in the literature to model haloes as ellipsoids. They all differ in details, but most methods model haloes using the eigenvectors from some form of the inertia tensor. The eigenvectors correspond to the direction of the major axes, and the eigenvalues to the length of the semimajor axes c≤b≤a, from which the axial ratios s≡c/a and q≡b/a are defined. The two forms of the inertia tensor used in the literature to determine shape are the unweighted

4

and the weighted (or reduced),

5

where

6

is the elliptical distance in the eigenvector coordinate system from the centre to the nth particle. In both the cases, the eigenvalues (λ_a≤λ_b≤λ_c) determine the axial ratios described at the beginning of Section 1 with ⁠. The orientation of the halo is determined by the corresponding eigenvectors.

One would like to recover the shape of an isodensity surface. The method used here begins by determining with s= 1 and q= 1, including all particles within some radius. Subsequently, new values for s and q are determined, and the volume of analysis is deformed along the eigenvectors in proportion to the eigenvalues. There are two options to choose from when deforming the volume. The volume within the ellipsoid can be kept constant, or one of the eigenvectors can be kept equal to the original radius of the spherical volume. In our analysis of shapes, the longest axis is kept equal to the original spherical radius. After the deformation of the original spherical region, is calculated once again, but now using the newly determined s and q and only including the particles found in the new ellipsoidal region. The iterative process is repeated until convergence is achieved. Convergence is achieved when the variance in both axial ratios, s and q, is less than a given tolerance.

The analysis presented here begins with a sphere of R= 0.3R_vir, unless otherwise stated, and keeps the largest axis fixed at this radius. For determining halo shapes accurately, we limit our analysis to isolated haloes with N_p≥ 7000 within R_vir. This corresponds to M_vir≥ 2.21 × 10¹²h⁻¹ M_⊙ for the 80 h⁻¹ Mpc box simulations, M_vir≥ 7.49 × 10¹²h⁻¹ M_⊙ for the 120 h⁻¹ Mpc box simulation and M_vir≥ 9.3 × 10¹¹h⁻¹ M_⊙ in the resimulated region of the 120 h⁻¹ Mpc box. For a discussion of our resolution tests, see Appendix A.

4 SHAPES AS A FUNCTION OF HALO MASS

The simulations examined here enable us to analyse haloes spanning a mass range from galaxy to cluster-sized objects. These data provide an opportunity to study the variation of shape and its intrinsic scatter with halo mass. Various statistics are used to derive robust estimates of the dependence of shape on halo-centric radius. Combining the detailed spatial and dynamical information from the simulations, we can relate quantities such as angular momentum or velocity anisotropy tensor to the shape and the orientation of the halo. In this section, we aim to present a comprehensive analysis of the properties of haloes at all redshifts. In Section 5, we will address evolutionary aspects of individual halo shapes.

4.1 Median relationships for distinct haloes

We begin by fitting the mass dependence of halo shape and find that the mean value of the axial ratio s≡c/a decreases monotonically with increasing halo mass as illustrated in Fig. 1. In other words, less massive haloes have a mean shape which is more spherical than that of more massive haloes. Since we use four different simulations (L80_0.9b, L120_0.9, L120_0.9r and L200_0.9) with varying mass and length-scales, we are able to determine 〈s〉(M_vir) over a wide mass range. We find that over the accessible mass range, the variation of the mean shape parameters with halo mass is well described by

7

with best fit values

8

The parameters, α and β, were determined by weighted χ² minimization on the best fit mean data points determined via Kolmogorov–Smirnov (KS) analysis assuming a Gaussian distribution within a given mass bin (see Section 4.4). M_*(z) is the characteristic non-linear mass at z such that the rms top-hat smoothed overdensity at scale σ(M_*, z) is δ_c= 1.68. The M_* for z= 0 is 8.0 × 10¹²h⁻¹ M_⊙ for the simulations with Ω_b= 0.045 and 8.6 × 10¹²h⁻¹ M_⊙ for the simulations with Ω_b= 0.03. Only bins containing haloes above our previously stated lower bound resolution limit were used and only mass bins with at least 20 haloes were included in the fit. This work extends the mass range of the similar relationships found by previous authors (Bullock 2002; Kasun & Evrard 2005; Springel et al. 2004; JS); we compare our results with these previous works in Section 6.

Figure 1

Mean axial ratios s=c/a for four simulations of different mass resolution are presented with a fit (solid black lines) given by equation (7) and dispersion of 0.1. The triangles, squares, solid circles and × symbols are the average s for a given mass bin. The solid circles have been shifted by 0.05 in log for clarity. The open circles and error bars are the best-fitting KS mean and 68 per cent confidence level assuming a Gaussian parent distribution. The dashed lines connect the raw dispersion for each point, and the coloured solid lines are the best-fitting (KS test) dispersion. (See the electronic edition for colour version of the figure).

4.2 Shapes of haloes at higher redshifts

The use of M_* in equation (7) alludes to the evolution of the 〈s〉(M_vir) relation. After examining the 〈s〉(M_vir) relation at higher redshifts, we find that the relation between 〈s〉 and M_vir is successfully described by equation (7) with the appropriate M_*(z). The M_* for z= 1.0, 2.0 and 3.0 are 3.5 × 10¹¹, 1.8 × 10¹⁰ and 1.3 × 10⁹h⁻¹ M_⊙, respectively, for the simulations with Ω_b= 0.045. We present our results for various redshifts in Fig. 2 from the L120_0.9r, L80_0.9, L120_0.9 and L200_0.9 simulations. We have also included data points provided by Springel (private communication) in Fig. 2 for comparison. These data are from a more complete sample than the data presented in Springel et al. (2004) and are for shapes measured at 0.4R_vir.

Figure 2

〈s〉(M) for z= 0.0, 1.0, 2.0, 3.0. The binning is the same as in Fig. 1, but now for many different redshifts. The solid line is the power-law relation set out in equation (7). The L120_0.9 points are shifted by 0.05 in log for clarity. The Springel data agree quite well with our data and model for z= 0.0, 1.0, 2.0.

4.3 Dependence on σ₈

Of the parameters in the ΛCDM cosmological model, the parameter which is the least constrained and the most uncertain is the normalization of the fluctuation spectrum, usually specified by σ₈. Therefore, it is of interest to understand the dependence of the 〈s〉(M_vir) relation on σ₈. Since M_* is dependent on σ₈, the scaling with M_* in equation (7) may already be sufficient to account for the σ₈ dependence. As stated in Section 2, L80_0.75 and L80_0.9a were produced with the same Gaussian random field but different values for normalization. Therefore, the differences between the two simulations can only be a result of the different values for σ₈. As Fig. 3 illustrates, the two simulations do indeed produce different relations. We find that the M_* dependence in equation (7) is sufficient to describe the differences between simulations of different σ₈. One should expect this from the result of the previous subsection, that the redshift evolution was also well described by the M_* dependence. The values of M_* for z= 0.1 are 5.99 × 10¹²h⁻¹ M_⊙ for σ₈= 0.9 and 2.22 × 10¹²h⁻¹ M_⊙ for σ₈= 0.75. The values of M_* for σ₈= 0.75 at z= 1 and 2 are 1.09 × 10¹¹ and 4.57 × 10⁹h⁻¹ M_⊙, respectively. A simple fit to the redshift dependence of M_* in these cosmologies is log (M_*) =A−Blog(1 +z) −C[log(1 +z)]², with A (B, C) = 12.9 (2.68, 5.96) for σ₈= 0.9 and A(B, C) = 12.5 (2.94, 6.28) for σ₈= 0.75, and is accurate to within 1.6 and 3.1 per cent, respectively, for z≤ 3.

Figure 3

〈s〉 versus M(h⁻¹ M_⊙) with different values of σ₈. Different values of σ₈ predict different values for the 〈s〉 versus M relationship. Here, one can see that a universe with a lower σ₈ produces haloes which are more elongated, although the power-law relationship (equation 7) remains valid, as shown by the agreement between the points and the lines representing this prediction.

4.4 Mean–dispersion relationship

In the previous subsections, we used the mean belonging to the best KS test fit, assuming a Gaussian parent distribution, as an estimate of the true mean of axial ratios within a given mass bin. In this subsection, we examine the validity of this assumption and test whether the dispersion has the mass dependence suggested by JS. In Fig. 4, we present the distribution of s in the six bins from Fig. 1 for L120_0.9. In each of the plots, we have also included the KS best-fitting Gaussian, from which the mean was used to determine the best-fitting power law in equation (7). The error bars on the mean indicated in Fig. 1 are the 68 per cent confidence limits of the KS probability. The limits are determined by varying the mean of the parent distribution until the KS probability drops below 16 per cent for values greater and less than the best fit dispersion in each mass bin. The values from this analysis corresponding to the distributions in Fig. 4 can be found in Table 2. In Fig. 4, the lowest mass bin, which also contains the most haloes, is well fit by a Gaussian. This is seen in Table 2 not only by the best-fit KS probability, but also by the tight constraint of the confidence limits. The higher mass bins are consistent with having Gaussian parent distributions though the parent distributions' values for the mean and dispersion are not as well constrained. There is no indication of a structured tail to lower values of s, but Table 3 indicates that the distributions have negative skewness. This arises from a small number of haloes with very low values of s, which are always determined to be ongoing major mergers with very close cores.

Figure 4

The distribution of s in the L120_0.9b simulation in the mass ranges indicated in units of h⁻¹ M_⊙. The number of haloes in each bin is also indicated. The Gaussian fit shown for each graph is the best fit, based on a KS test analysis.

Table 2

KS best-fitting values.

Mass (h−1 M⊙)	〈s〉	σs	KS prob
4.8 × 1012	0.583 (+0.003 −0.003)	0.108 (+0.006 −0.005)	0.80
9.6 × 1012	0.554 (+0.006 −0.006)	0.110 (+0.005 −0.007)	0.61
1.92 × 1013	0.518 (+0.009 −0.004)	0.094 (+0.007 −0.013)	0.72
3.84 × 1013	0.519 (+0.014 −0.013)	0.108 (+0.016 −0.030)	0.95
7.68 × 1013	0.486 (+0.005 −0.010)	0.082 (+0.020 −0.015)	0.78
1.54 × 1014	0.467 (+0.012 −0.014)	0.073 (+0.050 −0.033)	0.93

Table 3

Weighted versus unweighted shapes.

Method	Mean mass ( h−1 M⊙)	〈s〉	σs	Kurtosis	Skewness
Weighted	5 × 1011	0.694	0.094	−0.029	−0.132
Unweighted	5 × 1011	0.686	0.097	−0.302	−0.053
Weighted	5 × 1012	0.614	0.111	0.253	−0.188
Unweighted	5 × 1012	0.603	0.118	0.605	−0.287
Weighted	5 × 1013	0.518	0.117	1.756	−0.532
Unweighted	5 × 1013	0.524	0.126	1.628	−0.592
Weighted	5 × 1014	0.455	0.114	2.066	−0.515
Unweighted	5 × 1014	0.472	0.133	1.428	0.150

JS found that the distribution of s within a given mass bin is Gaussian. They found no indication of a tail or any low values of s. This is most likely due to their treatment of haloes with multiple cores (see Section 6). Bullock (2002) found a large tail to low values of s using R=R_vir. After repeating our analysis at R=R_vir, we find the exact opposite. We find even less indication of a tail than in the distributions shown in Fig. 4. The difference is most likely due to the centres of haloes determined by the different halo finders used. Bailin & Steinmetz (2005) find a more subtle but significant tail to lower values of s. This is most likely just a side effect of combining all mass bins into one histogram. If the histogram were divided into bins over smaller ranges in mass, this tail would most likely be seen as a consequence of the combination of Gaussian distributions. These distributions would have the property that mass bins with a lower number of haloes also would have lower mean values, as in Fig. 4. Kasun & Evrard (2005) find that the distribution about the mean ‘is well fit by a Gaussian’ and contains no haloes with low values of s. We find the same in our spherical window analysis.

Now we turn our attention to the relationship between 〈s〉 and the dispersion for each bin. At z= 0, we find a dependence of the dispersion on 〈s〉 (top plot in Fig. 5). The relationship is steeper than that of JS, who determined that σ_s= 0.21 〈s〉 for the mass range they studied. However, one can see in the bottom plot of Fig. 5 that at z= 1 this relationship is no longer visible. This may be due to the fact that the number of haloes in each bin at z= 1 is much lower and therefore dominated by systematics. It could also be that a large enough range in mass is not probed at z= 1 to see the relationship. We are unable to draw a conclusion similar to JS. We therefore assume a constant value of σ_s= 0.1, which is consistent with our results at all redshifts.

Figure 5

〈s〉 versus σ_s. Top: we find some evidence for a dependence of the dispersion on 〈s〉 at z= 0.0 if we perform a weighted linear least-squares fit to the best values (black circles) of 〈s〉 and σ_〈s〉 (black short dashed line). It is slightly stepper than that of JS (green long dashed line). Also shown are the raw average and dispersion points (red triangles) Bottom: by z= 1.0 this relationship seems to have disappeared. Due to the lack of a clear relationship between 〈s〉 and σ_s, we favour a constant value of σ_s= 0.1, which is consistent with all redshifts and is roughly consistent with JS for the values of 〈s〉 probed.

4.5 Middle axis relationship

The largest-to-smallest axial ratio, s, does not uniquely determine an ellipsoidal shape. There is still the determination of the relationship of the middle axis (b) to the smallest (c) or largest (a) axis. We find in our analysis, as did JS, that the function P(p≡c/b|s) exhibits a nice symmetric behaviour. More commonly examined is the distribution of q≡b/a which can be trivially obtained from P(p|s). Figure 6 contains six histograms of p for different ranges in s. The curves are a fit proposed by JS,

9

with for s < s_min= 0.55 and for s≥s_min. It should also be noted that P(p | s) = 0 below ⁠. JS found that a cut-off of s_min= 0.5 best fit their data, but otherwise we find agreement with their results.

Figure 6

The distribution of p=c/b in given bins of s shows a very similar behaviour to that found in JS. The fit line is from equation (9), originally found in JS.

4.6 Radial dependence of shape

The ellipsoidal shape measured for a halo is also found to be dependent on the radius at which it is determined. There is a systematic dependence of shape on radius with more massive haloes having a steeper gradient in s with radius than lower mass haloes. In order to study the radial dependence of shape, haloes in the L80_0.9b simulation were examined at five different fractions of their virial radius (Fig. 7). For all halo mass bins, there is a tendency for haloes to be more spherical at larger radii, with more massive haloes having a steeper change in 〈s〉 with radius.

Figure 7

〈s〉 as a function of radius at z= 0. M_vir is in units of h⁻¹ M_⊙.

We also examined the value of 〈p〉 with radius and found no radial dependence. Therefore, a and b have the same radial dependence, and the largest axis a becomes relatively shorter with radius. We examine the relationship of the radial dependence of s with mass and combine it with equation (7) to find a shape–radius relationship,

10

with

11

In addition, we examined the shape at different radii within a spherical window. The shape of the halo did not change very much on average with radius. This is consistent with the result of Bailin & Steinmetz (2005) (see Section 6).

4.7 Triaxiality

Often, ellipsoids are described in terms of their triaxiality (prolate, oblate, or triaxial). One way of expressing the triaxiality of an ellipsoid is by using the triaxiality parameter Franx, Illingworth & de Zeeuw (1991):

12

An ellipsoid is considered oblate if 0 < T < 1/3, triaxial if 1/3 < T < 2/3 and prolate if 2/3 < T < 1. In Fig. 8, we divide up the haloes into the same mass bins as in Fig. 7 and analyse the triaxiality at R= 0.3R_vir and R_vir. We find that most haloes are prolate in shape with very few oblate haloes, even at R_vir. The deficit of haloes with T very close to 1 is not physical. Due to the iterative process we use to define shapes, if any two of the axes become degenerate (same length) the process has trouble converging. In most cases, it does converge but with a large systematic error. Some authors have suggested that haloes become oblate at large radii. We find only a small trend to less prolateness at large radii, but no evidence of a shift to oblate. Fig. 8 also shows that larger haloes, mainly those above M_*, are almost entirely prolate. Because we expect haloes with masses above M_* to be undergoing a higher rate of merging than haloes with masses below M_*, and because it has been shown that this merging happens along preferred directions (Knebe et al. 2004; Faltenbacher et al. 2005; Zentner et al. 2005), the prolateness is most likely due to merging. This is in support of the idea that halo merging is responsible for the distribution of shapes. Another quantity related to the merger history is the internal velocity of a halo, and therefore one would expect a relationship between the velocity structure of haloes and their shape. We examine this in the next section.

Figure 8

Triaxiality of haloes at z= 0 at R= 0.3R_vir (solid) and R_vir (dashed). Beginning with the top left histogram and moving right, then down, the triaxiality of haloes is divided into the same mass bins as in Fig. 7.

4.8 Alignment with velocity and angular momentum

We have shown that the shape of a halo is related to the mass and have seen some indications that the shape is due to merging. Merging is also related to the angular momentum of haloes (Vitvitska et al. 2002) and their velocity dispersion. In order for a collisionless system such as a DM halo to sustain its shape after merging, there must be an internal pressure provided by the velocity dispersion. If this is the case, one would expect the internal velocities to be correlated with the shape. In order to investigate this, we examine the alignment of the angular momentum and the velocity anisotropy of the halo with the shape. The angular momentum used here is calculated using the same particles found in the final ellipsoidal volume from our iterative method for determining the shape, although the results do not have a large dependence on which subset of particles within the halo is used. We find, as was pointed out by Warren et al. (1992), Tormen (1997), and subsequently seen by others, that the angular momentum is highly correlated with the smallest axis of the halo. In Fig. 9, we show the absolute cosine of the angle between the indicated axis and the angular momentum vector. If the orientations were random, the plot would be of a straight line at a value of 0.5. A peak at |cos θ| = 1 means that the axes are most often aligned and a peak at |cos θ| = 0 means that the axes are most often perpendicular to the angular momentum. As one can see, the smallest axis is most often aligned and the largest axis is most often at an angle of π/2 from the angular momentum. Although the angular momentum is aligned with the smallest axis as would be expected for an object which is rotationally supported, DM haloes are found not to be rotationally supported. Therefore, the significance of this alignment points not to a cause and effect relationship but to a shared origin. It has been shown in previous studies that the angular momentum of haloes is largely determined by the last major merger (Vitvitska et al. 2002), and that, at least during very active periods, merging (both minor and major) happens along preferred directions (Knebe et al. 2004; Faltenbacher et al. 2005; Zentner et al. 2005). It would seem, based on this, that the shapes and orientations of DM haloes, at least during active merging periods, can be attributed to directional merging.

Figure 9

Top: the probability distribution of the cosine of the angles between the largest, middle and smallest axis and the angular momentum vector. If the angular momentum were randomly oriented, the graph would be a flat line at a value of 0.5. Bottom: the probability distribution of the cosine of the angles between the shape axes and velocity anisotropy axes. The velocity anisotropy is highly correlated with the shape.

In order to determine whether haloes are relaxed and self-supporting, we examine the relation of the velocity anisotropy to the shapes of haloes. The velocity anisotropy is defined in the same way as the inertia tensor used to measure the shape,

13

We do not use a weighted version of this and we do not iterate, because neither of these makes much physical sense. We calculate the velocity anisotropy tensor again using the particles found within the ellipsoidal shell defined by the shape analysis. As with the angular momentum, the alignment is very insensitive to the particles used. We then determine the angle between the respective axes (i.e. a, b and c). In Fig. 9, we plot the distribution of absolute cosines between a_shape and a_vel, b_shape and b_vel, and c_shape and c_vel. From this, one can see that all three of the axes are highly correlated. The strength of the alignment between the velocity anisotropy tensor and the shape suggests that the shape is supported by internal velocities. But how does the shape relate more directly to the velocity anisotropy?

In Fig. 10, we show the triaxiality of the velocity anisotropy and the density. The velocity anisotropy on average is more spherical in shape than the density. This is the expected trend from the Jeans equation for an ellipsoidal distribution. The velocity anisotropy is directly related to the potential which has the same orientation as the shape but is more spherical due to the fact that potential is related to the spatial integral of the density. It would therefore seem that the mass dependence of shape cannot be explained by different relaxation times. This is also supported by the fact that haloes are more aspherical in the centre where the relaxation time would be shorter than at the virial radius.

Figure 10

Point density plot of triaxiality of the velocity anisotropy tensor and the shape for R= 0.3R_vir (blue/solid) and R_vir (red/long dash). The contours are for 10, 20, 30 and 40 haloes within 0.05 squared bins of triaxiality.

5 MERGER HISTORY AND SHAPES

So far, we have investigated the evolution of halo shapes in fixed mass bins as a function of redshift and for two different values of σ₈. Additional insight into the origin of shapes and their dispersion can be gained by tracing the evolution of individual haloes. In order to quantify the evolution or mass accretion history (MAH) of the haloes, we have constructed merger trees. For more information on the merger trees, please see Allgood (2005). From these merger trees, we determine the MAH for each halo at z= 0 by following the evolution of its most massive progenitors. Wechsler et al. (2002) showed that the MAH of a halo can usually be well fit by a single parameter model,

14

where M_o is the mass of the halo at z= 0 and a_c is the scalefactor at which the log slope of the MAH is 2. Although in Wechsler et al. (2002), they only allowed a_c to be a free parameter, we find that by also allowing M_o to be a free parameter we are able to better recover a_c for haloes which had experienced a recent major merger. Haloes with lower values of a_c form earlier, and as shown by Wechsler et al. (2002), have a higher concentration. By means of equation (14), we assign an a_c to every halo found at z= 0.

In Fig. 11, we plot s versus a_c for the haloes in the L120_0.9 simulation split into separate mass bins. We find that haloes which formed earlier are on average more spherical with a dispersion of 0.08 − 0.1 (see Fig. 12) for all mass bins. This implies that the scatter in the 〈s〉(M_vir) relation cannot be completely attributed to the different values of a_c for that particular mass bin. However, the dependence on a_c is less for the higher mass haloes, and this would explain the mean–dispersion relationship explored in Section 4.4. Higher mass haloes are found to have a smaller dispersion than lower mass haloes at z= 0. This can also be seen in Fig. 12. It is very likely that the residual scatter is due, at least in part, to the pattern of infall. Since s is derived from an inherently three-dimensional quantity, namely the inertial tensor, a one-dimensional parametrization may not be sufficient to capture all of the physics involved. A more careful study of infall is needed to explain the dispersion in s.

Figure 11

〈s〉 versus characteristic formation epoch for different mass bins (mass quoted in units of 10¹²h⁻¹ M_⊙). Only bins that contain at least 10 haloes are shown (square points). There is a distinct trend of shape with a_c for the lower mass bins. At higher mass, there is still a trend but it is uncertain how strong the trend is due to the lower number statistics. Solid black line is a linear fit to the points, and dashed line is the 1σ scatter about the points.

Figure 12

Distributions of s for a given range in a_c and mass. The solid histogram is for the mass range: 3.2 × 10¹²h⁻¹ M_⊙ < M_vir < 6.4 × 10¹²h⁻¹ M_⊙ and the dotted histogram is for a mass range of 1.28 × 10¹³h⁻¹ M_⊙ < M_vir < 5.12 × 10¹³h⁻¹ M_⊙. The dispersion in s in a given mass bin can be explained in part by the different MAHs.

The above investigation makes clear that the dispersion of halo shapes cannot be explained by appealing to a single parameter description of the MAHs. However, an average evolutionary pattern for haloes which is dependent on both a_c and mass is seen. Fig. 13 displays the evolution of 〈s〉 (sorted by a_c) with scalefactor a for a particular mass bin at z= 0. Haloes that formed early (lower a_c) are more spherical today as was pointed out above. Moreover, they become spherical more rapidly (indicated by the increasing slopes for haloes of low a_c), although the transformation rate towards spherical shapes seems to slow for all values of a_c with increasing expansion factor. In Fig. 13, the results for the lowest mass bin (3.2 × 10¹²h⁻¹ M_⊙ < M_vir < 6.4 × 10¹²h⁻¹ M_⊙) are shown. Apart from a systematic shift to lower values of s, the corresponding plots for higher mass bins look very similar.

Figure 13

The evolution of haloes with different values of a_c in the mass bin 3.2 × 10¹²h⁻¹ M_⊙ < M_vir(z= 0) < 6.4 × 10¹²h⁻¹ M_⊙. Haloes become more spherical after a short period after a_c. The haloes which form earlier become spherical more rapidly. Log binning was chosen to even out the number of haloes in each bin.

Haloes that have early formation times (low a_c) at a fixed mass today have typically accreted more mass since a_c than haloes with higher values of a_c. The rapid transformation towards spherical shapes for early forming haloes implies either that lower mass haloes become spherical more rapidly after a_c, probably due to shorter dynamical times, or that mass accretion after a_c is more spherical, therefore causing the halo to become more spherical as well. By examining other mass bins, we find that haloes of different masses today, but with comparable values for a_c, show the same rate of change in s, but with different initial values of s. This finding suggests that the rate at which a halo becomes spherical depends on its a_c rather than on its mass. We find that we can approximate the dependence of s on the expansion factor a for a > a_c+ 0.1 by a simple power law

15

where ν has to be fitted for the particular halo. In Fig. 14, we display the values of ν versus a_c determined by χ² fitting for the L80_0.9b simulation. The L80_0.9b was divided into logarithmic bins of mass and a_c. The average MAH for each bin was fit by equation (15) using χ² minimization. All bins containing at least 20 haloes were used to determine the fit. We find a tight correlation between the ν and a_c which can be approximated by

16

Figure 14

The rate of change exponent ν (see equation 15) versus expansion factor at halo formation a_c. The solid line displays the fitting formula given by equation (16). Only bins containing at lease 20 haloes are displayed. The a_c value of each point is the average value for that respective bin and the error bars represent the variance determined by the χ² fitting.

This fit is represented by the solid line in Fig. 14. The remarkable success of equation (15) to fit the data supports the idea that the transformation from aspherical to spherical halo shapes is driven by mass accretion becoming more spherical after a_c. The physical reason for the observed behaviour merits further investigation.

6 COMPARISON WITH PREVIOUS DETERMINATIONS OF HALO SHAPE

In the previous sections, we have explored many aspects of halo shapes. Central to this discussion has been 〈s〉(M_vir). This relationship has been examined by many recent studies, all of which seem to determine different relationships. In this section, we address these discrepancies.

First, an examination of the difference in the inferred shape from the use of the weighted versus the unweighted inertia tensor is needed. Most recent authors prefer the weighted (or reduced) inertia tensor (equation 5) which is the method we have chosen to use. The motivation for the use of the weighted inertia tensor, ⁠, is due to the bias present in the unweighted method to particles at larger radii. By weighting the contribution from each particle in the sum by the distance to the particle squared, is less sensitive to large substructure in the outer regions of the analysis volume. To test the difference between the methods, we examined a sample of haloes from the L120_0.9 simulation using the iterative method with both versions of inertia tensor (Fig. 15). Both iterative methods give similar results for the mean quantities (inset in Fig. 15) as a function of mass with individual haloes differing by Δs≤ 0.15. The detailed distributions are different and have some interesting features (Table 3). We find that haloes which have lower s values for the weighted method over the unweighted method have larger substructure near the centre and haloes which have lower s values for the unweighted method have larger substructure near the outer edge of the analysis volume. At very small axial ratios, the unweighted method seems to always give larger axial ratios. As we have shown, haloes with s≲ 0.3 are late forming and are strongly contaminated by substructure. This leads to the unweighted method giving a lower value of s indicated by the high skewness shown in Table 3. For the two well-resolved mass ranges (10¹²h⁻¹ M_⊙ < M_vir < 10¹³h⁻¹ M_⊙ and 10¹³h⁻¹ M_⊙ < M_vir < 10¹⁴h⁻¹ M_⊙), the unweighted method has a larger negative skewness. Note that all mass bins except for the last unweighted bin have negative skewness (this was discussed in Section 4), but the most massive bin suffers from low statistics. The distribution of σ_s values is always broader for the unweighted method.

Figure 15

Axial ratios for haloes in a cosmological simulation are divided into mass bins, and the shapes are calculated using the weighted (s_w) and non-weighted (s_n) iterative inertia tensor methods. The two methods agree within ∼10 per cent and give the same value when averaged over a given mass bin (inset graph).

The analyses in the literature differ by more than just the form of the inertia tensor used. In order to compare our results to a selected number of previous results (Fig. 16), we have repeated the shape analysis using the methods described in the corresponding papers. The previous work which is most similar to the current work is that of Springel et al. (2004). Our findings are very similar to the Springel et al. (2004) results except we found that 〈s〉(M_vir) has a slightly higher normalization. Through private communication with Volker Springel, we were provided with an updated set of data points which come from a more complete sample and are in much better agreement with our results (open green squares in Fig. 16). Not only do our results agree at z= 0, but also at higher redshift (see Fig. 2).

Figure 16

Comparison of 〈s〉(M_vir) relation with previous studies. We attempt to reconcile the differences between our results and those of other authors. We present the results of a shape analysis of the L80_0.9b, L120_0.90r and L200_0.9 simulations using the iterative inertia methods at R= 0.3R_vir (black, pink, cyan) and kolonon-iterative spherical window analysis at R_vir (red,green,violet). In addition, we present the results of a shape analysis of the L80_0.9b and L120_0.90r using the iterative method at R_vir (blue). The black line is our proposed fit from equation (7), and this should be compared to the results of Springel (private communication) (green open squares) and JS (orange dot–dashed). The blue line is a fit to the blue points, which should be compared to the Bullock (2002) line (violet long dash). Finally, the red line is a renormalized version of the Kasun & Evrard (2005) fit which should be compared to their fit (brown small dash). The thin black dot–dashed line at 〈s〉∼ 0.7 is the spherical shell fit of Bailin & Steinmetz (2005). The bold portions of the lines indicate the mass ranges where the fit was compared to simulated data by the respective authors.

JS studied 12 high-resolution clusters with N∼ 10⁶ particles and five cosmological simulations with N= 512³ particles in a 100 h⁻¹ Mpc box with both an standard CDM (SCDM) and ΛCDM cosmology. The simulations were performed with a P³M code with fixed time-stepping and a spatial resolution of 10–20 h⁻¹ kpc. They used a friends-of-friends (FOF) halo finder and analysed the shapes of the high-resolution clusters in isodensity shells as a function of radius, finding that the haloes are more spherical at larger radii. After determining the relationship of shape with radius, they developed a generalized ellipsoidal NFW density profile (Fig. 17 shows that our method gives similar results). They applied this generalized fitting routine to the cosmological simulations and determined generalized NFW parameters and shape statistics. The shapes were determined using an isodensity shell at an overdensity of 2500ρ_c (where ρ_c is the critical density) which corresponds roughly to R= 0.3R_vir. The mass range analysed only covered one order of magnitude in mass (2.1 × 10¹³h⁻¹ M_⊙≤M_vir≲ 1 × 10¹⁴h⁻¹ M_⊙). They found a result very similar to ours, ⁠, and a dispersion which is well-fitted by a Gaussian distribution with ⁠. We do not find any evidence for a steepening of the exponent with redshift as they do.

Figure 17

Comparison of an isodensity shell (red) and a tensor ellipsoid (blue). Particles are selected by the JS (red) and inertia tensor (blue) methods and projected on the x–y plane of the simulation box. The shortest/longest axis ratio is s= 0.49 (0.48) for the isodensity shell (tensor ellipsoid). The semimajor axis of the isodensity shell 0.23R_vir, consistent with JS for 2500ρ_c isodensity.

The disagreement between our findings and those of JS regarding the scaling of 〈s〉 with mass is due to the procedures used. In the JS analysis, they determine the shape of an isodensity shell at 2500ρ_c± 3 per cent, completely ignoring the interior of the shell. JS analysed haloes with masses greater than 6.2 × 10¹²h⁻¹ M_⊙ which tend to be dynamically young and often have double cores. This can affect the shape a lot, but their analysis would not pick this up, due to the neglect of the central region. In our iterative inertia tensor analysis, we include the centres. In order to confirm that the difference is truly due to the shell versus the solid ellipsoid, we analysed haloes from the L200_0.9 simulation in the mass range 1 − 4 × 10¹⁴h⁻¹ M_⊙ using the technique presented in JS. We examined isodensity shells at 2500ρ_c with a thickness of ±30 per cent, instead of the ±3 per cent used by JS. We needed to examine thicker shells in order to obtain enough particles to do the analysis because the L200_0.9 has less mass resolution than the simulations studied by JS. We found that the inertia tensor method gives 〈s〉_tensor= 0.485 ± 0.008 and σ_s= 0.091 ± 0.006, and the JS method gives 〈s〉_JS= 0.515 ± 0.008 with the same scatter for the same haloes. The difference is due to the fact that s_JS is systematically larger at low s_tensor. This pattern is born out by a quantitative analysis. When we split the sample at s_tensor= 0.45 (roughly where agreement begins), we get that above s_tensor= 0.45 the samples agree quite well with 〈s〉_tensor= 0.571 ± 0.007 versus 〈s〉_JS= 0.576 ± 0.011. Whereas below s_tensor= 0.45, we get 〈s〉_tensor= 0.424 ± 0.007 versus 〈s〉_JS= 0.472 ± 0.008. The difference at the low end is due to the missing of the dynamically active cores by JS. If JS had extended their analysis to lower mass haloes where multiple cores are not as common their determination of 〈s〉 would converge with ours. Because X-ray observers normally do not choose to only analyse the outer shells of clusters due to the fact that the X-ray observations get noisier with the distance from the centre and because optical observers may not see the multiple cores when analysing cluster member velocities, we prefer the method which includes the effect of the multiple peaks. In Paper II, we show that using our method with some additional assumptions one can account for the observed X-ray ellipticity measurements.

Bullock (2002) analysed the shapes of haloes in a ΛCDM simulation with σ₈= 1.0 at three different redshifts (z= 0.0, 1.0, 3.0). The simulations were performed using the ART code in a 60 h⁻¹ Mpc box with 256³ particles and spatial resolution of 1.8 h⁻¹ kpc. The analysis of shape was done using the weighted inertia tensor in a spherical window with R=R_vir. The axial ratios were determined iteratively until convergence was obtained using a similar criterion as we have used, but with the window remaining spherical. The use of the weighted inertia tensor and iterative axial ratio determination seemed to almost eliminate the effect of using a spherical window (discussed below). Bullock (2002) found that 〈s〉(M_vir) ≃ 0.7 (M_vir/10¹²h⁻¹ M_⊙)^−0.05 (1 +z)^−0.2 was a good fit to the simulation. The empirical scaling of Bullock (2002) is similar to what we find, but the power law is steeper. This can be attributed to the lower resolution and possibly the use of a spherical window. Bullock's higher normalization is due to the higher σ₈.

Kasun & Evrard (2005) determined the shapes of cluster haloes (M₂₀₀ > 3 × 10¹⁴h⁻¹ M_⊙) in the Hubble Volume simulation. They calculated the axial ratios using the unweighted inertia tensor in a spherical window at R₂₀₀, the radius of the sphere within which the mean density is 200ρ_c(z), with ρ_c(z) being the critical density at redshift z. They determined a relationship of 〈s〉(M_vir) = 0.631[1 − 0.023 ln (M_vir/10¹⁵h⁻¹ M_⊙)](1 +z)^−0.086. We compare our analysis with theirs by performing the same spherical analysis at R=R_vir, which is slightly larger than R₂₀₀. We find that in our largest box simulation (L200_0.9) where we have good statistics on cluster mass haloes we find good agreement. In examining the other two simulations for lower mass haloes, we are unable to recover the extrapolation of the Kasun & Evrard (2005) relationship. In fact, we see a transition from the Kasun & Evrard (2005) relationship to the relationship of Bailin & Steinmetz (2005) (discussed below). We also find that the mean shape relationship has almost no dependence on radius when using a spherical window function.

Bailin & Steinmetz (2005) analyse the shapes of haloes at different radii in spherical shells in the mass range of 10¹¹h⁻¹ M_⊙ < M≲ 5 × 10¹³h⁻¹ M_⊙. After determining the axial ratios, they then apply an empirical correction of to correct for the use of a spherical window. They find that all haloes have an axial ratio of 〈s〉∼ 0.63 at R= 0.4R_vir with the scaling applied, which implies that they measure 〈s〉∼ 0.766 in the spherical window. This result is in very good agreement with our spherical analysis (green and red data points in Fig. 16). However, we do not find that haloes of different masses have the same mean axial ratio. There seems to be some evidence that the 〈s〉(M_vir) relationship flattens out below M_*, but it is definitely not constant with radius. Simulations with even higher mass resolution are needed to investigate for the possibility of flattening below M_*. For an extra check, we also analysed the haloes in a spherical shell between 0.25R_vir and 0.4R_vir, and measure a roughly flat value for all haloes of s= 0.77. The disagreement about the 〈s〉(M) relationship most likely lies in the determination of the empirical spherical window correction. The correction was determined using Monte Carlo (MC) haloes with no substructure, but we find that substructure plays a role in the determined shape of the halo.

7 COMPARISON WITH OBSERVATIONS

Since all of the differences between the shape statistics extracted from pure collisionless simulations by various authors can be reconciled by considering the different methods used to determine shapes, a comparison between observations and simulations is in order. Much of the attention halo shapes have received lately is due to the recent estimates of the shape of the MW's host halo. Most estimates find the MW's host halo to have an oblate shape with s≥ 0.8. This is in contrast with s≈ 0.6 ± 0.1 for 10¹²h⁻¹ M_⊙ haloes found in pure collisionless simulations, though there is some evidence that the haloes become more spherical when baryonic cooling is included (Kazantzidis et al. 2004; Bailin et al. 2005) and that some become oblate. The presence of gas cooling will invariably make the haloes more spherical, but the extent of the effect is not yet fully understood. Recently, there have, however, been studies of the M giants in the leading edge of the Sagittarius dwarf stream (Helmi 2004; Law et al. 2005), which concluded that the best-fitting shape of the host halo is a prolate ellipsoid with s= 0.6.

Another way of measuring the shape of DM haloes is through weak lensing. Hoekstra et al. (2004) and Mandelbaum et al. (2005) performed studies of galaxy–galaxy weak lensing using the Red-Sequence Cluster Survey and the Sloan Digital Sky Survey, respectively. Hoekstra et al. (2004) determine the average shapes of haloes by measuring the orientation of the galaxies, then stacking the galaxy images with the orientations aligned, and finally measuring the shear field around the stacked image. This measurement of the shear provides a rough estimate of the dark matter halo shapes at z≈ 0.33. They found an average projected ellipticity of 〈ε〉≡〈1 −q_2D〉= 0.20^+0.04_−0.05, corresponding to s= 0.66^+0.07_−0.06, for haloes with an average mass of 8 × 10¹¹h⁻¹ M_⊙. The Hoekstra et al. (2004) determination was hindered by the fact that the galaxies were stacked together regardless of morphological type. One would reasonably assume that morphological type is related to merger history and possibly orientation with the host halo which will in turn affect the measured halo shapes as we have shown.

In the Mandelbaum et al. (2005) analysis of SDSS, galaxies colour was used as a proxy for morphological type. Mandelbaum et al. (2005) studied 2 million lens galaxies with r > 19 and 31 million source galaxies dividing the lenses into bins of colour and luminosity. They find a suggestion that isophotal ellipticities of spiral (blue) galaxies may be anti-aligned with the halo ellipticities at the 2σ level and a suggestion that elliptical (red galaxy) host halo ellipticities are luminosity dependent. Since we and others have shown that halo angular momentum is highly aligned with the smallest axis, this finding would suggest that the angular momentum of spiral galaxies is not aligned with the angular momentum of the their host haloes. Recent theoretical work by Bailin et al. (2005) seems to support this idea as well, although they only study eight spiral galaxy simulations. If this were the case, one would assume that in the Hoekstra et al. (2004) analysis the measured signal would be diminished by this. It would, however, not be completely nullified because the Hoekstra et al. (2004) sample is dominated by elliptical galaxies. Indeed, Mandelbaum et al. (2005) show by combining the appropriate luminosity and colour bins that their findings are in agreement with those of Hoekstra et al. (2004). In Fig. 18, we have plotted the distribution of projected axial ratios for the L80_0.9a and L80_0.75 simulations for 1000 random lines of sight through the box for galaxy mass haloes. In projection the differences between the simulations become small hindering any sort of determination of σ₈ via lensing studies. The peaks of the distributions are both broadly in agreement with the findings of Hoekstra et al. (2004) and Mandelbaum et al. (2005). Further study of the galaxy/host halo alignment relationship and galaxy morphology/halo merger history relationship is in order to better understand and predict the galaxy–galaxy lensing measurements.

Figure 18

Projected ellipticity of galaxy mass haloes at z= 0.33 for σ₈= 0.9 (black solid) and σ₈= 0.75 (red dashed).

8 SUMMARY AND DISCUSSION

In this study, we examined the shapes of haloes characterized by ellipsoidal fits to the weighted inertia tensor, using six different simulations in order to adequately cover three orders of magnitude in mass, from galaxy to cluster scales. We examined the dependence of the shape parameters on halo mass and radius. We also explored the evolution of the shape and its connection to infall and velocity anisotropy. Our main conclusions can be summarized as follows.

We find that the mean largest-to-smallest axial ratio s=c/a at radius 0.3R_vir is well described by

17

• The distribution of s in each mass bin is well-fitted by a Gaussian with σ= 0.1. The relation found here is steeper than that of JS at z= 0, thus predicting less spherical cluster mass haloes and more spherical galaxy mass haloes.

• Within fixed mass bins the redshift dependence of 〈s〉(M_vir) is well characterized by the evolution of M_*, unlike the findings of JS who predict a much steeper relation of ∝M^−0.07_vir at high redshifts. We find that equation (17) works well for different values of σ₈ (a variation of σ₈ results in a variation of M_* which appears as a normalization parameter in the 〈s〉(M_vir) relation). Also worth noting is that at higher redshift the possible broken power-law behaviour disappears, as would be expected if it were truly due to M_*, because by z= 1 M_* is below our mass resolution.

• We find that the mean shape relation becomes shallower and more spherical at increased radius. We also find that higher mass haloes have a steeper relationship between shape and radius than smaller mass haloes. Since cluster-sized haloes are on average younger than galaxy-sized haloes, we are comparing dynamically different objects. The presence of an increased amount of massive substructure near the centre of dynamically young objects may be the reason for the steeper relation of shape with radius for cluster mass haloes than galaxy mass haloes.

• We find that the 〈s〉(M_vir) at z= 0 for galaxy mass haloes is ∼0.6 with a dispersion of 0.1. This result is in good agreement with only one estimate for the axial ratio of the MW halo. Helmi (2004) claims that a study of the M giants in the leading edge of the stream tidally stripped from the Sagittarius dwarf galaxy leads to a best-fitting prolate halo with s= 0.6. After analysis of the same data, Law et al. (2005) confirm this finding. Other studies (Ibata et al. 2001; Majewski et al. 2003; Martínez-Delgado et al. 2004) which examined different aspects of these streams concluded that the MW halo is oblate and nearly spherical with q≳ 0.8. If the shape of the MW halo is truly this spherical, it is either at least 2σ more spherical than the median, or else baryonic cooling has had a large effect on the shape of the dark matter halo (see e.g. Kazantzidis et al. 2004).

• Describing halo shapes by the triaxiality parameter T introduced by Franx et al. (1991), we find that the majority of haloes are prolate with the fraction of haloes being prolate increasing for haloes with M_vir >M_*. Since halo shapes are closely connected to their internal velocity structure, we compute the angular momentum and the velocity anisotropy tensor and relate them to both the orientation of the halo and the triaxiality. In agreement with previous studies, we find that the angular momentum is highly correlated with the smallest axis of the halo and that the principal axes of the velocity anisotropy tensor tend to be highly aligned with the principal axes of the halo. The strong alignment of all three axes of the two tensors is remarkable since the velocity tensor tends to be more spherical, thus the determination of its axes might be degenerate which would disturb the correlation with the spatial axes. If the accretion of matter determines the velocity tensor, the tight correlation between velocity tensor and density shape argues for a determination of the halo shape by directional accretion.

• Finally, we examine the evolution of shapes by following the merger trees of the individual haloes. We find that the different MAH of haloes cannot fully explain the observed dispersion about the mean s within fixed mass bins. It is likely that an analysis of the three-dimensional accretion is essential for the explanation of the dispersion at a fixed value of mass and a_c. However, haloes with earlier formation times (lower a_c) tend to be more spherical at z= 0. Furthermore, there is a pattern of haloes becoming spherical at a more rapid rate for haloes that formed earlier, and this rate appears to be independent of the final mass. The evolution of the shape for a >a_c+ 0.1 is well described by

• 

  18

• where ν= 1.74 ×a^−0.3_c. We detect a definite trend for the transformation from highly aspherical to more spherical halo shapes after a_c. The change of s seems to be less dependent on the total halo mass but strongly influenced by the relative mass increase since a_c which suggests that haloes are becoming more spherical with time due to a change in the accretion pattern after a_c from a directional to a more spherical mode.

We thank Anatoly Klypin for running some of the simulations which are used in this work and for help with the others. We also thank Volker Springel and Eric Hayashi for helpful private communications. The L120_0.9r and L80_0.9b were run on the Columbia machine at NASA Ames. The L80_0.9a, L80_0.75 and L120_0.9 were run on Seaborg at NERSC. BA and JRP were supported by a NASA grant (NAG5-12326) and a National Science Foundation (NSF) grant (AST-0205944). AVK was supported by the NSF under grants AST-0206216 and AST-0239759, by NASA through grant NAG5-13274, and by the Kavli Institute for Cosmological Physics at the University of Chicago. RHW was supported by NASA through Hubble Fellowship HF-01168.01-A awarded by Space Telescope Science Institute. JSB was supported by NSF grant AST-0507816 and by startup funds at UC Irvine. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

1

To calculate the density, we use the publicly available code smooth: http://www-hpcc.astro.washington.edu/tools/tools.html

Appendix

APPENDIX A: RESOLUTION TESTS

A potential source of systematic error in the determination of halo shapes is the limited number of particles involved for low-mass haloes. In Fig. A1, we show the result of a MC test to examine this. Our dark matter haloes are adequately described by an elliptical NFW density profile independent of cosmological model. Therefore, in the figure we show the result of applying the reduced tensor method of equation (5) to MC haloes built to have an ellipsoidal NFW profile. The axial ratios of the MC haloes are drawn from Gaussian distributions of mean 〈s〉= 0.7 and 〈q〉= 0.85, and dispersion σ= 0.1. Approximately 450 haloes were generated, each having ∼1000 particles, in order to have a catalogue comparable to the sample of haloes in the L80_0.9 box with mass in the range M_vir= 10^11.3− 10^11.7h⁻¹ M_⊙. The scatter plot shows that individual values of s determined at 0.3R_vir which contains roughly 300 of the 1000 particles in the halo. The recovered shape can be in error by up to ∼0.1. However, the scatter and mean of the distribution are very well determined by the inertia tensor method. The inset shows the histogram for the input values of the MC haloes (solid line), and the histogram for the output values (i.e. determined by the tensor method; dotted line). We conclude that the tensor method underestimates s by only 0.03 for haloes containing ∼1000 particles. For haloes with more than 2500 particles, the error falls to 0.01.

Figure A1

Results of applying our shape determination procedure at 0.3R_vir to 450 MC haloes produced with determined axial ratios. We found that the error in the recovered value of s could be as large as ∼0.1.

One shortcoming of the above test is that we are attempting to recover shapes from smooth haloes. In reality haloes have substructure, and the substructure plays a role in the shape determination of the halo. The presence of dense lumps close to the core will bias 〈s〉 to lower values relative to values determined from isodensity shells. In order to test our ability to recover the axial ratios in a cosmological simulation, we determine the mean shape of haloes down to very low masses for all of our simulations and then compare them to one another. The main difference between this test and the previous one is that the haloes in the cosmological simulations contain substructure, but we do not know a priori what the distribution of shapes should be. In Fig. A2, we have plotted the shapes of haloes down to very low particle numbers. There appears to be two resolution effects at work here. The first effect is an extension of the result we found in the MC test above. At small particle number, the recovered shape becomes very aspherical and all of the simulations turn over at the same number of particles. For all of the simulations, the turnover is detected at ∼3000 particles within R_vir (which is ∼1000 particles within 0.3R_vir, where we are determining the shape). This is consistent with MC test which were performed with comparable particle numbers. At higher particle numbers, there seems to be another effect driven by the particle number and not the mass. Haloes containing n_p < 7000 particles for any given simulation show a trend of becoming more spherical on average than the simulations of higher resolution for the same mass. In Fig. A3, we show the determined value for s for the most massive haloes shared by the L120_0.9 simulation and the resimulated subregion, L120_0.9r. We find that haloes shapes between the two simulations can differ by as much as 0.15 in s. This is not unexpected because in the resimulation process the identical haloes are captured at slightly different times. However, the main point here is that the average value of s is lower in the higher resolution simulation for the same haloes by ∼0.05. This behaviour can be explained by the more obstinate survival of substructure in the high-resolution simulations.

Figure A2

〈s〉 versus mass. This plot is a replica of Fig. 1 except we show mass bins below the determined resolution limit.

Figure A3

A direct comparison of the axial ratios of the most massive haloes in the L102_0.9r simulation to the corresponding haloes in the L120_0.9 simulation.

Author notes

†

Hubble Fellow, Enrico Fermi Fellow.