Virial Halo Mass Function in the Planck Cosmology

In spite of a strong dependence on the matter power spectrum, some analytical approaches have successfully predicted the number density of dark-matter halos (e.g., Press & Schechter 1974; Bond et al. 1991; Sheth & Tormen 2002). These approaches commonly relate the halos with their mass M to the linear density field smoothed at some scale $R={\left(3M/4\pi {\bar{\rho }}_{{\rm{m}}}\right)}^{1/3}$. To develop an analytic formula of the halo abundance, Press & Schechter (1974) assumed that the fraction of mass in halos of mass greater than M at redshift z is set to be twice the probability that smoothed Gaussian density fields exceed the critical threshold for spherical collapse, δ_c . In this ansatz, the number density of halos can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dM}}=\displaystyle \frac{{\bar{\rho }}_{{\rm{m}}}}{M}\displaystyle \frac{d\mathrm{ln}{\sigma }^{-1}}{{dM}}\,f(\sigma ),\end{eqnarray} \tag{ 1 }$$

where σ is the root mean square (rms) fluctuations of the linear density field smoothed with a filter encompassing this mass M. The rms is usually defined with a spherical top-hat filter,

$$\begin{eqnarray}&&{\sigma }^{2}(M,z)=\int \displaystyle \frac{{k}^{2}{dk}}{2{\pi }^{2}}\,{W}_{\mathrm{TH}}^{2}({kR}){P}_{{\rm{L}}}(k,z),\end{eqnarray} \tag{ 2 }$$

where P_L(k, z) is the linear matter power spectrum as a function of wavenumber k and redshift z, and W_TH is the Fourier transform of the real-space top-hat window function of radius R. To be specific, the top-hat window function is given by ${W}_{\mathrm{TH}}(x)=3[\sin x-x\cos x]/{x}^{3}$. Press & Schechter (1974) found that the function f, referred to as the multiplicity function, is expressed as

$$\begin{eqnarray}&&{f}_{\mathrm{PS}}=\sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{\delta }_{c}}{\sigma }\exp \left(-\displaystyle \frac{{\delta }_{c}^{2}}{2{\sigma }^{2}}\right).\end{eqnarray} \tag{ 3 }$$

The multiplicity function in Press & Schechter (1974) was later revised by introducing excursion set theory (Bond et al. 1991) and adopting the ellipsoidal collapse (Sheth & Tormen 2002). Note that previous analytic models predict that the multiplicity function is a universal function and any dependence on redshifts and cosmological models can be encapsulated in σ(M, z).

Motivated by these analytic predictions, we assume that the multiplicity function can be described by a four-parameter model,

$$\begin{eqnarray}\begin{array}{rcl}f(\sigma ,z) & = & A(z)\sqrt{\displaystyle \frac{2}{\pi }}\,\exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{a(z){\delta }_{c,z}^{2}}{{\sigma }^{2}}\right)\\ & & \times \ \left[1+{\left(\displaystyle \frac{\sqrt{a(z)}{\delta }_{c,z}}{\sigma }\right)}^{-2p(z)}\right]\\ & & \times \ {\left(\displaystyle \frac{\sqrt{a(z)}{\delta }_{c,z}}{\sigma }\right)}^{q(z)},\end{array}\end{eqnarray} \tag{ 4 }$$

where A(z), a(z), p(z), and q(z) are free parameters in the model. In Equation (4), we introduce the redshift-dependent critical overdensity for spherical collapse, δ_c,z . In a flat ΛCDM cosmology, this quantity is well approximated as (Kitayama & Suto 1996)

$$\begin{eqnarray}&&{\delta }_{c,z}=\displaystyle \frac{3\,{\left(12\pi \right)}^{2/3}}{20}\left[1+0.123\mathrm{log}{{\rm{\Omega }}}_{{\rm{m}}}(z)\right],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{m}}}(z)=\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{m}}0}{\left(1+z\right)}^{3}}{{{\rm{\Omega }}}_{{\rm{m}}0}{\left(1+z\right)}^{3}+(1-{{\rm{\Omega }}}_{{\rm{m}}0})}.\end{eqnarray} \tag{ 6 }$$

To find a best-fit model of f(σ, z) to the simulation data, we need to construct the multiplicity function from the halo catalogs. We start with a binned halo mass function, which is directly observable from the simulations. Let Δn_bin be the comoving number density of dark-matter halos in a bin size of ${\rm{\Delta }}\mathrm{log}M$ with mass ranges [M₁, M₂]. Using Equation (1), one finds

$$\begin{eqnarray}&&{\rm{\Delta }}{n}_{\mathrm{bin}}=\mathrm{ln}10\,{\displaystyle \int }_{\mathrm{log}{M}_{1}}^{\mathrm{log}{M}_{2}}\,d\mathrm{log}M\,\displaystyle \frac{{\bar{\rho }}_{{\rm{m}}}}{M}\,\displaystyle \frac{d\mathrm{log}{\sigma }^{-1}}{d\mathrm{log}M}\,f(\sigma ).\end{eqnarray} \tag{ 7 }$$

In the limit of ${\rm{\Delta }}\mathrm{log}M\to 0$, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{sim}}(\sigma ){| }_{M={M}_{\mathrm{bin}}} & = & \displaystyle \frac{1}{\mathrm{ln}10}\displaystyle \frac{{M}_{\mathrm{bin}}}{{\bar{\rho }}_{{\rm{m}}}}\displaystyle \frac{{\rm{\Delta }}{n}_{\mathrm{bin}}}{{\rm{\Delta }}\mathrm{log}M}\\ & & \times \ {\left(\displaystyle \frac{d\mathrm{log}{\sigma }^{-1}}{d\mathrm{log}M}\right)}^{-1}{| }_{M={M}_{\mathrm{bin}}},\end{array}\end{eqnarray} \tag{ 8 }$$

where M_bin is a center of the binned mass.

In Equation (8), we require the logarithmic derivative of $\mathrm{log}\sigma$ with respect to the halo mass M. This derivative is known to be affected by the size of simulation boxes (e.g., Lukić et al. 2007; Reed et al. 2007) because Fourier modes with scales beyond the box size are missed in the simulations. To account for the finite-volume effect on the estimate of f(σ), we follow an approach in Reed et al. (2007). The mass variance in a finite-volume simulation σ_loc is not equal to the global value in Equations (2). Nevertheless, we assume that the halo mass function in the finite-volume simulation, ${\left({dn}/{dM}\right)}_{\mathrm{loc}}$, can be written as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dn}}{{dM}}\right)}_{\mathrm{loc}}=\displaystyle \frac{{dn}}{{dM}}{\left.\right|}_{\sigma ={\sigma }_{\mathrm{loc}}}.\end{eqnarray} \tag{ 9 }$$

This ansatz is motivated by the consideration in Sheth & Tormen (2002). Here, we define the halo mass function in the finite-volume simulation as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dn}}{{dM}}\right)}_{\mathrm{loc}}=\displaystyle \frac{{\bar{\rho }}_{{\rm{m}}}}{M}\displaystyle \frac{d\mathrm{ln}{\sigma }_{\mathrm{loc}}^{-1}}{{dM}}\,{f}_{\mathrm{loc}}({\sigma }_{\mathrm{loc}}).\end{eqnarray} \tag{ 10 }$$

Using Equations (9) and (10), we can predict the global mass function with f(σ) = f_loc(σ) once we calibrate the functional form of f_loc as a function of σ_loc. In practice, Equation (8) provides an estimate of f_loc(σ_loc) (not f(σ_loc)). To evaluate the σ_loc–M relation, we directly compute the variance of smoothed density fields at the initial conditions of our simulations as varying smoothing scale of $R={\left(3M/4\pi {\bar{\rho }}_{{\rm{m}}}\right)}^{1/3}$. To do so, we grid N-body particles onto meshes with 512³ cells using the cloud-in-cell assignment scheme and apply the three-dimensional Fast Fourier Transform. Figure 2 shows the finite-volume effect of the mass variance measured in the phi1 and H2 runs at z = 127. We find that the finite-volume effect can be approximated as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{loc}}(M)=\sigma (M){\left(\displaystyle \frac{M}{{M}_{0}}\right)}^{\eta },\end{eqnarray} \tag{ 11 }$$

where M₀ = 2 × 10⁹ h⁻¹ M_⊙ and η = −0.02 give a reasonable fit to the phi1 run, while M₀ = 1 × 10¹⁰ h⁻¹ M_⊙ and η = −0.01 can explain σ_loc in the H2 run. For the L run, we assume no finite-volume effects on the mass variance and set σ_loc = σ.

Zoom In Zoom Out Reset image size

To estimate f_sim, we first measure the comoving number density of halos with 160 logarithmic bins in the range of M = [10⁸, 10¹⁶]. The bin size ${\rm{\Delta }}\mathrm{log}M$ is set to 0.05. We then compute the multiplicity function using Equation (8) with σ → σ_loc as in Equation (11).

For the calibration of our model with numerical simulations, we need a robust estimate of statistical errors in halo mass functions. In this paper, we adopt an analytic model of the sample variance of ${dn}/{dM}$ that was developed in Hu & Kravtsov (2003). Apart from a simple Poisson noise, the model takes into account the fluctuation of the number density of dark-matter halos in a finite volume.

Assuming that the fluctuation in ${dn}/{dM}$ is caused by underlying linear density modes at scales comparable to the simulation box size, one finds

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{Err}[{\rm{\Delta }}{n}_{\mathrm{bin}}]}{{\rm{\Delta }}{n}_{\mathrm{bin}}}={\left(\displaystyle \frac{1}{{\rm{\Delta }}{n}_{\mathrm{bin}}{V}_{\mathrm{sim}}}+S({M}_{\mathrm{bin}},z)\right)}^{1/2},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}\begin{array}{rcl}S(M,z) & = & {b}_{L}^{2}(M,z)\displaystyle \int \,\displaystyle \frac{{k}^{2}{dk}}{2{\pi }^{2}}\,{W}_{\mathrm{TH}}^{2}({{kR}}_{\mathrm{box}})\\ & & \times \,{P}_{{\rm{L}}}(k,z),\end{array}\end{eqnarray} \tag{ 13 }$$

where Err[Δn_bin] represents the statistical error in Equation (7), ${V}_{\mathrm{sim}}={L}_{\mathrm{box}}^{3}$, ${R}_{\mathrm{box}}={\left[3{V}_{\mathrm{sim}}/(4\pi )\right]}^{1/3}$, and b_L(M, z) is the linear bias at the halo mass of M and redshift z. The first term on the right-hand side in Equation (12) corresponds to the Poisson error, while the term of S(M, z) represents the sample variance. To compute Equation (13), we use the model of b_L in Tinker et al. (2010).

We validate the model of Equation (12) with the L run at z = 0. In the validation, we divide the L run into N_sub subvolumes, compute the mass function in each subvolume, and then estimate the standard deviation of the mass function over the N_sub subvolumes. We set each subvolume so that it has an equal volume. Figure 3 summarizes our validation. In this figure, we show the fractional error of the binned mass function and up-turns at higher masses indicate that the sample variance term becomes dominant. Because of Equation (8), the fractional error of f_sim is given by Equation (12).

Zoom In Zoom Out Reset image size

Before showing the results of our calibration, we perform a sanity check to see the redshift dependence of the multiplicity function in the simulation. Figure 4 shows the result of our sanity test. In this figure, we first find a best-fit model for the multiplicity function at z = 0, denoted as f_Best(σ, z = 0). We then compare the multiplicity function in the simulation at z > 0 and f_Best(σ, z = 0). If the function f(σ, z) is universal across different redshifts, then we should find a residual between the simulation results at z > 0 and f_Best(σ, z = 0), as small as the case for z = 0. In the top panel in Figure 4, the blue-dashed line represents the best-fit model f_Best(σ, z = 0), while blue circles show the simulation result at z = 0. Comparing between the blue-dashed line and blue circles, our fitting at z = 0 provides a representative model of the simulation result. The orange squares and green diamonds in the figure show the simulation results at z = 1.01 and z = 4.04, respectively. The bottom panel shows the residual between the simulation results and the model of f_Best(σ, z = 0), which highlights the prominent redshift evolution of the multiplicity function at δ_c,z /σ ≳ 3 from z = 0 to 4. In this figure, the error bars show the statistical uncertainties and we account for the sample variance caused by finite box effects in our simulation (see Section 3.3).

Zoom In Zoom Out Reset image size

This section will give the main results of this paper. Figure 10 in Appendix B summarizes the residual between the multiplicity functions in our simulations and the best-fit model at different redshifts. We find that our fitting works well across a wide range of redshifts (0 ≤ z ≤ 7). A typical difference between the simulation results and our best-fit model is an of order of 0.02 dex, except for bins at δ_c,z /σ(z) ≳ 2–3. Although the bins at δ_c,z /σ(z) ≳ 2 suffer from statistical fluctuations that are induced by the sample variance, we find that the residual is still within 0.06 dex even for such rare objects. The goodness-of-fit for our best-fit models ranges from 0.1 to 1.1. It would be worth noting that we include possible systematic errors due to modeling of subhalos in our fitting. These systematic errors can reduce the score of χ² for the best-fit model, making the best-fit χ² smaller than the number of degrees of freedom.

The redshift evolution in best-fit parameters is shown in Figure 5. The gray shaded region in each panel shows ±1σ errors of a given parameter, which are inferred from the Jacobian of Equation (14). In each panel, the circles at z ≤ 4.57 represent averages of the best-fit parameter within coarse bins of z, while those at z = 5.98 and 7.00 show the best-fit parameters. To compute the average, we set seven coarse redshift bins. The edge of coarse bins is set to [0.00, 0.30), [0.30, 0.61), [0.61, 1.01), [1.01, 1.49), [1.49, 2.44), [2.44, 3.36), and [3.36, 4.57]. We choose this binning of redshifts to ensure that the dynamical time (i.e., the ratio of the virial radius and virial circular velocity for halos) is comparable to the Hubble time at bin-centered redshifts.

Zoom In Zoom Out Reset image size

The blue-solid lines in Figure 5 present our calibrated model. We find that the following form can explain the redshift evolution of the best-fit parameters and return a similar level of the residual as in Figure 10 when used in Equation (4):

$$\begin{eqnarray}&&A(z)=0.325-0.017\,z,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\begin{array}{rcl}a(z) & = & 0.940\,{\left(1+z\right)}^{-0.02}{\left(\displaystyle \frac{1.686}{{\delta }_{c,z}}\right)}^{2}\\ & & \times \,{\left[1-0.015{\left(\displaystyle \frac{z}{1.5}\right)}^{0.01}\right]}^{-1},\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&p(z)=0.692,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&q(z)=1.611\,{\left(1+z\right)}^{0.12},\end{eqnarray} \tag{ 20 }$$

where δ_c,z is given by Equation (5).

Our model (Equations (17)–(20)) can be compared with previous models in the literature. The most popular model in Sheth & Tormen (1999) predicts a universal multiplicity function f(σ) with A = 0.322, a = 0.707, p = 0.3 and q = 1.0. It would be worth noting that the parameter of A in Sheth & Tormen (1999) is fixed by the normalization condition:

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }\,d\mathrm{ln}\sigma \,f(\sigma )=1,\end{eqnarray} \tag{ 21 }$$

where it means that all dark-matter particles reside in halos. Because our model has been calibrated by the data with a finite range of σ, it is not necessary to satisfy the condition of Equation (21). The integral of Equation (21) for our model can be well approximated as $1.37{\left(1+z/0.55\right)}^{-0.35}\,{\left(1-z/10.5\right)}^{0.21}$ at z ≤ 7 within a 0.5%-level accuracy. In reality, the upper limit in the integral of Equation (21) may be set by the free-streaming scale of dark-matter particles. If dark matter consists of Weakly Interacting Massive Particles (WIMPs), with a particle mass being ∼100 GeV, then the minimum halo mass is estimated to be 10⁻¹²–10⁻³ M_⊙ (e.g., Hofmann et al. 2001; Berezinsky et al. 2003; Green et al. 2004; Loeb & Zaldarriaga 2005; Bertschinger 2006; Profumo et al. 2006; Diamanti et al. 2015). When we set the upper limit to be σ_upp(z) =σ(M = 10⁻⁶ h⁻¹ M_⊙, z), the fraction of mass in field halos can be approximated as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{halo}} & = & {\displaystyle \int }_{-\infty }^{\mathrm{ln}{\sigma }_{\mathrm{upp}}(z)}\,d\mathrm{ln}\sigma \,f(\sigma ,z)\\ & & \simeq \,0.724{\left(1+\displaystyle \frac{z}{1.29}\right)}^{-0.28}{\left(1-\displaystyle \frac{z}{14.5}\right)}^{1.02},\end{array}\end{eqnarray} \tag{ 22 }$$

where the approximation is valid at z ≤ 7 within a 0.4%-level accuracy. Equation (22) is less sensitive to the choice of the minimum halo mass as long as we vary the minimum halo mass in the range of 10⁻¹²–10⁻³ h⁻¹ M_⊙. Our model predicts that about 72% of the mass density in the present-day universe resides in dark-matter halos.

For parameter a, our model (Equation (18)) shows a modest redshift evolution with a level of ∼20% from z = 7 to z = 0. Note that an effective critical density $\sqrt{a(z)}\,{\delta }_{c,z}$ becomes less dependent on z and only evolves by 2%–3% in the range of 0 ≤ z ≤ 7. The redshift dependence of a is mostly consistent with the model in Bhattacharya et al. (2011), but the overall amplitude differs by ∼20%. Note that the model in Bhattacharya et al. (2011) has been calibrated for the halo mass function when the mass is defined by the FoF algorithm. Because the FoF mass is expected to strongly depend on inner density profiles and substructures (More et al. 2011), our model does not need to match the one in Bhattacharya et al. (2011). The present-day values of a, p, and q are in good agreement with the result in Comparat et al. (2017), which presented the calibration of the halo mass function at z = 0 with the MultiDark simulation (Prada et al. 2012; Klypin et al. 2016). Note that Comparat et al. (2017) adopted same halo finder, mass definition and cosmology as ours.

Despali et al. (2016) argued that the multiplicity function can be expressed as a universal function once one includes the redshift dependence of the spherical critical density δ_c . Although we also adopt the redshift-dependent critical density, we find that the parameters in the multiplicity function depend on redshifts (also see Figure 4). The main difference between the analysis in Despali et al. (2016) and ours is halo identifications in N-body simulations. We define halos by the FoF algorithm in six-dimensional phase-space, while Despali et al. (2016) adopted a spherical-overdensity algorithm with a smoothing scale being the distance to the tenth nearest neighbor for a given N-body particle. Note that halo finders based on particle positions may not distinguish two merging halos. The ROCKSTAR algorithm utilizes the velocity information of N-body particles. This allows us to very efficiently determine particle-halo memberships, even in major mergers.

We next compare our model of the halo mass functions with previous models in the literature. For comparison, we consider six representative models (which are summarized in Table 2). Three of the models assume a universal functional form of the multiplicity function, while the other models include some redshift evolution. Among the previous studies, Despali et al. (2016) investigated the mass function when using the virial halo mass, and their results can be directly compared to ours. Tinker et al. (2008) and Watson et al. (2013) studied mass functions for various spherical-overdensity masses. Here we use their fitting formula, which takes into account the dependence of spherical-overdensity parameters, while the calibration for the virial halo mass has not been done in Tinker et al. (2008) and Watson et al. (2013). Bhattacharya et al. (2011) have calibrated the halo mass function when the mass is defined by the FoF algorithm, which indicates that a direct comparison with our results may not be appropriate (e.g., see More et al. 2011, for details).

The left-hand panel in Figure 6 shows the comparison of fitting formulas for the halo mass function at z = 0, while the right-hand panel presents the case at z = 6. Our model is shown by blue-solid lines in both panels and the gray shaded region in the figure represents the range of halo masses not explored by our simulations. At z = 0, our model is in good agreement with most of previous models in the range of 10⁹⁻¹³ h⁻¹ M_⊙, but there are 15%-level discrepancies at the high-mass end (M ≳ 10¹⁴ h⁻¹ M_⊙). This figure implies that previous fitting formulas are not sufficient to predict cosmology-dependence of the mass function for cluster-sized halos at z = 0, even within a concordance ΛCDM cosmology. A commonly adopted model by Tinker et al. (2008) may have a systematic error to predict the cluster abundance with a level of 20%–30%, but the exact amount of systematic errors should depend on the definition of spherical-overdensity halo masses. We need larger N-body simulations (e.g., see Ishiyama et al. 2021, for a relevant example) to precisely calibrate the mass function at high-mass ends and make a robust conclusion about systematic errors in cluster cosmology (e.g., McClintock et al. 2019; Bocquet et al. 2020; Klypin et al. 2021). We leave further investigation of cluster mass functions for future studies.

Zoom In Zoom Out Reset image size

At z = 6, our model shows an offset from the universal models by Sheth & Tormen (1999) and Despali et al. (2016). The difference reaches a 15%–20% level at M = 10⁹⁻¹⁰ h⁻¹ M_⊙ and becomes larger at higher masses. This is mainly caused by the redshift evolution of the amplitude in the multiplicity function, A(z). Our model predicts that the amplitude in the multiplicity function decreases at higher redshifts, and this trend is clearly found in the simulation results (see Figure 4).

An important implication of our calibration is that modeling of the halo mass function for ΛCDM cosmologies can affect constraints of nature of dark-matter particles by high-redshift galaxy number counts (e.g., Pacucci et al. 2013; Schultz et al. 2014; Menci et al. 2016; Corasaniti et al. 2017). Warm dark matter (WDM) is an alternative candidate of cosmic dark matter with free streaming due to its thermal motion. Some physically motivated extensions of the Standard model predict the existence of WDM with masses in keV range, such as sterile neutrino (e.g., Adhikari et al. 2017; Boyarsky et al. 2019). Structure formation with WDM particles can be suppressed at scales below free-streaming lengths λ_fs, while the bottom-up formation of dark-matter halos remains as in the standard CDM paradigm at scales larger than λ_fs. A characteristic mass scale for λ_fs has been estimated as ∼10¹⁰ M_⊙ for WDM with a particle mass of O(1) keV (e.g., Bode et al. 2001).

In a hierarchical structure formation, less massive galaxies form at higher redshifts. At a fixed cosmic mean density, WDM particles with larger masses are less effective at suppressing the growth of low-mass halos (e.g., Schneider et al. 2012). Assuming that high-z galaxies only form in collapsed halos, the observed abundance of high-z galaxies can thus provide a lower limit to the particle mass of WDM.

Recent high-resolution N-body simulations for WDM have indicated that the correction of the halo abundance due to free streaming can be expressed as a universal form of:

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dM}}(M,z){\left|{}_{\mathrm{WDM}}={ \mathcal R }(M)\displaystyle \frac{{dn}}{{dM}}(M,z)\right|}_{\mathrm{CDM}}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{ \mathcal R }(M)={\left[1+{\left(\displaystyle \frac{\alpha {M}_{\mathrm{hm}}}{M}\right)}^{\beta }\right]}^{\gamma },\end{eqnarray} \tag{ 24 }$$

where ${dn}/{dM}{| }_{\mathrm{WDM}}$ is the halo mass function for WDM cosmologies, ${dn}/{dM}{| }_{\mathrm{CDM}}$ is the counterpart of CDM, Lovell (2020) found that α = 2.3, β = 0.8 and γ = −1.0 provide a reasonable fit to simulation results. In Equation (24), M_hm is so-called "half-mode" mass, which is defined as the mass scale that corresponds to the power spectrum wavenumber at which the square root of the ratio of the WDM and CDM power spectra is 0.5 (Schneider et al. 2012). The mass M_hm depends on the particle mass of WDM m_WDM:

$$\begin{eqnarray}&&{M}_{\mathrm{hm}}=\displaystyle \frac{4\pi }{3}\,{\bar{\rho }}_{m}{\left[{s}_{\mathrm{WDM}}{\left({2}^{\mu /5}-1\right)}^{-\tfrac{1}{2\mu }}\right]}^{3},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}\begin{array}{rcl}{s}_{\mathrm{WDM}} & = & 0.153\,[{h}^{-1}\,\mathrm{Mpc}]{\left(\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{WDM}}}{0.25}\right)}^{0.11}\\ & & \times \,{\left(\displaystyle \frac{{m}_{\mathrm{WDM}}}{1\,\mathrm{keV}}\right)}^{-1.11}{\left(\displaystyle \frac{h}{0.7}\right)}^{1.22},\end{array}\end{eqnarray} \tag{ 26 }$$

where μ = 1.12 and Ω_WDM is the dimensionless density parameter of WDM. According to Equation (23), an accurate calibration of mass function for ΛCDM cosmologies is essential to predict the counterpart of WDM. We caution here that Equations (23) and (24) have been validated at z = 0 and z = 2 in Lovell (2020). We assume that these equations are valid at z ∼ 6. We leave a validation of Equations (23) and (24) at z ≳ 2 for future studies.

Figure 7 demonstrates the importance of the calibration of CDM halo mass function when one constrains WDM masses with high-redshift galaxy number counts. For given limiting magnitude and redshift, the cumulative galaxy number density should be smaller than the whole halo mass function within WDM cosmologies. This leads to

$$\begin{eqnarray}&&{\phi }_{\mathrm{obs}}\equiv {\int }_{{L}_{\mathrm{cut}}}^{\infty }{dL}\,\displaystyle \frac{{{dn}}_{\mathrm{gal}}}{{dL}}\leqslant {\int }_{{M}_{\min }}^{{M}_{\max }}\,{dM}\,\displaystyle \frac{{dn}}{{dM}}{\left.\right|}_{\mathrm{WDM}},\end{eqnarray} \tag{ 27 }$$

where L_cut is the luminosity corresponding to the limiting magnitude, ${{dn}}_{\mathrm{gal}}/{dL}$ is the galaxy luminosity function, and we set ${M}_{\min }=1\,{h}^{-1}{M}_{\odot }$ and ${M}_{\max }={10}^{16}\,{h}^{-1}{M}_{\odot }$. Note that our lowest mass ${M}_{\min }=1\,{h}^{-1}{M}_{\odot }$ is much smaller than the half-mode mass M_hm ≳ 10¹⁰ h⁻¹ M_⊙ for m_WDM ≳ 1 keV. For the UV luminosity function at z = 6 in the Hubble Frontier Fields with the limiting AB magnitude of −12.5 (Livermore et al. 2017), the lower bound of ϕ_obs has been estimated as $\mathrm{log}({\phi }_{\mathrm{obs}}\,[{\mathrm{Mpc}}^{-3}])\gt 0.01$ at a 2σ confidence level (Menci et al. 2016). Assuming the best-fit cosmological parameters in Planck Collaboration et al. (2016) and WDM is made of the whole abundance of dark matter, our model of the halo mass function with Equation (24) allows us to reject WDM with a mass smaller than 2.71 keV at the 2σ level. This lower limit is degraded to 2.27 keV and 1.96 keV for the commonly adopted models by Sheth & Tormen (1999) and Press & Schechter (1974), respectively. This simple example highlights that calibration of the mass function for ΛCDM cosmologies is essential to accurate predictions of the mass function for WDM cosmologies.

Zoom In Zoom Out Reset image size

Our lower limit of the WDM particle mass can be compared to other methods. For example, Baur et al. (2016) found a lower limit of 2.96 keV from observations of the Lyα forest, while Palanque-Delabrouille et al. (2020) placed a lower limit of 5.3 keV (a similar limit is found by Iršič et al. (2017)). Chatterjee et al. (2019) used high-redshift 21 cm data from EDGES to rule out WDM with m_WDM < 3 keV. Note that our limit is less sensitive to details of baryonic physics than the others because our analysis relies on the cumulative abundance of dark-matter halos.

For a conservative analysis, we consider that a significant small halos of M = 1 h⁻¹ M_⊙ can be responsible to the observed galaxy abundance at high redshifts. To further tighten the limit of WDM particle mass, it would be interesting to discuss more realistic halo mass scales to the faintest galaxy at z ∼ 6. In the Planck ΛCDM cosmology, we find that the minimum halo mass of ∼10⁷ h⁻¹ M_⊙ provides the cumulative halo abundance of 1.35 Mpc⁻³, which is close to the observed galaxy abundance at z = 6. When setting the minimum halo mass to 10⁷ h⁻¹ M_⊙ in Equation (27), we find a stringent 2σ limit of m_WDM > 14.1 keV. However, this stringent limit is very sensitive to the choice of the minimum halo mass. For the minimum halo mass of 10⁶ h⁻¹ M_⊙, the limit changes to m_WDM > 6.23 keV. This simple analysis implies that a more detailed modeling of galaxy-halo connections at M_vir ∼ 10⁶⁻⁷ h⁻¹ M_⊙ and z ∼ 6 would be worth pursuing in future work.

Our model of halo mass functions is calibrated against N-body simulations in the ΛCDM cosmology consistent with Planck16. In terms of studies of large-scale structure, Ω_m0 and σ₈ are the primary parameters and the simulations in this paper adopt Ω_m0 = 0.31 and σ₈ = 0.83. Therefore, our calibration of Equations (17)–(20) may be subject to an overfitting to the specific cosmological model. To examine the dependence of our model on cosmological models, we use another halo catalog from N-body simulations with a different ΛCDM model. For this purpose, we use the Bolshoi simulation in Klypin et al. (2011) and the first MultiDark-Run1 simulation performed in Prada et al. (2012). The Bolshoi simulation consists of 2048³ particles in a volume of ${250}^{3}\,{\left[{h}^{-1}\,\mathrm{Mpc}\right]}^{3}$ and assumes the cosmological parameters of Ω_m0 = 0.27, Ω_b0 = 0.0469, Ω_Λ = 1 −Ω_m0 = 0.73, h = 0.70, n_s = 0.95, and σ₈ = 0.82. These are consistent with the five-year observation of the cosmic microwave background obtained by the WMAP satellite (Komatsu et al. 2009) and we refer to them as the WMAP5 cosmology. The MultiDark-Run1 simulation adopted the same cosmological model and the number of particles as in the Bolshoi simulation, while the volume is set to $1\,{\left[{h}^{-1}\,\mathrm{Gpc}\right]}^{3}$. We use the ROCKSTAR halo catalog at z = 0, 0.53, 1.00, 1.96, 2.93 and 4.07 from the Bolshoi and MultiDark-Run1 simulations. ¹³ To compute our model prediction for the WMAP5 cosmology, we fix the functional form of Equation (4) and parameters in Equations (17)–(20) but include the cosmology-dependence of the critical density δ_c,z and the top-hat mass variance σ, accordingly.

Figure 8 summarizes the multiplicity function in the Bolshoi and MultiDark-Run1 simulations. In this figure, the dashed lines show the predictions by our model for the WMAP5 cosmology, while the different symbols represent the simulation results. We find that our model can reproduce simulation results within a 10% level even for the WMAP5 cosmology at 0.7 ≲ δ_c,z /σ ≲ 2.5. At high-mass ends (δ_c,z /σ ≳ 2.5), our model tends to underestimate the halo abundance by ∼30%–40% in a wide range of redshifts. It is worth noting that the residual between our model and the WMAP5-based simulation is less dependent on redshifts. In fact, we found a better matching between our models and the WMAP5-based simulations when reducing the overall amplitude in the parameter a(z) by 4%–5%. In summary, our model cannot predict the simulation results for the WMAP5 cosmology with the same level as in the Planck cosmology. The 10%-level difference in Ω_m0 can cause systematic uncertainties in our model predictions with a level of 10%, except for high-mass ends. Future studies would need to calibrate the cosmological dependence of a(z) for precision cosmology based on galaxy clusters. At low masses and high redshifts, our model can provide a reasonable fit to the simulations adopting the WMAP5 cosmology. According to this fact, we examine how much the WDM limit in Section 5.4 is affected by the choice of underlying cosmology. For the WMAP5 cosmology, we find a 2σ limit of m_WDM > 2.75 keV, which differs from our fiducial limit by only ∼1.4%.

Zoom In Zoom Out Reset image size

Our calibration of halo mass functions relies on dark-matter-only N-body simulations and ignores possible baryonic effects. Baryonic effects on halo mass functions have been studied with a set of hydrodynamical simulations (e.g., Stanek et al. 2009; Cui et al. 2012, 2014; Sawala et al. 2013; Cusworth et al. 2014; Martizzi et al. 2014; Bocquet et al. 2016; Beltz-Mohrmann & Berlind 2021). The evolution of cosmic baryons is governed not only by gravity but also complex processes associated with galaxy formation. Relevant processes include gas cooling, star formation and energy feedback from supernovae (SN) and Active Galactic Nuclei (AGNs). Adiabatic gas heated only by gravitational processes affects the halo mass function with a level of ≲2%–3% (Cui et al. 2012), while radiative cooling, star formation and SN feedback increase individual halo masses due to condensation of baryonic mass at the halo center, which changes the halo mass function (Stanek et al. 2009; Cui et al. 2012). Efficient SN feedback can decrease the halo mass function at M ≲ 10¹¹ M_⊙ by ∼20%–30% (Sawala et al. 2013). The mass function at M ≃ 10¹³⁻¹⁴ M_⊙ would be affected by AGN feedback, while current hydrodynamical simulations adopt a sub-grid model to include the AGN feedback. Because a variety of sub-grid models have been proposed, the impact of the AGN feedback on the mass function is still uncertain (Cui et al. 2014; Cusworth et al. 2014; Martizzi et al. 2014; Bocquet et al. 2016).

To account for the baryonic effects on the halo mass function, it is important to correct individual halo masses according to baryonic processes. Baryons do not change the abundance of dark-matter halos, but they do affect the internal structures and spherical masses of halos. Hence, one may be able to model the halo mass function in the presence of baryons by abundance matching between the gravity-only and hydrodynamical simulations for a given definition of the halo mass (e.g., Beltz-Mohrmann & Berlind 2021). In this sense, our calibration of the halo mass function provides a baseline model and still plays an important role in understanding the baryonic effects on large-scale structures.