The Splashback Radius of Halos from Particle Dynamics. II. Dependence on Mass, Accretion Rate, Redshift, and Cosmology

Our results are based on a suite of dissipationless Λ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.

Table 2.  N-Body Simulations

Name	L	N3	${m}_{{\rm{p}}}$	$\epsilon$	$\epsilon /(L/N)$	${z}_{\mathrm{initial}}$	${z}_{\mathrm{final}}$	${N}_{\mathrm{snaps}}$	${z}_{{\rm{f}} \mbox{-} \mathrm{snap}}$	Cosmology	Reference
L2000	2000	10243	5.6 × 1011	65	1/30	49	0	100	20	WMAP (Bolshoi)	DK15
L1000	1000	10243	7.0 × 1010	33	1/30	49	0	100	20	WMAP (Bolshoi)	DKM13
L0500	500	10243	8.7 × 109	14	1/35	49	0	100	20	WMAP (Bolshoi)	DK14
L0250	250	10243	1.1 × 109	5.8	1/42	49	0	100	20	WMAP (Bolshoi)	DK14
L0125	125	10243	1.4 × 108	2.4	1/51	49	0	100	20	WMAP (Bolshoi)	DK14
L0063	62.5	10243	1.7 × 107	1.0	1/60	49	0	100	20	WMAP (Bolshoi)	DK14
L0031	31.25	10243	2.1 × 106	0.25	1/122	49	2	64	20	WMAP (Bolshoi)	DK15
L0500-Planck	500	10243	1.0 × 1010	14	1/35	49	0	100	20	Planck	DK15
L0250-Planck	250	10243	1.3 × 109	5.8	1/42	49	0	100	20	Planck	DK15
L0125-Planck	125	10243	1.6 × 108	2.4	1/51	49	0	100	20	Planck	DK15
L0100-PL-1.0	100	10243	2.6 × 108	0.5	1/195	119	2	64	20	Self-similar, n = −1.0	DK15
L0100-PL-2.5	100	10243	2.6 × 108	1.0	1/98	49	0	100	20	Self-similar, n = −2.5	DK15

Note. The N-body simulations used in this paper. Here L denotes the box size in comoving ${h}^{-1}\,\mathrm{Mpc}$, N³ is the number of particles, ${m}_{{\rm{p}}}$ is the particle mass in ${h}^{-1}\,{M}_{\odot }$, $\epsilon$ is the force-softening length in physical ${h}^{-1}\,\mathrm{kpc}$, ${z}_{\mathrm{initial}}$ and ${z}_{\mathrm{final}}$ are the redshift ranges of the simulations, ${N}_{\mathrm{snaps}}$ is the number of snapshots written to disk, and ${z}_{{\rm{f}} \mbox{-} \mathrm{snap}}$ is the redshift of the first snapshot. The references correspond to Diemer et al. (2013a, DKM13), Diemer & Kravtsov (2014, DK14), and Diemer & Kravtsov (2015, DK15). More details on our logic for choosing force resolutions are given in DK14.

Download table as:  ASCIITypeset image

Table 3.  Cosmological Parameters

Cosmology	H0	Ωm	${{\rm{\Omega }}}_{{\rm{\Lambda }}}$	${{\rm{\Omega }}}_{{\rm{b}}}$	${{\rm{\Omega }}}_{{\rm{k}}}$	${{\rm{\Omega }}}_{\nu }$	${\sigma }_{8}$	${n}_{{\rm{s}}}$	P(k)	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(k)\propto {k}^{n}$	⋯

Note. Cosmological parameters of the N-body simulations listed in Table 2. The Bolshoi cosmology roughly corresponds to the WMAP7 cosmology of Komatsu et al. (2011). The Planck values correspond to the Planck-only best-fit values given in Table 2 of Planck Collaboration et al. (2014). Some of the parameters in both the Planck and Bolshoi cosmologies are rounded for convenience. The initial matter power spectrum for the Bolshoi and Planck cosmologies was computed using the Boltzmann code Camb (Lewis et al. 2000).

Download table as:  ASCIITypeset image

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}_{200{\rm{m}}}$ as the halo radius. This definition matters because halos whose centers lie inside ${R}_{200{\rm{m}}}$ of another, larger halo are considered subhalos and are treated rather differently (Section 2.2). Rockstar computes ${R}_{200{\rm{m}}}$ 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}_{200{\rm{m}}}$ 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.

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}_{\mathrm{sp}}$, the splashback radius r_sp, and the enclosed mass m_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_sp of such particles. From the remaining distribution, we compute various estimators of R_sp and M_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}_{200{\rm{m}}}\geqslant 1000$ particles (corresponding to ${M}_{200{\rm{m}}}\geqslant 1.4\times {10}^{11}\,{h}^{-1}\,{M}_{\odot }$) in a $30\,{h}^{-1}\,\mathrm{Mpc}$ slice through the L0125 simulation. The density field is visualized using the gotetra code (P. 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}_{\mathrm{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_sp is significantly larger than ${R}_{\mathrm{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).

Zoom In Zoom Out Reset image size

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}_{200{\rm{m}}}\geqslant 1000$, a limit that is applied to all halo samples hereafter. We combine all halos with valid R_sp and M_sp measurements into samples that are distinguished only by their redshift and cosmology. In order to compute Γ, we require halo masses at the current snapshot and 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_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 Γ (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_sp are used. Similarly, some halos are assigned very large values of Γ that are indicative of a merger or disruption event. Thus, we exclude halos with ${\rm{\Gamma }}\lt 0$ or ${\rm{\Gamma }}\gt 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 250,000 halos at z = 0, about 150,000 at z = 1, and about 3500 at z = 8. The Planck sample contains about 170,000 halos at z = 0 and about 120,000 at z = 1. Unless stated otherwise, we plot the median R_sp and M_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.

We begin by analyzing the distributions of R_sp and M_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}_{\mathrm{sp}}^{\mathrm{mn}}$, ${R}_{\mathrm{sp}}^{50 \% }$, and ${R}_{\mathrm{sp}}^{87 \% }$ (corresponding to the mean, median, and 87th percentile of the particle apocenter distribution). Generally, the distributions of R_sp and M_sp are reasonably well described by log-normal functions, except for a tail toward positive values. The tails are stronger in M_sp/M_200m than in R_sp/R_200m, presumably because the splashback mass can increase due to large subhalos that have crossed into the halo but not yet influenced R_sp. The tails are weakest for ${R}_{\mathrm{sp}}^{50 \% }$ and ${M}_{\mathrm{sp}}^{50 \% }$ and increase toward higher percentiles. For example, the M_sp/M_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 Δ_sp is much wider due to the combined scatter from R_sp and M_sp but shows no systematically discernible tails.

Zoom In Zoom Out Reset image size

As the residuals from the median values are nearly log-normal, we will hereafter quantify the distributions as the median R_sp or M_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 Γ, 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

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

where p is the percentile divided by 100 and ${\sigma }_{{\rm{p}}}$ is zero for ${R}_{\mathrm{sp}}^{\mathrm{mn}}$ and ${M}_{\mathrm{sp}}^{\mathrm{mn}}$. The parameters differ slightly for R_sp and M_sp and are given in Table 4. They were derived from a least-squares fit to the measured scatter in the Γ–R_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 Δ_sp is well approximated by the scatter in R_sp and M_sp added in quadrature,

$$\begin{eqnarray}&&{\sigma }_{{{\rm{\Delta }}}_{\mathrm{sp}}}=\sqrt{{\sigma }_{{M}_{\mathrm{sp}}}^{2}+3{\sigma }_{{R}_{\mathrm{sp}}}^{2}}.\end{eqnarray} \tag{ 4 }$$

For example, the scatter at intermediate masses ($\nu =1$) and accretion rates (${\rm{\Gamma }}=1$) is about 0.045 dex in both ${R}_{\mathrm{sp}}^{\mathrm{mn}}$ and ${M}_{\mathrm{sp}}^{\mathrm{mn}}$ and increases to about 0.055 dex for the 87th percentile. The lowest scatter of about 0.02 dex occurs at ${\rm{\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 (${\rm{\Gamma }}=10$), about 0.08 dex for ${R}_{\mathrm{sp}}^{\mathrm{mn}}$ and 0.1 dex for ${R}_{\mathrm{sp}}^{87 \% }$, resulting in a scatter of about 0.2 dex in ${{\rm{\Delta }}}_{\mathrm{sp}}^{87 \% }$.

Zoom In Zoom Out Reset image size

Table 4.  Best-fit Parameters

Parameter	${R}_{\mathrm{sp}}^{\mathrm{mn}}$	${R}_{\mathrm{sp}}^{ \% }$	${M}_{\mathrm{sp}}^{\mathrm{mn}}$	${M}_{\mathrm{sp}}^{ \% }$
Parameters for Rsp and Msp
a0	0.6498	0.3203	0.6792	0.2648
b0	0.6004	0.2674	0.4051	0.6660
${b}_{{\rm{\Omega }}}$	0.0920	0.1134	0.2919	0.1688
${b}_{\nu }$	0.0616	0.2080	0	0
c0	−0.8063	−0.9596	3.3659	4.7287
${c}_{{\rm{\Omega }}}$	17.5205	16.2459	1.4698	2.3889
${c}_{\nu }$	−0.2935	0	−0.0756	−0.0841
${c}_{{\rm{\Omega }}2}$	−9.6243	−9.4979	0	0
${c}_{\nu 2}$	0.0392	−0.0185	0	0
Meta-parameters for Dependence on Percentile
${a}_{{\rm{p}}}$	0	0.6148	0	0.8435
${b}_{{\rm{p}}}$	0	0.5452	0	−0.6392
${b}_{{\rm{\Omega }}{\rm{p}}}$	0	0	0	0.0032
${b}_{{\rm{\Omega }}{\rm{p}}2}$	0	0	0	4.9393
${b}_{\nu {\rm{p}}}$	0	−0.2233	0	0.2254
${c}_{{\rm{\Omega }}{\rm{p}}}$	0	0.0039	0	−0.7057
${c}_{{\rm{\Omega }}{\rm{p}}2}$	0	8.9691	0	−1.2419
${c}_{{\rm{\Omega }}2{\rm{p}}}$	0	−0.0005	0	0
${c}_{{\rm{\Omega }}2{\rm{p}}2}$	0	10.6132	0	0
${c}_{\nu {\rm{p}}}$	0	−0.4511	0	−0.3911
${c}_{\nu 2{\rm{p}}}$	0	0.0880	0	0.0742
Parameters for 68% Scatter (in dex)
${\sigma }_{0}$	0.0526	0.0445	0.0528	0.0276
${\sigma }_{{\rm{\Gamma }}}$	0.0038	0.0044	0.0025	0.0023
${\sigma }_{\nu }$	−0.0121	−0.0146	−0.0112	−0.0125
${\sigma }_{{\rm{p}}}$	0	0.0226	0	0.0473

Download table as:  ASCIITypeset image

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 Γ because the extrapolation in time is less reliable for rapidly evolving halos. In particular, the scatter is barely increased at low accretion rates (${\rm{\Gamma }}\,\lesssim 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σ (i.e., 95%) scatter can be slightly larger than twice the 1σ (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.

Before we discuss the various dependencies of R_sp and M_sp in detail, we summarize our results with a convenient fitting function. We find that the Γ–R_sp/R_200m and Γ–M_sp/M_200m relations are—at any redshift, cosmology, and peak height and for any R_sp definition—well fit by an expression similar to those suggested by Diemer & Kravtsov (2014) and More et al. (2015),

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

where ${X}_{\mathrm{sp}}$ can stand for either R_sp or M_sp and A, B, and C are free parameters. Those parameters are, in turn, functions of mass and redshift such that

$$\begin{eqnarray}A & = & {A}_{0}\\ B & = & ({B}_{0}+{B}_{{\rm{\Omega }}}{{\rm{\Omega }}}_{{\rm{m}}})\times (1+{B}_{\nu }\nu )\\ C & = & ({C}_{0}+{C}_{{\rm{\Omega }}}{{\rm{\Omega }}}_{{\rm{m}}}+{C}_{{\rm{\Omega }}2}{{\rm{\Omega }}}_{{\rm{m}}}^{2})\\ & & \times (1+{C}_{\nu }\nu +{C}_{\nu 2}{\nu }^{2}),\end{eqnarray} \tag{ 6 }$$

where we have introduced a total of nine free parameters. However, not all of these parameters are necessary to fit either R_sp or M_sp. In principle, there is no reason to expect that R_sp and M_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 (6) depend on the definition of R_sp. We fit ${R}_{\mathrm{sp}}^{\mathrm{mn}}$ and ${M}_{\mathrm{sp}}^{\mathrm{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

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

We constrain all free parameters simultaneously using a Levenberg–Marquart least-squares fit to the median Γ–R_sp relation at redshifts 0, 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 Γ values are weighted much more heavily than high Γ 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 Γ–R_sp relation. The fitting function should not be extrapolated beyond the 50th and 87th percentiles or the range of mass accretion rates for which it was constrained ($0\lt {\rm{\Gamma }}\lt 12$).

Figure 3 shows a summary of our main results and fitting function. Each column shows R_sp, M_sp, and Δ_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 Γ–Ω_m–ν parameter space where data are available. Figure 3 shows ${R}_{\mathrm{sp}}^{75 \% }$, but the same holds for ${R}_{\mathrm{sp}}^{\mathrm{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_sp and M_sp do not necessarily have to predict the correct median-enclosed-overdensity Δ_sp if the distributions are asymmetric. However, combining the fitting functions for R_sp and M_sp, we find that the result has roughly the expected agreement (bottom row of Figure 3; any error on R_sp is tripled; see Equation (4)). We have implemented our fitting function in the publicly available Python code Colossus.⁵

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_sp on halo mass, largely due to the limited accuracy of the R_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, ν). Interestingly, this mass dependence is weak for ${R}_{\mathrm{sp}}^{\mathrm{mn}}$ and ${R}_{\mathrm{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_sp/R_200m and M_sp/M_200m approach constant values at high Γ 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-Γ value with mass or redshift. The fits for R_sp and M_sp imply that Δ_sp also asymptotes to a particular value at high Γ, about 500 for ${R}_{\mathrm{sp}}^{\mathrm{mn}}$. At z = 0 and for our fiducial cosmology, ${{\rm{\Delta }}}_{{\rm{m}}}=500$ corresponds to ${{\rm{\Delta }}}_{{\rm{c}}}=135$, meaning that the average R_sp (using any definition) is always larger than ${R}_{200{\rm{c}}}$, even at very high accretion rates. At higher redshift, however, ${{\rm{\Delta }}}_{{\rm{c}}}$ becomes similar to ${{\rm{\Delta }}}_{{\rm{m}}}$, meaning that R_sp can reach radii smaller than ${R}_{200{\rm{c}}}$.

Another noteworthy feature of the fitting function presented in Equations (5)–(7) is that it encapsulates any dependence on redshift and cosmology into a dependence on ${{\rm{\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$, Ω_m barely changes, which manifests itself as a constant Γ–R_sp relation, even at z = 8. We find similar dependencies for other peak-height bins, as well as for higher percentiles.

Zoom In Zoom Out Reset image size

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 Γ–R_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_sp/R_200m at low redshifts due to its higher ${{\rm{\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 Ω_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.

Zoom In Zoom Out Reset image size

We further test the cosmology dependence using simulations that have ${{\rm{\Omega }}}_{{\rm{m}}}=1$ throughout, namely, self-similar cosmologies with a power-law power spectrum (Table 3). In such universes, we expect the Γ–R_sp relation not to evolve with redshift at all. We confirm this self-similar scaling, which constitutes further evidence that R_sp depends on Ω_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 Γ–R_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 ΛCDM cosmology (e.g., Figure 4 in Diemer & Kravtsov 2015). Our fitting function (with ${{\rm{\Omega }}}_{{\rm{m}}}=1$) matches the Γ–R_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 Ω_m is the variable that controls R_sp, not redshift.

By contrast, the self-similar simulation with a much shallower power spectrum slope, $n=-1$, leads to a rather different Γ–R_sp relation that is not described by our fitting function. We conclude that Ω_m is not the only parameter than influences R_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_sp. However, the impact of n in ΛCDM is degenerate with the effects of Ω_m and mass, making it difficult to disentangle the dependencies. Furthermore, it is not clear a priori how to define n in a Λ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.

Besides ν and Ω_m, our fitting function depends on the definition of R_sp and M_sp, i.e., whether we use the mean of the r_sp and m_sp distributions, their median, or higher percentiles. In Section 5, we demonstrate that different definitions can be useful for different purposes. Figure 6 shows the Γ–R_sp relation for different definitions of R_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 (7). 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.

Zoom In Zoom Out Reset image size

So far, we have considered R_sp primarily as a function of Γ and secondarily of other variables because Γ has the strongest effect. Unfortunately, Γ 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_sp in the absence of any knowledge about Γ, i.e., the median R_sp as a function of halo mass and redshift but marginalized over all mass accretion rates.

First, it is instructive to consider the average Γ 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 Γ 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 ν). Neglecting a mild dependence on redshift, we find that the logarithmic scatter in the Γ distributions is roughly

$$\begin{eqnarray}&&{\sigma }_{{\rm{\Gamma }}}\,\approx 0.41-0.07\nu ,\end{eqnarray} \tag{ 8 }$$

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

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

where

$$\begin{eqnarray}A & = & 1.2222+0.3515z\\ B & = & -0.2864+0.0778z-0.0562{z}^{2}+0.0041{z}^{3}.\end{eqnarray} \tag{ 10 }$$

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

Overall, the ∼0.07 dex scatter in the ν–R_sp relation is surprisingly low, considering the scatter in the Γ–R_sp relation and that the distribution of accretion rates is relatively extended. The low scatter allows for meaningful inferences regarding R_sp and M_sp in the absence of any knowledge about Γ.

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_sp for individual halos. Thus, More et al. (2015) used the stacked density profiles of Diemer & Kravtsov (2014) and defined R_sp as the radius where the logarithmic slope of the median profile is steepest. Based on this definition, they found that R_sp decreases as a function of Γ and increases slightly with increasing Ω_m—trends 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 Γ, 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 Γ–R_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 ${{\rm{\Gamma }}}_{\mathrm{DK}14}$ were computed based on ${M}_{\mathrm{vir}}$ in More et al. (2015) and ${M}_{200{\rm{m}}}$ in this work, but the difference is negligible given the accuracy of this comparison. We choose ${R}_{\mathrm{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 Γ. 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.

Zoom In Zoom Out Reset image size

More et al. (2015) also provided formulae for the ν–R_sp and ν–M_sp relations, again based on stacked density profiles. The dependence on ν was presumed to be due purely to the mass dependence of Γ. 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_sp at high peak height. Interestingly, Figure 4 of More et al. (2015) shows a hint of the reversed redshift evolution (lowered ν–R_sp relation at the highest redshifts), but the trend was not significant enough to be captured in their fitting function.

Zoom In Zoom Out Reset image size

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}_{\mathrm{sp}}^{75 \% }$ gives a good approximation, especially at high peak heights, where the results of More et al. (2015) were most constrained.

In Paper I, we undertook a halo-by-halo comparison of the results of Sparta and Shellfish. We found that ${R}_{\mathrm{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 Γ–R_sp relation measured by Sparta and Shellfish in Figure 11. The relations agree reasonably well, but, driven by systematic differences at low Γ, Shellfish prefers values of R_sp/R_200m that do not rise as sharply at low Γ. We further discuss the physical connection between splashback shells and apocenter distributions in Section 5.1.

Zoom In Zoom Out Reset image size

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_sp and M_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_sp is uniquely defined and that all particles reach their first apocenter exactly at R_sp. The values of r_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}_{\mathrm{sp}}^{\mathrm{mn}}$ as the definition for this comparison.

Zoom In Zoom Out Reset image size

We caution that ${{\rm{\Gamma }}}_{\mathrm{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, Γ is equal to s regardless of what interval Γ is averaged over, but realistic halo mass accretion rates tend to decrease with time at low redshifts, meaning that an average such as Γ 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}_{\mathrm{vir}}=1/2\,{r}_{\mathrm{turn} \mbox{-} \mathrm{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.⁶

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 ${{\rm{\Omega }}}_{{\rm{m}}}\,\ne 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 dashed lines in Figure 12 show fits to the numerical model predictions that were given for R_sp/R_200m and Δ_spand reconstructed for M_sp/M_200m. We note that the downturn in M_sp/M_200m at low Γ is present in the model but exaggerated due to the reconstruction from the fits to R_sp/R_200m and Δ_sp (X. Shi, private communication). The fitting functions were constrained for mass accretion rates in the range $0.5\lt s\lt 5$ (Shi 2016).

Both models correctly predict the general trend of decreasing R_sp/R_200m with mass accretion rate. In detail, however, the models do not match our results: the slope of the Γ–R_sp relation is steeper at low redshift and the normalization is higher at high redshift. The evolution of M_sp/M_200m with Γ is also steeper at low redshift and entirely different at high redshift. We discuss the physical reasons for these disagreements in Section 5.4.

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}_{\mathrm{sp}}^{\mathrm{mn}}$ should thus be the best definition of R_sp. However, it is not a priori clear that ${R}_{\mathrm{sp}}^{\mathrm{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 radius.⁷ 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}_{\mathrm{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_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_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}_{\mathrm{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_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.

Perhaps the most striking feature of the data presented in this paper, and thus the fitting function given in Equations (5)–(7), is their relative complexity, with significant dependencies on halo mass and Ω_m (at fixed Γ). 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_sp and conventional spherical overdensity radii. We note that R_sp/R_Δ would vary more strongly if spherical overdensity radii other than ${R}_{200{\rm{m}}}$ were used.

At a fixed R_sp, the ratio R_sp/R_200m depends on the mass profile around R_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 Γ in our fitting model.

Another open question is related to the dependence of R_sp/R_200m on cosmological parameters. We have shown that the scaling with Ω_m works for the WMAP and Planck cosmologies and is thus likely appropriate for any realistic cosmology. However, R_sp/R_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 Ω_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 ${{\rm{\Gamma }}}_{\mathrm{dyn}}$, but perhaps R_sp exhibits tighter correlations with other definitions that we have not yet considered (e.g., definitions based on shorter or longer timescales or relying directly on M_sp instead of ${M}_{200{\rm{m}}}$). 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_sp and conventional definitions is the meaning of the overdensity Δ. In the context of spherical overdensity radii, Δ is the fundamental quantity that determines R_Δ and ${M}_{{\rm{\Delta }}}$. In the splashback picture, Δ_sp is merely a consequence of independently determined R_sp and M_sp. Our results show that not only does Δ_sp vary systematically depending on a number of variables but also the scatter in Δ_sp is larger than the scatter in R_sp and M_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 Δ, R_Δ can change suddenly when a massive subhalo is accreted. In contrast, R_sp changes more slowly in this case, while M_sp and thus Δ_sp change rapidly.

In summary, we have little reason to expect a simple relation between R_sp and R_Δ. While we employ the quantities R_sp/R_200m and M_sp/M_200m to establish a connection between R_sp and R_Δ, R_sp is an independent definition of the halo boundary that cannot easily be expressed in terms of spherical overdensity radii.

The measurement of R_sp performed by More et al. (2016) and confirmed by Baxter et al. (2017) indicates a surprisingly small splashback radius, namely, ${R}_{\mathrm{sp}}/{R}_{200{\rm{m}}}=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}_{\mathrm{sp}}/{R}_{200{\rm{m}}}=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}_{\mathrm{sp}}^{75 \% }$ at the given peak height and redshift. Given the scatter in the ν–R_sp relation, the observed R_sp would represent a 2σ 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.

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 Λ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_sp/M_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\gt 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.