# The Graininess of Dark Matter Haloes

###### Abstract

We use the recently completed one billion particle Via Lactea II CDM simulation to investigate local properties like density, mean velocity, velocity dispersion, anisotropy, orientation and shape of the velocity dispersion ellipsoid, as well as structure in velocity space of dark matter haloes. We show that at the same radial distance from the halo centre, these properties can deviate by orders of magnitude from the canonical, spherically averaged values, a variation that can only be partly explained by triaxiality and the presence of subhaloes. The mass density appears smooth in the central relaxed regions but spans four orders of magnitude in the outskirts, both because of the presence of subhaloes as well as of underdense regions and holes in the matter distribution. In the inner regions the local velocity dispersion ellipsoid is aligned with the shape ellipsoid of the halo. This is not true in the outer parts where the orientation becomes more isotropic. The clumpy structure in local velocity space of the outer halo can not be well described by a smooth multivariate normal distribution. Via Lactea II also shows the presence of cold streams made visible by their high 6D phase space density. Generally, the structure of dark matter haloes shows a high degree of graininess in phase space that cannot be described by a smooth distribution function.

###### keywords:

Galaxies: haloes — Galaxies: structure — Galaxies: kinematics and dynamics — dark matter — methods: -body simulations — methods: numerical^{†}

^{†}pagerange: The Graininess of Dark Matter Haloes–B

^{†}

^{†}pubyear: 2008

## 1 Introduction

Spherical averaging is a commonly used method to describe the characteristics of dark matter haloes that form by gravitational collapse in a cosmological environment. For example one of the basic characteristics of a dark matter halo is its spherically-averaged density profile which can be well described by a smooth function within resolved scales (e.g. Navarro et al., 1996; Moore et al., 1998; Diemand et al., 2005; Diemand et al., 2007; Diemand et al., 2008; Stadel et al., 2008). Other examples include the velocity dispersion and anisotropy profiles.

Consider a set of observers at a given distance from the halo centre and capable of measuring local properties of the halo (e.g., density, velocity dispersion, velocity anisotropy). In a smooth, spherically symmetric halo, these measurements would yield identical values, up to statistical fluctuations. However, it is well known that dark matter haloes that form in pure dark matter simulations, which neglect possible baryonic effects in the centre, are in general triaxial: close to prolate in the central part and becoming rounder in the outskirts of the halo (e.g. Katz, 1991; Dubinski & Carlberg, 1991; Allgood et al., 2006; Bett et al., 2007; Kuhlen et al., 2007). Such a shape variation obviously leads to a significant variance in local properties compared to the spherically averaged value at a given radius (Knebe & Wießner, 2006). In addition, haloes have a high level of subhaloes which also effect results of local measurements. But local deviations from a smooth model go beyond these shape and subhalo driven variations, as we demonstrate in this paper. For example, overdensities are expected due to the presence of subhaloes, but we also find underdense regions, which are unexpected in a smooth triaxial background halo with subhaloes (see section 3.2).

One aim of this paper is to quantify the degree of variation between local and spherically averaged values for key quantities such as the density and velocity dispersion. A better understanding of the local variations of these properties is essential in order to obtain a better description of the phase space structure of dark matter haloes and their formation process.

Such local variations can have profound consequences. For example, in dynamical models of dark matter haloes and galaxies it is often assumed that most properties are just a function of radius (e.g. Dehnen & McLaughlin, 2005; Baes & van Hese, 2007; Van Hese et al., 2008). From the work presented here it becomes clear that this is not a valid assumption for the anisotropy (see section 3.5) and a lot of information is smeared out by spherically averaging. Models based on simple analytical forms of the distribution function, which are also used for generating -body realisations of dark matter haloes and galaxies (e.g. Kazantzidis et al., 2004; Zemp et al., 2008; Zhang & Magorrian, 2008), often have the property that the local bulk motion is zero and the velocity distribution function is similar to a smooth multivariate normal distribution (if it is not even deliberately set to a multivariate normal distribution as in Hernquist (1993)). Both of these characteristics are not valid in the outskirts of dark matter haloes that form in a hierarchical cosmological context (see sections 3.3 and 3.7). As we see, most of these effects are not taken into account in current dynamical models of dark matter haloes and galaxies and they will be a challenge to be modelled with more realistic distribution functions.

In addition, it is also important to quantify these variations locally at the Earth’s position within the Galaxy in order to interpret results from present and future direct dark matter detection experiments on Earth
(e.g. Moore
et al., 2001; Stiff
et al., 2001; Helmi
et al., 2002; Stiff &
Widrow, 2003; Kamionkowski & Koushiappas, 2008; Vogelsberger et al., 2008; Fairbairn &
Schwetz, 2008). While small subhalos contribute significantly to the indirect detection signal (e.g. Diemand
et al., 2006; Kuhlen
et al., 2008) and infall caustics have no significant effect (Diemand &
Kuhlen, 2008), it is worthwhile also to investigate whether additional lumpiness from tidal debris might enhance the annihilation rate (see section 3.2)^{1}^{1}1After submitting this work an independent, analytic study of this question (Afshordi et al., 2008) has appeared. Our simulation results agree that such an enhancement is rather small (see section 3.2)..

Also recent results from surveys like SDSS and RAVE need a better understanding of the structure of the outer halo. There is a lot of known structure in the Galactic stellar halo like for example the Sagittarius stream, the Monoceros stream or the Virgo overdensity (e.g. Martínez-Delgado et al., 2007; Cole et al., 2008; Casetti-Dinescu et al., 2008). Obviously, these features are baryonic, but stars are collisionless too, and it’s probable that these features originate from infalling dark matter dominated objects.

In general, a more detailed picture of the phase space structure of dark matter haloes is needed in order to understand and model their properties that are set by the hierarchical formation process. This work is a further important step towards that aim.

Here we present a study of the distribution of these properties as a function of galactocentric distance of a Milky Way size dark matter halo with data from the Via Lactea II project (Diemand et al., 2008). In section 2 we give a summary of the simulation data used in this study and present the results in section 3. We then discuss our results and present the conclusions in section 4. In appendix A we give a comparison of the results with lower resolution simulations and different definitions of locality.

## 2 Simulation and data

8 | 25 | 50 | 100 | 200 | 400 | ||

0.5 | 1.004 | 1.612 | 2.930 | 5.728 | 10.69 | ||

-0.2116 | -1.075 | -0.2053 | -4.207 | -2.975 | 9.285 | ||

0.7794 | 2.946 | -8.476 | -0.1219 | 0.6685 | 9.969 | ||

-0.3097 | -0.8743 | 0.7920 | -4.188 | 4.331 | -12.52 | ||

239.5 | 246.2 | 219.1 | 187.8 | 147.5 | 122.2 | ||

144.1 | 158.5 | 145.3 | 125.9 | 98.99 | 73.03 | ||

143.1 | 145.4 | 129.4 | 109.8 | 83.36 | 76.23 | ||

127.0 | 119.9 | 100.9 | 85.82 | 70.86 | 61.62 | ||

[1] | 0.1188 | 0.2938 | 0.3620 | 0.3872 | 0.3893 | 0.09919 | |

[] | 0.1480 | 4.633 | 3.164 | 4.893 | 16.18 | 37.31 | |

[] | 1.121 | 5.596 | 3.577 | 13.06 | 16.53 | 48.12 | |

[] | 1.117 | 3.166 | 3.532 | 12.80 | 5.159 | 31.24 | |

[1] | 0.6308 | 0.4893 | 0.5066 | 0.3894 | 0.6017 | 0.6091 | |

5.681 | 14.03 | 16.71 | 23.24 | 26.47 | 30.57 | ||

6.077 | 11.45 | 19.44 | 28.22 | 27.08 | 36.59 | ||

4.474 | 11.54 | 13.96 | 19.52 | 20.91 | 29.13 | ||

6.745 | 9.491 | 8.780 | 12.27 | 15.60 | 21.20 | ||

7.013 | 6.405 | 7.724 | 11.22 | 13.60 | 18.72 | ||

10.77 | 11.04 | 12.23 | 13.84 | 14.71 | 20.81 | ||

11.42 | 16.21 | 15.17 | 14.75 | 16.86 | 17.79 | ||

[1] | 0.2059 | 0.1509 | 0.1165 | 0.1863 | 0.2534 | 0.8370 | |

[] | 7.244 | 22.61 | 24.29 | 23.76 | 23.25 | 23.43 | |

[] | 7.990 | 22.65 | 24.84 | 24.12 | 19.91 | 22.37 | |

[] | 6.249 | 17.83 | 22.77 | 22.67 | 23.25 | 23.03 | |

[1] | 0.1692 | 0.2208 | 0.2075 | 0.1931 | 0.2038 | 0.2036 | |

[1] | 0.01060 | 0.05870 | 0.09630 | 0.1739 | 0.2492 | 0.2779 | |

[1] | 0 | 0 | |||||

[1] | 0.1415 | 0.2813 | 0.3384 | 0.3761 | 0.3828 | 0.6193 |

Rows are: galactocentric distance , spherically averaged density, radius of spheres, ensemble averaged mass in spheres, in shell spherically averaged value of density, mean radial velocity, mean -velocity, mean -velocity, total velocity dispersion, radial velocity dispersion, -velocity dispersion, -velocity dispersion, anisotropy parameter, angle between long axis of shape ellipsoid and long axis of velocity ellipsoid, angle between intermediate axis of shape ellipsoid and intermediate axis of velocity ellipsoid, angle between short axis of shape ellipsoid and short axis of velocity ellipsoid and triaxiality parameter of the velocity dispersion ellipsoid. Further we have the ensemble dispersions of the above mentioned quantities. The last three rows are the fraction of spheres that are affected by subhaloes, the fraction of empty spheres and the Gini coefficient for the probability density functions of the local mean density.

For the analysis presented here, we used data from the Via Lactea II simulation (Diemand et al., 2008) which simulates the assembly of a Milky Way size cold dark matter halo within a CDM universe. The initial conditions at a starting redshift of consist of a 40 comoving Mpc periodic box and were generated with a modified, parallel version of GRAFIC2 (Bertschinger, 2001). We use the traditional method of refining a region of interest with a large number of particles and leaving the rest on lower resolution so that we correctly account for the large scale tidal forces (Katz & White, 1993; Bertschinger, 2001). In total, the simulation consists of more than high resolution particles with a mass of and a gravitational softening length of 40 pc. We used the WMAP 3-year cosmological parameters (Spergel et al., 2007) with , and , and .

The time evolution until redshift was performed with the parallel tree-code PKDGRAV2 (Stadel, 2001). PKDGRAV2 uses a fast multipole expansion technique in order to calculate the forces with hexadecapole precision and a time-stepping scheme that is based on the true dynamical time of the particles with an accuracy parameter of (Zemp
et al., 2007; Stadel &
Zemp, 2008). The simulation used close to CPUh on the Jaguar Cray XT3 super-computer at the Oak Ridge National Laboratory^{2}^{2}2http://www.ornl.gov.

Throughout the paper, we only use data from the snapshot. We present all quantities in physical units and in a coordinate system centred on the particle with the deepest potential within the dark matter halo. A detailed discussion of the Via Lactea II simulation at can be found in Diemand et al. (2008) and we report here only the main characteristics. The radius , where the enclosed density is 200 times the mean matter density (background), is and contains a total mass of . The maximum circular velocity of the halo is and is reached at the radius . We determine the shape of the halo at 20 kpc by diagonalising the shape tensor^{3}^{3}3Often the tensor is incorrectly denoted as the inertia tensor. However, the tensor does not give the correct relation between the angular momentum vector L and the angular velocity vector which is given by where is the correct definition of the inertia tensor. Of course, the tensors and have the same eigenvectors since they differ only by a diagonal tensor with one eigenvalue of algebraic multiciply of 3. Though, the eigenvalues of and are obviously different and have entirely different meanings. with elements

(1) |

through an iterative procedure that adapts to the local shape as described in Katz (1991) so that the summation over the particles only contains those within the triaxial ellipsoid. We then rotate the whole halo so that the long axis is the -axis, the intermediate axis is the -axis, and the short axis is the -axis. The resulting axis ratios at 20 kpc are and where are the square roots of the eigenvalues of the shape tensor. Of course, a dark matter halo is not a rigid body and the axis ratios as well as their orientations can change with radius, with haloes typically becoming rounder in their outer part. The orientation of the principle axis is constant with radius (to within a few degrees) between 20 kpc and 400 kpc.

## 3 Local Properties

### 3.1 Procedure

In order to measure local properties of the dark matter halo, we chose the following procedure. We consider six distinct distances: 8, 25, 50, 100, 200 and 400 kpc from the halo centre. At each distance we randomly sample spheres with a radius depending only on the distance and chosen to include of particles. Specifically we set and calculate the radii at other distances by

(2) |

where is the measured spherically averaged density at radius . The number of sampling spheres is large enough so that even at 400 kpc the whole sky, as seen from the centre of the galaxy, is covered more than once. A summary of the characteristics of the shells at different distances from the galactic centre is given in Table 1. In order to save space, in most cases we only show the plots for 8, 25, 100 and 400 kpc in the following analysis.

In order to assess the influence of the triaxial shape of the halo, we also consider 3 sub-sets of spheres whose centres lie within a cone of 20 along the long (), the intermediate () respectively the short () axis, i.e. spheres that fulfil the condition

(3) |

where is the unit vector in galactocentric coordinates towards the centre of the sphere and is the unit vector along the -, - respectively -axis. It is expected that each sub-set contains approximately 6% of the total number of spheres and the actual numbers lie within the statistical range. The slightly twisted shape of the dark matter halo has only a small influence on the sub-sets along the different axes at different distances since the deviations of the eigenvector directions are only a few degrees which is small compared to the cone opening angle of 20.

Furthermore we consider the sub-sample of spheres that are affected by subhaloes. We take a group catalogue of the Via Lactea II simulation that contains all subhaloes with a peak circular velocity greater than , of which there are 26182 within 500 kpc. These groups are identified by a friends-of-friends method in phase space (6DFOF) as described in Diemand et al. (2006). By a density profile fitting algorithm, we determine a size for each subhalo. A sphere with centre at location is then deemed affected by a subhalo if the criterion

(4) |

is fulfilled for any subhalo with radius and centred at , i.e. the subhalo at least partially overlaps the sphere. We give the fraction of spheres that are affected by subhaloes in Table 1. This fraction increases with distance from the centre of the galaxy, reaching a value of approximately 28% at 400 kpc. The fraction of subhalo-affected spheres is a lower limit since more subhaloes would survive in higher resolution simulations at all distances from the centre of the galaxy.

Throughout this paper we use the following notation scheme: quantities calculated by averaging over particles within a sphere are denoted by , quantities calculated by averaging over particles within a shell from to are denoted by , and averages over the ensemble of spheres at radius are denoted by , where we also might use indices in order to denote ensemble averages over only sub-sets of the spheres, for any given quantity .

### 3.2 Local densities

In Fig. 1, we plot the probability density function of the local mean density within the spheres in logarithmic scale. The local density is simply defined by the mass within a sphere divided by the volume of the sphere. In the central region, the structure of the halo looks pretty smooth in density. At 8 kpc, we find a maximum total spread of a factor 2 in local density with a standard deviation of . This is due to the triaxiality of halo since the underdense spheres lie along the -axis and the overdense spheres along the -axis. In the outer part of the halo, however, the smooth picture does not apply. For example at 200 kpc and 400 kpc, underdense regions are found away from the short axis and the distributions of the sub-samples along the other axes are broader than in the centre of the halo. The possible range in densities at 400 kpc is close to 4 orders of magnitude with a peak probability at around , i.e. the typical value for the local densities is below the spherical average value. Beyond 100 kpc, the standard deviation becomes larger than the mean value, e.g. at 400 kpc .

Naively, one would expect a decrease in the density variation due to shape in the outskirts of the halo, since its shape becomes rounder further out. However, this is exceeded by the steeper fall off of the density profile in the outskirts of the dark matter halo. The actual density contrast between the mean local density of the -axis sample with respect to the mean local density of the -axis, , increases from to .

It is also striking that there are no high density spikes from subhaloes in the central region. This is just an artefact of our rather conservative definition of locality (i.e. the large size of the spheres). There are only very few surviving and well resolved subhaloes within 8 kpc in Via Lactea II. The subhaloes that survive close to the centre are small and compact due to tidal mass loss so that with our still relatively large definition of locality, their contribution to the mass of the spheres is small. By decreasing the size of the spheres the extremes of the distribution become more populated (see also appendix A). With higher resolution, one expects more subhaloes to survive which are even more compact. But one would also decrease the size of the spheres that define locality in our scheme. Hence, it is not clear what the outcome of these two competing effects would be. From simple analytical models, one expects some contribution from these small haloes (Kamionkowski & Koushiappas, 2008) though the details depend on their volume filling factor.

While high peaks from subhaloes are expected, underdense regions away from the short axis of the shape ellipsoid are surprising. In order to further investigate these underdensities, we repeated the measurement of local properties with spheres of 4 times smaller radii, i.e. . In these spheres, one would expect on average approximately 21 particles in the inner region. We then simply count the small spheres that contain no particles and calculate from that the fraction of empty spheres (given in Table 1). Assuming balls-in-bins statistics (see appendix B for more details), the probability of getting empty spheres is given by . Hence one would not expect to find any sphere of size that is completely empty. The fact that approximately 2% of all these smaller spheres at 400 kpc are empty shows clearly that the outskirts of dark matter haloes are far from smooth in position space. Of course, to some degree this is due to the triaxial shape as actually most of the empty spheres in the outskirts are in low density regions along the -axis. Taking the lower density into account, one would expect approximately 6 particles per sphere in the -axis sample at 400 kpc. The probability of getting empty spheres is . The measured fraction of empty spheres in the -sample is which is still clearly much higher.

We quantify this further with a statistical indicator that measures inequality: the Gini coefficient. The Gini coefficient is defined as where is the area under the Lorenz curve of the probability distribution (Lorenz, 1905; Gini, 1912). For the discrete probability density functions of the local mean density this results in

(5) |

with

(6) |

and

(7) |

with being the bin width in logarithmic scale. The range of the Gini coefficient is between 0 and 1. A Gini coefficient of 0 means that all the spheres have equal density whereas a Gini coefficient of 1 means that all the mass is in one sphere. As can be seen in Table 1, the Gini coefficient increases with galactocentric distance showing again that the outskirts of dark matter haloes have a clumpy structure.

This clumpy structure is due to the hierarchical build up of dark matter haloes by accretion of subhaloes. These subhaloes leave tidal streams along their orbits when they are falling into the host halo. In the outer parts of the halo we find many overlapping streams locally (see also section 3.7) that originate from the many subhaloes that passed a certain region. Combined with the long local dynamical time-scale in the outskirts of haloes which inhibits effective mixing, this leads to the clumpy structure and large spread in density we see on small local scales in the outer parts of dark matter haloes.

The tidal debris can have an effect on the local dark matter annihilation (Afshordi et al., 2008). We can quantify this by the boost factor . We can only say something about local boost factors since by averaging over a spherical shell, the boost factor would be mainly driven by the triaxiality. The biggest effect is expected in the outskirts of the halo since there we find the highest degree of clumpiness. In order just to measure the influence of the debris structure and exclude the effect of the subhaloes, we only include spheres that are not subhalo-affected and get values from to . This is in agreement with Afshordi et al. (2008) who get of for the boost factor originating in the structure of the debris. If we include the subhalo-affected sample we get values in the range from to locally, showing showing that even in the outer halo most of the clumpiness comes from the subhalos themselves and not from their more diffuse tidal debris.

Generally, all the probability density function plots in this paper are lower limits for a possible real variation of local properties in nature since numerical effects due to numerical under-resolving (even in Via Lactea II) lead to artificial heating, smoothing, and in general loss of structure. A more detailed discussion of numerical effects follows in appendix A.

### 3.3 Local mean velocities

In this section we investigate the properties of the local mean velocities, which we calculate by number weighted averages within the spheres. For this purpose, we split the velocity vector field in a radial, azimuthal ( defined as the angle in the -plane from the -axis) and a polar ( defined as the angle form the -axis) component so that the unit vectors denote a right handed system with . In Fig. 2 and 3, we plot the probability density functions for the radial respectively -component of the velocity field in linear scale; omitting the plot for the mean -velocity as it does not show qualitatively different features than the mean -velocity.

Often it is assumed that bulk velocities in dark matter haloes are zero. Again, this is approximately true in the central part. The standard deviation of the local mean radial velocity is about 4% of the shell averaged radial velocity dispersion. In the outer parts, not even the spherically averaged values are zero (see also Table 1) and the local bulk velocities reach up to twice the values of the local velocity dispersion. At 400 kpc, the dispersion of the mean radial velocity exceeds 40% of the shell averaged radial velocity dispersion.

At different radii one can see two populations of spheres: one that is infalling and another one that is moving outwards. This is best seen along the -axis but often the two populations are too much smeared out. Such a distribution naturally arises from a cosmological infall pattern and has been seen before in numerical simulations (Diemand & Kuhlen, 2008).

A similar picture emerges from the mean - and -velocities with a comparable spread of the probability density distributions at a given distance as for the mean radial velocities. We find that the sub-population of spheres along the long axis have a narrower and more centrally peaked distribution in mean - and -velocities than the sub-populations along the short and intermediate axes indicating some degree of coherent tangential flow around the short and intermediate axes.

The distribution of the subhalo-affected sample of spheres is always much broader than the total sample, and it contributes most to the extremes (both high and low) of the overall distribution. This indicates that the subhalo population as a whole is kinematically hotter than the background, in agreement with earlier studies by Diemand et al. (2004). Indeed we find subhaloes in the inner halo that have a speed of around .

In the inner region, we find that the subhalo-affected sample is skewed towards negative radial velocities (i.e. infalling spheres) and positive -velocities. The skew in radial velocity is a signature of the tidal disruption process: we see these subhaloes on their last infall before they either are completely disrupted or they loose so much mass that they do not contribute much to the sphere averaged properties any more on their way out since they are now much more compact and less massive.

All this conclusively shows that the traditional notion of a dark matter halo with no local bulk velocities is an inaccurate description for structures that form in cosmological simulations, especially in the outer parts of haloes.

### 3.4 Local velocity dispersions

In Fig. 4, 5 and 6 we plot the probability density functions of the local total, radial and -velocity dispersion. We omit the figure for the -velocity dispersion since it shows a similar behaviour as the -velocity dispersion. The local velocity dispersion is given by for and , where , respectively are the values of the dispersions in the eigencoordinate system of the local dispersion ellipsoid (see section 3.6 for more details).

We see that the local velocity dispersion has strong underlying variations. Whereas the central part only has variations on the few per cent level with respect to the spherically averaged value (e.g. at 8 kpc we have ), one can find regions in the outer part of the halo that deviate by approximately an order of magnitude in velocity dispersion and at 400 kpc we have for example .

The total velocity dispersion shows the most compact probability density function when measured with respect to the spherically averaged value, i.e. is smaller than for and the distributions peak generally around the spherically averaged value.

Interestingly, we find that regions in the centre along the long axis are colder (i.e. is smaller) than spheres along the intermediate or short axes. This trend is not observed in the outskirts of the halo. We also find that along the long axis, the radial velocity dispersions are higher and the tangential velocity dispersions are low in the inner region of the halo. This is due to the prolate shape and orientation of the local velocity dispersion ellipsoid and is discussed in more detail in section 3.6. In general, we find that the distributions of velocity dispersions for the two tangential velocity components are much broader along the short axis than along the long and intermediate axes. The subhalo-affected sample shows a less peaked distribution in velocity dispersions than the total sample. In the regions with enough subhaloes, the distribution is skewed towards lower dispersions for all velocity components. This is of course due to the lower internal velocity dispersion in subhaloes.

It remains to be seen if these empirical trends of velocity dispersion with orientation are universal for dark matter haloes or if they depend on the detailed hierarchical build-up history of every individual halo. Nevertheless, it is clear that the degree of chaotic motion can strongly vary locally and deviate by substantial factors from the spherically averaged value.

### 3.5 Local anisotropy parameter

In Fig. 7, we plot the probability density function of the local anisotropy parameter within the spheres. The local anisotropy parameter is given by

(8) |

where is the tangential velocity dispersion squared and where we neglect the correlation term.

This form of the anisotropy parameter is often used but has some problems since the following assumptions which are hidden in this expression are in general not fulfilled: i) , ii) and iii) . These three conditions are only approximately fulfilled in the central part of a halo and certainly not in the outskirts. One should therefore see the expression for as a definition for local anisotropy so that a comparison with previous work is possible.

We find that the radial anisotropy of the velocity dispersion tensor in our halo increases with radius (see also Table 1), in agreement with previous studies (see e.g. Cole & Lacey, 1996; Colín et al., 2000; Fukushige & Makino, 2001; Diemand et al., 2004; Hansen & Moore, 2006). In the central part we find that regions along the long axis are preferentially on radial orbits, whereas regions along the short axis are tendencially on tangential orbits. This effect disappears at larger radii. It can be understood as a direct consequence of the variation of the orientation of the local velocity dispersion ellipsoid discussed in the next section 3.6.

The probability densities peak in general close to the spherical average value except at 400 kpc where it is peaked towards highly radial anisotropies. This is especially true for the the sub-profiles of the - and -axis sample whereas the -axis sample peaks around . With increasing distance also the spread in local anisotropy becomes larger and can vary from close to perfectly radial to close to perfectly tangential. By inspecting the subhalo-affected sample, we find that this sample is skewed towards tangential values of at galactocentric distances with enough subhaloes.

### 3.6 Local velocity dispersion ellipsoid

In each sphere we also calculate the local velocity covariance tensor, also know as the velocity dispersion tensor, given by

(9) |

By diagonalising , the velocity dispersions in the eigencoordinate system are given by the square roots of the eigenvalues of the velocity covariance tensor. These velocity dispersions in the eigencoordinate system define the local velocity dispersion ellipsoid and are sorted by . We name the appropriate eigenvectors accordingly, e.g. etc.

First, we check the orientation of the local velocity dispersion ellipsoid. For this purpose we calculate the angles between the eigenvectors of the velocity covariance tensor and the appropriate shape axis given by

(10) | |||||

(11) | |||||

(12) |

In Fig. 8 respectively 9 we plot the probability density function of and . We do not show the figure for since it shows a qualitatively similar behaviour as for . We also plot the curve, which is proportional to , that corresponds to an isotropic distribution.

We observe that in the inner regions the local velocity dispersion ellipsoid is close to perfectly aligned with the shape ellipsoid since all angle probability density functions peak sharply at small angles. This effect is most prominent for the major axis sub-sample, which has a very tight and sharp distribution around small angles for all , and less so for the intermediate and short axes, which in some cases show a broader spread towards larger angles. The general alignment of the local velocity dispersion ellipsoid with the shape of the dark matter halo nicely explains the observed behaviour of the anisotropy parameter , for which we mainly found radial orbits along the -axis, isotropic orbits along the intermediated -axis and tangential orbits along the -axis. The further out we go, the more isotropic the different angle distributions become. However, at nearly all distances from the Galaxy centre we find that smaller angles are slightly more probable than in a perfectly isotropic distribution.

Interesting is the behaviour of the different sub-samples further out. The alignment with the shape seems to be best along the -axis. Only in the outskirts at 400 kpc, this sub-sample distribution becomes more and more isotropic for all angles. The sub-samples along the other two axes show even in the inner part a broader distribution that does not follow the overall distribution. For example, the -axis sample shows a distribution peaked around large angles for and for intermediate distances around 50 kpc and 100 kpc whereas the distribution is mildly peaked towards small angles at 50 kpc (not shown in figure). In other words: only the long axis seems to be slightly aligned with the long axis of the shape ellipsoid whereas the the intermediate and short axes are perpendicular to the corresponding shape axes. In general, also the subhalo sample shows a rather isotropic distribution although it also tends to be aligned with the shape ellipsoid in the inner part of the halo.

Globally, we expect such an alignment from the tensor virial theorem (Binney & Tremaine, 1987) and such a correlation has previously been found in cosmological -body simulations, e.g. in Allgood et al. (2006). If we calculate the velocity dispersion ellipsoid for all the particles in the shell, we also find a very tight alignment of this shell dispersion ellipsoid with the shape ellipsoid (see e.g. Fig. 8 respectively 9). Only in the outer parts at 400 kpc we get a deviation from that behaviour. But in principle, we do not expect that the tensor virial theorem holds locally. It seems that the alignment of the local velocity dispersion ellipsoid with the shape ellipsoid only holds in relaxed and well mixed regions - the central region of the Via Lactea II halo. Hence, this local alignment might just be a numerical artefact since the central region is too relaxed due to numerical under-resolving if compared to the true degree of relaxation in reality. But an alignment of the local dispersion ellipsoid with the shape is also found in observations. For example, galaxies from the SAURON survey that were dynamically modelled with the Schwarzschild technique also show such a correlation (Cappellari et al., 2007).

Furthermore, we calculate the shape of the local velocity dispersion ellipsoids which we measure with the triaxiality parameter (Franx et al., 1991) defined by

(13) |

Ellipsoids are called oblate if , triaxial if and prolate if . In Fig. 10 we plot the probability density function of the triaxiality parameter in linear scale.

The total probability density function shifts from a peak in the oblate region (e.g. at 25 kpc) via a peak in the triaxial region (e.g. at 100 kpc) to a more prolate shape at 400 kpc. Only the innermost region at 8 kpc does not follow this overall trend. Interesting is again the distribution of the different sub-samples. The shape of the velocity dispersion ellipsoids along the -axis generally peak in the prolate region although the distribution becomes broader towards oblate shapes further out. The -axis sample has in the inner region a rather oblate shape which drifts via triaxial (e.g. at 100 kpc) to a prolate distribution in the outskirts of the halo. Completely different is the behaviour of the -axis sample: at 8 kpc it’s distribution is peaked in the prolate region then drifts via the triaxial region at 25 kpc to a oblate distribution at 50 kpc (not shown in the figure). Then it swings back via a rather triaxial shape to a prolate shape again at 400 kpc. The shape of the shell averaged velocity dispersion tensor does not fully follow the trend of the overall distribution and is in general not close to the peak of the total probability density function.

At the moment it is not clear what the underlying cause for these shape and orientation variations with galactocentric distance is and if these trends are universal. Further numerical investigations are necessary.

### 3.7 Local velocity space

We now turn to a closer look at the local velocity space structure. For this purpose, we present 3D visualizations of the velocity structure within spheres cut out along the -axis, with radii given by from Table 1.

Fig. 11 shows the positions and velocity vectors for every particle within these spheres (1st and 3rd row) for all six galactocentric distances we investigated. We also plot the location of the particles in the local velocity space (2nd and 4th row). The colour (and length in position space) encodes the magnitude of the velocity vector. The big white arrow in the centre of the position space plots points towards the galactic centre whereas the black cube in the centre of the local velocity space plots marks the origin. For the velocity space plots, we split the velocity vector field in a radial- (-axis), - (-axis) and a -component (-axis).

These spheres have been selected to be free of subhaloes. Therefore, one would not expect to see clumpy structure in the phase space of these spheres. In the inner spheres no velocity space structure is apparent. The orientation of the velocity vectors appears random and the actual velocity distribution can be well fit by a multivariate normal distribution (the different velocity dispersion components are given in Table 1). But further out, there is evidence for locally coherent motion, visible as groups of vectors with the same colour pointing in the same direction. This is even clearer in the velocity space plots, which exhibit an increasing degree of clumpiness the further out one goes. At 400 kpc, for example, there are only a bit more than a dozen clumps in velocity space and no smooth component at all. Obviously a smooth multivariate normal distribution would not provide a good fit. This is evidence that the outer part of the halo is built up from a collection of large scale streams from the tidal disruption of infalling subhaloes. But so too is the smooth component in centre, in that it also likely consists of an overlap of many, many streams. The central limit theorem then guarantees that the resulting distribution in the centre closely resembles something like a multivariate normal distribution (Helmi et al., 2003) or a generalisation thereof (Hansen et al., 2006).

In Fig. 12 we show the central region of Via Lactea II: the sphere has a radius of 50 kpc. We plot the position space density (left) and the true phase space density (right) calculated with EnBiD^{4}^{4}4http://sourceforge.net/projects/enbid/ (Sharma &
Steinmetz, 2006). It is important here to calculate the true phase space density in 6D and not the commonly used pseudo phase space density (where is the velocity dispersion), since in the latter information in velocity space is lost by averaging over the particles in local position space, instead of calculating the density in local velocity space directly from the particles. We only show the top 5 to 6 orders of magnitude in position and phase space density, in order not to overload the two pictures.

Many cold streams are clearly visible in the central region. Although these cold streams only contribute a few percent to the local mass density, their velocity dispersion is just a few , resulting in a very high phase space density for these particles. These streams are not visible in traditional density or density squared pictures and can only be revealed by visualising the true phase space density. Via Lactea II seems to be the first structure formation simulation with sufficient resolution to reveal these streams in phase space, since previous lower resolution runs did not show such phase space features. Also subhaloes are better visible due to their higher contrast in phase space density. Approximately 2000 peaks are seen in the central 50 kpc sphere phase space density image. The high contrast makes the phase space density the ideal method of finding subhaloes.

From our findings it is obvious that also the central region shows a huge amount of additional structure aside from the expected subhalo density peaks. A more detailed discussion about streams and their properties will follow in a future publication.

## 4 Discussion and conclusions

In this paper we present a study of local properties of the Via Lactea II dark matter halo. Commonly, characteristics of dark matter haloes are described by spherically symmetric profiles. Here, we show that locally, at a fixed galactocentric distance, these properties can vary by orders of magnitude from their canonical, spherically averaged values. For example, while the density is smooth in the centre, in the outskirts the density range spans four orders of magnitude. This is due to the presence of both subhaloes and underdense regions or holes in the matter distribution. The widespread assumption that there are no bulk flows is only warranted in the central region, where we also find the local velocity dispersion ellipsoid to be aligned with the shape ellipsoid of the halo. Neither of these findings hold in the outer parts of the halo, where particles exhibit bulk motions and a more isotropic distribution of the velocity dispersion ellipsoid’s orientation is found. The local velocity space structure can, in general, not be well described by a smooth multivariate normal distribution, at least not in the outer parts of the halo. A qualitatively new feature in Via Lactea II is the detection of streams that are made visible only through their true 6D phase space density.

Such variations of local properties is due to both the triaxial shape of the dark matter halo and its generally clumpy structure in phase space. We find that the phase space structure of a dark matter halo shows a significant departure from the canonical picture of a smooth background density profile with subhaloes in position space and multivariate normal distributions in velocity space. Spherically averaged quantities do not adequately describe properties of the dark matter halo at a given galactocentric distance - especially not the wealth of structure we find locally - since a lot of information gets smeared out by spherically averaging. In general, dark matter haloes show a high degree of graininess - clumps in phase space - which will probably show up at even smaller scales that now seem to be smooth in future simulations with higher resolution. Therefore, a smooth and featureless distribution function does not accurately describe dark matter haloes that form in a cosmological structure formation simulation.

We find several correlations between the shape ellipsoid and the local velocity dispersion ellipsoid. It is unclear at the moment to what degree our findings are universal or how these correlations depend on the hierarchical build up history.

Knebe & Wießner (2006) earlier estimated analytically the effect of a triaxial halo shape on the variance of the local density. Their values of the dispersion normalised to the spherically averaged value for a halo with a similar triaxiality correspond quite good with our values in the centre of the halo but we get much higher dispersion values (larger than the mean) in the outskirts of the halo. This difference is due to their assumption of a smooth triaxial density profile and we get a higher scatter due to the presence of subhaloes and underdense regions.

For dark matter direct detection experiments it is essential to know the local dark matter properties at 8 kpc. The spherically averaged value for the density in Via Lactea II is , close to which is often used as a canonical value in the literature (Kamionkowski & Kinkhabwala, 1998; Particle Data Group, 2008). But as we have shown here, this value can vary locally. One of the large uncertainty factors is the missing information about the orientation of the disk with respect to the dark matter halo. There are different claims from observations and theory about a possible alignment of the angular momentum axis of the disk with the shape axes of the halo (see e.g. Navarro et al., 2004; Bailin et al., 2005; Sharma & Steinmetz, 2005, and references therein) so that a clear answer is not possible at the moment. But often it is claimed that the angular momentum axis of the disk and the short axis () tend to be aligned (Bailin et al., 2005) which would mean that the disk would lie preferentially in the -plane in our coordinate system. An additional problem is that we measure these local properties within a sphere of radius but for dark matter detection experiments more a scale of is relevant and it is not clear what the local properties of the dark matter distribution on a 1 AU scale are or how one could reasonably extrapolate that over 8 orders of magnitude - especially considering the highly non-linear numerical effects mentioned above and in the appendix that affect the local phase space structure. Also the missing baryonic physics in Via Lactea II like adiabatic contraction, stellar disk and bulge, inspiralling compact objects like black holes etc. can modify the central dark matter structure in either way. Therefore, it is still not clear what the detailed structure of the dark matter locally is.

## Acknowledgement

The Via Lactea II simulation was carried out at the Oak Ridge National Laboratory through an award from the Department of Energy’s Office of Science as part of the 2007 Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program^{5}^{5}5http://www.science.doe.gov. It is especially a pleasure to thank Bronson Messer and the Scientific Computing Group at the National Center for Computational Sciences^{6}^{6}6http://www.nccs.gov for their support and help. It is also a pleasure to thank Glenn van de Ven for helpful comments. Further supporting computations for generating the initial conditions, code optimisation and analysis were performed on MareNostrum at the Barcelona Supercomputing Centre^{7}^{7}7http://www.bsc.org.es, the zBox at University of Zurich, Columbia at NASA Ames^{8}^{8}8http://www.nas.nasa.gov and on Pleiades at University of California Santa Cruz^{9}^{9}9http://pleione.ucsc.edu/pleiades. M. Z. gratefully acknowledges support by the Swiss National Science Foundation. Additionally, support for this work was provided by NASA through grants HST-AR-11268.01-A1 and NNX08AV68G (P. M.) and Hubble Fellowship grant HST-HF-01194.01 (J. D.). M. K. acknowledges support from the William L. Loughlin Fellowship at the Institute for Advanced Study.

## References

- Afshordi et al. (2008) Afshordi N., Mohayaee R., Bertschinger E., 2008, preprint (arXiv:0811.1582)
- Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
- Baes & van Hese (2007) Baes M., van Hese E., 2007, A & A, 471, 419
- Bailin et al. (2005) Bailin J., Kawata D., Gibson B. K., Steinmetz M., Navarro J. F., Brook C. B., Gill S. P. D., Ibata R. A., Knebe A., Lewis G. F., Okamoto T., 2005, ApJ Lett., 627, L17
- Bertschinger (2001) Bertschinger E., 2001, ApJ Suppl., 137, 1
- Bett et al. (2007) Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215
- Binney & Tremaine (1987) Binney J., Tremaine S., 1987, Galactic Dynamics, 3. edn. Princeton University Press, Princeton, NJ USA
- Cappellari et al. (2007) Cappellari M., Emsellem E., Bacon R., Bureau M., Davies R. L., de Zeeuw P. T., Falcón-Barroso J., Krajnović D., Kuntschner H., McDermid R. M., Peletier R. F., Sarzi M., van den Bosch R. C. E., van de Ven G., 2007, MNRAS, 379, 418
- Casetti-Dinescu et al. (2008) Casetti-Dinescu D. I., Carlin J. L., Girard T. M., Majewski S. R., Peñarrubia J., Patterson R. J., 2008, AJ, 135, 2013
- Cole et al. (2008) Cole N., Newberg H. J., Magdon-Ismail M., Desell T., Dawsey K., Hayashi W., Liu X. F., Purnell J., Szymanski B., Varela C., Willett B., Wisniewski J., 2008, ApJ, 683, 750
- Cole & Lacey (1996) Cole S., Lacey C., 1996, MNRAS, 281, 716
- Colín et al. (2000) Colín P., Klypin A. A., Kravtsov A. V., 2000, ApJ, 539, 561
- Dehnen & McLaughlin (2005) Dehnen W., McLaughlin D. E., 2005, MNRAS, 363, 1057
- Diemand & Kuhlen (2008) Diemand J., Kuhlen M., 2008, ApJ Lett., 680, L25
- Diemand et al. (2006) Diemand J., Kuhlen M., Madau P., 2006, ApJ, 649, 1
- Diemand et al. (2007) Diemand J., Kuhlen M., Madau P., 2007, ApJ, 657, 262
- Diemand et al. (2008) Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J., 2008, Nature, 454, 735
- Diemand et al. (2005) Diemand J., Madau P., Moore B., 2005, MNRAS, 364, 367
- Diemand et al. (2004) Diemand J., Moore B., Stadel J., 2004, MNRAS, 352, 535
- Diemand et al. (2005) Diemand J., Zemp M., Moore B., Stadel J., Carollo C. M., 2005, MNRAS, 364, 665
- Dubinski & Carlberg (1991) Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
- Fairbairn & Schwetz (2008) Fairbairn M., Schwetz T., 2008, preprint (arXiv:0808.0704)
- Franx et al. (1991) Franx M., Illingworth G., de Zeeuw T., 1991, ApJ, 383, 112
- Fukushige & Makino (2001) Fukushige T., Makino J., 2001, ApJ, 557, 533
- Gini (1912) Gini C., 1912, Variabilità e mutabilità. Libreria Eredi Virgilio Veschi
- Hansen & Moore (2006) Hansen S. H., Moore B., 2006, New Astron., 11, 333
- Hansen et al. (2006) Hansen S. H., Moore B., Zemp M., Stadel J., 2006, JCAP, 1, 14
- Helmi et al. (2002) Helmi A., White S. D., Springel V., 2002, Phys. Rev. D, 66, 063502
- Helmi et al. (2003) Helmi A., White S. D. M., Springel V., 2003, MNRAS, 339, 834
- Hernquist (1993) Hernquist L., 1993, ApJ Suppl., 86, 389
- Ibata et al. (2002) Ibata R. A., Lewis G. F., Irwin M. J., Quinn T., 2002, MNRAS, 332, 915
- Kamionkowski & Kinkhabwala (1998) Kamionkowski M., Kinkhabwala A., 1998, Phys. Rev. D, 57, 3256
- Kamionkowski & Koushiappas (2008) Kamionkowski M., Koushiappas S. M., 2008, Phys. Rev. D, 77, 103509
- Katz (1991) Katz N., 1991, ApJ, 368, 325
- Katz & White (1993) Katz N., White S. D. M., 1993, ApJ, 412, 455
- Kazantzidis et al. (2004) Kazantzidis S., Magorrian J., Moore B., 2004, ApJ, 601, 37
- Knebe & Wießner (2006) Knebe A., Wießner V., 2006, Publ. Astron. Soc. Aust., 23, 125
- Kuhlen et al. (2007) Kuhlen M., Diemand J., Madau P., 2007, ApJ, 671, 1135
- Kuhlen et al. (2008) Kuhlen M., Diemand J., Madau P., 2008, ApJ, 686, 262
- Lorenz (1905) Lorenz M. O., 1905, Publ. Am. Stat. Assoc., 9, 209
- Martínez-Delgado et al. (2007) Martínez-Delgado D., Peñarrubia J., Jurić M., Alfaro E. J., Ivezić Z., 2007, ApJ, 660, 1264
- Moore et al. (2001) Moore B., Calcáneo-Roldán C., Stadel J., Quinn T., Lake G., Ghigna S., Governato F., 2001, Phys. Rev. D, 64, 063508
- Moore et al. (1998) Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ Lett., 499, L5
- Navarro et al. (2004) Navarro J. F., Abadi M. G., Steinmetz M., 2004, ApJ Lett., 613, L41
- Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
- Particle Data Group (2008) Particle Data Group 2008, Phys. Lett. B, 667, 212
- Sharma & Steinmetz (2005) Sharma S., Steinmetz M., 2005, ApJ, 628, 21
- Sharma & Steinmetz (2006) Sharma S., Steinmetz M., 2006, MNRAS, 373, 1293
- Siegal-Gaskins & Valluri (2008) Siegal-Gaskins J. M., Valluri M., 2008, ApJ, 681, 40
- Spergel et al. (2007) Spergel D. N., Bean R., Doré O., Nolta M. R., Bennett C. L., Dunkley J., Hinshaw G., Jarosik N., Komatsu E., Page L., Peiris H. V., Verde L., Halpern M., Hill R. S., Kogut A., Limon M., Meyer S. S., Odegard N., Tucker G. S., Weiland J. L., Wollack E., Wright E. L., 2007, ApJ Suppl., 170, 377
- Stadel et al. (2008) Stadel J., Potter D., Moore B., Diemand J., Madau P., Zemp M., Kuhlen M., Quilis V., 2008, preprint (arXiv:0808.2981)
- Stadel (2001) Stadel J. G., 2001, PhD Thesis, University of Washington
- Stadel & Zemp (2008) Stadel J. G., Zemp M., 2008, in preparation
- Stiff & Widrow (2003) Stiff D., Widrow L. M., 2003, Phys. Rev. Lett., 90, 211301
- Stiff et al. (2001) Stiff D., Widrow L. M., Frieman J., 2001, Phys. Rev. D, 64, 083516
- Van Hese et al. (2008) Van Hese E., Baes M., Dejonghe H., 2008, preprint (arXiv:0809.0901)
- Vogelsberger et al. (2008) Vogelsberger M., White S. D. M., Helmi A., Springel V., 2008, MNRAS, 385, 236
- Zemp et al. (2008) Zemp M., Moore B., Stadel J., Carollo C. M., Madau P., 2008, MNRAS, 386, 1543
- Zemp et al. (2007) Zemp M., Stadel J., Moore B., Carollo C. M., 2007, MNRAS, 376, 273
- Zhang & Magorrian (2008) Zhang M., Magorrian J., 2008, MNRAS, 387, 1719

## Appendix A Comparisons

In this section we investigate the influence of different definitions of locality and numerical resolution on our local estimates.

### a.1 Definition of locality

Our definition of locality is somewhat arbitrary. If the region is too big then local properties get smeared out by the averaging procedure and if the region is too small then statistical fluctuations become important. Hence, we chose the size of our spheres so that they contain particles as in similar previous work (Moore et al., 2001; Helmi et al., 2002). Exactly: an average density sphere of radius 1 contains 1356 particles. In order to see the effect of different choices for our local estimate, we repeated the exercise from section 3 with spheres of different radii of the following size: , , and , with a corresponding increase or decrease in the number of particles sampled per sphere. The resulting Poisson scatter ranges from 22% for the smallest spheres () to 1% for the largest spheres () if we take the spherical averaged density as a reference. By only looking at the lowest density along the -axis we get a range of 41% for the smallest spheres to 1.7% for the largest spheres.

In Fig. 13, we present the effect of the different definitions of locality on the local density distribution and see that they show the expected result. The regions that contain 8 times more particles are of course smoothed to a higher degree so that the resulting spread in the probability distribution function is a bit reduced. In a similar way the spread for the 8 times smaller spheres is increased. But the different definitions of locality mainly affect the rare outliers on the tails of the distribution. By reducing to 64 times less particles, the extremes of the distribution become much more populated, which is mostly due to the increased Poisson scatter and to a lesser extend due to increased graininess on the smaller probed scales.

Any size of a sphere that contains at least a few hundred particles would be fine for a definition of locality, indicating that our choice is a rather conservative choice and the resulting variances of the different distributions from the previous sections are lower limits. The probability density functions for the other characteristics show the same effects as function of locality definition as the one discussed here for the density.

### a.2 Influence of resolution

We compare the results from Via Lactea II (VL2) to a medium resolution simulation (VL2m) of the same halo. In the medium resolution run, the particle mass was a factor 64 higher than in the high resolution run. Therefore, in order to probe local properties with approximately the same number of particles spheres with times larger radii were used.

In Fig. 14 we show the probability density functions of the local density for the medium and high resolution run for different sizes of the spheres. In general the rarer peaks at the low and high end of the distribution are missing in the medium resolution run when one compares spheres with equal number of particles in both simulations (i.e. VL2 with 1 and VL2m with 4 ) since the low resolution run resolves less substructure and it is smoothed over a 4 times larger scale. When comparing spheres of equal physical size in both runs, then the low resolution run shows a broader distribution due to additional Poisson noise.

The lack of high density peaks in the low resolution simulation is due to numerical effects. Less subhaloes survive in the medium resolution run since the subhaloes are resolved with fewer particles and are more easily tidally disrupted. This leads to the effect that a larger region in the centre is smooth in lower resolution runs. The higher the resolution the smaller this apparently smooth region is.

Additionally, the higher degree of numerical artefacts (e.g. heating, relaxation etc.) in the medium resolution run smears out streams. For example, no dark matter streams are visible in the inner 50 kpc in phase space density maps of the medium resolution of Via Lactea II and the peak phase space densities of the subhaloes are much lower than the values from the high resolution run.

It is not clear at the moment to what degree the dark matter streams in the high resolution run are broadened by artificial heating in the simulations. Of course there are also real dynamical heating sources, such as dark matter subhaloes or baryonic structures like a stellar disk, a bulge, or gas clouds (see e.g. Ibata et al., 2002; Siegal-Gaskins & Valluri, 2008).

Material that ends up in the central regions of dark matter haloes originates from early forming halos (Diemand et al., 2005) and was accreted into the main halo early on (Helmi et al., 2002). One therefore expects a higher degree of phase mixing due to the early accretion and short time-scales in the centre. Nevertheless it is not clear if the smooth appearance of the central regions in Via Lactea II is entirely due to effcient phase mixing: The small, early progenitors that build up the central dark matter halo of Via Lactea II are under resolved and therefore too low numerical resolution might appear as an efficient phase mixing.

## Appendix B Balls-in-bins statistics

Let’s try to solve the following problem. We have bins and balls are thrown randomly into the bins. What is the probability of finding bins that contain balls after throwing balls?

The probability of a ball hitting a bin is given by and the probability of missing a bin is given by . Obviously, we have , for and for . After throwing balls we can write

(14) |

This is a recursion formula that allows us now to construct the general form of . By setting for all , we get for the non-zero values for the first few values of

(15) | |||||

(16) | |||||

(17) | |||||

(18) | |||||

(19) | |||||

(20) | |||||

(21) | |||||

(22) | |||||

(23) |

and so on. It becomes obvious that the general pattern is given by a binomial distribution

(24) | |||||

(25) |

We are specifically interested in the case of getting empty bins (or spheres in our case), i.e. . This results in

(26) |

where is the expectation value of balls per bin. We can simplify this in the large limit to

(27) |

which we use to estimate the expected fraction of empty bins respectively spheres.