Tracing the Dark Matter Sheet in Phase Space
MMon. Not. R. Astron. Soc. , 1–17 (2011) Printed 9 June 2018 (MN L A TEX style file v2.2)
Tracing the Dark Matter Sheet in Phase-Space
Tom Abel (cid:63) , Oliver Hahn † & Ralf Kaehler ‡ Kavli Institute for Particle Astrophysics and Cosmology, Stanford University,SLAC National Accelerator Laboratory, Menlo Park, CA 94025
MNRAS submitted
ABSTRACT
The primordial velocity dispersion of dark matter is small compared to the velocitiesattained during structure formation. The initial density distribution is close to uniformand it occupies an initial sheet in phase-space that is single valued in velocity space.Because of gravitational forces, this three dimensional manifold evolves in phase-spacewithout ever tearing, conserving phase-space volume and preserving the connectivity ofnearby points. N–body simulations already follow the motion of this sheet in phase-space.This fact can be used to extract full fine-grained phase-space structure information fromexisting cosmological N–body simulations. Particles are considered as the vertices of anunstructured three dimensional mesh moving in six dimensional phase-space. On thismesh, mass density and momentum are uniquely defined. We show how to obtain thespace density of the fluid, detect caustics, and count the number of streams as well astheir individual contributions to any point in configuration-space. We calculate the bulkvelocity, local velocity dispersions, and densities from the sheet – all without averagingover control volumes. This gives a wealth of new information about dark matter fluid flowwhich had previously been thought of as inaccessible to N–body simulations. We outlinehow this mapping may be used to create new accurate collisionless fluid simulation codesthat may be able to overcome the sparse sampling and unphysical two-body effects thatplague current N–body techniques.
Key words: cosmology: theory, dark matter, large-scale structure of Universe – galaxies:formation – methods: numerical
For the past 40 years, N–body simulations have allowed to nu-merically study the evolution of the distribution of matter inthe expanding Universe (cf. Peebles 1971; Bertschinger 1998;Springel et al. 2005). A significant number of simulation codeshave been developed for this purpose (e.g. Efstathiou et al.1985; Couchman 1991; Bryan & Norman 1997; Stadel 2001;Springel et al. 2001; Teyssier 2002; Wadsley et al. 2004, toname just a few). All such approaches to structure formationmodel the collisionless fluid of dark matter by a set of massiveparticles (typically of equal mass) and differ in how the gravi-tational forces are calculated at the positions of the particles.The forces are applied to update the velocities which in turnare used to update the positions. The system is then evolvedforward in time. From such simulations much has been learnedabout the formation and evolution of cosmological structuresand they have become a standard tool in physical cosmology.While three dimensional calculations have difficulty in sam-pling the six dimensional phase-space well (see e.g. Buchert &Bartelmann 1991) they have found a very large range of ap- (cid:63)
Email: [email protected] † Email: [email protected] ‡ Email: [email protected] plications and have driven much of the progress that has beenmade in the past decades of understanding structure forma-tion. In quite a range of these applications the space densityof the dark matter fluid is required and in many others thephase-space density is of great importance.Some open questions that require detailed informationabout the dark matter density and its velocity distributionare related to dark matter detection. For the indirect detectiontechniques, the predictions of the dark matter annihilation lu-minosity depend sensitively on the density of the dark matterstreams and the distribution and the relative velocities of theparticles. Assuming well mixed phase-space and assuming ashape of the velocity distribution function, this annihilationrate would scale with the square of the space density, ρ . Incurrent work, these estimates are typically carried out by fit-ting spherical profiles to the main dark matter halo and itssubhaloes and then assign annihilation luminosities by scalingthe square of the smoothed halo and subhalo profiles appro-priately. This smoothes the dark matter fluid sufficiently toavoid noisy estimates of the annihilation signal (Diemand et al.2008; Springel et al. 2008). However, in general the annihila-tion rates depend on the relative velocity of the dark matterparticles interacting. Here one can distinguish between darkmatter annihilation within individual streams of dark matteras well the contribution from stream-stream interactions which c (cid:13) a r X i v : . [ a s t r o - ph . C O ] J un T. Abel, O. Hahn & R. Kaehler differ strongly in the relative velocities (see e.g. Hogan 2001;Afshordi et al. 2009). This fine grained dark matter phase-space structure is equally important for considerations of darkmatter direct detection experiments where one applies veloc-ity cuts in the experimental analysis in order to reject certainbackgrounds. The annular modulation of the relative velocityand the main dark matter streams in the solar neighborhoodprovided by the earths motion around the sun allows one topotentially map the fine grained phase-space structure shouldthe experiments be able to detect dark matter.The best current approaches to probe phase-space struc-tures were surveyed by Maciejewski et al. (2009). These typ-ically start with some tessellation of phase-space such as aDelaunay triangulation or a Voronoi tesselation (Bernardeau& van de Weygaert 1996), or cartesian trees (Sharma & Stein-metz 2006; Ascasibar 2010). The mass of the particles foundwithin the cells give the local densities. When only the spacedensity is required, adaptive kernel smoothing is often em-ployed. In fact, most images of dark matter simulations shownare projections of kernel smoothed particle distributions. Con-figuration space density estimators also play a particularlyimportant role when studying the topology and character ofthe cosmic web (e.g Schaap & van de Weygaert 2000; Pelu-pessy et al. 2003; Arag´on-Calvo et al. 2007; Colberg et al.2008; Neyrinck 2008). However, in all cases, control volumesare defined such that they contain sufficient numbers of par-ticles to reduce sampling noise. Unfortunately, this will aver-age over large regions of configuration and phase-space andconsequently effectively degrade the spatial resolution of thecalculation.In any case, there is ample motivation to study the dis-tribution, evolution and current state of dark matter in theUniverse further by observations as well as simulations. In thiscontribution, we introduce a novel way to analyze N -body sim-ulations. Our approach naturally arises from considering howthe collisionless fluid of dark matter is expected to evolve inphase-space.As is well known, (e.g. Shandarin & Zeldovich 1989), atearly enough times, i.e. before shell crossing, the motion of thedark matter fluid is well described by the Zel’dovich approxi-mation (Zel’dovich 1970) as a potential flow x t = q + g t ∇ φ ( q ) , (1) v t = ˙ g t ∇ φ ( q ) . (2)Here φ is a potential field that is proportional to the initialgravitational potential field of perturbations, q are the initialparticle positions and g ( t ) is the growth factor of linear pertur-bations. At early times, i.e. for t →
0, also g →
0. These parti-cles occupy a three dimensional submanifold S of the entire sixdimensional phase-space with the time-dependent mapping: q (cid:55)→ ( q + g t ∇ φ ( q ) , ˙ g t ∇ φ ( q )) , (3)where for t → q (cid:55)→ ( q , ˙ g ∇ φ ( q )) , (4)the three-dimensional structure ( q ∈ R ) can be easily seen.The map between x t and q is bijective until shell crossingoccurs, i.e when more than one stream of dark matter existsat one spatial location. We will refer to this three dimensionalsubmanifold as the dark matter sheet (even when discussingit in one or two spatial dimensions). The volume of the sheetcontinues to grow as structure forms and evolves (Shandarin& Zeldovich 1989; Vogelsberger et al. 2008). At the same time, current N –body simulations of struc-ture formation do already follow individual dark matter par-ticles through phase-space (e.g. Bertschinger 1998). The N -body technique thus corresponds to sampling the sheet at afinite number of points q with the entire mass concentratedat their positions.Vogelsberger et al. (2008), White & Vogelsberger (2009)and Vogelsberger & White (2011) developed a powerful ap-proach to augment cosmological simulations to record moreknowledge about the evolution of this dark matter sheet. Theyderive an equation of motion for the distortion tensor aroundevery particle to linear order and then evolve it with every par-ticle during the simulation. They refer to this as the geodesicdeviation equation (GDE). This gives access to informationabout the evolution of the stream density along every parti-cle trajectory and allows to track the number of caustics aninfinitesimal fluid element surrounding a dark matter particlewill experience. This technique goes a long way in obtainingmore information about the fine grained phase-space structurein dark matter haloes (Vogelsberger & White 2011).More restricted calculations in this context have been car-ried out in fixed potentials (Stiff et al. 2001) or one dimensionsfor Newtonian gravity (Alard & Colombi 2005) and GeneralRelativity (Rasio et al. 1989). Also in the context of stellardynamics there is a large body of literature which exploresdetails of the phase-space structure of stellar system. From anumerical point of view the work of Cuperman et al. (1971)is particularly remarkable. Some 40 years ago these authorsrealized that in one dimension one can follow the phase-spaceboundary of a collisionless fluid and they give a beautiful im-plementation and calculations treating the phase-space fluidas a continuum. Studying the connection of their formalismto the N–body technique is revealing and what follows hereis in some ways the extension to three dimensions with theexception of our approach to velocity dimensions and the waythe Poisson equation is solved.We suggest that the three-dimensional manifold can bedecomposed by a space-filling grid that connects a finite num-ber of vertices q . The simplest version is to decompose thevolume into three dimensional simplices, i.e. tetrahedra, whichhave the nice topological property of being either convex ordegenerate. For any choice of a regular lattice of vertices q ,such a tetrahedral decomposition can be achieved by a Delau-nay triangulation. In contrast to the particle discretization, wecan now think of the dark matter mass being spread out overthe corresponding volume elements.This mesh traces the dark matter sheet as it subsequentlyevolves in phase-space. The motion of the mesh vertices areevolved using the Vlasov-Poisson equation of motion leadingto complex foldings of the submanifold (e.g. Arnold et al. 1982;Tremaine 1999). Any such folding is coincidental with a vol-ume inversion of a simplex. This volume inversion occurs whenthe simplex topologically evolves through a degenerate state(where the tetrahedron is planar because one vertex movesthrough the plane defined by the remaining three) which isequivalent to the emergence of a caustic. Note how this corre-sponds exactly to the sign changes of the distortion tensors ofVogelsberger et al. (2008) (their eq. 24).The motion of the vertices does not change the connec-tivity of the mesh so that at all times the simplex structurecan be constructed from knowledge of the q . For cosmological N -body simulations, there exists a unique mapping betweena particle ID and q , so that the phase-space structure of the c (cid:13) , 1–17 racing the dark matter sheet dark matter sheet can be reconstructed at all times. Projectingthe sheet onto configuration space gives then a volume fillingdensity field of the dark matter fluid that we propose to use asthe density field that should be used to solve Poisson’s equa-tion in future solvers for collisionless fluids. Current N -bodysolvers do not evolve the vertices consistently with a densityfield construed in the proposed way.As a first step towards this goal, we analyze the resultsof standard cosmological N -body simulations using this newdefinition of the dark matter sheet. The plan of the paperis as follows. First we will explain one and two dimensionalanalogues to introduce the relevant concepts. We then describethe details of our implementation before we apply the methodto analyze cosmological large-scale structure as well as thephase-space properties of a single galaxy cluster halo. The distribution function f ( x , p ; t ) describes the density of afluid in phase-space. It evolves via ∂f∂t = − p m · ∇ x f − ∇ x φ · ∇ p f, (5)where φ is the gravitational potential and m is the dark mat-ter particle mass. Fluid elements get stretched or compressedin coordinate space by advection p m · ∇ x f , and in the momen-tum coordinates by the gravitational forces p m · ∇ x f . Notethat in a Lagrangian frame, the first term on the right handside is zero. Furthermore, the second term describes how thefluid is stretched in momentum space and does not affectthe space density of the fluid parcel. This just states Liou-ville’s theorem (Gibbs 1902) that the volume in phase-spaceis conserved. Hence, any fluid volume (cid:52) x (cid:52) v will remain con-stant. We are interested here in the space density of the fluid,the projection of f into coordinate space. i.e. the integral ρ ( x ) = (cid:82) f ( x , v )d v . The contribution to the space densityof any stream of dark matter is only affected by the volumeit occupies in the space coordinates, i.e. (cid:52) x . Consequently,all that is necessary to follow the evolution of the dark matterdensity is to follow the Lagrangian evolution of fluid elements.The mass inside a volume element is conserved and its contri-bution to the space density of dark matter is described by thevolume it occupies. Conversely, for a given WIMP model oneknows the initial velocity dispersion at any point in space (e.g.Hogan 2001; Vogelsberger et al. 2008). Therefore, if one knowsthe spatial part of the phase-space density one has informa-tion about the density in velocity space. For a given shapeof the initial distribution function in the velocity directions(e.g. a Maxwellian) one has a reliable measure of the intrinsicvelocity density at all times.Using the Vlasov equation to describe DM is justified formost particle physics inspired models of dark matter. For aWIMP scenario with a 100 GeV particle e.g. there would be10 such particles in the Milky Way alone. So an elementin phase space that contains a million such WIMPS, giving awell defined phase space density, would have a spatial extentat mean density equivalent to a cube 500 meters on a side.Such a volume would only be a few meters on a side at theDM density expected in the solar neighbourhood. Clearly weare interested in scales much larger than this and the approx-imation of using a density in phase space is justified to a highdegree of confidence. ρ [ a r b i t r a r y ] neighbors7 point stencil−1−0.500.51 v x [ a r b i t r a r y ] Figure 1.
The one–dimensional plane wave collapse of Zel’dovich(Zel’dovich 1970; Binney 2004). The top panel gives the phase-space diagram showing the velocities of the particles at their loca-tions. The bottom panel gives the density of the dark matter insidethe stream, one computed with a seven point stencil (red squares),and the other computed from the volume between two neighboringpoints (solid line). Knowing the spatial volume between particlesalong one stream is sufficient to obtain accurate density estimatesat and between the points.
It is instructive to first describe a straightforward and wellknown one–dimensional example of the evolution of a collision-less fluid from which a number of lessons can be learned whichapply equally well in higher dimensions.
The phase-space diagram and the evolved density in aZel’dovich plane wave collapse is shown in Figure 1. The ini-tial sheet at very early times would be coincident with the x -axis as the initial velocity perturbation is small and the initialstate models a nearly homogeneous Universe. Sampling thisinitial state with particles of equal mass results in a grid ofuniformly placed particles. Their configuration space volumeis now simply related to their distances in the x –direction.Figure 1 shows the results of computing their local streamdensity from two approaches. In one, labelled “neighbour”, wetake the V i = x i +1 − x i as the volume between particle i and i + 1. One full particle mass is distributed in this volume andthe density at ( x i + x i +1 ) / ρ neigh = m p / | V | . Thevalues shown as “squares” in the same figure are computedincluding information from points further along the stream, ρ pt = 6 m p / | x i +3 − x i − | . It is defined at the particle position x i . A number of observations can be made. Volumes definedin this way may be positive or negative depending on whetherparticles have the same or opposite ordering that they hadinitially. Volume elements may also become zero. The den-sity involving more points along the stream gives rise to some c (cid:13)000
The phase-space diagram and the evolved density in aZel’dovich plane wave collapse is shown in Figure 1. The ini-tial sheet at very early times would be coincident with the x -axis as the initial velocity perturbation is small and the initialstate models a nearly homogeneous Universe. Sampling thisinitial state with particles of equal mass results in a grid ofuniformly placed particles. Their configuration space volumeis now simply related to their distances in the x –direction.Figure 1 shows the results of computing their local streamdensity from two approaches. In one, labelled “neighbour”, wetake the V i = x i +1 − x i as the volume between particle i and i + 1. One full particle mass is distributed in this volume andthe density at ( x i + x i +1 ) / ρ neigh = m p / | V | . Thevalues shown as “squares” in the same figure are computedincluding information from points further along the stream, ρ pt = 6 m p / | x i +3 − x i − | . It is defined at the particle position x i . A number of observations can be made. Volumes definedin this way may be positive or negative depending on whetherparticles have the same or opposite ordering that they hadinitially. Volume elements may also become zero. The den-sity involving more points along the stream gives rise to some c (cid:13)000 , 1–17 T. Abel, O. Hahn & R. Kaehler i i+1jj+1 oo oo oo oo oooooooo oo ooo ooo ooo ooo oo oo ooo oo oo ooo oo o ⦁ o a b ⦁⦁ c Figure 2.
There are a number of possible cell elements (and thatcan be further decomposed into simplices) one may consider as theunit cell. These unit cells can be used to tessellate the three dimen-sional dark matter sheet which connects the particles of an N–bodysimulation. smoothing and density extrema are clipped. The central highconfiguration space densities are reached for two reasons. Theprimordial stream densities along the sheet become larger andmany streams overlap adding their densities. The number ofstreams in space is always an odd number at any location inspace. Only at the caustics may one measure even numbers.The particle locations trace the sheet in phase-space. Anyunstructured space-filling grid that connects adjacent fluid ele-ments may be used to trace the dark matter sheet as it evolvesin phase-space. In fact, there is significant ambiguity here asillustrated in Figure 2. The two–dimensional analog shownthere has is based on triangles (the 2D simplex). The smallestpossible elements one may thus choose to follow would be theDelaunay triangulation of the points. However, these wouldgive two resolution elements per square initial cell (case c).It would seem unreasonable that the mass of the fluid wouldbe conserved exactly in each element, given that we only haveinformation at the vertices. This could only be true if themesh points were not distorted very much and the gradientsof the flow, both in configuration and momentum space hadlength scales much larger than the sides of the triangle. How-ever, cases a) depicted in the figure would likely be a betterchoice as a fundamental resolution element as its eight nodeson the surface would be able to more accurately describe thedeformations caused by the flow pattern. The mass inside thatboundary would be conserved to a better degree than choicesfor a fundamental resolution element with smaller area. Note,however, that we still can use triangles to calculate the volume(area in 2D) of the sheet.After considering some preliminaries in one and two di-mensions, let us proceed to the three–dimensional case.
The simplex in three dimensions is the tetrahedron. We willuse it to calculate volumes and tesselate our chosen funda-mental volume element. This is exactly equivalent to choosingline elements in one and triangles in two dimensions as the
310 265 74 x y z
Figure 3.
One of the possible decompositions of a cubical cell intotetrahedra. This choice results in 6 tetrahedra all of equal volumefraction. This particular case has two of the cube corners as verticesfor all tetrahedra and while the other corners are connected to twotetrahedra each. elements which are summed in volume calculations. Figure 3shows one of the choices we have employed to tesselate a cu-bical fundamental cell. The figure also gives a numbering ofvertices of the six tetrahedra which make up the cell. The con-nectivity of vertices is chosen such that in a regular uniformgrid all tetrahedron volumes are positive. For much of thecalculation we keep the sign, because, as we will see, it canbe a useful diagnostic of the flow. If we shift a tetrahedronsuch that one vertex coincides with the origin, the volumeis simply given by the determinant or, equivalently, from ascalar and a cross product involving its other three vertices, V tet = det | a , b , c | / a · ( b × c ) /
6. This implies that thevolume can be negative if the tetrahedron has been turnedinside out. This sign inversion occurs when one vertex movesthrough the plane defined by the other three vertices of thetetrahedron and is an efficient way to find caustics. Within anN-body calculation, tracking the number of volume inversionscould be used to trace the number of caustics crossed by a fluidelement. This volume can now be used straightforwardly to es-timate the local stream density of the fluid element describedby the tetrahedron. In our case, one tetrahedron contains onesixth of the mass of an N-body particle spread over its volumeso that ρ s = m P / V tet = m P | a · ( b × c ) | . (6)Alternatively one may choose a cubical region rather than asingle tetrahedron as the fundamental volume element to con-sider. This can be achieved e.g. by using the 24 (2 × × N -body simulation.This implies that the vertices are moving in a Lagrangian wayso that the spatial sampling is degrading over time in low-density regions and improving in high-density regions. Oneimplication of the Lagrangian motion of the mesh vertices is c (cid:13) , 1–17 racing the dark matter sheet Figure 4.
The need to resolve critical points of the flow: (left) ahalo-void-halo configuration that, when probed with a tetrahedralstructure leads to a final state of the structure as given in the mid-dle panel. The correct motion of the fluid element could have beenmore correctly been described by the right panel. The loss of spa-tial resolution can be circumvented by smoothing or by adaptiverefinement of the volume elements. that unresolved volume elements may not evolve accordingto the properties of the underlying flow. One particular suchexample is given in Figure 4 which illustrates that volumeelements covering a divergent flow can evolve as if they werein a convergent flow if only the vertices are convergent. Thisbehaviour implies that single volume elements, i.e. tetrahedra,of our space discretization are not to be trusted as perfect bagsof fluid. In the example from Figure 4, the single tetrahedronsuggests that shell-crossing occurred across a divergent flowwhich is unphysical – e.g. implying spurious links betweenhaloes across unresolved void regions in cosmological context.There are two possibilities to achieve a more robust den-sity estimate: (1) The volume estimates of several neighbour-ing tetrahedra are combined, so that density estimates arebased on a larger region of the flow, or (2) local refinementof the volume elements is performed whenever e.g. axis ratiosor curvature in phase-space suggest that the volume elementmight miss the flow properties. The second option is certainlythe more exciting possibility to achieve a density estimationmethod that is well consistent with flow properties. It requireshowever that the discretized dark matter sheet is refined whilethe N -body simulation is run, so that the newly inserted ver-tices are evolved with the flow. For this reason, we focus on anevaluation of the averaging approach in this paper and con-sider refinement strategies in a future publication. The choice of tetrahedra used for most of the plots in this pa-per were given with a slightly different tesselation than the oneshown in Figure 3. The connectivity list for the six tetrahedraspecifically is [4, 0, 3, 1], [7, 4, 3, 1], [ 7, 5, 4, 1], [7, 2, 5, 1], [7,3, 2, 1], [7, 6, 5, 2] where the unit cube vertices are labeled asin Figure 3. A natural way to obtain the dark matter streamdensity at the location of an N-body particle is to considerthe volume surrounding any vertex x as the union of all thetetrahedra that share this vertex. The volume of this unionis computed using the modulus of the volumes, i.e. invertedtetrahedra are not subtracted. Such a volume carries the massof four particles which allows us to estimate the local densityat a particle, ρ s,p = 4 m p / (cid:80) i =1 V itet . We call this local aver-age over adjacent tetrahedra the “primordial” stream densityat a particle’s position. It is a well defined quantity even af-ter shell-crossing. However, it is only defined for the vertices x of the tessellation. The simpler definition of equation (6),however, works very well and, from visual inspections of darkmatter renderings, appears perfectly justified. The next step in calculating a configuration space densityestimate, as well as in evaluating other stream properties is anintegration through velocity space. This is achieved by findingall intersections of tetrahedra with the point y at which thedensity or other properties are to be determined. We refer toall these intersections, which are not part of the primordial (orfundamental) stream, as the “secondary” streams. Because westart from a complete tessellation of the sheet, the number oftetrahedra enclosing any spatial point is the number of streamscontributing at this point.At the heart of a fast algorithm is thus a way to speed upthis search for point-tetrahedron-collisions. To find the phasespace properties e.g. at the locations of all the simulation par-ticles one at first sight expects to require 6 ∗ N searches. I.e.for each of the N –points check all 6 N tetrahedra whether theyoverlap that location. This indeed would be very numericallyinefficient. Fortunately, this can be radically improved usinga chaining mesh that organizes tetrahedra contained in meshcells or use tree structures which group tetrahedra into re-gions of space they occupy. One could use six trees, e.g., to or-ganize bounding boxes containing tetrahedra which is advan-tageous as they are generally smaller than the circumsphereof tetrahedra which would require less storage. We, however,implemented two other versions. One with a chaining meshand another that employs three bounding-box oct-trees offsetfrom each other containing disjoint sets of tetrahedra so thattetrahedra cutting along tree node boundaries only do so typi-cally for one of the trees. We then parallelized the search withOpenMP and arrived at a very practical tool to carry out thetessellation and measure the quantities we were interested in.For the most straightforward case where the stream den-sity is given by equation (6) the total density at a given loca-tions is now just the sum of these ρ s for all the intersectingtetrahedra. This simple approach is what we used for the ren-derings in the following section.Since y will typically lie inside a tetrahedron, we can alsointerpolate the primordial stream density to y if we definedthem by averaging over a cubical volume at the vertices of thistetrahedron as described above. A one-over-distance-squaredweighting (Shepard’s method) works well. Also in this case, thetotal configuration space density is simply obtained by findingall tetrahedra that contain the point under consideration andsumming over all their stream densities interpolated to thislocation. When evaluating velocities, we interpolate the veloc-ity field inside of a tetrahedron from the velocities of the fourparticles that constitute its vertices. Again, this is achieved byShepard’s method and no further averaging is needed.To be more explicit, let us emphasize here that all of whatfollows below is obtained solely by post-processing existingsimulations. No line of code has to be changed in the read-ers’ favourite cosmology code. As long as the code writes outthe particle IDs at every snapshot of the simulation, and theconnectivity of the initial particle distribution is known or canbe constructed, one can post-process simulations that alreadyexist and measure, visualise and analyse it in new ways.With these definitions and implementation details inhand, we now proceed to apply the method to a number ofN–body simulations of large scale structure formation. c (cid:13)000
The need to resolve critical points of the flow: (left) ahalo-void-halo configuration that, when probed with a tetrahedralstructure leads to a final state of the structure as given in the mid-dle panel. The correct motion of the fluid element could have beenmore correctly been described by the right panel. The loss of spa-tial resolution can be circumvented by smoothing or by adaptiverefinement of the volume elements. that unresolved volume elements may not evolve accordingto the properties of the underlying flow. One particular suchexample is given in Figure 4 which illustrates that volumeelements covering a divergent flow can evolve as if they werein a convergent flow if only the vertices are convergent. Thisbehaviour implies that single volume elements, i.e. tetrahedra,of our space discretization are not to be trusted as perfect bagsof fluid. In the example from Figure 4, the single tetrahedronsuggests that shell-crossing occurred across a divergent flowwhich is unphysical – e.g. implying spurious links betweenhaloes across unresolved void regions in cosmological context.There are two possibilities to achieve a more robust den-sity estimate: (1) The volume estimates of several neighbour-ing tetrahedra are combined, so that density estimates arebased on a larger region of the flow, or (2) local refinementof the volume elements is performed whenever e.g. axis ratiosor curvature in phase-space suggest that the volume elementmight miss the flow properties. The second option is certainlythe more exciting possibility to achieve a density estimationmethod that is well consistent with flow properties. It requireshowever that the discretized dark matter sheet is refined whilethe N -body simulation is run, so that the newly inserted ver-tices are evolved with the flow. For this reason, we focus on anevaluation of the averaging approach in this paper and con-sider refinement strategies in a future publication. The choice of tetrahedra used for most of the plots in this pa-per were given with a slightly different tesselation than the oneshown in Figure 3. The connectivity list for the six tetrahedraspecifically is [4, 0, 3, 1], [7, 4, 3, 1], [ 7, 5, 4, 1], [7, 2, 5, 1], [7,3, 2, 1], [7, 6, 5, 2] where the unit cube vertices are labeled asin Figure 3. A natural way to obtain the dark matter streamdensity at the location of an N-body particle is to considerthe volume surrounding any vertex x as the union of all thetetrahedra that share this vertex. The volume of this unionis computed using the modulus of the volumes, i.e. invertedtetrahedra are not subtracted. Such a volume carries the massof four particles which allows us to estimate the local densityat a particle, ρ s,p = 4 m p / (cid:80) i =1 V itet . We call this local aver-age over adjacent tetrahedra the “primordial” stream densityat a particle’s position. It is a well defined quantity even af-ter shell-crossing. However, it is only defined for the vertices x of the tessellation. The simpler definition of equation (6),however, works very well and, from visual inspections of darkmatter renderings, appears perfectly justified. The next step in calculating a configuration space densityestimate, as well as in evaluating other stream properties is anintegration through velocity space. This is achieved by findingall intersections of tetrahedra with the point y at which thedensity or other properties are to be determined. We refer toall these intersections, which are not part of the primordial (orfundamental) stream, as the “secondary” streams. Because westart from a complete tessellation of the sheet, the number oftetrahedra enclosing any spatial point is the number of streamscontributing at this point.At the heart of a fast algorithm is thus a way to speed upthis search for point-tetrahedron-collisions. To find the phasespace properties e.g. at the locations of all the simulation par-ticles one at first sight expects to require 6 ∗ N searches. I.e.for each of the N –points check all 6 N tetrahedra whether theyoverlap that location. This indeed would be very numericallyinefficient. Fortunately, this can be radically improved usinga chaining mesh that organizes tetrahedra contained in meshcells or use tree structures which group tetrahedra into re-gions of space they occupy. One could use six trees, e.g., to or-ganize bounding boxes containing tetrahedra which is advan-tageous as they are generally smaller than the circumsphereof tetrahedra which would require less storage. We, however,implemented two other versions. One with a chaining meshand another that employs three bounding-box oct-trees offsetfrom each other containing disjoint sets of tetrahedra so thattetrahedra cutting along tree node boundaries only do so typi-cally for one of the trees. We then parallelized the search withOpenMP and arrived at a very practical tool to carry out thetessellation and measure the quantities we were interested in.For the most straightforward case where the stream den-sity is given by equation (6) the total density at a given loca-tions is now just the sum of these ρ s for all the intersectingtetrahedra. This simple approach is what we used for the ren-derings in the following section.Since y will typically lie inside a tetrahedron, we can alsointerpolate the primordial stream density to y if we definedthem by averaging over a cubical volume at the vertices of thistetrahedron as described above. A one-over-distance-squaredweighting (Shepard’s method) works well. Also in this case, thetotal configuration space density is simply obtained by findingall tetrahedra that contain the point under consideration andsumming over all their stream densities interpolated to thislocation. When evaluating velocities, we interpolate the veloc-ity field inside of a tetrahedron from the velocities of the fourparticles that constitute its vertices. Again, this is achieved byShepard’s method and no further averaging is needed.To be more explicit, let us emphasize here that all of whatfollows below is obtained solely by post-processing existingsimulations. No line of code has to be changed in the read-ers’ favourite cosmology code. As long as the code writes outthe particle IDs at every snapshot of the simulation, and theconnectivity of the initial particle distribution is known or canbe constructed, one can post-process simulations that alreadyexist and measure, visualise and analyse it in new ways.With these definitions and implementation details inhand, we now proceed to apply the method to a number ofN–body simulations of large scale structure formation. c (cid:13)000 , 1–17 T. Abel, O. Hahn & R. Kaehler number particles m p /h − M (cid:12) (cid:15)/h − kpc32 . × . × . × . × Table 1.
The specifics of the suite of N -body simulations used inthis paper. All simulations are of a 40 h − Mpc cosmological volume, m p is the particle mass, (cid:15) the force softening. For our first applications, we chose to test the method on simu-lations of the same physical volume, differing only in numericalresolution. N -body Simulations We have carried out cosmological N -body simulations of avolume of 40 h − Mpc length, run with the tree-PM code
Gadget-2 (Springel 2005). The initial conditions for thesesingle-mass-resolution simulations were generated with the
Music code (Hahn & Abel 2011) keeping large-scale phasesidentical with changing mass and spatial resolution. We as-sume a concordance ΛCDM cosmological model with densityparameters Ω m = 0 . Λ = 0 . σ = 0 . H = 100 h kms − Mpc − with h = 0 .
703 and a spectral index n s = 0 . Our method allows to separate physically distinct structures.The number of streams at a given location can only ever bean odd number as any fold will add two more streams to anexisting one (e.g. Arnold et al. 1982; Shandarin & Zeldovich1989). The notable exception is at the location of causticswhere points may sit such as to only measure an odd numberof additional streams. For the number of streams defined at theparticle positions, we can use this fact to select the structuresthat are constituting the first caustic. We illustrate the mean-ing of the local number of streams in Figure 5. The particlesthat record that they are part of only their primordial streamclearly define the voids at the mass scale that is resolved inthe calculations. Particles that count two streams surround thesheets formed between voids. When they undergo their firstcaustic, they have already crossed through the sheet and areturning around on the side opposite to from where they fell infrom. The particles which measure three or more streams arealso shown and they trace the location of the collapsed objectswell. We still consider all particles which count two streams or
Figure 6.
The mass fraction distribution in streams. Many moreodd numbers of streams are found than even ones. We show them asseparate lines with the odd ones displaced down by a factor of tenfor clarity. They differ by a factor approximately 4 and that ratiodoes not depend much on the resolution of the underlying darkmatter simulation. An asymptotic slope of ∼ N − stream develops forthe higher resolutions for an intermediate number of streams. Atlow resolution steeper relations are inferred. For our most resolvedsimulations about 90 per cent of the mass is in collapsed structures.Some 50 per cent of the mass is in locations with 20 streams ormore. more as part of collapsed objects. This is the same decompo-sition that can be made in the GDE approach of Vogelsbergeret al. (2008) as shown for the environment of Aquarius haloesin Vera-Ciro et al. (2011) (their Figure B1) and Vogelsberger& White (2011) (their Figure 4).We can see the distributions of the mass fractions as afunction of the number of streams in Figure 6. We plot themfor particles recording odd and even counts separately. One c (cid:13) , 1–17 racing the dark matter sheet Figure 5.
Particles with different numbers of streams in a slice 0 . h − Mpc thick. Top left: particles whose primordial stream does notoverlap with other parts of the sheet. Top right: particles which are on their first caustic, i.e. they measure a number of streams of two attheir location. Bottom left: particles for which the number of streams is greater to or equal to three. Bottom right gives the average numberof streams on an infinitesimally thin slice. Distinct physical components become clearly visible and separated. may have expected that the fraction of particles that are oncaustics vs particles that have odd numbers of streams to de-crease, as the caustics are better resolved for the high resolu-tion runs. Instead, the offset between odd and even numberedmass fractions is approximately constant. This is just a featureof cold dark matter simulations that with increasing resolutionalso more smaller objects can be resolved. This is illustrateddirectly in the bottom panel of that Figure. The cumulativedistribution of mass above a given number of streams clearlydoes not converge. This just reflects that there are many smallscale density fluctuations that collapse even earlier when thesimulations can resolve them. It is quite plausible that the to-tal fraction of collapsed mass will ultimately approach the veryhigh value of 99 per cent that has been estimated analyticallyby Shen et al. (2006) from the ellipsoidal collapse model. For the volume averaged fraction as a function of streamsand the corresponding cumulative distribution in Figure 7,we observe a similar lack of convergence. Smaller and smallerpancakes are resolved as the resolution increases, making moreand more volume have had shell crossing in the past.It is clear from these results that questions about theshell-crossed mass and volume fractions in cold dark mattersimulations are intimately tied to a scale. Only when introduc-ing such a scale, e.g. through filtering of density perturbationsor a constant force softening across resolutions, we could hopeto obtain a meaningful measure of these quantities.This is compatible with previous results on mass and vol-ume fractions in the various parts of the cosmic web e.g. byHahn et al. (2007) using a fixed scale, or by Arag´on-Calvoet al. (2010) using adaptive filtering. In both cases the filter- c (cid:13)000
Particles with different numbers of streams in a slice 0 . h − Mpc thick. Top left: particles whose primordial stream does notoverlap with other parts of the sheet. Top right: particles which are on their first caustic, i.e. they measure a number of streams of two attheir location. Bottom left: particles for which the number of streams is greater to or equal to three. Bottom right gives the average numberof streams on an infinitesimally thin slice. Distinct physical components become clearly visible and separated. may have expected that the fraction of particles that are oncaustics vs particles that have odd numbers of streams to de-crease, as the caustics are better resolved for the high resolu-tion runs. Instead, the offset between odd and even numberedmass fractions is approximately constant. This is just a featureof cold dark matter simulations that with increasing resolutionalso more smaller objects can be resolved. This is illustrateddirectly in the bottom panel of that Figure. The cumulativedistribution of mass above a given number of streams clearlydoes not converge. This just reflects that there are many smallscale density fluctuations that collapse even earlier when thesimulations can resolve them. It is quite plausible that the to-tal fraction of collapsed mass will ultimately approach the veryhigh value of 99 per cent that has been estimated analyticallyby Shen et al. (2006) from the ellipsoidal collapse model. For the volume averaged fraction as a function of streamsand the corresponding cumulative distribution in Figure 7,we observe a similar lack of convergence. Smaller and smallerpancakes are resolved as the resolution increases, making moreand more volume have had shell crossing in the past.It is clear from these results that questions about theshell-crossed mass and volume fractions in cold dark mattersimulations are intimately tied to a scale. Only when introduc-ing such a scale, e.g. through filtering of density perturbationsor a constant force softening across resolutions, we could hopeto obtain a meaningful measure of these quantities.This is compatible with previous results on mass and vol-ume fractions in the various parts of the cosmic web e.g. byHahn et al. (2007) using a fixed scale, or by Arag´on-Calvoet al. (2010) using adaptive filtering. In both cases the filter- c (cid:13)000 , 1–17 T. Abel, O. Hahn & R. Kaehler ∝ N s t r ea m - V o l u m e f r ac t i o n −5 −4 stream V F ( > N s t r ea m ) Figure 7.
The volume fraction distribution in streams; a resolutionstudy. For our most resolved simulations about 85% of the volume isin voids, around 7 (14) per cent is in collapsed structures ( N streams larger or equal to three) for the 128 (256 ) run. As expected inCDM, the fraction of volume occupied by collapsed objects doesnot converge. ing scales are related to the non-linear scales today which isthe relevant scale for much of galaxy formation. The method presented here is also an ideal tool to visualizethe data from current N–body simulations and thus to furtherhelp extracting physical insights from the calculations. Fig-ure 8 gives an example that we have obtained with a custom-written OpenGL based renderer of the tetrahedral mesh. Theprimordial stream densities are averaged over abutting tetra-hedra as described in the implementation section above. Thentetrahedra are projected, taking advantage of the OpenGLprimitives designed specifically for polygonal meshes. Detailsof this algorithm and variations as well as efficient implemen-tations of current graphics hardware will be given in a forth-coming publication (Kaehler, Hahn and Abel, in prep.).We compare this to the visual impression one obtains byplotting individual points of a calculation vs. our new densitydefinition in Figure 9. One can clearly see how filaments andsheets in and surrounding voids can be distinguished easily −4 −5 −4 p r i m o r d i a l s h ee t ρ / ρ total sheet ρ / ρ Figure 10.
Two dimensional mass-weighted histogram comparingthe density of the primordial dark matter sheet to the total sheetdensity in the 128 simulation. Much of the mass is contained atconfiguration space densities about ten times the mean. The primor-dial stream density also scatters very much but its median densityis close to the average density of the Universe. now. The visual impression is commensurate with the statis-tics we present next. Figure 10 gives the relation of the density in the primordialstream to the total density evaluated at the locations of allthe particles, i.e. a mass-weighted histogram. At low densitiesthese are identical as this material is traced by the originalsheet and no folding has occurred. There is an enormous scat-ter at higher densities which we can quantify further. Figure 11gives the mass-weighted density distribution for all the simu-lations we have analyzed. The top panel summarizes the indi-vidual contributions for the 256 simulation. The primordialstream density distribution peaks slightly below mean densitywhile the total mass weighted density distribution is at muchhigher densities. The reason is that all the streams not insidethe primordial stream contributing to the density at the loca-tion of the particle contribute much more to the total densityat high densities. That distribution is given by the green linein the top panel labelled as “secondary”. There is an appar-ent power law part in the primordial stream densities visible.We will discuss this further when considering volume-weighteddistributions. These are given in Figure 12 where we show thevolume-weighted dark matter density. The total densities weestimate with our method are labelled as “Sheet” . We alsoindicate the resolution of the dark matter simulation used tocompute it. The median of the stream density is 1 . c (cid:13) , 1–17 racing the dark matter sheet Figure 8.
A rendering of the projected dark matter density in the 256 run using our density estimator and our custom GPU based renderer. Figure 9.
Comparison of the visual appearance of renderings of the dark matter density in the 256 run using our new density estimatorwith a simpler density estimate based on the log of the number of dark matter particles falling within given image pixels. While many of thewell sampled regions are clearly apparent in both, the detailed structure of filaments, sheets and how they connect to voids becomes onlyapparent in our new approach shown on the right.c (cid:13)000
Comparison of the visual appearance of renderings of the dark matter density in the 256 run using our new density estimatorwith a simpler density estimate based on the log of the number of dark matter particles falling within given image pixels. While many of thewell sampled regions are clearly apparent in both, the detailed structure of filaments, sheets and how they connect to voids becomes onlyapparent in our new approach shown on the right.c (cid:13)000 , 1–17 T. Abel, O. Hahn & R. Kaehler totalprimordialsecondaryvoronoi d M / d l og ρ −4 −3 total256 d M / d l og ρ −4 −3 primordial ∝ ρ - . d M / d l og ρ −4 −3 secondary256 ρ / ρ d M / d l og ρ −4 −3 Figure 11.
Mass-weighted density distributions. The top left panelshows the histogram for four different densities defined at the parti-cle locations for the 256 run. The density estimated from a zerothorder estimate of the Voronoi tessellation, the total sheet dark mat-ter density, the primordial stream density and the secondary streamdensity. The results of the resolution study is given in the otherthree panels. The total space density is given in the top right panel.The bottom left is the mass-weighted density pdf of sheet densityin the primordial stream. The bottom right panel gives the densitycontributed by material that is not in the primordial stream. ∝ ρ - d V / d l og ρ −4 −3 SheetVoronoiSheet 64 Sheet 128 Sheet 256 ∝ ρ - . ρ / ρ d V / d l og ρ −4 −3 Figure 12.
Volume-weighted density distributions. The top panelshows the histograms for the 256 run, the lower those for the 32 run. The zeroth-order density estimated from a Voronoi tessellationis shown with a dashed line, the total sheet dark matter densitywith a solid line. At both resolutions, both the Voronoi and thestream density approach a ρ − power-law at high densities. Also,the two methods produce different estimates at intermediate densi-ties of ρ/ ¯ ρ ∼
10. The bottom panel also shows in grey the densityhistograms from our method for all simulations to aid the compar-ison. the unique Voronoi cells around each particle. The density inthat volume is then simply defined as the mass of the en-closed particle divided by that volume element. Following vande Weygaert & Schaap (2009), we refer to this as the zeroth–order Voronoi density. Albeit, this is clearly more noisy thanthe DTFE density estimator developed by Schaap & van deWeygaert (2000) and Pelupessy et al. (2003) since the densityis defined for the smallest region. DTFE is more advancedand employs averaging over nearby tetrahedral cells. Theseauthors give comparisons of that estimator to the smooth ker-nel estimates obtained with the otherwise popular Smoothed-Particle Hydrodynamics estimator. Like DTFE, the simplezeroth–order Voronoi estimator we use tessellates the entirevolume and has no parameters and as such is a well suitedbenchmark for comparison.Quite strikingly, at low and high densities our density esti-mate is very similar to the zeroth–order Voronoi volume baseddensity estimator (Figure 13). Both methods do not convergeat the very low density tail when varying resolution. This is c (cid:13) , 1–17 racing the dark matter sheet −5 −4 V o r o n o i ρ / ρ total sheet ρ / ρ Figure 13.
Two dimensional histogram comparing the zeroth–order Voronoi density estimate vs. the total sheet density. The cor-respondence is quite good. The largest difference is observed forvalues between a third and thirty times the mean density of darkmatter. The zeroth–order Voronoi density estimators overestimatesthe volumes in regions around filaments and sheets. physical in that the simulations model smaller voids when theresolution is increased. These can achieve lower densities thantheir larger counterparts in lower resolution simulations. So,at higher resolutions, it is not surprising to see a tail grow-ing at those lowest densities. Also, the peak distribution shiftscontinuously to lower densities as more particles are employed.Significant differences in the volume fraction at a given den-sity of our method with the zeroth–order Voronoi method areseen at intermediate densities around the mean density. Thisis understandable as the volumes that the Voronoi tessella-tion provides, tend to connect particles in the voids to theparticles in the sheets and filaments. At that point it spreadsparticles in filaments into volumes that are larger and conse-quently estimates lower densities. In fact, these estimate arevery significantly different. So much so that, if we integrate thedensity from our method over the volumes computed from thezeroth–order Voronoi estimator, the total mass in the box isoverestimated by more than a factor of ten. It certainly wouldbe interesting to compare our density estimate not only tothis zeroth-order Voronoi density but also to other density es-timators such as DTFE, adaptive Kernel softening, etc.. Thisis, however, beyond the scope of this first exposition of ourapproach.Next, we will now apply our method to measuring a num-ber of well studied quantities in dark mater haloes.
Navarro et al. (1996) have discovered a universal radial pro-file of the dark matter density in virialized haloes. This is oneof the key findings of cosmological N–body simulations anda large body of literature has largely confirmed the finding.We will give these profiles next. To get the best possible es- timate, we chose 100,000 test positions and bin these in 50radii, spaced logarithmically in radius. Figure 14 summarizesour findings. The density profiles computed from the dark mat-ter sheet are somewhat shallower and have about 50% largercentral densities at all resolutions for the single halo we haveanalyzed. Physically, it is conceivable that volume elementsformed by particles on radial orbits oscillate such that thebounding regions have a higher probability to be found atlarge radii while still contributing to the density interior tosmall radii. Similarly, one can picture particles of the sheetorbiting the center at larger radii such that the volume ele-ment they span can contribute to the center. Our method hasa well defined density at all radii and it is bound to be a con-stant at the lowest of radii where one averages always overthe same tetrahedra. At large radii, the new density estimateand the Voronoi estimates all agree extraordinarily well. Thisis true even in the infall region. The Voronoi estimates at allradii are perfectly consistent with a simpler estimate based onthe particle mass binned in shells divided by the shell volume(not shown).The masses included within a radius do converge alsoquite well with our estimate being consistently slightly higherat all radii. This is understandable as mass from particles out-side a given radius can contribute if they span a volume ele-ment that has nodes inside the radius.The number of streams that contribute to the profile aregiven in the lower panel of Figure 14. Not surprisingly, theseincrease with increasing resolution. If one scales them by fac-tors of eight between increasing resolutions, some closer con-vergence is observed. Between the 128 and 256 simulations,there remains a difference of about 30% which is likely justdue to streams contributing from larger radii to the regionsinside the particles spanning the volume element.It certainly would be interesting that the dark matterdensity profile in the central parts of haloes could be differentthan one estimates by measuring dark matter particles inside agiven radius. It is also suggestive how well our density profilesconverge from 128 to 256 particles. However, if the massprofile were indeed different, also the forces contributing tothe particle motions would be changed. So even if our densityestimator were more accurate, one could not prove that theresult shown in Figure 14 is the correct physical one untilone has evolved the dark matter sheet consistently, i.e. usingaccelerations created by the density distribution of the actualsheet elements. We discuss some potential approaches to carryout such simulations in the discussion section. While our method gives access to the full fine-grained phase-space structure, we chose to only show moments to serve asan example of what the method is good for, and to be ableto compare to work done on this previously. While the colli-sionless fluid does not experience microphysical collisions, thescattering provided by the time-varying gravitational poten-tial leads to mixing in phase-space. The velocity dispersion ofthe particles is a measure of the effective pressure of the darkmatter, which is of relevance for understanding the dynamicalstructure of orbits, i.e. e.g. the expectations of how observ-able stars move in the DM potential. Figure 15 summarizesthe radial profiles of the velocity dispersion for the same halowe have analyzed above for density profiles. We again draw c (cid:13)000
Navarro et al. (1996) have discovered a universal radial pro-file of the dark matter density in virialized haloes. This is oneof the key findings of cosmological N–body simulations anda large body of literature has largely confirmed the finding.We will give these profiles next. To get the best possible es- timate, we chose 100,000 test positions and bin these in 50radii, spaced logarithmically in radius. Figure 14 summarizesour findings. The density profiles computed from the dark mat-ter sheet are somewhat shallower and have about 50% largercentral densities at all resolutions for the single halo we haveanalyzed. Physically, it is conceivable that volume elementsformed by particles on radial orbits oscillate such that thebounding regions have a higher probability to be found atlarge radii while still contributing to the density interior tosmall radii. Similarly, one can picture particles of the sheetorbiting the center at larger radii such that the volume ele-ment they span can contribute to the center. Our method hasa well defined density at all radii and it is bound to be a con-stant at the lowest of radii where one averages always overthe same tetrahedra. At large radii, the new density estimateand the Voronoi estimates all agree extraordinarily well. Thisis true even in the infall region. The Voronoi estimates at allradii are perfectly consistent with a simpler estimate based onthe particle mass binned in shells divided by the shell volume(not shown).The masses included within a radius do converge alsoquite well with our estimate being consistently slightly higherat all radii. This is understandable as mass from particles out-side a given radius can contribute if they span a volume ele-ment that has nodes inside the radius.The number of streams that contribute to the profile aregiven in the lower panel of Figure 14. Not surprisingly, theseincrease with increasing resolution. If one scales them by fac-tors of eight between increasing resolutions, some closer con-vergence is observed. Between the 128 and 256 simulations,there remains a difference of about 30% which is likely justdue to streams contributing from larger radii to the regionsinside the particles spanning the volume element.It certainly would be interesting that the dark matterdensity profile in the central parts of haloes could be differentthan one estimates by measuring dark matter particles inside agiven radius. It is also suggestive how well our density profilesconverge from 128 to 256 particles. However, if the massprofile were indeed different, also the forces contributing tothe particle motions would be changed. So even if our densityestimator were more accurate, one could not prove that theresult shown in Figure 14 is the correct physical one untilone has evolved the dark matter sheet consistently, i.e. usingaccelerations created by the density distribution of the actualsheet elements. We discuss some potential approaches to carryout such simulations in the discussion section. While our method gives access to the full fine-grained phase-space structure, we chose to only show moments to serve asan example of what the method is good for, and to be ableto compare to work done on this previously. While the colli-sionless fluid does not experience microphysical collisions, thescattering provided by the time-varying gravitational poten-tial leads to mixing in phase-space. The velocity dispersion ofthe particles is a measure of the effective pressure of the darkmatter, which is of relevance for understanding the dynamicalstructure of orbits, i.e. e.g. the expectations of how observ-able stars move in the DM potential. Figure 15 summarizesthe radial profiles of the velocity dispersion for the same halowe have analyzed above for density profiles. We again draw c (cid:13)000 , 1–17 T. Abel, O. Hahn & R. Kaehler voronoi 256 voronoi 128 voronoi 64 ρ / ρ enclosed 256 enclosed 128 enclosed 64 M ( r ) / h - M ๏ r / h -1 Mpc s t r ea m s Figure 14.
Dark Matter density profile in the most massive halo of2 × M (cid:12) with R vir ≈ particles at random positions in spherical shells for which wemeasure the stream-density-weighted bulk velocity, subtractit from the stream-density-weighted local velocity dispersion,before we finally average it to obtain the velocity dispersionin radial shells.These velocity dispersions differ quite significantly fromthe similarly termed quantities presented previously (e.g.Navarro et al. 2010, and references therein) . While this mayseem surprising, one has to keep in mind that we measure thedispersion at a single point, i.e. we do not carry out any aver-aging over volume. Dispersions quoted in the past measured a combination of a bulk and local dispersion. This will not onlyhave large sampling error but also confuse turbulent bulk mo-tions with actual microphysical velocity distributions. Indeed,we can see that our measured velocity dispersion does notconverge at scales of about one half the virial radius, as onlyabout one thousand streams contribute there for our highestresolution results. The distribution functions at this locationwill be quite anisotropic and a single temperature will be abad fit (see below). The halo is remarkably cold in the cen-ter – having less than half the velocity dispersion expectedfrom the virial velocity. In the same figure, we also show thepseudo phase-space density of Taylor & Navarro (2001) whichhas been found to be a perfect power-law entirely indepen-dent of resolution (Navarro et al. 2010). When we measurethe average of only the fine-grained quantities, as shown inFigure 15, this perfect power-law disappears. This may sug-gest that much of the measure is dominated by large-scale bulkflows. It is worthwhile to explore this further with higher res-olution simulations, where one can more confidently separatethermal motions from the bulk velocity dispersion.It is though just as interesting to check the actual dis-tribution function of dark matter velocities at a given point.The seminal work of Lynden-Bell (1967) discussed this in thecontext of stellar systems. The global distribution has beenmeasured from simulations many times (see e.g. Hoeft et al.2004; Wise & Abel 2007; Navarro et al. 2010; Vogelsbergeret al. 2009) but, to the best of our knowledge, this was neverdone at individual points in the simulations. Figure 16 sum-marizes the distributions found at three different points in themost massive halo. We show it at the center where a relativelyhot component overlays a colder one. At the center, the distri-butions of the individual velocity components have peaks thatalmost coincide and widths which are quite similar as well.They are not too far from an isotropic Maxwell-Boltzmanndistribution in their cores. As we step out in radius, the sit-uation changes rapidly and the microphysical flow structureclearly shows more and more anisotropy. Interestingly, a quitehot component is seen along the x direction. At the same thevelocity distribution in the y direction is the coldest at all radii.These distributions are consistent with the visual impressionobtained from the velocity dispersion slice in Figures 17 and 15which also shows that some of the hottest DM fluid elementsare found just inside the virial radius.It is remarkable how much physics can be learned fromeven these low resolution simulations analyzed here. For thehalo we just discussed, there are not even 600,000 dark matterparticles inside the virial radius of 1 . h − Mpc. At the sametime, there are already enough streams to compute meaningfulmeasures of the structure of phase-space. We are certainlylooking forward to carrying out a more detailed analysis onhigher resolution simulations. This point is born out by thevisual impression given by infinitesimal slices as shown in 17,which we will describe next.
To aid in the interpretation of the profiles we have just pre-sented, we also give two dimensional slices through the darkmatter density and “entropy”, which we define analogously tothe adiabats used when studying e.g. galaxy clusters hydro-dynamically simply as S DM = σ / ( ρ/ ¯ ρ ) / , which then hasunits of the square of a velocity. Also, the average numberof streams contributing to every point on the slice is given. c (cid:13) , 1–17 racing the dark matter sheet ρ / ρ h - M p c −3−2−10123 σ −1 0 1 2 3 log σ / ( ρ / ρ ) h -1 Mpc −3 −2 −1 0 1 2 3 h - M p c −3−2−10123 h -1 Mpc −3 −2 −1 0 1 2 3
Figure 17.
Infinitesimally thin slice through the 256 simulations for the most massive halo. We show the density in units of the meandensity (top left), the stream-density weighted velocity dispersion in kilometers per second (top right), the dark matter entropy [(km/s) ],computed from the density and stream weighted velocity dispersion σ / ( ρ/ ¯ ρ ) / , and the average number of streams (not stream densityweighted) at the bottom right. Clearly, our approach gives information so far thought inaccessible from current simulations. Material from the voids falls in perfectly cold. We can thinkof the velocity dispersion as a measure of the temperature ofthe fluid. It is ill defined in the single stream regions fallingfrom the voids. However, these carry very little mass. Then wesee a region that extends to about two Mpc from the centerwhich hosts multi-stream material of the order of about tenstreams. The virial radius, which is approximately at one Mpc,shows a marked increase in the velocity dispersion and a muchsmoother density structure. Even on this scale, we can see thecold central isothermal part of the object, both in the velocitydispersion and in the entropy. Substructure is easily seen as cold low entropy material embedded in the hot halo. Many ofthe structures seen here are very reminiscent of adiabatic hy-drodynamical simulations of galaxy clusters (e.g. Frenk et al.1999) and even first star formation (Abel et al. 2002) wheregas enters haloes predominantly through filaments and shockheats, resulting in a halo with rising entropy profiles with ra-dius. c (cid:13)000
Infinitesimally thin slice through the 256 simulations for the most massive halo. We show the density in units of the meandensity (top left), the stream-density weighted velocity dispersion in kilometers per second (top right), the dark matter entropy [(km/s) ],computed from the density and stream weighted velocity dispersion σ / ( ρ/ ¯ ρ ) / , and the average number of streams (not stream densityweighted) at the bottom right. Clearly, our approach gives information so far thought inaccessible from current simulations. Material from the voids falls in perfectly cold. We can thinkof the velocity dispersion as a measure of the temperature ofthe fluid. It is ill defined in the single stream regions fallingfrom the voids. However, these carry very little mass. Then wesee a region that extends to about two Mpc from the centerwhich hosts multi-stream material of the order of about tenstreams. The virial radius, which is approximately at one Mpc,shows a marked increase in the velocity dispersion and a muchsmoother density structure. Even on this scale, we can see thecold central isothermal part of the object, both in the velocitydispersion and in the entropy. Substructure is easily seen as cold low entropy material embedded in the hot halo. Many ofthe structures seen here are very reminiscent of adiabatic hy-drodynamical simulations of galaxy clusters (e.g. Frenk et al.1999) and even first star formation (Abel et al. 2002) wheregas enters haloes predominantly through filaments and shockheats, resulting in a halo with rising entropy profiles with ra-dius. c (cid:13)000 , 1–17 T. Abel, O. Hahn & R. Kaehler σ / k m s - l og σ / ( ρ / ρ ) / −1012345 r / h -1 Mpc ( ρ / ρ ) / σ −7 −6 −5 −3 Figure 15.
Radial spherically averaged profiles of the velocity dis-persion (top), the dark matter “entropy” σ / ( ρ/ ¯ ρ ) / (middle) andthe pseudo phase-space density (bottom) for the same halo as inFigure 14. The velocity dispersion is remarkably flat inside aboutone tenth of the virial radius. The dark matter “entropy” profilealso shows signs of already converging at the modest resolutionsemployed here. Using the microscopic velocity dispersion of our ap-proach which removes the bulk flows does not give the typical pow-erlaw behavior in the pseudo-phase-space density found when usingthe total dispersion of particles at those radii. There are large advantages to have well-defined grids whichallows one to interpolate to any point in space. This is a veryobvious observation of course, it is, however, a large step for-ward in understanding dark matter simulations. This has ledto a number of approaches being devised that allow such inter-polation, such as the methods discussed in the introduction.Here we discuss but a few approaches on how to use the tes-sellated dark matter sheet for interpolation.When probing the sheet at the particle locations, we find center d f ( v ) / d v −5 v x v y v z -1 Mpc d f ( v ) / d v −5 -1 Mpc|v| / km/s d f ( v ) / d v −5 v / km/s −1000 0 500 1000 Figure 16.
The velocity distribution function binned in velocityspace at three points at varying distance from the center of thehalo (left panels) and the individual distributions of the velocitycomponents (right). We used 200, 80 and 30 bins for the points atthe center, 0.4 h − Mpc and 1 h − Mpc from the center along oneaxis. These had 150,000, 13,000 and 400 individual streams thatcontributed. the primordial stream densities, total space densities, numberof streams, velocities etc.. Hence, we have sampled the fullvolume and have a non-uniform distribution of the fields weaim to interpolate. One may choose to achieve further interpo-lation by using a distance-weighted estimate from the nearestparticle locations. An efficient way to find two dimensionalslices is to take all tetrahedron edges and compute their in-tersections with the plane to be interpolated to. Along everyedge one can now linearly interpolate the values of the verticesto the plane. The resulting scattered data on the plane is thentriangulated again and interpolated with linear interpolationbetween nodes. As an example the slice of the total sheet den-sity is shown in Figure 18 which also gives a visual clue to howcloud in cell interpolation would sample the density field.Similarly, this allows us to extract one dimensional skew-ers at arbitrary resolution from N–body simulations. As anexample, Figure 19 compares the velocity along a random linethrough the volume for different resolutions. The large scalemodes are all consistent by design. It is interesting to see thatconvergence is quite slow and suggests to extend this analysisrigorously to much higher resolutions.
Vogelsberger et al. (2008), White & Vogelsberger (2009) andVogelsberger & White (2011) developed the GDE formalismto allow a calculation of the primordial stream density. Theirapproach modifies the simulation code to integrate an evo-lution equation of the tidal tensor along with every particletrajectory. In principle, this can be much more accurate sincethe local stretching of the dark matter sheet is calculated at c (cid:13) , 1–17 racing the dark matter sheet Figure 18.
The logarithm of the density in an infinitesimally thinslice in units of the mean density for the 128 simulation. The whitedots show the location of the particles which would contribute to acloud in cell interpolation on a grid with cells as large as the meanparticle spacing. The squares at the top left show the area to whichthese individual particles would contribute to at their locations. −1 Mpc]−100−75−50−250255075100 v x [ k m / s ] Figure 19.
The velocity field along a one dimensional line extractedusing tetrahedron edges to interpolate to a slice plane. The differ-ences in resolution are understandably quite large given the largerange of mass resolution of the simulations. The large scale fea-tures remain recognizable even at the lowest of resolutions. The“N-shaped” infall regions are seen for many structures. every time step of the calculation. It will be of great interestto compare our approach to theirs in detail. This will requireto carry out the calculations with both methods on an identi-cal simulation to facilitate a particle by particle comparison.This should be particularly interesting given that our methodcan, in addition to the primordial stream density, also pro-vide the total space density and number of streams at every location. Since our approach also gives full fine-grained phase-space information, it seems plausible one could combine bothapproaches to a hybrid which inherits the advantages of both.More generally, both the GDE and our approach sug-gest a number of possible approaches to improve the accu-racy of N–body calculations. Almost all current cosmologicalN–body solvers employ the particle mesh method at least forthe largest scales in the calculation. The cloud in cell approx-imation is used to interpolate the dark matter particles to agrid on which the gravitational potential will be evaluated be-fore differencing it to obtain the gravitational forces on theparticles. Since one integrates twice to get the potential fromthe density field and only differentiates once, this method givesreasonably smooth gravitational forces. However, it inherentlymodels a very noisy inaccurate density distribution obtainedfrom CIC which will have the largest relative errors in poorlysampled regions such as voids, pancakes and filaments (seeFigure 18).We have shown that our density estimator would have sig-nificantly more fidelity and reliability for these large regions.It is in principle quite straightforward to modify an existingparticle mesh code to make use of our density estimator andthen derive more accurate potentials and forces from it. It onlyinvolves the interpolation step to the grid. When interpolatingthe contribution back to the particle positions one could makeuse of the known analytical solutions to the Newtonian poten-tial of homogeneous polyhedra (Waldvogel 1979). Similarly,these analytic formulae could be applied in direct summationand tree-based codes. A priori it may seem difficult to imag-ine how to construct trees efficiently when considering thatthe tetrahedra may become exceedingly distorted and elon-gated and would cover many nodes of the tree. However, anynew code would most likely ever only be employed using localmesh refinement given that the tessellation we suggest givesmany opportunities to discover the regions of the flow whichmay be prone to errors. The local curvature of the flow com-pared to the tetrahedra edges is one measure but also theaxis ratio of individual tetrahedra provides an estimate wherethe flow would benefit from refinement. The key to any suchnew method will have to be to fully consider the dark matteras a fluid so that spurious particle-particle interactions maybe avoided and multi-mass resolution becomes feasible. Givena locally refined mesh, tree structures will remain useful inrapidly finding neighbours and retain n log n scaling.There are a number of improvements possible that willhelp to develop our GPU assisted volume rendering further.Using vertex values and some form of Shepard’s method tocarry out distance dependent weighting should still likely bevery fast on current GPUs even when drawing billions of tetra-hedra.Higher order interpolation in fact could be another av-enue to improve on the method suggested here. We have onlyimplemented the very simplest of ideas. Namely that volumeelements in phase-space are uniformly filled with the darkmatter fluid. This is similar in spirit to donor cell methodsused decades ago for hydrodynamical flows. We believe that itwill be possible to improve on our approach significantly. Up-winded schemes with linear reconstruction were a large gainin accuracy in numerical fluid dynamics and similar improve-ments are certainly possible here.Given that one now has a natural grid that can be used tointerpolate any state variables as well as the full fine-grainedphase-space structure, one can also define differentials on it. c (cid:13)000
The velocity field along a one dimensional line extractedusing tetrahedron edges to interpolate to a slice plane. The differ-ences in resolution are understandably quite large given the largerange of mass resolution of the simulations. The large scale fea-tures remain recognizable even at the lowest of resolutions. The“N-shaped” infall regions are seen for many structures. every time step of the calculation. It will be of great interestto compare our approach to theirs in detail. This will requireto carry out the calculations with both methods on an identi-cal simulation to facilitate a particle by particle comparison.This should be particularly interesting given that our methodcan, in addition to the primordial stream density, also pro-vide the total space density and number of streams at every location. Since our approach also gives full fine-grained phase-space information, it seems plausible one could combine bothapproaches to a hybrid which inherits the advantages of both.More generally, both the GDE and our approach sug-gest a number of possible approaches to improve the accu-racy of N–body calculations. Almost all current cosmologicalN–body solvers employ the particle mesh method at least forthe largest scales in the calculation. The cloud in cell approx-imation is used to interpolate the dark matter particles to agrid on which the gravitational potential will be evaluated be-fore differencing it to obtain the gravitational forces on theparticles. Since one integrates twice to get the potential fromthe density field and only differentiates once, this method givesreasonably smooth gravitational forces. However, it inherentlymodels a very noisy inaccurate density distribution obtainedfrom CIC which will have the largest relative errors in poorlysampled regions such as voids, pancakes and filaments (seeFigure 18).We have shown that our density estimator would have sig-nificantly more fidelity and reliability for these large regions.It is in principle quite straightforward to modify an existingparticle mesh code to make use of our density estimator andthen derive more accurate potentials and forces from it. It onlyinvolves the interpolation step to the grid. When interpolatingthe contribution back to the particle positions one could makeuse of the known analytical solutions to the Newtonian poten-tial of homogeneous polyhedra (Waldvogel 1979). Similarly,these analytic formulae could be applied in direct summationand tree-based codes. A priori it may seem difficult to imag-ine how to construct trees efficiently when considering thatthe tetrahedra may become exceedingly distorted and elon-gated and would cover many nodes of the tree. However, anynew code would most likely ever only be employed using localmesh refinement given that the tessellation we suggest givesmany opportunities to discover the regions of the flow whichmay be prone to errors. The local curvature of the flow com-pared to the tetrahedra edges is one measure but also theaxis ratio of individual tetrahedra provides an estimate wherethe flow would benefit from refinement. The key to any suchnew method will have to be to fully consider the dark matteras a fluid so that spurious particle-particle interactions maybe avoided and multi-mass resolution becomes feasible. Givena locally refined mesh, tree structures will remain useful inrapidly finding neighbours and retain n log n scaling.There are a number of improvements possible that willhelp to develop our GPU assisted volume rendering further.Using vertex values and some form of Shepard’s method tocarry out distance dependent weighting should still likely bevery fast on current GPUs even when drawing billions of tetra-hedra.Higher order interpolation in fact could be another av-enue to improve on the method suggested here. We have onlyimplemented the very simplest of ideas. Namely that volumeelements in phase-space are uniformly filled with the darkmatter fluid. This is similar in spirit to donor cell methodsused decades ago for hydrodynamical flows. We believe that itwill be possible to improve on our approach significantly. Up-winded schemes with linear reconstruction were a large gainin accuracy in numerical fluid dynamics and similar improve-ments are certainly possible here.Given that one now has a natural grid that can be used tointerpolate any state variables as well as the full fine-grainedphase-space structure, one can also define differentials on it. c (cid:13)000 , 1–17 T. Abel, O. Hahn & R. Kaehler
Consequently, it becomes possible to study vorticity, diver-gences as well as carry out the Cauchy-Stokes decompositionof the dark matter velocity fields. This way one can separatebulk, shear and rotational components of the velocity fieldswhich undoubtedly will make it possible to track down thephysical origin of dark matter density profiles as well as to bet-ter understand the internal structure of haloes and the cosmicweb.There is a remarkably large number of applications wherewe think our method can aid to gain new insights. Whether itis gravitational lensing to find more accurate lensing potentialsto studying the origin of the angular momentum profiles(seee.g. Fig. (12) in Bullock et al. 2001). Obviously the connec-tion between dark matter and the baryons they host can beexplored much more fully now as well.As we were preparing this manuscript, Shandarin et al.(2011) posted a paper on the electronic preprint server whichexplores the same basic idea as the one we present here. Theconcept of tessellating phase-space and tracking the dark mat-ter sheet is identical to ours. Details of the implementationand what to think of as fundamental parts of the approachare not the same. Their choice of tessellation is quite differ-ent. They pick the minimal combination of tetrahedra of theunit cube possible which has five elements where one of themis twice the volume of the other four and itself does not tesse-late the space uniformly. Consequently they alternately rotateadjacent cubes such that the edges of tetrahedra never cross.This is effective albeit likely more cumbersome for a prac-tical implementation. The powerlaw Shandarin et al. (2011)gives for the volume fraction as a function of the number ofstreams ( f V ( N stream ) = 0 . N − . ± . stream ) is to be comparedwith our Figure 7. Their power-law fits approximately our 32 run and likely just reflects the fact that the single simulationthey study had approximately two times worse mass resolu-tion than our 32 run with an effective gravitational softeninglength about five times larger than ours. So both approachesdo agree. We at this time would not attach a special meaningto this power-law as it clearly is strongly resolution depen-dent with our highest resolution run giving something close to N − stream . Our description also discusses the fine-grained struc-ture in the velocity directions of phase space, discusses haloproperties and profiles and gives visualizations of the darkmatter density not given by Shandarin et al. (2011). We presented a novel approach to better understand thedynamics of cold collisionless fluids. We apply it by post-processing cosmological N–body simulations and documentthe significant improvement it represents over previous at-tempts to quantify the macroscopic and microscopic aspects ofthe dark matter fluid flow. In particular, we show new resultsfor density estimates, a dark matter “entropy”, bulk velocities,velocity distribution functions – many of which are computedhere for the first time. We are confident that our approachto tracing the dark matter sheet in phase-space gives impor-tant physical insights which were inaccessible with previousapproaches.
ACKNOWLEDGEMENTS
T.A. is grateful for numerous conversations with Greg Bryanin the past ten years on how one might solve for dark matterdynamics directly in phase-space and acknowledges support bythe National Science foundation through award number AST-0808398 and the LDRD program at the SLAC National accel-erator laboratory as well as the Terman fellowship at StanfordUniversity. He also acknowledges help by Patrick Abel in con-structing tetrahedra from paper which helped considerably inunderstanding the many possible options of tessellations of thedark matter sheet.
REFERENCES
Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93Afshordi N., Mohayaee R., Bertschinger E., 2009,Phys. Rev. D, 79, 083526Alard C., Colombi S., 2005, MNRAS, 359, 123Arag´on-Calvo M. A., Jones B. J. T., van de Weygaert R.,van der Hulst J. M., 2007, A&A, 474, 315Arag´on-Calvo M. A., Platen E., van de Weygaert R., SzalayA. S., 2010, ApJ, 723, 364Arnold V. I., Shandarin S. F., Zeldovich I. B., 1982, Geo-physical and Astrophysical Fluid Dynamics, 20, 111Ascasibar Y., 2010, Computer Physics Communications, 181,1438Bernardeau F., van de Weygaert R., 1996, MNRAS, 279, 693Bertschinger E., 1998, ARA&A, 36, 599Binney J., 2004, MNRAS, 350, 939Bryan G. L., Norman M. L., 1997, in Clarke D. A., WestM. J., eds, Computational Astrophysics; 12th KingstonMeeting on Theoretical Astrophysics Vol. 123 of Astronom-ical Society of the Pacific Conference Series, Simulating X-Ray Clusters with Adaptive Mesh Refinement. pp 363–+Buchert T., Bartelmann M., 1991, A&A, 251, 389Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., KlypinA. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240Colberg J. M., Pearce F., Foster e. a., 2008, MNRAS, 387,933Couchman H. M. P., 1991, ApJ, 368, L23Cuperman S., Harten A., Lecar M., 1971, Ap&SS, 13, 411Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Pot-ter D., Stadel J., 2008, Nature, 454, 735Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985,ApJS, 57, 241Frenk C. S., White S. D. M., Bode P., Bond J. R., BryanG. L., Cen R., Couchman H. M. P., Evrard A. E., GnedinN., Jenkins A., Khokhlov A. M., Klypin A., Navarro J. F.,Norman M. L., Ostriker J. P., Owen J. M., Pearce F. R.,Pen U., Steinmetz M., Thomas P. A. e. a., 1999, ApJ, 525,554Gibbs J., 1902, Elementary principles in statistical mechan-ics, developed with especial reference to the rational founda-tion of thermodynamics. Library of American civilization,Scribner’s sonsHahn O., Abel T., 2011, MNRAS, 415, 2101Hahn O., Porciani C., Carollo C. M., Dekel A., 2007, MN-RAS, 375, 489Hoeft M., M¨ucket J. P., Gottl¨ober S., 2004, ApJ, 602, 162Hogan C. J., 2001, Phys. Rev. D, 64, 063515Lynden-Bell D., 1967, MNRAS, 136, 101 c (cid:13) , 1–17 racing the dark matter sheet Maciejewski M., Colombi S., Alard C., Bouchet F., PichonC., 2009, MNRAS, 393, 703Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462,563Navarro J. F., Ludlow A., Springel V., Wang J., VogelsbergerM., White S. D. M., Jenkins A., Frenk C. S., Helmi A., 2010,MNRAS, 402, 21Neyrinck M. C., 2008, MNRAS, 386, 2101Peebles P. J. E., 1971, A&A, 11, 377Pelupessy F. I., Schaap W. E., van de Weygaert R., 2003,A&A, 403, 389Rasio F. A., Shapiro S. L., Teukolsky S. A., 1989, ApJ, 344,146Schaap W. E., van de Weygaert R., 2000, A&A, 363, L29Shandarin S., Habib S., Heitmann K., 2011, ArXiv e-prints,1111.2366.Shandarin S. F., Zeldovich Y. B., 1989, Reviews of ModernPhysics, 61, 185Sharma S., Steinmetz M., 2006, MNRAS, 373, 1293Shen J., Abel T., Mo H. J., Sheth R. K., 2006, ApJ, 645, 783Springel V., 2005, MNRAS, 364, 1105Springel V., White S. D. M., Frenk C. S., Navarro J. F.,Jenkins A., Vogelsberger M., Wang J., Ludlow A., HelmiA., 2008, Nature, 456, 73Springel V., White S. D. M., Jenkins A., Frenk C. S., YoshidaN., Gao L., Navarro J., Thacker R., Croton D., Helly J.,Peacock J. A., Cole S., Thomas P., Couchman H., EvrardA., Colberg J., Pearce F., 2005, Nature, 435, 629Springel V., Yoshida N., White S. D. M., 2001, New Astron-omy, 6, 79Stadel J. G., 2001, PhD thesis, UNIVERSITY OF WASH-INGTONStiff D., Widrow L. M., Frieman J., 2001, Phys. Rev. D, 64,083516Taylor J. E., Navarro J. F., 2001, ApJ, 563, 483Teyssier R., 2002, A&A, 385, 337Tremaine S., 1999, MNRAS, 307, 877van de Weygaert R., Schaap W., 2009, in Mart´ınez V. J.,Saar E., Mart´ınez-Gonz´alez E., Pons-Border´ıa M.-J., eds,Data Analysis in Cosmology Vol. 665 of Lecture Notes inPhysics, Berlin Springer Verlag, The Cosmic Web: Geomet-ric Analysis. pp 291–413Vera-Ciro C. A., Sales L. V., Helmi A., Frenk C. S., NavarroJ. F., Springel V., Vogelsberger M., White S. D. M., 2011,MNRAS, 416, 1377Vogelsberger M., Helmi A., Springel V., White S. D. M.,Wang J., Frenk C. S., Jenkins A., Ludlow A., Navarro J. F.,2009, MNRAS, 395, 797Vogelsberger M., White S. D. M., 2011, MNRAS, 413, 1419Vogelsberger M., White S. D. M., Helmi A., Springel V., 2008,MNRAS, 385, 236Wadsley J. W., Stadel J., Quinn T., 2004, New Astronomy,9, 137Waldvogel J., 1979, Zeitschrift Angewandte Mathematik undPhysik, 30, 388White S. D. M., Vogelsberger M., 2009, MNRAS, 392, 281Wise J. H., Abel T., 2007, ApJ, 665, 899Zel’dovich Y. B., 1970, A&A, 5, 84 c (cid:13)000