DraftversionFebruary2,2008 PreprinttypesetusingLATEXstyleemulateapj MASS OF CLUSTERS IN SIMULATIONS Andrea V. Maccio`1,2, Giuseppe Murante3 & Silvio A. Bonometto1,2 Draft version February 2, 2008 ABSTRACT We show that dark matter halos, in n–body simulations, have a boundary layer (BL), which neatly separatesdynamicallyboundmassfromunboundmaterials. WedefineT(r)andW(r)asthe differential kinetic andpotentialenergyofhalosandevaluatethemin sphericalshells. We noticethat,insimulated 3 halos, such differential quantities fulfill the following properties: (i) The differential virial ratio = R 0 2T/W has at least one persistent (resolution independent) minimum r¯, such that, close to r¯, (ii) the − 0 function w = dlogW/dlogr has a maximum, while (iii) the relation (r¯) w(r¯) holds. BL’s are set 2 where these th−ree properties are fulfilled, in halos found in simulationsRof T≃CDM and ΛCDM models, n runad–hoc, usingthe ART andGADGET codes; theirpresenceis confirmedinlargersimulationsofthe a samemodels,withalowerlevelofresolution. Herewefindthat 97%ofthe 300clusters(permodel) J we have with M >4.2 1014h−1M , owns a BL. Those clusters∼which appear∼not to have a BL are seen ⊙ · 4 to be undergoing major merging processes and to grosslyviolate sphericalsymmetry. The radius r¯ rc ≡ 1 has significant properties. First of all, the mass M it encloses almost coincides with the mass M , c dyn evaluated from the velocities of all particles within r , according to the virial theorem. Also, materials c 1 at r > r are shown not be in virial equilibrium. Using r we can then determine an individual density c c v contrast ∆ for each virialized halo, that we compare with the ”virial” density contrast ∆ 178Ω0.45 7 c v ≃ m (whereΩ isthematterdensityparameter)obtainedassumingasphericallysymmetricandunperturbed 4 m fluctuationgrowth. Asexpected,foreachmassscale,∆ iswithintherangeof∆ ’s. Howeverthespread 2 v c in ∆ is wide, while the average ∆ is 25% smaller than the corresponding ∆ . We argue that the 1 c c ∼ v matching of properties derived under the assumption of spherical symmetry must be a consequence of 0 3 an approximate sphericity, after violent relaxation destroyed features related to ellipsoidal non–linear 0 growth. On the contrary, the spread of final ∆c is an imprint of the different initial 3–D geometries of / fluctuations and of the variable environment during their collapse, as suggested by a comparison of our h results with Sheth & Tormen analysis. p - Subject headings: galaxies: clusters: general — cosmology: theory — dark matter — large-scale o structure of the Universe — methods: N-body simulations r t s a 1. INTRODUCTION preserve such symmetries during their growth and feel no : influence fromtheir environment(Gott & Rees 1975,Pee- v The inner region of dark matter halos in n–body sim- bles 1980). However, even under such simplifying condi- i ulations is certainly virialized. A precise detection of the X tions,acosmologicalconstantisasignificantcomplication. boundaries of such a region, however, is a long–lasting Lahav et al (1991) studied the growth of uniform spher- r and not yet completely solved problem. If we try just to a ical fluctuations within such context. They also outlined identify a “collapsed” region, i.e. a region detached from that the most direct way through which observations can the cosmological expansion, the task is easier (see, e.g., distinguish among flat models with different Ω (matter Monaco & Murante 1999 and references therein). How- m densityparameter)iscomparingclusterabundancesatlow ever,ingeneral,acollapsedregionisnotyetinvirialequi- and high redshift. The dynamics of the spherically sym- librium. Thepointisthat,inacosmologicalcontext,each metricgrowthofauniformdensityfluctuation,inamodel structureisembeddedinamoreorlesshomogeneousback- withΩ +Ω =1(whereΩ isthevacuumdensityparam- ground, and we pass gradually from the inner halo to ex- m Λ Λ eter),wasfurther exploredby Ekeetal(1996). Brianand ternalmaterials. Theaimofthisworkistoshowaprecise Norman (1998) gave a polynomial expression for the ex- criterion able to find virialized regions in dark matter ha- pected density contrast of virialized clusters, which reads los. This criterionwillbe testedonsimulationsandfound ∆ =A+Bx+Cx2,withx=Ω 1,andA=18π2 178, to be more efficient than previously adopted rules. v m− ≃ B =82,C = 39,andholdsatanyredshiftz;thisrelation As a matter of fact, several global properties of galaxy − is well approximated by the expression: clusters, such as mass (M), radius (r) and, also, temperature (TX) or luminosity (LX) are strictly linked to the ∆v =AΩm0.45 . (1.1) achievementofequilibrium. Analyticalpredictionsonsuch globalquantitieswereperformedassumingthatstructures Under the same assumptions, using a Press & Schechter evolve from uniform, spherically symmetric fluctuations, (1974) approximation, one can predict the expected mass function of clusters. 1DepartmentofPhysicsG.Occhialini,Universita`diMilano–Bicocca,PiazzadellaScienza3,I20126Milano(Italy) 2I.N.F.N.,ViaCeloria16,I20133 Milano(Italy) 3I.N.A.F.,OsservatorioAstronomicodiTorino–Torino(Italy) 1 2 Mass of clusters in simulations The impact of this approximation on the estimate of aboundary layer(BL),whosedepth∆rrangesaround50– clusterglobalquantitiescouldbeappreciatedonlythrough 100h−1kpc, and whose properties fulfill precise analytical numerical simulations. Much work was therefore devoted prescriptions. Such a BL has been found in 97% of the ∼ to this aim (see, e.g., Lacey & Cole, 1993,1994; Cole & clustersinoursimulatedmodels andis oneofthe findings Lacey 1996; Carlberg et al., 1996; Eke et al. 1996; Eke of this work. et al. 1998; Brian & Norman, 1998; Gardini et al. 1999). The BL is not a physical confinement barrier, and par- However, in the very analysis of the numerical work, the ticles can traveloutwardsand inwards throughit. On the above approximations were often implicit. An example is contrary, it confines a volume, where material fulfills pre- theidentificationofdarkhalosbasedonsphericaloverden- ciseequilibriumconditions,whichceasetoholdoutsideit. sity(SO)algorithms,withreferencetothedensitycontrast Inside the BL itself, the formal condition for virial equi- ∆ givenineq.(1.1). InthisworkweshallalsouseanSO librium of an isolated system holds, as pressure forces at v algorithm(describedbelow)asaselectioncriterion. How- its upper and lower borders cancel each other; on the BL, ever, using eq. (1.1), to work out values for the cluster kineticenergyhasaminimumwithrespecttopotentialen- radius (r ) or mass (M ), implies a bias. The values (1.1) ergy,andmatterdensity israpidlydecreasing. Thesecon- v v areobtainedunderpreciserestrictionsandnumericalwork ditions follow from a precise analytical requirement (Rw should also verify their impact. In a sense, the problem is requirement) derived under the assumption of spherical even more severe when the values (1.1) are used to define symmetry. This Rw requirement will be only approxi- virialized halos without requiring explicitly that they are matelysatisfiedbyhalosinsimulations,wherethehaloge- spherical. ometriesapparentlyviolatesphericalsymmetry. Weshow, In this workwe plan to find a rule, suitable for defining however,thatnormaldeviationsfroma sphericalshapein virialized systems, taking into account that a variety of DM halos do not prevent them from having a BL. growth histories led to different final features of virialized We explicitly note that we do not force the BL to exist, structures. A number of recent papers (Sheth & Tormen but we apply the Rw requirement and find it. For this (1999, 2002), Sheth, Mo & Tormen, 2001) had a similar reason, finding the BL is different from finding a given motivation. They tried to study the growth of primeval spherical overdensity; the latter always exists when aver- fluctuations takingintoaccountthe effectofanellipsoidal age densities of halos decrease with increasing radii, and collapse. Then, in order to avoid having the collapse on this does not imply even an approximately sphericalsym- some axis going to zero, they set up a recipe to freeze it metry. Onthe contrary,wecouldhavebound andrelaxed outonceithasshrunkbysomecriticalfactor. Thisfreeze– regionswhicharenotconfinedbyaBL:ouranalysisofnu- out radiusis explicitly chosenso that the density contrast mericalsimulations shows that this is not the case. There at virialization for the whole halo is again ∆ as given in is alsoa majordifference between the Rw requirementwe v eq. (1.1), i.e. in spherical growth models. In contrast, in are defining here, and other prescriptions, such as SKID, thisworktheuseof∆ todefinevirializedsystemsiscom- which defines a halo as the set of all particles gravitation- v pletely overcome. The point is that the detailed growth ally bound to the halo itself. In fact, bound particles are history of individual halos, in which different ellipticities not guaranteed to belong to virialized zones and may be- and velocity anisotropies played an important role, is not long to regions which are bound but still unrelaxed. expected to leave major traces on final halo shapes af- In this paper we shall first discuss when the BL should terviolentrelaxationhasoccurredduringthevirialization exist and introduce the Rw requirement used to detect process (Navarro, Frenk & White 1997). On the other it. To have significant statistics based on this rule and hand, it may result in a significant spread of the equilib- toconfirmitsvalidityrequiressimulationsperformedwith rium density contrasts ∆ of individual halos. Therefore, a large dynamic range, as can only be done by parallel c when the virial radius of relaxed objects is sought, devi- computing. ations from spherical symmetry and velocity anisotropies Finding and situating BL’s allows a better evaluation may not be as important as deviations of ∆ from ∆ . of individual cluster masses, radii and density contrasts, c v Weexpectthatthevalue∆ ,givenbyeq.(1.1),iswithin as well as an analysis of their possible dependence onvar- v the range of the actual ∆ values, but we think that it is ious physical parameters. The study of the structure of c important to first test how extended this range is, and dark matter halos in numerical simulations can also be how ∆ is set with respectto actual values. Furthermore, put on a sounder footing. In previous analyses (see, e.g., v a different ∆ implies a different radius and mass (r and Cole&Lacey1996)ther dependence oftheintegralvirial c c M ) for each object. Accordingly, we expect that some ratio 2T(< r)/W(< r) was discussed and found not to c − individualobjectmassesincreaseandothersdecrease. We be “useful for defining the boundaries of the virialized re- will explore the extent ofthese changeson galaxyclusters gion”. One of the findings of our analysis is precisely the innumericalsimulationsofcriticalCDMmodelswithand possibility of detecting such a boundary, using the virial without a cosmological constant. Such changes may have ratioforshellsinsteadofspheres. Fig.1comparesintegral an impact onthe mass function (see section 5 and Fig. 13 anddifferentialvirialratios,showinghowmuchmoreinfor- and14),whilesomescalingrelationsmightbesignificantly mationappearstobecontainedinthelatter(seeSec.3for altered, for instance the relation between M and T . detailsonhow (r)pointsarecalculated). Previousanal- X R In order to implement our program, we need a simple yseshavealsostressedhalopropertiesapparentlydepend- and effective rule to identify virialized material in clus- ing “on how groups are identified” (Sheth et al., 2001). ters. In trying to find such rule we discovered a precise The existence of a border, with precise physical proper- regularity in the transition region between material be- ties,secludingclustermateriallets oneovercomeanysuch longing to clusters and the surrounding medium. In fact, ambiguity. the volume occupied by virializedmaterials is confinedby It may be also important to remember that according Macciò et al. 3 to the definitions given above, in this work ∆ , r , M length, is η 112h−1kpc, yielding a Plummer equivalent v v v are quantities worked out starting from the fixed density softening ǫ ≃ 40.6h−1kpc. These simulations were de- pl ≃ contrast given by (1.1), while ∆ , r , M are the same scribed in more detail in Gardini et al. (1999), and were c c c quantities when worked out from the setting of the BL also used in Macciò et al. (2002). The parameters of the around each cluster. Let us also point out that CHb the models are reported in detail in Table 1. procedure discussed in this work is not the basis for an In order to see how resolution affects the results of CHe alternative cluster–finding algorithm, such as, e.g., the Rw requirement, we performed further simulations Friends-of-Friends,SKIDorSO.Onthecontrary,ouranal- with the ART code (courtesy of A. Klypin). This al- ysisassumesthatclusterlocationsaregiven. Herewegive lows us to select regions inhabited by clusters and to re– adifferentprescriptionfordefiningthephysicaldimension run them with increasing particle mass resolution: from of the cluster, not its position. 7.7 1011Ω h−1M to 1.2 1010Ω h−1M . m ⊙ m ⊙ × × The planof the paper is as follows. In section2 we give TheARTcode(AdaptiveRefinementTree: Kravtsovet suitable information on the simulations used, which were al. 1997, Knebe et al. 2000) starts with a uniform grid, partially run ad–hoc, for this work. In section 3 and 4 which covers the whole computational box. This grid de- we define the boundary layer (BL) and discuss how it can fines the lowest (zeroth) level of resolution of the simu- be found in simulated clusters. In section 5 we show the lation. The standard Particles-Mesh algorithms are used results of our work, and in section 6 we discuss them and to compute the density and gravitationalpotential on the some future perspectives. zeroth-level mesh. The code then reaches high force resolution by refining all high density regions using an auto- mated refinement algorithm. The refinements are recur- 2. MODELSANDN–BODYSIMULATIONS sive: the refined regions can also be refined, each subse- Weusetwomainsetsofsimulations,performedwithdif- quent refinement having half of the previous level’s cell ferent codes: the parallel AP3M N–body code described size. This creates a hierarchy of refinement meshes of dif- in Gardini et al. (1999), the parallel PM ART code de- ferent resolution, size, and geometry covering regions of veloped by Kratshov & Klypin (1997), and the parallel interest. Becauseeachindividualcubiccellcanberefined, tree–code GADGET (Springel et al. 2001). the shape of the refinement mesh can be arbitrary and The code of Gardini et al. (1999) was developed from effectively match the geometry of the region of interest. theserialpublicAP3McodeofCouchman(1991),extend- This algorithm is well suited for simulations of a selected ing it to different cosmological models and allowing suit- region within a large computational box. ableflexibilityinparticlemasses. Usingthiscode,twosim- Thecriterionforrefinementisthelocaldensityofparti- ulations were performed. The first of them, dealing with cles: ifthenumberofparticlesinameshcell(asestimated a “tilted” Einstein-de Sitter model (hereafter TCDM), is bytheCloud-In-Cellmethod)exceedstheleveln ,the thresh the same already considered by Gardini et al. (1999); the cell is split (“refined”) into 8 cells of the next refinement normalization of the run was, however, rescaled to yield level. Therefinementthresholdmaydependontherefine- σ = 0.55 at the final step, instead of σ = 0.61, so that ment level. The code uses the expansion parameter a as 8 8 the final abundance of rich clusters is the same for both the time variable. Besides spatial refinements, during the modelsatthefinalepoch(asususal,σ isthem.s. density integration time refinements are also performed. 8 fluctuation on the scale of 8h−1Mpc, h being the Hubble The ART code can handle particles of different masses. parameter in units of 100 km/s/Mpc). The latter simu- This lets us increase the mass (and correspondingly the lation dealt with a ΛCDM model, i.e. a flat CDM uni- force) resolution of regions centered around the highest versewithnon–zerocosmologicalconstantwithσ =1.08. mass clusters. The multiple mass resolution is imple- 8 Owing to recent observationalresults (see, e.g., Schuecker mented in the following way. We first set up a realiza- et al, 2002), these normalizations are both rather high, tion of the initial spectrum of perturbations, so that ini- as a reasonable value for σ in ΛCDM models is 0.75; tial conditions for our largest number (5123) of particles 8 ∼ however, this has no impact when studying inner cluster can be generated in the simulation box. Coordinates and dynamical properties. velocities of all the particles are then calculated using all The above simulations (denoted A and B, respec- waves ranging from the fundamental mode k = 2π/L to tively) are performed in cubes with sides of 360h−1Mpc. the Nyquist frequency k = 2π/L N1/3/2, where L is CDM+baryons are represented by 1803 particles, whose the box size and N is the number o×f particles in the sim- individualmassis2.22 1012h−1M⊙ forTCDMand0.777 ulation. We obtained lower resolution regions by merg- 1012h−1M⊙ for ΛCDM·. We use a 2563 grid to compute· ing high resolution particles into particles of larger mass the FFT’s needed to evaluate the long range contribution where needed. This process can be repeated to get still totheforce(PM)andweallowformeshrefinementwhere lower resolution. The larger mass (merged) particles have theparticledensityattainsorexceeds 30timesthemean velocities and displacements equal to the average of the ∼ value. The starting redshifts are zin =10 for TCDM and velocitiesanddisplacements ofthe smaller-massparticles. zin =20 for ΛCDM. The particle sampling of the density Usingthistechniquewefirstlygeneratedinitialconditions fieldisobtainedbyapplyingtheZel’dovichapproximation forlow–resolutionsimulations. Thesesimulationswererun (Zel’dovich 1970; Doroshkevich et al. 1980) starting from from their initial redshift to z = 0 and they were used to a regulargrid. We adoptthe samerandomphasesinboth locate the Lagrangian regions forming the most massive A and B. clusters. Then, we were able to obtain initial conditions The number of steps were 1000 equal p–time steps (the athighresolution(both inforce andmass)in the selected time parameter is p a2/3, where a is the expansion fac- Lagrangian zones, with the same random Fourier phases. ∝ tor). Thecomovingforceresolution,givenbythesoftening 4 Mass of clusters in simulations Such zones were surrounded by lower and lower resolu- inmore detailandcompareit with other groupidentifica- tion shells taking into account the tidal field which acts tion algorithms, e.g. Governato et al. 1999). In this work on them. Finally, these ICs were used to run the high– we set φ = 0.2 and take N so to have a mass threshold f resolution cluster simulations. 3.0 1013h−1M . Above this mass threshold there are ⊙ × The simulations employed here were performed us- 10000halos in the simulations A,B and 1250halos in ing 1283 zeroth-level grid in a computational box of ∼the simulations C,D. ∼ 180h−1Mpc. The thresholdforcellrefinement(see above) Tests performed by Gardini et al. (1999) show that no was low on the zeroth level: n (0) = 2. Thus, ev- cluster is missed above 4.2 1014h−1M . We have 303 thresh ⊙ × ery zeroth-level cell containing two or more particles was suchclustersinA,316inB,40inC,36inD;6clustersof refined. The threshold was higher on deeper levels of re- C andD wherethenblownupto variousresolutionlevels. finement: n =3andn =4forthefirstleveland thresh thresh higher levels, respectively. 3. THEVIRIALTHEOREMAPPLIEDTOSINGLE For the low resolution runs the step in the expansion CLUSTERS:THEORY parameterwaschosento be ∆a =2 10−3 onthe zeroth Finding the volume wherecluster materialsareinvirial 0 × levelofresolution. Thisgivesabout500stepsforparticles equilibrium may seem a tough and somewhat ambiguous located in the zeroth level for an entire run to z = 0 and task. If kinetic and potential energies, within a radius r, 128.000 for particles at the highest level of resolution. are defined according to: Using the ART code we performed two simulations: a ΛCDMmodelandaTCDMmodel,thatwewillindicateas 2T( r . This simple condition, + − + − however, is not strictly necessary; it is sufficent that the (i) it is in virial equilibrium, difference w increases, becoming positive, at r > r . + R− This will be shown here below. Let us however antici- (ii) the r dependence of pressure is such that the r.h.s. pate that, while the Rw requirement can be fairly easily of eq. (3.4), yielding a term of the form 3(p V + + − tested, it is not easy to deal with such further require- p V ), can be neglected in the virial balance, − − ment. The tests performed on a large number of clusters (iii) inside ∆r, w is constant. anddescribed inthe next sections, however,let us conjec- ture that, when the Rw requirement is fulfilled, this fur- Further commentsonthe (ii) requirementcanbe foundin ther requirement is somehow inescapable. We now show AppendixA,whereweshowthathaloprofilesareexpected why, at r > r , the equilibrium conditions are violated if to approach it. Together with (i), it implies that + a boundary layer is found. Notice that, at r >r , we can + still assume that W = W(r )(r/r )−w, but w, starting d + + mv2 r Z(r )=0 , (3.9) from the value it had in ∆r, will now vary with r and, i − i dr i Xi (ri∈∆r) according to eq. (3.8), in the layer. Let us then take into account the condition that w is constantanduse eq.(3.8)with w′ =0 andthe virial ratio w =w′rln(r/r ) ; (3.13) + R− mv2 ∆r = Pi(ri∈∆r) i , (3.10) R (r ) Pi(ri∈∆r)Z i this equation cannot be fulfilled, if the difference w, R− to obtain that vanishing up to r , increases at r > r , becoming posi- + + =w , (3.11) tive. Infactwerequiredthatw ismaximumandtherefore R w′ <0. all along the interval ∆r. This equation coincides with Let us finally notice that if a layer is characterized by eq. (3.29) in Goldstein, Pole & Safko (2002) (in turn, higher w values, this means that it is emptier than inner the latter equation yields eq. (3.9) hereabove, in the limit or outer layers. A maximum of w, therefore, indicates a w const.). R≡equiring that w′ =0, then, also yields that low density layer. In turn, a minimum of R indicates a low kinetic energy, with respect to potential energy. BLs, d therefore, are partially depleted r–intervals, where parti- R =0 , (3.12) cles, on average,are particularly slow. dr 6 Mass of clusters in simulations 4. THEVIRIALTHEOREMAPPLIEDTOSINGLE expected to create severe problems. On the contrary, we CLUSTERS:PRACTICALUSE found that all clusters with mass M > 4.2 1014h−1M ⊙ ∼ · The mathematical definition must now be applied to have at least one value of rm fulfilling the above require- numerical simulations. A failure to find BLs could be at- ments, while some clusters have 2 or 3 points where the tributed to insufficient numericalresolutionorto badvio- fitting criteria are fulfilled. In the few cases, when the lations of sphericalsymmetry. However,in principle, even fitting criteria were satisfiedat more than one rm, we dis- if the resolution is enough and no substantial departures carded rm values yielding ∆c outside an interval 0.3∆v– from sphericity are visible, BLs could still be absent. 3∆v. This selection criterion left us with 97% of the ∼ Let us now describe our procedure. After finding clus- initialclusters. Whenmorethanonerm wasstillavailable, ters in simulations and their centers C , as described we selected the one corresponding to the smallest χ2; this W in sec. 2, we evaluated W(r) and (r), for spherical r– last criterion was needed for 10% of the clusters only. intervals containing a fixed numberR Let us recall again that the∼requirement of intersection betweenaminimum of andamaximumofw wasfound N =M /(24m ), (4.1) forsphericalsystems. WRhen we dealwiththe halosofour 4 p simulations, sphericity is likely to be an approximation. of particles (here M4 of the total mass within 4.5h−1Mpc When a single intersection exists, however, it is natural from CW and mp is the particle mass). However, succes- to guess that this is the point where the Rw requirement sive points along r are obtained by shifting outwards by is best approached. However, both in this case and when N/8 particles. In the relevantr range,the abovecriterion several r points are found, it is natural to perform a m yields N 102 for TCDM (simulation A) and N 3 102 closerinspection to make sure that no unexpected feature ∼ ∼ · for ΛCDM(simulationB),whilesuccessivepointslayata may bias our understanding of the physical situation. distance 20h−1kpc. Thisisbelow(butnotmuchbelow) This closer inspection requires an increased resolution. the force∼resolution ( 100h−1kpc), as desired. Asdescribedinsec.2,aTCDMmodelandaΛCDMmodel ∼ WedonotexpectBL’stobemuchthickerthanourforce were then reconsidered within a smaller box with sides of resolution;accordingly,wefirstseektheminimaof . Ow- 180h−1Mpc,makinguseoftheARTandGADGETcodes. R ing to the procedure by which points are worked out, This gave us 1/8 of the clusters provided by the AP3M a minimum is considered significRant only if it is such with simulations,b∼utenabledustofocusonparticularclusters, respectto8neighbors. Findingthemaximaofwisharder, making use of the ART package facilities to increase the as w(r), obtained by differentiating W(r), should then be resolutiontherewhilestillkeepingidenticalboundarycon- furtherdifferentiatedtofinditsmaxima. Inpractice,how- ditions. Inordertotestthattheregularitiesfoundarenot ever,thissecondorderdifferentiationisnotneeded. Quite linked to peculiarities of the ART package, some clusters in general,there will be just a few points rm where has were magnified using GADGET, instead of ART. Alto- R significantminima. Letthen m bethevirialratioatsuch gether we have 12 clusters magnified, 6 for each model. 4 R points and let us fit W(r) with the expression of the latter were obtained using ART and 3 using GAD- GET; hence, one cluster was magnified with both codes W˜ =A(r/rm)−Rm (4.2) showing that our results are independent of the n–body code used. In fact, the final positions of particles have in ν = 8 points > r , covering an interval ∆r 100– 150h−1kpc. The fit mis performed in two steps: w∼e first just minor differences, while R and W coincide. A number of subsequent tests also concerned the stability of the minimize the function φ = W/W˜ 1 just using A as a minima of . In general, the number of local minima r parameter; then, in eq. (4.2), we a−llow m to take val- decreases wRhen the resolution is increased, but the poinmt R ues different from the value of in rm, and perform a whereW fitstheexpression(4.1)withtheworstresolution R 2–parameter fit. is a minimum also with greater resolutions. In just a few Weconsideredthefittobesatisfactory,iftwoconditions cases, new minima arise when the resolution is improved; hold: (i) When performing the first fit, the residual but in no case is a minimum yielding the BL erasedwhen changingresolution. The“worst”casesforeachmodelare ν χ2 = φ2(r )<10−3ν , (4.2) shown in Fig. 2. Here we have a TCDM cluster, taken i from the simulationC, which passes from 3 minima (1283 X i=1 particles), to 1 minimum (2563 particles) and, again, to 3 (i.e., a mean square discrepancy 3%, between the nu- minima(5123 particles). TheBLisdisplacedoutwardsby merical values of W(r) and its fitt∼ing expression, was our 150h−1kpc, when we pass from 128 to 256; then, when ∼ limit of tolerance). This detects r–intervals where eqs. passing from 256 to 512, no appreciable shift occurs. Out (3.11)-(3.12) are both fulfilled. (ii) When performing the ofthe12clustersconsidered,onlythisonepresentsashift secondfit, the best–fitvalue of is within1–σ fromthe slightlyabovethe resolutionofthe initialsimulation. The m actual value. Here we test Rthat also the slope of W is ΛCDM cluster of simulation D shown in Fig. 2, shows 4 m nottooRfarfromwhatisrequired. Suchconstraintsshould minima with 1283 particles, 2 minima with 2563 particles therefore select points where is minimum and (nearly) and, again, to 3 minima with 5123 particles; but when intersects w. R changing resolution, the shifts are <40h−1kpc. As outlined above, one might reasonably expect that a For these 12 clusters we follow in detail the joint be- largefractionofclusters,attheirboundaries,isfarenough havior of and w. We show them in Fig. 3 for the best R from sphericity that the expected system of maxima and andworstcaseswefound. Ofcourse,aprecisecoincidence minima,aswellastheircoincidences,isdisrupted. Thenu- betweenaminimumof andamaximumofwwouldindi- R mericalnoiseincalculatingwisalsoquitelargeandcanbe catethat,withintheresolutionlimits,thesystemisspher- Macciò et al. 7 ical. Seekingthiskindofcoincidence,thereforeissomehow sity contrasts we found; the average ∆ , however,in both c fatuous. Furthermore, in order to work out w, we ought cases, is smaller than ∆ by 25%. v ∼ todifferentiateW,whichisobtainedfromadiscretepoint It is, however, clear that gravitationally evolved, sta- distribution. Hence, although the nice coincidence of the tionary objects are not uniquely identified with a given, best case is quite appealing,this case is not safer than the fixed density contrast. In further work, we plan to deal worst case we show. By tolerating a 3% disagreement in with the impact of variable ∆ on the physical character- c W, for r/r < 1.1, we allow for a tolerance 0.3 on istics ofDM halos,suchas concentration,averageangular m ∼ ∼ theexponentw. Intheworstcase,wefindaw– discrep- momentum, density and velocity profiles. It may also be R ancy of this order, but in most cases, w– turns out to significanttotestourresultsagainsttheeffectsofdifferent | R| be 0.1. We tentatively conclude that this resolution is particle selection criteria, e.g. using the SKID algorithm ∼ sufficient to show physical features and that the level of instead of the SO one to pre-select the DM halos. the w– (dis)agreement is a measure of the a–sphericity In Figs. 8 and 9 we plot r against M , for the simu- c c R of actual systems. The important issue is that we always lations A and B, respectively. Here, error bars account find a maximum of w to match the minimum of , quite for the variance of r around its average, which is wider c R close to the r selected at the initial resolution level. than the error of the averager , shown in Fig. 6. Contin- m c This point can be further illustrated by Fig. 4, which uous and dotted curves show the expected behavior of r c shows a plot similar to Fig. 3, for a typical case at the against M , when the density contrast is set either at the c initial resolution level. Both the and the w curves are average∆ value or at ∆ . c v R muchnoisierhere. ,inparticular,showsawidesetoflo- Itisalsosignificanttocomparetheplotofr againstM c c R calminima;theonlypersistentminima,however,arethose with the behaviour of the radii r against the masses M c v marked by thick filled circles. As expected, the behavior of the corresponding clusters. This comparison is shown of w is still noisier and most of the oscillations shown are in Fig. 10, for the simulation A only; in order to avoid clearlynumerical. Theroleofw,however,justamountsto confusion,we omiterrorbarswhich,however,areapprox- choosing among the points marked by a thick filled circle imately of the same size in both cases. The highest mass and selects the point marked by a vertical arrow. Even points, for which error bars cannot be evaluated, are also at this resolution level, the presence of a maximum of w, omitted. Filled circlesyield r vs. M , empty circlesyield c c closetosuchpoint,canbesuspected. Whentheresolution r vs. M . Filled circles nicely fit the line worked out c v is increased, this maximum is more clearly visible, while fromtheaveragedensitycontrast∆ ;onlyathighmasses, c no w maxima exist close to the other minima. where the statistics is poorer, some of the circles are far R A further preliminary test of our technique, concerning from this line. The same line is also a reasonable fit for galactic-sized and unvirialized DM halos, is given in Ap- the open circles. These, however, are systematically far- pendix B. ther at all mass scales. In the figure, however, we also Before concluding this section, let us draw the reader’s plot an horizontal line showing the average of the r val- c attentiononthesettingoftheBLaroundparticleagglom- uesofallclusterswithmass>4.2 1014h−1M . Eventhis ⊙ · eration, as shown in Fig. 5 for a ΛCDM cluster. This fig- averagecan be considered a reasonablefit for open circles ure shows that the BL is actually set outside of the main (anunweightedχ2 evaluationgivessimilarprobabilitiesto matter agglomeration; its spherical shape is somehow in both fitting lines). This plot shows that there are actual contrastwith the apparenta–sphericityofthe cluster ma- dangers, if quantities workedout assuming a constant ∆ v terialsandeasilyjustifiesminordiscrepanciesfromresults are compared with real data. In this case, one might ten- based on an assumption of full spherical symmetry. tatively suggestthat there is a typical, mass independent, cluster radius R 1.56h−1Mpc, quite close to the Abell 5. RESULTS value. Ourproced≃ure,instead,outlinesthatthisarisesbe- Once the sphere confining cluster materials is set, we causesomeclusterswereattributedabiasedmassandthe can evaluate the density contrast∆ and the mass M in- corresponding points were then set at a wrong abscissa. c c side it. In Figs. 6 and 7, points give ∆ and M for all Although, altogether, this abscissa reshuffling is substan- c c clusters in simulations A and B, respectively. They show tially casual, the danger that an characteristic scale erro- a rather wide spread of ∆ values. By subdividing the neously appears is far from absent. c M abscissain intervalsof constantlogarithmicwidth, we A critical result of our analysis is however shown in c evaluate the average density contrast in each of them, to Figs. 11 and 12 (for simulations A and B, respectively). seek systematic trends with mass. Averages are still sub- Insimulations,the “numerical”clustermassescanbeeas- jecttoasignificantuncertaintyowingtothespreadof∆ . ily evaluated by summing up particle masses. The actual c They are shown, at the 1–σ level, in the plots. At high approach, in the real world, ought to be quite different masses, some logarithmic intervals of Fig. 7 are empty or andisbased,firstofall,onananalysisofgalaxyvelocities contain a single object. In the latter cases no error bar is or X-ray emitting gas temperature. These analysis make given. use of the virial hypothesis; therefore, measured masses There seems to be no evidence of any peculiar trend of ofgalaxyclustersareestimatesoftheir dynamical masses. density contrasts with mass apart from, perhaps, a mod- It is therefore important to test how “numerical” cluster est indication of an increasing density contrast at very masses in simulations agree with the masses high scales in ΛCDM. It is therefore possible to consider v2 r the overallaverageamong ∆c’s. This averageis indicated Mdyn = h i{v,c}· {v,c} , (5.1) by the continuous horizontal line and compared with the G ”virial”densitycontrast∆ ,asgivenbyeq.(1.1). Inboth evaluated averaging over the velocities of the particles v cases, ∆ (dotted line) is well inside the range of the den- within r . This comparison has been often performed, v v,c 8 Mass of clusters in simulations showing a reasonable coincidence between M and M . arises when M is replaced by M . (ii) Variations how- v dyn v c The results of such fit are shown in Figs. 11 and 12 by ever exist, which are of the same order of the differences the dashed histograms, which confirm the slight excess of between the standard theoretical expression(5.3; Press & M vs. M already noticed by previous authors. In the Schechter) and the improved theoretical expression (5.4; dyn v same plots we also show the results of a fit between M Sheth&Tormen). (iii)Whenthelatterexpressionisused, dyn and M , obtained on the basis of the setting of the BL. the fits to the numerical mass functions are best if differ- c There seems to be little doubt that the latter fit is better. ent values of the parameter a are used for M and M v c The average values of M /M are both 0.97 0.03, ( 0.7 and 0.5, respectively). The differences between the dyn c ∼ ± ∼ quiteclosetounityandconsistentwithitatthe1–σlevel. mass functions n (>M ) and n (>M ) are more evident c v c c We have also compared the shapes of the cluster (inte- in the intermediate mass range (7 1014h−1M < M < ⊙ gral) mass functions n (> M ) and n (> M ). In Fig. 13 1015h−1M ). They arisefromawid·espread“re-shuffling” c c c v ⊙ and14,upperandlowerplotsshowtheirbehaviors,respec- among the masses of DM halos, when the two definitions tively. In the same plots, we also show Press & Schechter areused,causingtherebydifferentnumbersofhalosindif- and Sheth & Tormen theoretical curves. They can be ob- ferent mass intervals. This can also be seen from Fig. 10, tained from two different expressions of where the mean cluster radii (r , r ) of halos in this mass c v range show appreciable differences. M f(ν)νdlogν = nc(M)MdlogM ; (5.2) 6. CONCLUSIONS ρ¯ m Observationalworkongalaxyclustersmadeagiantleap here ρ¯m is matter density, the bias factor ν = δc/σM is forwardthankstoX–rayobservationsand,inturn,cluster theratioofthecriticaloverdensityinthesphericalgrowth data are providing more and more cosmological informa- model to the rms density fluctuation on the length scale tion. The ROSAT AllSky Survey(RASS) allowedthe ex- corresponding to the mass M, and nc(M) is the differen- traction of large samples of galaxy clusters, like REFLEX tialmassfunction,whichshallthenbeintegratedtoobtain (Boehringer et al. 2001) or NORAS (Viana et al. 2001) , nc(>M). tracing the large scale structure (LSS) of the Universe up In the Press & Schechter formulation, assuming spheri- to 1000h−1Mpcofdistance(histheHubbleconstantin cal unperturbed fluctuation growth, uni∼ts of 100 km/s/Mpc). These data are consistent with the spectrum of a low–density CDM model, with matter f(ν)ν = 2/π νexp( ν2/2) , (5.3) density parameter Ω 0.3. Also an analysis of the red- m − ≃ p shiftdependence ofthe X–rayclusterluminosity function, while Sheth & Tormen,taking into accountthe ellipsoidal based on a Press & Schechter approximation, sets a rela- nature of the collapse, argue that tion between Ω and σ , which requires again Ω < 0.5, m 8 m 0.6 (Ω : matter density parameter; σ : squared mean m 8 f(ν)ν =A(1+ν′−2q) 2/π ν′exp( ν′2/2) , (5.4) density contrast, on the scale of 8h−1Mpc) (Borgani et − p al. 2001). Such constraints on the model will become with ν′ =√aν; this expressionintroduces the parameters morestringent,intheforthcomingyears,thankstoXMM– A, q and a, however constrained by the requirement that Newton and Chandra data, providing high–z cluster sam- all matter is bound in objects of some mass M, however ples with tenths of thousands objects (see, e.g., Bartlett small, andthereforethe integraloff(ν)mustbe unity. In et al. 2001). our plots, A is normalized to yield the number of clusters In our opinion, such fast improvements of the observa- found in simulations, and this constraint is automatically tionalmaterialsdeservetobe matchedbyasuitableeffort satisfied. Sheth&Tormen(1999)showedthat,intheGIF of the theoretical analysis. In this work we try to con- simulations for anSCDM, 0CDMandΛCDM models (see tributetoit,byprovidinganimprovedpatterntomeasure Kauffmann et al 1999), a good fit to the numerical mass individual cluster densities, radiiand masses in numerical function was obtained if q =0.3, a=0.707. simulations, without reference to a standard density con- In Fig. 13 and 14 we report the mass function obtained trast. Although, both in real samples as well as in simu- fromeq.(5.4)for these values,forTCDM andΛCDM,relations, cluster members are directly observable together spectively. Thefitwiththen (>M )issurelybetterthan with their environs, a definition of cluster radii, masses, c v using eq. (5.3). We also tried to modify the parameter a, and density contrasts is needed both for statistical aims, that, according to Sheth & Tormen arguments, is related when treating a large set of halos, and for individual ob- to the density contrast of virialized structures. As a mat- jects,e.g.,tocomparetheobservedmassdistributionwith ter of fact n (> M ) is derived with a variable density analytical expressions of radial profiles. In the standard c c contrast, as explained in previous sections, and, however, approach, a value for the density contrast is usually bor- with an averagefinal density contrastwhichis 25%be- rowedfromanalyticalresults,whichareobtainedunderas- ∼ low ∆ . It is therefore reasonable to expect that a lower sumptionswhichcanbeincorrectintherealworld,aswell v value of a is needed to compare theoretical predictions asinsimulations. Suchassumptions,however,wereexplic- with the numerical behavior of n (> M ). In Fig. 13 and itlyconsideredtobethebestpossiblechoice. Inparticular, c c 14, we show the expression (5.4), both for a = 0.707 and Cole & Lacey (1996)had explored the possibility to make 0.5. It is also evident that, while a = 0.707 provides a use of the virial ratio, but considered its value in spheres, good fit for M –mass functions, the fit obtained provided instead of layers. Here we see that, passing from the inte- v by a=0.5 for the M –mass function is perhaps better. gral to the differential virial ratio (r) = 2T(r)/W(r), c R − The results of these comparisons, altogether, are that: allowsthedetectionofaneatboundary,whereatransition (i) No major variations in the numerical mass function from halo to surrounding material occurs. In fact, first of Macciò et al. 9 all, (i) this inspection allows the detection of minimum where we had a large statistics of clusters (> 300 with pointsof inallhalos,whichcouldbeeasilymissedinan M > 4.2 1014h−1M , for each model). We applied to c ⊙ R · integral inspection. Some of them are safely shown not to all of them an Rw–requirement, meant to detect the BL, be numericalfeatures, as they persistwithout appreciable and found it in 97% of the clusters considered. A de- shifts when we comparedifferentresolutionsimulationsof tailed inspection of these clusters shows that the mass M c a given halo. Furthermore, (ii) we found that, CHb quite within the BL substantially coincides with the mass M , dyn close to the radius r¯, where one of the persistent minima evaluated from the velocities of all particles within r =r¯, c occurs, there is always a maximum of the function w(r); according to the virial theorem. Furthermore, a direct in- finally, CHb (iii) at the radius r¯, we have also that spection of the clusters without a BL, showed they were exceptional structures, mostly due to undergoing major w(r¯) (r¯) . (6.1) ≃R merging,whichcertainlyjustifieda dynamicallyunsettled Let us recall that, in general, w(r) is defined with refer- condition and a badly a–spherical geometry. Let us how- ence to the behavior of the potential energy around r¯, by ever notice (see Fig. 5) that “standard” a–sphericity of setting: halostructuresdoesnotappeartobeaseriousobstacleto W(r)=W(r¯)(r/r¯)−w(r) (6.2) BL detection. At present, model clusters to fit X–ray data are mostly and, therefore, its r dependence accounts for the radial built using hydro–codes. This allows substantial improve- dependence of W(r). ments with respect to initial analysesbasedon pure New- In this paper we have shown that the above properties tonian dynamics. There are, however, some features in (i),(ii),(iii)areconsistentwiththerequirementsthat: (a) the interpretationofboth kinds ofsimulations,whichrely r¯sets a neat boundary between bound and unbound material. (b)Asignificantinterval(depth 50–100h−1kpc), on a sound definition of the virialization radius Rc. This ∼ radius, whose setting immediately follows from the detec- wherew(r)isconstant,liesaboutr¯. Wecalleditabound- tion of the BL position, corresponds to a wide spread of ary layer (BL). Accordingly, the r dependence of W(r), density contrasts ∆ , and, in general, is significantly dif- inside of the BL, is expressed by eq. (6.1), by replacing c ferent from the radius R , obtained using a given density the exponent w(r) with its constant value w¯. v contrast ∆ . Even keeping to models without hydrody- The relationbetween the observedfeatures(i), (ii), (iii) v namics,itisclearthatusingdifferentradiiR anddensity and the properties (a), (b), was shownunder the assump- c contrast ∆ could lead to different relations between L , tion of spherical symmetry. The fact that the conditions c X T and M . As a consequence the one–to–one correspon- (i), (ii), (iii) are not satisfied exactly, but with a (slight) X c dence between M and T could be lost and it would no approximation, can be easily attributed to deviations of c X longerbe grantedthatamassfunctionandatemperature numerical structures from sphericity. However, such devi- function are substantially equivalent quantities. ations ought to be quite mild, owing to the statistics we performed on a large number of halos, and this indicates that violent relaxation erased most previous a–spherical ACKNOWLEDGEMENTS features. Non–spherical collapses, however, leave a sub- stantialimprintonthefinaldensitycontrasts∆ foundfor We thank INAF for allowing us the CPU time to per- c virialized halos whose spread accounts for different initial form the ART simulation C and D at the CINECA con- 3–D geometries of fluctuations and different interactions sortium (grant cnami44a on the SGI Origin 3800 ma- with the environment during the collapsing stages. chine). GADGET simulations of high–resolution clusters The simulations used to focus the above properties, at have been run on the 16 Linux PC Beowulf cluster at the differentlevelsofresolution,inTCDMandΛCDMmodes, Osservatorio Astronomico di Torino. Alessandro Gardini were run both with the GADGET and the ART codes. is thanked for allowing us to make use of the simulations Results do not depend on the code. We could also ex- A and B; we also wish to thank him, Loris Colombo and tend the analysis, at the lowest resolution level, to sim- Giuseppe Tormen for useful discussions during the prepa- ulations run in quite large boxes (360h−1 Mpc aside), ration of this paper. REFERENCES Bartlett J., et al., 2001 “Galaxy Clusters and the High Redshift Gott,R.&Rees,M.1975,A&A,45,365G Universe Observed in X-rays” Proceedings of the XXI Moriond Governato F.,Babul A.,QuinnT.,Tozzi P.,BaughC.M.,KatzN., Conference,(preprintastro-ph/0106098) &LakeG.,1999, MNRAS,307,949 Boehringer,H.,etal.,2000, ApJS,129,435B Kauffmann G., Colberg J., Diaferio A. & White S.D.M. 1999, Borgani,S.,etal.,2001,ApJ,561,13B MNRAS,303,188 Brian,G.&Norman,M.,1998, ApJ,495,80 Klypin A., Kravtsov A., Bullock J. & Primack J., 2001. ApJ 554, Carlberg,R.G.,etal.,1996,ApJ,462,32 903 Cole,S.&Lacey, C.,1996,MNRAS,281,716 Knebe A., Kravtsov A., Gottlo¨ber S. & Klypin A., 2000 MNRAS, Couchman, H.M.P.,1991, ApJ,268,L23 317,630 Doroshkevivh, A.G., Kotok, E., Novikov, I.,Polyudov, A.N., KravtsovA.,KlypinA.&KhokhlovA.,1997ApJ,111,73K Shandarin,S.F.&Schctman, S.,1980,MNRAS,192,321 Lacey,C.&Cole,S.,1993,MNRAS,262,627 Eke,V.R.,Cole,S.,&Frenk,C.S.,1996,MNRAS,282,263 Lacey,C.&Cole,S.,1994,MNRAS,271,676 Eke,V.,Navarro,J.&Frenk,C.,1998,ApJ,503,569E Lahav, O., Lilje, P.R., Primack, J.R. & Rees, M., 1991, MNRAS, Gardini,A.,Murante,G.&Bonometto, S.A.,1999,ApJ,524,510 282,263E Ghigna S.,MooreB.,Governato F.,LakeG.,QuinnT.&Stadel J. MaccioÀ.V.,GardiniA.,GhignaS.&BonomettoS.A.,2002,ApJ, ,2000,ApJ544,616 564,1M GirardiM.,BorganiS.,GiuricinG.,MardirossianF.&MezzettiM., Monaco,P&Murante,G.,PhysicalReviewD,60/10,id.103502 , 1998,ApJ,506,45G 1999 Goldstein H., Poole C. & Safko, 2002, Classical Mechanics, 3rd ed. NavarroJ.,FrenkC.&WhiteS.D.M.1996, ApJ,462,563 –AddisonWesleypub. NavarroJ.,FrenkC.&WhiteS.D.M.1997, ApJ,490,493 10 Mass of clusters in simulations Peebles P.J.E., 1980, The Large Scale Structure of the Universe, ShethR.K.&TormenG.,1999MNRAS,308,119 PrincetonUniversityPress,Princeton ShethR.K.&TormenG.,2002MNRAS329,61 Power C., Navarro J., Jenkins A., Frenk C., White S.D.M., SchuekerP.,CaldwellR.R.,BoehringerH.,CollinsC.A.&GuzzoL., Springel V., Stadel J. & Quinn T., 2002, MNRAS (submitted) astro-ph/0211480 andA&A(inpress) astro–ph/0201544 SpringelV.,YoshidaN.&WhiteS.D.M.,2001NewA,6,79S PressW.H.&Schechter P.,1974, ApJ,187,425 VianaP.,NicholR.&LiddleA.,2002, Apj,569L,75,V ShethR.K.,MoH.J.&TormenG.,2001MNRAS,323,1 Zel’dovich,Ya.B.,1970A&A,5,84

Mass of clusters in simulations

Similar Mass of clusters in simulations