Table Of Content

LA-UR-06-0847 DRAFTVERSIONFEBRUARY5,2008 PreprinttypesetusingLATEXstyleemulateapjv.26/01/00 THEHALOMASSFUNCTION:HIGH–REDSHIFTEVOLUTIONANDUNIVERSALITY ZARIJALUKIC´1,KATRINHEITMANN2,SALMANHABIB3,SERGEIBASHINSKY3,ANDPAULM.RICKER1,4 1Dept.ofAstronomy,UniversityofIllinois,Urbana,IL61801 2ISR-1,ISRDivision,LosAlamosNationalLaboratory,LosAlamos,NM87545 3T-8,TheoreticalDivision,LosAlamosNationalLaboratory,LosAlamos,NM87545 4NationalCenterforSupercomputingApplications,Urbana,IL61801 DraftversionFebruary5,2008 ABSTRACT 8 WestudytheformationofdarkmatterhalosintheconcordanceΛCDMmodeloverawiderangeofredshifts, 0 fromz=20tothepresent.Ourprimaryfocusisthehalomassfunction,akeyprobeofcosmology.Byperforming 0 alargesuiteofnested-boxN-bodysimulationswithcarefulconvergenceanderrorcontrols(60simulationswith 2 box sizes from 4 to 256h- 1Mpc), we determine the mass function and its evolution with excellent statistical and systematic errors, reachinga few percentover most of the consideredredshiftand mass range. Across the n a studiedredshifts,thehalomassisprobedover6ordersofmagnitude(107–1013.5h- 1M⊙). Historically,therehas J beenconsiderablevariationinthehighredshiftmassfunctionasobtainedbydifferentgroups. We havemadea 4 concertedefforttoidentifyandcorrectpossiblesystematicerrorsincomputingthemassfunctionathigh–redshift 1 and to explain the discrepanciesbetween some of the previousresults. We discuss convergencecriteria for the required force resolution, simulation box size, halo mass range, initial and final redshifts, and time stepping. 2 Becauseofconservativecutsonthemassrangeprobedbyindividualboxes,ourresultsarerelativelyinsensitive v tosimulationvolume,theremainingsensitivitybeingconsistentwithextendedPress-Schechtertheory.Previously 0 obtainedmassfunctionfitsnearz=0,whenscaledbylineartheory,areingoodagreementwithourresultsatall 6 redshifts,althoughamildredshiftdependenceconsistentwiththatfoundbyReedetal. mayexistatlowredshifts. 3 Overall,ourresultsareconsistentwitha“universal”formforthemassfunctionathighredshifts. 2 0 Subjectheadings:methods:N-bodysimulations—cosmology:halomassfunction 7 0 1. INTRODUCTION differentfitting formulaeand their origin in §2). Heitmann et / h al.(2006a)showthatthePress-Schechterformcanbeseverely Abroadsuiteofastrophysicalandcosmologicalobservations p incorrectathighredshifts: atz 10,thepredictedmassfunc- - providescompellingevidencefortheexistenceofdarkmatter. tionsinksbelowthenumericalr≥esultsbyanorderofmagnitude o Although its ultimate nature is unknown, the large-scale dy- r namics of dark matter is essentially that of a self-gravitating atthe upperend ofthe relevantmassscale. Consequently,in- t correct, or at best imprecise, predictions for the reionization s collisionless fluid. In an expandinguniverse, gravitationalin- a historycanresultfromthefailureoffittingformulae. stability leads to the formation and growth of structure in the : Sincehaloformationisacomplicatednonlineargravitational v dark matter distribution. The existence of localized, highly process,thecurrenttheoreticalunderstandingofthemass,spa- i overdensedarkmatterclumps, orhalos, is a keypredictionof X tial distribution, and inner profiles of halos remains at a rela- cosmologicalnonlineargravitationalcollapse. Thedistribution r tivelycrudelevel. Numericalsimulationsarethereforecrucial a of dark matter halo masses is termed the halo mass function as driversof theoreticalprogress, havingbeen instrumentalin and constitutes one of the most important probes of cosmol- obtaining important results such as the Navarro-Frenk-White ogy. At low redshifts, z 2, the mass function at the high- ≤ (NFW)profile(Navarroetal. 1997)fordarkmatterhalosand massend(clusterscales)is verysensitivetovariationsin cos- an(approximate)universalformforthemassfunction(Jenkins mological parameters, such as the matter content of the Uni- verse Ω , the dark energy content along with its equation-of- et al. 2001, hereafter Jenkins). In order to better understand m the evolutionof the mass functionat high redshifts, a number state parameter,w(Holderetal. 2001), andthenormalization ofnumericalstudieshavebeencarriedout.High–redshiftsimu- of the primordial fluctuation power spectrum, σ . At higher 8 lations,however,sufferfromtheirownsetofsystematicissues, redshifts,thehalomassfunctionisimportantinprobingquasar andsimulationresultscanbeatconsiderablevariancewitheach abundanceandformationsites(Haiman&Loeb 2001),aswell other,differingon occasionbyasmuchasan orderofmagni- as the reionization history of the Universe (Furlanettoetal. tude! 2006). Motivatedby allof these reasonswe have carriedouta nu- Manyrecentlysuggestedreionizationscenariosarebasedon mericalinvestigationoftheevolutionofthemassfunctionwith theassumptionthatthemassfunctionisgivenreliablybymod- theaimofattaininggoodcontroloverbothstatisticaland,more ified Press-Schechter type fits (Press & Schechter 1974, here- importantly,possiblesystematicerrorsinN-bodysimulations. after PS; Bondet al. 1991). However,the theoreticalbasis of OurfirstresultshavebeenreportedincondensedforminHeit- thisapproachisatbestheuristicandcarefulnumericalstudies mann et al. (2006a). Here we provide a more detailed and arerequiredinordertoobtainaccurateresults. Twoexamples completeexpositionofourwork,includingseveralnewresults. serve to illustrate this statement. Reed et al. (2003) report a Wefirstpayattentiontosimulationcriteriaforobtainingac- discrepancywiththeSheth-Tormenfit(Sheth&Tormen1999, curate mass functionswith the aim of reducingsystematic ef- hereafter ST) of 50% at a redshiftof z=15 (we explain the ∼ fects.Ourtwomostsignificantpointsarethatsimulationsmust 1 2 TheHaloMassFunction be started early enoughto obtain accurate results and that the cusstheimportanceofpost-processingcorrectionssuchasFOF box sizes must be large enough to suppress finite-volume ar- particlesamplingcompensationandfinite-volumeeffects. We tifacts. As in most recent work following that of Jenkins, we discussourresultsandconcludein§6. definehalomassesusingafriends-of-friends(FOF)halofinder withlinkinglengthb=0.2. Thischoiceintroducessystematic 2. DEFINITIONSANDPREVIOUSWORK issuesofitsown(e.g.,connectiontosphericaloverdensitymass Themassfunctiondescribesthenumberdensityofhalosof asafunctionofredshift),whichwetouchonasrelevantbelow. a givenmass. Inorderto determinethemassfunctionin sim- Asitisnotquantitativelysignificantinthecontextofthispaper, ulationsonehastofirstidentifythehalosandthendefinetheir weleaveadetaileddiscussiontolaterwork(Z.Lukic´ etal.,in mass. No precise theoreticalbasis exists for these operations. preparation;seealsoReedetal. 2007). Nevertheless,dependingonthe situationathand,the observa- Themoredetailedresultsinthispaperenableustostudythe tional and numerical communities have adopted a few “stan- massfunctionatstatisticalandsystematicaccuraciesreaching dard”waysofdefininghalosandtheirassociatedmasses. For afewpercentovermostofourredshiftrange,asubstantialim- arecentreviewoftheseissueswithregardtoobservations,see, provementovermostpreviouswork. Atthislevelwefinddis- e.g., Voit (2005),butforamoretheoreticallyorientedreview, crepancies with the “universal” fit of Jenkins at low redshifts see,e.g.,White (2001). (z<5),butitmustbekeptinmindthattheuniversalityofthe originalfitwasonlymeanttobeatthe 20%level. Recently, 2.1. HaloMass ± Reedetal. (2007)havefoundviolationofuniversalityathigh There are basically two ways to find halos in a simulation. redshifts(uptoz=30).Tofitthemassfunctiontheyhaveincor- One,theoverdensitymethod,isbasedonidentifyingoverdense porated an additionalfree parameter, the effective spectral in- regionsaboveacertainthreshold.Thethresholdcanbesetwith dexn ,withtheaimofunderstandingandtakingintoaccount eff respecttothecriticaldensityρ =3H2/8πG(orthebackground theextraredshiftdependencemissingfromconventionalmass– c densityρ =Ω ρ , whereΩ isthematterdensityoftheUni- function–fittingformulae.Oursimulationresultsareconsistent b m c m with the trends found by Reedetal. (2007) at low redshifts verse including dark matter and baryons). The mass M∆ of a halo identified this way is defined as the mass enclosed in (z 5),butathigherredshiftswedonotobserveastatistically sig≤nificantviolationoftheuniversalformofthemassfunction. a sphere of radius r∆ whose mean density is ∆ρc. Common values for ∆ range from 100 to 500 (or even higher). As ex- Resultsfromsomeprevioussimulationshavereportedgood plained in Voit (2005), cluster observers prefer higher values agreementwiththePress-Schechtermassfunctionathighred- for∆. Propertiesofclustersareeasiertoobserveinhigherden- shifts. Since the Press-Schechter fit has been found signifi- sity regionsand these regionsare more relaxed than the outer cantly discrepantwith low–redshiftresults (z<5), thiswould partswhicharesubjecttotheeffectsofinflowandincomplete implyastrongdisagreementwithextendingthewell-validated mixing. The disadvantageofdefininga haloin thismanneris low–redshiftnotionof(approximate)massfunctionuniversal- that sphericity of halos is implied, an assumption which may ity to high z. Our conclusionis thatthe simulationson which beeasilyviolated,e.g.,inthecaseofhalosthatformedinare- thesefindingswerebasedviolatedoneormoreofthecriteriato centmergereventorhalosathighredshifts.Athigherredshifts, bediscussedbelow. thenonlinearmassscaleM decreasesrapidly,andtheratioof As simulationsare perforcerestricted to finite volumes, the ∗ the consideredhalomass M to M can becomelarge. This obtainedmassfunctionclearlycannotrepresentthatofaninfi- halo ∗ translates into producinglarge-scale structures roughly analo- nitebox.Notonlyissamplingakeyissue,butalsothefactthat goustosuperclusterstructurestoday.Whilethesestructuresare simulationswith periodicboundaryconditionshavenofluctu- gravitationallybound, theyare often notvirialized, norspher- ationsonscaleslargerthantheboxsize. Tominimizeandtest ical. Even the much smaller structures (which are considered for these effects we were conservative in our choices of box inthispaper)arenotvirializedathighredshifts,andtherefore, size and the mass range probed in each individual box. We assumptions about sphericity are most likely violated. Hence also used nested-volumesimulationsto directlytest for finite- the sphericaloverdensitymethoddoesnotsuggest itself as an volume effects. Because we used multiple boxes and aver- obviouswaytoidentifyhalosathighredshift. aged mass function results over the box ensemble, extended The other method, the FOF algorithm, is based on finding Press-Schechtertheorycanbeusedtocorrectforresidualfinite neighborsofparticlesandneighborsofneighborsasdefinedby volume–effects (Mo&White 1996; Barkana&Loeb 2004); agivenseparationdistance(see,e.g.,Einastoetal.1984;Davis this approach is different from the individual box corrections et al. 1985). The FOF algorithmleads to halos with arbitrary appliedbyReedetal. (2007). Detailsaregivenin§5.3. shapessince nopriorsymmetryassumptionshavebeenmade. The paper is organized as follows. In §2 we give a brief Thehalomassisdefinedsimplyasthesumofparticleswhich overviewofthemassfunctionandpopularfittingformulae,dis- aremembersofthehalo. Whilethisdefinitioniseasytoapply cussingaswellpreviousnumericalworkonthehalomassfunc- tosimulations,theconnectiontoobservationsisdifficulttoes- tion athighredshifts. In §3we givea shortdescriptionof the N-bodycodeMC2(Mesh-basedCosmologyCode)andasum- tablish directly. (For an investigationof connections between differentdefinitions of halos masses and approximate conver- maryoftheperformedsimulations.In§4wederiveanddiscuss sionsbetweenthem,seeWhite2001). somesimple criteriaforthestartingredshiftandconsidersys- Itisimportanttokeepinmindthatthedefinitionofahalois tematic errorsrelatedto the numericalevolutionsuch as mass essentiallytheadoptionofsomesortofconventionforthehalo and force resolution and time stepping. These considerations boundary. In reality, a sharp distinction between the particles inturnspecifytheinputparametersforthesimulationsinorder ina haloandparticlesinthe simulation“field”doesnotexist. to span the desired mass and redshift range for our investiga- Jenkinsshowedthatthe choiceofa FOF finderwith a linking tion. In§5wepresentresultsforthemassfunctionatdifferent lengthb=0.2to definehalo massesprovidesthe bestfitfora redshiftsaswellasthehalogrowthfunction.Herewealsodis- universalformofthemassfunction.Thischoicehassincebeen Lukic´,Heitmann,Habib,Bashinsky,Ricker 3 adoptedbymanynumericalpractitionersasastandardconven- Even in linear theory, equation (8) is only an approximation. Ausefuldiscussionofthevarioushalodefinitionscanbe tion because baryons began their gravitational collapse with foundinWhite (2002). velocities different from those of CDM particles. Until re- InthispaperweusetheFOFalgorithmtoidentifyhalosand combination at z 1100, well into the matter era with non- ∼ theirmasses. ItwasrecentlypointedoutbyWarrenetal.(2006, negligiblegrowthofCDMinhomogeneities,thebaryonswere hereafter Warren) that FOF masses suffer from a systematic heldagainstcollapsebythepressureoftheCMBphotons(see, problem when halos are sampled by relatively small numbers e.g. Hu&Sugiyama (1996)). While thereafter the relative ofparticles. Althoughhaloscanberobustlyidentifiedwith as baryon-CDMvelocitydecayedas1/a,theresidualvelocitydif- few as 20 particles, if a given halo has too few particles, its ferencewassufficienttoaffectthegrowthfunctiond(z)atz=50 FOFmassturnsouttobesystematicallytoohigh. Wedescribe by morethan 1%and at z=10 by about0.2%(Yoshidaetal. howwecompensateforthiseffectin§5.2. Inthecurrentpaper, 2003;Naoz&Barkana 2007). all results for the mass function are displayed at a fixed FOF linkinglengthofb=0.2,usingtheWarrencorrection. 2.3. FittingFunctions Overthelastthreedecadesseveraldifferentfittingformsfor 2.2. DefiningtheMassFunction the mass function have been suggested. The mass function is Theexactdefinitionofthemassfunction,e.g.,integratedver- notonlyasensitivemeasureofcosmologicalparametersbyit- sus differential form or count versus number density, varies selfbutalsoakeyingredientinanalyticandsemianalyticmod- widely in the literature. To characterize different fits, Jenk- eling of the dark matter distribution, as well as of several as- ins introduced the scaled differential mass function f(σ,z) as pectsof the formation,evolution, anddistributionof galaxies. afractionofthetotalmassperlnσ- 1thatbelongstohalos: Therefore, if a reliable and accurate fit for the mass function applicabletoawiderangeofcosmologiesandredshiftswereto dρ/ρ M dn(M,z) f(σ,z) b = . (1) exist, it would be of obviousutility. In this section we briefly ≡ dlnσ- 1 ρ (z)dln[σ- 1(M,z)] b reviewthecommonfittingfunctionsandcomparethematdif- Heren(M,z)isthenumberdensityofhaloswithmassM,ρ (z) b ferentredshifts. isthebackgrounddensityatredshiftz,andσ(M,z)isthevari- Thefirstanalyticmodelforthemassfunctionwasdeveloped anceofthelineardensityfield. AspointedoutbyJenkins,this by PS. Their theoryaccountsfor a sphericaloverdenseregion definitionofthemassfunctionhastheadvantagethattoagood in an otherwise smooth background density field, which then accuracyitdoesnotexplicitlydependonredshift,powerspec- evolvesasaFriedmannuniversewithapositivecurvature. Ini- trum, or cosmology; all of these are encapsulated in σ(M,z). tially, the overdensity expands, but at a slower rate than the Forthemostpart,wewilldisplaythemassfunction backgrounduniverse(thusenhancingthedensitycontrast),un- dn tilitreachesthe‘turnaround’density,afterwhichcollapsebe- F(M,z) (2) ≡ dlogM gins. Althoughfromapurelygravitationalstandpointthiscol- lapseendswithasingularity,itisassumedthatinreality–due as a function of logM itself. [In §5 we include results for to the sphericalsymmetrynotbeing exact– the overdensere- f(σ,z).] gionwillvirialize. ForanEinstein-deSitteruniverse,theden- To compute σ(M,z), the power spectrum P(k) is smoothed sity of suchan overdenseregionatthe virializationredshiftis with a spherical top-hat filter function of radius R, which on averageenclosesamassM(R=[3M/4πρ (z)]1/3): z 180ρc(z). Atthispoint,thedensitycontrastfromthelinear b ≈ theory of perturbation growth [δ(~x,z) = d(z)δ(~x,0)] would be σ2(M,z)= d2(z) ∞k2P(k)W2(k,M)dk, (3) δc(z) 1.686inanEinstein-deSittercosmology. ForΩm<1, 2π2 Z0 theva≈lueofthe thresholdparameterδc canvary(seeLacey& whereW(k,M)isthetop-hatfilter: Cole1993),butthedependenceoncosmologyhaslittle quantitative significance (see, e.g., Jenkins). Thus, throughoutthis 3 ,r<R W(r)= 4πR3 (4) paperweadoptδ =1.686. 0,r>R c (cid:26) Followingtheabovereasoningandwiththeassumptionthat W(k)= 3 [sin(kR)- kRcos(kR)]. (5) the initial density perturbations are described by a homoge- (kR)3 neousandisotropicGaussianrandomfield, thePSmassfunc- Theredshiftdependenceentersonlythroughthegrowthfactor tionisspecifiedby d(z),normalizedsothatd(0)=1: 2δ δ2 f (σ)= cexp - c . (9) σ(M,z)=σ(M,0)d(z). (6) PS π σ 2σ2 r (cid:18) (cid:19) In theapproximationofnegligibledifferencein the CDM and The PS approachassumes that all mass is inside halos, as en- baryonpeculiarvelocities,thegrowthfunctioninaΛCDMuni- forcedbytheconstraint verseisgivenby(Peebles1980) +∞ D+(a) fPS(σ)dlnσ- 1=1. (10) d(a)= D+(a=1), (7) Z- ∞ While as a first rough approximation the PS mass function where we consider d as a function of the cosmological scale agreeswithsimulationsatz=0reasonablywell,itoverpredicts factora=1/(1+z),and thenumberoflow–masshalosandunderpredictsthenumberof 5Ω H(a) a da′ massive halos at the currentepoch. Furthermore, it is signifi- D+(a)= 2m H [a′H(a′)/H ]3 (8) cantlyinerrorathighredshifts(see,e.g.,Springeletal.2005; 0 Z0 0 Heitmannetal.2006a;§5.4). withH(a)/H = Ω /a3+(1- Ω ) 1/2. Inparticular,forz 1, AfterPS,severalsuggestionsweremadeinordertoimprove 0 m m whenmatterdominatesthecosmologicalconstant,D+(a) ≫a. the mass function fit. These suggestions were based on more (cid:2) (cid:3) ≃ 4 TheHaloMassFunction TABLE1 MASSFUNCTIONFITSFOR f(σ) Reference FittingFunction f(σ) MassRange Redshiftrange ST,Sheth&Tormen(2001) 0.3222 2aδcexp - aδc2 1+ σ2 p unspecified unspecified π σ 2σ2 aδ2 c Jenkins 0.315eqxp - lnσ- 1h+0.6i1h3.8 (cid:16) (cid:17) i - 1.2 lnσ- 1 1.05 z=0- 5 | | ≤ ≤ Reedetal.(2003) fST(σ)exp(cid:2) - 0.7/ σ(cosh(2(cid:3)σ))5 - 1.7 lnσ- 1 0.9 z=0- 15 ≤ ≤ Warren 0.7234 σ- 1(cid:8).625+0.(cid:2)2538 exp - 1(cid:3).1σ(cid:9)9282 (1010- 1015)h- 1M⊙ z=0 p Reedetal.(2007) A 2πa (cid:0)1+ aσδ2c2 +0.6G(cid:1)1(σ)(cid:2)+0.4G2(cid:3)(σ) - 0.5≤lnσ- 1≤1.2 z=0- 30 qδcexhp - c(cid:16)aδc2 -(cid:17) 0.03 δc 0.6 i ×σ 2σ2 (neff+3)2 σ h (cid:0) (cid:1) i Note. —Shownareexamplesofcommonlyusedfittingfunctions. STuseda=0.707andp=0.3,whileSheth&Tormen (2002)suggestthata=0.75leadstoa betterfit. TheWarrenfitrepresentsbyfarthelargestuniformsetofsimulationsbasedonmultipleboxeswiththesamecosmologyrunwiththesamecode. Weuse itasareferencestandardthroughoutthispaper. Reedetal. (2003)suggestanempiricaladjustementoftheSTfit,whichisslightlymodifiedinReedetal. (2007). Forthelatter,G1(σ)andG2(σ)aregivenbyeqs.(16)and(17),respectively,c=1.08,ca=0.764,andA=0.3222. refined dynamicalmodeling, direct fitting to simulations, or a 1.2 Warren et al. combinationofthetwo. 1 Jenkins et al. Press-Schechter 0.8 z=0 Using empiricalargumentsST proposedan improvedmass Sheth-Tormen 0.6 functionfitoftheform: 9 10 11 12 13 14 15 f (σ)=0.3222 2aδcexp - aδc2 1+ σ2 p , (11) on 1.25 ST π σ 2σ2 aδ2 cti 1 r (cid:18) (cid:19)(cid:20) (cid:18) c(cid:19) (cid:21) un 0.5 z=5 awi=th0.a75=a0s.7a0n7imanpdropve=d0v.a3l.ue(.S)hSehtheth&eTtoarlm. (e2n020010)2resduegrigveesdt Mass F 049 10 11 12 13 this fit theoreticallyby extendingthe PS approachto an ellip- en 3 tical collapse model. In this model, the collapse of a region Warr 2 z=10 droeupnednidnsgnsohteoanrlfiyeoldn.iTtsheinditeipalenodveenrdceenissitcyhobsuetnalssuochonthtahteitsurer-- nction/ 017 8 9 10 11 12 u 6 lcionveearrsrtehgeimZee.l’dAovhicahloaipspcrooxnismidaetrieodnv(Zireial’ldizoevdicwhh1e9n7t0h)einthtihrde Mass F 435 z=15 axiscollapses(seealsoLee&Shandarin(1998)foranearlier, X 21 0 differentapproachtothesameidea). 7 8 9 10 11 Jenkinscombinedhighresolutionsimulationsforfourdiffer- 6 5 ent CDM cosmologies(τCDM, SCDM, ΛCDM, and OCDM) 4 3 z=20 spanningamassrangeofover3ordersofmagnitude( (1012- 2 1015)h- 1M ), and including several redshifts betwe∼en z = 5 01 ⊙ 7 7.5 8 8.5 9 9.5 10 and 0. Independentof the underlying cosmology, the follow- log(M [M./h]) O ingfitprovidedagoodrepresentationoftheirnumericalresults FIG.1.—RatiooftheJenkins,PS,andSTmassfunctionfitswithrespectto (within 20%): theWarrenfitforfivedifferentredshiftsoverarangeofhalomasses. Topto ± f (σ)=0.315exp - lnσ- 1+0.613.8 . (12) bottom:Redshiftsz=0,5,10,15,and20.Notethattherangesoftheaxesare Jenkins | | differentinthedifferentpanels. WedonotshowtheJenkinsfitbelowmasses Theaboveformulaisveryclose(cid:0)totheSheth-Torm(cid:1)enfit,lead- of1011h- 1M⊙atz=0,sinceitisnotvalidforsuchlowmassesatthatredshift. ingtosomeimprovementatthehigh-massend. Thedisadvan- tage is thatit cannotbe simply extrapolatedbeyondthe range ofthefit,sinceitwastunedtoaspecificsetofsimulations. respecttotheWarrenfitinFigure1. WedonotshowtheJenk- By performing 16 nested-volume dark matter simulations, insfitbelow1011h- 1M⊙atz=0sinceitdivergesinthisregime. Warren was able to obtain significant halo statistics spanning The originalST fit, the Jenkins fit, and the Warren fit all give amassrangeof5ordersofmagnitude( (1010- 1015)h- 1M ). similarpredictions. ThediscrepancybetweenPSandtheother ⊙ ∼ Becausethisrepresentsbyfarthelargestuniformsetofsimulations– fitsbecomesmoresevereforhighermassesathighredshifts.PS basedonmultipleboxeswiththesamecosmologyrunwiththe dramaticallyunderpredictshalosinthehigh-massrangeathigh samecode–weuseitasareferencestandardthroughoutthispa- redshifts(assumingthattheotherfitsleadtoreasonableresults per. UsingafunctionalformsimilartoST,Warrendetermined inthisregime). Forlow-masshalosthedisagreementbecomes thebestmassfunctionfittobe less severe. For z=0 the Warren fit agrees, especially in the fWarren(σ)=0.7234 σ- 1.625+0.2538 exp - 1.1σ9282 . (13) SloTwfi-mt.aAssttrhaenhgiegbhe-mloawss1e0n1d3hth- 1eMdi⊙ff,etroenbceettienrctrheaanse5s%upwtoit2h0t%he. (cid:18) (cid:19) TheJenkinsfitleadstosimilarresultsovertheconsideredmass For a quantitative com(cid:0) parison of the(cid:1)different fits at different range.Athigherredshiftsandintermediate-massrangesaround redshifts,weshowtheratioofthePS,Jenkins,andSTfitswith Lukic´,Heitmann,Habib,Bashinsky,Ricker 5 109h- 1M ,theWarrenandSTfitdisagreebyroughlyafactor 1000 ⊙ of2. (107 - 108) M./h O SeveralothergroupshavesuggestedmodificationsoftheST fit. In§5wecompareourresultswithtwoofthem.Reedetal. 1 (108 - 109) M./h (2003)suggestanempiricaladjustmenttotheST fitbymulti- O plyingitwithanexponentialfunction,leadingto (109 - 1010) M./h f (σ)= f (σ)exp - 0.7/ σ(cosh(2σ))5 , (14) 0.001 O Reed03 ST ] 3 tvoalaidsuopvperretshseiornanogfeth-e1S.7T≤fitl(cid:8)antσla- 1rg≤e(cid:2)σ0-.91..ITnhRiseaeddjue(cid:3)st(cid:9)taml.en(t2l0ea0d7s) -pc/h) 1e-06 (1010 - 1011) MO./h the adjustmenttotheST fitisslightlymodifiedagain,leading M tothefollowingnewfit: [( n 2a σ2 p 1e-09 fReed07(σ)=A π 1+ aδ2 +0.6G1+0.4G2 (1011 - 1012) MO./h r (cid:20) (cid:18) c(cid:19) (cid:21) ×δσcexp"- c2aσδ2c2- (ne0ff.0+33)2(cid:18)δσc(cid:19)0.6#, (15) 1e-12 (1014 - 1015) M.(/1h013 - 1014()1 0M12O. /-h 1013) MO./h ln(σ- 1- 0.4)2 O G =exp - , (16) 1e-15 1 2(0.6)2 0 5 10 15 20 (cid:20) (cid:21) Redshift ln(σ- 1- 0.75)2 G =exp - , (17) FIG. 2.—HalogrowthfunctionbasedontheWarrenmassfunctionfitfor 2 2(0.2)2 differentmassbins. Thecurvesforthelowermassbinshaveamaximumat (cid:20) (cid:21) z>0whichreflectsacrossoverofthemassfunctionsatdifferentredshifts. withc=1.08,ca=0.764,andA=0.3222. Theadjustmenthas very similar effects to that of Reedetal. (2003), as we show in §5. Reedetal. (2007) notethatthe(small)suppressionof eightdifferentmassbins,coveringthemassrangeinvestigated the mass functionrelativeto ST as a functionof redshiftseen inthispaper,asafunctionofredshiftz. Asexpectedfromthe insimulations(seealsoHeitmannetal. 2006a)canbetreated paradigmof hierarchicalstructure formationin a ΛCDM cos- by adding an extra parameter, the power spectral slope at the mology, small halos form much earlier than larger ones. An scaleofthehaloradius,n (formallydefinedbyequation(42) eff interesting feature in the lower mass bins is that they have a below). Wereturntothisissuewhenwediscussournumerical maximumatdifferentredshifts.Thenumberofthesmallestha- results in §5. We summarize the described, most commonly losgrowsuntilaredshiftof z=2andthendeclineswhenhalos usedfittingfunctionsinTable1. startmergingandformingmuchmoremassivehalos. Thisfea- Although fitting functionsmay be a useful way to approxi- tureisreflectedinacrossingofthemassfunctionsatdifferent matelyencapsulateresultsfromsimulations,meaningfulcom- redshiftsforsmallhalos. parisonstoobservationsrequireovercomingmanyhurdles,e.g., anoperationalunderstandingofthedefinitionofhalomass(see, 2.5. MassFunctionatHighRedshift:PreviousWork e.g., White 2001), how it relates to various observations, and error control in N-body codes (see, e.g., O’Shea et al. 2005; Mostofthe effortto characterize,fit, andevaluatethemass Heitmannetal.2005). Inthispaper,ourfocusisfirstoniden- functionfromsimulationshasbeenfocusedonornearthecur- tifyingpossiblesystematicproblemsintheN-bodysimulations rentcosmologicalepoch,z 0.Thisismainlyfortworeasons: ∼ themselvesandhowtheycanbeavoidedandcontrolled. (1)sofarmostobservationalconstraintshavebeenderivedfrom low-redshiftobjects(z<1);(2)theaccuratenumericalevalua- 2.4. HaloGrowthFunction tionofthemassfunctionathighredshiftsisanontrivialtask. The increasing reach of telescopes on the ground and in Ausefulwaytostudythestatisticalevolutionofhalomasses space,suchastheupcomingJamesWebbSpaceTelescope,al- in simulations is to transform the mass function into the halo growthfunction,n(M ,M ,z) M2FdlogM(Heitmannetal. lows us to study the Universe at higher and higher redshifts. 1 2 ≡ M1 Recentdiscoveriesinclude970galaxiesatredshiftsbetweenz= 2006a),whichmeasuresthemass-binnednumberdensityofha- R 1.5andz=5fromtheVIMOSVLTDeepSurvey(LeFevreetal. losasafunctionofredshift. Thehalogrowthfunction,plotted 2005),andtherecentobservationofagalaxyatz=6.5(Mobasheretal. versusredshiftinFigure2,showsataglancehowmanyhalos 2005). The epoch of reionization (EOR) is of central impor- in a particular mass bin and box volume are expectedto exist tance to the formationof cosmic structure. Althoughourcur- atacertainredshift. Thishelpssettherequiredmassandforce rentobservationalknowledgeoftheEORisratherlimited,fu- resolutionin a simulation whichaimsto capturehalosathigh ture 21 cm experiments have the potential for revolutionizing redshifts.Foragivensimulationvolume,thehalogrowthfunc- thefield.Proposedlow-frequencyradiotelescopesincludeLO- tion directly predicts the formationtime of the first halos in a FAR(LowFrequencyArray)1, the MileuraWide Field Array givenmassrange. (MWA)(Bowmanetal. 2006)2, andthenext-generationSKA In orderto derivethis quantityapproximately,we first con- (Square Kilometer Array) 3. The observationalprogressis an vertan accuratemassfunctionfit(we use the Warrenfithere) importantdriverforhigh-redshiftmassfunctionstudies. into a functionof redshiftz. Ithasbeenshownrecentlybyus (Heitmann et al. 2006a) that mass function fits work reliably 1Seehttp://www.lofar.org enoughouttoatleastz=20,andcanthereforebeusedtoesti- 2Seehttp://haystack.mit.edu/arrays/MWA/ matethehalogrowthfunction.Figure2showstheevolutionof 3Seehttp://www.skatelescope.org 6 TheHaloMassFunction TABLE2 SUMMARYOFTHEPERFORMEDRUNS BoxSize Resolution ParticleMass SmallestHalo Mesh z z No.ofRealizations (h- 1Mpc) (h- 1kpc) in final (h- 1M ) (h- 1M ) ⊙ ⊙ 10243 256 250 100 0 8.35 1010 3.34 1012 5 10243 128 125 200 0 1.04×1010 4.18×1011 5 10243 64 62.5 200 0 1.31× 109 5.22×1010 5 10243 32 31.25 150 5 1.63×108 6.52× 109 5 10243 16 15.63 200 5 2.04×107 8.16×108 5 10243 8 7.81 250 10 2.55×106 1.02×108 20 10243 4 3.91 500 10 3.19×105 1.27×107 15 × × Note.—Massandforceresolutionsofthedifferentruns.Thesmallesthalosweconsidercontain40particles.Allsimulationshave2563particles. Theoreticalstudiesofthemassfunctionathighredshiftsare lyzealargesuiteofN-bodysimulationswithvaryingboxsizes challengingduetothesmallmassesofthehalosatearlytimes. between4and256h- 1Mpc,includingmanyrealizationsofthe Inordertocapturethesesmall-masshalos,highmassandforce smallboxes,tostudythemassfunctionatredshiftsuptoz=20 resolutionarebothrequired. Forthelargesimulationvolumes andtocoveralargemassrangebetween107and1013.5h- 1M . ⊙ typical in cosmological studies, this necessitates a very large Withrespecttoourpreviouswork,thenumberofsmall-boxre- numberofparticles,aswellasveryhighforceresolution.Such alizations has been increased to improve the statistics at high simulationsareverycostly,andonlyaverylimitednumbercan redshifts. Our results categorically rule out the PS fit as be- beperformed,disallowingexplorationofawiderangeofpossi- ing more accurate than any of the more modern forms at any blesimulationparameters.Alternatively,manysmallervolume redshiftuptoz=20,thediscrepancyincreasingwithredshift. simulationboxes, each with moderateparticle loading,can be 3. THECODEANDTHESIMULATIONS employed.Thisleadsautomaticallytohighforceandmassres- olutioningridcodes(suchasparticle-mesh[PM])andalsore- Allsimulationsinthispaperarecarriedoutwiththeparallel duces the costs for achieving sufficient resolution for particle PMcodeMC2. ThiscodesolvestheVlasov-Poissonequations codes (such as tree codes) or hybridcodes (such as TreePM). foranexpandinguniverse.Itusesstandardmassdepositionand The disadvantages of this strategy are the limited statistics in force interpolationmethods allowing periodic or open bound- individualrealizations(because fewer halos form in a smaller aryconditionswithsecond-order(global)symplectictimestep- box)andtheunreliabilityofsimulationsbelowanintermediate pingandfastfouriertransformbasedPoissonsolves. Particles redshiftatwhichthelargestmodeintheboxisstill(accurately) are depositedon the grid using the cloud-in-cellmethod. The linear. In addition, results from small boxes may be biased, overall computational scheme has proven to be accurate and sincetheyonlyfocusonasmallregionandvolume.Therefore, efficient: relatively large time steps are possible with excep- onemustshowthatthesimulationsarefreefromfinite-volume tionalenergyconservationbeing achieved. MC2 has been ex- artifacts, e.g.missing tidalforces, andruna sufficientnumber tensivelytestedagainststate-of-the-artcosmologicalsimulation of statistically independent simulations to reduce the sample codes(Heitmannetal.2005,2007). variance. Both strategies, employing large volume or multi- Weusethefollowingcosmologyforallsimulations: ple small-volume simulations, have been followed in the past Ω =1.0, Ω =0.253, Ω =0.048, tot CDM baryon in orderto obtainresultsat highredshifts. The differentmass σ =0.9, H =70kms- 1Mpc- 1, n=1, (18) rangesinvestigatedbydifferentgroupsare shownin Figure3. 8 0 Thefitsareshownforredshiftsz=10and20. IntheAppendix in concordancewith cosmic microwavebackgroundand large we providea very detailed discussion on previousfindings as scalestructureobservations(MacTavishetal. 2006)(thethird- organizedbysimulationvolume. yearWilkinsonMicrowaveAnisotropyProbeobservationssug- In summary, there is considerable variation in the high- gest a lower value of σ8; Spergeletal. (2007)). The transfer redshift (z>10) mass function as found by different groups, functionsaregeneratedwithCMBFAST(Seljak&Zaldarriaga independent of box size and simulation algorithm. Broadly 1996). We summarizethe differentruns, includingtheir force speaking,theresultsfallintotwoclasses: eitherconsistentwith and mass resolution, in Table 2. As mentioned earlier, we lineartheoryscalingofauniversalform(Jenkins,Reed,ST,or identify halos with a standard FOF halo finder with a linking Warren) at low redshift (Reed et al. 2003, 2007; Springel et lengthofb=0.2.DespiteseveralshortcomingsoftheFOFhalo al. 2005; Heitmann et al. 2006a; Maio et al. 2006; Zahn et finder, e.g., the tendencyto link up two haloswhich are close al. 2007) or more consistent with the PS fit (Jang-Condell& toeachother(see,e.g.,Gelb&Bertschinger1994,Summerset Hernquist2001;Yoshidaetal. 2003a,2003b,2003c;Cenetal. al.1995)orstatisticalbiases(Warren),theFOFalgorithmitself 2004;Ilievetal. 2006;Trac&Cen2006). iswelldefinedandveryfast. Asdiscussedin§2.1,weadoptthe Ouraimhereistodeterminetheevolutionofthemassfunc- correctionforsamplingbiasgivenbyWarrenwhenpresenting tion accurately, at the few percentlevel, and at the same time ourresults. characterize many of the numerical and physical factors that 4. INITIALCONDITIONSANDTIMEEVOLUTION controlthe errorinthe massfunction(detailsbelow). We fol- lowuponourpreviouswork(Heitmannetal. 2006a)andana- Inanear-idealsimulationwithveryhighmassandforceres- olution, the first halos would form very early. By z= 50, a Lukic´,Heitmann,Habib,Bashinsky,Ricker 7 1e+06 126 Mpc/h 32 Mpc/h 8e+05 8 Mpc/h es cl6e+05 arti p of er b m4e+05 u N 2e+05 0 0 10 20 30 40 50 60 70 80 90 100 110 120 |Gradient of f (q)|/D p FIG.4.—Probabilitydistributionof|∇φ|inunitsoftheinterparticlespacing ∆p.Allcurvesshownaredrawnfrom2563particlesimulationsfromaninitial density grid of 2563 zones. Thephysical box sizes are 126h- 1Mpc (black line),32h- 1Mpc(redline),and8h- 1Mpc(greenline). Asexpected, h|∇φ|i increases with decreasing box size (which is equivalent to increasing force resolution).Therefore,zinandzcrossarehigherforthesmallerboxes. havingmovedonlytherelativelysmalldistanceassignedbythe initialZel’dovichstep. Onlyaftertheparticleshavetraveleda sufficientdistanceandcomeclosetogethercantheyinteractlo- cally to form the first halos. In the followingwe estimate the redshiftwhentheZel’dovichgriddistortionequalstheinterpar- ticlespacing,leadingtothemostconservativeestimateforthe redshiftofpossiblefirsthaloformation.Fromthisestimate,we derivethenecessarycriterionforthestartingredshiftforagiven boxsizeandparticlenumber. 4.1. InitialRedshift In order to capture halos at high redshifts, we have found thatitisveryimportanttostartthesimulationsufficientlyearly. Weconsidertwocriteriaforsettingthestartingredshift:(1)en- FIG. 3.—Summaryofrecentworkonthemassfunctionathighredshift. suringthelinearityofallthemodesintheboxusedtosample Themassfunctionfitsareshownatz=10(top)andz=20(bottom)forthe theinitialmatterpowerspectrum,and(2)restrictingtheinitial cosmologyusedthroughoutthispaper(theothergroupsusedslightlydifferent particlemoveto preventinterparticlecrossingandto keep the parameters).Atz=10,Jang-Condell&Hernquist(2001)(grayshadedregion) covertheverylowmassrangeusingaverysmallbox,asdoCenetal.(2004) particle grid distortion relatively small. The first criterion is (greenshadedregion).ThelargerboxesofReedetal. (2007)andSpringelet commonlyusedtoidentifythestartingredshiftinsimulations. al. (2005)(redshadedregion)leadtoresultsathigherhalomasses.Notethat However,asshownbelow,itfailstoprovidesufficientaccuracy inthisregimethePSfitdeviatessubstantiallyfromtheotherfits,whileatthe of the mass functions, accuracy which can be obtained when verylowmassendallfitstendtomerge.Oursuiteofvariableboxsizescovers amassrangeof107to1013.5h- 1M⊙betweenz=0and20,amuchlargerrange a second (much more restrictive) control is applied. Further- thanpreviouslycoveredbyanygroupwithauniformsetofsimulations. At more,itisimportanttoallowasufficientnumberofexpansion z=20Yoshidaetal. (2003a,2003b,2003c,2003d)covertheverylowmass factorsbetweenthestartingredshiftz andthehighestredshift end of the mass function, while Zahn etal. (2007) investigate larger mass in ofphysicalsignificance. Thisisneededtomakesurethatarti- halos. Oursimulationsoverlapwithbothofthemattheedges. Bycombining aheterogeneoussetofsimulations,Reedetal. (2007)coverawiderangein factsfromtheZel’dovichapproximationarenegligibleandthat massandredshift.FigurequalityreducedforthearXivversionofthepaper. thememoryoftheartificialparticledistributionimposedatz in (gridorglass)islostbythe timeanyhalophysicsistobe ex- tractedfromthesimulationresults. Althoughnotstudiedhere,it isimportantto notethathigh- redshift commonly used to start cosmological simulations, a redshift starts do require the correct treatment of baryons as largenumberofsmallhaloswouldalreadybepresent(see,e.g., notedin §2.2. In addition,redshiftstartsthataretoohighcan Reed et al. [2005] for a discussion of the first generation of lead to forceerrorsfora varietyofreasons, e.g., interpolation star-forminghalos). Ina morerealistic situation,however,the systematics,round-off,andcorrelatederrorsintreecodes. initialconditionsatz=50haveofcoursenohalos,theparticles 8 TheHaloMassFunction TABLE3 140 INITIALREDSHIFTESTIMATESFROMTHELINEARITYOF 5 realizations for 2563 particles ∆2(kNy) Average for 2563 particles 120 (Bho- 1xMSpizce) (hMkNpyc- 1) T(z=0,kNy) zin ossing100 5A vreearlaigzea tfioorn s1 2fo8r3 1p2a8rt3i cplaersticles 13128266 16250.5030 1410...783.0··10110000--2- 566 43555305 hift of first cr 80 s · ed 60 e r g Note.—Thenumberofparticlesis2563,thesameinallsimulations. era v 40 A 4.1.1. InitialPerturbationAmplitude 20 The initial redshift in simulations is often determined from the requirement that all mode amplitudes in the box below 0 the particle Nyquist wavenumbercharacterizedby k /2 with 1 10 100 Ny k =2π/∆ , where ∆ is the mean interparticle spacing, be Box length [Mpc/h] Ny p p sufficientlylinear.Thesmallertheboxsizechosen(keepingthe 400 numberofparticlesfixed),thelargerthelargestk-value.There- fore,inordertoensurethatthesmallestinitialmodeinthebox 5 realizations for 2563 particles is wellin thelinearregime,thestartingredshiftmustincrease astheboxsizedecreases. Inthefollowingwegiveanestimate Average for 2563 particles basedonthiscriterionfortheinitialredshiftfordifferentsim- g300 5 realizations for 1283 particles n ulation boxes. We (conservatively) require the dimensionless ssi Average for 1283 particles powerspectrum∆2=k3P(k)/2π2tobesmallerthan0.01atthe cro initi∆al2r(ekdNsyh,izfitn.)T=hke3iPn(i2tkiπNaly2,pzoinw)e∼rsBpekcn2t+rπ3uT2m(2z(iksN+gyi1,vz)e2=nb0y), (19) hift of first 200 in s d w(sheee,ree.Bg.i,sBthuennno&rmWalhiiztaeti[o1n99o7f]thfeorparimfitotirndgiaflupnocwtioernsfpoercBtruinm- est re h cluding COBE results) and T(k) is the transfer function. We Hig100 assumethespectralindextoben=1,whichissufficienttoob- tain an estimate for the initialredshift. Fora ΛCDM universe the normalizationis roughlyB 3.4 106(h- 1Mpc)4. There- ∼ × fore,z issimplydeterminedby in zin≃4150kN2yT(z=0,kNy). (20) 0 1 Bo1x0 length [Mpc/h] 100 We present some estimates for different box sizes in Table 3. Forthesmallerboxes(<8h- 1Mpc),theestimatesfortheinitial redshiftsareataroundz =50. FIG. 5.—Average redshift offirst crossing (top) and highest redshift of in firstcrossing(bottom)asafunctionofboxsize. Theinitialconditions(five It is clear that this criterion simply sets a minimal require- different realizations) areshownforboxes between 1and512h- 1Mpcwith mentforzin andneglectsthe factthattheinitial particlemove 1283and2563particles. Foreachinitialcondition,zficrrostssandzrcmrosssareshown shouldbesmallenoughtomaintainthedynamicalaccuracyof bythecrosses. Thesolidlines showtheaveragefromthefiverealizations. perturbationtheory(linearorhigherorder)usedtosettheini- Asexpected,scatterfromthedifferentrealizationsislargerforsmallerboxes. Theseplotsprovideestimatesoftherequiredinitialredshiftforasimulation tial conditions. Also, this criterion certainly does not tell us since|∇φ|/∆pisz-independentintheZel’dovichapproximation(seetext). thatif,e.g.,z =50,thenwemayalreadytrustthemassfunc- in tionat,say,z=30.Anexampleofthisisprovidedbytheresults of Reedetal. (2003), who find that their high-redshiftresults 1970). Initiallyeachparticleisplacedonauniformgridorina betweenz=7and15havenotconvergeiftheystarttheirsim- glassconfigurationandisthengivenadisplacementdetermined ulations at z =69. (A value of z =139 was claimed to be in in bytherelation sufficientintheircase.) x=q- d(z) φ, (21) We now consider another criterion – ostensibly similar in ∇ spirit–thatparticlesshouldnotmovemorethanacertainfrac- whereqistheLagrangiancoordinateofeachparticle.Thegra- tion of the interparticlespacingin the initialization step. This dient of the potential φ is independent of the redshift z. The secondcriteriondemandsmuchhigherredshiftstarts. Zel’dovichapproximationholdsinthemildlynonlinearregime, aslongasparticletrajectoriesdonotcrosseachother(nocaus- 4.1.2. FirstCrossingTime ticshaveformed).Studyingthemagnitudeof φ allowsusto |∇ | estimatetwoimportantredshiftvalues: first,theinitialredshift In cosmologicalsimulations, initial conditions are most of- z at which the particles should not have moved on average tengeneratedusingtheZel’dovichapproximation(Zel’dovich in Lukic´,Heitmann,Habib,Bashinsky,Ricker 9 morethanafractionoftheinterparticlespacing∆ =L /n , TheZel’dovichapproximationmatchestheexactdensityand p box p whereL isthephysicalboxsizeandn thenumberofparti- velocity fields to linear order in Lagrangian perturbation the- box p cles in the simulation; second, the redshift at which particles ory. Therefore, there is in principle an error arising from the first move more than the interparticle spacing, z , i.e., at resultingdiscrepancywiththedensityandvelocityfieldsgiven cross whichtheyhavetraveledonaverageadistancegreaterthan∆ . bytheexactgrowingmodeinitializedinthefarpast. p Foragivenrealizationofthepowerspectrum,themagnitude This error is linear in the number of expansion factors be- of φ dependsontwoparameters: thephysicalboxsizeand tweenz andtheredshiftofinterestz . Ithasbeenexplored in phys |∇ | theinterparticlespacing. Togethertheseparametersdetermine in the context of simulation error by Valageas (2002) and by the range of scales under consideration. The smaller the box, Crocceetal. (2006). Depending on the quantity being cal- the smaller the scales; therefore, φ increases and both z culated,thenumberofexpansionfactorsbetweenz andz in in phys |∇ | and z increase. Increasing the resolution has the same ef- requiredtolimittheerrortosomegivenvaluemayormaynot cross fect. InFigure4weshowtheprobabilitydistributionfunction beeasytoestimate. Forexample,unlikequantitiessuchasthe for φ for three differentbox sizes, 8, 32, and 126h- 1Mpc, skewness of the density field, there is no analytical result for |∇ | representingvaluesstudied by other groups, as well as in this howthiserrorimpactsthedeterminationofthemassfunction. paper. Tomakethecomparisonbetweenthedifferentboxsizes Neither does there exist any independent means of validating morestraightforward,wehavescaled φ withrespecttothe the result aside from convergence studies. Nevertheless, it is interparticle spacing ∆ . All curves a|∇re d|rawn from simula- clearthattobeconservative,oneshouldaimforafactorof 20 p ∼ tionswith2563particlesona2563grid,inaccordancewiththe in expansion factor in order to anticipate errors at the several setupofourinitialconditions. Thebehavioroftheprobability percentlevel, aruleofthumbthathasbeenfollowedbymany function follows our expectations: the smaller the box, or the N-bodypractitioners(andoftenviolatedbyothers!). Thisrule highertheforceresolution,thelargertheinitialdisplacements ofthumbgivesredshiftstartsthatareroughlyinagreementwith oftheparticlesonaverage. Fromthemeanandmaximumval- the estimates in the previous subsection. Convergence tests uesofsuchadistributionwecandetermineappropriatevalues doneforoursimulationsshowthatthesuppressioninthemass forz andz . Forourestimatesweassumed(z) 1/(1+z), functionisverysmall(lessthan1%)forsimulationswhoseevo- in cross whichisvalidforhighredshifts. Themaximumand≃rmsinitial lutioncoversafactorof15intheexpansionfactorandcanbe displacementsoftheparticlescanthenbeeasilycalculated: up to 20% for simulations that evolved by only 5 expansion factors.However,duetomodestparticleloads,wewereunable max( φ/∆ ) δmax |∇ | p , (22) todistinguishbetweentheerrorinducedbytoofewexpansion in ≃ 1+zin factorsandthebreakdownoftheZel’dovichapproximation. δrms rms(|∇φ|/∆p). (23) Another possible problem, independent of the accuracy of in ≃ 1+z the Zel’dovich approximation, is the initial particle distribu- in tion itself. Whether based on a grid or a glass, the small- Theveryfirst“gridcrossing”ofaparticleoccurswhenδmax= in distance (k>k ) mass distribution is clearly not sampled at 1; onaveragetheparticleshavemovedmorethanoneparticle Ny all by the initial condition. Therefore, unlike the situation spacingwhenδrms=1. Thisleadstothefollowingestimates: in thatwouldariseifafullydynamicallycorrectinitialcondition zfirst max( φ/∆ )- 1, (24) were given, some time must elapse before the correct small- cross ≃ ∇ p zrms rms( φ/∆ )- 1. (25) separationstatisticscanbeestablishedinthesimulation. Thus, cross ≃ ∇ p all other things being equal, for the correct mass function to We show these two redshifts in Figure 5 for 10 different box exist in the box, one must run the simulation forward by an sizesrangingfrom1to512h- 1Mpcandfor2563and1283parti- amount sufficiently greater than the time taken to establish cles. Thetoppanelshowstheaverageredshiftofthefirstcross- the correct small-scale power on first-halo scales while eras- ing as a function of box size (which corresponds to the max- ingmemoryon these scales ofthe initialconditions. If thisis imum in Fig. 4). The bottom panel shows the redshift where not done, structure formationwill be suppressed, leading to a the first “gridcrossing”occurs(correspondingto the righttail loweringofthehalomassfunction. in Fig. 4). To estimate the scatter in the results, we havegen- Because thereisnofullysatisfactorywaytocalculatez in in eratedfivedifferentrealizationsforeachbox. Asexpected,the ordertocomputethemassfunctionatagivenaccuracy,wesub- small boxesshow much morescatter. The averageredshiftof jected every simulation box to convergence tests in the mass the first crossing in the 1h- 1Mpc box varies between z = 63 functionwhilevaryingz . Theresultsshowninthispaperare and83,whilethereisalmostnoscatterinthe512h- 1Mpcbox. in allconvergedtothesub-percentlevelinthemassfunction. We Since φ/∆ isindependentofredshiftintheZel’dovichap- p giveanexampleofonesuchconvergencetestbelow. |∇ | proximation,asimplescalingdeterminestheappropriateinitial redshiftfromtheseplots. Forexample,ifaparticleshouldnot 4.2.1. InitialRedshiftConvergenceStudy havemovedmorethan0.3∆ onaverageattheinitialredshift, p the average redshiftof first crossing has to be multiplied by a As mentionedabove,we havetested andvalidatedouresti- factor 1/0.3=3.¯3. For an 8h- 1Mpc box this leads to a mini- matesfor the initialredshiftfor all the boxesused in the sim- mum starting redshift of z=230, while for a 126h- 1Mpc box ulation suite via convergence studies. Here, we show results this suggests a starting redshift of zin =50. The 1283 particle foran8h- 1Mpcboxwithinitialredshiftszin=50,150,and250 curve can be scaled to the 2563 particle curve by multiplying in Figure 6, where the mass functionsat z=10 are displayed. by a factor of 2. Curves for differentparticle loadingscan be For the lowestinitial redshift, zin=50, the averageinitialpar- obtainedsimilarly. ticlemovementis1.87∆p,whilesomeparticlestravelasmuch as5.03∆ . Thisclearlyviolatestherequirementthattheinitial p particle grid distortion be kept sufficiently below 1 grid cell. 4.2. TransientsandMixing Thestartingredshiftz =150leadstoanaveragedisplacement in 10 TheHaloMassFunction FIG. 6.—Dependenceofthemassfunctionontheinitialredshift. Theresultsareatz=10fromthree8h- 1Mpcboxsimulationswithzin=50(left),zin=150 (middle),andatzin=250(right).Themassfunctionintheleftpanelissystematicallylowerthantheothertwobyroughly15%.Poissonerrorbarsareshown. of0.63∆ andamaximumdisplacementof1.71∆ ,andthere- Ω , including baryons and dark matter, the physical box size p p m forejustbarelyfulfillstherequirements. Forz =250wefind L ,andthenumberofsimulationparticlesn3: in box p anaveragedisplacementinthisparticularrealizationof0.37∆ andThaembaoxttiommumplodtisipnleaaccehmoefntthoeft1h.r0e0e∆papn.elsofFigure6showps mparticle=2.775×1011Ωm n hL-b1oMx pc 3h- 1M⊙. (26) (cid:18) p (cid:19) the ratio of the mass functions with respect to the Warren fit. Therequiredforceresolutiontoresolvethechosensmallest In the middle and right panels the ratio for the largest halo is haloscanbeestimatedverysimply.Supposeweaimtoresolve outsidethedisplayedrange. Themassfunctionfromthesimu- avirializedhalowithcomovingradiusr∆ atagivenredshiftz, lationstartedatzin=50(leftpanel)isnoticeablylower, 15%, where∆istheoverdensityparameterwithrespecttothecritical ∼ thanfortheothertwosimulations.Themassfunctionsfromthe two higherredshiftstartsareingoodagreement,showingthat densityρc. Thecomovingradiusr∆ isgivenby tighseencoecrnhasoleicrcvoeanftcoivlruesa,ivoaennrdailgltuhesagttrraoidtneeddicbsatyonFrtsiiagofunerleoyf6uaispsept(hr0oa.xt5ii–mf0aa.6tse)iml∆yup0.la.3tTi∆ohnep r∆=9.51×10- 5(cid:20)ΩΩ(mz)(cid:21)1/3(cid:18)∆1 hM- 1M∆c⊙(cid:19)1/3h- 1Mpc, (27) isstartedtoolate,halosarefoundtobemissingovertheentire where Ω(z)=Ωm(1+z)3/[Ωm(1+z)3+ΩΛ] and the halo mass massrange.Withthelatestart,thereislesstimetoformbound M∆c=mpartnh,wherenh isthenumberofparticlesinthehalo. objects. Also,someparticlesthatarestillstreamingtowardsa Wemeasuretheforceresolutionintermsof halodonothaveenoughtimetojoinit. Bothoftheseartifacts δ = Lbox. (28) leadtoanoveralldownshiftofthemassfunction. f n g Tosummarize,requiringalimitoninitialdisplacementssets In the case of a grid code, n is literally the number of grid the starting redshift much higher than simply demanding that g pointsperlineardimension;foranyothercode,n standsforthe all modes in the box stay linear. Indeed, the commonly used g number of “effective softening lengths” per linear dimension. lattercriterion(withδrms 0.1)isnotadequateforcomputing the halo massfunctionat∼highredshifts. Onemustverifythat To resolve halos of mass M∆c, a minimal requirement is that the code resolution be smaller than the radius of the halo we the chosen z sets an early enough start as shown here. We in wishtoresolve: commentonpreviousresultsfromothergroupswithrespectto thisfindingbelowin§6. δf<r∆. (29) Note that this minimal resolution requirement is aimed only 4.3. ForceandMassResolution at capturing halos of a certain mass, not at resolving their in- terior profile. Next, inserting the expression for the particle Wenowtakeupaninvestigationofthemassandforcereso- mass (eq. [26]) and the comoving radius (eq. [27]) into the lutionrequirements.Thefirstusefulpieceofinformationisthe requirement(eq [29]) and employingthe relation between the size ofthesimulationbox: fromFigure2wecaneasilytrans- interparticlespacing∆ andtheboxsize∆ =L /n ,theres- late the numberdensityintowhen thefirst halois expectedto p p box p olutionrequirementreads appearinaboxofvolumeV. Forexample,ahorizontallineat n=10- 6wouldtellusatwhatredshiftwewouldexpectonaver- δ n Ω(z) 1/3 agetofind1haloofacertainmassina(100h- 1Mpc)3box.The ∆f <0.62 h∆ . (30) first halo of mass 1011- 1012h- 1M⊙ will appear at z 15.5, p (cid:20) (cid:21) and the first cluster-like object of mass 1014- 1015h-≃1M at We now illustrate the use of this simple relation with an ex- ⊙ z 2. Of course, these statements only hold if the mass and ample. Let ∆=200 and consider a ΛCDM cosmology with fo≃rceresolutionaresufficienttoresolvethesehalos. Themass Ωm =0.3. Then for PM codes for which δf/∆p =np/ng, we ofaparticleinasimulation,andhencethehalomass,isdeter- havethe followingconclusions. If the numberof meshpoints minedbythreeparameters: themattercontentoftheUniverse isthesameasthenumberofparticles(np=ng),haloswithless