The splashback radius of halos from particle dynamics. II. Dependence on mass, accretion rate, redshift, and cosmology

TL;DR

This study investigates the splashback radius of dark matter halos using cosmological simulations, revealing complex dependencies on accretion rate, mass, redshift, and cosmology, and providing new fitting functions for these relationships.

Contribution

It introduces detailed analysis of the splashback radius dependence on multiple parameters and offers new fitting functions, highlighting its independence from traditional spherical overdensity definitions.

Findings

R_sp/R_200m decreases with accretion rate Gamma

Fitting functions achieve 5% accuracy across parameter space

R_sp as the steepest density slope encloses ~75% of apocenters

Abstract

The splashback radius $R_{\rm sp}$, the apocentric radius of particles on their first orbit after falling into a dark matter halo, has recently been suggested as a physically motivated halo boundary that separates accreting from orbiting material. Using the SPARTA code presented in Paper I, we analyze the orbits of billions of particles in cosmological simulations of structure formation and measure $R_{\rm sp}$ for a large sample of halos that span a mass range from dwarf galaxy to massive cluster halos, reach redshift 8, and include WMAP, Planck, and self-similar cosmologies. We analyze the dependence of $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ on the mass accretion rate $\Gamma$, halo mass, redshift, and cosmology. The scatter in these relations varies between 0.02 and 0.1 dex. While we confirm the known trend that $R_{\rm sp}/R_{\rm 200m}$ decreases with $\Gamma$, the relationships turn out to be more complex than previously thought, demonstrating that $R_{\rm sp}$ is an independent definition of the halo boundary that cannot trivially be reconstructed from spherical overdensity definitions. We present fitting functions for $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ as a function of accretion rate, peak height, and redshift, achieving an accuracy of 5% or better everywhere in the parameter space explored. We discuss the physical meaning of the distribution of particle apocenters and show that the previously proposed definition of $R_{\rm sp}$ as the radius of the steepest logarithmic density slope encloses roughly three-quarters of the apocenters. Finally, we conclude that no analytical model presented thus far can fully explain our results.

Figures

Tables

Table 1. Table 1 Definitions of the Symbols Used in This Paper

Symbol	Meaning
${\rho }_m$	Mean matter density of the universe
${\rho }_c$	Critical density of the universe
${\Omega }_m$	Fractional matter density, ${\Omega }_m\equiv {\rho }_m/{\rho }_c$
$r$	Some radius in physical units, measured from the halo center
$R$	A particular definition of the halo boundary
$\Delta$	An overdensity with respect to either ${\rho }_m$ or ${\rho }_c$
$R_{200m}$	Radius enclosing an overdensity of $200\times {\rho }_m$
$R_{200c}$	Radius enclosing an overdensity of $200\times {\rho }_c$
$R_{\mathrm{vir}}$	$R_{\Delta }$ with varying overdensity (Bryan & Norman, 1998)
$M_{200m}$	Mass inside $R_{200m}$
$N_{\Delta }$	Number of particles inside $R_{\Delta }$, e.g. $N_{200m}$
$\nu$	Peak height, $\nu \equiv {\nu }_{200m}={\delta }_c/\sigma (M_{200m},z)$
$r_s$	Scale radius of an NFW profile
$c_{\Delta }$	Concentration, $c_{\Delta }\equiv R_{\Delta }/r_s$
$r_{\mathrm{sp}}$	Splashback radius of a particle
$m_{\mathrm{sp}}$	Splashback mass of a particle, i.e.
$R_{\mathrm{sp}}$	Splashback radius of a halo
$M_{\mathrm{sp}}$	Splashback mass of a halo
${\Delta }_{\mathrm{sp}}$	Splashback overdensity wrt. ${\rho }_m$, ${\Delta }_{\mathrm{sp}}\equiv 3M_{\mathrm{sp}}/\left(4\pi R_{\mathrm{sp}}^3\right)/{\rho }_m$
$R_{\mathrm{sp}}^{\mathrm{mn}}$	$R_{\mathrm{sp}}$ defined as the mean of the particle $r_{\mathrm{sp}}$
$R_{\mathrm{sp}}^{50\%}$	$R_{\mathrm{sp}}$ defined as the median of the particle $r_{\mathrm{sp}}$
$R_{\mathrm{sp}}^{75\%}$	$R_{\mathrm{sp}}$ defined as the 75th percentile of the particle $r_{\mathrm{sp}}$
$R_{\mathrm{sp}}^{\ast }$	Summary symbol for multiple definitions of $R_{\mathrm{sp}}$
$v_{\Delta }$	Circular velocity, $v_{\Delta }\equiv \sqrt{G{M}_{\Delta }/R_{\Delta }}$
$t_{\mathrm{dyn}}$	Dynamical time or crossing time, $t_{\mathrm{dyn}}\equiv 2R_{200m}/v_{200m}$
$s$	Instantaneous mass accretion rate, $d\log \left(M\right)/d\log \left(a\right)$
${\Gamma }_{\mathrm{dyn}}$	Mass accretion rate over one $t_{\mathrm{dyn}}$, $\Delta \log \left(M\right)/\Delta \log \left(a\right)$
${\Gamma }_{\mathrm{DK14}}$	Mass accretion rate as defined in Diemer & Kravtsov (2014)

Table 2. Table 2 N 𝑁 N -Body Simulations

Name	$L$	$N^3$	$m_p$	$ϵ$	$ϵ/\left(L/N\right)$	$z_{\mathrm{initial}}$	$z_{\mathrm{final}}$	$N_{\mathrm{snaps}}$	$z_{f-\mathrm{snap}}$	Cosmology	Reference
L2000	$2000$	${1024}^3$	$5.6\times {10}^{11}$	$65$	$1/30$	$49$	$0$	$100$	$20$	$WMAP$ (Bolshoi)	DK15
L1000	$1000$	${1024}^3$	$7.0\times {10}^{10}$	$33$	$1/30$	$49$	$0$	$100$	$20$	$WMAP$ (Bolshoi)	DKM13
L0500	$500$	${1024}^3$	$8.7\times {10}^9$	$14$	$1/35$	$49$	$0$	$100$	$20$	$WMAP$ (Bolshoi)	DK14
L0250	$250$	${1024}^3$	$1.1\times {10}^9$	$5.8$	$1/42$	$49$	$0$	$100$	$20$	$WMAP$ (Bolshoi)	DK14
L0125	$125$	${1024}^3$	$1.4\times {10}^8$	$2.4$	$1/51$	$49$	$0$	$100$	$20$	$WMAP$ (Bolshoi)	DK14
L0063	$62.5$	${1024}^3$	$1.7\times {10}^7$	$1.0$	$1/60$	$49$	$0$	$100$	$20$	$WMAP$ (Bolshoi)	DK14
L0031	$31.25$	${1024}^3$	$2.1\times {10}^6$	$0.25$	$1/122$	$49$	$2$	$64$	$20$	$WMAP$ (Bolshoi)	DK15
L0500-Planck	$500$	${1024}^3$	$1.0\times {10}^{10}$	$14$	$1/35$	$49$	$0$	$100$	$20$	$Planck$	DK15
L0250-Planck	$250$	${1024}^3$	$1.3\times {10}^9$	$5.8$	$1/42$	$49$	$0$	$100$	$20$	$Planck$	DK15
L0125-Planck	$125$	${1024}^3$	$1.6\times {10}^8$	$2.4$	$1/51$	$49$	$0$	$100$	$20$	$Planck$	DK15
L0100-PL-1.0	$100$	${1024}^3$	$2.6\times {10}^8$	$0.5$	$1/195$	$119$	$2$	$64$	$20$	Self-similar, $n=-1.0$	DK15
L0100-PL-2.5	$100$	${1024}^3$	$2.6\times {10}^8$	$1.0$	$1/98$	$49$	$0$	$100$	$20$	Self-similar, $n=-2.5$	DK15

Table 3. Table 3 Cosmological Parameters

Cosmology	$H_0$	${\Omega }_m$	${\Omega }_{\Lambda }$	${\Omega }_b$	${\Omega }_k$	${\Omega }_{\nu }$	${\sigma }_8$	$n_s$	$P\left(k\right)$	Reference
$WMAP$ (Bolshoi)	$70$	$0.27$	$0.73$	$0.0469$	$0$	$0$	$0.82$	$0.95$	CAMB	Klypin et al. (2011), Komatsu et al. (2011)
$Planck$	$67$	$0.32$	$0.68$	$0.0491$	$0$	$0$	$0.834$	$0.9624$	CAMB	Planck Collaboration et al. (2014)
Self-similar	$70$	$1$	$0$	$0$	$0$	$0$	$0.82$	…	$P\left(k\right)\propto k^n$	…

Table 4. Table 4 Best-Fit parameters

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}}$
$a_0$	$0.6498$	$0.3203$	$0.6792$	$0.2648$
$b_0$	$0.6004$	$0.2674$	$0.4051$	$0.6660$
$b_{\Omega }$	$0.0920$	$0.1134$	$0.2919$	$0.1688$
$b_{\nu }$	$0.0616$	$0.2080$	$0$	$0$
$c_0$	$-0.8063$	$-0.9596$	$3.3659$	$4.7287$
$c_{\Omega }$	$17.5205$	$16.2459$	$1.4698$	$2.3889$
$c_{\nu }$	$-0.2935$	$0$	$-0.0756$	$-0.0841$
$c_{\Omega 2}$	$-9.6243$	$-9.4979$	$0$	$0$
$c_{\nu 2}$	$0.0392$	$-0.0185$	$0$	$0$
Meta-Parameters for Dependence on Percentile
$a_p$	$0$	$0.6148$	$0$	$0.8435$
$b_p$	$0$	$0.5452$	$0$	$-0.6392$
$b_{\Omega p}$	$0$	$0$	$0$	$0.0032$
$b_{\Omega \mathrm{p2}}$	$0$	$0$	$0$	$4.9393$
$b_{\nu p}$	$0$	$-0.2233$	$0$	$0.2254$
$c_{\Omega p}$	$0$	$0.0039$	$0$	$-0.7057$
$c_{\Omega \mathrm{p2}}$	$0$	$8.9691$	$0$	$-1.2419$
$c_{\Omega 2p}$	$0$	$-0.0005$	$0$	$0$
$c_{\Omega 2p2}$	$0$	$10.6132$	$0$	$0$
$c_{\nu p}$	$0$	$-0.4511$	$0$	$-0.3911$
$c_{\nu 2p}$	$0$	$0.0880$	$0$	$0.0742$
Parameters for 68% Scatter (in dex)
${\sigma }_0$	$0.0526$	$0.0445$	$0.0528$	$0.0276$
${\sigma }_{\Gamma }$	$0.0038$	$0.0044$	$0.0025$	$0.0023$
${\sigma }_{\nu }$	$-0.0121$	$-0.0146$	$-0.0112$	$-0.0125$
${\sigma }_p$	$0$	$0.0226$	$0$	$0.0473$

Equations

$$M_{\Delta }=\frac{4\pi }{3}\Delta {\rho }_{\mathrm{ref}}{R}_{\Delta }^3$$

$$\Gamma (a_1)\equiv \frac{\Delta \log (M)}{\Delta \log (a)}=\frac{\log (M_1)-\log (M_0)}{\log (a_1)-\log (a_0)},$$

$${\sigma }_{\mathrm{sp}}={\sigma }_0+{\sigma }_{\Gamma }\Gamma +{\sigma }_{\nu }\nu +{\sigma }_{\mathrm{p}}p$$

$${\sigma }_{{\Delta }_{\mathrm{sp}}}=\sqrt{{\sigma }_{M_{\mathrm{sp}}}^2+3{\sigma }_{R_{\mathrm{sp}}}^2}.$$

$$X_{\mathrm{sp}}=A+B{e}^{-\Gamma /C}$$

$$A$$

$$B$$

$$C$$

$$A_0$$

$$B_0$$

$$B_{\Omega }$$

$$B_{\nu }$$

$$C_0$$

$$C_{\Omega }$$

$$C_{\Omega 2}$$

$$C_{\nu }$$

$$C_{\nu 2}$$

$${\sigma }_{\Gamma }\approx 0.41-0.07\nu$$

$${\Gamma }_{\mathrm{dyn}}=A\nu +B{\nu }^{3/2}$$

$$A$$

$$B$$

Full text

††journal: The Astrophysical Journal, 843:140 (16pp), 2014 July 10

\SIthousandsep

,

1. Introduction

According to our current understanding of structure formation, cold dark matter hierarchically collapses into condensations called halos. Baryons follow this collapse on large scales and cool at the centers of halos to form galaxies (Rees & Ostriker, 1977; Silk, 1977; White & Rees, 1978). There is a tight connection between the masses and evolutionary histories of galaxies and their halos, as demonstrated by the success of various classes of models for the galaxy-halo connection that have been put forward over the past decades. For example, subhalo abundance matching (Kravtsov et al., 2004; Tasitsiomi et al., 2004; Vale & Ostriker, 2004; Conroy et al., 2006; Conroy & Wechsler, 2009; Moster et al., 2010; Behroozi et al., 2013a) assigns galaxies to halos based on rank orderings of stellar mass and some halo property, such as mass. Halo occupation distributions (Peacock & Smith, 2000; Seljak, 2000; Berlind & Weinberg, 2002; Cooray & Sheth, 2002) assign one or multiple galaxies to a halo based on its mass or other properties (Hearin et al., 2016). Finally, semi-analytical models (Kauffmann et al., 1993; Somerville & Primack, 1999; Bower et al., 2006; Guo et al., 2010) attempt to describe the sophisticated mechanisms of galaxy formation but are ultimately based on merger trees that represent the evolutionary histories of halos. Even in hydrodynamic simulations designed to follow the formation of galaxies from first principles, certain parameters are sometimes explicitly tied to the halo mass or radius — for example the seeding of black holes (Vogelsberger et al., 2013) or stellar wind velocities (Davé et al., 2016).

Thus, the models described above share one important caveat: they depend upon a particular definition of the halo boundary and mass. The most widely accepted definition is for the halo radius to enclose some overdensity $\Delta$ such that

[TABLE]

where $\rho_{\rm ref}$ is either the critical or mean matter density of the universe (e.g. Cole & Lacey, 1996). This spherical overdensity definition has a number of manifest advantages: radius and mass are trivially related via the reference density, $M_{\Delta}$ can be measured in both simulations and observations by counting the mass included in shells of increasing radius, and re-scaling halo radii by $R_{\Delta}$ leads to a self-similar form of the density profile that can approximately be described as a function of only mass and a concentration parameter (Navarro et al., 1995, 1996, 1997, 2004; Burkert, 1995; Cole & Lacey, 1996).

By contrast, spherical overdensity radii and masses suffer from a number of issues. First, the extent to which the profiles are self-similar at different masses and redshifts depends on the somewhat arbitrarily chosen overdensity threshold (Diemer & Kravtsov, 2014, 2015). One can derive a so-called virial overdensity of $\Delta_{\rm vir}=178$ from arguments based on the collapse of an isolated top-hat overdensity in an $\Omega_{\rm m}=1$ universe (Gunn & Gott, 1972; Peebles, 1980; Lacey & Cole, 1993), where the overdensity evolves with time in $\Lambda$CDM cosmologies (e.g., Lahav et al., 1991). However, the peaks in the initial Gaussian random field are not in a top-hat shape (Dalal et al., 2008, 2010), the particles do not instantaneously virialize as assumed in the model (e.g. Shaw et al., 2006; Sánchez-Conde et al., 2007; Ludlow et al., 2012), and halos do not form in isolation, creating complicated density fields that extend well past $R_{\rm vir}$ (Prada et al., 2006; Hayashi & White, 2008; Oguri & Hamana, 2011; Diemer & Kravtsov, 2014). One manifestation of this extended structure is that subhalos falling into a more massive host begin to lose mass long before they cross the host’s $R_{\rm vir}$ (Behroozi et al., 2014; Peñarrubia & Fattahi, 2017). Finally, spherical overdensity masses can grow unphysically despite a constant halo density profile because the reference density decreases with cosmic time, an effect called pseudo-evolution (Diemand et al., 2005; Cuesta et al., 2008; Diemer et al., 2013b; Zemp, 2014; More et al., 2015).

In order to mitigate these issues, a number of alternative mass definitions have been put forward. The most popular of these is the friends-of-friends (FOF) mass (Davis et al., 1985; Jenkins et al., 2001). Although appealingly simple, this algorithm relies on a somewhat arbitrarily chosen linking length parameter, and, for common choices of this parameter, FOF groups can include neighboring halos (White, 2001). Furthermore, FOF masses have been shown to suffer from dependencies on mass resolution and halo concentration (More et al., 2011; Benson, 2017). Another alternative was suggested by Cuesta et al. (2008), who argued for the radius where the average radial velocity changes from outflowing to infalling. This radius, however, is not clearly defined in some low-mass halos and encloses a large amount of matter falling toward the halo for the first time that arguably should not be included (Diemer & Kravtsov, 2014). Anderhalden & Diemand (2011) suggested counting all particles that ever entered the halo, a definition that suffers from similar theoretical issues. The ORIGAMI algorithm (Falck et al., 2012; Neyrinck, 2012) defines halos by identifying particles that have switched positions with other particles along three orthogonal axes. Theoretical considerations aside, the most important issue with all of these mass definitions is that they cannot be measured in the real universe.

Recently, it has been argued that a more natural halo boundary is provided by the splashback radius, $R_{\rm sp}$, the radius where particles reach the apocenter of their first orbit after infall (Diemer & Kravtsov, 2014; Adhikari et al., 2014; More et al., 2015; Mansfield et al., 2017). The theoretical inspiration for this definition is provided by the spherical collapse model, where spherically symmetric shells of matter successively fall onto an initial power-law density perturbation, creating a power-law inner density profile (Fillmore & Goldreich, 1984; Bertschinger, 1985; Mohayaee & Shandarin, 2006; Ascasibar et al., 2007; Diemand & Kuhlen, 2008; Vogelsberger et al., 2011; Lithwick & Dalal, 2011; Adhikari et al., 2014; Shi, 2016). Particles at the apocenter of their first orbit pile up due to their low radial velocity, creating a caustic that manifests itself as a sharp drop in the density profile. This so-called splashback radius represents a clear boundary between matter orbiting in the halo and matter on a first infall toward the halo.

The sharp drop in stacked halo density profiles at the splashback radius was recently detected in cosmological simulations (Diemer & Kravtsov, 2014), and its location was shown to primarily depend on mass accretion rate. Adhikari et al. (2014) reproduced this dependence with a simple theoretical model, making a convincing case for the connection between the splashback radius and the density drop (see also Shi, 2016). More et al. (2015) adopted the definition of $R_{\rm sp}$ as the radius where the density profile reaches its steepest slope and investigated its dependence on mass accretion rate and redshift. Finally, More et al. (2016, see also and ) detected a sharp drop in the stacked density profiles of galaxy cluster members at the splashback radius, though such observations are complicated by the systematics of the cluster identification method (Zu et al., 2016; Busch & White, 2017, see also Rines et al. 2013, Tully 2015, Patej & Loeb 2016, and Umetsu & Diemer 2017 for hints of the splashback radius in observations of individual clusters and weak lensing signals).

All of the theoretical and observational work discussed above has been based on the definition of $R_{\rm sp}$ as the radius where the logarithmic slope of the density profile is steepest. While this definition is intuitive, the radius of the steepest slope is affected by a trade-off between the sharply falling inner profile and the outer infall region. Thus, it is not clear what fraction of particles actually reach their orbital apocenter inside the radius of the steepest slope and whether this fraction is universal across halo masses, redshifts, and cosmologies. Moreover, Mansfield et al. (2017) showed that substructure can wipe out the signature of splashback in simulated density profiles and leads to a significant bias in the splashback radius measured from stacked density profiles. Finally, the scatter in the $R_{\rm sp}$ distribution cannot be determined from stacked density profiles.

For all of these reasons, it is desirable to measure $R_{\rm sp}$ in individual simulated halos using a method that does not rely on spherically averaged density profiles. Mansfield et al. (2017) performed such measurements using the full three-dimensional density information to obtain nonspherical splashback shells. While their method relies only on the density field at a given time, it demands relatively well-resolved halos with more than \num50000 particles and can fail for the slowest accreting fraction of halos (Mansfield et al., 2017). In order to measure $R_{\rm sp}$ in less well-resolved systems, Diemer (2017, hereafter Paper I) suggested an algorithm based on the apocenter passages of individual particles. This method necessarily uses all of the snapshots of a simulation but was shown to converge for halos resolved by as few as $1000$ particles.

In this second paper in the series, we investigate the relation between the $R_{\rm sp}$ and $M_{\rm sp}$ of individual halos and their spherical overdensity mass, accretion rate, redshift, and cosmology. In order to facilitate the use of $R_{\rm sp}$ as a practical definition of the halo boundary, we provide accurate fitting functions for these dependencies. While the particle apocenters are not directly observable, we discuss the connection of our new $R_{\rm sp}$ measurements to results based on stacked density profiles.

The paper is organized as follows. In Section 2, we briefly summarize our simulations and algorithm, referring the reader to Paper I for details. We show our results in Section 3, and compare them to previous work and theoretical models in Section 4. We further discuss the implications of our results in Section 5 and summarize our conclusions in Section 6.

Throughout the paper, we adopt the same symbols as in Paper I which are summarized in Table 1. Any input quantities to the fitting functions for $R_{\rm sp}$, such as peak height and mass accretion rate, are defined in terms of conventional masses and radii that can be measured by a standard halo finder. While theoretical models typically refer to the instantaneous accretion rate $s\equiv d\log(M)/d\log(a)$, this quantity cannot be measured in simulation data due to the noisy nature of mass accretion histories. Thus, Diemer & Kravtsov (2014) defined the mass accretion rate over a finite range of time,

[TABLE]

where $M=M_{\rm vir}$ and the $a_{0}$-$a_{1}$ pairs were chosen manually to correspond to roughly a crossing time (see also Lau et al., 2015; More et al., 2015; Mansfield et al., 2017). As in Paper I, we instead choose $M=M_{\rm 200m}$ and measure the accretion rate over one dynamical time, $a_{1}\equiv a(t-t_{\rm dyn})$. The dynamical time used in this definition depends only on the chosen overdensity and cosmology, not on the properties of individual halos (Paper I). We emphasize that one has to be careful when interpreting mass accretion rates in terms of the growth of the physical density profile. Pseudo-evolution, i.e. spurious growth due to the changing definition of the overdensity with redshift, contributes to changes in $M_{\rm 200m}$. For example, for a static, non-evolving Navarro-Frenk-White (NFW) profile, the “accretion rate” is $\Gamma\approx 0.5$ (regardless of redshift or halo mass). Thus, halos with $\Gamma<0.5$ are maintaining the same physical mass profile within $R_{\rm 200m}$, or are even losing mass (e.g. due to tidal disruption as they approach another halo).

2. Methods

In this section, we describe the $N$-body simulations used for this project and give a brief summary of the Sparta algorithm, referring the reader to Paper I for details.

2.1. $N$-Body Simulations

Our results are based on a suite of dissipationless $\Lambda$CDM simulations of different box sizes and cosmologies (Table 2). Our fiducial cosmology is the same as that of the Bolshoi simulation (Klypin et al., 2011), but we also use simulations of the $Planck$ cosmology in order to investigate the cosmology dependence of the splashback radius (see Table 3). The initial conditions for the simulations were generated using the second-order Lagrangian perturbation theory code 2LPTic (Crocce et al., 2006). The simulations were started at redshift $z=49$, sufficiently high to avoid transient effects (Crocce et al., 2006), and were run with the publicly available code Gadget2 (Springel, 2005). In Paper I, we showed that the number of saved snapshots (generally $100$) is sufficient for the Sparta algorithm to give reliable results.

We used the phase-space halo finder Rockstar (Behroozi et al., 2013b) to extract halos and subhalos from the snapshots of each simulation and the Consistent-Trees code (Behroozi et al., 2013c) to establish subhalo relations and assemble merger trees. We note that the halo catalogs and merger trees used in this paper are based on $R_{\rm 200m}$ as the halo radius. This definition matters because halos whose centers lie inside $R_{\rm 200m}$ of another, larger halo are considered subhalos and are treated rather differently (Section 2.2). Rockstar computes $R_{\rm 200m}$ using only bound particles in order to avoid spurious contributions from their hosts. While the merger trees are based on these bound-only radii, we generally use $R_{\rm 200m}$ as computed from all particles, bound and unbound, and explicitly state when we are using bound-only masses and radii. For the vast majority of host halos, the difference between the two masses is small.

2.2. The Sparta algorithm

In each host halo, we track all particles as they fall into the halo for the first time and record whether a particle was part of a subhalo at infall. Thereafter, we follow the particle’s trajectory and, at the apocenter of its first orbit, record the time $t_{\rm sp}$, the splashback radius $r_{\rm sp}$, and the enclosed mass $m_{\rm sp}$. We exclude particles that were part of a subhalo larger than $0.01$ times the host mass at infall, because dynamical friction biases the $r_{\rm sp}$ of such particles. From the remaining distribution, we compute various estimators of $R_{\rm sp}$ and $M_{\rm sp}$, namely, the mean, median, and higher percentiles of the distribution (Paper I).

For the results presented in this paper, Sparta analyzed between $38$ and $640$ million particle apocenter passages per simulation, a total of $4.4$ billion splashbacks. Figure 1 shows a visualization of the conventional virial and splashback radii of halos with $N_{\rm 200m}\geq 1000$ particles (corresponding to $M_{\rm 200m}\geq 1.4\times 10^{11}\>h^{-1}M_{\odot}$) in a $30\>h^{-1}{\rm{Mpc}}$ slice through the L0125 simulation. The density field is visualized using the gotetra code (Mansfield et al. 2017, in preparation), which is based on a tetrahedron density estimator (Abel et al., 2012; Hahn et al., 2013, 2015). The splashback radii shown correspond to $R_{\rm sp}^{87\%}$, i.e., the radius enclosing $87\%$ of the particle apocenters (Table 1), the definition that most closely matches the results of Shellfish (Paper I). Generally speaking, $R_{\rm sp}$ is significantly larger than $R_{\rm vir}$. A few halos were not assigned a splashback radius because they had recently been subhalos, but this fraction is relatively small (about $5\%$; see Paper I).

3. Results

In this section, we analyze the distribution of $R_{\rm sp}$, $M_{\rm sp}$, and $\Delta_{\rm sp}$ as a function of halo mass, accretion rate, redshift, and cosmology. As shown in previous work (Diemer & Kravtsov, 2014; More et al., 2015; Mansfield et al., 2017), the parameter that has the strongest influence on $R_{\rm sp}/R_{\rm 200m}$ is the mass accretion rate. Thus, the majority of the section focuses on the $\Gamma$–$R_{\rm sp}$ relation. However, we also discuss the distribution of $R_{\rm sp}$ marginalized over $\Gamma$, partly because the accretion rate of individual halos is difficult to measure observationally.

3.1. Halo Sample

As shown in Paper I, it does not matter which simulation a halo originated from because our results are insensitive to mass resolution as long as $N_{\rm 200m}\geq 1000$, a limit that is applied to all halo samples hereafter. We combine all halos with valid $R_{\rm sp}$ and $M_{\rm sp}$ measurements into samples that are distinguished only by their redshift and cosmology. In order to compute $\Gamma$, we require halo masses at the current snapshot and at a particular time in the past. We exclude any halos that were not host halos at the current or past snapshot but include halos that temporarily became subhalos at intermediate times (so-called backsplash halos). We confirmed that excluding these halos makes a negligible difference to our results. We note that virtually all halos without a valid $R_{\rm sp}$ measurement had recently been subhalos and might thus be excluded anyway.

Finally, we exclude the most extreme mass accretion rates from consideration. As discussed in Paper I, the lowest values of $\Gamma$ (in particular negative values) correspond to halos that are being disrupted because they are falling into or passing close by another halo. Due to the resulting tidal disruption, their radius and mass undergo drastic changes and are not particularly well defined, regardless of whether conventional definitions or $R_{\rm sp}$ are used. Similarly, some halos are assigned very large values of $\Gamma$ that are indicative of a merger or disruption event. Thus, we exclude halos with $\Gamma<0$ or $\Gamma>12$ from our samples and do not include them when deriving our fitting function. This cut affects less than 1% of halos at $z=0$, and about 2% at higher redshifts.

After all cuts, the sample for the fiducial cosmology includes about \num250000 halos at $z=0$, about \num150000 at $z=1$, and about $3500$ at $z=8$. The $Planck$ sample contains about \num170000 halos at $z=0$ and about \num120000 at $z=1$. Unless stated otherwise, we plot the median $R_{\rm sp}$ and $M_{\rm sp}$ of a halo sample because the mean is more sensitive to outliers. We compute the statistical uncertainty in each bin from the standard deviation and omit bins with fewer than 30 halos.

3.2. Distribution and Scatter

We begin by analyzing the distributions of $R_{\rm sp}$ and $M_{\rm sp}$ at a fixed mass accretion rate, mass, redshift, and cosmology. Figure 2 shows examples of the distribution of residuals for three definitions of the splashback radius, mass, and overdensity, namely, $R_{\rm sp}^{\rm mn}$, $R_{\rm sp}^{50\%}$, and $R_{\rm sp}^{87\%}$ (corresponding to the mean, median, and 87th percentile of the particle apocenter distribution). Generally, the distributions of $R_{\rm sp}$ and $M_{\rm sp}$ are reasonably well described by log-normal functions, except for a tail toward positive values. The tails are stronger in $M_{\rm sp}/M_{\rm 200m}$ than in $R_{\rm sp}/R_{\rm 200m}$, presumably because the splashback mass can increase due to large subhalos that have crossed into the halo but not yet influenced $R_{\rm sp}$. The tails are weakest for $R_{\rm sp}^{50\%}$ and $M_{\rm sp}^{50\%}$, and increase toward higher percentiles. For example, the $M_{\rm sp}/M_{\rm 200m}$ distribution of the 87th percentile (orange lines in Figure 2) has a peak that is slightly shifted off the median value. The distribution of the enclosed overdensity $\Delta_{\rm sp}$ is much wider due to the combined scatter from $R_{\rm sp}$ and $M_{\rm sp}$, but shows no systematically discernible tails.

As the residuals from the median values are nearly log-normal, we will hereafter quantify the distributions as the median $R_{\rm sp}$ or $M_{\rm sp}$ and the logarithmic 68% scatter in dex. Figure 2 hints at some of the most important trends: the scatter is smallest for low percentiles, low $\Gamma$, and large halo masses. In contrast, redshift does not have a major impact on the scatter (not shown in Figure 2). We find that the scatter, expressed in units of dex, can be approximated as

[TABLE]

where $p$ is the percentile divided by 100, and $\sigma_{\rm p}$ is zero for $R_{\rm sp}^{\rm mn}$ and $M_{\rm sp}^{\rm mn}$. The parameters differ slightly for $R_{\rm sp}$ and $M_{\rm sp}$, and are given in Table 4. They were derived from a least-squares fit to the measured scatter in the $\Gamma$–$R_{\rm sp}$ relation of the fiducial and $Planck$ samples at redshifts $0.2$, $0.5$, $1$, $2$, $4$, and $8$ and in the peak-height bins shown in Figure 3 (we ignore the scatter at $z=0$ which is artificially increased; see Paper I). The scatter in the enclosed overdensity $\Delta_{\rm sp}$ is well approximated by the scatter in $R_{\rm sp}$ and $M_{\rm sp}$ added in quadrature,

[TABLE]

For example, the scatter at intermediate masses ($\nu=1$) and accretion rates ($\Gamma=1$) is about $0.045$ dex in both $R_{\rm sp}^{\rm mn}$ and $M_{\rm sp}^{\rm mn}$, and increases to about $0.055$ dex for the 87th percentile. The lowest scatter of about $0.02$ dex occurs at $\Gamma\approx 0.5$ and $\nu\approx 3$. We note that Equation (3) extrapolates to lower (and even negative) scatter but should not be taken seriously below $\sigma=0.02$. The highest scatter occurs at low masses ($\nu=0.5$) and high accretion rates ($\Gamma=10$), about $0.08$ dex for $R_{\rm sp}^{\rm mn}$ and $0.1$ dex for $R_{\rm sp}^{87\%}$, resulting in a scatter of about $0.2$ dex in $\Delta_{\rm sp}^{87\%}$.

We note that Equation (3) does not describe the scatter at $z=0$, or, more generally, at the final redshift of a simulation. At those snapshots, the scatter is increased significantly by the correction term introduced to balance the asymmetric time distribution of particle splashbacks (Paper I). This term de-biases the results, on average, but induces additional scatter that strongly depends on $\Gamma$ because the extrapolation in time is less reliable for rapidly evolving halos. In particular, the scatter is barely increased at low accretion rates (\Gamma\lower 3.01385pt\hbox{;\buildrel<\over{\sim};}1) but increased by up to a factor of $2$ at high accretion rates. Finally, we caution that (due to the tails in the distributions) the 2$\sigma$ (i.e., 95%) scatter can be slightly larger than twice the 1$\sigma$ (i.e. 68%) scatter. The difference exhibits a rather complex dependence on mass and redshift, and we refrain from adding further complexity to our fitting function.

3.3. Fitting Function

Before we discuss the various dependencies of $R_{\rm sp}$ and $M_{\rm sp}$ in detail, we summarize our results with a convenient fitting function. We find that the $\Gamma$–$R_{\rm sp}/R_{\rm 200m}$ and $\Gamma$–$M_{\rm sp}/M_{\rm 200m}$ relations are — at any redshift, cosmology, peak height, and for any $R_{\rm sp}$ definition — well fit by an expression similar to those suggested by Diemer & Kravtsov (2014) and More et al. (2015),

[TABLE]

where $X_{\rm sp}$ can stand for either $R_{\rm sp}$ or $M_{\rm sp}$, and $A$, $B$, and $C$ are free parameters. Those parameters are, in turn, functions of mass and redshift such that

[TABLE]

where we have introduced a total of nine free parameters. However, not all of these parameters are necessary to fit either $R_{\rm sp}$ or $M_{\rm sp}$. In principle, there is no reason to expect that $R_{\rm sp}$ and $M_{\rm sp}$ should be fit by exactly the same functional form. Thus, it is not surprising that slightly different parameters are used in the two fits. Furthermore, the fit parameters in Equation (3.3) depend on the definition of $R_{\rm sp}$. We fit $R_{\rm sp}^{\rm mn}$ and $M_{\rm sp}^{\rm mn}$ separately and opt to further parameterize the dependence of the fit parameters on the percentile. For this purpose, we introduce $p$, the percentile value divided by $100$ (e.g., $0.5$ for the median). The fit parameters depend on $p$ in a nontrivial manner, where

[TABLE]

We constrain all free parameters simultaneously using a Levenberg–Marquart least-squares fit to the median $\Gamma$–$R_{\rm sp}$ relation at redshifts [math], $0.2$, $0.5$, $1$, $2$, $4$, and $8$, at the same peak-height bins as those shown in Figure 3, and in both the fiducial and $Planck$ cosmologies. The statistical uncertainty is used as an inverse weight in the fit. However, in this scheme, low $\Gamma$ values are weighted much more heavily than high $\Gamma$ values where the halo sample is less populated, and the $\chi^{2}$ values are much greater than one, indicating that the error bars are underestimated. Thus, we add a $1\%$ systematic error in quadrature with the statistical error, which balances the weights and leads to more reasonable $\chi^{2}$ values between $1$ and $2$. We constrain the dependence of the parameters on $p$ by simultaneously fitting the 50th, 63rd, 75th, and 87th percentiles. Based on the results of Paper I, the highest percentiles are known to be unreliable, and we thus do not attempt to fit their $\Gamma$–$R_{\rm sp}$ relation. The fitting function should not be extrapolated beyond the 50th and 87th percentiles or beyond the range of mass accretion rates for which it was constrained ($0<\Gamma<12$).

Figure 3 shows a summary of our main results and fitting function. Each column shows $R_{\rm sp}$, $M_{\rm sp}$, and $\Delta_{\rm sp}$ for a different redshift, and the colors indicate different halo masses. Given the statistical uncertainties, the fits agree with the median relations to 5% or better everywhere in $\Gamma$–$\Omega_{\rm m}$–$\nu$ parameter space where data are available. Figure 3 shows $R_{\rm sp}^{75\%}$, but the same holds for $R_{\rm sp}^{\rm mn}$ and up to the 87th percentile. Similarly, Figure 3 shows the results for our fiducial cosmology, but the $Planck$ results are fit equally well. We note that the relations for the median $R_{\rm sp}$ and $M_{\rm sp}$ do not necessarily have to predict the correct median enclosed overdensity $\Delta_{\rm sp}$ if the distributions are asymmetric. However, combining the fitting functions for $R_{\rm sp}$ and $M_{\rm sp}$, we find that the result has roughly the expected agreement (bottom row of Figure 3; any error on $R_{\rm sp}$ is tripled; see Equation (4)). We have implemented our fitting function in the publicly available Python code Colossus 111Colossus is a Python module for computations related to cosmology, large-scale structure, and dark matter halos (Diemer, 2015; Diemer & Kravtsov, 2015). In addition to all fitting functions proposed in this paper, we have also implemented the More et al. (2015) fit as well as the Adhikari et al. (2014) and Shi (2016) semi-analytical models (see Section 4)..

Our fitting function highlights a number of interesting features in the data. First, More et al. (2015) found no significant evidence for a dependence of $R_{\rm sp}$ on halo mass, largely due to the limited accuracy of the $R_{\rm sp}$ determination from stacked density profiles. Here, we qualitatively confirm the results of Mansfield et al. (2017) in that we find a slight but significant dependence on mass (or, equivalently, $\nu$). Interestingly, this mass dependence is weak for $R_{\rm sp}^{\rm mn}$ and $R_{\rm sp}^{50\%}$, but becomes more significant at higher percentiles. We discuss the mass dependence further in Section 5.2.

According to our fitting function, $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ approach constant values at high $\Gamma$ that do not depend on redshift, mass, or cosmology (they do, however, depend on the percentile definition). We caution that our data do not unambiguously require such behavior. But the data also do not show a statistically significant preference for an evolution of the high-$\Gamma$ value with mass or redshift. The fits for $R_{\rm sp}$ and $M_{\rm sp}$ imply that $\Delta_{\rm sp}$ also asymptotes to a particular value at high $\Gamma$, about $500$ for $R_{\rm sp}^{\rm mn}$. At $z=0$ and for our fiducial cosmology, $\Delta_{\rm m}=500$ corresponds to $\Delta_{\rm c}=135$, meaning that the average $R_{\rm sp}$ (using any definition) is always larger than $R_{\rm 200c}$, even at very high accretion rates. At higher redshift, however, $\Delta_{\rm c}$ becomes similar to $\Delta_{\rm m}$, meaning that $R_{\rm sp}$ can reach radii smaller than $R_{\rm 200c}$.

3.4. Dependence on Redshift and Cosmology

Another noteworthy feature of the fitting function presented in Equations (5–3.3) is that it encapsulates any dependence on redshift and cosmology into a dependence on $\Omega_{\rm m}(z)$. Figure 4 demonstrates that such a scaling is strongly suggested by the data. Instead of comparing different peak heights within each panel as in Figure 3, the columns show different peak heights, and the lines in each panel correspond to different redshifts. At redshifts higher than $z\approx 2$, $\Omega_{\rm m}$ barely changes, which manifests itself as a constant $\Gamma$–$R_{\rm sp}$ relation, even at $z=8$. We find similar dependencies for other peak-height bins, as well as for higher percentiles.

Although similar redshift scalings were predicted in analytical models (Adhikari et al., 2014; Shi, 2016) and seemed to work well for the fitting functions of More et al. (2015) and Mansfield et al. (2017), the dependence on cosmology has not been explicitly tested in simulations. Figure 5 compares the $\Gamma$–$R_{\rm sp}$ relation for a particular peak-height bin in a number of cosmologies. First, we focus on the fiducial and $Planck$ cosmologies (dark and light blue lines, respectively). As expected, the $Planck$ cosmology produces slightly higher values of $R_{\rm sp}/R_{\rm 200m}$ at low redshifts due to its higher $\Omega_{\rm m,0}=0.32$ (compared to $0.27$ in the fiducial cosmology). At high redshift, the difference vanishes, indicating that a scaling with $\Omega_{\rm m}$ captures the difference between the cosmologies. We note that the $Planck$ cosmology also has a $\sigma_{8}$ that is 2% higher than that in the fiducial cosmology, but our data are not constraining enough to definitively exclude dependencies on $\sigma_{8}$ or other cosmological parameters.

We further test the cosmology dependence using simulations that have $\Omega_{\rm m}=1$ throughout, namely, self-similar cosmologies with a power-law power spectrum (Table 3). In such universes, we expect the $\Gamma$–$R_{\rm sp}$ relation not to evolve with redshift at all. We confirm this self-similar scaling, which constitutes further evidence that $R_{\rm sp}$ depends on $\Omega_{\rm m}$ rather than $z$. Furthermore, the self-similarity allows us to combine halos at different redshifts into one sample per simulation, ensuring coverage over a wide range of peak heights. The samples from each self-similar simulation differ only by the slope of the initial power spectrum, $n$ (see Diemer & Kravtsov 2015 for more details on the simulations and the technique of combining redshifts). Figure 5 shows the $\Gamma$–$R_{\rm sp}$ relation in two self-similar cosmologies, namely, those with $n=-2.5$ (purple) and $n=-1$ (red). A slope of $-2.5$ is close to the slope on scales relevant for halo formation in a $\Lambda$CDM cosmology (e.g., Figure 4 in Diemer & Kravtsov, 2015). Our fitting function (with $\Omega_{\rm m}=1$) matches the $\Gamma$–$R_{\rm sp}$ relation from this simulation well, even though the self-similar models were not used when constraining the fit parameters. Nevertheless, the self-similar cosmology with $n=-2.5$ exhibits essentially the same relations as our fiducial cosmology at high redshift. This match is yet another confirmation that $\Omega_{\rm m}$ is the variable that controls $R_{\rm sp}$, not redshift.

By contrast, the self-similar simulation with a much shallower power spectrum slope, $n=-1$, leads to a rather different $\Gamma$–$R_{\rm sp}$ relation that is not described by our fitting function. We conclude that $\Omega_{\rm m}$ is not the only parameter than influences $R_{\rm sp}$; the power-spectrum slope clearly has an impact too. In principle, we could introduce $n$ as an extra parameter into our fitting function and use the self-similar models to constrain the dependence of $R_{\rm sp}$. However, the impact of $n$ in $\Lambda$CDM is degenerate with the effects of $\Omega_{\rm m}$ and mass, making it difficult to disentangle the dependencies. Furthermore, it is not clear a priori how to define $n$ in a $\Lambda$CDM cosmology where the slope is scale-dependent (and thus halo mass–dependent). As the self-similar models have little practical application, we leave an investigation of the effect of $n$ for future work.

3.5. Definitions of $R_{\rm sp}$

Besides $\nu$ and $\Omega_{\rm m}$, our fitting function depends on the definition of $R_{\rm sp}$ and $M_{\rm sp}$, i.e. whether we use the mean of the $r_{\rm sp}$ and $m_{\rm sp}$ distributions, their median, or higher percentiles. In Section 5, we will demonstrate that different definitions can be useful for different purposes. Figure 6 shows the $\Gamma$–$R_{\rm sp}$ relation for different definitions of $R_{\rm sp}$. At first sight, it appears that the main difference is the normalization of the curves. However, there are also nontrivial changes in the shape of the relations, demanding the relatively large number of free parameters introduced in Equation (3.3). At the highest percentiles (greater than the 87th), the evolution of the normalization becomes superlinear, and the shape changes in a complex manner with percentile. As the measurement of the highest percentiles is relatively uncertain anyway (Paper I), we limit the applicability of our fitting function to the range between the 50th and 87th percentiles and omit it for the 99th percentile in Figure 6.

3.6. Constraints on $R_{\rm sp}$ without Knowledge of the Mass Accretion Rate

So far, we have considered $R_{\rm sp}$ primarily as a function of $\Gamma$ and secondarily of other variables because $\Gamma$ has the strongest effect. Unfortunately, $\Gamma$ is also the variable that is hardest to measure observationally: in practice, we almost never know the mass accretion rate of individual halos. Thus, we also give expressions for $R_{\rm sp}$ in the absence of any knowledge about $\Gamma$, i.e. the median $R_{\rm sp}$ as a function of halo mass and redshift but marginalized over all mass accretion rates.

First, it is instructive to consider the average $\Gamma$ as a function of halo mass and redshift, as shown in Figure 7. The accretion rate increases with peak height (because larger halos are more likely to be actively forming at any given time in hierarchical structure formation), but there is also a strong trend with redshift, with much higher $\Gamma$ at high redshift. Moreover, the distribution at fixed mass and redshift is broad and not particularly well described by a normal or log-normal distribution (partially because a few percent of halos have negative accretion rates, especially at low $\nu$). Neglecting a mild dependence on redshift, we find that the logarithmic scatter in the $\Gamma$ distributions is roughly

[TABLE]

where $\sigma_{\Gamma}$ is in units of dex. Despite the large scatter, there are clear trends in the median accretion rate $\Gamma(\nu,z)$ that can be approximated with the expression

[TABLE]

where

[TABLE]

This fitting function is shown with dashed lines in Figure 7, and fits the median $\Gamma$ to better than 5% at all $\nu$ and $z$ and for both the fiducial and $Planck$ cosmologies. It is clear that a dependence on $\Omega_{\rm m}$ instead of $z$ would not work in this case, as $\Gamma$ strongly increases at high redshift even though $\Omega_{\rm m}\approx 1$ (e.g., Zhao et al., 2009).

Given the trends in $\Gamma(\nu,z)$, we expect the $\nu$–$R_{\rm sp}$ relation to experience a conflation of multiple competing effects: the $\Gamma$–$R_{\rm sp}$ relation increases with redshift due to an increasing $\Omega_{\rm m}$, but the increasing $\Gamma$ at high $z$ leads to lower $R_{\rm sp}$. As the distribution of $\Gamma$ is nonsymmetric and the $\Gamma$–$R_{\rm sp}$ relation is nonlinear, there is no guarantee that we can construct a $\nu$–$R_{\rm sp}$ relation from our $\nu$-$\Gamma$ and $\Gamma$–$R_{\rm sp}$ relations. However, we find that such a procedure does, in fact, work surprisingly well. Figure 8 shows the $\nu$–$R_{\rm sp}$ relation for a number of redshift bins and $R_{\rm sp}$ definitions. The dashed lines show the fit obtained from our $\Gamma$-based fitting function (Equations (5–3.3)), evaluated using the $\Gamma$ fit of Equation (3.6). The fits are accurate to 5% for $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ at all peak heights, redshifts, and $R_{\rm sp}$ definitions and for both the fiducial and $Planck$ cosmologies. As expected, the corresponding maximum error in $\Delta_{\rm sp}$ is about 15%.

Naturally, the scatter in the $\nu$–$R_{\rm sp}$ relation is increased compared to that of the $\Gamma$–$R_{\rm sp}$ relation because we are averaging over all mass accretion rates. In particular, the 68% scatter is about $0.07$ dex in $R_{\rm sp}/R_{\rm 200m}$ regardless of peak height or redshift and between $0.04$ dex and $0.1$ dex in $M_{\rm sp}/M_{\rm 200m}$, where the highest scatter occurs at high redshift and low peak height. The distributions in $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ combine to a more or less constant scatter of $0.15$ dex in $\Delta_{\rm sp}$. As with the $\Gamma$–$R_{\rm sp}$ relation, the distribution of $R_{\rm sp}$ and $M_{\rm sp}$ values is reasonably described by a log-normal, but the tails toward high and low values are enhanced when marginalizing over $\Gamma$. As a result, the 2$\sigma$ (i.e., 95%) scatter can be slightly larger than twice the 1$\sigma$ (i.e., 68%) scatter. The exact distribution shows complex dependencies on peak height and redshift, and we refrain from quantifying it further.

Overall, the $\sim 0.07$ dex scatter in the $\nu$–$R_{\rm sp}$ relation is surprisingly low, considering the scatter in the $\Gamma$–$R_{\rm sp}$ relation and that the distribution of accretion rates is relatively extended. The low scatter allows for meaningful inferences regarding $R_{\rm sp}$ and $M_{\rm sp}$ in the absence of any knowledge about $\Gamma$.

4. Comparison to Previous Work

The Sparta algorithm operates based on an entirely different principle than any of the previous estimates of $R_{\rm sp}$ used in either simulations or observations. Thus, we expect that our findings might disagree with other measurements and models. In this section, we compare our data with previous simulation results based on density profiles and the Shellfish algorithm, as well as with theoretically motivated models.

4.1. Comparison with Results Based on Density Profiles

Most work on the splashback radius thus far has been based on spherically averaged density profiles, both in simulations and in observations. As those profiles suffer from noise due to resolution and substructure effects, they cannot generally be used to measure $R_{\rm sp}$ for individual halos. Thus, More et al. (2015) used the stacked density profiles of Diemer & Kravtsov (2014) and defined $R_{\rm sp}$ as the radius where the logarithmic slope of the median profile is steepest. Based on this definition, they found that $R_{\rm sp}$ decreases as a function of $\Gamma$ and increases slightly with increasing $\Omega_{\rm m}$ — trends which we confirm in this work. Presumably owing to the relatively poor accuracy and restricted range of peak height used by More et al. (2015), they did not detect any dependence on mass at fixed $\Gamma$, as found in this investigation and in Mansfield et al. (2017).

With the data at hand, we can for the first time elucidate the relation between the radius of the steepest slope and the apocenter passages of particles. Figure 9 compares the $\Gamma$–$R_{\rm sp}$ relation as derived by Sparta to the More et al. (2015) fitting function based on stacked density profiles. We note that the mass accretion rates $\Gamma_{\rm DK14}$ were computed based on $M_{\rm vir}$ in More et al. (2015) and based on $M_{\rm 200m}$ in this work, but the difference is negligible given the accuracy of this comparison. We choose $R_{\rm sp}^{75\%}$ as this definition matches the More et al. (2015) results most closely, indicating that the steepest profile slope occurs at a radius that encloses about 75% of the particle apocenters. The More et al. (2015) function matches the overall shape and redshift evolution of the relations relatively well, particularly for high-mass halos. Choosing a lower percentile might bring the overall normalizations into better agreement at low mass, but the Sparta relations become almost mass-independent at low $\Gamma$. Thus, there is no percentile for which the low-mass relation is matched well by the fit of More et al. (2015). We further discuss the connection between the density and splashback profiles in Section 5.2.

More et al. (2015) also provided formulae for the $\nu$–$R_{\rm sp}$ and $\nu$–$M_{\rm sp}$ relations, again based on stacked density profiles. The dependence on $\nu$ was presumed to be due purely to the mass dependence of $\Gamma$. We compare their function to our results in Figure 10. As the More et al. (2015) function does not depend on redshift, it cannot match the trends found in Section 3.6, but it once again coincides more or less with the 75th percentile $R_{\rm sp}$ at high peak height. Interestingly, Figure 4 of More et al. (2015) shows a hint of the reversed redshift evolution (lowered $\nu$–$R_{\rm sp}$ relation at the highest redshifts), but the trend was not significant enough to be captured in their fitting function.

We conclude that there is no exact one-to-one match between the radius of the steepest slope and the definitions used in this paper. However, $R_{\rm sp}^{75\%}$ gives a good approximation, especially at high peak heights, where the results of More et al. (2015) were most constrained.

4.2. Comparison with Shellfish

In Paper I, we undertook a halo-by-halo comparison of the results of Sparta and Shellfish. We found that $R_{\rm sp}^{87\%}$ agrees with Shellfish to a few percent, though with about 15% scatter. The small overall difference, however, does not mean that there could not be systematic trends with mass or accretion rate. Thus, we compare the $\Gamma$–$R_{\rm sp}$ relation measured by Sparta and Shellfish in Figure 11. The relations agree reasonably well, but, driven by systematic differences at low $\Gamma$, Shellfish prefers values of $R_{\rm sp}/R_{\rm 200m}$ that do not rise as sharply at low $\Gamma$. We further discuss the physical connection between splashback shells and apocenter distributions in Section 5.1.

4.3. Comparison with Theoretical Models

Given that the spherical collapse model provides the theoretical foundation for the splashback radius as a halo boundary, it is natural to use this type of model to predict $R_{\rm sp}$ and $M_{\rm sp}$. In Figure 12, we compare the models of Adhikari et al. (2014) and Shi (2016) to our results. Both models assume spherical symmetry, meaning that $R_{\rm sp}$ is uniquely defined and that all particles reach their first apocenter exactly at $R_{\rm sp}$. The values of $r_{\rm sp}$ measured by Sparta represent a distribution around this radius that is scattered due to the complexities of realistic structure formation. Assuming that this scatter does not bias the mean of the distribution, we use $R_{\rm sp}^{\rm mn}$ as the definition for this comparison.

We caution that $\Gamma_{\rm dyn}$ is not equivalent to the instantaneous mass accretion rate $s$ used in the models, which assume $s=\mathrm{const}$ such that $M\propto a^{s}$. In this case, $\Gamma$ is equal to $s$ regardless of what interval $\Gamma$ is averaged over, but realistic halo mass accretion rates tend to decrease with time at low redshifts, meaning that an average such as $\Gamma$ likely overestimates the instantaneous accretion rate. Thus, the different definitions complicate the interpretation of the comparison shown in Figure 12.

Adhikari et al. (2014) used a spherical shell collapse model (e.g., Gunn & Gott, 1972), assuming an NFW profile for the mass inside a given radius. The concentration is set by matching the slope of the profile at $R_{\rm vir}=1/2\ r_{\rm turn-around}$ to the spherical infall prediction, resulting in concentrations that do not necessarily match those observed in cosmological simulations. Due to this definition of the concentration, the prediction cannot be extended past $s=6$. Given the NFW mass profile, Adhikari et al. (2014) numerically computed the radius at which shells reach their first apocenter. The dotted lines in Figure 12 show this prediction.222We do not show the fitting function given in Adhikari et al. (2014), but rather the results of an improved numerical calculation of $\Delta_{\rm sp}$ (S. Adhikari, private communication). We compute $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ from an NFW profile with the concentration used in the model. We have implemented both this method and the model of Shi (2016) in the Python module Colossus.

Instead of assuming a particular function for the density profile, Shi (2016) performed a self-consistent calculation of shell collapse. As in the predecessor models of Fillmore & Goldreich (1984) and Bertschinger (1985), this calculation results in a power-law density profile with a sharp drop-off at the splashback radius. The power-law slope depends on the slope of the initial perturbation, which also sets the accretion rate. The model takes $\Omega_{\rm m}\neq 1$ into account and thus makes redshift-dependent predictions. However, due to the self-similarity of the problem at fixed redshift, the model predicts no mass dependence. The black dashed lines in Figure 12 show fits to the numerical model predictions that were given for $R_{\rm sp}/R_{\rm 200m}$ and $\Delta_{\rm sp}$, and reconstructed for $M_{\rm sp}/M_{\rm 200m}$. We note that the downturn in $M_{\rm sp}/M_{\rm 200m}$ at low $\Gamma$ is present in the model but exaggerated due to the reconstruction from the fits to $R_{\rm sp}/R_{\rm 200m}$ and $\Delta_{\rm sp}$ (X. Shi, private communication). The fitting functions were constrained for mass accretion rates in the range $0.5<s<5$ (Shi, 2016).

Both models correctly predict the general trend of a decreasing $R_{\rm sp}/R_{\rm 200m}$ with mass accretion rate. In detail, however, the models do not match our results: the slope of the $\Gamma$–$R_{\rm sp}$ relation is steeper at low redshift and the normalization is higher at high redshift. The evolution of $M_{\rm sp}/M_{\rm 200m}$ with $\Gamma$ is also steeper at low redshift and entirely different at high redshift. We discuss the physical reasons for these disagreements in Section 5.4.

5. Discussion

We have analyzed our results for the splashback radius, mass, and overdensity as a function of halo mass, accretion rate, redshift, and cosmology and expressed them in a convenient fitting function. In this section, we physically interpret some of the results and discuss the theoretical implications.

5.1. On the Physical Meaning of the Apocenter Distribution

The Sparta algorithm provides a new, independent way to measure the splashback radius, adding to two previous definitions (the radius where the density profile is steepest and the nonspherical splashback shell determined by Shellfish). In the spherical collapse model, all of these definitions are equivalent, because a spherical shell reaches apocenter at a fixed time, and all particles splash back at the same time and radius. This pileup causes the sharp density drop in the model (Fillmore & Goldreich, 1984; Bertschinger, 1985).

In reality, the situation is more complicated. First, nonsphericity causes an intrinsic scatter in the apocenter radii (Adhikari et al., 2014). Second, the energy and angular momentum of particles at infall slightly influences their apocenter (Paper I). One could argue that these effects can be seen as adding scatter but not shifting the mean apocenter, and that $R_{\rm sp}^{\rm mn}$ should thus be the best definition of $R_{\rm sp}$. However, it is not a priori clear that $R_{\rm sp}^{\rm mn}$ has any signatures that are observable in the real universe, meaning we need to establish a connection to the properties of the density profile.

The density profile does not carry a unique signature of the “true” splashback radius either. While the inner profile falls steeply near the splashback radius, the density due to nonlinearly infalling shells becomes increasingly important with radius333We note that this density contribution is not well described by the so-called 2-halo term, i.e. the statistical contribution due to the clustering of halos (e.g., Smith et al., 2003; Hayashi & White, 2008; Tinker et al., 2010). This contribution begins to dominate only at much larger radii (Diemer & Kravtsov, 2014).. Thus, the location of the steepest density slope represents a trade-off between the inner and outer profiles and cannot trivially be interpreted as the splashback radius. In observations, however, this radius is the most accessible quantity, and we have shown that it is reasonably approximated by $R_{\rm sp}^{75\%}$ as measured by Sparta. This connection will be investigated in more detail in future work, ideally using hydrodynamical simulations to directly connect the apocenter distribution to observables such as the density of satellite galaxies in clusters.

The fact that the radius of the steepest slope is merely one possible definition of $R_{\rm sp}$ was illustrated by the results of Mansfield et al. (2017), who found that massive substructures bias the radius of the steepest slope by about 30% compared to measurements where substructure has been removed. As a result, they find a somewhat larger $R_{\rm sp}$ on average, even though their measurements are also based on the density field. Another effect contributing to this difference may be the nonspherical nature of their shells (which is converted to a volume-equivalent radius). We have identified $R_{\rm sp}^{87\%}$ as the best proxy for the Shellfish results (Paper I), but the two measurements differ significantly in some halos, which remains to be investigated in more detail.

In summary, there is no one definition of $R_{\rm sp}$ that clearly corresponds to the density drop in the spherical collapse model and that can be measured in both simulations and observations. In the future, we hope to establish a tighter connection between the density drop measured by Shellfish and the apocenter distribution by considering the three-dimensional distribution of apocenters rather than only their radii.

5.2. On the Relationship between Splashback and Spherical Overdensity

Perhaps the most striking feature of the data in presented in this paper, and thus the fitting function given in Equations (5)–(3.3), is their relative complexity, with significant dependencies on halo mass and $\Omega_{\rm m}$ (at fixed $\Gamma$). While the latter dependence was expected from theoretical considerations (Adhikari et al., 2014; Shi, 2016), the halo mass dependence was not (though it was recently found in simulations by Mansfield et al. 2017). These complexities raise the question of whether we expect there to be a simple relation between $R_{\rm sp}$ and conventional spherical overdensity radii. We note that $R_{\rm sp}/R_{\Delta}$ would vary more strongly if spherical overdensity radii other than $R_{\rm 200m}$ were used.

At a fixed $R_{\rm sp}$, the ratio $R_{\rm sp}/R_{\rm 200m}$ depends on the mass profile around $R_{\rm sp}$ which depends on mass, accretion rate, and redshift in a nontrivial way (e.g., Figures 3 and 10 in Diemer & Kravtsov, 2014). These dependencies are expected from the fact that concentration depends on mass, redshift, and cosmology in a complex fashion (e.g., Bullock et al., 2001). We expect a significant correlation between concentration and mass accretion rate (Wechsler et al., 2002), raising the question of whether $c$ could be substituted for $\Gamma$ in our fitting model.

Another open question is related to the dependence of $R_{\rm sp}/R_{\rm 200m}$ on cosmological parameters. We have shown that the scaling with $\Omega_{\rm m}$ works for the $WMAP$ and $Planck$ cosmologies, and is thus likely appropriate for any realistic cosmology. However, $R_{\rm sp}/R_{\rm 200m}$ varies significantly between self-similar simulations that are distinguished only by their power-spectrum slope $n$, meaning that the splashback radius is sensitive to cosmological parameters beyond $\Omega_{\rm m}$. An alternative way to frame such issues could be to ask whether we are considering the optimal variables. For example, we quantify the mass accretion rate as $\Gamma_{\rm dyn}$, but perhaps $R_{\rm sp}$ exhibits tighter correlations with other definitions that we have not yet considered (e.g. definitions based on shorter or longer timescales, or definitions relying directly on $M_{\rm sp}$ instead of $M_{\rm 200m}$). We will systematically explore the correlations with other parameters (such as concentration) in future work.

On a theoretical level, one of the most important differences between $R_{\rm sp}$ and conventional definitions is the meaning of the overdensity $\Delta$. In the context of spherical overdensity radii, $\Delta$ is the fundamental quantity that determines $R_{\Delta}$ and $M_{\Delta}$. In the splashback picture, $\Delta_{\rm sp}$ is merely a consequence of independently determined $R_{\rm sp}$ and $M_{\rm sp}$. Our results show that not only does $\Delta_{\rm sp}$ vary systematically depending on a number of variables but also that the scatter in $\Delta_{\rm sp}$ is larger than the scatter in $R_{\rm sp}$ and $M_{\rm sp}$, highlighting that there is nothing fundamental about the splashback overdensity. This difference in interpretation has a bearing on the physical interpretation of halo growth. For example, due to the constant $\Delta$, $R_{\Delta}$ can change suddenly when a massive subhalo is accreted. In contrast, $R_{\rm sp}$ changes more slowly in this case, while $M_{\rm sp}$ and thus $\Delta_{\rm sp}$ change rapidly.

In summary, we have little reason to expect a simple relation between $R_{\rm sp}$ and $R_{\Delta}$. While we employ the quantities $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ to establish a connection between $R_{\rm sp}$ and $R_{\Delta}$, $R_{\rm sp}$ is an independent definition of the halo boundary that cannot easily be expressed in terms of spherical overdensity radii.

5.3. Compatibility with Observations

The measurement of $R_{\rm sp}$ performed by More et al. (2016) and confirmed by Baxter et al. (2017) indicates a surprisingly small splashback radius, namely, $R_{\rm sp}/R_{\rm 200m}=0.837\pm 0.031$ for a cluster sample with $\nu=2.4$ at $z=0.24$. While we cannot directly measure the mass accretion rate of the clusters, Figure 3 clearly shows that the theoretically expected value is higher.

In particular, the fitting function of More et al. (2015) predicts $R_{\rm sp}/R_{\rm 200m}=1.1$ for this sample, 32% higher than the observed value. The fitting function presented in this paper predicts almost exactly the same value for $R_{\rm sp}^{75\%}$ at the given peak height and redshift. Given the scatter in the $\nu$–$R_{\rm sp}$ relation, the observed $R_{\rm sp}$ would represent a 2$\sigma$ fluctuation even for an individual halo, whereas the More et al. (2016) result was derived by stacking the density profiles of thousands of clusters.

In contrast, Busch & White (2017) showed that the observed location of the steepest slope of the density profile can be sensitive to the details of the cluster identification algorithm. In addition, the assumption of spherical symmetry can affect the inferred radius of the steepest slope in three dimensions. Thus, the significance of the disagreement between simulations and observational results will remain unclear until the cluster-finding algorithm can be tested with realistic mock catalogs.

5.4. The Status of Analytical Models

In Section 4.3, we found that none of the semi-analytical models that have been proposed to date predict our results in detail. The model of Shi (2016) corresponds to the prediction of the spherical collapse model in a $\Lambda$CDM universe. Due to the self-similarity of the setup, this model predicts a power-law inner density profile, $\rho\propto r^{\alpha}$, whereas the density profiles in cosmological simulations steepen with radius. Given that the Adhikari et al. (2014) model is based on an NFW profile instead of a power law, the relative similarity of the predictions for $M_{\rm sp}/M_{\rm 200m}$ might seem surprising. The agreement can be explained by the Adhikari et al. (2014) procedure for setting the NFW concentration: the slope of the mass profile is set to the spherical collapse model prediction of $3s/(3+s)$, even at $s>3/2$, leading to density slopes of $\alpha\to-3$ when $s\to 0$ (and thus $c\to\infty$) and $\alpha\to-1$ when $s\to 6$ (and thus $c\to 0$). These extreme values of concentration mean that the NFW profile in the Adhikari et al. (2014) model approaches a power-law shape for both low and high accretion rates.

In summary, spherically symmetric models of shell collapse are a promising class of models for predicting the splashback radius. However, these models need to be coupled with realistic halo density profiles in order to match the simulation results in detail. One possible avenue toward more accurate models would be to set the concentration of the inner profile according to a numerically calibrated concentration–mass relation, automatically introducing a dependence on halo mass that is not present in self-similar collapse models.

6. Conclusions

Using the Sparta algorithm described in Paper I, we have quantified the splashback radii and masses of a large sample of halos from $N$-body simulations of different $\Lambda$CDM cosmologies. We have investigated the dependence of those quantities on halo mass (expressed as peak height, $\nu$), accretion rate $\Gamma$, redshift (expressed as $\Omega_{\rm m}$), and cosmology. The relatively complex dependencies indicate that the splashback radius represents an independent definition of the halo boundary that cannot simply be reconstructed from conventional spherical overdensity definitions. Our main conclusions are as follows.

At fixed $\Gamma$, mass, and redshift, $R_{\rm sp}$ and $M_{\rm sp}$ are distributed roughly log-normally, with some tails toward high values. The 68% scatter in both $R_{\rm sp}$ and $M_{\rm sp}$ varies between about $0.02$ and $0.1$ dex, where scatter decreases with $\nu$ and increases with $\Gamma$. If the accretion rate is unknown and we average over all $\Gamma$, the distribution is still close to log-normal, but the scatter in the relations increases to about $0.07$ dex in $R_{\rm sp}$ and between $0.04$ and $0.1$ dex in $M_{\rm sp}$. 2. 2.

In agreement with previous work, we find that $R_{\rm sp}/R_{\rm 200m}$ and $M_{\rm sp}/M_{\rm 200m}$ decrease with accretion rate and increase with $\Omega_{\rm m}$, but we also find a significant dependence on $\nu$. We do not find any dependence on cosmological parameters (beyond the dependence on $\Omega_{\rm m}$) within different $\Lambda$CDM cosmologies. We parameterize the median $R_{\rm sp}$ and $M_{\rm sp}$ as a function of $\Gamma$, $\nu$, and $\Omega_{\rm m}$, which is accurate to 5% or better for all masses $M_{\rm 200m}>1.7\times 10^{7}\>h^{-1}M_{\odot}$, up to $z=8$, and for the $WMAP$ and $Planck$ cosmologies that span the currently favored range of cosmological parameters. This function is implemented in the publicly available Python code Colossus. 3. 3.

We give a fitting function for the accretion rate as a function of $\nu$ and $z$. Using the fitted $\Gamma$ as input to the fitting functions for $R_{\rm sp}$ and $M_{\rm sp}$, we obtain predictions for the $\nu$–$R_{\rm sp}$ and $\nu$–$M_{\rm sp}$ relations that are accurate to 5% or better. 4. 4.

We compare our results to measurements of $R_{\rm sp}$ from stacked halo density profiles and find that the radius of the steepest slope as measured by More et al. (2015) corresponds roughly to the 75th percentile of the splashback distribution of particles. Similarly, we find that the 87th percentile results in a $\Gamma$–$R_{\rm sp}$ relation that is roughly compatible with that measured using the Shellfish code of Mansfield et al. (2017). 5. 5.

We compare our results to several semi-analytical models of $R_{\rm sp}$. While these models reproduce the general trends of $R_{\rm sp}$ and $\Delta_{\rm sp}$, they do not predict any mass dependence and cannot explain our results in detail.

Rather than a definitive statement, the results presented in this paper represent a starting point in our investigation of the splashback radius as a viable alternative to conventional radius definitions. Our goal is to create self-consistent halo catalogs with $R_{\rm sp}$ measurements or estimates for all halos above a certain mass threshold and subhalo relations based on $R_{\rm sp}$. With such catalogs, a number of classical topics in structure formation can be re-visited, namely, semi-analytical models of galaxy formation, assembly bias, or the growth of halos and its connection to concentration. Furthermore, our theoretical understanding of $R_{\rm sp}$ and $M_{\rm sp}$ is still lacking, as an accurate analytical description of our results from first principles remains elusive.

We are grateful to Susmita Adhikari, Neal Dalal, and Xun Shi for enlightening discussions. BD gratefully acknowledges the financial support of an Institute for Theory and Computation Fellowship. This work was completed using the Midway computing cluster provided by the University of Chicago Research Computing Center.