The Splashback Radius of Halos from Particle Dynamics. III. Halo Catalogs, Merger Trees, and Host–Subhalo Relations

The original purpose of the Sparta code was to calculate the splashback radius from the dynamics of individual particles. Sparta tracks particles as they approach a halo and records the time and location of their infall and their first apocenter (or "splashback event"). We tag particles at infall if they are deemed to belong to a subhalo. Subhalos whose mass is greater than 1/100 of the host mass may suffer from significant dynamical friction and are thus biased in their splashback radii. After removing particles from such subhalos from consideration, we smooth the distribution of the remaining particle splashbacks in time and define the halo's splashback radius as their mean, median, or higher percentiles (Figure 1). The same procedure is repeated for the mass enclosed within the radii of the particle splashbacks (meaning that the splashback mass can slightly deviate from the mass within the splashback radius).

Since the publication of Paper II, we have made a number of improvements to the Sparta algorithm. First, we have fixed a bug due to which the Hubble drag term was underestimated by a factor of h (Croton 2013), leading to earlier splashback times and splashback radii that were a few percent too large (because halos grow over time). Furthermore, the Paper I version of Sparta could not compute ${R}_{\mathrm{sp}}$ for up to 5% of halos depending on the mass range. We have increased the completeness (the fraction of halos with valid splashback measurement) of the ${R}_{\mathrm{sp}}$ catalogs in a number of ways. First, the previous version of Sparta stopped tracking all particles as soon as a halo became a subhalo. The new version keeps tracking particles if the subhalo epoch lasts for only one snapshot because such fast flyby events do not necessarily mean that the particle trajectories are interrupted. Similarly, we do not categorically prohibit computing ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ for subhalos any more, as long as there are sufficient particle splashback events near the time in question (Paper I). We do, however, still abort all particle trajectories when a halo becomes a subhalo for more than one snapshot. Subsequently, the subhalo will quickly run out of past splashback events, and a determination of ${R}_{\mathrm{sp}}$ will not be possible any longer. This behavior is physically sensible because the splashback radius is ill-defined for subhalos. We have also updated our algorithm to determine which particles belong to a subhalo at infall (Section 3.4). Altogether, the algorithmic changes shift the average ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ by a few percent. To incorporate these differences, we recalibrate the model of Paper II in Section 5.6.

SO radii and masses are routinely calculated by halo finders such as Rockstar, but a number of subtleties can make it difficult to obtain exactly the desired definition. First, many halo finders offer only a limited number of definitions or particularly common ones such as ${R}_{200{\rm{c}}}$ or ${R}_{\mathrm{vir}}$. Second, we can include all particles (strict SO) or only gravitationally bound particles. The resulting masses tend to be close for host halos, but there are some exceptions. For subhalos, strict SO masses are often ill-defined, as they can include the entire host halo in some cases. To facilitate a thorough exploration of mass definitions, Sparta allows the user to compute any number of SO radii and masses with arbitrary density thresholds.

Sparta computes SO masses by considering a sorted list of particle radii within a user-defined factor times ${R}_{200{\rm{m}}}$. This factor generally needs to be somewhat larger than unity because the common ${{\rm{\Delta }}}_{\mathrm{vir}}$ definition tends to overdensities of about 178 at high redshift, rendering ${R}_{\mathrm{vir}}$ larger than ${R}_{200{\rm{m}}}$. For subhalos, we start at $4{R}_{200{\rm{m}},\mathrm{tcr}}$, where the halo radius is computed from tracked particles only (Section 3.4 and B. Diemer & P. S. Behroozi 2020, in preparation). Given a list of particles inside the search radius, sorted by radius, we start from the outermost particle and move inward until a density threshold is achieved for the first time. In some rare cases, the threshold can be crossed multiple times, and we wish to obtain the outermost occurrence. We set ${M}_{{\rm{\Delta }}}=i\times {m}_{{\rm{p}}}$, where i is the index of the particle where the threshold was crossed. We note that the internally used ${R}_{200{\rm{m}}}$ definition in Sparta interpolates the density between particles, but this difference is negligible in all but the smallest halos.

The calculation of an SO radius can fail for two reasons. First, the density may not reach the required threshold at any radius because even the very core of the halo is not dense enough. This issue occurs mostly for high-density definitions such as ${R}_{500{\rm{c}}}$ and can arise even in hosts halos. Second, the density may never fall below the threshold within the available particle distribution or within any reasonable radius. This issue cannot occur for hosts in Sparta because halos are defined via their ${R}_{200{\rm{m}}}$ internally. For subhalos, however, the entire host halo may be included if the subhalo is close enough to its center. This issue predominantly affects low-density definitions such as ${R}_{\mathrm{vir}}$ and ${R}_{200{\rm{m}}}$. If the calculation fails, we record an error code in the catalogs. We investigate the completeness of our SO measurements in Section 5.2.

In addition to strict SO masses, we also wish to include bound-only masses in our catalogs, particularly for subhalos. Our base catalogs from Rockstar already contain bound-only masses that are computed as follows. Rockstar considers the particles in a friends-of-friends (FOF) group (or subgroup for subhalos) and unbinds particles if their kinetic energy exceeds their potential energy. The process is not repeated iteratively. If a halo contains more than some fraction (half, in our case) of unbound particles, it is abandoned as a transient feature and not included in the catalog.

While this unbinding procedure is sensible in theory and gives excellent results in practice, all unbinding techniques turn out to be highly subjective and code dependent because they strongly depend on the initial particle distribution considered. For example, if we wanted to be conservative and be sure to include all bound halo particles, we could start from all particles inside $2{R}_{200{\rm{m}}}$—but with this much mass in the overall distribution, virtually all particles within ${R}_{200{\rm{m}}}$ would be bound. Conversely, if we considered only particles within ${R}_{500{\rm{c}}}$, very large fractions would be unbound. Since no SO definition is special, there is no "correct" distribution of initial particles. FOF-based halo finders such as Rockstar or SubFind use FOF groups as the initial particle set, which is sensible but depends on the linking length and the algorithm for discerning subclumps (e.g., Springel et al. 2001; Han et al. 2012; Behroozi et al. 2013a). This issue is exacerbated in subhalos, where the position of the subhalo within the host will influence its overall density and thus how many particles would be considered bound. After much experimentation, we provide a fairly general unbinding algorithm in Sparta. We check the inequality

$$\begin{eqnarray}&&-{{\rm{\Phi }}}_{{\rm{i}}}=G\displaystyle \sum _{j\ne i}^{{N}_{\mathrm{ptl}}}\displaystyle \frac{1}{{r}_{\mathrm{ij}}+\epsilon }\geqslant {f}_{\mathrm{bnd}}\displaystyle \frac{{v}_{{\rm{i}}}^{2}}{2},\end{eqnarray} \tag{ 1 }$$

where ${N}_{\mathrm{ptl}}$ is the number of particles to consider, f_bnd is a threshold for the binding-to-kinetic energy ratio, and is the force softening scale of the simulation. To compute the sum over many particles efficiently, we have adapted the tree algorithm of Rockstar for Sparta (see B. Diemer & P. S. Behroozi 2020, in preparation, for details). The user can set the initial radius for particles, the boundness factor f_bnd in Equation (1), and whether the procedure is iterated until it converges. We set the initial radius for host halos to ${R}_{200{\rm{m}}}$ and for subhalos to the tracer-only radius (Section 3.4), which eliminates the issue of unphysical all-particle radii in subhalos; however, it does not eliminate the issue of host material contributing to the potential.

We compare our bound-only masses to those from Rockstar in Section 5.4. Broadly speaking, we find that our results agree reasonably well for the majority of subhalos, although with strong outliers where the host contributes a significant amount of material. These tests demonstrate that the potential calculation in Sparta works as expected, but, given the issue of host contamination, it is not clear how useful the bound-only masses computed by Sparta would be. We have thus left them out of the catalogs and instead included the Rockstar bound-only masses for four common overdensity thresholds (Section 5.1).

The algorithms described in the previous sections apply mostly to host halos, where we can meaningfully compute splashback and SO radii and masses without appealing to FOF algorithms or other ways to decide on the membership of particles. We have, however, also introduced entirely new algorithms for dealing with subhalos in Sparta. These methods will be described in detail in B. Diemer & P. S. Behroozi (2020, in preparation), but we briefly review them here because the results are included in our catalogs.

Unlike the Sparta version of Paper I, the code now tracks all particles in subhalos. When the halo first becomes a subhalo (when it crosses ${R}_{200{\rm{m}}}$), we identify the particles that truly belong to the subhalo using an updated algorithm. Specifically, we apply user-defined thresholds to exclude particles that first joined the subhalo close to the host (as they were probably swept up during infall) and that are not strongly gravitationally bound. Thereafter, we assume that subhalos only ever lose particles and never accrete new ones. When particles stray more than $2{R}_{200{\rm{m}},\mathrm{tcr}}$ from the subhalo, they are removed. The "tracer masses," including ${R}_{200{\rm{m}},\mathrm{tcr}}$, are computed like normal SO masses but only including the remaining distribution of tracked particles. This definition gives similar but distinct results to Rockstar's unbinding algorithm (B. Diemer & P. S. Behroozi 2020, in preparation).

We also use the tracking algorithm to extend the lives of subhalos after they are lost by the halo finder. When a subhalo disappears from the input catalogs, we keep following its particles, at which point we term the halo a "ghost." We track the ghost's particles until there are only 10 particles left or until it has sunk to the center of the host halo. We compute the position and velocity of the ghost based on its most bound particles. Conceptually, ghosts are similar to "orphans" and "cores," where subhalos are represented as a single particle or a number of the most bound particles (e.g., Wang et al. 2006; Guo et al. 2010; Heitmann et al. 2019, 2020). Tracking all particles, however, allows us to compute physically meaningful tracer masses for ghosts. In B. Diemer & P. S. Behroozi (2020, in preparation), we show that tracking ghosts significantly increases the completeness of our merger trees.

Moria works in postprocessing, meaning that it can be run multiple times on the same Sparta output (Figure 2). The main reason for this choice is that the CPU time consumed by Moria is orders of magnitude smaller than that of Sparta because the former processes only halo data but not particles. For example, it would be highly inconvenient to rerun Sparta only to create a new halo catalog with a few extra fields from the original catalogs. Nevertheless, Moria performs I/O operations on large catalog and hdf5 files and uses tree searches to find host–subhalo relations, making it too computationally demanding for scripting languages. Thus, Moria is, like Sparta, written in pure C, but it runs on a single process.

The inputs to Moria are a Sparta output file and the original catalogs used by Sparta to create that file. The output format is flexible: the user can choose to write catalogs in Moria's native hdf5 format, the original catalog format (e.g., ASCII files for Rockstar), or a merger tree format that is essentially the same hdf5 format as the catalogs. The content of those files is entirely up to the user, who can choose any number of SO and splashback mass definitions that are present in the Sparta output file, as well as any number of fields from the original catalog. Mass definitions that exist in both the halo finder and Sparta outputs will be named accordingly, allowing for comparisons between different algorithms. Finally, the user can choose to compute an arbitrary number of host–subhalo relations based on any of the available radius definitions (Section 4.4). Virtually all parameters mentioned in the following sections can be changed by the user at run-time.

In Paper I, we showed that the Sparta algorithm for determining ${R}_{\mathrm{sp}}$ is at least about 95% complete for host halos with ${N}_{200{\rm{m}}}\geqslant 1000$ (and often better depending on halo mass and redshift). This completeness is somewhat unsatisfactory for a number of reasons. First, we would like to push below 1000 particles to include the large fraction of halos that are smaller than this limit. Second, even 5% of halos without splashback values can lead to inconvenient systematics, e.g., in mass functions. Third, many of the interrupted splashback histories occur because the halo temporarily becomes a subhalo, meaning that it enters within ${R}_{200{\rm{m}},\mathrm{bnd}}$ of a larger halo. When we compute host–subhalo relations for the splashback definitions, some of those subhalos may become hosts but not have ${R}_{\mathrm{sp}}$ measurements. This case occurs relatively rarely because the splashback radius is typically larger than ${R}_{200{\rm{m}}}$, but, depending on the exact definition, there is a significant fraction of halos for which it is slightly smaller (Paper II). We have partially addressed this issue in Section 3.1, but we would ideally like to achieve a completeness of 100% for host halos.

Moria addresses this issue by interpolating ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ across time. In many cases, splashback histories are interrupted by flyby events with more than one snapshot as a subhalo or a temporary lack of particles at first apocenter. Such gaps can be bridged because the overall splashback radius and mass histories tend to be relatively smooth, particularly when expressed as ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$. As a first guess, Moria looks for an earlier and later epoch within one dynamical time where the halo was a host and where ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ were successfully measured by Sparta. Furthermore, we check that ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$ (for all splashback definitions) are not too dissimilar from our expectation for a halo of the given mass and mass accretion rate, as quantified by the model from Paper II (see Section 5.6 and Equation (3)). We require the value to be within a factor of two or a logarithmic factor of $5\sigma$, where the standard deviation is quantified by Equation (6). These limits may seem overly generous, but the nominal standard deviation from the model can be as small as 0.02 dex or 5%. This deviation does not capture the significant tails in the distribution, which we do want to allow to avoid biasing our results. If at least one valid past or future epoch is found, we assume that ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$ vary slowly over a dynamical time and accept their values for the missing epoch. If both a past and future epoch are found, we linearly interpolate ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$ in time. Another potential source of error are the values of ${R}_{200{\rm{m}}}$ and ${M}_{200{\rm{m}}}$ used to convert the ratios to ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$. If ${R}_{200{\rm{m}}}$ includes a large amount of material from a nearby halo (Section 3.2), keeping the ratio fixed can lead to extreme values of ${R}_{\mathrm{sp}}$. Thus, we check the ratio between ${R}_{200{\rm{m}},\mathrm{all}}$ and ${R}_{200{\rm{m}},\mathrm{bnd}}$ and use the latter if the ratio exceeds two. This choice is consistent with the fitting function on Equation (3), where we exclude halos with extreme all-to-bound ratios.

The accuracy of the interpolation depends on a number of factors, including the time spacing of snapshots, how many snapshots lack ${R}_{\mathrm{sp}}$ determinations, and the lowest particle number considered. We have tested the interpolation from past snapshots at z = 0 in TestSim100. For halos with ${N}_{200{\rm{m}}}\gt 200$, we find that the guess based on previous snapshots is biased low by about 3% in radius and 5% in mass on average, with about 15% scatter. The low estimate is expected as ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ typically increases with time. These results are likely to be a worst-case scenario because ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ increases rapidly at z = 0 as mass accretion rates are falling quickly; this is particularly relevant in a small box such as TestSim100, where the majority of halos are low-mass galaxy halos. We find similar values for guesses based on future snapshots and much better results for interpolated values, depending on the time interval over which we interpolate. For short intervals of one or a few snapshots, the typical error is less than 1%. We have also experimented with extrapolating the evolution of ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ from past snapshots but found that, while it can slightly reduce the bias, it leads to somewhat erratic results and large scatter.

If no epoch with valid ${R}_{\mathrm{sp}}$ determination is found within a dynamical time, we use the mass and mass accretion rate of the halo to compute the fitting model of Section 5.6. This procedure is not circular because the model is calibrated using only halos where ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ were computed by Sparta. The fitting function is problematic in that we are imposing a median relation with significant scatter, but this guess is preferable over having no estimate at all (which could lead to significant differences in the subhalo statistics; see Section 4.4). The catalogs contain a status value that indicates how the splashback data for each halo and epoch were computed, allowing the user to easily exclude interpolated or predicted values. We test the completeness of our catalogs in Section 5.2.

Over the past few years, it has become clear that the mass accretion rate of halos is an important property that determines their density profiles, splashback radius, and certain baryonic properties (Diemer & Kravtsov 2014; Lau et al. 2015; More et al. 2015; Green et al. 2020). However, the optimal definition of the accretion rate is not obvious: while we would ideally like to define an infinitesimal derivative of mass akin to theoretical models (Adhikari et al. 2014; Shi 2016), this quantity has no meaning in simulations where it is entirely dominated by shot noise and mergers.

In Moria, we keep the definition of Paper I, the logarithmic change in all-particle ${M}_{200{\rm{m}}}$ per scale factor over one dynamical time, or ${{\rm{\Gamma }}}_{\mathrm{dyn}}$. However, Xhakaj et al. (2019b) recently reported that our calculation of ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ differs from Rockstar's owing to three seemingly subtle differences: the mass definition (bound-only ${M}_{\mathrm{vir}}$ in Rockstar, all-particle ${M}_{200{\rm{m}}}$ in Moria), interpolating between two past snapshots (which is done in Rockstar but not in Moria), and the definition of the dynamical time (${t}_{\mathrm{dyn},\mathrm{vir}}$ in Rockstar, ${t}_{\mathrm{dyn},200{\rm{m}}}$ in Moria). The latter is by far the most important factor, highlighting how sensitive accretion rates can be to the time interval over which they are computed. For comparison, the Rockstar mass accretion rates are included in our catalogs (Section 5.1). We note that the Behroozi et al. (2013a) definition of dynamical time is a factor of two shorter than ours, meaning that their rate measured over $2{t}_{\mathrm{dyn}}$ corresponds most closely to our definition.

Once again, a complication arises for subhalos and backsplash halos that were subhalos one dynamical time ago. We could replace the unphysical all-particle masses of subhalos with bound-only masses from Rockstar, but this mixture would lead to an apples-to-oranges comparison and ignore the physical reality that there is no clear equivalency between the SO masses of hosts and subhalos. Moreover, the splashback radius is influenced by the total change of mass within the orbits of particles, which includes all matter regardless of whether it "truly" belongs to the halo. Thus, we use ${M}_{200{\rm{m}},\mathrm{all}}$ but adaptively reduce the interval over which we compute ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ to avoid subhalo epochs. We stop if the interval becomes shorter than $0.25{t}_{\mathrm{dyn}}$, which occurs if the halo was a subhalo until very recently. In this rare case, we use the median value of ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ as predicted by the fitting function of Section 5.6, ignoring the substantial scatter of (0.41–0.07ν) dex. Physically speaking, we have to accept that there are halos for which it is difficult to define a meaningful accretion rate. For example, the evolution of backsplash halos can be dominated by stripping due to the flyby event. For subhalos, there is no sensible way to define an accretion rate at all. For completeness, the catalogs give ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ at the epoch when the subhalo was accreted, which we compute as for host halos. The catalogs contain a flag that indicates how the accretion rate was computed.

At z ≈ 0, ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ is computed by the default procedure for about 99% of host halos, over a reduced interval for 1%, and from the fitting function for a negligible fraction. In the smallest box at z = 0, L0063-WMAP7, those numbers rise to 4% and 0.2%. For subhalos, between 1% and 17% have reduced-interval accretion rates (at their infall redshift), and between 1% and 4% have only the fitting function estimate. These numbers vary only slightly with redshift. Moreover, Mansfield & Avestruz (2020) find that mass accretion rates in the Erebos simulations are well converged even at 100 particles per halo. In conclusion, we are able to measure mass accretion rates for virtually all halos for which this quantity is physically sensible.

Our process of assigning host–subhalo relations is illustrated in Figure 3. We find all halos whose centers lie within another, larger halo (see García & Rozo 2019, for alternative percolation algorithms). Here we need to define what "larger" means, since we will be dealing with a range of mass definitions. This ordering definition is a user-defined parameter in Moria. Given that none of the mass definitions discussed so far can be computed and are meaningful for all halos, we follow Behroozi et al. (2013a) in using ${V}_{\max }$ (as measured by Rockstar) to sort the halos. Starting with the largest halo, we use a tree search algorithm to look for all halos within its radius in a given mass definition (respecting the periodic boundary conditions). We mark all such halos as subhalos of the current halo unless they are already a subhalo.

Zoom In Zoom Out Reset image size

Some subhalos are inside multiple hosts, such as halo 3 in Figure 3. In this case, we assign the largest parent's ID (i.e., 1 in our example; this logic is equivalent to the "upid" field in Consistent-Trees catalogs). Another choice arises for subhalos whose only host is also a subhalo (halo 5 in Figure 3). Moria lets the user decide which parent ID to use; in the catalogs presented in this paper we assign the direct parent's ID (4 in this example). We have tested our assignment against Consistent-Trees and find perfect agreement when using the same radius definition.

One subtle question is how to handle missing radius estimates, for instance, when SO definitions cannot be computed (Section 3.2). We set those radii to zero, meaning that those halos cannot have subhalos. This solution seems sensible regardless of whether the SO could not be computed because the the halo never reaches the required central density at the center or because it is so close to a host that it includes its mass. Both issues apply predominantly to small halos that would not contain many subhalos regardless. Finally, we note that not all parent halos are necessarily part of the catalogs, depending on the chosen resolution cut (Section 5.1).

The term "merger tree" can refer to a format of outputting time-ordered halo information or to the history of one halo and its progenitors; here, we mean the former. The Moria merger trees contain basically the same information as the collection of catalogs from all snapshots, but the information is ordered differently. A typical ordering is "depth-first," where the subsequent entries represent the most massive progenitor branch of one halo going back in time. Once that branch ends, the next progenitor (starting at some earlier time) is listed, and so on. In this way, one can construct a logical scheme where halos are listed sequentially as in the catalog entries, except that a redshift or snapshot number must be assigned to each entry (e.g., Springel et al. 2001; Behroozi et al. 2013b; Rodriguez-Gomez et al. 2015; Elahi et al. 2019; see Srisawat et al. 2013 for a code comparison).

The Moria merger tree format is designed in an entirely different fashion as illustrated in Figure 4. We store halos in a 2D array spanning all halo histories and the number of snapshots in the simulation (a similar system is used in the HBT+ halo finder of Han et al. 2018). This format has the obvious advantage that it is easy to extract a catalog-like data set (by selecting a time slice) or a history-like data set (by selecting a halo index). At first sight, the format seems wasteful since about one-third of the array is occupied by halos that did not yet or do no longer exist. These empty fields, however, occupy very little additional disk space owing to the hdf5 compression algorithm. The only difference between our catalogs and a time slice of the tree data is that the tree includes the entire histories of all halos that are included in the catalogs at any snapshot, meaning that, at a given redshift, the tree will include some halos that would not pass the threshold in the corresponding catalog (Section 4.6). If so, a flag in the tree data indicates that a halo is present only in the tree file. We provide Python code to load the catalog and tree data, which takes care of this cut automatically if desired.

Zoom In Zoom Out Reset image size

The halo histories in the trees correspond to progenitor–descendant relations that we take directly from the original Consistent-Trees catalogs. These relations uniquely determine the next epoch in a halo history (horizontal gray arrows in Figure 4) or the tree branch into which a halo merges at the end of its life if it does not survive until the end of the simulation (vertical arrows). Moria does not attempt to improve the progenitor–descendant relations except for ghosts (blue fields in Figure 4).

We are now free to sort the tree file along the halo axis in any convenient fashion. We split the halos into subtrees, where each subtree contains a host halo at the final snapshot, its subhalos, and all histories that merged into them at previous redshifts. The subtrees are sorted by the same quantity that is used to impose a mass cut on the catalogs and trees, the peak ${M}_{200{\rm{m}}}$ along a halo history in our case. Thus, the first "line" in a tree file is the halo with the largest ${M}_{200{\rm{m}},\mathrm{peak}}$ that survived until the final snapshot. This row in the 2D array will be followed by the halos that merged into the top halo at the second-to-last snapshot, again ordered by their ${M}_{200{\rm{m}},\mathrm{peak}}$. Thus, a halo history is not necessarily followed by its own progenitors, leading to upward-pointing gray arrows that cross multiple lines in Figure 4.

Moria makes it a priority to apply a well-defined lower-mass cutoff to the catalogs, with the twofold purpose of dissuading the user from accidentally using poorly resolved halos and saving disk space. The latter effect is substantial because the steep mass function of halos means that halos near the resolution limit account for a significant fraction of a catalog's size. The cutoff can be given in units of halo mass according to some definition, as a number of particles corresponding to that mass, ${V}_{\max }$, or ${V}_{\max }$ corresponding to a number of particles. In the latter case, we estimate the limiting value of ${V}_{\max }$ using the empirical formula of Klypin et al. (2011),

$$\begin{eqnarray}&&{V}_{\max }\approx 2.8\times {10}^{-2}{M}_{\mathrm{vir}}^{0.316}.\end{eqnarray} \tag{ 2 }$$

Generally, it is preferable to set limits in units of particle number since they directly connect to the mass resolution of each simulation. We discuss the mass limits chosen for our catalogs in Section 5.1.

Furthermore, the user can decide how to treat phantoms and ghosts. Phantoms are rare cases where Consistent-Trees infers the existence of a halo from past and future snapshots even though Rockstar's FOF algorithm did not detect it. This can happen to both host and subhalos, although the latter make up for the vast majority of phantoms. Sparta uses the tracked particles in phantoms to compute their positions and velocities, which will not exactly agree with those from the halo finder (Section 3.4). The user can choose to replace the Rockstar values with those from Sparta, but we have not used this feature in our catalogs because it causes the host–subhalo relations to slightly differ from Consistent-Trees even if the same mass definition is used. We do not include ghosts in our catalogs because they do not have any catalog-defined properties, by definition; they are, however, included in the merger trees.

Our catalogs and trees are based on Rockstar and Consistent-Trees catalogs run with ${R}_{200{\rm{m}},\mathrm{bnd}}$ as the main definition, with the "strict SO" setting switched off and the "bound-props" setting on, meaning that all properties, including SO masses, are calculated from bound particles only. While the distance units used in the Rockstar catalogs are comoving $\,{h}^{-1}\,\mathrm{Mpc}$ and $\,{h}^{-1}\,\mathrm{kpc}$, all halo radii are given in physical $\,{h}^{-1}\,\mathrm{kpc}$ in the Moria catalogs.

We have chosen a mass limit of ${N}_{200{\rm{m}},\mathrm{bnd},\mathrm{peak}}\geqslant 200$ in Moria, meaning that a halo is included in the catalogs if it had at least 200 particles within ${R}_{200{\rm{m}},\mathrm{bnd}}$ at some time in the past. The trees contain the entire history of a halo if it exceeds the limit at any snapshot. There is no one ideal limit for all purposes because different halo properties converge at different particle numbers (Mansfield & Avestruz 2020). SO masses tend to be relatively well converged down to a few hundred particles, even for subhalos (e.g., Leroy et al. 2020), although their mass is subject to significant statistical variance (Trenti et al. 2010; Benson 2017). For splashback radii, we would ideally enforce a higher limit but instead caution the user that splashback masses and radii below 1000 particles become less complete (Paper I and Section 5.2).

We have included a large number of mass definitions, accretion rates, and status fields, some of which we summarize in Table 3. We refer the reader to the code documentation for units and further details. Most importantly, we include SO masses and radii in the ${M}_{500{\rm{c}}}$, ${M}_{200{\rm{c}}}$, ${M}_{\mathrm{vir}}$, and ${M}_{200{\rm{m}}}$ definitions, each calculated with all particles (from Sparta), bound-only (from Rockstar), and from tracers only (for subhalos, from Sparta). We also include splashback radii and masses corresponding to the mean, median, 70%, 75%, 80%, 85%, and 90% of the particle splashback distribution. For each of these radius definitions, we provide host–subhalo relations, that is, an array of host IDs for each halo. We also include the original host–subhalo assignment from Consistent-Trees. We include a number of flags that indicate the status of the halo in Sparta, which definitions were successfully computed, and whether fields were modified or added by Moria. In addition to the fields listed in Table 3, we have added almost all fields available in the Rockstar/Consistent-Trees catalogs, including scale radius, velocity dispersion, angular momentum, spin parameters, axis ratios, position and velocity offsets, half-mass radii and times, and tidal forces (see Behroozi et al. 2013a, 2013b, although some fields have been added since those original publications).

Table 3.  Partial List of Fields Included in the Moria Catalogs and Merger Trees

Field	From	Explanation
Organizational and merger information:
id	C,S	Halo ID from the original catalogs (except ghosts, where the ID is determined by Sparta)
descendant_id	C,S	ID of the halo at the next snapshot, or of the halo it merged into if it is ending (tree only)
descendant_index	M	Tree index of descendant; same as current halo unless merging (tree only)
mask_alive	M	True if there is an alive halo at this index and snapshot, false if not (tree only)
mask_cut	M	True if halo is part of catalog as well as tree (tree only)
num_prog	C	Number of progenitors (not counting ghosts)
Halo status:
status_sparta	S	Indicates host/sub/ghost status, switches in host, becoming a subhalo, becoming host, etc.
phantom	C	Nonzero if halo was interpolated by Consistent-Trees
status_acc_rate	M	Status indicating whether accretion rate was computed or guessed
status_sparta_rsp	S	Status of ${R}_{\mathrm{sp}}$ analysis in Sparta at this snapshot (success, not enough particles, etc.)
status_moria_rsp	M	Status of ${R}_{\mathrm{sp}}$ analysis in Moria (taken from Sparta, interpolated, guessed from fitting function, etc.)
status_moria_hps_M<definition>	M	Status of each SO mass (success, not dense enough, never below threshold, etc.)
status_moria_hps_R<definition>	M	Status of each SO radius definition (probably the same as for the SO mass of the same definition)
Important times in halo history and merger information:
Mpeak_Scale	C	Scale where highest ${M}_{200{\rm{m}},\mathrm{bnd}}$ was reached during halo's history
Acc_Scale	C	Scale where subhalo last became a subhalo, if it ever did
First_Acc_Scale	C	Scale where this halo first became a subhalo, if ever
scale_of_last_MM	C	Scale where last merger with mass ratio greater than 0.3 occurred, if ever
Time_to_future_merger	C	Time (in Gyr) until the halo merges into a larger halo, if ever
Future_merger_MMP_ID	C	ID of most massive progenitor into which halo merges (−1 if it does not exist at this time)
Mass and radius definitions:
R200m_all_spa_internal	S	Internally used ${R}_{200{\rm{m}}}$ from Sparta (${R}_{200{\rm{m}},\mathrm{all}}$ for hosts, ${R}_{200{\rm{m}},\mathrm{tcr}}$ for subhalos)
M200m_all_spa_internal	S	Internally used ${M}_{200{\rm{m}}}$ from Sparta (corresponding to R200m_all_spa_internal)
nu200m_internal	M	Peak height, ${\nu }_{200{\rm{m}}}$, calculated from internal ${M}_{200{\rm{m}}}$ or bound mass if ${M}_{200{\rm{m}},\mathrm{all}}\gt 2{M}_{200{\rm{m}},\mathrm{bnd}}$
R<delta>_all_spa	S	SO radius computed from all particles (with Δ = 500c, 200c, vir, and 200m)
R<delta>_tcr_spa	S	SO radius computed from tracked particles for subhalos, same as all-particle radius for hosts
R<delta>_bnd_cat	R	SO radius computed from bound particles only
M<delta>_all_spa	S	SO mass computed from all particles
M<delta>_tcr_spa	S	SO mass computed from tracked particles for subhalos, same as all-particle mass for hosts
M<delta>_bnd_cat	R	SO mass computed from bound particles only
Rsp-apr-mn	S	${R}_{\mathrm{sp}}$ computed from mean of the particle splashback distribution
Rsp-apr-p<percentile>	S	${R}_{\mathrm{sp}}$ computed from percentiles of the particle splashback distribution (50, 70, 75, 80, 85, 90)
Msp-apr-mn	S	${M}_{\mathrm{sp}}$ computed from mean of the particle splashback distribution
Msp-apr-p<percentile>	S	Msp computed from percentiles of the particle splashback distribution (50, 70, 75, 80, 85, 90)
Vmax	R	Maximum circular velocity
Macc	C	${M}_{200{\rm{m}},\mathrm{bnd}}$ at accretion (for subhalos)
Vacc	C	${V}_{\max }$ at accretion (for subhalos)
M200m_peak_cat	C	Highest ${M}_{200{\rm{m}},\mathrm{bnd}}$ attained during halo's history
Vpeak	C	Highest ${V}_{\max }$ attained during halo's history
Vmax@Mpeak	C	${V}_{\max }$ at scale Mpeak_Scale
First_Acc_Mvir	C	${M}_{\mathrm{vir}}$ when first becoming subhalo (if ever)
First_Acc_Vmax	C	${V}_{\max }$ when first becoming subhalo (if ever)
Host–subhalo relations:
parent_id_cat	C	Parent ID according to original catalog; equivalent to upid (most massive host)
parent_id_R<delta>_all_spa	M	Parent ID for all-particle SO definitions
parent_id_R<delta>_bnd_cat	M	Parent ID for bound-only SO definitions
parent_id_R<delta>_tcr_spa	M	Parent ID for tracer SO definitions (experimental and very similar to all-particle assignment)
parent_id_Rsp-apr-mn	M	Parent ID for ${R}_{\mathrm{sp}}$ from mean of the particle splashback distribution
parent_id_Rsp-apr-p<percentile>	M	Parent ID for ${R}_{\mathrm{sp}}$ from percentiles of the particle splashback distribution (50, 70, 75, 80, 85, 90)
Accretion rates:
acc_rate_200m_dyn	M	Fiducial accretion rate ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ as computed by Moria, see status_acc_rate; dimensionless units
Acc_Rate_1*Tdyn	C	Acc. rate over $0.5{t}_{\mathrm{dyn}}$ as defined here but using ${t}_{\mathrm{dyn},\mathrm{vir}};$ in absolute units of $\,{M}_{\odot }/h/\mathrm{yr}$
Acc_Rate_2*Tdyn	C	Acc. rate in ${M}_{200{\rm{m}},\mathrm{bnd}}$ over ${t}_{\mathrm{dyn}}$ as defined here but using ${t}_{\mathrm{dyn},\mathrm{vir}}$
Acc_Rate_Inst	C	Acc. rate in ${M}_{200{\rm{m}},\mathrm{bnd}}$ over one snapshot (very noisy)
Acc_Rate_100Myr	C	Acc. rate in ${M}_{200{\rm{m}},\mathrm{bnd}}$ over 100 Myr (typically short compared to ${t}_{\mathrm{dyn}}$)
Acc_Rate_Mpeak	C	Acc. rate in ${M}_{200{\rm{m}},\mathrm{bnd},\mathrm{peak}}$ from current z to $z+0.5$
Acc_Log_Vmax_1*Tdyn	C	Difference in ${\mathrm{log}}_{10}{V}_{\max }$ over $0.5{t}_{\mathrm{dyn},\mathrm{vir}}$
Acc_Log_Vmax_Inst	C	Difference in ${\mathrm{log}}_{10}{V}_{\max }$ over one snapshot
Log_(VmaxVmax_max(Tdyn;Tmpeak))	C	${\mathrm{log}}_{10}{V}_{\max }$ over ${V}_{\max }(t-{t}_{\mathrm{dyn}})$ or ${V}_{\max }({t}_{\mathrm{peak}})$ (the latter if ${M}_{\mathrm{peak}}$ happened more than ${t}_{\mathrm{dyn}}$ ago)
Other halo properties:
x	R,S	Halo position from Rockstar, except for ghosts (and optionally phantoms) where computed by Sparta
v	R,S	Halo velocity from Rockstar, except for ghosts (and optionally phantoms) where computed by Sparta

Note. The letters in the "From" column indicate the possible sources of a field, namely, R (Rockstar), C (Consistent-Trees), S (Sparta), or M (Moria). This list is incomplete; see the Sparta documentation for additional fields and units. In identifiers, "delta" can be replaced by "200m," "vir," "200c," or "500c." The "percentile" for ${R}_{\mathrm{sp}}$ definitions can be 50, 70, 75, 80, 85, or 90.

Download table as:  ASCIITypeset images: 1 2

We have set the compression level of the hdf5 files to the lowest level, 1, which gives a good compromise between file size and speed. The files are readable with any up-to-date hdf5 distribution. In total, our catalogs for the 14 Erebos simulations contain about 7 million halo histories and a total of 342 million halo epochs. On average, 64% of the tree arrays contain halos that are alive; the remaining 36% are filled with zeros. Out of the active halo epochs, about 88% are occupied by hosts, 8.6% by subhalos, and 3.8% by ghosts. These fractions vary from simulation to simulation, with smaller box sizes containing higher fractions of subhalos and ghosts, as well as an overall greater number of halos due to the higher variance of the density field on smaller scales.

As discussed in Sections 3.1–3.3, the calculation of SO and splashback radii can fail in certain situations. Figure 5 shows the corresponding incompleteness of our catalogs for a variety of simulations and redshifts. The top four panels address SO definitions, where the completeness varies strongly between all-particle and bound-only definitions and between hosts and subhalos (different line styles in Figure 5).

Zoom In Zoom Out Reset image size

The calculation of bound-only radii can fail if the halo never reaches the threshold density at its center. By construction, there are no halos for which the halo finder fails to calculate the main definition, ${R}_{200{\rm{m}},\mathrm{bnd}}$. As the threshold increases to 500c, an increasing number of subhalos fail to reach the threshold, but the number remains well below a percent (dotted–dashed lines). The same issue can occur even for host halos, although with a frequency of less than 10⁻³ (dashed lines). The results for bound-only definitions in host halos are similar to those for all particles, which is not surprising since virtually all particles near the halo center are bound.

We also attempt to determine all-particle SO masses for subhalos even though they are often intrinsically ill-defined because the density never falls below the threshold (unless the entire host halo is included). The resulting incompleteness can reach up to 50% for ${R}_{200{\rm{m}}}$ but decreases significantly toward ${R}_{500{\rm{c}}}$, which can be computed for the vast majority of subhalos (solid lines in Figure 5). The incompleteness in all-particle subhalo masses is a function of mass rather than particle number, leading to the large differences between simulations. This is not the case for bound-only definitions, host halos, or splashback definitions though. Either way, the large incompleteness of all-particle subhalo masses is not important because they are given for comparison rather than as a physically meaningful mass definition.

In some very rare cases, a host halo is close enough to another halo that Sparta cannot find a radius where the density profile decreases below $200{\rho }_{{\rm{m}}}$. We observe this case in the 10⁻⁴ of host halos in L0063-WMAP7 (dashed line in top left panel of Figure 5). We emphasize that this incompleteness captures only cases where ${R}_{200{\rm{m}},\mathrm{all}}$ could not be computed; there are many more cases where it is computed but includes significant material from neighboring halos (Section 5.3).

We now turn to the splashback radius (bottom row of Figure 5). The issue of completeness was discussed at length in Paper I, with the conclusion that Sparta produces splashback data for about 95% of halos with ${N}_{200{\rm{m}}}\geqslant 1000$. Here, the situation has changed because we have included smaller halos down to 200 particles, selected them by peak mass rather than current mass, and computed a value for ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ regardless of whether they could be computed by Sparta. As discussed in Section 4.2, we either interpolate the splashback radius from past and/or future values (solid lines) or estimate it using the fitting function (dashed lines). For most applications, the model estimates should not be used because the scatter in the true ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ ratio is large (Paper II). At z ≈ 0, we find that model estimates account for just under 10% of halos with N_peak = 1000 in the smallest boxes and much less than a percent in the larger boxes. This trend is driven by the larger number of mergers and close flybys in the smaller boxes, which lead to interrupted splashback histories (Diemer 2020a). The second and third panels show redshifts 2 and 4 (including a pink line for L0031, which runs only up to z = 2). By z = 4, the fraction of estimated values falls below a percent in all simulations. Finally, we consider the self-similar simulations in the bottom right panel. At N = 1000, the modeled fraction varies between 10% for n = −1 and around 1% for the other simulations.

In summary, Figure 5 confirms that pushing far below the canonical value of 1000 particles per halo leads to possibly problematic fractions of interpolated or estimated splashback values. We find very similar results for the Planck cosmology at the same box sizes. For SO definitions, we do not find a large evolution in the fractions with redshift, meaning that the general trends shown in Figure 5 persist. Given that we will largely use bound-only and tracer masses for subhalos, we can think of the SO definitions as complete. The splashback definitions are nearly complete for hosts as long as we use reasonably well resolved halos. One exception are the last few snapshots of a simulation: due to the missing future particle splashback events, the completeness and reliability of Sparta's splashback results decrease for the last $\approx 0.2{t}_{\mathrm{dyn}}$ (Paper I). Thus, it is generally advisable to consider the splashback results at z ≈ 0.1 rather than at z = 0. In the following sections, we combine data from all simulations of a given cosmology.

We now consider the ratios between different mass and radius definitions quantitatively, starting with SO definitions. Both all-particle and bound-only definitions are commonly used, for instance, in studies of the halo mass function and merger trees, respectively. While it is widely known that unbinding and the density threshold can have a significant effect, these differences are rarely quantified.

In Figure 6, we show histograms of the ratio between all-particle and bound-only SO definitions for halos from all simulations in our WMAP7 sample for which both were measured and for which the bound-only radius contains at least 200 particles. We use the ${R}_{200{\rm{m}},\mathrm{bnd}}$ definition to split hosts and subhalos. Furthermore, we crudely split the sample by peak height into low-mass (0.5 < ν < 1) and higher-mass halos (1.5 < ν < 2). Similarly, we plot the ratio at z = 0 and z = 2 to highlight the general redshift trends. When referring to SO definitions, we use radii and masses interchangeably as they are uniquely coupled via M ∝ R³. The logarithmic mass ratios translate into three times smaller logarithmic radius ratios.

Zoom In Zoom Out Reset image size

For host halos (top row of Figure 6), we notice a very small fraction where ${R}_{\mathrm{all}}\lt {R}_{\mathrm{bnd}}$ owing to numerical effects; in general, the bound-only mass is always smaller. We have checked that all-particle masses measured by Sparta and Rockstar agree perfectly (if Rockstar's strict SO option is activated), except in some rare cases where the FOF groups used by Rockstar are slightly incomplete at large radii. However, in TestSim100 we find that only about 1% of halos exhibit more than 1% difference in ${R}_{200{\rm{m}}}$. Returning to the all-bound comparison, Figure 6 shows that unbinding plays almost no role for the vast majority of host halos. The median ratios are unity in all definitions, and the mean ratios deviate from unity by no more than 3%. The strength of the tail depends on the density threshold, with ${R}_{500{\rm{c}}}$ leading to fewer large ratios. For example, at z = 0, the all-bound ratio for ${M}_{500{\rm{c}}}$ is higher than 10% for only 0.5% of halos but for about 9% when using ${M}_{200{\rm{m}}}$. These differences are caused by halos that are close to another, larger halo whose material is partially included in the all-particle mass. High density thresholds are useful in avoiding this proximity issue, but they also lead to more failures of the calculation.

While it is generally safe to use all-particle and bound-only masses interchangeably for hosts, there are noticeable differences for the SO radii in the lower-threshold definitions. For example, there is a standard deviation of 14% in ${R}_{\mathrm{vir}}$ in the z = 2 low-mass example. The corresponding tail toward large radius ratios highlights the importance of unbinding in the context of subhalo relations: even few halos with spuriously large radii could, in principle, create a larger number of subhalos. We investigate this difference in Diemer (2020a) and find that subhalo relations based on all-particle definitions lead to up to 5% more subhalos at the low-mass end.

We now turn to subhalos (bottom panels of Figure 6). Here, the ratios are peaked near unity, but a significant fraction of halos have larger ratios. The distributions begin to fall off around a mass ratio of 100, but this cutoff is not necessarily physical because Sparta stops looking for a solution when the all-particle radius gets too large. On the one hand, the main conclusion from Figure 6 is that all-particle masses are ill-defined for subhalos and should simply not be used. The median mass ratios vary between 1.2 and 5, increasing with redshift and decreasing with mass and overdensity threshold. Considering radii, on the other hand, a surprisingly large fraction of subhalos do have well-defined all-particle radii: the median radius ratio varies between 1.1 and 1.8. The large standard deviations of up to unity reflect that the result of the calculation will strongly depend on how close to the host center the subhalo is located. In summary, the bound-to-all ratio appears to increase in cases where halos are close to neighbors (for hosts) or to their host center (for subhalos).

We now explicitly confirm this picture by considering the mass evolution of subhalos in the WMAP7 cosmology in Figure 7. We use our new merger tree format to pick out all epochs where a halo becomes a subhalo and to extract its mass evolution before that time (while it was a host) and afterward (while it is a subhalo, as opposed to a ghost or backsplash halo). We define the time of infall as the last snapshot before each halo became a subhalo according to the ${R}_{200{\rm{m}}}$ definition (vertical dashed line in Figure 7). We require halos to contain at least 200 particles at infall, but this limit does not qualitatively change the results. Similarly, the mass evolution does not change much if we limit the infall redshift or peak height; the figure includes infalls at any redshift. We omit epochs where a mass definition cannot be measured, meaning that the median all-particle masses would be even larger if we could include the most extreme cases.

Zoom In Zoom Out Reset image size

Figure 7 shows that, before infall, the all-particle and bound-only masses evolve more or less in unison. After infall, the median all-particle mass rapidly increases by a factor of up to 10. The large scatter highlights that the mass evolution strongly depends on the orbital parameters and on the internal structure of host and subhalo. Despite the scatter, we observe a dip in the median all-particle mass after one dynamical time, when many subhalos have reached their first apocenter in the less dense outskirts of their hosts. The corresponding dips from subsequent orbits are smoothed out by phase mixing. Focusing on the evolution around infall, the all-particle mass begins to deviate from the bound-only mass about half a dynamical time before infall. As expected, the increase over the bound-only mass is larger for low-threshold definitions such as ${M}_{200{\rm{m}}}$ because their larger radius includes more host material. While ${M}_{500{\rm{c}}}$ increases by only about 10% at infall, ${M}_{200{\rm{m}}}$ has increased about 40% on average. This mass increase is largely responsible for the tails in the host halo distribution in Figure 6.

As mentioned in Section 3.3, we found that bound-only masses depend strongly on the chosen algorithm, most notably on the initial set of particles. In Figure 8, we compare bound-only masses computed by Sparta to those from Rockstar. Here, we consider all particles within ${R}_{200{\rm{m}}}$, defined as the all-particle radius for hosts and the tracer radius for subhalos. We do not iteratively unbind to match Rockstar's procedure and because iterating does not significantly alter our conclusions. The results in Figure 8 refer to TestSim100 because we have not included Sparta's bound-only masses in our fiducial catalogs.

Zoom In Zoom Out Reset image size

As expected, the results agree very well for host halos, with median ratios of unity in all definitions and standard deviations between 2% in ${M}_{500{\rm{c}}}$ and 8% in ${M}_{200{\rm{m}}}$. This agreement reconfirms that we can safely neglect boundedness for host halos except in some extreme merger scenarios. For subhalos, the estimates can differ significantly, with a long tail toward large ratios. The shaded areas in Figure 8 show the 68% and 95% intervals for subhalos. While the logarithmic scale of the figure highlights the tails, the median ratios vary between only 1.1 for ${M}_{500{\rm{c}}}$ and 1.35 for ${M}_{200{\rm{m}}}$. For ${M}_{500{\rm{c}}}$, the estimates agree to 10% or better for about half the halos; for ${M}_{200{\rm{m}}}$ that fraction goes down to 30%. However, the tails and standard deviations are extremely large in all definitions. The reason for these disagreements are halos where ${R}_{200{\rm{m}}}$ includes numerous host particles, rendering the entire distribution more bound than it would be if we considered only "true" subhalo particles. The 6D-FOF groups used for subhalos in Rockstar do not constitute a "correct" answer either, but they contain less host material because the velocity of the particles is considered. We thus provide the Rockstar bound-only masses in our catalogs. We continue the search for physically motivated definitions of subhalo masses in B. Diemer & P. S. Behroozi (2020, in preparation), where we further explore the concept of tracer masses.

In Paper II, we investigated the splashback radii and masses computed by Sparta as a function of ${R}_{200{\rm{m}}}$ and ${M}_{200{\rm{m}}}$, accretion rate, redshift, and cosmology. However, we did not show the distribution of ratios and did not consider subhalos. Figure 9 shows ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$ but using the bound-only SO radii and masses from Rockstar to make the ratios more comparable to Figure 6 and to avoid the issue of spuriously large all-particle radii (Section 5.3). The distributions are clearly peaked around the median (shaded areas), but they do exhibit broad tails. Some of the most extreme ratios are caused by flyby encounters that can result in negative accretion rates and large splashback radii. The cutoff above a radius ratio of three is partly caused by Sparta's algorithm, which does not trace particles to arbitrary distances (Paper I). The distribution of ratios does not strongly depend on redshift, although the average splashback radius and mass slightly decrease with z (in agreement with Paper II). We choose ${R}_{\mathrm{sp},75 \% }$ as a representative definition in Figure 9; the distributions for the other percentiles look similar (although slightly shifted in radius and mass). One important difference between Paper II and our new catalogs is that Moria reconstructs splashback quantities by interpolation, extrapolation, and estimates from our fitting function. Those results are shown as dashed lines in Figure 9. The good agreement with the measured splashback quantities is somewhat by construction since we interpolate in ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$ space, but it demonstrates that the interpolation does not introduce any unexpected features.

Zoom In Zoom Out Reset image size

For subhalos, Sparta does not compute particle splashbacks because orbits reflect the combined potential of host and subhalo and do not typically exhibit well-defined apocenters. However, if a halo becomes a subhalo only briefly, there may be enough past and future particle splashbacks to compute ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$. The corresponding distributions are similar to those of host halos (solid light-blue lines in Figure 9). The estimated subhalo splashback radii (dashed light-blue lines) refer to the time of infall and are thus not directly comparable to the current-time radius and mass according to the halo finder; they are given mostly for completeness.

Based on the mass evolution of SO halos shown in Figure 7, we might worry about the evolution of splashback radii and masses close to merger events. For example, if the radii were artificially inflated as a result of particles being stripped from the halo, that could lead to other halos in the vicinity falsely being identified as subhalos. However, we find no evidence for such an issue. Figure 10 shows the average radius and mass evolution before infall for all (future) subhalos with more than 1000 particles in the WMAP7 cosmology. As in Figure 7, we observe a smooth mass evolution in bound-only ${M}_{200{\rm{m}}}$ (blue lines) and a sudden growth of all-particle ${M}_{200{\rm{m}}}$ shortly before infall (purple lines); ${R}_{200{\rm{m}},\mathrm{all}}$ follows this growth by construction. The splashback radii (orange and red lines), on the other hand, do not exhibit any sudden growth. The splashback mass is independent of the radius and does include material in the vicinity of the host, leading to a gradual increase in ${M}_{\mathrm{sp}}$ during the last dynamical time before the merger. The mean (dashed lines) grows much more sharply than the median, indicating a strongly asymmetric distribution with a tail toward large mass increases. We conclude that our splashback radii and masses behave as expected in the vicinity of mergers, at least on average.

Zoom In Zoom Out Reset image size

As discussed in Section 3.1, we have made improvements to the Sparta code that change the splashback radii and masses by a few percent. Thus, we recalibrate the Paper II model including the following changes. We calculate mass accretion rates using Moria as described in Section 4.3. Including halos that were subhalos or did not exist one dynamical time ago helps to extend the high-z samples, especially at high accretion rates. We have confirmed that excluding backsplash halos makes a small difference and thus leave them in the sample. Since we are constraining ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$, we should use only halos for which ${R}_{200{\rm{m}}}$ and ${M}_{200{\rm{m}}}$ are reliably measured. In particular, we wish to avoid halos that include a significant contribution from a neighboring halo, which is relatively common at low masses (Sections 5.2 and 5.3). We could sidestep this issue by using bound-only masses, but there are good reasons to calibrate the relationship against all-particle quantities. First, bound-only masses cannot be measured observationally and depend on the unbinding algorithm; second, particle orbits (and thus ${R}_{\mathrm{sp}}$) react to the total mass inside their radius rather than the bound mass; and third, all-particle masses are more clearly related to the halo density profile and the steepening feature. Thus, we use all-particle radii and masses but exclude halos for which ${M}_{200{\rm{m}},\mathrm{all}}$ is more than 20% greater than ${M}_{200{\rm{m}},\mathrm{bnd}}$. While this cut may seem strict, including halos with mass increases of up to 50% visibly increases the scatter in the lowest peak height bins, indicating that the splashback–SO relation is erratic in such halos. The higher peak height bins are barely affected by the cut. As in Paper II, we use only halos with ${N}_{200{\rm{m}}}\gt 1000$ and with successful ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ measurements from Sparta (as opposed to interpolated values).

With this halo sample in hand, we repeat the fitting procedure of Paper II. In particular, we fit the splashback results from both the WMAP7 and Planck cosmologies at redshifts 0.13, 0.3, 0.5, 1, 2, 4, and 8, avoiding z = 0 owing to the corrections at the final snapshots (Section 3.1). We bin ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ as a function of mass accretion rate ($0\lt {{\rm{\Gamma }}}_{\mathrm{dyn}}\lt 12$) and peak height (with bin edges at 0.5, 1, 1.5, 2, and 3). We use a least-squares fit to minimize the difference between the binned relations and our fit. We combine the statistical uncertainty in each bin (due to the finite number of halos) with an additional 1% systematic error in quadrature. This systematic error prevents high-occupation bins from dominating the fit and captures a number of inaccuracies in our methodology. The fitting function remains the same as in Equations (5)–(7) in Paper II. In particular, the relationships between ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ and the splashback quantities ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$ (which we summarily call ${X}_{\mathrm{sp}}$) are fit with the general function

$$\begin{eqnarray}&&{X}_{\mathrm{sp}}={A}_{0}+{{Be}}^{-{{\rm{\Gamma }}}_{\mathrm{dyn}}/C},\end{eqnarray} \tag{ 3 }$$

where A₀, B, and C are free parameters. The latter two are functions of peak height and ${{\rm{\Omega }}}_{{\rm{m}}}(z)$,

$$\begin{eqnarray}\begin{array}{rcl}B & = & ({B}_{0}+{B}_{{\rm{\Omega }}}{{\rm{\Omega }}}_{{\rm{m}}})\times (1+{B}_{\nu }\nu )\\ C & = & ({c}_{0}+{C}_{{\rm{\Omega }}}{{\rm{\Omega }}}_{{\rm{m}}}+{C}_{{\rm{\Omega }}2}{{\rm{\Omega }}}_{{\rm{m}}}^{2})\times (1+{C}_{\nu }\nu +{C}_{\nu 2}{\nu }^{2}).\end{array}\end{eqnarray} \tag{ 4 }$$

We separately fit ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$, as well as the mean-based definition and percentiles, leading to four sets of best-fit parameters. The dependence of the parameters on the percentile is parameterized in terms of p, the percentile divided by 100,

$$\begin{eqnarray}\begin{array}{rcl}{A}_{0} & = & {a}_{0}+{a}_{{\rm{p}}}\times p\\ {B}_{0} & = & {b}_{0}+{b}_{{\rm{p}}}\times p\\ {B}_{{\rm{\Omega }}} & = & {b}_{{\rm{\Omega }}}+{b}_{{\rm{\Omega }}{\rm{p}}}\times \exp ({b}_{{\rm{\Omega }}{\rm{p}}2}\times p)\\ {B}_{\nu } & = & {b}_{\nu }+{b}_{\nu {\rm{p}}}\times p\\ {C}_{{\rm{\Omega }}} & = & {c}_{{\rm{\Omega }}}+{c}_{{\rm{\Omega }}{\rm{p}}}\times \exp ({c}_{{\rm{\Omega }}{\rm{p}}2}\times p)\\ {C}_{{\rm{\Omega }}2} & = & {c}_{{\rm{\Omega }}2}+{c}_{{\rm{\Omega }}2{\rm{p}}}\times \exp ({c}_{{\rm{\Omega }}2{\rm{p}}2}\times p)\\ {C}_{\nu } & = & {c}_{\nu }+{c}_{\nu {\rm{p}}}\times p\\ {C}_{\nu 2} & = & {c}_{\nu 2}+{c}_{\nu 2{\rm{p}}}\times p.\end{array}\end{eqnarray} \tag{ 5 }$$

We also fit for the logarithmic scatter in the relation,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{sp}}={\sigma }_{0}+{\sigma }_{{\rm{\Gamma }}}{{\rm{\Gamma }}}_{\mathrm{dyn}}+{\sigma }_{\nu }\nu +{\sigma }_{{\rm{p}}}p.\end{eqnarray} \tag{ 6 }$$

The new best-fit parameters are given in Table 4; some parameters are not necessary in all four fits as indicated by zeros. The new model is implemented in the publicly available Colossus code. The fit quality is very similar to that reported in Paper II, namely, better than 5% essentially everywhere in the Γ−ν−${{\rm{\Omega }}}_{{\rm{m}}}$ parameter space where we have data. ${{\rm{\Delta }}}_{\mathrm{sp}}$ is derived from ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$, with a combined uncertainty of 15% or better. For a visual impression of these accuracies, we refer the reader to Figure 3 of Paper II. The model is valid for $0\lt {{\rm{\Gamma }}}_{\mathrm{dyn}}\lt 12$ and any percentile between 50 and 90 (higher percentiles become extremely noisy and were thus not computed by Sparta). We have checked the model against both the WMAP7 and Planck cosmologies and find similarly good agreement. The model does not describe self-similar simulations, although it comes close for slopes of n ≈ −2.5 that resemble ΛCDM (Paper II).

Table 4.  Best-fit Parameters for the Splashback–SO Model

Parameter	${R}_{\mathrm{sp},\mathrm{mn}}$	${R}_{\mathrm{sp}, \% }$	${M}_{\mathrm{sp},\mathrm{mn}}$	${M}_{\mathrm{sp}, \% }$
Parameters for ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$
a0	0.6597	0.3071	0.6962	0.2874
b0	0.5562	0.2508	0.3736	0.6616
${b}_{{\rm{\Omega }}}$	0.1141	0.1527	0.3005	0.1321
${b}_{\nu }$	0.0698	0.1956	0	0
c0	−0.8508	−1.2214	3.3445	4.5913
${c}_{{\rm{\Omega }}}$	18.4464	17.5374	1.3718	3.0928
${c}_{\nu }$	−0.3332	0	−0.0825	−0.1155
${c}_{{\rm{\Omega }}2}$	−10.0596	−10.3158	0	0
${c}_{\nu 2}$	0.0474	−0.0189	0	0
Meta-parameters for Dependence on Percentile
${a}_{{\rm{p}}}$	0	0.6428	0	0.8228
${b}_{{\rm{p}}}$	0	0.5074	0	−0.6567
${b}_{{\rm{\Omega }}{\rm{p}}}$	0	0	0	0.0032
${b}_{{\rm{\Omega }}{\rm{p}}2}$	0	0	0	4.9536
${b}_{\nu {\rm{p}}}$	0	−0.2128	0	0.2882
${c}_{{\rm{\Omega }}{\rm{p}}}$	0	0.0024	0	−0.6609
${c}_{{\rm{\Omega }}{\rm{p}}2}$	0	9.7115	0	−1.1050
${c}_{{\rm{\Omega }}2{\rm{p}}}$	0	−0.0005	0	0
${c}_{{\rm{\Omega }}2{\rm{p}}2}$	0	10.7626	0	0
${c}_{\nu {\rm{p}}}$	0	−0.4735	0	−0.3761
${c}_{\nu 2{\rm{p}}}$	0	0.0940	0	0.0784
Parameters for 68% Scatter (in dex)
${\sigma }_{0}$	0.0501	0.0419	0.0456	0.0224
${\sigma }_{{\rm{\Gamma }}}$	0.0035	0.0043	0.0017	0.0009
${\sigma }_{\nu }$	−0.0108	−0.0141	−0.0079	−0.0091
${\sigma }_{{\rm{p}}}$	0	0.0235	0	0.0454

Download table as:  ASCIITypeset image

Figure 11 shows a visual comparison of the old and new fitting functions. The differences in ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{sp}}/{M}_{200{\rm{m}}}$ are below 5% for all redshifts, masses, and cosmologies that we tested and below 2% for most of the parameter space. The differences in ${{\rm{\Delta }}}_{\mathrm{sp}}$ (which is derived from ${R}_{\mathrm{sp}}$ and ${M}_{\mathrm{sp}}$) multiply to a maximum of 15% but are typically around 5%. There is no systematic upward or downward trend in the results. These small differences confirm that, despite numerous changes and bug fixes, Sparta's algorithm to determine ${R}_{\mathrm{sp}}$ is robust.

Zoom In Zoom Out Reset image size

We have also updated our fit to the mass accretion rate as a function of peak height and redshift. This formula is convenient when computing splashback properties without any knowledge of the mass accretion rate,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{dyn}}=A\nu +B{\nu }^{3/2},\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}A & = & 1.1721+0.3255z\\ B & = & -0.2565+0.0932z-0.0571{z}^{2}+0.0042{z}^{3}.\end{array}\end{eqnarray} \tag{ 8 }$$

The differences to the fit from Paper II are minuscule, but the new fit is more consistent with the parameters in Table 4.

While we have extensively investigated the relationship between Sparta's dynamically determined ${R}_{\mathrm{sp}}$ and SO radii, we have not yet established a connection to the most commonly used definition of the splashback radius: the radius where the logarithmic slope of the spherically averaged density profile is steepest, R_steep (Figure 1). The splashback radius was first discovered as a result of the sharp drop in density (Diemer & Kravtsov 2014), first quantified according to this definition (More et al. 2015), and all observational detections have hitherto relied on the density profile (More et al. 2016, and following work). Moreover, in theoretical models, there is no ambiguity between R_steep and the radius where the most recently accreted particles reach their first apocenter (Adhikari et al. 2014; Shi 2016). In realistic halos, however, ${R}_{\mathrm{sp}}$ and ${R}_{\mathrm{steep}}$ are not equivalent because the apocenters are spread out over a significant radial range owing to nonsphericity, complex accretion histories, and interactions with other halos (Paper I). Similarly, the density drop represents a trade-off between the radially decreasing density of orbiting particles (or galaxies) and the stream of infalling material at large radii. As a result, it is not clear that there is a one-to-one correspondence between features measured in density profiles (or correlation functions) and the dynamically determined splashback radii from Sparta (Xhakaj et al. 2019a; Aung et al. 2020; Garcia et al. 2020).

We attempt to quantify the connection in Figure 12, where we show the percentiles that most closely match R_steep as approximated by the formula of More et al. (2015). This fitting function was derived from the averaged (median) density profiles of halos binned in accretion rate and peak height. To be comparable to observations, the profiles included subhalo particles. More et al. (2015) used the Diemer & Kravtsov (2014) definition of accretion rate, ${{\rm{\Gamma }}}_{\mathrm{DK}14}$, rather than ${{\rm{\Gamma }}}_{\mathrm{dyn}}$. There, the time interval was defined by fixed redshift intervals instead of the dynamical time. To evaluate the fitting function, we crudely convert ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ to ${{\rm{\Gamma }}}_{\mathrm{DK}14}$ based on linear fits to the median relation in the WMAP7 sample, which is well approximated as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{DK}14}={{\rm{\Gamma }}}_{0}+\gamma {{\rm{\Gamma }}}_{\mathrm{dyn}}.\end{eqnarray} \tag{ 9 }$$

At redshifts 0 and 0.5, the accretion rates are virtually identical because both methods choose the same redshift to compare the current masses to. For the higher redshifts defined in Diemer & Kravtsov (2014), z = [1, 2, 4], we find Γ₀ = [0.51, 1.43, 0.88] and γ = [0.74, 0.53, 0.78], respectively. We now obtain ${R}_{\mathrm{steep}}/{R}_{200{\rm{m}}}$ and ${M}_{\mathrm{steep}}/{M}_{200{\rm{m}}}$ from the More et al. (2015) fit at each z and ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ and find the percentile where the fitting function of Equation (3) most closely matches these values (if such a value exists between the 50th and 90th percentiles). We also evaluate the scatter using Equation (6) and find the range of percentiles where R_steep and M_steep lie within the 1σ scatter (shaded areas in Figure 12).

Zoom In Zoom Out Reset image size

Clearly, the relationship between dynamical ${R}_{\mathrm{sp}}$ and ${R}_{\mathrm{steep}}$ depends on accretion rate, mass, and redshift in a complex fashion. For example, Figure 12 explains the trend we noticed in the example of Figure 1, where higher ${{\rm{\Gamma }}}_{\mathrm{dyn}}$ drives ${R}_{\mathrm{steep}}$ toward a lower percentile. The mass dependence is driven solely by Sparta's measurements because the More et al. (2015) formula does not depend on mass, a conclusion that should be revisited with improved simulation data. The ranges where the scatter around the medians matches ${R}_{\mathrm{steep}}$ are large, which reflects the fact that the differences between the percentiles are comparable to the scatter in the relation. Nevertheless, the median relations are statistically well defined.

We conclude that there is no one percentile value that can be used to mimic ${R}_{\mathrm{steep}}$ and ${M}_{\mathrm{steep}}$. The exact conversion will depend not only on accretion rate, mass, and redshift but also on the exact fitting procedure used to derive R_steep and on whether the profile of dark matter particles or subhalos is considered (Xhakaj et al. 2019a). We leave a more detailed calibration of the conversion for future work. For practical purposes, percentiles such as those in Figure 12 can be obtained by numerically matching the fitting functions of More et al. (2015) and Section 5.6 for a given accretion rate, mass, and redshift, for example, using the implementations in Colossus.