17.4. The halo mass function¶
We end this chapter with a brief discussion of the overall mass distribution of dark-matter halos, but we will revisit this topic in more detail in Chapter 19.2. In this chapter, we have learned that dark-matter overdensities grow in a scale-independent manner during the epochs of matter- and dark-energy domination at \(z\lesssim 3400\) during which galaxy formation occurs. Initial small overdensities grow until they gravitationally collapse into a dark matter halo. In the spherical-collapse approximation of Section 17.2, this collapse finishes when the linearly-extrapolated overdensity reaches the critical value of \(\delta_c \approx 1.686\). At that point, a virialized halo with the properties that we discussed in Section 17.3.1 has formed. We have also discussed the overall statistics of density fluctuations at different epochs in Section 17.1.5. In the linear approximation, the fluctuations are given by an isotropic Gaussian random field with a power spectrum given, at \(z=0\), by Equation (17.79). In this final section, we will see how by combining the statistics of density fluctuations with the simple spherical-collapse approximation for the formation of halos, we can derive the halo mass function—the number of dark-matter halos of a given mass—at different epochs.
Because dark-matter halos are extended regions that are overdense with respect to the mean background density, we commence by smoothing the density field \(\delta(\vec{x},t)\) using a window function \(W(\vec{x};R)\) where the radius \(R\) is related to the mass \(M\) of the dark-matter halo that forms; we will always assume below that this window function is spherically symmetric. To determine the relation between \(R\) and \(M\), recall that the “\(\vec{x}\)” in the density field \(\delta(\vec{x},t)\) is the comoving position, that is, the position with the overall effect of the expansion of the Universe removed. The radius \(R\) is therefore also a comoving distance. Thus, we can compute it at a very early time \(t_i\), when \(\delta(\vec{x},t_i) \ll 1\) everywhere and the density \(\rho(\vec{x},t_i)\) was indistinguishable from the mean matter density \(\bar{\rho}_m(t_i)\). Considering then, for example, a spherical top-hat window function \(W(\vec{x};R)\), we have that the enclosed mass within the comoving radius \(R\) is \(M = 4\pi \bar{\rho}_{0,m} R^3/3 = 4\pi \Omega_{0,m}\rho_c R^3/3\), where \(\rho_c\) is the critical density. In this expression, the mean matter and critical density are evaluated today, because the calculation is performed in comoving coordinates. Using Equation (17.3) for the critical density, we have that \begin{align}\label{eq-dmhaloform-PS-radius-for-mass} R & = \left({2\,GM \over \Omega_{0,m}\,H_0^2}\right)^{1/3} = 1.82\,\mathrm{Mpc}\,\left({M \over 10^{12}\,M_\odot}\right)^{1/3}\,. \end{align} Thus, the matter that forms a dark-matter halo like that of the Milky Way, with a mass of \(\approx 10^{12}\,M_\odot\), comes from a region in the very early Universe within a comoving radius of roughly 2 Mpc. This is about an order of magnitude larger than the virial radius of the dark-matter halo that forms, reflecting the higher average density of the collapsed dark-matter halo when compared to the mean matter density. Through the mass-to-radius mapping of Equation (17.114), we can express the smoothing scale equivalently as a radius \(R\) or a mass \(M\). Below we will use both \(R\) and \(M\) depending on the context.
We then convolve the linear density field \(\delta(\vec{x},t)\) with the window function \(W(\vec{x};R)\) to calculate the smoothed field \(\delta_S(\vec{x},t;M)\) as \begin{equation} \delta_S(\vec{x},t;M) = \int \mathrm{d}\vec{x}'\,\delta(\vec{x'},t)\,W(\vec{x}-\vec{x}';R)\,. \end{equation} Because the linear density field \(\delta(\vec{x},t)\) follows an isotropic Gaussian probability distribution with mean zero, so does the smoothed density field \(\delta_S(\vec{x},t;M)\) and its statistics is therefore fully determined by its variance \(\sigma^2(M,t)\), which is given by \begin{equation}\label{eq-dmhaloform-sigmaM} \sigma^2(M,t) = {1\over 2\pi^2}\int_0^\infty\mathrm{d}k\,k^2\,P_m(k,a)\,W^2_k(R)\,, \end{equation} where \(P_m(k,a)\) is the linear matter power spectrum from Equation (17.79) evaluated at the scale factor \(a\) that corresponds to the time \(t\); as usual, we use \(t\), \(a\), and \(z\) interchangeably to denote the time at which cosmological quantities are evaluated and we therefore also denote this variance as \(\sigma^2(M,a)\) or \(\sigma^2(M,z)\). The function \(W_k(R)\) is the Fourier transform of the window function \(W(\vec{x};R)\), which for a spherical top-hat is given by \begin{equation} W_k(R) = 3\,{j_1(kR) \over kR} = 3\,{\sin kR - kR\,\cos kR \over (kR)^3}\,, \end{equation} where \(j_1(x)\) is a spherical Bessel function (see Appendix B.3.5). Because of the complexity of the transfer function, \(\sigma^2(M,a)\) needs to be computed numerically, which is tricky due to the oscillatory nature of the Bessel function. However, because of the scale-independent growth of density perturbations, we can compute \(\sigma^2(M,a=1)\) and then scale using the relative growth factor \(D^2_+(a)/D_+^2(a=1)\) to obtain \(\sigma^2(M,a)\) at another scale factor.
In the spherical-collapse picture, dark-matter halos with mass \(M\) will have formed abundantly at redshift \(z\) if \(\sigma(M,z) \gtrsim \delta_c\), because overdensities in regions that encompass a mass \(M\) will at that redshift typically have grown to have \(\delta > \delta_c\). Computing \(\sigma(M,z)\) is straightforward using Equation (17.116) and the linear matter power spectrum from Equation (17.79).
Below, we implement the calculation of \(\sigma(M,z)\) using the BBKS transfer function from Equation (17.74) with the Sugiyama (1995) prescription for the variable \(q\) (Equation 17.76). To be able to efficiently compute \(\sigma(M,z)\) at different redshifts, which only differ by an overall relative growth factor \(D_+(a)/D_+(a=1)\), we factor out the dependence of the matter power spectrum on the growth factor, such that we directly compute \(\sigma(M,z=0)\) and we can obtain \(\sigma(M,z)\) as \(\sigma(M,z) = \sigma(M,z=0)D_+(a)/D_+(a=1)\). Because an important cosmological parameter that characterizes the amplitude of density fluctuations is the value of \(\sigma(M,z=0)\) evaluated at a comoving radius of \(R = 8\,h^{-1}\,\mathrm{Mpc}\)—this parameter is known as \(\sigma_8\)—we first use the code to compute \(\sigma_8\) for our adopted cosmological parameters.
[11]:
# Constants
Om0= 0.3111
Ode0= 0.6889
H0= 67.70*u.km/u.s/u.Mpc
h= (H0/(100.*u.km/u.s/u.Mpc)).to_value(u.dimensionless_unscaled)
Ob0= 0.0490
ns= 0.9665
As= numpy.exp(3.047)*1e-10
kp= 0.05/u.Mpc
from scipy import integrate
from astropy.constants import c
# Transfer function approximation
T0CMB= 2.7255*u.K
def transfer_bbks(k,model='BBKS',Om0=Om0,h=h,Ob0=Ob0,T0CMB=T0CMB):
if model.lower() == 'bbks':
q= (k/(Om0*h**2./u.Mpc)).to_value(u.dimensionless_unscaled)
elif model.lower() == 'sugiyama':
q= (k*(T0CMB/2.7/u.K)**2./(Om0*h**2./u.Mpc)/numpy.exp(-Ob0-numpy.sqrt(2.*h)*Ob0/Om0))\
.to_value(u.dimensionless_unscaled)
return numpy.log(1.+2.34*q)/2.34/q*(1.+3.89*q+(16.1*q)**2.
+(5.46*q)**3.+(6.71*q)**4.)**-0.25# Growth factor stuff
def H(a):
# Assumes flat Universe
return H0*numpy.sqrt(Om0/a**3.+Ode0)
def growth_factor(a):
return (5.*Om0/2.*H(a)/H0\
*integrate.quad(lambda ap: 1./ap**3./(H(ap)/H0)**3.,
0.,a)[0]).to_value(u.dimensionless_unscaled)
# Matter power spectrum without the growth factor, so we can apply it later
def matter_power_spectrum_nogrowth(k,model='BBKS'):
return 8.*numpy.pi**2./25.*c**4.*As*kp/Om0**2./H0**4.\
*transfer_bbks(k,model=model)**2.*(k/kp)**ns
def window(x):
# Top-hat window function in Fourier space
return 3.*(numpy.sin(x)-x*numpy.cos(x))/x**3.
def sigma2_R(R,a=1.):
# sigma2 computed as a function of R at scale factor a
return integrate.quad(lambda k: k**2.\
*matter_power_spectrum_nogrowth(k/u.Mpc,
model='Sugiyama').to_value(u.Mpc**3)\
*window((k*R/u.Mpc).to_value(u.dimensionless_unscaled))**2,
0,numpy.inf)[0]/2./numpy.pi**2.*growth_factor(a)**2.
def sigma2(M,a=1.):
# sigma2 computed as a function of mass M
R= 1.82*u.Mpc*(M/1e12/u.Msun)**(1./3.)
return sigma2_R(R,a=a)
def rescale_sigma(sig,a=1.,aref=1.):
# rescale sigma to another scale factor
return sig*growth_factor(a=a)/growth_factor(a=aref)
print(f"The value of sigma_8 is {numpy.sqrt(sigma2_R(8./h*u.Mpc)):.3f}")
The value of sigma_8 is 0.846
This value is similar to the exact value computed using the exact, numerically-calculated transfer function reported by Planck Collaboration et al. (2020b): they find \(\sigma_8 = 0.8102\).
Next, we can determine how \(\sigma(M,z)\) depends on mass at different redshifts and which masses are above the critical density threshold is shown in Figure 17.8.
[13]:
ms= numpy.geomspace(1e5,1e16,21)*u.Msun
sigma_z0= numpy.array([numpy.sqrt(sigma2(m,a=1.)) for m in ms])
figure(figsize=(8,5))
semilogx(ms,sigma_z0,label=r'$z=0$')
semilogx(ms,rescale_sigma(sigma_z0,a=1./(1.+2.)),label=r'$z=2$')
semilogx(ms,rescale_sigma(sigma_z0,a=1./(1.+5.)),label=r'$z=5$')
xlabel(r'$M/M_\odot$')
ylabel(r'$\sigma$')
axhline(1.686,color='k',lw=1.5)
text(5e14,2,r'$\sigma=\delta_c$',fontsize=18.)
legend(frameon=False,fontsize=18.);
Figure 17.8: The variance of the smoothed linear matter density field.
Looking at the \(z=0\) curve, we see that \(\sigma(M,z)\) monotonically decreases with increasing mass. Low-mass density fluctuations are therefore always typically larger than higher-mass density fluctuations. This means that low-mass dark-matter halos will generally form at earlier times than high-mass dark-matter halos, because they have a much higher chance of exceeding the critical density threshold. This implies that structure formation in the standard cold-dark-matter cosmological model is hierarchical where structure forms from the bottom up. This did not have to be the case: if dark matter is hot or very warm in the early Universe, meaning that it moves at relativistic or close-to relativistic speeds, free-streaming of dark matter away from overdensities means that small-scale overdensities cannot directly form through gravitational collapse in the way that high-mass dark-matter halos can. In such models, \(\sigma(M,z)\) decreases at low masses and the only way to create low-mass dark-matter halos is through the fragmentation of high-mass dark-matter halos. Such a structure formation scenarios is known as a top-down model.
The mass \(M^*(z)\) at which \(\sigma(M^*,z) = \delta_c\) is the mass above which dark-matter halos are rare at the given redshift, because density fluctuations at higher masses \(M \gtrsim M^*(z)\) generally do not exceed \(\delta_c\). We see from the figure above that this mass is \(M^* \approx 10^{13}\,M_\odot\) at \(z=0\). Because \(\sigma(M,z)\) decreases by an overall factor at higher redshifts, \(M^*(z)\) decreases and we see, for example, that it is only \(M^* \approx 10^{10}\,M_\odot\) at \(z=2\) and \(M^* \approx 10^{6}\,M_\odot\) at \(z=5\). Because galaxies typically inhabit dark-matter halos with masses \(\gtrsim 10^8\,M_\odot\), all galaxies are therefore quite rare at \(z=5\) compared to their abundance today and massive galaxies are rare at \(z=2\). Groups and clusters of galaxies with \(M \gtrsim 10^{13}\,M_\odot\) are rare still today.
Now that we have determined the statistics of overdensities corresponding to dark-matter halos with mass \(M\), all that remains to be done is using this statistics to determine the halo mass function. A simple physical picture, albeit with some ad-hoc ingredients, was proposed by Press & Schechter (1974). Because, in the spherical-collapse picture, a dark-matter halo with a mass \(M\) has formed at time \(t\) at position \(\vec{x}\) when \(\delta_S(\vec{x},t) > \delta_c \approx 1.686\), Press & Schechter (1974) posited that the mass fraction \(f(>M,t)\) of dark-matter halos with masses larger than \(M\) is twice the probability that \(\delta_S(\vec{x},t) > \delta_c\) (we will explain the factor of two shortly). Because \(\delta_S(\vec{x},t)\) follows a Gaussian probability distribution function with mean zero and variance \(\sigma^2(M,t)\) given by Equation (17.116), the probability that \(\delta_S(\vec{x},t) > \delta_c\) is given by (see Equation B.19) \begin{equation}\label{eq-dmhaloform-hmf-pdeltas} P[\delta_S(\vec{x},t) > \delta_c] = {1 \over 2}\,\mathrm{erfc}\left[{\delta_c \over \sqrt{2}\,\sigma(M,t)}\right]\,. \end{equation} This expression demonstrates the need for the factor of two in converting this to the mass fraction of halos with masses larger than \(M\). Because \(\sigma(M,t)\) increases without bound as \(M \rightarrow 0\) in cold-dark-matter models, \(P[\delta_S(\vec{x},t) > \delta_c] \rightarrow 1/2\) when \(M \rightarrow 0\) and without the factor of two, only half of the mass of the Universe would be part of some bound object. This makes sense in the linear-theory treatment: only regions that are initially denser than the mean matter density can ever collapse, so the half of the mass that is in initially underdense regions never reaches \(\delta = \delta_c\). However, mass that is initially part of an underdense region can be accreted into an overdense region, essentially if it is part of a larger region that is itself overdense. Press & Schechter (1974) assumed that the result of this accretion process is simply that the abundance of halos with mass \(M\) is twice that predicted by \(P[\delta_S(\vec{x},t) > \delta_c]\) in Equation (17.118) and a more rigorous derivation of the halo mass function in the spherical-collapse picture that we will discuss in Chapter 19.2 shows that this assumption is in fact correct. But for now we will simply assume the factor of two and move on. The discussion in Chapter 19.2 will also make it clear that in detail the Press-Schechter halo mass function is for main halos only, that is, it does not include subhalos orbiting within a larger halo.
Given Equation (17.118), the cumulative mass per unit volume in halos with masses larger than \(M\) is then given by \begin{equation} \mathcal{M}(>M,t) = \bar{\rho}_{0,m}\,f(>M,t) = \bar{\rho}_{0,m}\,\mathrm{erfc}\left[{\delta_c \over \sqrt{2}\,\sigma(M,t)}\right]\,, \end{equation} where we denote the cumulative mass as \(\mathcal{M}\) and this is the cumulative mass function in a comoving volume (hence, we use \(\bar{\rho}_{0,m}\) rather than \(\bar{\rho}_{m}\)(t)). The halo mass function in a comoving volume itself is the derivative of this cumulative distribution with respect to \(M\) and then divided by \(M\) to convert to a number density of halos; thus, it is given by \begin{equation}\label{eq-dmhaloform-hmf-ps} {\mathrm{d} n \over \mathrm{d} \ln M} = \sqrt{{2\over \pi}}\,{\bar{\rho}_{0,m} \over M}\,{\delta_c \over \sigma}\,\exp\left(-{\delta_c^2 \over 2\sigma^2}\right)\,{\mathrm{d} \ln \sigma^{-1} \over \mathrm{d} \ln M}\,, \end{equation} where we have abbreviated \(\sigma \equiv \sigma(M,t)\) and, as is usual, we have expressed the halo mass functions in logarithmic mass bins. This is the Press-Schechter halo mass function. Before implementing it, we note that we can write this in the form \begin{equation}\label{eq-dmhaloform-hmf-ps-alt} {\mathrm{d} n \over \mathrm{d} \ln M} = {\bar{\rho}_{0,m} \over M}\,f_{\mathrm{PS}}(\nu)\,{\mathrm{d} \ln \sigma^{-1} \over \mathrm{d} \ln M}\,, \end{equation} where the peak height \(\nu = \delta_c/\sigma(M,t)\) and \begin{equation}\label{eq-dmhaloform-hmf-fps} f_{\mathrm{PS}}(\nu) = \sqrt{{2\over \pi}}\,\nu\,e^{-\nu^2/2}\,, \end{equation} is the multiplicity function. To compute the halo mass function we need \(\sigma(M,t)\) from Equation (17.116) and its logarithmic derivative \(\mathrm{d} \ln \sigma^{-1} / \mathrm{d} \ln M = -\mathrm{d}\ln \sigma / \mathrm{d} \ln M\) with respect to \(M\). This derivative is straightforward to compute numerically using a similar integral as for \(\sigma(M,t)\), but an important point to note is that the derivative does not depend on time, because \(\sigma(M,t)\) only changes by an overall factor at different times.
The following functions implement \(\mathrm{d} \ln \sigma / \mathrm{d} \ln M\):
[17]:
def dwindowdx(x):
return 3.*((3.*x*numpy.cos(x)+(x**2.-3.)*numpy.sin(x))/x**4.)
def dlnsigmadlnM(M):
# constant with a, so set a=1.
R= 1.82*u.Mpc*(M/1e12/u.Msun)**(1./3.)
return integrate.quad(lambda k: k**2.\
*(k*R/u.Mpc).to_value(u.dimensionless_unscaled)\
*matter_power_spectrum_nogrowth(k/u.Mpc,
model='Sugiyama').to_value(u.Mpc**3)\
*window((k*R/u.Mpc).to_value(u.dimensionless_unscaled))\
*dwindowdx((k*R/u.Mpc).to_value(u.dimensionless_unscaled)),
0,numpy.inf)[0]\
/6./numpy.pi**2.*growth_factor(1.)**2./sigma2_R(R,a=1.)
Then it is straightforward to implement the halo mass function using the form from Equation (17.121):
[18]:
from astropy.constants import G
rhom= Om0*3.*H0**2./8./numpy.pi/G
def halo_mass_function(M,a=1.,dlnsigdlnM=None,sig=None):
# d n / d ln M
# Use sigma and d ln sigma / d ln M if provided, otherwise compute
if sig is None:
sig= numpy.sqrt(sigma2(M,a=a))
if dlnsigdlnM is None:
dlnsigdlnM= dlnsigmadlnM(M)
# f function
def fnu(nu):
return numpy.sqrt(2./numpy.pi)*nu*numpy.exp(-nu**2./2.)
nu= 1.686/sig
return -rhom/M*fnu(nu)*dlnsigdlnM
We can then compute the halo mass function at different redshifts:
[19]:
ms= numpy.geomspace(1e5,1e16,51)*u.Msun
# Only compute sigma and dln sigma / d ln M once,
# rescale sigma to other redshifts as necessary
sig_z0= numpy.sqrt(numpy.array([sigma2(m,a=1.) for m in ms]))
dlnsigdlnM= numpy.array([dlnsigmadlnM(m) for m in ms])
figure(figsize=(8,5))
loglog(ms,
u.Quantity(halo_mass_function(ms,a=1.,sig=sig_z0,
dlnsigdlnM=dlnsigdlnM)).to_value(1./u.Mpc**3),
label=r'$z=0$')
loglog(ms,
u.Quantity(halo_mass_function(ms,a=1./(1.+2.),
sig=rescale_sigma(sig_z0,a=1./(1.+2.)),
dlnsigdlnM=dlnsigdlnM)).to_value(1./u.Mpc**3),
label=r'$z=2$')
loglog(ms,
u.Quantity(halo_mass_function(ms,a=1./(1.+5.),
sig=rescale_sigma(sig_z0,a=1./(1.+5.)),
dlnsigdlnM=dlnsigdlnM)).to_value(1./u.Mpc**3),
label=r'$z=5$')
loglog(ms,
u.Quantity(halo_mass_function(ms,a=1./(1.+10.),
sig=rescale_sigma(sig_z0,a=1./(1.+10.)),
dlnsigdlnM=dlnsigdlnM)).to_value(1./u.Mpc**3),
label=r'$z=10$')
xlabel(r'$M/M_\odot$')
ylabel(r'$\mathrm{d} n/\mathrm{d} \ln M\,(\mathrm{Mpc}^{-3})$')
ylim(1e-10,4e4)
legend(frameon=False,fontsize=18.);
Figure 17.9: The Press-Schechter halo mass function.
At low masses \(M \lesssim M^*\), the halo mass function is approximately a power-law \(\propto M^{-1}\), but above \(M \gtrsim M^*\), this power-law is exponentially cut-off by the \(e^{-\nu^2/2}\) factor in Equation (17.122). The power-law behavior results from the fact that at \(\nu < 1\) (or \(\sigma > \delta_c\)), \(f_{\mathrm{PS}}(\nu)\) is approximately constant with \(M\), such that the only mass dependence comes from the \(\bar{\rho}_{0,m} / M\) factor in Equation (17.121). Physically, this makes sense: while we can compute the variance of \(\delta_S(\vec{x},t;M)\) in linear theory when \(\sigma(M,t) > \delta_c\), in reality such regions have gravitationally (and non-linearly!) collapsed into dark-matter halos and, therefore, their abundance does not depend on how much larger \(\sigma(M,t)\) is than \(\delta_c\). Once \(\sigma(M,t) \gg \delta_c\), the only evolution in the halo mass function comes from the merging of smaller dark-matter halos into larger ones. This causes the lowest-mass end of the mass function to decline as time goes on (and redshift decreases). But the overall similarity of the low-mass power-law behavior at low masses demonstrates that merging does not have a huge effect on the overall halo mass function. At \(M \gtrsim M^*(t)\), the halo mass function does sensitively depend on \(\nu\) or, equivalently, on how much smaller \(\sigma(M,t)\) is than \(\delta_c\).
The halo mass function derived by Press & Schechter (1974) is in reasonable agreement with halo mass functions computed directly using \(N\)-body simulations of structure formation, but in detail does not agree. Over the last few decades, there have therefore been many studies that use \(N\)-body simulations to empirically determine the halo mass function from the number of dark-matter halos that form in such simulations. Generally, such halo mass functions are defined through an alternative form for the function \(f_{\mathrm{PS}}(\nu)\) that appears in Equation (17.121), because simulations demonstrate that this function is remarkably independent of the background cosmology. The most useful one of these alternative halo mass functions is that of Sheth & Tormen (1999), who proposed the alternative multiplicity function \begin{equation}\label{eq-dmhaloform-hmf-fst} f_{\mathrm{ST}}(\nu) = 0.3222\,\sqrt{{2\over \pi}}\,\nu'\,e^{-\nu'^2/2}\,\left(1+\nu'^{-0.6}\right)\,, \end{equation} where \(\nu' = \sqrt{0.707}\nu\). While this function was empirically derived using \(N\)-body simulations, Sheth et al. (2001) later showed that this general form of \(f(\nu)\) naturally arises when taking into account deviations from spherical symmetry in the modeling of the collapse process (that is, through generalizing the spherical-collapse model to ellipsoidal collapse), although the derived equation is not quite that of Equation (17.123) (instead it is Equation 17.123 with \(\nu'=\nu\)). We implement the Sheth & Tormen (1999) halo mass function and then compare it to the Press & Schechter (1974) one at \(z=0\) in Figure 17.10.
[20]:
from astropy.constants import G
rhom= Om0*3.*H0**2./8./numpy.pi/G
def halo_mass_function(M,a=1.,model='PS',dlnsigdlnM=None,sig=None):
# d n / d ln M
# Use sigma and d ln sigma / d ln M if provided, otherwise compute
if sig is None:
sig= numpy.sqrt(sigma2(M,a=a))
if dlnsigdlnM is None:
dlnsigdlnM= dlnsigmadlnM(M)
# f function
if model.lower() == 'ps':
def fnu(nu):
return numpy.sqrt(2./numpy.pi)*nu*numpy.exp(-nu**2./2.)
elif model.lower() == 'st':
def fnu(nu):
nu= numpy.sqrt(0.707)*nu
return 0.3222*numpy.sqrt(2./numpy.pi)*nu*numpy.exp(-nu**2./2.)\
*(1.+nu**-0.6)
nu= 1.686/sig
return -rhom/M*fnu(nu)*dlnsigdlnM
figure(figsize=(8,5))
loglog(ms,
u.Quantity(halo_mass_function(ms,a=1.,model='ps',
sig=sig_z0,dlnsigdlnM=dlnsigdlnM))\
.to_value(1./u.Mpc**3),
label=r'$\mathrm{Press\ \&\ Schechter\ (1974)}$')
loglog(ms,u.Quantity(halo_mass_function(ms,model='st',
sig=sig_z0,
dlnsigdlnM=dlnsigdlnM))\
.to_value(1./u.Mpc**3),
label=r'$\mathrm{Sheth\ \&\ Tormen\ (1999)}$')
xlabel(r'$M/M_\odot$')
ylabel(r'$\mathrm{d} n/\mathrm{d} \ln M\,(\mathrm{Mpc}^{-3})$')
ylim(1e-10,4e4)
legend(frameon=False,fontsize=18.);
Figure 17.10: The Sheth-Tormen halo mass function.
While overall the two halo mass functions are very similar, the large range on the \(y\) axis obscures the fact that the differences are in fact relatively large and important for quantitative use of the halo mass function. The Sheth & Tormen (1999) mass function, however, remains a workhorse in practical research use of the halo mass function.
Because the number of dark-matter halos of different masses is an essential ingredient of any cosmological investigation regarding galaxies, the Press-Schechter halo mass function has had an enormous impact on the field of galaxy formation and evolution. To give just a few examples, the Press-Schechter halo mass function allowed for studies of the abundance and properties of rare luminous quasars (e.g., Efstathiou & Rees 1988), investigations of the importance of hierarchical clustering for galaxy formation (e.g., White & Frenk 1991), models for the formation of galactic disks (e.g., Mo et al. 1998), and, through the extension of the Press-Schechter formalism to also include mergers (Bond et al. 1991), studies of the merger rate (e.g., Lacey & Cole 1993).
.
.
.