SUBMITTEDTOAPJDECEMBER11,2012 SLAC-PUB-15331 PreprinttypesetusingLATEXstyleemulateapjv.11/10/09 THEROCKSTARPHASE-SPACETEMPORALHALOFINDERANDTHEVELOCITYOFFSETSOFCLUSTERCORES PETERS.BEHROOZI,RISAH.WECHSLER,HAO-YIWU KavliInstituteforParticleAstrophysicsandCosmology,DepartmentofPhysics,StanfordUniversity; DepartmentofParticlePhysicsandAstrophysics,SLACNationalAcceleratorLaboratory;Stanford,CA94305 SubmittedtoApJDecember11,2012 ABSTRACT Wepresentanewalgorithmforidentifyingdarkmatterhalos,substructure,andtidalfeatures.Theapproach isbasedonadaptivehierarchicalrefinementoffriends-of-friendsgroupsinsixphase-spacedimensionsandone timedimension,whichallowsforrobust(grid-independent,shape-independent,andnoise-resilient)trackingof substructure; as such, it is named ROCKSTAR (RobustOverdensityCalculation using K-Space Topologically Adaptive Refinement). Our method is massively parallel (up to 105 CPUs) and runs on the largest current simulations (>1010 particles) with high efficiency (10 CPU hours and 60 gigabytes of memory required per billion particles analyzed). A previous paper (Knebe et al. 2011) has shown ROCKSTAR to have excellent recoveryofhaloproperties;weexpandonthesecomparisonswithmoretestsandhigher-resolutionsimulations. Weshowasignificantimprovementinsubstructurerecoverycomparedtoseveralotherhalofindersanddiscuss thetheoreticalandpracticallimitsofsimulationsinthisregard. Finally,wepresentresultswhichdemonstrate conclusivelythatdarkmatterhalocoresarenotatrestrelativetothehalobulkorsubstructureaveragevelocities andhavecoherentvelocityoffsetsacrossawiderangeofhalomassesandredshifts.Formassiveclusters,these offsets can be up to 350 km s- 1 at z=0 and even higher at high redshifts. Our implementation is publicly availableathttp://code.google.com/p/rockstar. Subjectheadings:darkmatter—galaxies:abundances—galaxies:evolution—methods:N-bodysimulations 1. INTRODUCTION cretionhistories,seee.g.Wechsleretal.2002;Springeletal. In the paradigm of Lambda Cold Dark Matter (ΛCDM), 2005; Faltenbacher et al. 2005; Allgood et al. 2006; Harker etal.2006;Tweedetal.2009;Wetzeletal.2009). the majority of the matter density of the universe does not This level of uncertainty is unacceptable for current and couple to electromagnetic fields, leaving it detectable only future surveys, including data expected to come from the throughitsgravitationalandpossiblyweakforceinteractions. Baryon Oscillation SpectroscopicSurvey(BOSS), Dark En- Nonetheless,theeffectsofdarkmatteronthevisibleuniverse ergy Survey (DES), BigBOSS, Panoramic Survey Tele- arespectacular,asthesteepgravitationalpotentialsofbound scope and Rapid Response System (Pan-STARRS), Ex- darkmatterhaloschannelbaryonstogether,formingthebirth- tended Roentgen Survey with an Imaging Telescope Ar- placesofvisiblegalaxies. Inthismodel,thelocations,sizes, ray (eROSITA), Herschel, Planck, James Webb Space Tele- and merging histories for galaxies are thus intricately con- scope(JWST),andLargeSynopticSurveyTelescope(LSST) nectedtothegrowthofbounddarkmatterstructures. (Schlegel et al. 2009; The Dark Energy Survey Collabora- Testing this model in detail requires extensive computer tion 2005; Schlegel et al. 2011; Chambers 2006; Predehl simulations,asthecomplicatednonlinearevolutionofstruc- et al. 2010; Pilbratt et al. 2010; The Planck Collaboration ture growth cannotbe fully evaluated by hand. The simula- 2006;Gardneretal.2006;LSSTScienceCollaborationsetal. tions generally follow the evolution of a set of dark matter 2009),inordertofullyrealizetheconstrainingpowerofthese particlesand outputthe positionsand velocitiesof the parti- observations.Indeed,derivedquantitiessuchasthehalomass cles at several discrete timesteps. These outputs must then functionandautocorrelationfunctionsmustbeunderstoodat bepostprocessedtodeterminethelocationsandpropertiesof the one percentlevel in order for theoreticaluncertaintiesto bounddarkmatterstructuresalsoknownas“halos”—namely, beatthesamelevelasstatisticaluncertaintiesinconstraining, thelocationsandpropertieswhichinfluencetheformationof e.g.,darkenergy(Wuetal.2010;Cunha&Evrard2010,for visiblegalaxies. Thispostprocessing(“halo-finding”)neces- example).Althoughtheremaybealimittotheaccuracypos- sarily involves both ambiguity and imprecision—ambiguity siblewithdarkmattersimulationsalone,giventheimpactof in the definitions (e.g., the center of a bound halo) and im- baryonsondarkmatterhaloprofiles(Staneketal.2009)and precisionindetermininghalopropertiesduetolimitedinfor- the power spectrum (Rudd et al. 2008), different halo find- mation (e.g., for halos consisting of only a few dark matter ers running on the same simulation demonstrate that a large particles,orfordeterminingparticlemembershipintwoover- fractionof this uncertaintyis still due to the processof halo lappinghalos). Currently,as a consequence,essential statis- findingitself(Knebeetal.2011). ticalmeasures(e.g.,thenumberdensityofhaloswithagiven As previously mentioned, some of these uncertainties are mass)areknowntoatbest5-10%evenforaspecificcosmol- duetolimiteduseofinformation. Considerableprogresshas ogy(Tinkeretal.2008). Moreover,halopropertiesarerarely been made since the first generation of position-space halo checked for consistency over multiple timesteps, which can finders(Davis et al. 1985;Lacey & Cole 1994), which used beaseriousproblemforrobustmodelingofgalaxyformation only particle locations to determine bound halos. Currently, theories(seehoweverBehroozietal.2011;thishasalsobeen themostadvancedalgorithmsareadaptivephase-spacefind- addressedinspecificcontextswhencreatingmergertrees,for ers (e.g., Maciejewski et al. 2009; Ascasibar et al. 2010), exampleimprovingsubhalo tracking and smoothermass ac- whichmakeuseofthefullsixdimensionsofparticlepositions Work supported in part by US Department of Energy contract DE-AC02-76SF00515. 2 BEHROOZI,WECHSLER,&WU and velocities. At the same time, an often-overlooked(with 2.1. SummaryofPreviousApproaches the possible exception of Tweed et al. 2009) aspect of halo Previously published approaches to halo finding may be findingistheextrainformationstoredinthetemporalevolu- classified,withfewexceptions,intotwolargegroups. Spher- tionofboundparticles. Inthis paper,we detailan advanced ical overdensity (SO) finders, such as Amiga’s Halo Finder adaptive phase-space and temporal finder which is designed (AHF;Knollmann&Knebe2009),AdaptiveSphericalOver- tomaximizeconsistencyofhalopropertiesacrosstimesteps, density Halo Finder (ASOHF; Planelles & Quilis 2010), ratherthanjustwithinasinglesimulationtimestep. Together Bound Density Maxima (BDM; Klypin et al. 1999), SO withacompanionpaper(Behroozietal.2011),whichdetails (Lacey&Cole1994;Jenkinsetal.2001;Evrardetal.2002), theprocessofcomparingandvalidatinghalocatalogsacross parallel SO (pSO; Sutter & Ricker 2010), Voronoi Bound timestepstocreategravitationallyconsistentmergertrees,our Zones(VOBOZ;Neyrincketal.2005),andSplineKernelIn- combinedapproachis the first to use particle informationin terpolativeDenmax(SKID;Stadel2001)proceedbyidentify- seven dimensions to determine halo catalogs, allowing un- ingdensitypeaksintheparticledistributionandthenadding precedentedaccuracyandprecisionindetermininghaloprop- particles in spheres of increasing size around the peaks un- erties. til the enclosed mass falls below a predetermined density Furthermore, in contrast to previous grid-based adaptive threshold (a top-down approach). Friends-of-friends (FOF) phase-spaceapproaches,oursisthefirstgrid-independentand andHOP-based(Eisenstein&Hut1998)halofinders,suchas orientation-independentapproach;itisalsothefirstpublicly- FOF,SUBFIND,theLANLhalofinder,parallelFOF(pFOF), availableadaptivephase-spacecodedesignedtorunonmas- Ntropy-fofsv, and AdaptaHOP (Davis et al. 1985; Springel sivelyparallelsystemsontheverylargesimulations(>1010 et al. 2001; Habib et al. 2009; Rasera et al. 2010; Gard- particles)whicharenecessarytoconstrainstructureformation neretal.2007;Tweedetal. 2009),collectparticlestogether acrosstherangeofscalesprobedbycurrentandfuturegalaxy which fall above a certain density threshold and then, if so observations. Finally, we remark that our halo finder is the designed, search for substructure inside these particle col- first which is successfully able to probesubstructuremasses lections (a bottom-up approach). Phase-space finders typi- downto theverycentersofhosthalos(seealso Knebeetal. cally extend these two approachesto include particle veloc- 2011),whichisessentialforafullmodelingofgalaxystatis- ity information, either by calculating a phase-space density, ticsandwillenablefuturestudiesoftheexpectedbreakdown such as the Hierarchical Structure Finder (HSF; Maciejew- betweenhalopositionsandgalaxypositionsduetotheeffects ski et al. 2009) and the Six-DimensionalHierarchical Over- ofbaryoninteractionsattheverycentersofgalaxies. density Tree (HOT6D; Ascasibar et al. 2010) or by using a Throughout,wehavepaidcarefulattentionnotonlytothe phase-spacelinkinglength,asdoesSix-DimensionalFriends- basictaskofassigningparticlestohalos,butalsotothepro- of-Friends(6DFOF;Diemandetal.2006). cess of estimating useful properties from them to compare There are three notable exceptions to these algorithms, with observations. While in many cases galaxy surveys are namelytheORIGAMIhalofinder(discussedinKnebeetal. not able to probe halo properties to the same precision as 2011),theHierarchicalBound-Tracingalgorithm(HBT;Han halo finders in simulations, one significant counterexample etal.2011),andSURV(Giocolietal.2010). ORIGAMIop- exists. It is a common practice especially for halo finders erates by examiningphase-space shell crossingsfor the cur- based on the friends-of-friends algorithm (e.g., Davis et al. rent particle distribution as compared to the initial particle 1985; Springel et al. 2001; Habib et al. 2009; Rasera et al. distribution; shells which have crossed along three dimen- 2010; Gardner et al. 2007; see also §2.1) to calculate halo sionsare consideredto be halos(as opposedto shells which velocitiesbyaveragingallhaloparticlevelocitiestogetherto have crossed along one or two dimensions, which would be find a bulk velocity. Examinationof the differencebetween consideredas sheets and filaments, respectively). HBT uses velocitiesintheinnerregionsofhalosandthebulkaveraged a friends-of-friendsapproach to find distinct halos and uses velocity suggests that the bulk average velocity may be off- particle lists fromdistinct halosat previoustimesteps to test set by several hundred km s- 1 from the velocity at the ra- for the presence of subhalos. SURV is a very similar algo- diuswheregalaxiesareexpectedtoreside;differencesatthis rithm,exceptthatdistincthalosareidentifiedusingspherical scaleareeasilydetectableinclusterredshiftsurveysandmay overdensities. These algorithmsallrely heavilyon temporal alsofactorininterpretingobservationsofthekineticSunyaev- informationintheirapproachtohalofinding. Zel’dovicheffect. Asthisdifferencehasanimportantimpact onthe usefulnessof derivedhalo properties,weadditionally 2.2. LimitationsofPreviousAlgorithms performan investigationof the core-bulkvelocitydifference inhalosacrossawiderangeofredshiftsandmasses. In order to develop an improved halo finder, it is impor- Webeginthispaperwithasurveyofpreviousworkinhalo tantto understandsome ofthe shortcomingsofpreviousap- findingaswellaspreviouslimitationsin§2. Wediscussour proaches. Forthevastmajorityofhalos,eventhemostbasic improvedmethodologyin§3andconductdetailedtestsofour of algorithms (FOF and SO) do an acceptable job of deter- approachin§4. Wepresentananalysisofthetheoreticaland mininghalopropertiesto10%accuracy(Knebeetal.2011). practicallimitationsof simulationsin termsof trackingsub- However, recent interest in the detailed properties and his- structure in §5. Finally, our results concerning the velocity tories of halos — e.g., precision mass functions and merger offsets of cluster cores are presented in §6; we summarize trees and the shape of tidal structures — requires improve- our conclusionsin §7. Multiple simulationsincludingslight mentstoolderapproaches;thishasresultedinaproliferation variationsof cosmologicalparametersare considered in this ofnewcodesinthepastdecade(summarizedinKnebeetal. paper; all simulations model a flat ΛCDM universe and we 2011). always take the HubbleconstantH to be 70 km s- 1 Mpc- 1; The most significant improvements to halo finding have 0 equivalently,h=0.7. come from using the informationfrom six-dimensional(po- sitionandvelocity)phasespace. Twotraditionalweakpoints 2. PREVIOUSHALOFINDERS for 3D (position-space) halo finders have been major merg- TheROCKSTARHaloFinder 3 ers and subhalos close to the centers of their host halos. In halo concentrationsat a given halo mass. At the same time, both cases, the density contrast alone is not enough to dis- whileTweedetal.(2009)successfullyresolvesthisproblem, tinguishwhichparticlesbelongtowhichhalo: whentwoha- itnonethelessonlyfindshalosinpositionspace,andthushas losarecloseenough,theassignmentofparticlestohalosbe- thesameweaknessesinidentifyingsubhaloandmajormerger comesessentiallyrandomintheoverlapregion. However,as properties. long as the two halos have relative motion, six-dimensional TheHBT and SURV algorithms(Han et al. 2011;Giocoli halo finders can use particle velocity space information to etal.2010)usetheingeniousapproachoftracingsubhalosby very effectively determine particle-halo membership. This, using the particles found in previously distinct halos, which coupledwiththe abilityof6Dhalofindersto findtidalrem- couldpotentiallyalsosolvemanyoftheseproblems.Yet,they nants(whicharecondensedinphasespacebutnotinposition bothalsoincludeassumptionsaboutaccretionontosubhalos space), means that phase-space capabilities are required for (e.g., that subhalos cannot accrete background matter from themostaccurateandinterestingstudiesofdarkmatterhalos. the host) which are untrue in HBT’s case for large linking- Atthesametime,phasespacepresentsauniquechallenge. lengths (as halos will be identified as satellites far from the While position space and velocity space have well-defined actualvirialradiusofthehost)andforbothhalofinderswith distance metrics, there is not a single, unique way to com- majormergers(wheresatelliteandhostaremoreambiguous; bine position and velocity distance into a single metric. For they can in any case easily trade particles with each other). a useful definition of phase-space distance, one needs to be Theseassumptionsvastlysimplifythecodeatsomeexpense abletodecide,e.g.,whetheranobjectjustpassingbytheori- tothecompletenessandaccuracyofthemassfunction.More gin at 1 km s- 1 is closer or farther than an object at rest 1 seriously, the design of the algorithms requires temporal in- kpcaway. Oneapproach,usedby6DFOF,istochooseinad- formation to find subhalos; in cases where simulation out- vanceastaticconversionbetweenvelocityandpositionspace. puts are spaced very far apart or when only one timestep is Whilesimple,thisapproachseemssomewhatself-limiting:if available,theycannoteffectivelyfindsubstructure. Theseis- too short a linking-length is chosen, the full extent of sub- suesareinprinciplefixable:futureversionsofthealgorithms structures cannot be found; if too large a linking-length is could easily combine advanced single-timestep substructure chosen, then otherwisedistinct substructureswillbe merged findingwithchecksagainstprevioustimesteps’particlelists. together. A demonstrablysuperiorapproach,atleast in termsof re- 3. ANIMPROVEDAPPROACH:ROCKSTAR coveringhaloproperties(Knebeetal.2011),istoadaptively Ourprimarymotivationindevelopingahalofinderwasto define a phase-space metric. Both HSF and HOT6D subdi- improve the accuracy of halo merger trees that are required vide the simulation space into six-dimensional hyperboxes foran understandingofgalaxyevolution. Thedesignof our containing (at the maximum refinement level) as little as a halo finder was thus motivated by a requirementfor consis- single particle each. For a given particle, the enclosing hy- tent accuracy across multiple timesteps. This interest led to perboxgivesalocalestimateofthephase-spacemetric,based the developmentof a unique, adaptivephase-spacetemporal ontherelativesizesofthehyperbox’sdimensionsinposition algorithmwhich is provablyindependentof halo orientation and velocity space. The usefulness of this estimate depends and velocity relative to the simulation axes, and which also heavilyonthemethodforpartitioningspaceintohyperboxes; attempts to be highly preserving of particle-halo and halo- HSF uses, for example, an entropy-basedapproachto deter- subhalo membership across timesteps. In addition, we paid mine whether more information is contained in the particle special attention to the algorithm’sefficiency and paralleliz- locationsforpositionspaceorvelocityspace. ability,toallowittorunonthelargestcurrentdatasetsandso These algorithms all give excellent results for identifying thatitcouldeasilyscaletothenextgenerationofsimulations. halo centers at a single timestep. However, consistent halo Thusfar,wehaverunthehalofinder(andinmanycasesthe catalogsacrosstimestepsareoftencompromisedbyafunda- partner merger tree code) on the Bolshoi ( 1010 particles, mentalambiguityin the definitionof a hosthalo. Formajor ∼ Klypinetal.2011)andLasDamassimulations(200boxesof mergers, it is often unclear which halo center represents the larger “host” or central halo, and which represents the sub- 1- 4 109particleseach, McBrideetal. inpreparation),on halo. Phase-spacehalofindinghelpswhenthetwohalocen- sever×al20483 simulationsrun to create Dark EnergySurvey tersarerelativelyfarapart(i.e.,weaklyinteracting),because (DES)simulatedskycatalogs,on onehundredhighresolu- ∼ thereexistsastrongcorrelationbetweenthevelocitiesofpar- tionhalossimulatedaspartoftheRHAPSODY project(Wu ticles in the halo outskirtsand halo centers. However,when etal.2012a,b),andonhalosA-1throughA-5oftheAquarius the centers come close enough to interact strongly, this cor- Simulation(upto1.4 109particlesinasinglehalo;Springel × relation is weakened, and it becomesmuchmore difficultto etal.2008). accuratelyassignparticlestothehalos.Asaresult,itismuch As a first step, our algorithm performs a rapid variant of more difficultto determinewhich of the halo centers should the3Dfriends-of-friends(FOF)methodtofindoverdensere- beconsideredthehosthalo.Sincethedefinitionofhalomass gions which are then distributed among processorsfor anal- often includes the mass of subhalos, this problemcan result ysis (§3.1). Then, it builds a hierarchy of FOF subgroups inlargemassfluctuationsacrosstimestepsformerginghalos. in phase space by progressivelyand adaptively reducingthe A number of solutions to this problem have been pro- 6D linking length, so that a tunable fraction of particles are posedandexaminedwiththeAdaptaHOPhalofinder(Tweed capturedateachsubgroupascomparedtotheimmediatepar- et al. 2009). Tweed et al. (2009) found that a temporal ap- ent group (§3.2). Next, it converts this hierarchy of FOF proach(examiningthehostvs. subhaloassignmentatearlier groups into a list of particle memberships for halos (§3.3). timesteps)wasmostsuccessfulatfixingthisproblem. Other It then computes the host halo/subhalo relationships among methods, such as choosing the densest halo center to be the halos, usinginformationfromtheprevioustimestepif avail- hosthalo,haveinherentinstabilitiesbecauseofthespreadin able(§3.4). Finally,itremovesunboundparticlesfromhalos andcalculates halo properties, beforeautomaticallygenerat- 4 BEHROOZI,WECHSLER,&WU 1. The simulationvolumeisdivided Asthe 3D FOF groupsare used in ourmethodonly to di- into 3D Friends-of-Friends groups vide up the simulation volume into manageable work units, weinsteadmakeuseofamodifiedalgorithmwhichisanor- for easy parallelization. derofmagnitudefaster. Asisusual,neighboringparticlesare assignedtobeinthesamegroupiftheirdistanceiswithinthe linkinglength. However,ifaparticlehasmorethanacertain 2. For each group, particlepositions numberofneighborswithinthelinkinglength(16,inourver- and velocities are divided (normal- sion),thentheneighbor-findingprocessforthoseneighboring ized) by the group position and ve- particlesisskipped. Instead,neighborsfortheoriginalparti- cle outto twice the originallinkinglength are calculated. If locity dispersions, giving a natural anyofthoseparticlesbelongtoanotherFOFgroup,thatcor- phase-space metric. respondinggroup is joined with that of the originalparticle. Thus,althoughtheneighbor-findingprocesshasbeenskipped 3. A phase-space linking length is forthenearestparticles,groupswhichwouldhavebeenlinked throughthoseintermediateneighborsarestilljoinedtogether. adaptively chosen such that 70% of Thisprocessthereforelinkstogetheratminimumthesame the group’s particles are linked to- particles as in the standard FOF algorithm—finefor our de- getherin subgroups. sired purpose—but does so much faster: neighbors must be calculated over a larger distance, but many fewer of those 4. The process repeats for each calculations must be performed. Indeed, rather than slow- ing down as the linking length is increased, our variation of subgroup: renormalization, a new FOF becomesfaster. Because of this, we are free to choose linking-length, and a new level of anexceptionallylargevalueforthelinkinglength.Avalueof substructurecalculated. b=0.2istoosmall,asitdoesnotincludeallparticlesoutto thevirialradius(Moreetal.2011);afterevaluatingdifferent 5.Oncealllevelsofsubstructureare choicesforb(see§4.5),wechoseb=0.28,whichguarantees thatvirialsphericaloverdensitiescanbedeterminedforeven found, seed halos are placed at the themostellipsoidalhalos. lowest substructure levels and par- The most important parallelization work occurs at this ticles are assigned hierarchically to stage.Separatereadertasksloadparticlesfromsnapshotfiles. Depending on the number of available CPUs for analysis, a theclosestseedhaloinphasespace. masterprocessdividesthesimulationregionintorectangular boundaries, and it directs the reader tasks to send particles 6. Once particles have been as- withinthoseboundariestotheappropriateanalysistasks. signed to halos, unbound particles Eachanalysistask firstcalculates3DFOFsin itsassigned are removed and halo properties analysis region, and FOFs which span processor boundaries areautomaticallystitchedtogether. TheFOFgroupsarethen (positions, velocities, etc.) are distributed for further phase-space analysis according to in- calculated. dividualprocessorload. Theload-balancingprocedureisde- scribedinfurtherdetailinAppendixA. Currently,single3D FIG.1.—Avisualsummaryoftheparticle-haloassignmentalgorithm. FOF groups are analyzed by at most one processor. Also, in the currentimplementation, there is no supportfor multi- ple particle masses, although support could easily be added ingparticle-basedmergertrees(§3.5). A visualsummaryof by varying the linking length depending on particle mass. thesestepsisshowninFig.1. Providedenoughinterest,supportformultiprocessoranalysis of single large halos as well as support for multiple particle 3.1. Efficient,ParallelFOFGroupCalculation massesmaybeaddedinafutureversionofROCKSTAR. The 3D friends-of-friends (FOF) algorithm has existed since at least 1985 (Davis et al. 1985). In principle, imple- 3.2. ThePhase-SpaceFOFHierarchy mentation of the algorithm is straightforward: one attaches For each 3D FOF group which is created in the previous two particlesto the samegroupif theyarewithin a prespec- step, the algorithmproceedsbybuildinga hierarchyof FOF ified linking length. Typically, this linking length is chosen subgroupsinphasespace. Deeperlevelsofsubgroupshavea intermsofafractionbofthemeaninterparticledistance(in tighter linking-lengthcriterion in phase space, which means ourcode,asforothers,thecuberootofthemeanparticlevol- that deeper levels correspond to increasingly tighter isoden- ume);commonvaluesforgeneratinghalocatalogsrangefrom sity contoursaroundpeaksin the phase-spacedensity distri- b=0.15tob=0.2(Moreetal.2011). bution. Thisenablesaneasywaytodistinguishseparatesub- In practice, this meansthat one mustdetermine neighbors structures—abovesomethresholdphase-spacedensity,their foreveryparticlewithinasphereofradiusequaltothelink- particle distributions must be distinct in phase space; other- ing length. Even with an efficient tree code (we use a cus- wise,itwouldbedifficulttojustifytheseparationintodiffer- tombinaryspacepartitioningtree),thisrepresentsagreatdeal entstructures.1 of wasted computation,especially in dense cluster cores. In Beginning with a base FOF group, ROCKSTAR adaptively such cases, particles might have tens of thousandsof neigh- borswithinalinkinglength,allofwhichwilleventuallyend 1Thiswouldnotbetrueif,e.g.,haloshadPlummerprofilesorotherwise upinthesameFOFgroup. flatdensityprofilesintheircenters. TheROCKSTARHaloFinder 5 chooses a phase-space linking length based on the standard in the autocorrelation function close to the simulation force deviationsoftheparticledistributioninpositionandvelocity resolution. space. Thatis,fortwoparticles p and p inthebasegroup, Foraparentgroupwhichcontainsonlyasingleseedhalo, 1 2 thephase-spacedistancemetricisdefinedas: all the particles in the groupare assigned to the single seed. Foraparentgroupwhichcontainsmultipleseedhalos,how- ~x - ~x 2 ~v - ~v 2 1/2 ever, particles in the group are assigned to the closest seed d(p1,p2)= | 1σ22| +| 1σ22| , (1) halo in phase space. In this case, the phase-space metric is (cid:18) x v (cid:19) setbytheseedhaloproperties,sothatthedistancebetweena whereσ andσ aretheparticlepositionandvelocitydisper- halohandaparticle pisdefinedas: x v sionsforthegivenFOFgroup;thisisidenticaltothemetricof 1/2 Gottlöber(1998). Foreachparticle,thedistancetothenear- ~x - ~x 2 ~v - ~v 2 est neighbor is computed; the phase-space linking length is d(h,p)= |rh2 p| +| hσ2p| (3) thenchosensuchthataconstantfraction f oftheparticlesare dyn,vir v ! v linkedtogetherwithatleastoneotherparticle.Inlargegroups r =v t = max (4) dyn,vir maxdyn,vir (>10,000particles),wherecomputingthenearestneighborfor 4πGρ allparticlescanbeverycostly,thenearestneighborsareonly 3 vir calculated for a random10,000-particlesubset of the group, q whereσ istheseedhalo’scurrentvelocitydispersion,v is v max asthisissufficienttodeterminethelinkinglengthtoreason- its current maximum circular velocity (see §3.5), and “vir” ableprecision. specifies the virial overdensity, using the definition of ρ vir The proper choice of f is constrained by two considera- from Bryan & Norman (1998), which corresponds to 360 tions.Ifonechoosestoolargeavalue(f >0.9),thealgorithm timesthebackgrounddensityatz=0,however,otherchoices willtake muchlonger,anditcanalso findspurious(notsta- ofthisdensitycaneasilybeapplied. tistically significant)subgroups. Ifonechoosestoolow ofa Usingtheradiusr astheposition-spacedistancenor- value(f <0.5),thealgorithmmaynotfindsmallersubstruc- dyn,vir malization may seem unusual at first, but the natural alter- tures. Assuch, we usean intermediatevalue(f =0.7);with native (using σ ) gives unstable and nonintuitive results. At x therecommendedminimumthresholdforhaloparticles(20). fixedphase-spacedensity,subhalosandtidalstreams(which In our tests of the mock NFW (Navarro et al. 1997) halos havelowervelocitydispersionsthanthehosthalo)willhave described in Knebe et al. (2011), this results in a false posi- largerposition-spacedispersionsthanthehosthalo. Thus,if tiverateof10%for20-particlegroupscomparedtoacosmo- σ wereused,particlesintheoutskirtsofahalocouldbeeas- x logicalsubhalodistribution,whichdeclinesexponentiallyfor ilymis-assignedtoasubhaloinsteadofthehosthalo. Using larger group sizes. These false positives are easily removed r , on the otherhand, preventsthis problemby ensuring by the significance and boundedness tests described in §3.3 den,vir thatparticlesassignedtosubhaloscannotbetoofarfromthe and§3.5.3. maindensitypeakeveniftheyarecloseinvelocityspace.2 In- Once subgroups have been found in the base FOF group, tuitively,the largesteffectof usingr isthat velocity-space vir thisprocessisrepeated. Foreachsubgroup,thephase-space informationbecomesthedominantmethodofdistinguishing metric is recalculated, and a new linking-length is selected particlemembershipwhentwo halosare within each other’s such that a fraction f of the subgroup’s particles are linked virialradii.3 togetherinto sub-subgroups. Groupfindingproceedshierar- This process of particle assignment assures that substruc- chicallyinphasespaceuntilapredeterminedminimumnum- ture masses are calculated correctly independently of the ber of particles remain at the deepest level of the hierarchy. choiceof f,thefractionofparticlespresentineachsubgroup Here we set this minimum number to 10 particles, although relativetoitsparentgroup.Inaddition,forasubhalocloseto halopropertiesarenotrobustapproachingthisminimum. thecenterofitshosthalo,itassuresthathostparticlesarenot mis-assignedtothesubhalo—thecentralparticlesofthehost 3.3. ConvertingFOFSubgroupstoHalos willnaturallybecloserinphasespacetothetruehostcenter ForeachofthesubgroupsatthedeepestleveloftheFOFhi- thantheyaretothesubhalo’scenter. erarchy(correspondingtothelocalphase-spacedensitymax- ima),aseedhaloisgenerated.Thealgorithmthenrecursively 3.4. CalculatingSubstructureMembership analyzes higher levels of the hierarchyto assign particles to In addition to calculating particle-halo membership, it is these seed halosuntilall particlesin the originalFOFgroup alsonecessarytodeterminewhichhalosaresubstructuresof havebeenassigned.Topreventcaseswherenoisegivesriseto other halos. The most common definition of substructure is duplicatedseedhalos,weautomaticallycalculatethePoisson a bound halo contained within another, larger halo. Yet, as uncertainty in seed halo positions and velocities, and merge halo masses are commonly defined to include substructure, thetwoseedhalosiftheirpositionsandvelocitiesarewithin the question of which of two halos is the largest (and thus, 10σ of the uncertainties. Specifically, the uncertainties are whichshouldbecalledasatelliteoftheother)canchangede- calculated as µx =σx/√n and µv =σv/√n, where σx and σv pendingonwhichsubstructureshavebeenassignedtothem. arethepositionandvelocitydispersions,andnisthenumber Thisisoneofthelargestsourcesofambiguitybetweenspher- ofparticles,allforthesmallerofthetwoseedhalos. Thetwo ical overdensity halo finders, even those which limit them- halosaremergedif selvestodistincthalos. (x - x )2µ- 2+(v - v )2µ- 2<10√2. (2) 2Fordeterminationoftidalstreams,this“problem”becomesa“feature,” 1 2 x 1 2 v anduseofσxmaybepreferabletorvir. Inourtests,qthisthresholdyieldsanear-featurelesshaloauto- 3 An alternate radius (e.g., r200b or r500c) could be used instead, but it wouldhaveaneffectonlyonasmallfractionofparticlesinasmallfraction correlationfunction;lowervaluesresultinaspuriousupturn ofhalos(majormergers). 6 BEHROOZI,WECHSLER,&WU Webreaktheself-circularitybyassigningsatellitemember- notedin§6,thehalocorecanhaveasubstantialvelocityoff- ship based on phase-space distances before calculating halo set from the halo bulk. Since the galaxy hosted by the halo masses. Treatingeachhalocenterlikea particle, we use the will presumably best track the halo core, we calculate the same metric as Eq. 3 and calculate the distance to all other mainvelocityforthehalousingtheaverageparticlevelocity haloswithlargernumbersofassignedparticles. Thesatellite withintheinnermost10%ofthehaloradius. Forcalculating haloinquestionisthenassignedtobea subhalooftheclos- thebound/unboundmassofthehalo(see3.5.2),however,we estlargerhalowithinthesame3Dfriends-of-friendsgroup,if usethemoreappropriateaveragedhalobulkvelocityinclud- oneexists. If the halo catalogatan earlier timestep is avail- ingsubstructure. able, this assignment is modified to include temporal infor- mation.Halocoresatthecurrenttimestepareassociatedwith 3.5.2. HaloMasses halosattheprevioustimestepbyfindingthehaloattheprevi- Forhalomasses, ROCKSTAR calculatessphericaloverden- oustimestepwiththelargestcontributiontothecurrenthalo sitiesaccordingtomultipleuser-specifieddensitythresholds: core’sparticlemembership. Then,host-subhalorelationships e.g.,thevirialthreshold,fromBryan&Norman(1998),ora are checked against the previous timestep; if necessary, the densitythresholdrelativetothebackgroundortothecritical choiceofwhichhaloisthehostmaybeswitchedsoastopre- density. Asisusual, these overdensitiesarecalculatedusing servethehost-subhalorelationshipoftheprevioustimestep. alltheparticlesforallthesubstructurecontainedinahalo.On As explained above, these host-subhalo relationships are theotherhand,subhalomasseshavetraditionallybeenama- onlyusedinternallyforcalculatingmasses:particlesassigned jorpointofambiguity(fordensity-spacehalofinders). With tothehostarenotcountedwithinthemassofthesubhalo,but a phase-space halo finder, such as ROCKSTAR, the particles particles within the subhalo are counted as part of the mass belongingto the subhalo can be more reliably isolated from ofthe hosthalo. Thischoiceassuresthatthe massof a halo thehost,andthusless ambiguityexists: thesamemethodof won’tsuddenlychangeasitcrossesthevirialradiusofalarger calculatingsphericaloverdensitiesmaybeappliedtojustthe halo, and it provides more stable mass definitions in major particlesbelongingtothesubhalo.Intermsofthedefinitionof mergers.Oncehalomasseshavebeencalculated,thesubhalo wherethesubhalo“ends,”Eq.3impliesthatthesubhaloedge membershipisrecalulatedaccordingtothestandarddefinition is effectively where the distribution of its particles in phase (subhalosarewithinr∆ofmoremassivehosthalos)whenthe spacebecomesequidistantfromthesubhaloanditshosthalo. mergertreesareconstructed. Ifalternatemassdefinitionsarenecessary,thehalofindercan Forclarity,itshouldbenotedthateverydensitypeakwithin output the full phase-space particle-halo assignments; these the original FOF analysis group will correspond to either a maythenbepost-processedbytheusertoobtainthedesired hosthaloorasubhalointhefinalcatalog.Ithasbeenobserved masses. thatFOFgroupswill“bridge”or“premerge”longbeforetheir correspondingSOhalocounterparts(e.g.,Klypinetal.2011). 3.5.3. UnbindingParticles However, as we calculate full SO propertiesassociated with By default, ROCKSTAR performs an unbinding procedure eachdensitypeak,asingleFOFgroupisnaturallyallowedto beforecalculating halo mass and v , althoughthis may be max containmultipleSOhosthalos;thusthebridgingorpremerg- switched off for studies of e.g., tidal remnants. Because the ingofFOFgroupsdoesnotaffectthefinalhalocatalogs. algorithm operates in phase space, the vast majority of halo particles assigned to central halos are actually bound. We 3.5. CalculatingHaloPropertiesandMergerTrees findtypicalboundednessvaluesof98%atz=0;see§4.5and Typically, several properties of interest are generated for Behroozietal.(2012a).Evenforsubstructure,unboundparti- halocatalogs.Regardlessofthequalityofparticleassignment clestypicallycorrespondtotidalstreamsattheoutskirtsofthe in the halo finder, carefulattention to halo propertycalcula- subhalo,makingacomplicatedunbindingalgorithmunneces- tionisessentialforconsistent,unbiasedresults. sary. For this reason, as well as to improve consistency of halomassesacrosstimesteps, we performonlya single-pass 3.5.1. HaloPositionsandVelocities unbindingprocedureusingamodifiedBarnes-Hutmethodto accurately calculate particle potentials (see Appendix B for For positions, Knebe et al. (2011) demonstrated that halo details).4 Sincemanyuserswillbeinterestedinclassicalhalo finders which calculated halo locations based on the maxi- findingonlyasopposedtorecoveringtidalstreams,thecode mumdensity peakwere more accuratethanFOF-based halo by default does not output halos where fewer than 50% of finders which use the averaged location of all halo particles the particlesare bound; this thresholdis user-adjustable,but (see also Gao & White 2006). The reason for this may be changing it does not produce statistically significant effects simply understood: as particle density rapidly drops in the on the clustering or mass function until halos with a bound outerreachesofahalo,thecorrespondingdispersionofparti- fractionoflessthan15%areincluded(see§4.5). clepositionsclimbsprecipitously. Consequently,ratherthan Wenotethat,inmajormergers,amorecarefulapproachto increasingthe statistical accuracyof the halo centercalcula- unbinding must be used. In many cases where merging ha- tion,includingtheparticlesatthehaloboundaryactuallyre- los initially have large velocity offsets, particles on the out- ducesit. The highest accuracyis instead achievedwhen the skirts of the halos can mix in phase space before the halo expected Poisson error (σ /√N) is minimized. As our halo x cores themselves merge. This results in many particles be- finderhasaccess(viathehierarchyofFOFsubgroups)tothe ingunboundwithrespecttoeitherofthetwohalocores,even innerregionsofthehalodensitydistribution,anaccuratecal- thoughtheyareboundtotheoverallmergingsystem.Assuch, culationofthecenterispossiblebyaveragingtheparticlelo- anaiveunbindingoftheparticleswouldleadtothe merging cationsfortheinnersubgroupwhichbestminimizesthePois- halos’massesdroppingsharplyforseveraltimestepspriorto sonerror.Typically,fora106particlehalo,thisestimatorends upaveragingthepositionsoftheinnermost103particles. 4Providedenoughinterest,wemayaddtheoptionofamulti-passunbind- The picture for halo velocities is not quite as simple. As ingprocedureinfutureversionsofthehalofinder. TheROCKSTARHaloFinder 7 the final merger.5 To counter this effect in major mergers, To help with cluster studies, we calculate several halo re- ROCKSTAR calculatesthe gravitationalpotentialofthe com- laxationparameters. Theseincludethecentralpositionoffset binedmergingsystem,ratherthanforindividualhalos,when (X , defined as the distance between the halo density peak off determiningwhethertounbindparticles.Thisbehavioristrig- andthehalocenter-of-mass),thecentralvelocityoffset(V , off geredbydefaultwhentwohaloshaveamergerratioof1:3or definedas the differencebetween the halo core velocity and larger;thisvalueisuser-adjustable,buthaslittleeffectonthe bulkvelocity),andtheratioofkinetictopotentialenergy( T ) |U| recoveredmassfunctionorclustering(see§4.5). forparticleswithinthehaloradius. Wereferinterestedread- ersto Netoetal. (2007)foradescriptionandcomparisonof 3.5.4. AdditionalHaloPropertiesandMergerTrees methodsfordetermininghalorelaxedness. Two more common outputs of halo finders are vmax, the We also calculate ellipsoidal shape parameters for halos. maximum circular velocity and Rs, the scale radius. vmax is Following the recommendations of Zemp et al. (2011), we simplytakenasthemaximumofthequantity GM(r)r- 1;it calculatethe massdistributiontensorforparticleswithinthe shouldbenotedthatthisquantityisrobustevenforthesmall- haloradius,excludingsubstructure: esthalosbecauseoftheextremelyshallowdepepndenceofv max onradius.ForR ,wedividehaloparticlesintoupto50radial 1 s M = xx (11) equal-massbins(withaminimumof15particlesperbin)and ij N i j directlyfitanNFW(Navarroetal.1997)profiletodetermine XN themaximum-likelihoodfit. The sorted eigenvalues of this matrix correspond to the We also calculate the Klypin scale radius for comparison squares of the principal ellipsoid axes (a2 > b2 > c2). We (Klypinetal.2011),whichusesv andM tocalculateR max vir s output both the axis ratios (b and c) as well as the largest undertheassumptionofanNFW(Navarroetal.1997)profile. a a Inparticular,forNFWprofiles,theradiusR atwhichv ellipsoidaxisvector,A~. max max is measured is a constantmultiple b of the radius R , and is Finally,wementionthatourhalofinderautomaticallycre- s givenby: ates particle-based merger trees. For a given halo, its de- d 1 scendant is assigned as the halo in the next timestep which db b- 1ln(1+b)- 1+b =0 (5) has the maximum number of particles in common (exclud- (cid:20) (cid:21) ingparticlesfromsubhalos). Whileitispossibletousethese Thismaybesolvednumerically,withtheresultthat merger trees directly, we recommend instead to use the ad- vanced merger tree algorithm discussed in Behroozi et al. R =2.1626R (6) max s (2011). This algorithm detects and corrects inconsistencies Instead of using R directly (which is ill-determined for acrosstimesteps(e.g.,haloswhichdisappearandreappearas max small halos), we make use of the ratio of R /R to relate theycrossthedetectionthreshold)tofurtherimprovethetem- max s themassenclosedwithin2.1626R tov ,M ,andthecon- poralconsistencyofthemergertrees. s max vir centrationc=R /R : vir s 4. TESTS&COMPARISONS c R 2.1626 f(c) =v2maxGMvir f(2.1626) (7) TheROCKSTARalgorithmhasalreadyundergoneextensive vir testing and comparison to other halo finders in Knebe et al. x f(x) ln(1+x)- (8) (2011). Intestswithgeneratedmockhalos, ROCKSTAR suc- ≡ 1+x cessfullyrecoveredhalopropertiesforhalosdownto20parti- Thus, by numericallyinvertingthe functionon the left-hand cles,inmanycases(e.g.,forhalomassandbulkvelocity)bet- sideofEq.7,cmaybefoundasafunctionofM andv , terthanallseventeenotherparticipatinghalofinders.Incases vir max and the Klypin scale radius R can be derived. The Klypin whereitdidnotperformbest,itwasoftenonlymarginallybe- s scaleradiusismorerobustthanthefittedscaleradiusforha- hindone ofthe otherphase-spacehalofinders. Notably, out los with less than 100 particles; this is due not only to shot ofallthehalofinders,itwastheonlyonetofullysuccessfully noise, butalsoduetothefactthathaloprofilesarenotwell- recoverallhaloproperties(mass,location,position,velocity, determined at distances comparable to the simulation force andv )forasubhalocoincidingwiththecenterofitshost max resolution. halo. In addition, Knebe et al. (2011)compared mass func- We additionally calculate the angular momentum of the tions, v functions, correlation functions (for r>2 Mpc), max halo(usingboundparticlesouttothedesiredhaloradius)and and halo bulk velocities for a cosmological simulation with thehalospinparameter(λ),asintroducedbyPeebles(1969): 10243 particles; ROCKSTAR had results comparableto other halofindersinalltheseresults,althoughonlytheotherphase- J E spacehalofindersshareditslowmasscompletenesslimit( λ= | | (9) GM2.5 25particlesforM ). ∼ pvir 200c We thus focus on comparisonsbeyond those already cov- whereJ isthemagnitudeofthehaloangularmomentumand ered in Knebe et al. (2011). Our comparisons cover results E is the total energyof the halo (potentialand kinetic). For forseveraldifferentdarkmattersimulations, brieflysumma- comparison, we also calculate the Bullock spin parameter rizedin§4.1. Weprovideavisualdemonstrationofthealgo- (Bullocketal.2001),definedas rithm’s performancein §4.2, a detailed comparisonwith the J J massandcorrelationfunctionsforotherhalofindersin§4.3, λB= = (10) anevaluationofthedynamicalaccuracyofhalopropertiesin √2MvirVvirRvir 2GMv3irRvir §4.4,andwepresentjustificationforourchoiceofthedefault 5 Withthe BDMhalofinder(Klypin etaql.2011), forexample, wehave parameters in §4.5. Finally, we show figures demonstrating observedhalomasseswhichdropbyafactorofthreeonaccountofthisaffect. theexcellentperformanceofthehalofinderin§4.6. 8 BEHROOZI,WECHSLER,&WU FIG.2.—ROCKSTARallowsrecoveryofevenveryclosemajormergers. Thisfigureshowsanexampleofamajormergerinvolving1013M⊙halosfromthe Bolshoisimulation(Klypinetal.2011).Thetoppanelshowsthecompleteparticledistributionaroundthemerginghalos.Inthesecondrow,theleftpanelshows thehostparticledistribution,andtherightpanelshowsthesubhaloparticledistribution,withparticlescoloredaccordingtosubhalomembership. (Theparticle plottingsizehasbeenincreasedtoshowmoreclearlytheextentofthesmallsubstructuresintheright-handpanel). Thetwodifferentcolorsintheleft-hand panelhintatthefactthatthereareindeedthreehalosinvolvedinthemajormerger,twoofwhichareextremelyclosetomerging.Theuniformsubhaloshapesin theright-handpanelsuggestthatsubhaloparticlescanbedistinguishedwithoutbiasdespiteextremevariationsinthehostparticledensitybetweenthesubhalo centersandthesubhalooutskirts. Thebottomrowshowsmoreclearlytheextremelyclosemajormerger. Thebottomleft-handpanelshowsthefullparticle distributioninpositionspaceinasmallregionclosetothemerginghalocores;here,thebimodaldistributionisevident,butdistinguishingparticlemembership isimpossiblebeyondtheimmediatevicinityofthecores. Ontheotherhand,thebottomright-handpanelshowsthesameparticlesinvelocityspace,wherethe bimodalityfromthetwocoresshowsaclearvelocityseparation,allowingparticlestobereasonablyassignedeveninthehalobulk. TheROCKSTARHaloFinder 9 10-1 Tinker et al. (2008) 10-2 420 Mpc h-1 10-3 1.1 640 Mpc h-1 -3-1c dex]1100--45 Tinker 1 12040000 MMppcc hh--11 Mp Tinker et al. (2008) f/ [10-6 420 Mpc h-1 f f 10-7 640 Mpc h-1 0.9 1000 Mpc h-1 10-8 2400 Mpc h-1 10-9 0.8 1011 1012 1013 1014 1015 1011 1012 1013 1014 1015 Halo Mass [M ] Halo Mass [M ] O• O• 10-1 BDM 10-2 1.1 Rockstar Tinker et al. (2008) -3-1c dex]1100--43 Rockstar 1 Mp BDM f / f [10-5 RTionckkesrt aert al. (2008) f 0.9 10-6 10-7 0.8 109 1010 1011 1012 1013 1014 1015 109 1010 1011 1012 1013 1014 1015 Halo Mass [M ] Halo Mass [M ] O• O• FIG.3.—Thehalomassfunctionfordistincthalosisverysimilartopreviouslypublishedresults. ThetoprowshowscomparisonstotheTinkeretal.(2008) centralhalomassfunctionforeachofthefourLasDamasboxes(seeTable1).Goodagreementisseenabove100particles.Thebottomrowshowsacomparison betweentheROCKSTARandBDM(Klypinetal.2011)halofindersontheBolshoisimulation(20483particles,250Mpch- 1).Theleft-handplotsshowthefull massfunctions,andtheright-handplotsshowtheresiduals,withPoissonerrorsshownfortheBolshoisimulation.AsnotedinTinkeretal.(2008),thecalibrated massrangedoesnotextendbelow1011h- 1M⊙;furthermore,theauthorsstate“thebehaviorofthefittingfunctionatlowermassesisarbitrary.”Wethereforedo notextrapolatetheTinkeretal.(2008)massfunctioninourcomparisons. 100 10-1 10-1 10-2 -1ex]10-2 -1ex]10-3 d d -3c 10-3 -3c p p M BDM (all) M10-4 BDM (sats) [10-4 Rockstar [ Rockstar f f 10-5 10-5 10-6 10-6 100 1000 100 1000 v [km s-1] v [km s-1] max max FIG.4.—Thehalovelocityfunctionisalsoverysimilartopreviouslypublishedresults. Thisfigureshowscomparisonsbetweenvelocity(vmax)functionsfor theROCKSTARhalofinderandBDMontheBolshoisimulation(20483particles,250Mpch- 1). Theleft-handplotshowsallhalos;theright-handplotshows satellitehalosonly. 10 BEHROOZI,WECHSLER,&WU 4.1. SimulationParameters We use five sets of simulations for this paper. The first of these, Bolshoi, has already been extensively detailed in 105 Klypinet al. (2011). To summarizethe relevantparameters, BDM (v > 300 km s-1) Mitipscah2- 10,4r8u3npuasritnicglethseimAuRlaTtisoinmouflactoiomnocvoindge(sKidreavletsnogvthet25al0. 104 Rockstarmax 1997) on the NASA Ames Pleiades supercluster. The as- sumedcosmologyisΩm=0.27,ΩΛ=0.73,h=0.7,ns=0.95, 103 and σ = 0.82, similar to WMAP7 results (Komatsu et al. x 2011); the effective force-softeninglength is 1 kpc h- 1, and 102 theparticlemassis1.36 108M h- 1. ⊙ × We also use four simulations of different sizes from the 101 LasDamas project (McBride et al, in preparation).6 These have11203 to14003 particlesincomovingregionsfrom420 100 Mpc h- 1 to 2400 Mpc h- 1 on a side, and were run using the 0.01 0.1 1 10 GADGET-2 code (Springel2005). Again, the assumed cos- r [Mpc h-1] mology is similar to WMAP7, with Ωm =0.25, ΩΛ =0.75, h=0.7, n =1.0, and σ =0.8. The effective force-softening s lengths range from 8 to 53 kpc h- 1, and the particle masses rangefrom1.87 109M⊙ h- 1 to 4.57 1011h- 1M⊙. Details 106 × × ofallthesimulationsaresummarizedinTable1. 105 BDM (vmax > 150 km s-1) Rockstar 4.2. VisualDemonstration Inordertodemonstratehowthealgorithmperformsonha- 104 losinarealsimulation,Fig.2showsanexampleofhowpar- ticles are assigned to halos in a major merger; this example 103 x hasbeentakenfromtheBolshoisimulationatz=0.Fromthe top image in the figure, one might expectthat two large ha- 102 losseparatedbyabout200kpch- 1 aremergingtogether,but 101 careful analysis reveals that the larger halo actually consists ofanothermajormergerwhereinthehalocoresareseparated by only 15 kpc h- 1. As shown in the bottom-left panel of 100 0.01 0.1 1 10 thefigure,theexistenceofthethirdmassivehaloisvisibleat r [Mpc h-1] a moderate significance level in position space—however, a position-spacehalofinderwouldhavenowaytocorrectlyas- signparticlesbeyondtheimmediatelocalityofthetwocores. Yet,inthebottom-rightpanel,theseparationofthetwohalo 106 coresisclearlydistinguishableinvelocityspace.Assuch,not onlycan the distinct existenceof the close-to-merginghalos 105 BDM (vmax > 100 km s-1) bereliablyconfirmed,butparticleassignmenttothetwohalos Rockstar basedonparticlevelocitycoordinatescanbereliablycarried 104 Subfind outaswell. We also remark that phase-space halo-finding helps im- 103 prove the accuracy of subhalo shapes by removingthe need x 102 to perform a position-space cut to distinguish host particles fromsubstructureparticles. Satellite halosare usuallyoffset 101 invelocityspacefromtheirhosts,butjustasimportantly,they usuallyalsohaveamuchlowervelocitydispersionthantheir 100 hosts. Thisimpliesthatsatellitesmaybereliablyfoundeven inthedensecoresofhalos—althoughthepositionspaceden- 10-1 0.01 0.1 1 10 sityisveryhighforhostparticles,theaveragevelocitydisper- r [Mpc h-1] sion is large enoughthat the lower-dispersionsubhalo parti- clescanbereliablydistinguishedfromhostparticles. Conse- FIG.5.—Two-pointcorrelation functions areverysimilartopreviously- quently,asshowninthemiddle-rightpanelofFig.2,satellites publishedresultsexceptrelativelyclosetothecentersofhalos. Thisfigure arepickedouteveninthedensehalocenterswithoutbiasas showscorrelationfunctionsfortheROCKSTARhalofinderandBDMonthe regardtoshapeorsize. Bolshoi simulation (20483 particles, 250 Mpch- 1). Thetop panel shows thecorrelation function forvmax>300kms- 1 (10,000halos), themiddle 4.3. MassandCorrelationFunctionComparisons apnadnetlhsehboowtstotmhepcaonrreellsahtioownsfuthnectcioonrrefolartivomnaxfu>nc1ti5o0nkfmorsv-m1a(x1>001,00000khmalos-s1),. Anextensivecomparisonofthemassfunctionandcorrela- ThebottompanelincludesacomparisonwithSubfind. Subfindhaloswere tion functionbetween ROCKSTAR and otherhalo findershas onlyavailablefor1/125thofBolshoi(atotalof∼1000halosforthisvmax threshold). Forafaircomparison,bothROCKSTARandBDMresultsforthe bottompanelwerecomputedonthesameregionofBolshoi. 6http://lss.phy.vanderbilt.edu/lasdamas/