The Splashback Mass Function in the Presence of Massive Neutrinos

The dark matter halos form through gravitational collapse of the initially overdense sites. As the gravitational collapse is a highly nonlinear process occurring after the shell crossing, the linear or even higher-order perturbation theory is not eligible for the analytical description of the halo formation. Only for the extreme case that the collapse of an isolated overdense site occurs completely in a spherically symmetry way, the halo formation can be analytically described by the top-hat spherical dynamics, which basically predicts that the virialized halos form from the sites whose linearly extrapolated density contrast exceeds a certain value, called the virial density threshold (Gunn & Gott, 1972). A unique aspect of the top-hat spherical dynamics is that the virial density threshold is independent of halo mass and quite insensitive to the initial conditions of the universe (e.g., Press & Schechter, 1974; Lahav et al., 1991; Eke et al., 1996; Pace et al., 2010).

In realistic cases where the overdense sites are usually not isolated and the gravitational collapse does not proceed in a spherically symmetric way, the simple analytical description of the halo formation in terms of a constant value of the virial density threshold is no longer attainable (Bond & Myers, 1996; Sheth & Tormen, 1999; Chiueh & Lee, 2001; Sheth et al., 2001; Tinker et al., 2008). Nevertheless, several theoretical works based on $N$-body simulations indicated that the condition for the realistic formation of dark matter halos can still be expressed in terms of the virial density threshold by treating it as a mass-dependent stochastic variable rather than a deterministic constant value (e.g., Robertson et al., 2009; Maggiore & Riotto, 2010a, b; Corasaniti & Achitouv, 2011a; Corasaniti, & Achitouv, 2011b). The mass-dependence of the mean value of the virial density threshold stems from the departure of the gravitational collapse from the top-hat spherical dynamics, while its stochastic nature reflects the ambiguity in the identification of non-isolated halos.

In the high-mass limit corresponding to the massive cluster halos, the virial density threshold becomes less non-spherical but more stochastic (e.g., Robertson et al., 2009; Maggiore & Riotto, 2010a, b). Since the gravitational collapse of highest density peaks occur more or less in a spherically symmetric way (Bernardeau, 1994), the mean value of the virial density threshold converges to the spherical value, $\delta_{sc}=1.686$, (Gunn & Gott, 1972). On the other hand, the massive cluster halos often reside at the nodes of cosmic filaments and thus are more vulnerable to the disturbance from the surroundings. In consequence, it becomes more ambiguous to identify their clear-cut boundaries, which results in a higher degree of stochasticity of the virial density threshold in this limit (Maggiore & Riotto, 2010b).

Very recently, it has been suggested that the ambiguity in the cluster identification can be considerably reduced by using the splashback radius as a physical boundary instead of the virial counterpart (Ryu & Lee, 2021). The splashback radius of a DM halos was originally found in $N$-body simulations as its outer limit where the density profile exhibits an abrupt drop (Wang et al., 2009; Wetzel et al., 2014; Diemer & Kravtsov, 2014; Adhikari et al., 2014; Diemer, 2021). It corresponds to the radial distance to the apocenter of the elliptical orbit of the latest infalling particle (Adhikari et al., 2014; More et al., 2015), distinguishing between infalling and orbiting particles, and naturally taking into account the fly-by objects (More et al., 2015). If the halo boundary is defined by its splashback radius, which is in fact a solution to the self-similar spherical infall model (Fillmore & Goldreich, 1984; Bertschinger, 1985), the unrealistic assumption of isolated overdense region is no longer required to describe the halo formation, since the splashback boundary generically incorporates the pseudo-evolution of halo mass (Diemer et al., 2013).

Noting that the self-similar spherical infall model gives a lower mean value of the splashback density threshold (Shapiro et al., 1999), Ryu & Lee (2021) substituted it for the virial counterpart in the the generalized excursion set formalism (Maggiore & Riotto, 2010a, b; Corasaniti & Achitouv, 2011a; Corasaniti, & Achitouv, 2011b) and proposed an analytic formula for the splashback mass function characterized by a single parameter, diffusion coefficient, which measures the degree of the stochasticity of the splashback density threshold. Testing the analytically derived splashback mass functions against the $N$-body results for the Planck and WMAP7 cosmologies (Planck Collaboration et al., 2014; Komatsu et al., 2011), Ryu & Lee (2021) verified its validity in a wide mass range up to the redshift $z\sim 3$. Showing that the redshift evolution of the diffusion coefficient sensitively depends on the matter density parameter $\Omega_{m}$, Ryu & Lee (2021) suggested that the diffusion coefficient of the splashback mass function of dark matter halos be a sensitive probe of the initial conditions of the universe.

In this paper, we explore if the diffusion coefficient of the splashback mass function can effectively constrain the total neutrino mass, $\sum m_{\nu}$, speculating that the degree of the stochasticity of the splashback density threshold may be sensitive not only to the amount of dark matter but also to its nature. The upcoming Sections contain the following contents: a brief review of the analytic formula for the splashback mass function (Section 2); a comparison between the analytic formula and the numerically obtained splashback mass function in the presence of neutrinos and a difference in the redshift evolution of the diffusion coefficient between the massless and massive neutrino cases (Section 3); a summary and conclusion (Section 4). Throughout this paper, we will assume a cosmology where the neutrino ($\nu$) is present, and the cosmological constant ($\Lambda$) and cold dark matter (CDM) are dominant, i.e., $\nu\Lambda$CDM cosmology.

In the original excursion set formalism (Press & Schechter, 1974; Bond et al., 1991), the virial mass function, $dN(M,z)/d\ln M$, defined as the number densities of DM halos as a function of its logarithmic virial mass, $\ln M$, at redshift $z$, can be analytically written as

where $\bar{\rho}$ is the mean mass density of the universe, $\sigma(M,z)$ is the rms fluctuation of the linear density field on the mass scale of $M$ at redshift $z$, and $f(\sigma)$ is the multiplicity function that counts the number of the initial density peaks whose linearly extrapolated density contrast exceeds some virial density threshold when the rms density fluctuation has the value of $\sigma$. The excursion set theory relates $dN/d\ln M$ to the background cosmology through $\sigma(M,z)$ that is often computed as $\sigma^{2}(M,z)\equiv(2\pi^{2})^{-1}\int\,dk\,k^{2}P(k,z)W^{2}(k,M)$ where $P(k,z)$ is the linear density power spectrum and $W(k,M)$ is the top-hat window function on the virial mass scale $M$ (Press & Schechter, 1974; Bond et al., 1991).

This original excursion set theory was generalized by Maggiore & Riotto (2010a, 2010b) (hereafter, Maggiore-Riotto) who treated the virial density threshold as a stochastic random variable whose mean value is $\delta_{sc}=1.686$ according to the top-hat spherical model (Gunn & Gott, 1972). The resulting Maggiore-Riotto model turned out to describe very well the virial mass function of cluster halos, despite that it has only one fitting parameter, called the diffusion coefficient, which characterizes the degree of the stochasticity of the virial density threshold. Noting the limited validity of the Maggiore-Riotto formula to the cluster mass scale, Corasaniti & Achitouv (2011a, 2011b) modified it by taking into account possible variation of the mean value of the virial density threshold with mass scale. The modified formula, which has one additional parameter called the drifting average coefficient, was shown to improve the agreements with the numerical results on the galactic mass scale.

It was Ryu & Lee (2021) who for the first time proposed that the generalized excursion set theory can also be used for the evaluation of the splashback mass function. From here on, we let $M$ denote not the virial but the splashback mass of a DM halo. They put forth an idea that the splashback density threshold must be much less stochastic than the virial counterpart and that its mean value would not change with the mass scale. Then, they claimed that the splashback multiplicity function could be described by the following single-parameter formula, which is the same as the Maggiore-Riotto model but with $\delta_{sc}=1.686$ replaced by $\delta_{sp}=1.52$, the mean splashback density threshold predicted in the self-similar spherical infall model (Shapiro et al., 1999):

Here, $\Gamma$ is the incomplete gamma function, $D_{\rm B}$ is the diffusion coefficient defined as the ratio of the variance of the splashback density threshold to the mass variance, and $\kappa$ is related to the cross correlation of the linear density field on two different splashback mass scales $M$ and $M^{\prime}$ as¹¹1In the work of Ryu & Lee (2021), they adopted the same approximation of $\kappa\approx 0.475$ as used in Corasaniti, & Achitouv (2011b), assuming a $\Lambda$CDM cosmology. However, in the current work, we rigorously calculate $\kappa$ by Equation (3).:

Comparing Equations (1)-(3) with the numerical results from the Erebos $N$-body experiment (Diemer & Kravtsov, 2014, 2015) and determining $D_{\rm B}$ as a function of redshift, Ryu & Lee (2021) found the following. First, despite that it has only one parameter, the analytic formula excellently works in the splashback mass range of $5\times 10^{12}\leq M/(h^{-1}\,M_{\odot})\leq 10^{15}$ up to $z\sim 3$ for two different $\Lambda$CDM cosmologies. Second, the splashback mass function is characterized by a much lower value of $D_{\rm B}$ than the virial mass function, which justifies the assumption that the splashback density threshold is less stochastic than the virial counterpart.

Third, the diffusion coefficient almost monotonically decreases with $z$, well approximated by a linear fit, and vanishing at some critical redshift, $z_{c}$. This result implied that the splashback mass function at $z\geq z_{c}$ can be described by the following purely analytic formula with no free parameter:

It was also found by Ryu & Lee (2021) that the critical redshift, $z_{c}$, appears to sensitively depend on the initial conditions, exhibiting a tendency to have a lower value for the case of the less amount of dark matter.

We are going to numerically tackle the following two issues. First, is the analytic model reviewed in Section 2 still valid for the $\nu\Lambda{\rm CDM}$ cosmology? Second, what effect does the presence of neutrinos has on the evolution of the diffusion coefficient, $D_{\rm B}(z)$? To address these issues, we utilize the snapshot data from the Cosmological Massive Neutrino Simulations (MassiveNuS) (Liu et al., 2018) for which the effect of neutrinos was incorporated into the background with the help of the analytic linear response approximation (Ali-Haïmoud & Bird, 2013). In the periodic box of linear size $512\,h^{-1}$Mpc, the MassiveNuS contained a total of $N_{\rm dm}=1024^{3}$ DM particles with individual mass $10^{10}\,h^{-1}\,M_{\odot}$ (Liu et al., 2018).

The snapshot data from the MassiveNuS are available only for two $\nu\Lambda{\rm CDM}$ models which share the same matter density parameter ($\Omega_{m}$) and same large-scale amplitude of the linear density power spectrum ($A_{s}$) but different total neutrino mass, $\sum m_{\nu}/\,{\rm eV}=0.0$ and $0.1$, respectively. Table 1 lists the values of $\Omega_{m}$, $A_{s}$ and $\sigma_{8}$, for the two models. Note that $\sigma_{8}$ has a slightly lower value in the presence of massive neutrinos due to the suppression of the small-scale powers caused by the neutrino free streaming (Lesgourgues & Pastor, 2014). The catalogs of distinct halos and their subhalos identified by the Rockstar algorithm (Behroozi et al., 2013) in the snapshot data are also publicly available at the MassiveNuS website ²²2http://astronomy.nmsu.edu/aklypin/SUsimulations/MassiveNuS/.

We apply the publicly available SPARTA (Subhalo and PARticle Trajectory Analysis) algorithm (Diemer, 2017) to a total of $45$ snapshots and rockstar catalogs in the redshift range of $0\leq z\leq 5.285$ from the MassiveNuS for the two $\nu\Lambda{\rm CDM}$ cosmologies. The redshift interval, $\Delta z$, between the adjacent snapshots in the redshift range of $z\leq 1$ is approximately $\delta z\approx 0.04$. The SPARTA basically inspects the orbits of all particles in each rockstar halo over redshifts to compute the distances to the apocenters of their first orbits, $r_{sp}$, and takes the average over them, $R^{mn}_{sp}$, to define the splashback radius of each halo. The SPARTA algorithm also provides several different definitions of the splashback radius, $R^{50}_{sp}$, $R^{70}_{sp}$, and $R^{90}_{sp}$, which correspond to the $50\%$, $70\%$, and $90\%$ percentiles of $r_{sp}$, respectively. The analytic formula reviewed in Section 2 was proven to describe well the splashback mass function for the case of the $\Lambda$CDM cosmology, no matter which definition of $R_{sp}$ was used. In the current work, we exclusively use $R^{mn}_{sp}$ for the halo splashback radius, since the linear fit approximation to the evolution of the diffusion coefficient turned out to work best for the case of $R^{mn}_{sp}$ (Ryu & Lee, 2021).

To determine the splashback mass function, we consider only the distinct halos that are not embedded within the splashback radii of any larger halos. Note that the criterion for a distinct halo depends on the halo finder algorithm. For example, some rockstar object classified as a distinct halo by the spherical overdensity (SO) code could be classified as a subhalo by the MORIA code, a subrouine contained in the SPARTA algorithm, as its splashback radius can be overlapped with that of a neighbor larger halo. Selecting the distinct halos with splashback masses $M\geq 5\times 10^{12}\,\,h^{-1}\!M_{\odot}$, we split the range of $\ln M$ into short bins with equal length of $d\ln M$ and count the number of the selected halos, $dN$, whose logarithmic masses are in the range of $\left[\ln M,\ln M+d\ln M\right]$. Then, the splashback mass function is numerically determined as the ratio, $dN/d\ln M$, at each redshift, to which the analytic model, Equations (1)-(3), is adjusted by varying $D_{\rm B}$. Dividing the halo sample into eight subsamples of equal size, we also repeat the same calculation separately for each of the eight subsamples to obtain eight different splashback mass functions, $\{dN_{i}/d\ln M\}_{i=1}^{8}$. The one standard deviation scatter among $\{dN_{i}/d\ln M\}_{i=1}^{8}$ is then determined as the Jackknife errors in the original splashback mass function. As done in Ryu & Lee (2021), the standard $\chi^{2}$-statistics is employed to find the best-fit value of $D_{\rm B}$, while the error in $D_{\rm B}$ is computed as the associated Fisher information.

Figure 1 plots the numerically obtained splashback mass function (black closed circles) with the Jackknife errors at $z=0$, and compares them with the analytic formula (red solid line) with the best-fit value of $D_{\rm B}$, for the two $\nu\Lambda{\rm CDM}$ models in the first and third from the top panels. The ratios of the numerical result to the analytical formula for each model (red solid line) with one standard deviation scatter (gray area) are also shown in the third and first from the bottom panels. As can be seen, for both of the $\nu\Lambda{\rm CDM}$ models, despite that it has only one parameter, the analytic formula excellently matches the numerical results in the wide mass range of $5\times 10^{12}\leq M/(h^{-1}\,M_{\odot})\leq 10^{15}$ at $z=0$.

Repeating the same calculation but at higher redshifts, we trace $D_{\rm B}(z)$. Figures 2-3 plot the same as Figure 1 but at $z=0.47$ and $1.05$, respectively, demonstrating how well the analytic formula, Equations (1)-(3) with fixed $\delta_{sp}=1.52$, describes the splashback mass function even in the presence of neutrinos over a broad range of $z$. Figure 4 shows how $D_{\rm B}(z)$ (black filled circles) changes with redshifts for the two cases of $\sum m_{\nu}$. As can be seen, for both of the cases, $D_{\rm B}(z)$ decreases almost linearly with redshifts at $z\geq 0.2$, showing a trend of converging to zero at some critical redshift, $z_{c}$. If the diffusion coefficient vanishes, the analytic formula of the splashback mass function becomes parameter free. Figure 5 compares the numerical results at three edshifts beyond $z_{c}$ with the parameter-free model, Equation (4). Although the numerical results suffer from large Jackknife errors at these high redshifts due to the low number densities of massive cluster halos, the splashback mass functions agree quite well with the purely analytic model for both of the $\nu\Lambda{\rm CDM}$ models at $z>z_{c}$. These result are consistent with what Ryu & Lee (2021) found for the case of the $\Lambda$CDM cosmology, confirming the validity of Equations (1)-(4) even in the presence of neutrinos.

As done in Ryu & Lee (2021), approximating $D_{\rm B}(z)$ in the range of $0.2\leq z\leq z_{c}$ to a linear fit, $B(z_{c}-z)$, we determine the best values of $B$ and $z_{c}$ via the standard $\chi^{2}$-statistics, which are listed in Table 1 . As can be read, the two models differ in the best-fit values of $B$ and $z_{c}$, the massive neutrino case yielding a substantially larger value of $z_{c}$ and significantly lower value of $B$ than the massless neutrino case. This result indicates that the presence of massive neutrinos has an effect of making $D_{\rm B}(z)$ more diffusive in this redshift range, slowing down the decrease of $D_{\rm B}(z)$ with $z$. In other words, in the presence of hot DM particles like massive neutrinos that have large free streaming lengths, the DM halos experience larger amount of disturbance from the surrounding, which make it more ambiguous to precisely identify their physical boundaries.

Unlike the result of Ryu & Lee (2021) that $D_{\rm B}(z)$ linearly decreases with redshifts in the whole range from $z=0$ to $z_{c}$, however, we find an anomalous behavior of $D_{\rm B}(z)$ at $z<0.2$ in the presence of neutrinos: it increases with $z$, deviating from the linear fit at $z<0.2$, as can be seen in Figure 4. We also note that the massive neutrino case exhibits a much more conspicuous deviation of $D_{\rm B}(z)$ from the linear fit than the massless case at these low redshifts. We suspect that it should be caused or at least linked closely with the inherent limitation of the SPARTA code in locating the splashback boundaries of DM halos at $z\leq 0.13$ where the unknown future orbits of infalling particles make it impossible to determine particle apocenters (private communication with B. Diemer 2022).

We have numerically investigated the effect of neutrinos on the splashback mass functions of DM halos in the mass range of $0.5\leq M/(10^{13}\,h^{-1}\,M_{\odot})\leq 10^{2}$ by applying the SPARTA code (Diemer, 2017) to the snapshots of MassiveNuS (Liu et al., 2018) for two different $\nu\Lambda{\rm CDM}$ cosmologies with total neutrino mass of $\sum m_{\nu}=0.0\,{\rm eV}$ and $0.1\,{\rm eV}$. Comparing the numerically obtained splashback mass functions to the analytic single parameter formula proposed by Ryu & Lee (2021) at various redshifts, we have trailed the evolution of the single parameter, dubbed the diffusion coefficient $D_{\rm B}(z)$, in the range of $0\leq z\leq 3$. It has turned out that the analytic formula with fixed splashback density threshold $\delta_{sp}=1.52$ validly describes the splashback mass function for both of the $\nu\Lambda{\rm CDM}$ models at redshifts up to $z\approx 3$ and that the diffusion coefficient decreases almost linearly with redshifts at $z\geq 0.2$, evanescing at some critical redshift, which are all in consistent with the claim of Ryu & Lee (2021).

Fitting $D_{\rm B}(z)$ at $z\geq 0.2$ to a linear scaling relation of $B(z-z_{c})$ for both of the $\nu\Lambda{\rm CDM}$ models, we have determined the best-fit values of the slope, $B$, and critical redshift, $z_{c}$, and find that they sensitively depend on the total neutrino mass. The $\nu\Lambda{\rm CDM}$ model with $\sum m_{\nu}=0.1\,{\rm eV}$ yields a significantly lower value of $B$ and a substantially higher value of $z_{c}$ than the other model with $\sum m_{\nu}=0.1\,{\rm eV}$, which indicates that the diffusion coefficient decreases more slowly with $z$ in the presence of massive neutrinos. Since the diffusion coefficient of the analytic formula is a measure of the stochasticity of the splashback density threshold, this result implies that the presence of hot DM particles like massive neutrinos has an effect of disturbing more severely the splashback boundaries and making the halo identification more ambiguous. The dependence of $B$ and $z_{c}$ on $\sum m_{\nu}$ also implies that the diffusion coefficient of the splashback mass function can in principle be used as a probe of the total neutrino mass.

A failure of the linear fit approximation to $D_{\rm B}(z)$ has been witnessed for both of the models at $z<0.2$ when the diffusion coefficient exhibit increment rather than decrement with $z$, an anomalous tendency never detected in the absence of neutrinos for which case the linear fit always matches the diffusion coefficient in the entire range of $0\leq z<3$. Given the fundamental limitation of the SPARTA code at $z\leq 0.13$ that it is incapable of locating the apocenters of infalling DM particles whose future orbits are known (private communication with B. Diemer 2022), we suspect that the anomalous behavior of $D_{\rm B}(z)$ witnessed at $z<0.2$ should not be a real one caused by the presence of massive neutrinos but a spurious one caused by the inherent flaw of the SPARTA code at low redshifts.

It is, however, worth discussing a caveat that our conclusion is subject to. The MassiveNuS employed the analytic linear response approximation (Ali-Haïmoud & Bird, 2013) to mimic the presence of massive neutrinos rather than directly including the massive neutrinos as component particles. Although this approximation has been known to work quite well provided that $\sum m_{\nu}\leq 0.3\,{\rm eV}$ (Bird et al., 2018), it is not guaranteed to accommodate all possible effects that the massive neutrinos would have on the nonlinear evolution of cosmic web. It will require a more rigorous and accurate treatment of the neutrino effects to investigate the true behavior $D_{\rm B}(z)$ in the nonlinear regime and its dependence on the total neutrino mass. It will be also quite desirable to have an analytic model for $D_{\rm B}(z)$ derived from the first principles, with which one can quantitatively probe the nature of dark matter including massive neutrinos. Our future work is in this direction.