The Diffusion Coefficient of the Splashback Mass Function as a Probe of Cosmology

The concept of a splashback radius of a dark matter (DM) halo, defined as the distance from the halo center to the apocenter of the first orbit of the newest infalling objects (Adhikari et al., 2014; More et al., 2015) has recently risen up to limelight. Unlike the virial radius, it properly includes the flyby subhalos (More et al., 2015), does not suffer from pseudo-evolution (Diemer et al., 2013), and coincides well with the location on the halo periphery where the stacked density profile exhibits an abrupt decline (Wang et al., 2009; Wetzel et al., 2014; Diemer & Kravtsov, 2014; Adhikari et al., 2014; Diemer, 2020c). Given these favorable properties, the splashback radius is often dubbed ”a physically motivated boundary”, which naturally sunders objects orbiting in the inner region of a halo from still infalling ones in its outer region (More et al., 2015).

Very recently, Diemer (2020a) found another auspicious aspect of the splashback boundary from a N-body experiment. Determining the splashback mass of each DM halo as the total mass of the DM particles enclosed by its splashback radius, he constructed an analytic formula for the splashback mass function by employing the excursion set ansatz and claimed that it exhibits a substantially higher degree of universality than the virial mass function. In the excursion theory, the mass function of DM halos is referred to as being universal, if its multiplicity (see Section 2.1 for a definition) does not have explicit dependence on the redshift and shape of the linear density power spectrum (Jenkins et al., 2001). The success of the nearly universal form of the splashback mass function provided by Diemer (2020a) was demonstrated by its good agreements with the numerical results at all redshifts for the case of the self-similar cosmologies as well as for the standard cosmology where the cosmological constant ($\Lambda$) and cold dark matter (CDM) are the most dominant energy and matter contents of the universe, respectively, whether the initial conditions are described by the key cosmological parameters from the Planck experiment (Planck Collaboration et al., 2014) or from the Seven-Year Wilkinson Microwave Anisotropy Probes (hereafter, WMAP7 Komatsu et al., 2011).

The excursion set theory describes the growth of the initially overdense regions into DM halos as a random walking process bounded by a prescribed absorbing barrier, and evaluates the virial mass function of DM halos by enumerating the random walks that hit for the first time the absorbing barrier on a given virial mass scale (Bond et al., 1991). The degree of the universality of the virial mass function is contingent upon the shape, height and statistical nature of the absorbing barrier, which in turn is determined by the underlying collapse dynamics, background cosmology, and halo-identification procedure (Chiueh & Lee, 2001; Sheth et al., 2001; Maggiore & Riotto, 2010a, b). For instance, the absorbing barrier would have a flat shape and redshift-independent constant height, yielding a completely universal mass function, if the DM halos are formed via spherically symmetric collapse process and identified without any ambiguity in an Einstein de Sitter universe where the matter density parameter equals unity ($\Omega_{\rm m}=1$) (Press & Schechter, 1974; Bond et al., 1991).

The departure of the true gravitational collapse from spherical symmetry and the presence of dark energy lead the absorbing barrier to have a curved shape with weakly $z$-dependent height (Lahav et al., 1991; Eke et al., 1996; Kitayama & Suto, 1996; Sheth et al., 2001), which in consequence slightly decreases the degree of the universality of the excursion set mass function of DM halos. Moreover, in realistic situations, the virial radius of a DM halo is susceptible to the pseudo-evolution (Diemer et al., 2013) as well as to the tidal disturbance from the surrounding cosmic web, which implies that considerable amount of ambiguity should be involved in the determination of the virial mass (Robertson et al., 2009; Maggiore & Riotto, 2010a, b). To take into account this ambiguity, the absorbing barrier degrades into a stochastic variable whose standard deviation varies with redshift and key cosmological parameters (Maggiore & Riotto, 2010b; Corasaniti & Achitouv, 2011a; Corasaniti, & Achitouv, 2011b). Thus, the degree of the universality of the excursion set virial mass function is most closely linked with that of the stochasticity of the absorbing barrier.

Meanwhile, if the boundary of a halo is set at its splashback radius which does not experience pseudo-evolution and incorporates properly the orbits of the flyby objects ejected by the tidal disturbance (More et al., 2015), the identification of a DM halo and determination of its mass are expected to suffer much less from ambiguity, rendering the absorbing barrier to become less stochastic. The detected higher degree of universality of the splashback mass function in the seminal work of Diemer (2020a) can be naturally explained by this reduced degree of stochasticity of the absorbing barrier. What Diemer (2020a) focused on, however, was not providing a physical explanation for the nearly universal form of the spashback mass function in terms of the less stochastic absorbing barrier, but finding a practical formula for it characterized by redshift independent multiple parameters. In other words, even though it used the excursion set ansatz and turned out to be universal and robust against the variations of the background cosmologies, the formula of Diemer (2020a) characterized by multiple free parameters is not a physical model but a mere fitting formula.

Here, we will attempt to find a physical model for the splashback mass function in the generalized excursion set formalism. Instead of preoccupying with expressing the splashback mass function in an universal form, we will focus on deriving it from the first principles and reducing the number of free parameters that characterize it. The plan of this paper is as follows. In Section 2.1, we review the excursion set theory for the virial mass function. In Section 2.2 we present how the slashback mass function can be analytically obtained in the generalized excursion set theory. In Section 3, we present the results of numerical testing and model selection. In Section 4 we summarize the implication of our findings and discuss its prospect as a cosmological probe as well.

Diemer (2017) developed an elaborate code, called SPARTA, to efficiently locate the splashback boundaries of DM halos. This code basically follows the orbits of all DM particles in and around a given halo identified by the Rockstar halo-finder algorithm (Behroozi et al., 2013) and determines the distances to the apocenters of their first orbits, $r_{\rm sp}$. To derive the splashback radius of each halo, $R_{\rm sp}$, from the distribution of $r_{\rm sp}$, Diemer (2017) considered several different definitions. For instance, $R_{\rm sp}^{\rm min}$ is the splashback radius of a halo defined as the average of the $r_{\rm sp}$ values, while $R_{\rm sp}^{Q\%}$ is defined as the Q-th percentile of the distribution of $r_{\rm sp}$. For our analysis, we adopt mainly the former definition, $R_{\rm sp}^{\rm min}$, and refer to the mass enclosed by $R_{\rm sp}^{\rm min}$ as the splashback mass of a halo, unless otherwise stated.

Diemer et al. (2017) applied this code to the Erebos suite of 14 DM only N-body simulations for three different (Planck, WMAP7 and self-similar) cosmologies and produced the halo catalogs ²²2They are available at the website of http://erebos.astro.umd.edu/erebos/ . that provided additional information on the splashback radius and mass of each halo (see also Diemer, 2020b). From the Erebos suite, we utilize two simulations of $1024^{3}$ DM particles run by Diemer & Kravtsov (2014) and Diemer & Kravtsov (2015) with the GADGET2 code (Springel, 2005) on a periodic box of linear size $500\,h^{-1}$Mpc for the Planck and WMAP7 cosmologies, respectively. Table 1 lists the mass of individual DM particle ($m_{p}$) and value of $\Omega_{m}$ adopted for each simulation in the first three columns.

From the halo catalog of each simulation at $z=0$, we select only those DM halos whose splashback masses, $M$, are equal to or larger than $5\times 10^{12}\,h^{-1}\,M_{\odot}$ (i.e., galaxy group and cluster scales), corresponding to the particle number cut of $500$ and $575$ for the Planck and WMAP7 cosmologies, respectively. Splitting the range of $\ln M$ into short intervals of length $d\ln M$, we determine the number densities of the selected DM halos, $dN/d\ln M$, at each bin. We first fit the modified Corasaniti-Achitouv formula to this numerically obtained $dN/d\ln M$ by adjusting the values of $D_{B}$ and $\beta$ via the $\chi^{2}$-minimization:

where $n^{\rm num}(m_{i})$ and $n^{\rm ana}(m_{i};D_{B},\beta)$ denote the numerical and analytical values of $dN/d\ln M$, respectively, at the $i$th logarithmic mass bin, $m_{i}\equiv\ln M_{i}$. Here, the error $\sigma_{n}$ at each bin is obtained through the following Jackknife analysis. Splitting the halo sample into eight subsamples according to the halo positions, we evaluate $dN/d\ln M$ separately for each subsample. Then, the Jackknife error of the splashback mass function at each logarithmic mass bin, $\sigma(m_{i})$, is obtained as one standard deviation scatter of $dN/d\ln M$ among the eight subsamples.

Figure 1 plots the results (black filled circles) with the Jackknife errors as well as the modified Corasaniti-Achitouv  mass function (blue solid lines), Equation (11), with the best-fit values of $D_{B}$ and $\beta$ at $z=0$ for the Planck and WMAP7 cosmologies (in the first and third panels from the top, respectively). The ratio of the analytical to the numerical splashback mass functions, $n^{\rm ana}/n^{\rm num}$, are also shown for the two cosmologies (in the second and fourth panels from the top, respectively). As can be seen, the modified Corasaniti-Achitouv mass function indeed agrees very well with the numerical results, exhibiting, $1\leq n^{\rm ana}/n^{\rm num}\leq 1.1$ in the mass range of $0.5\leq M/(10^{13}\,h^{-1}\,M_{\odot})\leq 50$ for both of the cosmologies. In the higher mass range of $M/(10^{13}\,h^{-1}\,M_{\odot})>50$ where the numerical results suffer from larger errors, we note a substantial discrepancy, $n^{\rm ana}/n^{\rm num}\geq 1.5$.

We investigate if this analytical splashback mass function based on the generalized excursion set theory is valid even for the case that different definitions of $R_{\rm sp}$ are used to determine the splashback masses. Figures 2 and 3 plot the same as Figure 1 but for four different cases that the splashback radii are defined by $R_{\rm sp}^{50\%}$ (top-left panel), $R_{\rm sp}^{70\%}$ (top-right panel), $R_{\rm sp}^{80\%}$ (bottom-left panel) and $R_{\rm sp}^{90\%}$ (bottom-right panel) for the Planck and WMAP7 cosmologies, respectively. As can be seen, the validity of the Corasaniti-Achitouv  formula in the mass range of $0.5\leq M/(10^{13}\,h^{-1}\,M_{\odot})\leq 50$ is fairly robust against the variation of the definition of $R_{\rm sp}$ for both of the cosmologies. It is, however, worth mentioning that for these plots, we newly determine the best-fit values of $D_{B}$ and $\beta$ for each case via Equation (14). That is, the two parameters, $D_{B}$ and $\beta$, are dependent on the percentiles, $Q$. For the case of $R_{\rm sp}^{50\%}$ the best-fit values are found to be almost identical to those for the case of $R^{\rm min}_{\rm sp}$, while the other three cases yield substantially different best-fit parameters. We do not pursue here an analytical modeling of the percentile dependence of the two parameters and consistently use $R_{\rm sp}^{\rm min}$ as a definition of the splashback radius.

We investigate if and how the splashback mass functions differ in the best-fit values of $D_{B}$ and $\beta$ from the virial mass functions, for which we numerically determine $dN/d\ln M$ as a function of their virial masses and analytically evaluate the original Corasaniti-Achitouv formula with $\delta_{sc}=1.686$. Figure 4 plots the $68\%,\ 95\%$ and $99\%$ contours of $\chi^{2}$ (red, blue and green lines, respectively) in the $\beta$-$D_{B}$ configuration space at six different redshifts in the range of $0\leq z\leq 2$ for the two cosmologies. As can be seen, the values of $\beta$ and $D_{B}$ for the splashback mass functions (solid lines) are progressively lower than those for the virial mass functions (dotted lines), which implies the following: For the case that the splashback rather than the virial radius is used as the halo boundary, the formation process of a DM halo follows much better the spherical collapse dynamics and thus its identification becomes much less ambiguous.

Noting that the best-fit values of the drifting average coefficient, $\beta$, shown in Figure 4 are quite close to zero at all of the redshifts considered, we set it at zero and suggest that the modified Maggiore-Riotto  formula, Equation (12), with only one parameter $D_{B}$ should be a more efficient model for the splashback mass function. Figure 1 also plots the modified Maggiore-Riotto  formula with the best-fit value of $D_{B}$ (red solid lines) and its ratio to the numerical results. The comparisons of the modified Maggiore-Riotto and Corasaniti-Achitouv formulae with the numerical results at three higher redshifts in the range of $0\leq z\leq 2$ can be witnessed in Figure 5. The results shown in Figures 1 and 5 demonstrate that even though it has only one free parameter, the modified Maggiore-Riotto  formula can describe the numerically obtained splashback mass function as validly as the modified Corasaniti-Achitouv formula in the redshift range of $0\leq z\leq 2$, which justifies our adoption of $\beta=0$. We also test and confirm the robustness of the Maggiore-Riotto  formula against the definitions of $R_{\rm sp}$, the results of which are shown (red solid lines) in Figures 2-3.

It is worth recalling the numerical finding of García et al. (2021) that when the halo boundaries were defined as the radii that optimize the halo fits, the Press-Schechter like formula with a modified spherical barrier height of $\delta_{sc}=1.45$ lower than the top-hat solution of $\delta_{sc}=1.68$ matched quite well the numerically obtained halo mass function. Given that the halo boundaries defined by García et al. (2021) are slightly larger than the splashback radii, the success of the modified Maggiore-Riotto formula with $\delta_{sp}=1.52$ shown in Figures 1-5 is in line with their findings, providing a physical explanation to them.

Figure 6 plots the resulting best-fit values of $D_{B}$ (filled circles) as a function of $z$ for both of the cosmologies. Here, the error in the best-fit value of $D_{B}$ is found as the inverse of the square root of the Fisher information. As can be seen, the best-fit value of $D_{B}$ appears to decrease almost linearly with the increment of $z$ toward zero, which motivates us to model $D_{B}(z)$ by the following linear function of $z$,

characterized by two parameters, the slope, $A_{D}$, and threshold redshift, $z_{c}$, satisfying $A_{D}\equiv|dD_{B}(z)/dz|$ and $D_{B}(z_{c})=0$, respectively. Employing the $\chi^{2}$-minimization method, we statistically compare Equation (15) with the numerically obtained $D_{B}(z)$ to find the best-fit values of $A_{D}$ and $z_{c}$ as well as their errors for the two cosmologies, which are listed in the fourth and fifth columns of Table 1. Note that a significant difference in $z_{c}$ but no difference in $A_{D}$ are found between the two $\Lambda$CDM cosmologies, which hints that $z_{c}$ may sensitively depend on $\Omega_{m}$.

The result shown in Figure 6 also indicates that the absorbing barrier for the splashback mass function becomes less and less stochastic as $z$ increases, until it eventually turns to being deterministic at $z=z_{c}$. Figure 7 compares the numerical results (black filled circles) with the parameter free formula given in Equation (13) at $z=z_{c}$ for the two cosmologies. As can be seen, good agreements are found between them ($1\leq n^{\rm ana}/n^{\rm num}<1.1$) in the mass range of $0.5\leq M/(10^{13}\,h^{-1}\,M_{\odot})\leq 2$ for both of the cosmology, even though no empirical adjustment is made at all. Yet, in the higher mass section, the large errors make it difficult to make a proper comparison.

We have so far considered a rather restricted mass-range of $0.5\leq M/(10^{13}\,h^{-1}\,M_{\odot})\leq 10^{2}$, mainly because of the limited mass-resolution of the simulation box of $L_{\rm box}=500\,h^{-1}$Mpc and partially because the splashback mass functions are most useful as a cosmological diagnostics on the galaxy group and cluster scales. To examine if the validity of the modified Maggiore-Riotto formula may also be extended to a wider mass range, we repeat the same analysis but with the data obtained from the simulation boxes of two different linear sizes, $L_{\rm box}=125$ and $250\,h^{-1}$Mpc for both of the cosmologies. Two additional datasets from the larger simulation boxes with $L_{\rm box}=1000$ and $2000\,h^{-1}$Mpc are also available for the WMAP7 cosmologies. For a detailed information on these simulations with various $L_{\rm box}$, see Diemer (2017).

Figure 8 and 9 compare the modified Maggiore-Riotto  formula with the numerical results at six different redshifts from the three and five different simulation boxes for the Planck and WMAP7 cosmologies, respectively. For the analytic formula, we do not newly determine the best-fit value of $D_{B}$ for each case of $L_{\rm box}$ but use the same best-fit value obtained for the case of $L_{\rm box}=500\,h^{-1}$Mpc. As can be seen, the modified Maggiore-Riotto  formula with the same best-fit value of $D_{B}$ still agrees quite well with the numerical results in the low-mass section of $0.1\leq M/(10^{12}\,h^{-1}\,M_{\odot})\leq 1$ from the smaller simulation boxes for both of the cosmologies, exhibiting $0.9<n^{\rm ana}/n^{\rm num}<1.2$, even though no new adjustment of $D_{B}$ is made. That is, the validity of the modified Maggiore-Riotto formula for the splashback mass function is not limited to the high-mass section, unlike the case of the virial mass function for which the Maggiore-Riotto formula breaks down even in the intermediate mass section of $10^{13}\leq(M/h^{-1}\,M_{\odot})<10^{14}$. Meanwhile, for the WMAP7 cosmology, in the higher mass section of $1\leq M/(10^{15}\,h^{-1}\,M_{\odot})\leq 5$ where the numerical results carry relatively larger errors, the agreements between the modified Maggiore-Riotto  formula with the fixed value of $D_{B}$ and the numerical results from the larger simulation boxes turn out to be not so good as in the lower-mass section.

Using the data from this smaller simulation box, we also examine whether or not at $z\geq z_{c}$ the parameter free formula, Equation (13), is capable of matching the numerical results in a wider mass range. Figure 10 plots the parameter free formula and the numerical results from the simulation boxes of $L_{\rm box}=125$ and $500\,h^{-1}$Mpc at three redshifts higher than $z_{c}$ in the range of $2.5\leq z\leq 3.2$. As can be seen, in the mass range of $0.1\leq M/(10^{12}\,h^{-1}\,M_{\odot})<10$ the parameter free analytic formula agrees quite well with the numerical results, having $1<n^{\rm ana}/n^{\rm num}<1.15$ at all of the three redshifts. This result is consistent with the aforementioned claim that the collapse barrier for the halos identified by the splashback radii becomes deterministic at $z\geq z_{c}$.

Now that the modified Maggiore-Riotto  formula for the splashback mass function has been found to work very well in the wide ranges of $M$ and $z$, we would like to statistically decide which model for the splashback mass function is most preferred by the numerical data among the modified Corasaniti-Achitouv, Maggiore-Riotto, Sheth-Tormen, and Tinker et al.  formulae. Here, the modified Sheth-Tormen  mass function is obtained by replacing $\delta_{sc}$ in Equation (6) by $\delta_{sp}$ as done for the modified Corasaniti-Achitouv  and Maggiore-Riotto  formulae and newly determining the best-fit values of its three parameters for the splashback mass function. As for the modified Tinker et al.  mass function in which no term contains $\delta_{sc}$, we obtain it by newly determining the best-fit values of the four parameters by fitting Equations (1), (4) and (10) to the numerically obtained splashback mass function.

To select the model most preferred by the numerical results, we calculate the Bayesian information criterion (BIC) of each model, defined as (Schwarz, 1978)

where $N_{\rm pa}$ and $N_{\rm pt}$ denote the numbers of the free parameters of each model and the number of the experimental data points, respectively, and $\chi^{2}_{\rm min}$ denotes the minimum value of $\chi^{2}$ for each model. The model with lower value of ${\rm BIC}$ is more preferred by the numerical results.

Figure 11 plots the BIC values of the modified Corasaniti-Achitouv, Maggiore-Riotto, Sheth-Tormen, Tinker et al.  (blue, red, green, violet sold lines, respectively) formulae as a function of $z$ for the two $\Lambda$CDM cosmologies. As can be seen, the modified Maggiore-Riotto  formula characterized by only one single parameter, $D_{B}$, is found to have the smallest BIC value in the redshift range of $0.3<z\leq 3$, which proves that the modified Maggiore-Riotto formula is indeed the most preferred analytic model for the splashback mass function among the four in this redshift range.

Yet, it should be noted that the modified Maggiore-Riotto and Corasaniti-Achitouv formulae exhibit sharp increases with the decrement of $z$ in the lower redshift range of $z\leq 0.3$, where the modified Sheth-Tormen formula yields the smallest BIC values. To understand the reason for this abrupt increases of the BIC values at $z\leq 0.3$, we examine the $\chi^{2}$ values and find that both of the Corasaniti-Achitouv and Maggiore-Riotto formulae yield relatively high values of $\chi^{2}$ in the low mass section $M\leq 10^{13}\,h^{-1}\,M_{\odot}$ at $z\leq 0.3$, which in turn leads to produce larger $BIC$ values at these low redshifts. Given that the more nonlinearly evolved tidal fields at lower redshifts (Bond et al., 1996; Lee & Springel, 2010) are expected to more severely disturb the lower-mass DM halos (Bernardeau, 1994), we suspect that the single diffusion coefficient $D_{B}$ introduced by the generalized excursion set theory (Maggiore & Riotto, 2010a, b) is no longer capable of accounting for the stochasticity involved in the halo identifications due to the severe tidal disturbance at $z\leq 0.3$, even when the splashback radii are used as the halo boundaries.

We also want to examine which one is preferred by the numerical results between the modified Maggiore-Riotto  with one parameter and the universal formula with four parameters proposed by Diemer (2020a). Although the same functional form of the Tinker et al. multiplicity function, i.e., Equation (10), was adopted by Diemer (2020a) to construct his universal formula, the best-fit parameters found by fitting the modified Tinker et al. formula to the numerically obtained splashback mass function $dN/d\ln M$ at a fixed redshift turn out to be different from those found by fitting the universal formula to the numerically obtained $f(\nu)$ with $\nu\equiv\delta_{sc}/\sigma(M,z)$ over a broad range of $M$ and $z$. For this examination, the BIC method is not applicable since it is impossible to specify the number of data points, $N_{\rm pt}$, at a fixed redshift for the case of the universal formula of Diemer (2020a). Instead, we calculate the alternative Akaike Information Criterion (AIC) defined as (Cavanaugh & Neath, 2019)

and known to complement the BIC method in case that the fitting procedure is performed over a different range of an independent variable between the two target models.

Figure 12 plots the AIC values of the modified Maggiore-Riotto (red solid lines) and the universal formula of Diemer (2020a) (blue solid lines) in the redshift range of $0\leq z\leq 3$ for the two $\Lambda$CDM cosmologies. For the latter, we adopt the same best-fit values of the four parameters as given in Diemer (2020a). As can be seen, the modified Maggiore-Riotto formula yields a much lower value of AIC than the universal formula of Diemer (2020a) in this redshift range, indicating that the former is more preferred by the numerical results at least in this redshift range for the two $\Lambda$CDM cosmologies.

We have constructed an efficient physical model for the splashback mass function in the framework of the generalized excursion set theory. Noting that for the cases of the halos identified by their splashback radii the absorbing barrier of the excursion set formalism has a lower degree of stochasticity without drifting on average from the spherical height (Figure 4), we have modified the Maggiore-Riotto  formula (Maggiore & Riotto, 2010a, b) by replacing the top-hat solution of the spherical height, $\delta_{sc}=1.686$, by the self-similar spherical infall solution, $\delta_{sp}=1.52$ (Shapiro et al., 1999). The modified Maggiore-Riotto  formula, Equation (12), has been found to be in good accord with the numerically obtained splashback mass function from the Erbos simulations of linear size $500\,h^{-1}$Mpc at $0\leq z\leq 3$ on the galaxy group and cluster scales (e.g., $0.5\leq M/[10^{13}\,h^{-1}\,M_{\odot}]<10^{2}$ at $z=0$) for both of the Planck and WMAP7 cosmologies (Figure 5).

We have also modified the Sheth-Tormen, Corasaniti-Achitouvand Tinker et al.  formulae in the same manner to accommodate the splashback mass function and compared each of them with the numerical results to determine their best-fit parameters. Applying the BIC method to the four analytic formulae, we have found that among them the numerically obtained splashback mass functions prefer most the modified Sheth-Tormen  and Maggiore-Riotto  formulae at $z\leq 0.3$ and $0.3<z\leq 3$, respectively (Figure 11). With the help of the AIC method, it has been also shown that the numerical results prefer the modified Maggiore-Riotto  formula to the universal formula of Diemer (2020a) at $0\leq z<3$.

Despite that the modified Maggiore-Riotto formula is not so universal as that of Diemer (2020a) and that it turns out to be not preferred most by the BIC method in the redshift range of $0\leq z\leq 0.3$, we believe that it has a couple of advantages over the previous formulae. First, since it is a physical model derived analytically from the generalized excursion set theory, taking into proper account of the merit of identifying the DM halos with their splashback boundaries, it gives us a much better insight into the underlying connection between the initial conditions and the halo abundance. Second of all, its single parameter, the diffusion coefficient, $D_{B}$, has turned out to vary almost linearly with redshift, unlike the cases of the other formulae whose multiple parameters show almost desultory fluctuations with redshifts. Using this linear relation of $D_{B}(z)$, we have determined the critical redshift, $z_{c}$, at which $D_{B}$ vanishes (Figure 6). A crucial implication of this result is that the splashback mass functions at $z\geq z_{c}$ can be described by the purely analytical parameter-free formula (Figure 7). Besides, it has turned out that the Planck and the WMAP7 cosmologies yield different values of $z_{c}$, which implies that $z_{c}$ defined as $D_{B}(z_{c})=0$ might be useful as a complementary probe of $\Omega_{m}$ (Table 1).

It should be worth discussing here why the diffusion coefficient of the modified Maggiore-Riotto formula exhibits a linear decrement with redshift and why the value of the critical redshift $z_{c}$ depends sensitively on $\Omega_{m}$. One of the primary factors that make the halo identification ambiguous is the tidal disturbance from the surrounding cosmic web (Maggiore & Riotto, 2010a, b). At higher redshifts, $z\sim z_{c}$, when the tidal field is still at its quasi-linear stage, the DM halos can be identified unambiguously by the splashback boundary. At lower redshifts, $z<z_{c}$, however, even the splashback boundaries would be disturbed by the more nonlinearly evolved tidal field (Bond et al., 1996; Lee & Springel, 2010), which leads to a higher degree of stochasticity of the absorbing barrier and accordingly to a larger value of $D_{B}$.

The above logic also allows us to understand why $z_{c}$ has a larger value in the Planck cosmology than in the WMAP7 cosmology. In the Planck cosmology with a higher value of $\Omega_{m}$ a more rapidly evolving tidal fields would disturb the halo identification at earlier times, which in turn deviates $D_{B}$ from zero at higher redshifts than in the WMAP7 cosmology. Yet, we have no clue to understanding why the decrease of $D_{B}$ follows a linear scaling with $z$. A related important question is whether or not this linear scaling of $D_{B}(z)$ is universal throughout the background cosmologies including the non-standard ones. If it turns out to be the case, then the critical density, $z_{c}$, would become a powerful probe not only of $\Omega_{m}$ but also of the background cosmology. Our future work will be in the direction of physically modeling $D_{B}(z)$, extending our analysis to the non-standard cosmologies and testing the power of $z_{c}$ as a probe of cosmology. We hope to report the results elsewhere in the near future.