Constraining the Neutrino Mass with the Drifting Coefficient of the Field Cluster Mass Function

The galaxy clusters are often divided into two categories, the wall and the field clusters to which the members of the superclusters and the rest correspond, respectively. Presuming that the field clusters are more isolated and thus less susceptible to the disturbing effects of the surrounding cosmic web than their wall counterparts, Lee (2012) claimed that the field clusters should provide a more sensitive indicator of the background cosmology, developing an analytic model for the field clusters in the generalized excursion set framework (Maggiore & Riotto, 2010a, b; Corasaniti & Achitouv, 2011a, b). The analytic model of Lee (2012) was indeed a practical and theoretical success. On the practical side, its success was demonstrated by the excellent agreements with the N-body results obtained assuming a base cosmology, where the late-time cosmic acceleration is caused by the cosmological constant ($\Lambda$) and the structure formation is predominantly driven by the gravity of the cold dark matter (CDM).

On the theoretical side, its success resonates with the fact that it is a physical model having only one deterministic parameter. Unlike the conventionally used empirical formulae with physically meaningless multiple stochastic parameters for the cluster mass function (e.g., Sheth & Tormen, 1999; Warren et al., 2006; Tinker et al., 2008), the analytic model of Lee (2012) for the field cluster mass function imparts physical substance to its single parameter. Although it is inevitable to resort to the numerical experiments for the determination of its exact value, this single parameter dubbed drifting coefficient (Corasaniti & Achitouv, 2011a) is a physical measure of how far the density threshold for the realistic non-spherical collapse, $\delta_{c}$, drifts from that for the idealistic spherical collapse, $\delta_{sc}$, at a given mass scale.

Moreover, the recent work of Ryu & Lee (2020) discovered that the drifting coefficient of the field cluster mass function in fact carries significance beyond a physical parameter quantifying the deviation of $\delta_{c}$ from $\delta_{sc}$. Confirming that the validity of the analytic model of Lee (2012) for the field cluster mass functions is robust even against the variations of the key cosmological parameters including the dark energy (DE) equation of state, Ryu & Lee (2020) found that the drifting coefficient evolves differently even among those degenerate dynamical DE cosmologies which yield almost the same linear growth factors or the same cluster mass functions. In the light of the results of Ryu & Lee (2020), we speculate that the evolution of the drifting coefficient might be also useful to constrain the total mass of neutrino species, $M_{\nu}\equiv\sum m_{\nu}$, one of the utmost missions entrusted to the cosmological physics (Lesgourgues & Pastor, 2012).

The latest Planck analysis of the Cosmic Microwave Background (CMB) temperature power spectra combined with the priors from the weak gravitational lensing (WL) and Baryonic Acoustic Oscillations (BAO) concluded $M_{\nu}\leq 0.12\,{\rm eV}$ (Planck Collaboration et al., 2018), assuming the base flat $\Lambda$CDM cosmology (see also Vagnozzi et al., 2017). A higher value of $M_{\nu}$ above $0.12\,{\rm eV}$, however, can still be accommodated by the Planck data, if the assumption about the background cosmology is released (see Choudhury & Choubey, 2018; Choudhury & Hannestad, 2019, and references therein) or if different priors are used to complement the CMB probe (e.g., Giusarma et al., 2016). For the past decade, the cluster mass function has been prevalently promoted as an useful complementary probe of $M_{\nu}$ (e.g., Marulli et al., 2011; Ichiki & Takada, 2012; Costanzi et al., 2013; Villaescusa-Navarro et al., 2013; Castorina et al., 2014; Biswas et al., 2019; Hagstotz et al., 2019). Although the cluster mass function is only indirectly linked to $M_{\nu}$ through its dependence on the linear density power spectrum, it has a practical advantage as a probe of $M_{\nu}$, being more readily observable than the linear density power spectrum, the measurements of which are often plagued by the systematics stemmed from the existence of nonlinear galaxy bias (Giusarma et al., 2018, and references therein).

Due to the inherent non-sphericity and stochastic aspect of the cluster formation process that defies purely analytic modeling from the first principle, a theoretical prediction for the cluster abundance and its dependence on $M_{\nu}$ was conventionally made in the empirically modified excursion set formalism (e.g., Costanzi et al., 2013; Villaescusa-Navarro et al., 2013; Biswas et al., 2019). While a link between the cluster abundance and $M_{\nu}$ through the linear power spectrum is provided by the excursion set theory, the required accuracy and precision was achieved by the empirical modification of the theory, i.e., deteriorating of a physical model into a fitting formula with multiple free parameters (Warren et al., 2006; Tinker et al., 2008). Lack of a physical model for the cluster abundance undermines its power as a probe of $M_{\nu}$. To make matters worse, the notorious $\sigma_{8}$-$M_{\nu}$ degeneracy of the initial density power spectrum translates into the relative low sensitivity of the cluster mass function to $M_{\nu}$.

Given the aforementioned difficulties in constraining $M_{\nu}$ with the cluster abundance, what may be desirable to have is a new probe, well described by a physical model, free from the $\sigma_{8}$-$M_{\nu}$ degeneracy, and highly sensitive to the variation of $M_{\nu}$. Our goal here is to prove that the drifting coefficient of the field cluster mass function fulfills this expectation. Section 2 will be consigned to a brief review of the works of Lee (2012) and Ryu & Lee (2020). Section 3 will present a procedure through which the power and efficacy of the drifting coefficient as a new probe of $M_{\nu}$ is numerically appraised. Section 4 will be devoted to discussing a physical implication of the final results and a prospect for constraining $M_{\nu}$ with this new probe in practice, as well.

The differential mass function of the field clusters, $dN_{\rm I}/d\ln M$, gives their number densities in each logarithmic interval of $[\ln M,\ln M+d\ln M]$. Suggesting for the first time that the mass function of the field clusters should be more sensitive to the background cosmology than that of all clusters, Lee (2012) modified the generalized excursion set mass function theory (Maggiore & Riotto, 2010a, b; Corasaniti & Achitouv, 2011a, b) to derive the following single parameter model for $dN_{\rm I}/d\ln M$, which will be adopted for our analysis.

where $\kappa=0.475$ and $\beta$ is the drifting coefficient that quantifies how much the non-sphericity of the gravitational collapse drifts the density threshold $\delta_{c}$ away from the Einstein de Sitter spherical collapse threshold of $\delta_{sc}=1.686$ (Gunn, & Gott, 1972; Peebles, 1993) at a given mass scale. Since the non-spherical gravitational collapse process is too complicated for $\delta_{c}$ to be theoretically predicted from the first principle (Bond & Myers, 1996), $\beta$ has to be treated as a free adjustable parameter, as in the generalized excursion set formalism (Corasaniti & Achitouv, 2011a, b). Nevertheless, as Lee (2012) and Ryu & Lee (2020) explained, the density threshold $\delta_{c}$ (or equivalently $\beta$) is deterministic for the field clusters, while it is stochastic fo their wall counterparts (Robertson et al., 2009; Maggiore & Riotto, 2010a, b), which allows the field cluster mass function to have only one free parameter in the generalized excursion set formalism.

This analytical single parameter model, Equations (1)-(5), connects $dN_{\rm I}/d\ln M$ to the initial conditions of the universe through two different routes. The rms density fluctuation of the initial density field, $\sigma(M,z)$, expressed in terms of the linear density power spectrum, $P(k,z)$, is the usual route envisaged by the original excursion set mass function theory (Press & Schechter, 1974; Bond et al., 1991). While, the drifting coefficient, $\beta(z)$, is another independent route induced by the cosmology dependence of $\delta_{c}$ (Ryu & Lee, 2020). To effectively describe different behaviors of $\beta(z)$ among different DE cosmologies, the following fitting formula was proposed by Ryu & Lee (2020),

with three fitting parameters, $\beta_{A},\ q_{z}$ and $z_{c}$.

As mentioned in Section 1, Ryu & Lee (2020) tested this fitting formula for $\beta(z)$ as well as the above analytical single parameter model for $dN_{\rm I}/d\ln M$ against the large N-body simulations for various DE cosmologies including the $\Lambda$CDM and confirmed that it is quite valid regardless of the DE equation of states. Moreover, it was also shown by Ryu & Lee (2020) that $\beta(z)$, via Equation (6), allows us to distinguish even among those degenerate DE cosmologies which produce almost the same linear density power spectra and cluster mass functions. In Section 3, we are going to numerically test if Equations (1)-(5) are also valid for the $\nu\Lambda$CDM cosmology and to examine its power as a complementary probe of $M_{\nu}$.

Before proceeding further, it is worth mentioning that there is actually another route other than $P(k,z)$ and $\delta_{c}$ that connects $dN_{\rm I}/d\ln M$ to the background cosmology. This third route is nothing but the spherical collapse density threshold, $\delta_{sc}$, for which even purely theoretical predictions from the first principles can be made thanks to the spherical symmetry. The cosmology dependence of $\delta_{sc}$, however, was found too weak to stand out over those of $P(k,z)$ and $\delta_{c}$ (e.g., see Eke et al., 1996; Pace et al., 2010), which is why the connection between $dN_{\rm I}/d\ln M$ and the initial conditions can be almost entirely attributed to $P(k,z)$ and $\delta_{c}$. The same argument applies to the $M_{\nu}$-dependence of $\delta_{sc}$, which was already shown to be not so strong as that of $P(k,z)$ (e.g., LoVerde, 2014). Throughout this Letter as in Ryu & Lee (2020), we set the spherical density threshold $\delta_{sc}$ at the fixed Einstein-de Sitter value, $1.686$ (see also Maggiore & Riotto, 2010a, b; Lee, 2012).

We make an extensive use of the publicly available data from the Cosmological Massive Neutrinos Simulations (MassiveNuS) run by Liu et al. (2018) on a periodic box of comoving volume $512^{3}\,h^{-3}\,$Mpc³, containing $1024^{3}$ particles, each of which is as massive as $10^{10}\,h^{-1}\,M_{\odot}$. The MassiveNuS was recurringly performed for one $\Lambda$CDM cosmology with massless neutrinos and for $100$ different $\nu\Lambda$CDM cosmologies with massive neutrinos, whose initial conditions were described by the six key cosmological parameters as well as $M_{\nu}$. For the study of the sole effect of the massive neutrinos on $dN_{\rm I}/d\ln M$ and $\beta(z)$, we consider only those cosmologies which have identical initial conditions other than $M_{\nu}$ with one another. Among the $101$ cosmologies are found only $3$ to meet this selection criterion, which have the same matter density parameter, $\Omega_{m}=0.3$, and same amplitude of the primordial density power spectrum, $A_{s}=2.1\times 10^{9}$, but different total neutrino mass, $M_{\nu}=0.0,\ 0.1$ and $0.6$ eV, respectively.

The MassiveNuS engaged the Rockstar algorithm (Behroozi et al., 2013) to find the DM halos at various redshifts and recorded such key properties of each Rockstar halo as its virial mass ($M$), virial radius, comoving position vector, peculiar velocity vector and so forth. From the catalog of the Rockstar halos resolved at each redshift for each of the three cosmologies, we first exclude the subhalos embedded in larger parent halos and then set the cutoff mass at $3\times 10^{13}\,h^{-1}\,M_{\odot}$ to sort out the distinct cluster halos. Following the same procedure arranged in Ryu & Lee (2020), we apply the friends-of-friends (FoF) algorithm with the linkage length parameter of $l_{\rm c}=0.33$ to the distinct cluster halos for the identification of the superclusters composed of two or more members. Eliminating the wall clusters belonging to the identified superclusters, we end up having a sample of the distinct field cluster halos with $M\geq 3\times 10^{13}\,h^{-1}\,M_{\odot}$. Then, we reckon the field clusters at each logarithmic mass bin to numerically determine $dN_{I}/d\ln M$ to which the analytical single parameter model, Equations (1)-(5), is fitted by adjusting the value of $\beta$.

In the procedure of evaluating the analytic mass functions of the field clusters, the CAMB code (Lewis et al., 2000) is exclusively used for $P(k,z)$, while the standard $\chi^{2}$-statistics is employed for the best-fit value of $\beta$. Note that since the linear growth factor, $D(z)$, acquires a scale dependence in the presence of massive neutrinos, the rms density fluctuation $\sigma(M,z)$, is no longer equal to $D(z)\sigma(M,z=0)$. Instead, we calculate it as $\sigma(M,z)=\left[(2\pi^{2})^{-1}\int dk\,k^{2}P(k,z)W^{2}_{\rm th}(k,M)\right]^{1/2}$ where $W_{\rm th}$ is the spherical top-hat filter on the mass scale of $M$.

Figure 1 plots the linear density power spectra, $P(k,z)$, for the three different cases of $M_{\nu}$ at three different redshifts, computed by the CAMB code. As expected, the more massive neutrinos suppress more severely the linear density powers on the small scales ($k>0.02\,h\,{\rm Mpc}^{-1}$). Note the small differences in $P(k,z)$ between the cases of $M_{\nu}=0.0$ eV and $M_{\nu}=0.1$ eV at all of the three redshifts. Given that the large uncertainties in the high-mass tails of the cluster mass functions caused by poor-number statistics and cosmic variance are likely to exceed this small differences in $P(k,z)$, the cluster mass functions would be unable to discriminate the two $\nu\Lambda$CDM cosmologies from each other.

Figure 2 displays both of the numerical field cluster mass functions from the MassiveNuS (filled circles) and the analytic model with the best-fit value of $\beta$ (red solid lines) at $z=0$ for the three different cases of $M_{\nu}$. The errors in the numerical determination of $dN_{\rm I}/d\ln M$ is calculated as one standard deviation from the mean averaged over eight Jackknife resamples, (Ryu & Lee, 2020). The black dotted lines in the middle and right panels conform to the red solid line in the left panel. Figures 3-4 show the same as Figure 2 but at $z=0.42$ and $0.83$, respectively. As can be seen, the analytical single parameter model for $dN_{\rm I}/d\ln M$ agrees excellently well with the numerical results at all redshifts for all of the three cases of $M_{\nu}$, confirming its validity even in the presence of massive neutrinos and proving its robustness as a physical model.

Figures 2-4 clearly show that $dN_{\rm I}/d\ln M$ has a significantly lower amplitude for the case of $M_{\nu}=0.6\,{\rm eV}$ than for the other two cases of $M_{\nu}=0.0\,{\rm eV}$ and $M_{\nu}=0.1\,{\rm eV}$, between which almost no difference is found in $dN_{\rm I}/d\ln M$, no matter at what redshifts they are compared with each other. Although the difference in $dN/d\ln M$ between the two cases of $M_{\nu}=0.0\,{\rm eV}$ and $M_{\nu}=0.1\,{\rm eV}$ tends to slightly increase with $z$, the larger errors in the measurement of $dN_{\rm I}/d\ln M$ at higher redshifts weigh down their statistical significances. The comparison of Figures 2-4 with Figure 1 indicates that the $M_{\nu}$-dependence of the field cluster mass function is almost entirely dictated by the $M_{\nu}$-dependence of $P(k,z)$. As mentioned in Section 2, the cosmology-dependence of the field cluster abundance (including its $M_{\nu}$-dependence) has two different sources, $P(k,z)$ and $\beta$. The results shown in Figures 2-4, however, imply that the former overwhelms the latter in shaping the $M_{\nu}$-dependence of the field cluster mass function, which in turn warns that the field cluster mass function would fail not only in constraining $M_{\nu}$ below the Planck constraint but also in breaking the $\sigma_{8}$-$M_{\nu}$ degeneracy.

Figure 5 plots the numerically determined values of $\beta(z)$ at twenty different redshifts in the range of $0\leq z\leq 1$ for the three different cases of $M_{\nu}$, revealing that $\beta(z)$ evolves differently among the three cases. Here, the errors, $\sigma_{\beta}$, are obtained through the Fisher information analysis, as done in Ryu & Lee (2020). As can be seen, at $z\lesssim 0.3$ the drifting coefficient $\beta(z)$ has higher values for the case of $M_{\nu}=0.6\,{\rm eV}$ than for the other two cases. Whereas at $z\gtrsim 0.3$, the tendency is reversed. The most massive neutrinos case yields the lowest values of $\beta(z)$, while its highest values are found for the massless neutrinos case. In addition, we find that the slope of $\beta(z)$ substantially differs even between the two cases of $M_{\nu}=0.0\,{\rm eV}$ and $M_{\nu}=0.1\,{\rm eV}$, while no difference found in $\beta(z=0)$ between them.

We speculate that this redshift-dependence of the effect of massive neutrinos on $\beta(z)$ might help break the $\sigma_{8}$-$M_{\nu}$ degeneracy. Recall that the effect of massive neutrinos on the linear density power spectra and cluster mass function is consistent in its direction, regardless of the redshifts, as witnessed in Figures 1-4. The more massive neutrinos always reduce more severely the amplitudes of $P(k,z)$ and $dN_{\rm I}/d\ln M$ at all redshifts, which is why the two diagnostics suffer from the $\sigma_{8}$-$M_{\nu}$ degeneracy. In other words, the lower value of $\sigma_{8}$ has the same effect on $P(k,z)$ (and $dN/d\ln M$ as well) as the higher value of $M_{\nu}$. Meanwhile, our result shown in Figure 5 implies that the effect of the higher value of $M_{\nu}$ on $\beta(z)$ might be differentiated from that of the lower value of $\sigma_{8}$ on $\beta(z)$. The latter lowers the amplitude of $\beta(z)$ without changing its slope, while the former heightens its amplitude and concurrently steepens its slope. Yet, the possibility of breaking the $\sigma_{8}$-$\Omega_{m}$ degeneracy with $\beta(z)$ is only a speculation, since we have yet to demonstrate its feasibility in practice.

As done in Ryu & Lee (2020), to effectively quantify the differences in the evolution of the drifting coefficient among the three cosmologies, we fit Equation (6) to the numerically determined $\beta(z)$ by adjusting the values of $\beta_{A},\ q_{z}$ and $z_{c}$ to yield the minimum $\chi^{2}$. Figure 6 demonstrates how well the simple formula (red solid lines), Equation (6), suggested by Ryu & Lee (2020), agrees with the numerically obtained $\beta(z)$ (black filled circles) for all of the three cases of $M_{\nu}$. Figure 7 shows the best-fit values of $-\beta_{A},\ q_{z},\ z_{c}$ with their errors $\sigma_{\beta_{A}},\ \sigma_{qz},\ \sigma_{zc}$, which are all obtained through the $\chi^{2}$ fitting after due consideration of the uncertainties in $\beta(z)$ shown in Figure 6.

The most significant differences among the three cases are found in the values of $z_{c}$, which is consistent with the result of Ryu & Lee (2020) that $z_{c}$ was found to vary most sensitively with the dark energy equation of state. Assessing the statistical significances of the differences in $z_{c}$ among the three cases of $M_{\nu}$ by estimating the errors of their mutual differences, $\sigma_{\Delta(zc)}$, propagated from $\sigma_{zc}$, as done in Ryu & Lee (2020), we find the difference in $z_{c}$ between the two cases of $M_{\nu}=0.0\,{\rm eV}$ and $M_{\nu}=0.1\,{\rm eV}$ ($M_{\nu}=0.6\,{\rm eV}$) to exceed $4\sigma_{\Delta({zc})}$ ($10\sigma_{\Delta(zc)}$). Whereas, the differences in the other two parameters, $\beta_{A}$ and $q_{z}$, between the two cases of $M_{\nu}=0.0\,{\rm eV}$ and $M_{\nu}=0.1\,{\rm eV}$ ($M_{\nu}=0.6\,{\rm eV}$) are found to be statistically insignificant (not so significant as that in $z_{c}$).

It should be worth explaining here why $z_{c}$ is the most sensitive to the variation of $M_{\nu}$. Given the definition $z_{c}$ as a critical redshift at which $\delta_{c}=1.686$ (i.e, $\beta(z_{c})=0$), its value should be determined by two factors, both of which sensitively depend on $M_{\nu}$. The first factor is how fast the matter density parameter $\Omega_{m}$ approaches unity (i.e, the Einstein-de Sitter value) at high redshifts, while the second one is how rare the field clusters are in a given universe, since the gravitational collapse of the rarer objects proceeds in a more spherically symmetrical way (Bernardeau, 1994). Meanwhile, the other two parameters, $\beta_{A}$ and $q_{z}$, depend mainly on either of the two factors: $\beta_{A}$ on the second, while $q_{z}$ on the first.

Conducting a numerical analysis of the MassiveNuS data (Liu et al., 2018), we have found that the massive neutrinos have a unique redshift-dependent effect on the drifting coefficient of the field cluster mass function, $\beta(z)$, which measures the difference between the density thresholds for the realistic nonspherical and the idealistic EdS spherical collapse at a given mass scale. In our previous work (Ryu & Lee, 2020), we already found that $\beta$ vanishes to zero at a certain critical redshift $z_{c}$ but increases as the universe evolves from $z_{c}$ to $z=0$ like an inverse sine hyperbolic function of $z$. We have newly found here that the presence of more massive neutrinos lowers $z_{c}$ and induce a faster increase of $\beta(z)$ with the decrement of $z$ below $z_{c}$. The $\nu\Lambda$CDM cosmology with total neutrino mass of $M_{\nu}=0.6\,{\rm eV}$ has been found to yield higher (lower) values of $\beta$ at $0\leq z\lesssim z_{\rm th}$ ($z_{\rm th}\lesssim z\leq z_{c}$) than the $\Lambda$CDM cosmology with massless neutrinos with $z_{\rm th}\sim 0.3$. Noting that this redshift-dependent effect of massive neutrinos on $\beta$ is quite unique and distinct especially from the redshift-independent effect of $\sigma_{8}$ on $\beta(z)$, we suggest that the drifting coefficient of the field cluster mass function should allow us to break the notorious $\sigma_{8}$-$M_{\nu}$ degeneracy, which has haunted for long the conventional probes of $M_{\nu}$ based on the linear density power spectrum.

Our physical explanation for this distinct redshift-dependent effect of $M_{\nu}$ on $\beta(z)$ is that it is generated by a competition between the suppressed small-scale powers and the increased degree of the anisotropy of the cosmic web in the presence of massive neutrinos. As shown by Bernardeau (1994), the formation of a rare event like a massive cluster (or a field cluster) is well approximated by a spherical collapse process. The rarer an object is, the more spherically its gravitational collapse proceeds. In the presence of more massive neutrinos which suppress more severely the small-scale powers, a field cluster corresponds to an even rarer object since it originates from a more extreme local maximum in the initial density field. Therefore, it is naturally expected that in the presence of more massive neutrinos the collapse density threshold $\delta_{c}$ for the field clusters would become closer to the spherical threshold $\delta_{sc}$ (or equivalently, $\beta$ closer to zero).

The free streaming of massive neutrinos, however, has another effect of rendering the cosmic web more anisotropic in the deeply nonlinear stage. According to the previous works (e.g., Shim et al., 2014; Ho et al., 2018) which found the degree of the anisotropy of the cosmic web to depend on the background cosmology, the stronger gravity at a given scale pulls down the anisotropic feature of the cosmic web in the nonlinear stage. The free streaming of massive neutrinos plays a role along with DE in weakening the gravitational clustering on the cluster scale, which in consequence increases the degree of the anisotropy of the cosmic web. The stronger tidal influences from the more anisotropic cosmic web (Bond et al., 1996) deviate the collapse process further from the spherical symmetry, elevating $\beta$ above zero.

At high redshifts ($z_{\rm th}\lesssim z\leq z_{c}$), the first effect of massive neutrinos overwhelms the second, lowering $\beta$ close to zero, since the high-$z$ field clusters correspond to the rarest events formed through the collapses of the highest density peaks which proceed in almost perfectly spherically. However, at lower redshifts ($0\leq z\lesssim z_{\rm th}$) after the onset of the nonlinear evolution of the cosmic web, the second effect wins over the first, deviating $\beta$ further from zero. Our result shown in Figure 6 reveals the $M_{\nu}$-dependence of the threshold redshift, $z_{\rm th}$, at which the second effect becomes more dominant than the first. It is around $0.3$ for the case of $M_{\nu}=0.6\,{\rm eV}$, while it becomes around zero for the case of $M_{\nu}=0.1\,{\rm eV}$. The more massive neutrinos induce the turn-over of the second effect to occur earlier. Our future work is in the direction of constructing a more theoretical model for $\beta(z)$, within which the $M_{\nu}$-dependences of $z_{\rm th}$ and $z_{c}$ could be predicted.

Another important hint of this work is that the sensitivity of $\beta(z)$ to $M_{\nu}$ might be high enough to detect the effect of massive neutrinos on it, even in case that $M_{\nu}$ is as low as $0.1\,{\rm eV}$ below the Planck constraint (Planck Collaboration et al., 2018). The signal of the difference in $z_{c}$ between the $\Lambda$CDM and $\nu\Lambda$CDM with $M_{\nu}=0.1\,{\rm eV}$ ($M_{\nu}=0.6\,{\rm eV}$) cosmologies has been found to be approximately four (ten) times higher than the propagated errors. Given that the observational data from much larger volumes than that of the MassiveNuS are already in the pipeline (e.g., Euclid Collaboration et al., 2019), we conclude that the drifting coefficient of the field cluster mass function, $\beta(z)$, has a good prospect for providing a very powerful complementary probe of $M_{\nu}$ in practice.