ABSTRACT
We present a new algorithm for identifying dark matter halos, substructure, and tidal features. The approach is based on adaptive hierarchical refinement of friends-of-friends groups in six phase-space dimensions and one time dimension, which allows for robust (grid-independent, shape-independent, and noise-resilient) tracking of substructure; as such, it is named rockstar (Robust Overdensity Calculation 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 has shown rockstar to have excellent recovery of halo properties; we expand on these comparisons with more tests and higher-resolution simulations. We show a significant improvement in substructure recovery compared to several other halo finders and discuss the theoretical and practical limits of simulations in this regard. Finally, we present results that demonstrate conclusively that dark matter halo cores are not at rest relative to the halo bulk or substructure average velocities and have coherent velocity offsets across a wide range of halo masses and redshifts. For massive clusters, these offsets can be up to 350 km s−1 at z = 0 and even higher at high redshifts. Our implementation is publicly available at http://code.google.com/p/rockstar.
Export citation and abstract BibTeX RIS
1. INTRODUCTION
In the paradigm of lambda cold dark matter (ΛCDM), the majority of the matter density of the universe does not couple to electromagnetic fields, leaving it detectable only through its gravitational and possibly weak force interactions. Nonetheless, the effects of dark matter on the visible universe are spectacular, as the steep gravitational potentials of bound dark matter halos channel baryons together, forming the birthplaces of visible galaxies. In this model, the locations, sizes, and merging histories for galaxies are thus intricately connected to the growth of bound dark matter structures.
Testing this model in detail requires extensive computer simulations, as the complicated nonlinear evolution of structure growth cannot be fully evaluated by hand. The simulations generally follow the evolution of a set of dark matter particles and output the positions and velocities of the particles at several discrete time steps. These outputs must then be post-processed to determine the locations and properties of bound dark matter structures also known as "halos"—namely, the locations and properties that influence the formation of visible galaxies. This post-processing ("halo finding") necessarily involves both ambiguity and imprecision—ambiguity in the definitions (e.g., the center of a bound halo) and imprecision in determining halo properties due to limited information (e.g., for halos consisting of only a few dark matter particles, or for determining particle membership in two overlapping halos). Currently, as a consequence, essential statistical measures (e.g., the number density of halos with a given mass) are known to at best 5%–10% even for a specific cosmology (Tinker et al. 2008). Moreover, halo properties are rarely checked for consistency over multiple time steps, which can be a serious problem for robust modeling of galaxy formation theories (see, however, Behroozi et al. 2011; this has also been addressed in specific contexts when creating merger trees, for example, improving subhalo tracking and smoother mass accretion histories; see, e.g., Wechsler et al. 2002; Springel et al. 2005; Faltenbacher et al. 2005; Allgood et al. 2006; Harker et al. 2006; Tweed et al. 2009; Wetzel et al. 2009).
This level of uncertainty is unacceptable for current and future surveys, including data expected to come from the Baryon Oscillation Spectroscopic Survey (BOSS), Dark Energy Survey (DES), BigBOSS, Panoramic Survey Telescope and Rapid Response System (Pan-STARRS), Extended Roentgen Survey with an Imaging Telescope Array (eROSITA), Herschel, Planck, James Webb Space Telescope, and Large Synoptic Survey Telescope (LSST; Schlegel et al. 2009, 2011; The Dark Energy Survey Collaboration 2005; Chambers 2006; Predehl et al. 2010; Pilbratt et al. 2010; The Planck Collaboration 2006; Gardner et al. 2006; LSST Science Collaborations et al. 2009), in order to fully realize the constraining power of these observations. Indeed, derived quantities such as the halo mass function and autocorrelation functions must be understood at the one percent level in order for theoretical uncertainties to be at the same level as statistical uncertainties in constraining, e.g., dark energy (for example, Wu et al. 2010; Cunha & Evrard 2010). Although there may be a limit to the accuracy possible with dark matter simulations alone, given the impact of baryons on dark matter halo profiles (Stanek et al. 2009) and the power spectrum (Rudd et al. 2008), different halo finders running on the same simulation demonstrate that a large fraction of this uncertainty is still due to the process of halo finding itself (Knebe et al. 2011).
As previously mentioned, some of these uncertainties are due to limited use of information. Considerable progress has been made since the first generation of position-space halo finders (Davis et al. 1985; Lacey & Cole 1994), which used only particle locations to determine bound halos. Currently, the most advanced algorithms are adaptive phase-space finders (e.g., Maciejewski et al. 2009; Y. Ascasibar et al. 2012, in preparation), which make use of the full six dimensions of particle positions and velocities. At the same time, an often-overlooked (with the possible exception of Tweed et al. 2009) aspect of halo finding is the extra information stored in the temporal evolution of bound particles. In this paper, we detail an advanced adaptive phase-space and temporal finder that is designed to maximize consistency of halo properties across time steps, rather than just within a single simulation time step. Together with a companion paper (Behroozi et al. 2011), which details the process of comparing and validating halo catalogs across time steps to create gravitationally consistent merger trees, our combined approach is the first to use particle information in seven dimensions to determine halo catalogs, allowing unprecedented accuracy and precision in determining halo properties.
Furthermore, in contrast to previous grid-based adaptive phase-space approaches, ours is the first grid-independent and orientation-independent approach; it is also the first publicly available adaptive phase-space code designed to run on massively parallel systems on the very large simulations (>1010 particles) that are necessary to constrain structure formation across the range of scales probed by current and future galaxy observations. Finally, we remark that our halo finder is the first that is successfully able to probe substructure masses down to the very centers of host halos (see also Knebe et al. 2011), which is essential for a full modeling of galaxy statistics and will enable future studies of the expected breakdown between halo positions and galaxy positions due to the effects of baryon interactions at the very centers of galaxies.
Throughout, we have paid careful attention not only to the basic task of assigning particles to halos, but also to the process of estimating useful properties from them to compare with observations. While in many cases galaxy surveys are not able to probe halo properties to the same precision as halo finders in simulations, one significant counterexample exists. It is a common practice especially for halo finders based on the friends-of-friends (FOF) algorithm (e.g., Davis et al. 1985; Springel et al. 2001; Habib et al. 2009; Rasera et al. 2010; Gardner et al. 2007; see also Section 2.1) to calculate halo velocities by averaging all halo particle velocities together to find a bulk velocity. Examination of the difference between velocities in the inner regions of halos and the bulk averaged velocity suggests that the bulk average velocity may be offset by several hundred km s−1 from the velocity at the radius where galaxies are expected to reside; differences at this scale are easily detectable in cluster redshift surveys and may also factor in interpreting observations of the kinetic Sunyaev–Zel'dovich effect. As this difference has an important impact on the usefulness of derived halo properties, we additionally perform an investigation of the core-bulk velocity difference in halos across a wide range of redshifts and masses.
We begin this paper with a survey of previous work in halo finding, as well as previous limitations, in Section 2. We discuss our improved methodology in Section 3 and conduct detailed tests of our approach in Section 4. We present an analysis of the theoretical and practical limitations of simulations in terms of tracking substructure in Section 5. Finally, our results concerning the velocity offsets of cluster cores are presented in Section 6; we summarize our conclusions in Section 7. Multiple simulations including slight variations of cosmological parameters are considered in this paper; all simulations model a flat ΛCDM universe, and we always take the Hubble constant H0 to be 70 km s−1 Mpc−1; equivalently, h = 0.7.
2. PREVIOUS HALO FINDERS
2.1. Summary of Previous Approaches
Previously published approaches to halo finding may be classified, with few exceptions, into two large groups. Spherical overdensity (SO) finders, such as Amiga's Halo Finder (AHF; Knollmann & Knebe 2009), Adaptive Spherical Overdensity Halo Finder (ASOHF; Planelles & Quilis 2010), Bound Density Maxima (BDM; Klypin et al. 1999), SO (Lacey & Cole 1994; Jenkins et al. 2001; Evrard et al. 2002), parallel SO (pSO; Sutter & Ricker 2010), Voronoi Bound Zones (VOBOZ; Neyrinck et al. 2005), and Spline Kernel Interpolative Denmax (SKID; Stadel 2001), proceed by identifying density peaks in the particle distribution and then adding particles in spheres of increasing size around the peaks until the enclosed mass falls below a predetermined density threshold (a top-down approach). FOF and HOP-based (Eisenstein & Hut 1998) halo finders, such as FOF, SUBFIND, the LANL halo finder, parallel FOF (pFOF), Ntropy-fofsv, and AdaptaHOP (Davis et al. 1985; Springel et al. 2001; Habib et al. 2009; Rasera et al. 2010; Gardner et al. 2007; Tweed et al. 2009), collect particles together that fall above a certain density threshold and then, if so designed, search for substructure inside these particle collections (a bottom-up approach). Phase-space finders typically extend these two approaches to include particle velocity information, either by calculating a phase-space density, such as the Hierarchical Structure Finder (HSF; Maciejewski et al. 2009) and the Six-Dimensional Hierarchical Overdensity Tree (HOT6D; Y. Ascasibar et al. 2012, in preparation), or by using a phase-space linking length, as does Six-Dimensional Friends-of-Friends (6DFOF; Diemand et al. 2006).
There are three notable exceptions to these algorithms, namely, the ORIGAMI halo finder (discussed in Knebe et al. 2011), the Hierarchical Bound-Tracing algorithm (HBT; Han et al. 2011), and SURV (Giocoli et al. 2010). ORIGAMI operates by examining phase-space shell crossings for the current particle distribution as compared to the initial particle distribution; shells that have crossed along three dimensions are considered to be halos (as opposed to shells that have crossed along one or two dimensions, which would be considered as sheets and filaments, respectively). HBT uses an FOF approach to find distinct halos and uses particle lists from distinct halos at previous time steps to test for the presence of subhalos. SURV is a very similar algorithm, except that distinct halos are identified using spherical overdensities. These algorithms all rely heavily on temporal information in their approach to halo finding.
2.2. Limitations of Previous Algorithms
In order to develop an improved halo finder, it is important to understand some of the shortcomings of previous approaches. For the vast majority of halos, even the most basic of algorithms (FOF and SO) do an acceptable job of determining halo properties to 10% accuracy (Knebe et al. 2011). However, recent interest in the detailed properties and histories of halos—e.g., precision mass functions and merger trees and the shape of tidal structures—requires improvements to older approaches; this has resulted in a proliferation of new codes in the past decade (summarized in Knebe et al. 2011).
The most significant improvements to halo finding have come from using the information from six-dimensional (6D, position and velocity) phase space. Two traditional weak points for three-dimensional (3D, position-space) halo finders have been major mergers and subhalos close to the centers of their host halos. In both cases, the density contrast alone is not enough to distinguish which particles belong to which halo: when two halos are close enough, the assignment of particles to halos becomes essentially random in the overlap region. However, as long as the two halos have relative motion, 6D halo finders can use particle velocity-space information to very effectively determine particle–halo membership. This, coupled with the ability of 6D halo finders to find tidal remnants (which are condensed in phase space but not in position space), means that phase-space capabilities are required for the most accurate and interesting studies of dark matter halos.
At the same time, phase space presents a unique challenge. While position space and velocity space have well-defined distance metrics, there is not a single, unique way to combine position and velocity distance into a single metric. For a useful definition of phase-space distance, one needs to be able to decide, e.g., whether an object just passing by the origin at 1 km s−1 is closer to or farther than an object at rest 1 kpc away. One approach, used by 6DFOF, is to choose in advance a static conversion between velocity and position space. While simple, this approach seems somewhat self-limiting: if too short a linking length is chosen, the full extent of substructures cannot be found; if too large a linking length is chosen, then otherwise distinct substructures will be merged together.
A demonstrably superior approach, at least in terms of recovering halo properties (Knebe et al. 2011), is to adaptively define a phase-space metric. Both HSF and HOT6D subdivide the simulation space into 6D hyperboxes containing (at the maximum refinement level) as little as a single particle each. For a given particle, the enclosing hyperbox gives a local estimate of the phase-space metric, based on the relative sizes of the hyperbox's dimensions in position and velocity space. The usefulness of this estimate depends heavily on the method for partitioning space into hyperboxes; HSF uses, for example, an entropy-based approach to determine whether more information is contained in the particle locations for position space or velocity space.
These algorithms all give excellent results for identifying halo centers at a single time step. However, consistent halo catalogs across time steps are often compromised by a fundamental ambiguity in the definition of a host halo. For major mergers, it is often unclear which halo center represents the larger "host" or central halo and which represents the subhalo. Phase-space halo finding helps when the two halo centers are relatively far apart (i.e., weakly interacting) because there exists a strong correlation between the velocities of particles in the halo outskirts and halo centers. However, when the centers come close enough to interact strongly, this correlation is weakened, and it becomes much more difficult to accurately assign particles to the halos. As a result, it is much more difficult to determine which of the halo centers should be considered the host halo. Since the definition of halo mass often includes the mass of subhalos, this problem can result in large mass fluctuations across time steps for merging halos.
A number of solutions to this problem have been proposed and examined with the AdaptaHOP halo finder (Tweed et al. 2009). Tweed et al. (2009) found that a temporal approach (examining the host versus subhalo assignment at earlier time steps) was most successful at fixing this problem. Other methods, such as choosing the densest halo center to be the host halo, have inherent instabilities because of the spread in halo concentrations at a given halo mass. At the same time, while the Tweed et al. (2009) method successfully resolves this problem, it nonetheless only finds halos in position space and thus has the same weaknesses in identifying subhalo and major merger properties.
The HBT and SURV algorithms (Han et al. 2011; Giocoli et al. 2010) use the ingenious approach of tracing subhalos by using the particles found in previously distinct halos, which could potentially also solve many of these problems. Yet they both also include assumptions about accretion onto subhalos (e.g., that subhalos cannot accrete background matter from the host) that are untrue in HBT's case for large linking lengths (as halos will be identified as satellites far from the actual virial radius of the host) and for both halo finders with major mergers (where satellite and host are more ambiguous; they can in any case easily trade particles with each other). These assumptions vastly simplify the code at some expense to the completeness and accuracy of the mass function. More seriously, the design of the algorithms requires temporal information to find subhalos; in cases where simulation outputs are spaced very far apart or when only one time step is available, they cannot effectively find substructure. These issues are in principle fixable: future versions of the algorithms could easily combine advanced single-time step substructure finding with checks against previous time steps' particle lists.
3. AN IMPROVED APPROACH: ROCKSTAR
Our primary motivation in developing a halo finder was to improve the accuracy of halo merger trees that are required for an understanding of galaxy evolution. The design of our halo finder was thus motivated by a requirement for consistent accuracy across multiple time steps. This interest led to the development of a unique, adaptive phase-space temporal algorithm that is provably independent of halo orientation and velocity relative to the simulation axes and also attempts to be highly preserving of particle–halo and halo–subhalo membership across time steps. In addition, we paid special attention to the algorithm's efficiency and parallelizability, to allow it to run on the largest current data sets and so that it could easily scale to the next generation of simulations. Thus far, we have run the halo finder (and in many cases the partner merger tree code) on the Bolshoi (∼1010 particles; Klypin et al. 2011) and LasDamas simulations (200 boxes of (1–4) × 109 particles each; C. McBride et al. 2012, in preparation), on several 20483 simulations run to create DES simulated sky catalogs, on ∼100 high-resolution halos simulated as part of the RHAPSODY project (Wu et al. 2012a, 2012b), and on halos A-1 through A-5 of the Aquarius Simulation (up to 1.4 × 109 particles in a single halo; Springel et al. 2008).
As a first step, our algorithm performs a rapid variant of the 3D FOF method to find overdense regions, which are then distributed among processors for analysis (Section 3.1). Then, it builds a hierarchy of FOF subgroups in phase space by progressively and adaptively reducing the 6D linking length, so that a tunable fraction of particles are captured at each subgroup as compared to the immediate parent group (Section 3.2). Next, it converts this hierarchy of FOF groups into a list of particle memberships for halos (Section 3.3). It then computes the host halo/subhalo relationships among halos, using information from the previous time step if available (Section 3.4). Finally, it removes unbound particles from halos and calculates halo properties, before automatically generating particle-based merger trees (Section 3.5). A visual summary of these steps is shown in Figure 1.
3.1. Efficient, Parallel FOF Group Calculation
The 3D FOF algorithm has existed since at least 1985 (Davis et al. 1985). In principle, implementation of the algorithm is straightforward: one attaches two particles to the same group if they are within a prespecified linking length. Typically, this linking length is chosen in terms of a fraction b of the mean interparticle distance (in our code, as for others, the cube root of the mean particle volume); common values for generating halo catalogs range from b = 0.15 to b = 0.2 (More et al. 2011).
In practice, this means that one must determine neighbors for every particle within a sphere of radius equal to the linking length. Even with an efficient tree code (we use a custom binary space partitioning tree), this represents a great deal of wasted computation, especially in dense cluster cores. In such cases, particles might have tens of thousands of neighbors within a linking length, all of which will eventually end up in the same FOF group.
As the 3D FOF groups are used in our method only to divide up the simulation volume into manageable work units, we instead make use of a modified algorithm that is an order of magnitude faster. As is usual, neighboring particles are assigned to be in the same group if their distance is within the linking length. However, if a particle has more than a certain number of neighbors within the linking length (16, in our version), then the neighbor-finding process for those neighboring particles is skipped. Instead, neighbors for the original particle out to twice the original linking length are calculated. If any of those particles belong to another FOF group, that corresponding group is joined with that of the original particle. Thus, although the neighbor-finding process has been skipped for the nearest particles, groups that would have been linked through those intermediate neighbors are still joined together.
This process therefore links together at minimum the same particles as in the standard FOF algorithm—fine for our desired purpose—but does so much faster: neighbors must be calculated over a larger distance, but many fewer of those calculations must be performed. Indeed, rather than slowing down as the linking length is increased, our variation of FOF becomes faster. Because of this, we are free to choose an exceptionally large value for the linking length. A value of b = 0.2 is too small, as it does not include all particles out to the virial radius (More et al. 2011); after evaluating different choices for b (see Section 4.5), we chose b = 0.28, which guarantees that virial spherical overdensities can be determined for even the most ellipsoidal halos.
The most important parallelization work occurs at this stage. Separate reader tasks load particles from snapshot files. Depending on the number of available CPUs for analysis, a master process divides the simulation region into rectangular boundaries, and it directs the reader tasks to send particles within those boundaries to the appropriate analysis tasks.
Each analysis task first calculates 3D FOFs in its assigned analysis region, and FOFs that span processor boundaries are automatically stitched together. The FOF groups are then distributed for further phase-space analysis according to individual processor load. The load-balancing procedure is described in further detail in Appendix A. Currently, single 3D FOF groups are analyzed by at most one processor. Also, in the current implementation, there is no support for multiple particle masses, although support could easily be added by varying the linking length depending on particle mass. Provided enough interest, support for multiprocessor analysis of single large halos, as well as support for multiple particle masses, may be added in a future version of rockstar (Robust Overdensity Calculation using K-Space Topologically Adaptive Refinement).
3.2. The Phase-space FOF Hierarchy
For each 3D FOF group that is created in the previous step, the algorithm proceeds by building a hierarchy of FOF subgroups in phase space. Deeper levels of subgroups have a tighter linking-length criterion in phase space, which means that deeper levels correspond to increasingly tighter isodensity contours around peaks in the phase-space density distribution. This enables an easy way to distinguish separate substructures—above some threshold phase-space density, their particle distributions must be distinct in phase space; otherwise, it would be difficult to justify the separation into different structures.3
Beginning with a base FOF group, rockstar adaptively chooses a phase-space linking length based on the standard deviations of the particle distribution in position and velocity space. That is, for two particles p1 and p2 in the base group, the phase-space distance metric is defined as
where σx and σv are the particle position and velocity dispersions for the given FOF group; this is identical to the metric of Gottlöber (1998). For each particle, the distance to the nearest neighbor is computed; the phase-space linking length is then chosen such that a constant fraction f of the particles are linked together with at least one other particle. In large groups (>10,000 particles), where computing the nearest neighbor for all particles can be very costly, the nearest neighbors are only calculated for a random 10,000-particle subset of the group, as this is sufficient to determine the linking length to reasonable precision.
The proper choice of f is constrained by two considerations. If one chooses too large a value (f > 0.9), the algorithm will take much longer, and it can also find spurious (not statistically significant) subgroups. If one chooses too low of a value (f < 0.5), the algorithm may not find smaller substructures. As such, we use an intermediate value (f = 0.7), with the recommended minimum threshold for halo particles (20). In our tests of the mock NFW (Navarro et al. 1997) halos described in Knebe et al. (2011), this results in a false positive rate of 10% for 20-particle groups compared to a cosmological subhalo distribution, which declines exponentially for larger group sizes. These false positives are easily removed by the significance and boundedness tests described in Sections 3.3 and 3.5.3.
Once subgroups have been found in the base FOF group, this process is repeated. For each subgroup, the phase-space metric is recalculated, and a new linking length is selected such that a fraction f of the subgroup's particles are linked together into sub-subgroups. Group finding proceeds hierarchically in phase space until a predetermined minimum number of particles remain at the deepest level of the hierarchy. Here we set this minimum number to 10 particles, although halo properties are not robust approaching this minimum.
3.3. Converting FOF Subgroups to Halos
For each of the subgroups at the deepest level of the FOF hierarchy (corresponding to the local phase-space density maxima), a seed halo is generated. The algorithm then recursively analyzes higher levels of the hierarchy to assign particles to these seed halos until all particles in the original FOF group have been assigned. To prevent cases where noise gives rise to duplicated seed halos, we automatically calculate the Poisson uncertainty in seed halo positions and velocities and merge the two seed halos if their positions and velocities are within 10σ of the uncertainties. Specifically, the uncertainties are calculated as and , where σx and σv are the position and velocity dispersions and n is the number of particles, all for the smaller of the two seed halos. The two halos are merged if
In our tests, this threshold yields a near-featureless halo autocorrelation function; lower values result in a spurious upturn in the autocorrelation function close to the simulation force resolution.
For a parent group that contains only a single seed halo, all the particles in the group are assigned to the single seed. For a parent group that contains multiple seed halos, however, particles in the group are assigned to the closest seed halo in phase space. In this case, the phase-space metric is set by the seed halo properties, so that the distance between a halo h and a particle p is defined as
where σv is the seed halo's current velocity dispersion, vmax is its current maximum circular velocity (see Section 3.5), and "vir" specifies the virial overdensity, using the definition of ρvir from Bryan & Norman (1998), which corresponds to 360 times the background density at z = 0; however, other choices of this density can easily be applied.
Using the radius rdyn, vir as the position-space distance normalization may seem unusual at first, but the natural alternative (using σx) gives unstable and nonintuitive results. At fixed phase-space density, subhalos and tidal streams (which have lower velocity dispersions than the host halo) will have larger position-space dispersions than the host halo. Thus, if σx were used, particles in the outskirts of a halo could be easily misassigned to a subhalo instead of the host halo. Using rden, vir, on the other hand, prevents this problem by ensuring that particles assigned to subhalos cannot be too far from the main density peak even if they are close in velocity space.4 Intuitively, the largest effect of using rvir is that velocity-space information becomes the dominant method of distinguishing particle membership when two halos are within each other's virial radii.5
This process of particle assignment assures that substructure masses are calculated correctly independently of the choice of f, the fraction of particles present in each subgroup relative to its parent group. In addition, for a subhalo close to the center of its host halo, it assures that host particles are not misassigned to the subhalo—the central particles of the host will naturally be closer in phase space to the true host center than they are to the subhalo's center.
3.4. Calculating Substructure Membership
In addition to calculating particle–halo membership, it is also necessary to determine which halos are substructures of other halos. The most common definition of substructure is a bound halo contained within another, larger halo. Yet, as halo masses are commonly defined to include substructure, the question of which of two halos is the largest (and thus, which should be called a satellite of the other) can change depending on which substructures have been assigned to them. This is one of the largest sources of ambiguity between SO halo finders, even those that limit themselves to distinct halos.
We break the self-circularity by assigning satellite membership based on phase-space distances before calculating halo masses. Treating each halo center like a particle, we use the same metric as Equation (3) and calculate the distance to all other halos with larger numbers of assigned particles. The satellite halo in question is then assigned to be a subhalo of the closest larger halo within the same 3D FOF group, if one exists. If the halo catalog at an earlier time step is available, this assignment is modified to include temporal information. Halo cores at the current time step are associated with halos at the previous time step by finding the halo at the previous time step with the largest contribution to the current halo core's particle membership. Then, host–subhalo relationships are checked against the previous time step; if necessary, the choice of which halo is the host may be switched so as to preserve the host–subhalo relationship of the previous time step.
As explained above, these host–subhalo relationships are only used internally for calculating masses: particles assigned to the host are not counted within the mass of the subhalo, but particles within the subhalo are counted as part of the mass of the host halo. This choice assures that the mass of a halo will not suddenly change as it crosses the virial radius of a larger halo, and it provides more stable mass definitions in major mergers. Once halo masses have been calculated, the subhalo membership is recalculated according to the standard definition (subhalos are within rΔ of more massive host halos) when the merger trees are constructed.
For clarity, it should be noted that every density peak within the original FOF analysis group will correspond to either a host halo or a subhalo in the final catalog. It has been observed that FOF groups will "bridge" or "premerge" long before their corresponding SO halo counterparts (e.g., Klypin et al. 2011). However, as we calculate full SO properties associated with each density peak, a single FOF group is naturally allowed to contain multiple SO host halos; thus, the bridging or premerging of FOF groups does not affect the final halo catalogs.
3.5. Calculating Halo Properties and Merger Trees
Typically, several properties of interest are generated for halo catalogs. Regardless of the quality of particle assignment in the halo finder, careful attention to halo property calculation is essential for consistent, unbiased results.
3.5.1. Halo Positions and Velocities
For positions, Knebe et al. (2011) demonstrated that halo finders that calculated halo locations based on the maximum density peak were more accurate than FOF-based halo finders that use the averaged location of all halo particles (see also Gao & White 2006). The reason for this may be simply understood: as particle density rapidly drops in the outer reaches of a halo, the corresponding dispersion of particle positions climbs precipitously. Consequently, rather than increasing the statistical accuracy of the halo center calculation, including the particles at the halo boundary actually reduces it. The highest accuracy is instead achieved when the expected Poisson error () is minimized. As our halo finder has access (via the hierarchy of FOF subgroups) to the inner regions of the halo density distribution, an accurate calculation of the center is possible by averaging the particle locations for the inner subgroup that best minimizes the Poisson error. Typically, for a 106 particle halo, this estimator ends up averaging the positions of the innermost 103 particles.
The picture for halo velocities is not quite as simple. As noted in Section 6, the halo core can have a substantial velocity offset from the halo bulk. Since the galaxy hosted by the halo will presumably best track the halo core, we calculate the main velocity for the halo using the average particle velocity within the innermost 10% of the halo radius. For calculating the bound/unbound mass of the halo (see Section 3.5.2), however, we use the more appropriate averaged halo bulk velocity including substructure.
3.5.2. Halo Masses
For halo masses, rockstar calculates spherical overdensities according to multiple user-specified density thresholds: e.g., the virial threshold, from Bryan & Norman (1998), or a density threshold relative to the background or to the critical density. As is usual, these overdensities are calculated using all the particles for all the substructure contained in a halo. On the other hand, subhalo masses have traditionally been a major point of ambiguity (for density-space halo finders). With a phase-space halo finder, such as rockstar, the particles belonging to the subhalo can be more reliably isolated from the host, and thus less ambiguity exists: the same method of calculating spherical overdensities may be applied to just the particles belonging to the subhalo. In terms of the definition of where the subhalo "ends," Equation (3) implies that the subhalo edge is effectively where the distribution of its particles in phase space becomes equidistant from the subhalo and its host halo. If alternate mass definitions are necessary, the halo finder can output the full phase-space particle–halo assignments; these may then be post-processed by the user to obtain the desired masses.
3.5.3. Unbinding Particles
By default, rockstar performs an unbinding procedure before calculating halo mass and vmax, although this may be switched off for studies of, e.g., tidal remnants. Because the algorithm operates in phase space, the vast majority of halo particles assigned to central halos are actually bound. We find typical boundedness values of 98% at z = 0; see Section 4.5 and Behroozi et al. (2012a). Even for substructure, unbound particles typically correspond to tidal streams at the outskirts of the subhalo, making a complicated unbinding algorithm unnecessary. For this reason, as well as to improve consistency of halo masses across time steps, we perform only a single-pass unbinding procedure using a modified Barnes–Hut method to accurately calculate particle potentials (see Appendix B for details).6 Since many users will be interested in classical halo finding only as opposed to recovering tidal streams, the code by default does not output halos where fewer than 50% of the particles are bound; this threshold is user adjustable, but changing it does not produce statistically significant effects on the clustering or mass function until halos with a bound fraction of less than 15% are included (see Section 4.5).
We note that, in major mergers, a more careful approach to unbinding must be used. In many cases where merging halos initially have large velocity offsets, particles on the outskirts of the halos can mix in phase space before the halo cores themselves merge. This results in many particles being unbound with respect to either of the two halo cores, even though they are bound to the overall merging system. As such, a naive unbinding of the particles would lead to the merging halos' masses dropping sharply for several time steps prior to the final merger.7 To counter this effect in major mergers, rockstar calculates the gravitational potential of the combined merging system, rather than for individual halos, when determining whether to unbind particles. This behavior is triggered by default when two halos have a merger ratio of 1:3 or larger; this value is user adjustable but has little effect on the recovered mass function or clustering (see Section 4.5).
3.5.4. Additional Halo Properties and Merger Trees
Two more common outputs of halo finders are vmax, the maximum circular velocity, and Rs, the scale radius. vmax is simply taken as the maximum of the quantity ; it should be noted that this quantity is robust even for the smallest halos because of the extremely shallow dependence of vmax on radius. For Rs, we divide halo particles into up to 50 radial equal-mass bins (with a minimum of 15 particles per bin) and directly fit an NFW (Navarro et al. 1997) profile to determine the maximum-likelihood fit.
We also calculate the Klypin scale radius for comparison (Klypin et al. 2011), which uses vmax and Mvir to calculate Rs under the assumption of an NFW (Navarro et al. 1997) profile. In particular, for NFW profiles, the radius Rmax at which vmax is measured is a constant multiple b of the radius Rs and is given by
This may be solved numerically, with the result that
Instead of using Rmax directly (which is ill determined for small halos), we make use of the ratio of Rmax/Rs to relate the mass enclosed within 2.1626 Rs to vmax, Mvir, and the concentration c = Rvir/Rs:
Thus, by numerically inverting the function on the left-hand side of Equation (7), c may be found as a function of Mvir and vmax, and the Klypin scale radius Rs can be derived. The Klypin scale radius is more robust than the fitted scale radius for halos with less than 100 particles; this is due not only to shot noise but also to the fact that halo profiles are not well determined at distances comparable to the simulation force resolution.
We additionally calculate the angular momentum of the halo (using bound particles out to the desired halo radius) and the halo spin parameter (λ), as introduced by Peebles (1969):
where J is the magnitude of the halo angular momentum and E is the total energy of the halo (potential and kinetic). For comparison, we also calculate the Bullock spin parameter (Bullock et al. 2001), defined as
To help with cluster studies, we calculate several halo relaxation parameters. These include the central position offset (Xoff, defined as the distance between the halo density peak and the halo center of mass), the central velocity offset (Voff, defined as the difference between the halo core velocity and bulk velocity), and the ratio of kinetic to potential energy (T/|U|) for particles within the halo radius. We refer interested readers to Neto et al. (2007) for a description and comparison of methods for determining halo relaxedness.
We also calculate ellipsoidal shape parameters for halos. Following the recommendations of Zemp et al. (2011), we calculate the mass distribution tensor for particles within the halo radius, excluding substructure:
The sorted eigenvalues of this matrix correspond to the squares of the principal ellipsoid axes (a2 > b2 > c2). We output both the axis ratios (b/a and c/a) and the largest ellipsoid axis vector, A.
Finally, we mention that our halo finder automatically creates particle-based merger trees. For a given halo, its descendant is assigned as the halo in the next time step that has the maximum number of particles in common (excluding particles from subhalos). While it is possible to use these merger trees directly, we recommend instead to use the advanced merger tree algorithm discussed in Behroozi et al. (2011). This algorithm detects and corrects inconsistencies across time steps (e.g., halos that disappear and reappear as they cross the detection threshold) to further improve the temporal consistency of the merger trees.
4. TESTS AND COMPARISONS
The rockstar algorithm has already undergone extensive testing and comparison to other halo finders in Knebe et al. (2011). In tests with generated mock halos, rockstar successfully recovered halo properties for halos down to 20 particles, in many cases (e.g., for halo mass and bulk velocity) better than all 17 other participating halo finders. In cases where it did not perform best, it was often only marginally behind one of the other phase-space halo finders. Notably, out of all the halo finders, it was the only one to fully successfully recover all halo properties (mass, location, position, velocity, and vmax) for a subhalo coinciding with the center of its host halo. In addition, Knebe et al. (2011) compared mass functions, vmax functions, correlation functions (for r > 2 Mpc), and halo bulk velocities for a cosmological simulation with 10243 particles; rockstar had results comparable to other halo finders in all these results, although only the other phase-space halo finders shared its low-mass completeness limit (∼25 particles for M200c).
We thus focus on comparisons beyond those already covered in Knebe et al. (2011). Our comparisons cover results for several different dark matter simulations, briefly summarized in Section 4.1. We provide a visual demonstration of the algorithm's performance in Section 4.2, a detailed comparison with the mass and correlation functions for other halo finders in Section 4.3, an evaluation of the dynamical accuracy of halo properties in Section 4.4, and justification for our choice of the default parameters in Section 4.5. Finally, we show figures demonstrating the excellent performance of the halo finder in Section 4.6.
4.1. Simulation Parameters
We use five sets of simulations for this paper. The first of these, Bolshoi, has already been extensively detailed in Klypin et al. (2011). To summarize the relevant parameters, it is a 20483 particle simulation of comoving side length 250 Mpc h−1, run using the ART simulation code (Kravtsov et al. 1997) on the NASA Ames Pleiades supercluster. The assumed cosmology is Ωm = 0.27, ΩΛ = 0.73, h = 0.7, ns = 0.95, and σ = 0.82, similar to WMAP7 results (Komatsu et al. 2011); the effective force-softening length is 1 kpc h−1, and the particle mass is 1.36 × 108 M☉ h−1.
We also use four simulations of different sizes from the LasDamas project (C. McBride et al. 2012, in preparation).8 These have 11203 to 14003 particles in comoving regions from 420 Mpc h−1 to 2400 Mpc h−1 on a side and were run using the GADGET-2 code (Springel 2005). Again, the assumed cosmology is similar to WMAP7, with Ωm = 0.25, ΩΛ = 0.75, h = 0.7, ns = 1.0, and σ = 0.8. The effective force-softening lengths range from 8 to 53 kpc h−1, and the particle masses range from 1.87 × 109 M☉ h−1 to 4.57 × 1011 h−1 M☉. Details of all the simulations are summarized in Table 1.
Table 1. Simulation Parameters
Name | Particles | Size | Particle Mass | Softening | Ωm | ΩΛ | h | σ8 | ns | Code |
---|---|---|---|---|---|---|---|---|---|---|
Bolshoi | 20483 | 250 h−1 Mpc | 1.36 × 108 h−1 M☉ | 1 h−1 kpc | 0.27 | 0.73 | 0.7 | 0.82 | 0.95 | ART |
Consuelo | 14003 | 420 h−1 Mpc | 1.87 × 109 h−1 M☉ | 8 h−1 kpc | 0.25 | 0.75 | 0.7 | 0.8 | 1.0 | GADGET-2 |
Esmeralda | 12503 | 640 h−1 Mpc | 9.31 × 109 h−1 M☉ | 15 h−1 kpc | 0.25 | 0.75 | 0.7 | 0.8 | 1.0 | GADGET-2 |
Carmen | 11203 | 1000 h−1 Mpc | 4.94 × 1010 h−1 M☉ | 25 h−1 kpc | 0.25 | 0.75 | 0.7 | 0.8 | 1.0 | GADGET-2 |
Oriana | 12803 | 2400 h−1 Mpc | 4.57 × 1011 h−1 M☉ | 53 h−1 kpc | 0.25 | 0.75 | 0.7 | 0.8 | 1.0 | GADGET-2 |
Note. Descriptions are in Section 4.1.
Download table as: ASCIITypeset image
4.2. Visual Demonstration
In order to demonstrate how the algorithm performs on halos in a real simulation, Figure 2 shows an example of how particles are assigned to halos in a major merger; this example has been taken from the Bolshoi simulation at z = 0. From the top image in the figure, one might expect that two large halos separated by about 200 kpc h−1 are merging together, but careful analysis reveals that the larger halo actually consists of another major merger wherein the halo cores are separated by only 15 kpc h−1. As shown in the bottom-left panel of the figure, the existence of the third massive halo is visible at a moderate-significance level in position space—however, a position-space halo finder would have no way to correctly assign particles beyond the immediate locality of the two cores. Yet, in the bottom-right panel, the separation of the two halo cores is clearly distinguishable in velocity space. As such, not only can the distinct existence of the close-to-merging halos be reliably confirmed, but particle assignment to the two halos based on particle velocity coordinates can be reliably carried out as well.
Download figure:
Standard image High-resolution imageWe also remark that phase-space halo finding helps improve the accuracy of subhalo shapes by removing the need to perform a position-space cut to distinguish host particles from substructure particles. Satellite halos are usually offset in velocity space from their hosts, but just as importantly, they usually also have a much lower velocity dispersion than their hosts. This implies that satellites may be reliably found even in the dense cores of halos—although the position-space density is very high for host particles, the average velocity dispersion is large enough that the lower-dispersion subhalo particles can be reliably distinguished from host particles. Consequently, as shown in the middle-right panel of Figure 2, satellites are picked out even in the dense halo centers without bias as regard to shape or size.
4.3. Mass and Correlation Function Comparisons
An extensive comparison of the mass function and correlation function between rockstar and other halo finders has already been conducted in Knebe et al. (2011). Our algorithm has changed somewhat since that paper was published, and the mass and vmax function comparison in that paper did not separate out the effects of central halos from satellite halos. In addition, the correlation function comparison only compared the outskirts of halos beyond r = 2 Mpc h−1 and not smaller scales dominated by substructure.
In Figure 3, we present a comparison of the mass function for central halos in the Bolshoi and Las Damas simulations. The results from rockstar agree with those from Tinker et al. (2008) at all masses for halos with more than 100 particles on the level of a few percent, well within the calibration specification (5%) except when Poisson errors dominate. Detailed comparisons of the mass function from all of the Las Damas boxes, including 50 times the volume shown here, will be presented by C. McBride et al. (2012, in preparation). In the Bolshoi simulation, the mass function for central halos returned by rockstar is virtually identical to bdm (Klypin et al. 2011), differing by at most 5% over the entire mass range until Poisson errors dominate; these differences stem mainly from different conventions for unbinding particles, especially in major mergers.
Download figure:
Standard image High-resolution imageFigures 4 and 5 show comparisons between the vmax function for all halos and for satellite halos and the correlation function for vmax-selected halos between rockstar and the bdm algorithm for the Bolshoi simulation (Klypin et al. 2011). The vmax function for all halos is also mostly identical, differing by only 5% on average for halos above 100 km s−1. For satellites, the deviations are slightly more significant, especially for very massive halos. Most notably, as a phase-space halo finder, rockstar is able to track subhalos into the extreme inner reaches of halos. This enables it to track very massive subhalos even when their cores are very close to their hosts in position space, as long as they are sufficiently separated in velocity space. This is evident both in the increased number density of subhalos with large vmax as compared to bdm (Figure 4), but also in the increased one-halo contribution to the correlation function for massive halos (Figure 5). We also show comparisons for a small region of Bolshoi (1/125th of the total volume) where Subfind (Springel et al. 2001) halos were also available in Figure 5. For the halos considered in the latter comparison (vmax > 100 km s−1), bdm may overpredict the number of major mergers within 20 kpc h−1, whereas Subfind may underpredict the number of major mergers within 30 kpc h−1, given the deviations seen in the correlation function as compared to larger scales.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageFigure 6 shows the relationship between Mvir and vmax for both satellite halos and centrals using the rockstar halo finder on Bolshoi; as may be expected, satellites have a very similar relation as compared to centrals. Due to dynamical stripping, however, the satellite relation has more scatter and is biased toward lower halo masses at fixed vmax.
Download figure:
Standard image High-resolution image4.4. Dynamical Tests
In Knebe et al. (2011), the authors performed a series of tests on mock halos in order to assess the accuracy of halo property recovery. While rockstar performed extremely well in these tests, they nonetheless are not representative of the halos that one would expect to find in a real simulation (the tests considered only spherical, NFW/Plummer halos with very little substructure). This is understandable—in a real simulation, there is no a priori "correct" answer for the recovery of halo properties. Nonetheless, it is still possible to assess the precision of halo property recovery in a simulation by comparing halo properties between time steps.
To this end, we have used code from Behroozi et al. (2011) to check the halo property consistency for two position-space halo finders (bdm and Subfind; Klypin et al. 2011; Springel et al. 2001) and rockstar on the Bolshoi simulation. This code simulates the gravitational evolution of halos purely based on their recovered properties (mass, scale radius, position, and velocity). Given the positions and velocities at one time step, one can thus predict their positions and velocities at the next time step with high accuracy; the difference between the predicted and actual positions and velocities is one measure of the uncertainty in halo property recovery for the halo finder.
Figure 7 shows the results of this analysis. Both bdm and rockstar recover halo positions very precisely across time steps; although rockstar recovers positions better by a factor of ∼2 for lower halo masses, both perform close to the force (softening) resolution of the simulation. Subfind, where halo positions are averaged over all particles, performs poorly in comparison. rockstar recovers halo velocities much better (2–3 times as precisely) as compared to bdm, due largely to the fact that bdm uses only the innermost 100 particles to compute halo velocities (Klypin et al. 2011). By comparison, Subfind appears to do extremely well; however, because it averages particle velocities over the entire halo instead of over the central region, it suffers from substantial systematic errors (as discussed in Section 6), which make its actual performance in recovering estimates of galaxy velocities always worse than rockstar.
Download figure:
Standard image High-resolution image4.5. Evaluation of Default Parameters
In Figure 8, we show the residual mass and correlation functions for changes in the default parameter settings for rockstar. All mass and correlation functions were calculated at z = 0 for the Consuelo simulation. A summary of the main tunable parameters is shown in Table 2.
Download figure:
Standard image High-resolution imageTable 2. Summary of Algorithm Parameters
Variable | Default | Description | Section |
---|---|---|---|
b | 0.28 | Friends-of-friends linking length. | Section 3.1 |
f | 0.7 | Fraction of particles in each hierarchical refinement level. | Section 3.2 |
Δ(z) | Virial | Spherical overdensity definition | Section 3.3 |
BP | 1 | (Bound Properties) Whether halo properties are calculated only using bound particles. | Section 3.5.3 |
UT | 0.5 | (Unbound Threshold) The minimum ratio of bound to total mass required for halos. | Section 3.5.3 |
MP | 0.3 | (Major-merger Potential) Mass ratio for mergers at which particles are evaluated for boundedness relative to the merging system, as opposed to individual halos/subhalos. | Section 3.5.3 |
Download table as: ASCIITypeset image
In the top row of figures, it is evident that choosing a base 3D linking length of b = 0.2 does not capture particles all the way out to the virial radius for massive halos. Larger values of b (0.25–0.35) result in very similar mass functions; the standard value of b = 0.28 thus allows some degree of safety even for specific resimulations with unusual halo shapes. No effect on vmax results from choosing larger linking lengths, and thus the correlation function is unaffected. For smaller values of b, FOF groups become more fragmented; at low particle numbers, this can result in a single halo at high b being detected as multiple halos for a lower b.
In the second row, the effects of varying the refinement threshold for the 6D FOF hierarchy are shown. The default behavior is to retain 70% of particles from the next higher refinement level. If one retains more particles, then one is more likely to find low-significance objects—however, one is also more likely to find coincidental particle groups that do not correspond to actual halos. If one had a 90% particle retainment threshold, then approximately 4% of 30-particle halos (1011 M☉ for Consuelo) would be false positives; this is exactly the residual difference between the mass functions for a 90% threshold versus the standard threshold in Figure 8. However, a 50% threshold is somewhat too low, as several percent fewer halos are detected as compared with the standard threshold. For halos with masses above 100 particles, where percent-level comparisons are more trustworthy (Tinker et al. 2008), the mass functions are almost indistinguishable. Some effect is seen in the autocorrelation function for vmax > 240 km s−1 halos within 40 kpc. However, Section 5 suggests that rockstar with default parameters is already substantially complete for major mergers down to 10 kpc; as such, the additional halos found with f = 0.9 as opposed to the default f = 0.7 may represent false detections.
In the third row, the effects of calculating halo properties using all particles (as opposed to the default behavior of only using bound particles) are shown. When using all particles, halo masses increase by 1%–2% on average, implying that 98% of the initially assigned particles were bound. Including unbound particles affects satellite halos more than centrals, as evidenced by an increase in the correlation function in the regime where the one-halo term is dominant.
In the fourth row, the effects of changing the boundedness threshold requirement are shown. Lowering the threshold substantially (to UT = 0.1) does not result in any more halos being found; however, raising the threshold to UT = 0.7 results in fewer satellite halos being found and a slight decrement in the mass function.
In the fifth row, the effects of changing the major-merger potential threshold are shown; this has extremely little (<1%) effect on the mass function, and only a minor effect on the correlation function due to somewhat higher satellite masses for lower values of the threshold.
4.6. Performance
Figure 9 shows the excellent performance scaling of rockstar. The costliest part of our algorithm is the process of calculating the hierarchical refinement levels, which takes more time for halos with more substructure. For that reason, structure finding takes more time both for higher-resolution simulations and for late redshifts (the runtime per snapshot scales approximately as T∝a0.84). Nonetheless, the halo finder is still so efficient (5–10 CPU hours per billion particles at z = 0) that high-throughput access to particle data is important to avoid wasting CPU time. Indeed, despite its advanced capabilities, it is roughly 4–5 times as efficient as other halo finders that include substructure with published performance data (Knollmann & Knebe 2009; Springel et al. 2010), as shown in Figure 9. Some caution is necessary in comparing these timings directly, as the systems used were not the same. Nonetheless, the systems used for calculating rockstar's performance (Sun X2200 systems with 2 × 2.3 GHz AMD Opteron 2376 processors, 32 GB of memory, and Cisco DDR HCA Infiniband adapters accessing a Lustre Filesystem) are over five years old, negating any advantage from newer technology. The most modern FOF halo finders (e.g., Woodring et al. 2011) are faster by roughly a factor of two than rockstar on large data sets; however, these cannot identify substructure, nor can they produce complete catalogs of halos with SO masses.
Download figure:
Standard image High-resolution image5. ESTIMATING SATELLITE HALO COMPLETENESS LIMITS
For many decades, early computational simulations suffered from the so-called overmerging problem, where satellite halos were found to disappear almost as soon as they passed within the boundary of a larger cluster (see Klypin et al. 1999 for a discussion). This situation has improved dramatically on account of increased mass and force resolution in simulations; however, even in modern simulations there is often a need to add "orphan" satellite galaxies (galaxies that exist without a corresponding subhalo) in order to match the small-scale clustering seen in observations (Kitzbichler & White 2008). In detail, this depends on both simulation resolution and halo finding; for example, using the halo finder presented here, the clustering on scales larger than ∼100 kpc for galaxies in the Sloan Digital Sky Survey can be reproduced using only resolved halos and subhalos in the Bolshoi simulation (Reddick et al. 2012). With recent observations pushing the previous small-scale limit of clustering measurements from 100 kpc to 10 kpc (Tinker et al. 2012; Jiang et al. 2011), and with increased accuracy needs for modeling galaxy populations in clusters, it is important to understand the limitations of cosmological simulations' abilities to recover substructure at this level.
In observations, the distribution of satellite galaxies around such massive hosts is completely self-similar (density proportional to a power law of radius, with an exponent between −1.7 and −1.5; Tinker et al. 2012). However, Figure 10 acutely shows the tension between this observational result and what is found in the Bolshoi simulation for cluster-sized halos (Rvir ∼ 1 Mpc); in the Bolshoi simulation, there is a clear radial incompleteness scale in the recovered volume density for merger ratios below 1:30.
Download figure:
Standard image High-resolution imageMaking a quantitative estimate of the radial incompleteness scale is difficult because the true satellite halo distribution, ρ(r), is unknown. We can nonetheless connect the luminosity constraints in Tinker et al. (2012) to constraints on dark matter halos via the assumption that galaxy luminosity is tightly correlated with halo peak mass (i.e., the largest mass that progenitors of a given halo or subhalo have ever had). This assumption has been tested in a large number of studies and been found to accurately reproduce both the halo mass–stellar mass relation and the luminosity-dependent clustering of galaxies (Behroozi et al. 2010, 2012b; Reddick et al. 2012; Moster et al. 2010; Guo et al. 2010; Conroy & Wechsler 2009; Wang & Jing 2010; Conroy et al. 2006). As concerns halos, this assumption would imply that the radial dependence of satellite halos with a given peak mass should also follow a power law with an exponent between −1.7 and −1.5.
We thus conservatively define the satellite halo incompleteness scale as the radius at which the logarithmic slope of the satellite halo density ρ(r) becomes shallower than −1.5; steeper numbers could overestimate the true logarithmic slope of the satellite halo radial distribution. However, we include the possibility that the true cutoff slope should be −1.7 in our treatment of the errors. We show the results for the radial completeness scale as a function of satellite mass ratio (satellite halo peak mass compared to host mass) for massive hosts (1013.5 M☉ < M < 1014 M☉) in Figure 10. We also include results for the bdm halo finder in Figure 10 to show the comparison with non-phase-space halo finding. The benefit of phase-space finding is striking for satellite halo mass ratios above 1:30, with rockstar able to find major mergers consistently down to a tiny fraction of the virial radius. However, for satellites with smaller mass ratios, rockstar is unable to perform any better than bdm. This finding is especially unexpected considering the results in Knebe et al. (2011), wherein rockstar was the only halo finder able to accurately recover all the halo properties of a satellite (mass ratio 1:100) placed directly at the center of a million-particle 1014 M☉ halo.
In order to explain such a drastic reduction in ability to find satellites, it is instructive to consider the effects of tidal stripping on a satellite in a simulation with limited force and mass resolution. In particular, due to gravitational softening, the maximum density of a halo cannot continue increasing indefinitely toward its center; instead, it will threshold at some distance fres from the center. For halos with a low enough mass, there may not be any particles within fres of the center; for these halos, the effective maximum density is degraded even further. Due to this effect, satellites are vulnerable to tidal disruption much earlier for larger values of fres and for larger particle masses. In particular, the fluid Roche limit (under the assumption of a significantly more massive host) is given by
where dRoche is the radius within which tidal disruption occurs, RH is the host radius, ρH is the enclosed averaged host density, and ρs is the average satellite density. In terms of halos, this means that tidal disruption will occur at a distance R from the host center where the enclosed host halo average density is 2.4−3 ∼ 7% of the satellite halo peak density ρs.
A naive application of the Roche limit under the assumption of a spherical NFW profile, circular satellite halo orbits, the mean mass–concentration relation from Bolshoi, and the mass and force resolutions of Bolshoi yields an estimate of the theoretical disintegration limit shown in Figure 10. A number of significant limitations apply to this calculation, because it does not include any mass stripping for satellites after they enter the host, and it does not include any effects from eccentric orbits, where satellites at a given radius would be expected to have passed through closer radii many times. The latter bias is stronger for low-mass satellites, which experience less dynamical friction and so are expected to orbit many more times before destruction than larger halos; as such, low-mass satellites at a given radius are on average more stripped and have had more exposure to high tidal fields than high-mass satellites at the same radius. This calculation also does not include numerical simulation artifacts wherein forces on one part of a satellite may be calculated differently from forces on another part (e.g., direct summation versus multipole expansion). Nonetheless, many of these biases go in the direction of earlier destruction of satellites; as such, the calculation represents a useful lower limit on the radius where satellites may exist in massive clusters. For a much more detailed assessment of subhalo completeness, we refer the reader to H. Wu et al. (2012, in preparation).
This calculation and the results in Figure 10 would suggest that, regardless of the halo-finding technique used, simulations of clusters still suffer from an overmerging at the very innermost radii due to limited force and mass resolution. Depending on the scientific questions being addressed, these incomplete halos may be more or less relevant. For a high-resolution simulation like Bolshoi, the incompleteness does not impact the correlation function on scales larger than ∼100 kpc and is present only in the inner regions of massive halos (see, e.g., Reddick et al. 2012). In general, when determining the desired mass and force resolution required, the requirements to resolve satellites in the inner regions of massive clusters will be much more stringent. In these regions, tracking galaxies after their halos are destroyed with "orphans" may represent an attractive solution to reduce computational requirements, but the need for these galaxies is somewhat mitigated with increased force resolution and with the effective halo finding in inner regions available with rockstar.
6. HALO CORE VELOCITY OFFSETS
In the ΛCDM paradigm, the main method by which halos reach a relaxed state is dynamical friction, namely, energy in bulk motions relative to the halo center of gravity is transformed into increased halo velocity dispersion. Dynamical friction depends not only on the background density that a satellite halo is passing through, but also on the satellite mass and velocity. For a massive host halo, the high background density means that incoming satellites will initially transfer momentum to the host with high efficiency. However, massive host halos also have very strong tidal fields, and once satellites are disrupted into cold velocity streams, dynamical friction becomes much less effective at transferring momentum into the inner halo regions. This combination of dynamical friction and tidal forces suggests that momentum transfer may be more efficient in the outer regions of a host halo and less efficient in the inner regions, leading to an offset between the mean velocity of the central density peak and the bulk velocity of the halo.
We examine the presence of this effect by calculating radially averaged halo particle velocities (i.e., particle velocities averaged in spherical shells) and comparing to the halo bulk velocity (averaged over all halo particles) for a large number of central halos in both the Bolshoi and Consuelo simulations over a range of halo masses (1012–1014 M☉) and redshifts (z = 0–1.4). We additionally consider the effects of excluding particles belonging to substructure from calculations of the radially binned velocities; results for all of the median offsets (both including and excluding substructure) are shown in Figure 11.
Download figure:
Standard image High-resolution imageFigure 11 demonstrates that the core-bulk velocity difference can be quite dramatic at high redshift for massive halos: halos with 1014 M☉ < M < 3 × 1014 M☉ at z = 1.4 have a median offset of over 450 km s−1 between the velocity averaged within 0.1 Rvir (excluding substructure) and the bulk velocity of the halo. This effect is lower for both smaller halos and later times, presumably due to lower velocity offsets in incoming substructure for smaller halos and a reduced merger rate at later times. Nonetheless, a clear difference between the inner and bulk velocities is present in all halos tested even down to z = 0. In all cases, there is a definite radial trend, with the velocity in the inner region of the halo being farther away from the bulk velocity than the halo outskirts, consistent with the hypothesis of reduced efficiency of momentum transfer in the inner regions of the halo.
In the calculations that exclude substructure, a clear plateau is evident below 0.1–0.2 Rvir where the core velocity offset stabilizes, at least for halos with M > 1013 M☉; we have insufficient mass resolution to test this effect robustly for smaller halos. This transition is extremely consistent with expectations for where most of the satellite mass will be stripped in such halos. For a satellite halo with a concentration of c = 10 falling into a massive host with c = 5, Equation (12) would suggest that 90% of the satellite mass (and momentum) will be stripped from the main satellite density peak at radii greater than 0.3 Rvir, host, with another 9% stripped by 0.1 Rvir, host. As tidal forces increase dramatically toward the center of a halo, they should disrupt efficient momentum transfer within 0.1 Rvir, host even for satellites on highly radial orbits.
In the calculations that include substructure, the median velocity offsets are less than the results that exclude substructure because substructure is included when calculating the bulk halo velocity. In the very innermost regions of halos, however, both radial velocity averages converge, suggesting that our results should be robust to the subhalo-finding technique used. On the outskirts of halos, the presence of non-disrupted subhalos results in an upturn in the offset between the radially averaged velocity (including substructure) and the bulk halo velocity.
Due to significantly higher mass resolution, the Bolshoi simulation is able to probe the velocity offsets to significantly smaller radii as compared to the Consuelo simulation. Furthermore, the Bolshoi simulation is able to better resolve substructure than Consuelo; as such, the difference between radial velocity offset profiles that include and exclude substructure in Consuelo is less than the difference for Bolshoi. Except for this caveat, the median velocity offsets are in excellent agreement between the two simulations.
One more interesting feature of Figure 11 is that the median velocity offset never reaches zero for any of the halos in either simulation. That is to say, there is no single radius where the bulk-averaged halo velocity corresponds to the actual average motion of particles in a radial shell. We examine this issue further by considering 10 halos randomly drawn from the Bolshoi simulation at z = 0 in the mass range from 1014 to 3 × 1014 M☉. As before, we compute the average velocity in radial bins for each halo, but instead of the absolute offset from the bulk velocity, we show just the offsets between the x-components of the two velocities in Figure 12.
Download figure:
Standard image High-resolution imageThe results for individual halos in Figure 12 are just as striking as for the median offsets in Figure 11. Almost every halo in Figure 11 exhibits the same velocity offset plateau within 0.1 Rvir, suggesting that the velocity within 0.1 Rvir has a physical interpretation corresponding to the actual bulk motion of those particles. If substructure is excluded, the radially averaged x-velocity often never matches the bulk x-velocity, implying an overall offset between the main halo motion and satellites. If substructure is included, the radially averaged x-velocity can briefly match the bulk x-velocity (and it must do so, by the mean value theorem), but there is no reason why the radially averaged velocities along other axes must match the bulk velocity at the same time; in fact, as shown for the z-velocities in the bottom panel of Figure 12, they do not. As such, the bulk velocity of the particles in the halos shown here does not actually correspond to the velocity averaged over any of the radial shells.
To investigate the latter point in more detail, we show conditional density plots of the absolute offsets between radially averaged velocities and the bulk velocity for the complete sample of 1014 to 3 × 1014 M☉ (again at z = 0 in Bolshoi) in Figure 13. While there is significant scatter in the absolute offsets (on the order of 0.3 dex, after correcting for the variation over the bin width), only a tiny fraction (≪1%) of the halos ever have a radial bin in which the averaged velocity matches the bulk velocity to better than 20 km s−1, and most have significant offsets (∼100 km s−1) at all radii from the halo center.
Download figure:
Standard image High-resolution imageThese results support the use of the Poisson minimization method of Section 3.5.1 as a velocity estimator for rockstar, as the method primarily probes the inner regions of halos with velocity profiles like those shown in Figure 12. Indeed, for both the Bolshoi and Consuelo simulations, we find an extremely consistent mean velocity offset between the core velocity (as discussed in Section 3.5.1) and the halo bulk velocity, as shown in Figure 14. Significantly, these offsets are much larger than either the direct Poisson error estimates in determining velocities (Section 3.5.1) or the cross-time step velocity error estimates in Figure 7 (Section 4.4). The main differences between Bolshoi and Consuelo come from reduced mass resolution and substructure resolution in the latter simulation, resulting in a velocity estimate slightly closer to the bulk average for Consuelo as compared to Bolshoi.
Download figure:
Standard image High-resolution imageThese velocity offsets display a power-law-like behavior across more than four orders of magnitude in halo mass; however, this power law is not the same as the power law for the average velocity dispersion σv (which is proportional to M1/3); instead, the velocity offset relative to the velocity dispersion increases for increasing halo masses. At z = 0, the average core-velocity offset may be fit by the empirical relation
This relation implies that, for the most massive M ∼ 2 × 1015 M☉ clusters at z = 0, the average core-bulk velocity offset is of order 350 km s−1.
We note that Gao & White (2006) have previously studied the question of core-bulk velocity offsets in the Milky Way to cluster-size halos in the Millennium Simulation. In their work, the core velocity is defined as the average velocity of the most-bound 100–1000 particles as opposed to being defined directly by the halo finder. As may be expected from our finding that the core-bulk velocity offset is relatively constant within 0.1 Rvir, the choice of the number of particles to use in the definition of core velocity only impacts their results for the smallest halos in their analysis, which have only ∼2000–4000 particles to begin with. They find that the core-bulk velocity offsets are approximately 10%–15% of V200 (defined as ), increasing with increasing halo mass. As V200 and σv both scale as M1/3 and are approximately the same magnitude across a wide range of halo masses at z = 0, this compares well with our result in Figure 14 that the core-bulk offsets are on the order of 10%–20% of σv, also increasing with increasing halo mass.
Interestingly, Gao & White (2006) also explore the accretion histories of halo cores; they find that 75% of cluster-scale halos have accreted less than 10% of their core particles since z = 0.5. By contrast, cluster-scale halos have on average accreted more than 40% of their bulk mass in that same time period (Wechsler et al. 2002). While this does not by itself support our hypothesis that momentum transfer to the core becomes inefficient due to tidal stripping, it leads to the same conclusion that accreted material makes its way only very slowly to the halo core.
7. CONCLUSIONS
We have presented a novel method of finding halos in phase space, the rockstar halo finder. This method has several key strengths that, in combination, improve on previous approaches to robustly identifying halos and their substructures:
- 1.High accuracy of recovered halo properties (as discussed in Knebe et al. 2011).
- 2.The ability to consistently and usefully define a subhalo mass (Section 3.5.2).
- 3.
- 4.Explicit grid independence, orientation independence, and halo shape independence for recovered halos (Section 3).
- 5.Calibrated estimates of the accuracy of recovered halo properties in full cosmological simulations (Section 4.4).
- 6.Extremely efficient resource consumption in both CPU time and memory (Section 4.6).
- 7.Massively parallel implementation, supporting large future simulations (Appendix A).
- 8.A publicly available codebase:http://code.google.com/p/rockstar.
In addition, the halo finder complements the method of Behroozi et al. (2011) for creating merger trees, making it part of the first halo-finding system to recover halos using seven dimensions of information (phase space plus time). These properties make it ideal for a large number of scientific applications, especially for studies of subhalos and subhalo particle distributions. We have demonstrated the ability to probe substructure into the very inner regions of halos, at least as far as simulations are capable of resolving. This is critical for comparison with galaxy statistics in dense regimes and paves the way for direct comparisons of clustering between simulations and observations in the regions where the unknown effects of baryons are the most significant.
We have also investigated the radial dependence of halo velocity. We find that there exists a well-defined core velocity (for r < 0.1 Rvir) that is always distinct from the bulk-averaged halo velocity. Using the bulk velocity to estimate galaxy peculiar velocities thus carries significant systematic error, which can be up to 350 km s−1 in the most massive clusters at z = 0 and is larger at higher redshifts. Furthermore (as discussed in Section 6), the halo velocity averaged in radial bins never matches the bulk average velocity; as such, the bulk average velocity is only an ensemble property of the entire halo, and the local dynamics have a much more complicated and interesting structure than can be described by the bulk velocity and the velocity dispersion alone.
We thank Anatoly Klypin and Joel Primack for providing access to the Bolshoi simulation, which was run on the NASA Advanced Supercomputing Pleiades computer at NASA Ames Research Center, Michael Busha for running Subfind on the Bolshoi simulation, and Michael Busha, Cameron McBride, and the rest of the LasDamas collaboration for providing access to the Las Damas simulations, which were run on the Orange cluster at SLAC and on the NSF Teragrid machine Ranger (PI: Andreas Berlind). We are grateful to Anatoly Klypin, Joel Primack, Michael Busha, Houjun Mo, Alexie Leauthaud, Andrey Kravtsov, James Bullock, Alex Knebe, Rachel Reddick, Yu Lu, Oliver Hahn, Matt Becker, Louie Strigari, Jeremy Tinker, Tom Abel, Steven Rieder, Markus Haider, Shea Garrison-Kimmel, Gus Evrard, Marcelo Alvarez, Joanne Cohn, Martin White, Simon White, Manodeep Sinha, and Annika Peters for interesting discussions, helpful suggestions, and testing. We gratefully acknowledge the support of Stuart Marshall and the SLAC computational team, as well as the computational resources at SLAC. Finally, we thank our anonymous referee for a very careful and thorough review, which resulted in a great many improvements to the clarity of this paper. P.S.B. and R.H.W. received support from NASA HST Theory grant HST-AR-12159.01-A; R.H.W. was also supported by the National Science Foundation under grant NSF AST-0908883. This work was additionally supported by the U.S. Department of Energy under contract number DE-AC02-76SF00515.
APPENDIX A: LOAD BALANCING
The simplest approach to load balancing—equal volumes—works well on scales where the universe is homogenous, namely, when individual analysis regions are at least 100 Mpc on a side. However, cosmic variance on smaller scales requires a more sophisticated approach for algorithms designed to run on small simulations or very large numbers of processors.
The approach we use is to divide the simulation volume into chunks requiring approximately equal memory. We divide the number of available processors, N, into three factors (Nx, Ny, and Nz), which control the number of divisions along each of the principal axes of the simulation. Naturally, Nx, Ny, and Nz are chosen to be as close as possible to N1/3 so as to minimize the surface-to-volume ratio of the divisions. Initially, particles are binned into B = 104 bins according to their x-coordinate, and the simulation is divided into Nx regions of equal particle count along the x-axis. Then, within each of these regions, particles are binned according to their y-coordinate, and each region is split accordingly along the y-axis into Ny subregions of equal particle count. Finally, each subregion is split in the same way into Nz sub-subregions along the z-axis. This process is accurate and efficient, requiring a fixed total data transfer size of , an amount that is completely independent of the total number of particles.
3D FOF group calculation is performed in parallel on these regions without any inter-process communication. However, for phase-space analysis, the number of recovered groups can vary considerably depending on the overall density of the analysis region (e.g., voids will have lower fractions of bound particles). For that reason, FOFs are distributed to processors in small work units; as soon as a processor finishes one work unit, it receives another to work on (possibly from a completely different part of the simulation), so that maximum concurrency can be maintained at all times. However, a single FOF group is never (in the current implementation) split up for analysis among multiple processors. Thus, the analysis time for the largest FOF group sets the lower limit for the wall clock time taken to analyze the simulation snapshot, regardless of the number of processors available for use.
APPENDIX B: UNBINDING
To calculate which particles are bound to halos, we use a modified Barnes–Hut algorithm (Barnes & Hut 1986) to calculate particle potential energies with a binary space partitioning (BSP) tree.9 As noted in Salmon & Warren (1994), the cell refinement criterion suggested by Barnes and Hut gives unacceptably high maximum errors (as opposed to average errors); in addition, our use of a BSP tree (which gives fewer wasted refinement levels than an octree) renders the Barnes and Hut criterion somewhat inapplicable. Instead, given a cell containing n particles with mass center x0, we calculate a threshold monopole approximation distance—that is, the distance from x0 at which a monopole would be an appropriate approximation of the potential. Following Salmon & Warren (1994), we can estimate the relative error θ in the calculated potential as a function of distance d from the center x0 as being bounded by
where xi and mi are the locations and masses of the n particles in the cell, respectively, and rmax is the maximum distance between x0 and any of the particles in the cell. Thus, for a given choice of the relative error θ, this equation may be inverted to give a minimum distance for acceptable use of the monopole approximation:
This equation is valid regardless of the boundary sizes of the cell; hence, it is appropriate for use even with a BSP tree. We choose θ to give potentials accurate to 4%; from our tests, this is sufficient to correctly estimate the boundedness of 99.8% of all halo particles.
Footnotes
- 3
This would not be true if, e.g., halos had Plummer profiles or otherwise flat density profiles in their centers.
- 4
For determination of tidal streams, this "problem" becomes a "feature," and use of σx may be preferable to rvir.
- 5
An alternate radius (e.g., r200b or r500c) could be used instead, but it would have an effect only on a small fraction of particles in a small fraction of halos (major mergers).
- 6
Provided enough interest, we may add the option of a multi-pass unbinding procedure in future versions of the halo finder.
- 7
With the bdm halo finder (Klypin et al. 2011), for example, we have observed halo masses that drop by a factor of three on account of this affect.
- 8
- 9
Several other halo finders, such as bdm and AHF, use the assumption of spherical symmetry to calculate particle potentials; we find that, while this works well for the vast majority of halos, the potential calculated in this way can be dramatically wrong for halos undergoing major mergers, especially if the halo center is incorrectly identified.