key: cord-0634308-5k7xvkkr
authors: Diemer, Benedikt
title: Universal at last? The splashback mass function of dark matter halos
date: 2020-07-20
journal: nan
DOI: nan
sha: ec845e11ab861b1a0a0b143178bb5b44f801034a
doc_id: 634308
cord_uid: 5k7xvkkr

The mass function of dark matter halos is one of the most fundamental statistics in structure formation. Many theoretical models (such as Press-Schechter theory) are based on the notion that it could be universal, meaning independent of redshift and cosmology, when expressed in the appropriate variables. However, simulations exhibit persistent non-universalities in the mass functions of the virial mass and other commonly used spherical overdensity definitions. We systematically study the universality of mass functions over a wide range of mass definitions, for the first time including the recently proposed splashback mass, Msp. We confirm that, in LambdaCDM cosmologies, all mass definitions exhibit varying levels of non-universality that increase with peak height and reach between 20% and 500% at the highest masses we can test. Mvir, M200m, and Msp exhibit similar levels of non-universality. There are, however, two regimes where the splashback mass functions are significantly more universal. First, they are universal to 10% at z<2, whereas spherical overdensity definitions experience an evolution due to dark energy. Second, when additionally considering self-similar cosmologies with extreme power spectra, splashback mass functions are remarkably universal (to between 40% and 60%) whereas their spherical overdensity counterparts reach non-universalities between 180% and 450%. These results strongly support the notion that the splashback radius is a physically motivated definition of the halo boundary. We present a simple, universal fitting formula for splashback mass functions that accurately reproduces our simulation data.

In a universe where galaxies form at the center of dark matter halos, the number of halos as a function of their mass is one of the most basic statistics. This mass function serves as a crucial connecting point between cosmological models and observations, for instance when counts of clusters are translated into constraints on cosmological parameters (e.g., Allen et al. 2011; Kravtsov & Borgani 2012) . The seminal work of Press & Schechter (1974) presented a simple model where peaks in a Gaussian random field collapse if their overdensity passes a certain threshold (see also Bond et al. 1991) . This model assumes that the collapse threshold, and thus the mass function, are universal, meaning independent of redshift and cosmology. Here, the abundance of halos is expressed as the multiplicity function, f (σ), the fraction of matter that has collapsed into halos within a logarithmic interval in the variance of the underlying Gaussian density field. If true, this universality would allow us to find powerful functional forms for f (σ) that would describe the abundance of halos in all cosmologies, without the need to simulate each one separately.

While the specific form of f (σ) suggested by Press & Schechter (1974) was shown to deviate from simulation results (e.g., Gross et al. 1998; Tormen 1998) , the quest for a universal mass function continued. When halo masses are measured using a friends-of-friends (FOF) algorithm (Huchra & Geller 1982; Davis et al. 1985) , f (σ) has been found to be roughly independent of redshift (Jenkins et al. 2001; White 2002; Reed et al. 2003; Warren et al. 2006; Lukić et al. 2007; Heitmann et al. 2019 ). On closer inspection, however, the FOF mass function does exhibit deviations from universality at high redshifts and between cosmologies (Reed et al. 2007; Crocce et al. 2010; Bhattacharya et al. 2011) , with some dependence on the linking length parameter (More et al. 2011; Courtin et al. 2011 ). More fundamentally though, FOF masses are a theoretical construct that does not place the halo boundary at a particular density and are thus particularly difficult to connect to observations. Spherical overdensity (SO) masses are preferable for this purpose, but the relation between FOF and SO masses depends on redshift and cosmology and exhibits large, non-Gaussian scatter (White 2001; Tinker et al. 2008; Lukić et al. 2009; More et al. 2011; Courtin et al. 2011) . For SO masses, the question of universality depends on the overdensity threshold, a fixed or varying multiple of the critical or mean density of the Universe that results in definitions such as M 200c , M vir , and M 200m (Lacey & Cole 1994; White 2001) . The SO mass function has been shown to vary with redshift and cosmology (Tinker et al. 2008; Watson et al. 2013; Bocquet et al. 2016 Bocquet et al. , 2020 , although some claims of universality with redshift have been made for particular definitions (e.g., Despali et al. 2016) . At this juncture, one possible way forward is to predict the mass function with increasingly complex models that may eventually capture some nonuniversality, for instance based on excursion sets (e.g., Maggiore & Riotto 2010; Musso & Sheth 2012) or on peak statistics (e.g. Bond & Myers 1996; Monaco et al. 2002) In this work, we ask a different question: could the nonuniversality be an artifact of the SO mass definition? Should we expect universality when using such definitions, and are there alternatives? To this end, we undertake a first exploration of the splashback mass function of halos. The splashback radius, R sp , and its enclosed mass, M sp , have been proposed as a more physically motivated halo definition (Diemer & Kravtsov 2014; Adhikari et al. 2014; More et al. 2015) . The splashback radius is located at the first apocenter radius of the most recently infalling dark matter, which is marked by a sharp drop in the density profile. R sp depends on the mass accretion rate of a halo and should, by analogy with the spherical collapse model, separate infalling from orbit-2 Splashback Mass Function ing material (Fillmore & Goldreich 1984; Bertschinger 1985; Shi 2016) . Given this physical picture, we might hope that the splashback mass corresponds more tightly to the universal collapse threshold assumed in many mass function models. On the other hand, measuring the splashback radii of individual simulated halos is more complicated than for SO definitions, which could introduce artificial non-universalities. To test the universality of the mass function, we need measurements of SO and splashback masses over a wide range of halo mass, redshift, and cosmology, including universes with radically different power spectra. We have recently presented such a dataset, namely publicly available halo catalogs and merger trees with splashback data computed using the Sparta code (Diemer 2017; Diemer et al. 2017; Diemer 2020a, hereafter Paper I, Paper II, and Paper III, respectively) .

This paper is organized as follows. We discuss the theoretical underpinnings of mass functions in Section 2 and describe our simulation data and algorithms in Section 3. We analyze the universality of the mass function in Section 4. We further discuss the notion of universality in Section 5 before summarizing our results in Section 6.

2. THEORY In this section, we review the universal form of the mass function and the necessary numerical calculations, which can be reproduced with the publicly available code Colossus (Diemer 2018 ).

The mass function, by definition, represents the comoving number density of halos as a function of mass, n(M), in a given cosmology, at a given redshift z, and (in practice) in a simulation with a finite box size L, particle number N 3 , and other numerical limitations (Table 1) . Here, M is understood to represent the mass according to some definition. This basic count n(M), however, is not directly comparable between cosmic times or cosmologies because masses grow over time and because their number depends on cosmological parameters. Instead, we convert to the multiplicity function,

which is natural in the Press & Schechter (1974) picture (see e.g. Kravtsov & Borgani 2012 , for a derivation). Here, ρ m,0 is the mean matter density of the Universe at z = 0 and σ = σ(R L , z) is the variance of the linear density field on the Lagrangian scale of a halo. While mass has commonly been expressed as σ −1 , we choose to convert it to peak height, which compares σ to the overdensity at which we expect halos to collapse,

where δ c (z) is the critical overdensity for collapse and D + (z) is the linear growth factor. We also convert f (σ) to peak height,

where we have used that d ln(σ −1 ) = d ln(ν) at fixed redshift. The multiplicity function can now be interpreted as the fraction of mass that has collapsed into halos in a logarithmic unit interval in peak height. To investigate the universality of f (ν), we need to carefully consider the terms in Equation 2. The first ingredient is the variance of the linear power spectrum, P(k),

We discuss the integration limits in Section 2.2. We use a Boltzmann code to compute the power spectra, which are also used to generate the initial conditions of our simulations (Section 3.1). Colossus interpolates the tabulated spectra with a cubic spline (Diemer 2018) . One important ingredient in Equation 4 is W, the Fourier transform of the top-hat filter in real space,

We discuss alternative estimates of P(k) and other filter choices in Section 5.1. The Lagrangian scale of a halo with mass M is the comoving radius of a sphere that encompasses the halo's mass at the mean density of the Universe,

The final ingredient is the critical overdensity for collapse, which is often set to a constant, δ c,EdS 1.68647, derived from the spherical top-hat collapse model in an Einstein-de Sitter universe (Gunn & Gott 1972) . If δ c (z) ≡ δ c,EdS , peak height and 1/σ would be related by a constant factor. However, due to dark energy, the collapse density evolves slightly with redshift,

This expression (from Mo et al. 2010 ) is almost identical to that of Kitayama & Suto (1996) . For the cosmologies considered in this paper, δ c evolves by less than 1%, but even small changes in ν can lead to large changes in the universality of f (ν). We discuss the impact of this evolution in Section 5.1.

The σ integral in Equation 4 should ideally be taken from zero to infinity because all k-scales contribute to the top-hat filter. In simulations, the mass function is affected by their finite box size, which cuts off the power spectrum at a scale k min = 2π/L (Tormen & Bertschinger 1996; Yoshida et al. 2003; Bagla & Ray 2005; Power & Knebe 2006; Reed et al. 2007; Lukić et al. 2007; Gnedin et al. 2011) . We note that k min is not taken into account when normalizing the power spectrum to a given σ 8 because that correct value was used when setting up the initial conditions. The reduction in variance is due to power spectrum modes that were cut off in the initial conditions and is not recovered through renormalization. With k min > 0, the f (ν) form of the multiplicity function has the advantage that it avoids converting n(M) to f (σ) via the d ln σ/d ln M term in Equation 1. This term is susceptible to numerical inaccuracies due to the combination of a sharp cut-off at k min and the oscillatory shape of the top-hat filter (Equation 5 ). Nevertheless, our finite-volume correction does not perfectly account for the real effects of a limited box size. For example, Power & Knebe (2006) show that the mass function can be enhanced at intermediate masses due to halos whose mass is reduced by the missing spectral power. We quantitatively analyze the effects of our finite-volume correction in Section 5.1. Note. -The N-body simulations used in this paper. L denotes the box size in comoving units, N 3 the number of particles, m p the particle mass, and the force softening length in physical units. The redshift range of each simulation is determined by the first and last redshifts z initial and z final , but snapshots were output only between z f−snap and z final . The structure in the earlier snapshots of some simulations is not developed enough to yield any halos yet, the first catalog with halos is output at z f−cat . The references correspond to Diemer et al. (2013, DKM13) , Diemer & Kravtsov (2014, DK14) , and Diemer & Kravtsov (2015, DK15) . Our system for choosing force resolutions is discussed in DK14.

The reader may wonder why we did not impose an upper limit on the σ integral due to the limited resolution of the simulations; the answer is not obvious. The highest wave mode represented is the Nyquist frequency, k max = πN/L. This frequency corresponds roughly to a smallest Lagrangian scale that we can resolve, R ny = 2L/N. The scales we wish to resolve will be limited by a minimum number of particles in a halo, N min (Sections 3.3 and 3.4), resulting in a minimum Lagrangian mass M min = N min ρ m,0 L 3 /N 3 and a minimum Lagrangian radius

We have inserted 500 particles as a typical minimum mass because that is the threshold chosen in Section 3.3. Regardless of cosmology, redshift, or box size, the unresolved scale is a factor of a few smaller than the smallest halos we wish to resolve. We have confirmed that the upper cutoff does not change σ appreciably (for N min = 500), even in cosmologies with shallow power spectra. However, we caution that, due to the large k-extent of the top-hat filter, small changes become noticeable for resolution limits of N min < ∼ 100.

The field of fitting mass functions began with the Press-Schechter prediction that f (ν) ∝ ν×exp(−ν 2 /2). This function provides an impressive fit for a theoretical Ansatz without free parameters but had to be amended to fit simulations in detail. One particularly useful four-parameter extension is

which was used by Warren et al. (2006) and Tinker et al. (2008) and is itself based on expressions in Sheth & Tormen (1999 , 2002 . As previously discussed, σ and ν are not related by a constant factor in this work, making it more natural to express the fitting function in terms of peak height,

with the simple parameter transformations b = δ c /b and c = c /δ 2 c .

3. SIMULATIONS AND ALGORITHMS In this section, we briefly review our N-body simulations. We describe our procedure for computing mass functions and show how to combine the results from multiple simulations. This work is based on the halo catalogs presented in Paper III, we refer the reader to that paper for details.

We use the Erebos suite of N-body simulations whose most important properties we review in Table 1 . This collection of simulations allows us to explore a WMAP7 ΛCDM cosmology based on that of the Bolshoi simulation (Komatsu et al. 2011; Klypin et al. 2011 , Ω m = 0.27, Ω b = 0.0469, h = 0.7, σ 8 = 0.82, and n s = 0.95) as well as a Planck-like cosmology (Planck Collaboration et al. 2014 , Ω m = 0.32, Ω b = 0.0491, h = 0.67, σ 8 = 0.834, and n s = 0.9624). These cosmologies span the likely ranges of those parameters that are most important for structure formation, such as Ω m . We also consider self-similar Einstein-de Sitter universes with power-law initial spectra of slopes −1, −1.5, −2, and −2.5. The power spectra for the ΛCDM simulations were computing using Camb (Lewis et al. 2000) , and the initial conditions for all simulations were generated by 2LPTic (Crocce et al. 2006 ). The simulations were run with Gadget2 (Springel 2005) .

Given that we are concerned with the question of universality, the self-similar simulations play an important role in our investigation because they test power spectra than differ drastically from ΛCDM. In these universes, both the power spectrum and time are independent of physical scales, meaning that properties such as the mass function must be universal with redshift except for resolution effects (Efstathiou et al. 1988; Lacey & Cole 1994; Lee & Shandarin 1999; Knollmann et al. 2008; Elahi et al. 2009; Leroy et al. 2020; Joyce et al. 2020) . As self-similar simulations are rarely visualized, Figure 1 shows density slices through our four simulations at their final redshifts (see Efstathiou et al. 1988 and Elahi 2009 for older examples). The scales at which matter clusters visi-4 Splashback Mass Function n = −1 n = −1.5 n = −2 n = −2.5 Figure 1 . Density slices through self-similar universes with power spectrum slopes of n = −1, −1.5, −2, and −2.5. The simulations are shown at their final redshifts of 2, 1, 0.5, and 0, respectively. The slices are 5 h −1 Mpc thick, the color scale corresponds to the projected density on a logarithmic scale between 3 × 10 −3 and 100 with respect to the cosmic mean. The simulations were normalized to the same σ 8 at z = 0 but exhibit very different distributions of structure: while shallow power spectra lead to lots of small-scale and little large-scale structure, the structure in the n = −2.5 simulation collapses at almost the same time for large and small scales (a slope of −3 would correspond to all scales collapsing at the same time and can thus not be simulated with a finite box). The largest structures look similar in shape because the same random phases were used in all simulations. The visualizations were created using the gotetra code by P. Mansfield (https://github.com/phil-mansfield/gotetra), which is based on a tetrahedron-based estimate of the density field Abel et al. 2012; Hahn et al. 2013) . Higher-resolution images of the self-similar and ΛCDM simulations are available at benediktdiemer.com/visualization. bly evolve from the shallowest to the steepest power spectrum slopes.

We use the halo catalogs presented in Paper III for this work, which we briefly review here. We start with catalogs created by the Rockstar and Consistent-Trees codes (Behroozi et al. 2013a,b) . Rockstar finds the particles in friends-of-friends groups in six-dimensional phase space but outputs SO masses and radii, as well as numerous other halo properties. We use Rockstar to compute bound-only masses, that is, SO masses where gravitationally unbound particles are excluded. We have run the Sparta code (Paper I) on all simu-lations to compute the splashback and all-particle SO masses that will be the focus of this paper. We then combined the Sparta data with the original catalogs, creating enhanced catalogs with a large number of mass definitions. When calculating mass functions, we wish to consider only isolated (host) halos. Thus, the radius (or mass) definition impacts the mass function in two ways: it shifts halos in mass and changes the number of host halos. We have computed separate hostsubhalo relations for each mass definition (Paper III). . If converged, the mass functions are expected to overlap, either between different box sizes in ΛCDM (top four panels for WMAP7, bottom-left panel for Planck) or between different redshifts in a single self-similar universe (three bottom right panels). The smaller bottom panels show the fractional difference to a fit (using Equation 10, gray dashed lines). We show only bins where the mass function is constrained by at least 100 halos to avoid crowding the plot with noisy measurements. We have picked representative examples of mass definitions and redshifts, the results are similar for other definitions. Generally, the mass functions are converged to within the statistically expected uncertainty, partly because we have imposed the appropriate mass cuts at the low-mass end. We do notice slight convergence issues at high redshift (top right panel) and for the n = −1 self-similar simulation (bottom center panels), although it is not clear that the deviations are statistically significant or systematic. In the following figures, we combine the data in each panel into a single sample per cosmology, mass definition, and redshift (the latter only for ΛCDM).

or 500 times the critical or mean matter density of the Universe. M vir corresponds to an evolving threshold as approximated by Bryan & Norman (1998) . Before computing their mass functions, we need to consider the resolution limits of our simulations, the completeness of our catalogs, and the accuracy to which masses are measured. The latter is particularly important when measuring steep mass functions, where scatter can cause Eddington bias (Eddington 1913 ). We impose a lower limit of 500 particles regardless of the overdensity threshold. While this choice leads to slightly different cutoffs compared to a fixed mass definition, it limits the statistical variance due to particle noise (Trenti et al. 2010; Benson 2017) . We have derived this limit by overlaying the mass functions from simulations with different resolutions and find that they are converged to within the statistical uncertainties for almost all SO definitions, redshifts, and cosmologies ( Figure 2) , with the possible exception of n = −1 and n = −2.5 which exhibit noticeable, although not highly significant, deviations. A number of analyses have found that mass functions can be converged for more ambitious limits (e.g., Klypin et al. 2015; Leroy et al. 2020 ), but we find significant non-convergence at the low-mass end when including halos with, say, 200 particles. Moreover, there is no need to push to the smallest halos in each simulation because our staggered box sizes ensure that the mass range is covered relatively evenly (Table 1) .

When counting halos to compute mass functions (as opposed to taking averages of halo properties), we need to worry about any incompleteness. In rare cases, the computation of SO masses fails because the density never reaches the threshold at the halo center (most commonly for M 500c ) or because the density never falls below the threshold at large radii (most common in low-density definitions such as M 200m ). The latter effect predominantly occurs for all-particle masses, the former also affects bound-only masses. However, both effects are much more important for subhalos and thus do not pose 6 Splashback Mass Function a significant problem for this investigation. We find that less than 10 −4 of our halos are typically affected, with the fraction rising to 10 −3 for R 500c and in the self-similar simulations. Given the much larger systematic uncertainties, we make no attempt to correct for this issue.

Not all SO measurements are physically meaningful. For instance, the all-particle masses of a small fraction of halos include a large amount of material from another, nearby halo. In Paper III, we found that this issue is most common in small boxes and at high redshift, which could introduce slight redshift-dependent shifts into our mass functions. We identify the affected halos by comparing the all-particle and boundonly masses of the same definition and use the bound-only mass if the ratio is greater than two. This limit is arbitrary, but there is no clear boundary between "normal" mass ratios and pathological cases. Our choice echoes Rockstar's algorithm, which discards halos as unphysical if more than half of their FOF group is unbound. Our procedure shifts the mass functions by up to a few percent, but these shifts are mostly uniform with redshift or improve the universality of SO definitions. Excessively large all-particle radii also include some spurious subhalos that would otherwise be hosts, a small effect that we cannot correct for here.

Despite these issues, all-particle SO masses are preferred over bound-only masses because they are uniquely defined. In contrast, bound-only masses depend intimately on the halo finder and unbinding algorithm (Paper III), making it harder to establish the connection to theoretical models and observations. We have checked that our conclusions hold for boundonly definitions as well. The figures in the paper show allparticle SO mass functions; labels such as M 200c refer to the all-particle version unless otherwise specified.

Sparta measures splashback masses by following the orbits of all particles in a halo and finding their first apocenter. From the collection of the apocenter times and radii, the splashback radius and mass of a halo are derived by slightly smoothing the distribution in time, taking into account both past and future splashbacks. The splashback radius can then be defined as the mean of the particle distribution, R sp,mn , the median, R sp,50% , or higher percentiles such as R sp,75% and R sp,90% . The latter is the largest percentile computed by Sparta and closest to the non-spherical splashback shells of Mansfield et al. (2017) . At the final snapshots of a simulation, the Sparta estimates become less accurate because there are no future particle orbits that can be taken into account. Sparta corrects for this bias, but the algorithm does incur increased scatter (Paper I). Thus, we restrict ourselves to z ≥ 0.13 for our fits and convergence tests.

We impose a lower limit of 1000 particles on our splashback masses, as suggested in Paper I. Again, we have tested this limit by comparing simulations of different resolutions and find that they are converged, with the previously mentioned caveat regarding high redshifts and the most extreme self-similar universes (Figure 2 ). The splashback calculation can fail for a number of reasons, most commonly insufficient particle splashback events. As described in Paper III, we increase the completeness by interpolating and extrapolating splashback masses across time if possible. If this interpolation fails, we guess a value based on the fitting function of Paper III. Out of the halos that pass the particle number cut, between 0.1% and 1% have masses not directly computed by Sparta, with the fraction increasing towards small box sizes (low halo masses). However, the vast majority of those masses were computed from extrapolation of nearby past or future values, which gives fairly accurate estimates (Paper III). Less than 0.1% of halos had their masses estimated by the fitting function, whose predictions suffer from large scatter; these numbers do not strongly depend on redshift. In principle, even a small number of poor estimates could up-scatter masses into higher bins with intrinsically low counts. We have confirmed that removing all halos without Sparta-computed splashback masses has no appreciable effect on our results, even for the n = −1 self-similar simulation that suffers from the worst incompleteness effects (Paper III). We conclude that the incompleteness of our catalogs does not constitute a serious source of error. Figure 2 demonstrates that our mass functions for individual simulations are converged at fixed peak height. For the remainder of the paper, we will combine the ΛCDM mass functions from multiple simulations into one sample per redshift and cosmology. These combined mass functions cannot be converted between mass and peak height any longer due to the simulation-specific corrections of Section 2.2. Thus, we combine different simulations in νf (ν) space. For each comparison, we define a logarithmically spaced set of peak height bins, ν i . For a given simulation and mass definition, we mask out bins that are not fully above the lower particle number limit. We convert all halo masses to peak height using Equation 2 and compute

where n i, j is the number of halos in simulation j and bin i, V j is the volume of simulation j, and ν i+1 and ν i are the upper and lower edges of the bin. The index j is understood to run only over simulations that contribute to the mass bin (or over only one simulation in the case of Figure 2 ). We estimate the fractional error on f (ν) from the uncertainty on the total number count in each bin,

This estimate of the uncertainty approaches 1/ √ n at large n but does not reach zero for n = 0 and reproduces the asymmetry of the Poisson distribution (Lukić et al. 2007 ). It does not, however, include other systematic errors or sample variance (Hu & Kravtsov 2003) .

For the self-similar simulations, we wish to combine different redshifts of the same simulation rather than different simulations at fixed redshift (bottom right panels of Figure 2 ; see, e.g., Lacey & Cole 1994; Lee & Shandarin 1999) . Thus, the volume weighting in Equation 11 will not work, but we need to weight the different redshifts because they contribute very differently at different peak heights (with larger peak heights collapsing earlier). Thus, we replace V j with the number of halos, n i, j (Lacey & Cole 1994) . This procedure is biased because bins that randomly fluctuate upwards in n i, j will receive systematically higher weight than bins that fluctuate downwards. In practice, this bias is negligible because most bins are dominated by one or two redshifts with large number counts and small statistical uncertainty. The uncertainty of the Splashback Mass Function They are compared to a fit to the z = 0 data, which has no particular significance but serves as a smooth reference for comparison. All SO definitions exhibit significant non-universalities with redshift. The second row shows four splashback definitions, which we compare to our universal fitting function (dashed gray lines, Section 4.3). The splashback mass functions exhibit non-universalities at high peak height and redshift but are universal to about 10% at z ≤ 2.

weighted average is then

where N z is the number of redshifts we have averaged over (see Leroy et al. 2020 , for a more sophisticated method to test the convergence of self-similar mass functions). In principle, we could consider all snapshots of a self-similar simulation, but that would lead to highly correlated measurements. Instead, we assume that halo masses change significantly over one dynamical time (a crossing time defined as in Paper I).

We avoid the final 0.3 dynamical times due to the late-time correction to the splashback masses (Section 3.4) and move backwards in steps of one dynamical time. Figure 2 ). In the following figures and fits, we omit bins where the mass function is constrained by fewer than 200 halos to avoid noisy measurements. We plot each ν-bin at its logarithmic center (as opposed to the mean peak height of halos in the bin).

We are now in a position to investigate the universality of mass functions by comparing them between different redshifts and cosmologies (Sections 4.1 and 4.2). We encapsulate the splashback mass functions in a simple, universal fitting function (Section 4.3).

We begin by comparing mass functions at different redshifts in the WMAP7 cosmology (Figure 3) . In each panel, we compare multiplicity functions for a given mass definition. The bottom panels show the fractional deviation from a reference function, for which we use a fit of Equation 10 (dashed lines) because the binned data do not span the entire ν range and because they suffer from shot noise. Any disagreements between the ratios shown in the bottom panels indicate nonuniversality. All cosmologies Figure 4 . Non-universality of the multiplicity function for different mass definitions, taking into account the WMAP7 and Planck cosmologies (left panel) and all cosmologies including self-similar universes (right panel). The relative differences appear slightly larger than in Figures 3 and 6 because they are defined with respect to the minimum f (ν) rather than to some fitting function and because we include the Planck cosmology. In ΛCDM, the non-universality grows with ν for all definitions, partly because high-redshift data are only available at higher peak height. M 200m is non-universal to 20% over a wide range of peak heights whereas M vir and the splashback definitions are more universal up to ν ≈ 2 and then rise to non-universalities between 30% and 60%. When we take into account self-similar cosmologies, the splashback definitions are clearly more universal than their SO counterparts.

The top row shows the all-particle versions of four SO definitions. The comparisons are very similar for the bound-only results or for the Planck cosmology (for which we have fewer simulations and thus less coverage of ν-z space). We compare the mass function at each redshift to a fit to the z = 0 mass function (blue dashed lines). This fit has no physical significance, it is merely a convenient point of comparison. Clearly, all SO mass functions suffer from significant non-universality with redshift. For all but M 200m , f (ν) increases with redshift, especially at the high-ν end. The trend is particularly strong between z = 0.5 and z = 0, indicating that it is caused by the accelerated decline of Ω m (z) due to dark energy. This trend holds, to a lesser extent, even for the evolving ∆ vir threshold. M 200m shows a distinctly different redshift evolution, where the normalization of the mass function slowly decreases with redshift, reaching an offset of 20% at z = 6 (in agreement with Tinker et al. 2008) . The M 200c and M 200m mass functions are identical at high redshift where ρ c ≈ ρ m , the difference is in their evolution towards low z.

However, it can be misleading to make quantitative inferences about universality based on Figure 3 because the eye is drawn to noisy outlier bins and deviations from the arbitrary fitting functions. To accurately quantify the differences, we define the non-universality η as the fractional difference between the highest and lowest f (ν),

where f min and f max are the lowest and highest multiplicity functions in a given ν-bin and σ min and σ max are their uncertainties. The error term ensures that only statistically significant differences are considered, i.e., that η = 0 if f min and f max are compatible within one standard deviation. Given that the uncertainties do not include systematic errors, η should be seen as an upper limit. The left panel of Figure 4 shows η as a function of peak height for our ΛCDM cosmologies. We have included redshifts 0-6 in the WMAP7 cosmology and 0-4 in Planck, the range where we have reliable data. We use the same peak height bins as in Figure 3 . In a given bin, we remove mass functions measured with fewer than 100 halos as in Figure 2 ; compared to the stricter 200 halo limit, this choice extends the curves slightly at high ν. Moreover, we have removed a spurious drop of η at the highest peak heights that can occur because the mass functions of certain redshifts are not constrained any longer. To avoid noisy curves, we smooth the results with a 4th-order Savitzky & Golay (1964) filter with a window size of 11 bins. We have confirmed that this smoothing does not change any of the noteworthy features in Figure 4 . Most importantly, it does not systematically change the end points of the curves where they reach the largest non-universalities. Nevertheless, our definition of η is far from unique. With higher-resolution simulations, for instance, we could push to higher redshift and might infer larger non-universalities. Figure 4 is thus intended as a relative comparison of the different mass definitions. The η curves in the left panel of Figure 4 reflect our previous conclusions regarding SO mass functions but further illuminate their behavior at high ν: while M vir is universal to about 10% up to ν ≈ 2.5, its non-universality quickly increases with ν and reaches about 50% (in agreement with Figure 2 of Despali et al. 2016) . Figure 4 also allows us to quantify the non-universality of M 200c and M 500c , which is apparent even at low ν and reaches factors of more than unity.

We now turn to the splashback definitions (second row of Figure 3 ). Here, we compare to the universal fitting function introduced in Section 4.3 (gray dashed lines). In contrast to the SO definitions, this fit is performed at z ≤ 2 rather than z = 0, but we are interested in the relative differences between mass functions, not in the fit quality (which we consider in Section 4.3). For a fair comparison with the SO definitions, we use the z = 0 splashback masses rather than z = 0.13 despite the caveat mentioned in Section 3.4; we have checked Splashback Mass Function that the z = 0.13 results are very similar. All splashback mass functions exhibit excellent universality up to z ≈ 2, with the redshifts agreeing to about 10%. At z ≥ 3 and ν > 2, we notice a relative decline in the mass function. The nature of this high-z decrease is fundamentally different from the redshift evolution of the SO mass functions in that it clearly does not depend on dark energy. Instead, it could indicate an issue in Sparta's algorithm at high z or a genuine non-universality in the splashback mass functions, perhaps caused by the rapidly changing mass accretion rates at high z. Our conclusions are confirmed in Figure 4 , where the splashback definitions top out between 30% and 45% (though their differences at high ν are not particularly significant). Overall, M sp,90% exhibits the highest level of universality over most of the peak height range. The values of η are a little larger than one might infer from Figure 3 because we are dividing the difference by the lowest f (ν) now instead of the the universal fitting function. In summary, the mass functions of no definition are truly universal across all redshifts in ΛCDM. While high-threshold SO definitions such as M 200c and M 500c exhibit a large redshift evolution, M 200m , M vir , and M sp lead to comparable levels of non-universality with somewhat different dependencies on peak height. However, the M sp functions are noticeably more universal at z ≤ 2.

While universality with redshift is an important consideration, it is by no means sufficient: for f (ν) to be truly universal, it must hold across cosmologies with very different power spectra. ΛCDM cosmologies such as WMAP7 and Planck do not differ nearly enough in their spectra to claim such universality (e.g., Brown et al. 2020 ). Thus, we now consider the mass functions in our self-similar simulations ( Figure 5 ). For M 200m (left panel), we arbitrarily use a fit to the n = −1 case as a reference (dashed yellow line), for the splashback definitions we use the universal fitting function (dashed gray lines). All SO definitions exhibit significant non-universality with n; we show M 200m because we found it to be most universal in the ΛCDM case. Steeper power spectra result in shallower mass functions and thus a relative increase of the mass function at high masses. The splashback definitions fare much better, especially M sp,90% where the different cosmologies agree to about 15% or better at fixed ν. This agreement is remarkable given how different the formed structures are in our self-similar universes (Figure 1) .

Having established the basic trends, we now attempt a more ambitious test of universality: directly comparing the different cosmologies and redshifts to each other. To avoid crowding the plots with too many lines, Figure 6 shows a selection of the most extreme cases, namely low and high redshifts from WMAP7 and Planck and the self-similar simulations with the shallowest and steepest slopes. As before, we compare the SO mass functions to a fit to the z ≈ 0 WMAP7 data and the splashback mass functions to the universal fitting function. This comparison highlights that the non-universality of the SO mass functions is more significant than suggested by Figures 3 and 5 individually: in all SO definitions, the mass functions in n = −1 and n = −2.5 behave qualitatively and quantitatively differently from the low-z ΛCDM case. Even for M 200m , they strongly disagree at the high-ν end, and yet both disagree with ΛCDM. Again, we note that the M 200c and M 200m curves are the same for the self-similar simulations, they only appear different in relation to ΛCDM. By comparison, the splashback definitions are much more universal with cosmology: the differences between ΛCDM and self-similar cosmologies are no larger than the differences between redshifts in ΛCDM. The n = −1 mass functions tend to be slightly higher than ΛCDM, especially at the low-ν end. This effect could be caused by their sparse snapshot sampling of only about 3.5 snapshots per dynamical time, which was shown to cause a small bias (Paper I), or by residual convergence issues in this particular simulation (Figure 2) .

The conclusions based on Figure 6 are quantified by the η curves shown in the right panel of Figure 3 but combining different cosmologies and redshifts, namely the WMAP7 and Planck cosmologies at z ≈ 0 and z = 4 as well as the self-similar cosmologies with the most extreme power spectra. This selection essentially brackets the parameter space we have investigated. The comparison conclusively demonstrates that none of the SO mass functions are universal; they either exhibit strong redshift evolution, are incompatible with the self-similar simulations, or both. The splashback mass functions are much closer to universality, especially for the highest percentiles.

behaving qualitatively similarly to M 500c and M 200c . We conclude that the mass functions of any SO definition are not universal with cosmology. The splashback mass functions remain within η < ∼ 15% below ν = 2, while the non-universality increases towards high ν and reaches 40% for M sp,90% and up to 60% for the other percentiles. This peak non-universality is driven by the redshift differences in ΛCDM, not by differences to the self-similar cosmologies.

The universality of the splashback mass functions is remarkable given that they are computed with a fairly complicated algorithm. We emphasize that Sparta was not tuned to result in a universal mass function; the halo catalogs of Paper III were created before measuring the mass functions for this paper. The slight improvement in universality towards the largest splashback radii suggests that the physically most meaningful halo boundary resides near the sharp drop in density in the outskirts Xhakaj et al. 2019 ).

We now construct a universal fitting function for the splashback mass function that is valid regardless of cosmology and redshift (to the accuracy implied by Figures 3, 5, and 6) . We use Equation 10 with parameters that vary depending on the definition of the splashback radius. Following Paper II and Paper III, we provide one set of parameters for M sp,mn and one for all higher percentiles. We parameterize the latter dependence as p, the percentile divided by 100,

The best-fit values of the seven free parameters are given in Table 2 . They are derived from a fit to the binned mass functions at redshifts 0.13, 0.5, 1, and 2. We exclude z = 0 due to the late-time corrections (Section 3.4) and higher redshifts due to their slight non-universality. We include both the WMAP7 and Planck samples with equal weight per bin, but we omit the self-similar simulations since we are most interested in a good fit to ΛCDM. For the percentile-based definitions, we fit to all percentiles included in our halo catalogs (50, 70, 75, 80, 85, and 90) . We have experimented with up-weighting the 50th percentile to account for the gap at low percentiles, but we find that it has little impact. With our dataset thus defined, we constrain all parameters simultaneously using a Levenberg-Marquart least-squares minimiza- tion. We add a systematic uncertainty of 3% in quadrature to the statistical uncertainties shown in the figures. This extra error prevents the most statistically significant bins from dominating the fit and (artificially) results in a roughly unity χ 2 per degree of freedom. The accuracy of our fitting function can be read off the bottom panels in Figures 3, 5, and 6 . In ΛCDM and below ν ≈ 2, the formula is accurate to about 5% or better for all splashback definitions. At the highest peak heights, the nonuniversality with redshift means that the accuracy is degraded to about 10% for z < ∼ 2 and about 20% for z > 2. For the self-similar simulations, the accuracy depends more strongly on the definition, with M sp,mn fit to about 20% accuracy and M sp,90% to 15% or better regardless of the power spectrum slope. We have implemented our function in the publicly available python code Colossus (Diemer 2018) .

We note that the integral over our fit diverges because Equation 10 approaches A as ν → 0. It is possible to reparameterize the function to avoid this divergence (Tinker et al. 2008 ), but we do not attempt this conversion. Integrating the mass function from ν = 0.5 to infinity, we find that between 50% and 60% of dark matter resides within the splashback radius, depending on the percentile.

We have systematically investigated the universality of the multiplicity function and found it to strongly depend on the mass definition. As our results rely on an accurate conversion from mass to peak height, we investigate the subtleties of this calculation in Section 5.1. We discuss the theoretical and practical implications of our findings in Section 5.2 and consider possible directions for future research in Section 5.3.

Throughout the paper, we have defined non-universality as differences in the multiplicity function at fixed peak height. Given the steep fall-off at high ν, however, relative changes can equally be thought of as differences in ν at fixed number density. Thus, it is of paramount importance to calculate ν(M) in a physically meaningful and numerically accurate manner. Moreover, the details of the calculation vary across the literature, making it difficult to cross-compare results. In this section, we thus discuss the corrections we applied when computing σ and δ c , the filter function, and different strategies for computing the underlying power spectrum. Figure 7 shows a graphical impression of their impact on the universality of M sp,90% ; the effect on other definitions is similar.

The first panel shows the multiplicity functions at different redshifts in the WMAP7 cosmology, essentially the same as in the bottom-right panel of Figure 3 . Here, however, we have not corrected for finite box size and the collapse threshold is constant, δ c = δ c,EdS . In the second panel, we have applied the lower limit to the σ integral (Section 2.2), which accounts for the decreased variance on large scales due to missing low-k modes. This correction increases ν for the most massive halos in a given box, and brings the redshifts into better agreement.

In the third panel, we additionally evolve δ c with redshift to account for dark energy (Equation 7). This correction makes the mass functions in all definitions moderately more universal at the high-ν end. The correction operates only at z < ∼ 2 in ΛCDM because Ω m (z) approaches unity at higher redshifts. Our prescription for δ c (z) is not unique; numerous works have experimented with different values of δ c , its evolution, and dependencies on factors such as the density environment (e.g., Monaco 1995; Gross et al. 1998; Tormen 1998; Lacey & Cole 1994; Desjacques 2008; Courtin et al. 2011; Paranjape et al. 2013) . In this work, however, we do not treat δ c (z) as a free parameter. First, changing its normalization has no impact on the universality. Second, we are interested in how well different mass definitions capture the total mass that has collapsed into halos, not in modeling the physics that lead to particular functional forms of δ c (z). We thus impose the evolution predicted by the standard top-hat collapse collapse model in the presence of dark energy (Gunn & Gott 1972; Eke et al. 1996) . While this calculation can be complicated (e.g., Courtin et al. 2011; Corasaniti & Achitouv 2011; Benson et al. 2013) , the basic evolution of δ c (z) with Ω m (z) does not change even when taking into account minute algorithmic details (Pace et al. 2017) . Thus, alternative δ c (z) calculations could change our conclusions only if they include additional physics beyond the top-hat collapse model.

In the fourth and fifth panels of Figure 7 , we explore the effect of different power spectrum measurements. Reed et al. (2007) advocated summing over the measured power spectrum of the initial conditions of each simulation to capture the random deviations of modes from the input power spectrum. Comparing the third and fourth panels, we find that this procedure makes virtually no difference; we thus use the Camb spectra throughout the paper (see Despali et al. 2016, for similar findings) . This conclusion, however, does not imply that the exact shape of the power spectrum is not important: in the fifth panel, we have approximated P(k) using the Eisenstein & Hu (1998) transfer function, which is sufficiently accurate for most other applications but causes visible differences in the universality.

One final aspect of the calculation that has received renewed attention is the filter function in Equation 4. This filter ensures that σ(M) refers to those spatial scales that actually contribute to the progenitor density fields of halos. By using the top-hat filter (Equation 5), we implicitly assume that the material in halos originated from spherical regions with sharp boundaries in real space. As a first check, we have computed f (ν) with a Gaussian filter, W ∝ exp[−(kR) 2 /2], and a sharp-k filter (a top-hat in frequency space). Both lead to drastically non-universal mass functions (bottom panel of Figure 7 ). This result is somewhat expected because the initial peaks of simulated halos were found to have an intermediate shape somewhat closer to a top-hat than a Gaussian (Dalal et al. 2010; Chan et al. 2017) . Ideally, the filter should match this shape, leading Chan et al. (2017) to propose a filter that transitions between a top-hat and Gaussian. Their function has free parameters that depend on ν, making the predictions somewhat difficult to compute. Nevertheless, we have experimented with the Chan et al. (2017) filter using the range of parameter values they suggest. We find that the mass functions move away from universality as the filter moves away Gaussian filter Figure 7 . Impact of details of the peak height calculation on the universality of M sp,90% ; the conclusions are similar for other mass definitions. Each panel shows the ratio of mass functions in the WMAP7 cosmology compared to a fit at z = 0.13, the colors have the same meaning as in Figure 3 (redshifts 0.13, 0.5, 1, 2, 4, and 6 from blue to yellow). In the top panel, we use a fixed δ c and no finite volume correction. In the second panel, we impose a lower limit on the σ integral (Section 2.2), which has a visible impact. In the third panel, we also take into account the evolution of δ c in ΛCDM cosmologies (Equation 7), which causes a slight relative shift between redshifts. In the fourth and fifth panels we explore different ways to calculate the power spectrum. First, we directly sum over the modes of the measured power spectrum of the initial conditions rather than the Camb input spectrum; the differences are negligible. In the fifth panel, we approximate the power spectrum with the transfer function of Eisenstein & Hu (1998) . Although this fit is accurate to about 5%, it leads to significant relative changes in the mass functions. In the bottom panel, we replace the top-hat filter with a Gaussian, leading to strong non-universality (note the different scales on both axes). See Section 5.1 for a detailed discussion of these effects.

from a top-hat. Similarly, Leo et al. (2018) propose a smooth analogue to the sharp-k filter. We have tested their filter with a range of parameters recommended in their paper but find it to degrade the universality. In summary, the top-hat filter does not exactly capture the shape of peaks in the early Universe, but it is difficult to improve it, not least because the shape of peaks depends on their peak height (Bardeen et al. 1986 ). The preceding discussion highlights why self-similar universes are so useful in testing universality: many of the effects we have discussed do not apply to them. In particular, the input power spectrum is unambiguous, δ c is a constant, and finite volume corrections barely affect universes with n > ∼ − 2 because the overall power is dominated by small scales. As a result, the comparisons of self-similar simulations in Figures 5 and 6 are highly robust to the details of the ν calculation. The filter function, however, matters even in self-similar universes. While their variance is always σ 2 ∝ R −n−3 , its normalization depends on n. As in ΛCDM, using a filter other than top-hat leads to much larger differences between different power spectrum slopes, confirming that the top-hat filter is a reasonable choice (see Lacey & Cole 1994, for similar experiments) . Finally, the relative universality of the splashback mass functions in self-similar universes demonstrates that there is no need to include n as a parameter in our fitting function (cf. Reed et al. 2007 ).

Our key result is that the splashback mass function is surprisingly universal, especially for definitions that enclose a large radius such as M sp,90% . In this section, we consider the implications of this finding. In particular, we ask whether we have any reason to expect such universality, whether it could be replicated with other mass definitions, and what the practical uses of a universal mass function might be.

Throughout the paper, we have highlighted numerous reasons why we might not expect the mass function to be universal: the shape of the initial peaks influences their collapse and depends on peak properties beyond ν ( Bardeen et al. 1986; Bond & Myers 1996) ; the large extent of the filter function in k-space introduces a dependence on the power spectrum; the collapse of halos is more complicated than suggested by the spherical top-hat collapse model; and, not least, mass definitions may not capture the total mass collapsed into halos. Thus, even with a "perfect" definition of the halo boundary (in the sense of the spherical collapse model), we have no reason to expect a perfectly universal mass function -which makes the agreement of f (ν) between different redshifts and cosmologies even more remarkable.

Could this success be replicated with other definitions of the halo boundary, such as a low-threshold SO mass? Such an equivalency seems unlikely because the overdensity enclosed within R sp , ∆ sp , varies between about 50 and 400ρ m and depends strongly on mass accretion rate, redshift, and halo mass (Paper II). Low SO thresholds are also problematic in practice because they can merge neighboring halos, whereas overlapping splashback radii pose no issue. Another possibility is a threshold that varies with redshift like ∆ vir , but its evolution would not affect the universality in the self-similar cosmologies. Similarly, the agreement between FOF mass functions can be improved by changing the linking length, but the chosen value would need to depend on cosmology, mass, and redshift (Courtin et al. 2011; More et al. 2011) .

What are the practical implications of a universal mass function? We might hope to compare the observed mass function of galaxy clusters to a universal form to avoid running simulations for numerous different cosmologies. In practice, however, this application seems unlikely. First, the level of non-universality in the splashback mass functions is still too large for precision cosmology (e.g., McClintock et al. 2019 ), although the M sp functions, unlike those of SO masses, are universal to better than 10% in the regime relevant for cluster observations (z < ∼ 2 and ν > ∼ 3 in ΛCDM). The accuracy of observational comparisons, however, is limited by systematic uncertainties: while we can quantify simulated mass functions to high accuracy (Murray et al. 2013) , we cannot actually measure the masses (SO or otherwise) of real clusters. Instead, they are inferred from proxies such as their X-ray signal or richness, with uncertainties that will likely dominate the comparison for the foreseeable future (e.g., Kravtsov & Borgani 2012; DES Collaboration et al. 2020) .

Leaving aside observational inferences, a universal mass function would have significant theoretical ramifications. The key assumption underlying most models of the mass function is that the physics of halo collapse is insensitive to the power spectrum and that any further dependencies (e.g., on the shape of peaks) can be expressed as a function of only peak height.

Our results indicate that real halo collapse is, to a surprising degree, captured by this simplistic picture. However, they do not necessarily support the specific Ansatz of Press & Schechter (1974) , where the mass in halos above a certain peak height is estimated as the fraction of the Lagrangian volume where the smoothed density field exceeds δ c . While the Press & Schechter (1974) prediction differs from our M sp,90% fit by only 20% at low peak heights (as opposed to the infamous 50% for FOF masses), it still drastically underpredicts the abundance of halos at high ν. In related work, Garcia et al. (2020) recently investigated a generic halo boundary optimized to fit the halo-matter correlation function, which appears to roughly match the splashback radius. They found that the corresponding mass function is closely approximated by Press & Schechter (1974) , although for a non-standard value of δ c . While this match is intriguing, our results demonstrate that the mass included within R sp of large halos is much larger than suggested by the Press & Schechter (1974) model.

The universality of the splashback mass function at low redshift motivates us to consider possible causes of the remaining non-universality, besides those discussed in Section 5.2. One such effect is non-sphericity: while real halos have long been known to be triaxial (Jing & Suto 2002) , their masses are still routinely measured within spherical apertures. Numerous works have attempted to include non-sphericity via the collapse model calculation (e.g., Eisenstein & Loeb 1995; Monaco 1995; Bond & Myers 1996; Audit et al. 1997; Lee & Shandarin 1998; Sheth et al. 2001; Sheth & Tormen 2002) , typically expressed as a dependence of δ c on ν. Such corrections are justified by simulation results (Robertson et al. 2009 ), but they affect the predicted shape of f (ν) rather than its universality. If, in reality, the average non-sphericity depended on redshift or cosmology as well, it could introduce a systematic shift in the mass measurements. Such a bias is likely to affect M sp in particular because the splashback shells around halos can be strongly non-spherical and are poorly fit by ellipsoids . While the corresponding change in M sp may be modest due to the low densities in the outskirts, non-spherical boundaries can drastically change the number of subhalos, which also influences the mass function (Mansfield & Kravtsov 2020) . We plan to update our catalogs with non-spherical halo boundaries in the future.

Similarly, it remains to be seen whether the universality of the splashback mass function is affected by baryonic effects (Rudd et al. 2008; Cui et al. 2012; Velliscig et al. 2014) , although the large extent of the splashback radius should make it less sensitive to baryons than SO definitions. Baryonic simulations could also bridge the gap to observations of the galaxy cluster mass function (Section 5.2). In particular, it would be instructive to compare the scatter in scaling relations between observables and M sp to those for SO masses. Since the scatter between splashback and SO masses is much larger than that in the scaling relations (Paper II), it is not obvious that tight relations will emerge. For example, one might expect that small radii such as R 500c will correlate more tightly with X-ray properties measured near the cluster center. Richness, on the other hand, might correlate better with the splashback radius, given that it is designed to include all subhalos (Diemer 2020b) .

Finally, we have not included simulations of more exotic universes with non-Gaussian random fields, modified gravity, or alternative dark matter models (e.g., Łokas et al. 2004; Dalal et al. 2008; Pillepich et al. 2010; von Braun-Bates & Devriendt 2018; Lovell 2020) . We leave investigations of the splashback mass function in such cosmologies for future work.

6. CONCLUSIONS We have systematically investigated the universality of halo mass functions in N-body simulations over a wide range of halo mass, redshift, and cosmology. For the first time, we have considered splashback masses as well as conventional SO definitions and find that their mass functions are surprisingly universal. Our main conclusions are as follows.

1. We confirm that SO mass functions in ΛCDM are generally not universal across redshifts. Out of the commonly used definitions, M 200m is most universal but exhibits deviations of about 20% across a large range of masses. M vir is most universal at low peak heights but reaches more than 50% non-universality at the highest masses. In contrast, the M 200c and M 500c mass functions are non-universal at all masses. The splashback mass functions are universal to better than 10% at z < ∼ 2 but exhibit non-universalities around 40% at the highest masses and redshifts. It is unclear whether this remaining non-universality is fundamental or due to measurement errors.

I am deeply grateful to Susmita Adhikari, Shaun Brown, Neal Dalal, Giulia Despali, Michael Joyce, Andrey Kravtsov, Philip Mansfield, Joseph Mohr, and Darren Reed for illuminating discussions and detailed feedback on a draft. This work was partially completed during the Coronavirus lockdown and would not have been possible without the essential workers who did not enjoy the privilege of working from the safety of their homes. All computations were run on the Midway computing cluster provided by the University of Chicago Research Computing Center. This research made extensive use of the python packages NumPy (Oliphant 2006) , SciPy (Virtanen et al. 2019) , Matplotlib (Hunter 2007) , and Colossus (Diemer 2018) . Support for Program number HST-HF2-51406.001-A was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

The splashback radius of halos from particle dynamics: III. Halo Catalogs, Merger Trees, and host-subhalo relations

Galaxy Formation and Evolution

Guide to NumPy