19.2. Dark-matter halo growth by mergers and accretion¶
As we discussed in Chapter 17, in the \(\Lambda\)CDM cosmological framework, the dark-matter halos that host galaxies form when initial overdensities in the primordial density field grow so large that they gravitationally collapse and form a virialized object. These dark-matter halos are not static, but continue to grow through mergers with other virialized dark-matter halos and through accretion of diffuse dark matter. Galaxies grow within dark-matter halos from the gas that was part of the initial collapse and from stars and gas brought in through mergers and accretion. Thus, studying the growth of dark-matter halos after their initial formation is crucial to understanding galaxy evolution.
Because dark matter only interacts through gravity—at least on cosmological scales—determining the evolution of dark-matter halos through merging and accretion from an initial density field is in principle straightforward. All one needs to do is evolve the initial density field up to the present time, identify halos at different times along the way, and record if and when they merge. For a given cosmological model and initial density field, evolving the density field under the effect of gravity in an expanding Universe—typically using \(N\)-body codes—has no free parameters and, therefore, has an unambiguous, correct solution (however, numerical algorithms used to compute the evolution, do have free parameters; see Chapter 12.4). Real dark-matter halos are not simple objects with hard boundaries, so identifying halos in simulations is not unambiguous and a large number of algorithms for finding halos has been developed (e.g., friends-of-friends-type algorithms: e.g., Davis et al. 1985; algorithms that look for spherical overdensities around density peaks: e.g., Klypin et al. 1999; and complex algorithms that consider the full six-dimensional phase-space and temporal evolution of particles in a simulation: e.g., Behroozi et al. 2013b). But aside from details important for precision-cosmology applications (but less so for galaxy evolution), we can identify halos in \(N\)-body simulations of structure formation and extract their merger and accretion history. Because advances in computational infrastructure have made it such that running large, dark-matter-only \(N\)-body simulations of structure formation can be done quickly and easily, \(N\)-body simulations are the gold standard for studying the evolution of dark-matter halos nowadays.
Much about the evolution of dark-matter halos can be learned, however, from an analytical approach first introduced by Bond et al. (1991): the extended Press-Schechter (EPS) formalism, also known as the excursion-set formalism. As the name implies, this formalism is an extension of the basic Press-Schechter formalism for determining the halo mass function that we discussed in Chapter 17.4; the extension allows it to determine the full merger and accretion history of a target halo. It does this by re-formulating the Press-Schechter formalism in terms of excursion sets of Gaussian fluctuations (hence the alternative name). The advantage of the EPS formalism over the \(N\)-body approach is that it is blazingly fast: rather than spending weeks of supercomputer time to run a single cosmological \(N\)-body simulation using a complex, highly-optimized code, the EPS formalism can determine the full merger history of a dark-matter halo in mere seconds with only a few dozen lines of code. Because it is much faster, while still being accurate enough for most applications in galaxy evolution, the EPS formalism is useful for experimenting with deviations from the fiducial cosmology used to run large \(N\)-body simulations (usually, the cosmology specified by the cosmological parameters from the latest CMB experiment; e.g., Planck Collaboration et al. (2020b) at the time of writing). Because the EPS formalism allows one to focus on the history of a single halo, it is also helpful for quickly determining the statistical properties of its merger history, e.g., the expected number of major and minor mergers. Because of its good accuracy and speed, we discuss the EPS formalism in detail below to understand the important properties of the merger and accretion history of dark-matter halos.
Before we dive into the EPS formalism, it is important to be clear about its limitations. The EPS formalism is most useful for determining how much mass a dark-matter halo has gained through mergers and accretion, especially over limited redshift ranges where it compares well with \(N\)-body simulations. However, it does not directly tell us much more about these mergers, for instance, the orbit of the merging galaxies, which determines where any accreted material will be deposited. For gas accreted from the environment, it does not tell us what happens to this gas as it makes its way from the intergalactic medium to the central galaxy. Obtaining this information and more from the EPS formalism is the domain of semi-analytic galaxy formation (White & Frenk 1991; Kauffmann et al. 1993; Cole et al. 2000; Benson 2012), which adds prescriptions for gas cooling, star formation, merger orbits, etc. to the EPS formalism. The development of large-volume, cosmological hydrodynamical simulations of galaxy formation in the last decade (e.g., Illustris: Vogelsberger et al. 2014; EAGLE: Crain et al. 2015) has largely made semi-analytic models obsolete. Hydrodynamical simulations can provide the same information as semi-analytic galaxy formation models using modeling that is more directly connected to the underlying physics, although it still requires parameterizing processes like star formation and supernova feedback that cannot be directly self-consistently modeled.
19.2.1. The extended Press-Schechter formalism and the halo mass function¶
The EPS formalism is based on a rigorous derivation of the Press-Schechter approach from Chapter 17.4 by Bond et al. (1991). This derivation starts with the following insight. When we consider the density field \(\delta_S(\vec{x},t;R)\) at a given position \(\vec{x}\) smoothed using a window function of size \(R\), then for a generic window function, the variance of the difference between \(\delta_S(\vec{x},t;R)\) and the density field \(\delta_S(\vec{x},t;R-\Delta R)\) smoothed on a smaller scale \(R-\Delta R\), depends on the entire power spectrum (see Equation 17.116): \begin{align} \langle & \left[\delta_S(\vec{x},t;R-\Delta R)-\delta_S(\vec{x},t;R)\right]^2\rangle \nonumber\\ &\ = \langle \delta^2_S(\vec{x},t;R-\Delta R)\rangle+\langle \delta^2_S(\vec{x},t;R)\rangle-2\,\langle \delta_S(\vec{x},t;R-\Delta R)\,\delta_S(\vec{x},t;R)\rangle\,,\\ &\ = \sigma^2(R-\Delta R,t)+\sigma^2(R,t)-2\,{1\over 2\pi^2}\int_0^\infty\mathrm{d}k\,k^2\,P_m(k,a)\,W_k(R)\,W_k(R-\Delta R)\,.\label{eq-mergers-eps-dS-in-terms-of-powerspec} \end{align} However, when we use a window function \(W_k(R)\) that is a sharp \(k\)-space filter—constant up to \(1/R\) and zero otherwise, that is, \(W_k(R) = \Theta(1-kR)\), where \(\Theta\) is the Heaviside function (Equation B.31)—then this simplifies to \begin{align} \langle & \left[\delta_S(\vec{x},t;R-\Delta R)-\delta_S(\vec{x},t;R)\right]^2\rangle = \sigma^2(R-\Delta R,t)-\sigma^2(R,t) = {1\over 2\pi^2}{\large\int}_{1/R}^{1/(R-\Delta R)}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\mathrm{d}k\,k^2\,P_m(k,a)\,,\label{eq-mergers-eps-dS-in-terms-of-powerspec-sharpk} \end{align} because we can replace \(W_k(R-\Delta R)\) with \(W_k(R)\) in Equation (19.2) as the integrand is zero where they differ. Thus, for a sharp \(k\)-space filter, the variance of the difference between smoothing the density field in a given location on scales \(R\) and \(R-\Delta R\) only depends on the statistics of density fluctuations on the intermediate scales \(1/(R-\Delta R) < k < 1/R\). Because of this property, we can draw realizations for the smoothed linear density field \(\delta_S(\vec{x},t;R-\Delta R)\) starting at a given \(\delta_S(\vec{x},t;R)\) by adding an increment drawn from a Gaussian with zero mean and variance given by Equation (19.3). For a generic window function, the distribution of this increment additionally depends on the value of the density field smoothed on all scales larger than \(R\).
Because dark-matter halos in the Press-Schechter formalism are identified with regions whose density, when smoothed over some scale \(R\), exceed the critical density \(\delta_c\) for collapse in the spherical-collapse picture (see Chapter 17.2), the question of whether a mass element at location \(\vec{x}\) is part of a dark-matter halo with size larger than \(R\) becomes whether \(\delta_S(\vec{x},t;R) \geq \delta_c(t)\). As is standard in this field, we denote \(S \equiv \sigma^2(R,t)\), where \(S\) is understood to depend on scale, but this dependence is kept implicit because it need not be known for much of what follows. In hierarchical models such as the \(\Lambda\)CDM paradigm, there is a one-to-one mapping between \(R\) and \(S\) (see Figure 17.8) and we can, thus, denote the smoothed density field as \(\delta(\vec{x},S) \equiv \delta_S(\vec{x},t;R[S])\). For the sharp \(k\)-space filter, we can easily draw values for \(\delta(\vec{x},S)\) from its expected distribution by simulation the following trajectories as a function of \(S\).
When smoothed on the scale \(R = \infty\), the overdensity \(\delta_S(\vec{x},t;\infty) = 0\) for all \(\vec{x}\). Because \(S = \sigma^2(\infty,t) = 0\) (see Equation 17.116), every trajectory, therefore, starts at \(\delta(\vec{x},S) = 0\) at \(S=0\).
To move from \(\delta(\vec{x},S)\) to \(\delta(\vec{x},S+\Delta S)\), we draw an increment \(\Delta \delta(\Delta S)\) from a zero-mean Gaussian with variance \(\Delta S\) and obtain \(\delta(\vec{x},S+\Delta S) = \delta(\vec{x},S) + \Delta \delta(\Delta S)\).
This exactly corresponds to a random walk, also known as Brownian motion. Thus, to simulate values of \(\delta(\vec{x},\tilde{S})\), we draw random walks from \(S=0\) to \(S=\tilde{S}\). Figure 19.4 demonstrates three example realizations.
[6]:
from numpy.random import default_rng
rng= default_rng(17171)
# First we implement a straightforward random walk (because we use it below)
# then the walk with an absorbing barrier
def random_S_walk(nS,dS=0.1):
S= numpy.linspace(dS,(nS-1)*dS,nS-1)
out= numpy.cumsum(numpy.hstack(\
([0.],rng.standard_normal(size=len(S))*numpy.sqrt(dS))))
return numpy.hstack(([0.],S)),out
def random_S_walk_absorb(nS,dS=0.1,absorb=1.686,return_unabsorbed=False):
S= numpy.linspace(dS,(nS-1)*dS,nS-1)
out= numpy.cumsum(numpy.hstack(\
([0.],rng.standard_normal(size=len(S))*numpy.sqrt(dS))))
if return_unabsorbed: # in this case, we also return the unabsorbed walk and the index where it starts to differ
raw= out.copy()
# argmax stops at first match
out[numpy.argmax(numpy.hstack((out,[absorb]))>=absorb):]= absorb
if return_unabsorbed:
return numpy.hstack(([0.],S)),out,raw,numpy.argmax(numpy.hstack((out,[absorb]))>=absorb)
else:
return numpy.hstack(([0.],S)),out
figure(figsize=(8,5))
for ii in range(3):
s,d,rd,indx= random_S_walk_absorb(100,return_unabsorbed=True)
line = plot(s,d,zorder=-ii)
plot(s[indx:],rd[indx:],color=line[0].get_color(),ls='--',zorder=-ii)
axhline(1.686,color='k',lw=1.5,zorder=-100)
ylim(-3.7,3.7)
text(0.15,1.9,r'$\delta(\boldsymbol{{x}},S)=\delta_c$',fontsize=18.)
xlabel(r"$S$")
ylabel(r"$\delta(\boldsymbol{{x}},S)$");
Figure 19.4: Brownian motion trajectories with (solid) and without (solid \(\rightarrow\) dashed) an absorbing barrier.
We also indicate the line \(\delta(\vec{x},S) = \delta_c\) where the smoothed linear density field exceeds the threshold for collapse in the spherical-collapse framework. Following the Press-Schechter approach, dark-matter halos form when \(\delta(\vec{x},S) \geq \delta_c\) and, thus, to simulate halo formation using the random walks shown above, we need an additional ingredient: that linear density evolution ceases when \(\delta(\vec{x},S) \geq \delta_c\). Once a dark-matter halo collapses, its overdensity on smaller scales (larger \(S\)) far exceeds \(\delta_c\) and it cannot, in particular, become smaller than \(\delta_c\). Thus, trajectories such as the dashed one in Figure 19.4 cannot occur, but instead we can simply set \(\delta(\vec{x},S) = \delta_c\) for any \(S\) larger than the smallest \(S\) where \(\delta(\vec{x},S) \geq \delta_c\). This is illustrated in Figure 19.4 as the solid plateaus for the two uppermost trajectories (the third one does not exceed \(\delta_c\) for \(S<10\)). Thus, we see that the two uppermost trajectories collapse and form halos, the first one at \(S \approx 4.5\) and the second one at \(S \approx 9.5\). Note that we do not need any knowledge of the cosmological model for generating these trajectories, the cosmological model only enters in how \(S\) translates to \(R\) (and subsequently to halo mass \(M\)). The line \(\delta(\vec{x},S) = \delta_c\) acts as a barrier that absorbs any trajectory that crosses it and this stochastic process is known as a random walk with an absorbing barrier. For a random walk, it is guaranteed that \(\delta(\vec{x},S)\geq \delta_c\) for some \(S\) when we let \(S \rightarrow \infty\), so all mass elements in the EPS formalism are automatically part of a halo. Remember that in the original Press-Schechter formalism in Chapter 17.4, we needed to introduce a factor of two, because only half of the mass elements are naturally part of a halo in that formalism.
In the EPS formalism, the Press-Schechter ansatz is that the fraction of mass elements contained in halos with sizes less than \(R\) is equal to the fraction of trajectories with \(\delta(\vec{x},S[R]) < \delta_c\). In the absence of the absorbing barrier, the distribution of \(\delta(\vec{x},S)\) is simply a Gaussian \begin{equation} p(\delta[\vec{x},S]) = {1 \over \sqrt{2\pi\,S}}\,\exp\left(-{\delta^2[\vec{x},S]\over 2S}\right)\,, \end{equation} but in the presence of an absorbing barrier at \(\delta_c\), this is modified to (Chandrasekhar 1943b) \begin{equation} p(\delta[\vec{x},S]) = {1 \over \sqrt{2\pi\,S}}\,\left[\exp\left(-{\delta^2[\vec{x},S]\over 2S}\right)-\exp\left(-{\{2\delta_c-\delta[\vec{x},S]\}^2\over 2S}\right)\right]\,, \end{equation} The fraction of mass elements in halos with sizes smaller than \(R\) is then given by the fraction of trajectories that do not exceed \(\delta_c\) at \(S(R)\) \begin{align} P[\delta(\vec{x},S[R]) < \delta_c] & = \int_{-\infty}^{\delta_c}\mathrm{d} \left(\delta[\vec{x},S] \right)\,p(\delta[\vec{x},S])= \mathrm{erf}\left({\delta_c \over \sqrt{2S}}\right)\,, \end{align} where \(\mathrm{erf}(\cdot)\) is the standard error function from Equation (B.17). Because all mass elements are contained in halos, the fraction of mass elements in halos with sizes larger than \(R\) is given by \begin{align}\label{eq-mergers-eps-deltaS-le-deltac} P[\delta(\vec{x},S[R]) \geq \delta_c] & =\mathrm{erfc}\left({\delta_c \over \sqrt{2S}}\right)\,, \end{align} where \(\mathrm{erfc}(\cdot)\) is the complementary error function (see Equation B.19). This fraction is exactly the same as that obtained in the original Press-Schechter formalism (Equation 17.118), after accounting for the magical factor of two in the original Press-Schechter formalism. The remainder of the derivation of the halo mass function is then the same as before and we end up at Equation (17.120) again. Because trajectories are absorbed once they reach \(\delta(\vec{x},S) = \delta_c\), it is also clear that the halo mass function determined in this way applies to main halos only and does not include subhalos within larger halos.
Thus, the EPS formalism demonstrates that the Press-Schechter halo mass function of Equation (17.120) is exact in the case of the sharp \(k\)-space filter. It does not prove that it is correct for any other filter and, in fact, shows the opposite, because the derivation above depends crucially on the random-walk nature of \(\delta(\vec{x},S)\) trajectories for a sharp \(k\)-space filter and this assumption is broken for all other filters. For all other filters, the trajectories are not even Markovian, meaning that the increment at \(\delta(\vec{x},S)\) depends on values of \(\delta(\vec{x},S')\) for \(S' < S\). While such trajectories can be simulated using a diffusion-equation approach to the above derivation (e.g., Bond et al. 1991) and we can then numerically find \(P[\delta(\vec{x},S[R]) < \delta_c]\) in the presence of an absorbing barrier (e.g., Percival 2001), this is only of limited interested, because in the end any result from the EPS formalism needs to be validated and calibrated using \(N\)-body simulations. In practice, rigorous results from the EPS formalism derived under the assumption of a sharp \(k\)-space filter are typically applied to other window functions that are more suited to defining dark-matter halos (e.g., the top-hat window function that we used in Chapter 17.4). Comparing the results from this not-entirely-rigorous procedure to \(N\)-body simulations demonstrates that this works well in practice (e.g., we saw in Chapter 17.4 that the Press-Schechter halo mass function computed using a top-hat windows is close to the Sheth & Tormen (1999) one that was fit to \(N\)-body simulations).
We can compare the halo mass function computed using the sharp \(k\)-space filter to that obtained using the top-hat filter in Chapter 17.4, but in order to do that, we first need to discuss how we relate the size \(R\) of the filter to the mass \(M\) of the halo. For the top-hat window in Chapter 17.4, this correspondence was simple, because we could set \(M = 4\pi \bar{\rho}_{0,m} R^3/3\) where \(\bar{\rho}_{0,m}\) is the mean matter density at redshift zero and this led to Equation (17.114). For general windows, \(M\) should still scale like \(R^3\), but with a different prefactor and this is usually denoted as \(M = \gamma R^3\). The factor \(\gamma\) can be set to the volume of the \(R=1\) window in position space—this leads to \(\gamma = 4\pi/3\) for the top-hat window—but the sharp \(k\)-space filter has infinite volume, so there is no obvious correct setting for \(\gamma\). Instead, \(\gamma\) is often set to the inverse of the real-space window function evaluated at the origin—again equal to \(4\pi/3\) for the top hat—which gives \(\gamma = 6\pi^2\) for the sharp \(k\)-space filter. With that choice, we would find \(R = \left(4\,GM / [9\pi\,\Omega_{0,m}\,H_0^2]\right)^{1/3} = 0.75\,\mathrm{Mpc}\,\left(M / 10^{12}\,M_\odot\right)^{1/3}\). However, we are free to choose a different pre-factor, because \(\gamma\) is not predicted by the EPS formalism, so instead we will use \begin{align}\label{eq-merger-PS-radius-for-mass-sharpk-adopted} R & = 0.9\,\mathrm{Mpc}\,\left({M \over 10^{12}\,M_\odot}\right)^{1/3}\,, \end{align} because this leads to approximately the same halo mass function.
To compute \(\sigma(M,t)\) and the resulting halo mass functions, we use the code written in Chapter 17.4, but edited to allow both the top-hat and sharp \(k\)-space window functions. First we calculate \(\sigma(M,t)\) using Equation (17.116):
[7]:
# 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,kind='tophat'):
if kind.lower() == 'tophat':
# Top-hat window function in Fourier space
return 3.*(numpy.sin(x)-x*numpy.cos(x))/x**3.
elif kind.lower() == 'sharpk':
# Sharp k-space filter
return float(x <= 1.)
def sigma2_R(R,a=1.,window_kind='tophat'):
# sigma2 computed as a function of R at scale factor a
kmax= (1./R).to_value(1./u.Mpc) \
if window_kind.lower() == 'sharpk' \
else numpy.inf
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),
kind=window_kind)**2,
0,kmax)[0]/2./numpy.pi**2.*growth_factor(a)**2.
def M_to_R(M,window_kind='tophat'):
if window_kind.lower() == 'tophat':
R= 1.82*u.Mpc*(M/1e12/u.Msun)**(1./3.)
elif window_kind.lower() == 'sharpk':
R= 0.9*u.Mpc*(M/1e12/u.Msun)**(1./3.)
return R
def sigma2(M,a=1.,window_kind='tophat'):
# sigma2 computed as a function of mass M
return sigma2_R(M_to_R(M,window_kind=window_kind),
a=a,window_kind=window_kind)
def rescale_sigma(sig,a=1.,aref=1.):
# rescale sigma to another scale factor
return sig*growth_factor(a=a)/growth_factor(a=aref)
To compute the halo mass function, we use Equation (17.120). The calculation of \(\mathrm{d} \sigma^{-1} /\mathrm{d} \ln M\) in that equation is slightly complicated, because it involves the derivative of the window function, which is a Heaviside function for the sharp \(k\)-space filter and its derivative is a delta function (Equation B.32). The relevant integral then has the product of a Heaviside function and a delta function, which can be computed using Equation (B.33). The result is \begin{equation} {\mathrm{d} \ln \sigma(M,t) \over \mathrm{d} \ln M} = -{1\over 12\pi^2}\,R^{-3}\,P_m(R^{-1})\,\sigma^{-2}(M)\,, \end{equation} where all quantities on the right-hand side are evaluated at redshift zero, because this derivative is independent of redshift.
Putting this all together, the following code computes the halo mass function for the top-hat and sharp \(k\)-space window functions:
[8]:
from astropy.constants import G
rhom= Om0*3.*H0**2./8./numpy.pi/G
def dwindowdx(x,kind='tophat'):
if kind.lower() == 'tophat':
return 3.*((3.*x*numpy.cos(x)+(x**2.-3.)*numpy.sin(x))/x**4.)
elif kind.lower() == 'sharpk':
return 0. # Never actually used
def dlnsigmadlnM(M,window_kind='tophat'):
# constant with a, so set a=1.
R= M_to_R(M,window_kind=window_kind)
if window_kind.lower() == 'sharpk':
return -(1./R.to_value(u.Mpc))**3.\
*matter_power_spectrum_nogrowth(1./R,
model='Sugiyama').to_value(u.Mpc**3)\
/12./numpy.pi**2.*growth_factor(1.)**2.\
/sigma2_R(R,a=1.,window_kind=window_kind)
else:
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),
kind=window_kind)\
*dwindowdx((k*R/u.Mpc).to_value(u.dimensionless_unscaled),
kind=window_kind),
0,numpy.inf)[0]\
/6./numpy.pi**2.*growth_factor(1.)**2.\
/sigma2_R(R,a=1.,window_kind=window_kind)
def halo_mass_function(M,a=1.,window_kind='tophat',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,window_kind=window_kind))
if dlnsigdlnM is None:
dlnsigdlnM= dlnsigmadlnM(M,window_kind=window_kind)
# 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
The resulting redshift zero halo mass function is then shown in Figure 19.5:
[9]:
ms= numpy.geomspace(1e5,1e16,51)*u.Msun
figure(figsize=(8,5))
loglog(ms,
u.Quantity([halo_mass_function(m,a=1.,window_kind='tophat')
for m in ms]).to_value(1./u.Mpc**3),
label=r'$\mathrm{Top\!-\!hat}$')
loglog(ms,
u.Quantity([halo_mass_function(m,a=1.,window_kind='sharpk')
for m in ms]).to_value(1./u.Mpc**3),
label=r'$\mathrm{Sharp}\,k$')
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 19.5: The Press-Schechter halo mass function for different window functions.
Thus, we see that the halo mass functions computed using the top-hat and sharp \(k\)-space window functions agree very well. This helps justify using the results from the EPS formalism with \(S\) computed using a top-hat window function rather than the sharp \(k\)-space one used to derive the results. This is the approach we will take below.
19.2.2. Dark-matter halo accretion histories and merger trees¶
While the EPS formalism is useful for deriving the results from Press-Schechter theory in a (more) rigorous way, where it really shines is in determining the full merger and accretion history of a dark-matter halo. This is possible, because the random-walk formalism from above can be re-interpreted as giving the trajectory of a given mass element over time. We can then use these trajectories to determine whether a mass element that is part of a halo of mass \(M_1\) at, say, the present time \(t_1\) was part of a halo with mass \(M_2 < M_1\) at an earlier time \(t_2\). To achieve this, the first insight that we need is that while we use \(\delta_c\) and \(\sigma(M,t)\), all of the relevant equations in the EPS formalism only depend on their ratio \(\delta_c/\sigma(M,t)\). Because the time dependence of \(\sigma(M,t)\) is simply that it is proportional to the growth-factor \(D_+(a)\) (see Equations 17.79 and 17.116), we can factor out the growth factor and absorb it into a re-definition of \(\delta_c \rightarrow \delta_c(a) = \delta_c / [D_+(a)/D_+(1)]\) and then use \(\sigma(M) \equiv \sigma(M,a=1)\) and \(\delta_c(a)\) instead of \(\sigma(M,t)\) and \(\delta_c\). Note that the original \(\delta_c\) already depends on time in our \(\Lambda\)CDM Universe (Equation 17.93), but this time dependence is so weak that we ignore it. The only time dependence of the transformed \(\delta_c(a)\) then comes from the growth factor; to clearly distinguish whether we mean the original \(\delta_c\) or the transformed \(\delta_c(a)\), we will always explicitly include the scale-factor dependence of the latter.
To determine the largest halo that a mass element is part of at redshift \(z\) (scale factor \(a\)), we simulate random-walk trajectories exactly as we did above, but we move the absorbing barrier up to \(\delta_c(a) = \delta_c/[D_+(a)/D_+(1)]\). With the re-definition of \(\delta_c\), we have separated the mass and time dependence of the EPS formalism in the sense that the mass dependence is entirely contained in the variable \(S = \sigma^2(M)\) evaluated at redshift zero and only \(\delta_c(a)\) depends on time. Below, it will often be useful to use \(S\) as the mass variable and \(\delta_c(a)\) as the time variable.
To then reinterpret the random-walk trajectories that we simulated above as merger histories for mass elements, we go back to the original random-walk trajectory without the barrier and we now interpret any upward movement of \(\delta(\vec{x},S)\) when it is above \(\delta_c\) as meaning that the mass element resides in a halo with a mass corresponding to \(S\) at a scale factor set by \(\delta(\vec{x},S)=\delta_c(a)\). To illustrate this, we simulate another random-walk trajectory. To easily translate between \(S\) and \(M\) and \(\delta_c(a)\) and \(a\) (or redshift \(z\)), we first set up interpolating functions \(M(S)\) and \(z(\delta_c[a])\):
[10]:
from scipy import interpolate
ms_sigma= numpy.geomspace(1e5,1e16,21)*u.Msun
sigma_sharpk= numpy.array([numpy.sqrt(sigma2(m,a=1.,window_kind='sharpk'))
for m in ms_sigma])
_M_from_S= interpolate.InterpolatedUnivariateSpline(\
sigma_sharpk[::-1]**2.,
numpy.log(ms_sigma[::-1].to_value(u.Msun)),
k=3,ext=0)
M_from_S= lambda S: numpy.exp(_M_from_S(S))*u.Msun
Dpluszero= growth_factor(1.)
as_for_growthfactor_interp= numpy.geomspace(1e-2,1.,101)[::-1]
z_from_deltac= interpolate.InterpolatedUnivariateSpline(\
[1.686/(growth_factor(a)/Dpluszero)
for a in as_for_growthfactor_interp],
1./as_for_growthfactor_interp-1.,
k=3,ext=0)
Then we look at the following trajectory shown as the solid line:
[11]:
from scipy.signal import argrelextrema
def walk_to_growth(walk,return_plateaus=False):
S, deltaS= walk
# Find all local maxima
indx= argrelextrema(deltaS,numpy.greater)[0]
# Remove those below the z=0 threshold
indx= indx[deltaS[indx] >= 1.686]
# Remove those that are smaller than the previous one
delta_locmax= (numpy.hstack((deltaS[indx],0.))
-numpy.roll(numpy.hstack((deltaS[indx],0.)),1))[:-1]
while numpy.any(delta_locmax < 0.):
indx= indx[delta_locmax >= 0.]
delta_locmax= (numpy.hstack((deltaS[indx],0.))
-numpy.roll(numpy.hstack((deltaS[indx],0.)),1))[:-1]
if numpy.sum(indx) == 0:
# no local maxima above the 1.686 peak
if return_plateaus:
return S, deltaS, [], []
else:
return S, deltaS
# Now replace any values below the previous local maximum with the value
# of the local maximum
# First we create an array that has the value of the previous local maximum
locmax_freq= (numpy.roll(numpy.hstack((indx,len(S))),-1)
-numpy.hstack((indx,len(S))))[:-1]
locmax_arr= numpy.hstack(([[lmax for ii in range(lmax_freq)]
for lmax,lmax_freq in zip(deltaS[indx],locmax_freq)]))
# Add deltas before the first local maximum at the start
locmax_arr= numpy.hstack((deltaS[:indx[0]],locmax_arr))
# Now we can easily replace any values below the previous local maximum
# with the value of the local maximum
out= deltaS.copy()
out[deltaS < locmax_arr]= locmax_arr[deltaS < locmax_arr]
if not return_plateaus:
return S,out
# Also return the boundaries of the plateaus (=mergers) if requested
tot_indx_arr= numpy.arange(len(S))
plateau_start_indx= tot_indx_arr[deltaS < locmax_arr]\
[(tot_indx_arr[deltaS < locmax_arr]
-numpy.roll(tot_indx_arr[deltaS < locmax_arr],1)) > 1]
# Add the first one
plateau_start_indx= numpy.hstack(([tot_indx_arr[deltaS < locmax_arr][0]],
plateau_start_indx))
plateau_end_indx= tot_indx_arr[deltaS < locmax_arr]\
[(numpy.roll(tot_indx_arr[deltaS < locmax_arr],-1)
-tot_indx_arr[deltaS < locmax_arr]) > 1]
#Add the last one
plateau_end_indx= numpy.hstack((plateau_end_indx,
[tot_indx_arr[deltaS < locmax_arr][-1]]))
return S,out,plateau_start_indx,plateau_end_indx
rng= default_rng(82828)
walk= random_S_walk(300)
figure(figsize=(8,5))
plot(*walk,zorder=10)
line_growth = plot(*walk_to_growth(walk),ls='--',lw=3.)
line_growth[0].set_dashes([2,2,2,2])
# Also label mergers
S,deltaS,plateau_start,plateau_end= walk_to_growth(walk,return_plateaus=True)
for pl_start,pl_end in zip(plateau_start,plateau_end):
if pl_end-pl_start < 3: continue
axvspan(S[pl_start],S[pl_end],color='0.8')
axhline(1.686,color='k',lw=1.5)
t= text(20.,2.35,r'$\delta(\boldsymbol{x},S)=\delta_c(a=1)$',fontsize=18.)
t.set_bbox(dict(facecolor='w',alpha=0.75,edgecolor='none'))
xlabel(r"$S$")
ylabel(r"$\delta(\boldsymbol{x},S)$")
# Alternative scales
from matplotlib.ticker import FuncFormatter
ax_orig= gca()
# Mass on alt x
ax_twiny= ax_orig.twiny()
old_ticks= ax_orig.get_xticks()
ax_twiny.set_xticks(old_ticks[old_ticks > 0.])
def mass_label(x,pos):
M= M_from_S(x).to_value(u.Msun)
expo= numpy.abs(numpy.floor(numpy.log10(numpy.fabs(M))))
return rf"${M/10**expo:.0f}\times 10^{{{expo:.0f}}}$"
ax_twiny.xaxis.set_major_formatter(FuncFormatter(mass_label))
ax_twiny.set_xlim(ax_orig.get_xlim()) # has to be ~last
ax_twiny.set_xlabel(r"$M/M_\odot$")
# z on alt y
ax_twinx= ax_orig.twinx()
old_ticks= ax_orig.get_yticks()
dold_ticks= old_ticks[1]-old_ticks[0]
nticks= numpy.ceil((old_ticks[-1]-1.7)/dold_ticks)
ax_twinx.set_yticks(numpy.arange(1.7,(nticks+1)*dold_ticks,dold_ticks))
def z_label(x,pos):
return rf"${z_from_deltac(x):.1f}$"
ax_twinx.yaxis.set_major_formatter(FuncFormatter(z_label))
ax_twinx.set_ylim(ax_orig.get_ylim()) # has to be ~last
ax_twinx.set_ylabel(r"$\mathrm{redshift}$");
Figure 19.6: A Brownian-motion trajectory interpreted as the merger history of a mass element in the Universe.
To allow one to easily read off the history of the mass element, we also include the redshift corresponding to each \(\delta_c(a) \geq \delta_c\) and the mass corresponding to \(S\) as alternative \(y\) and \(x\) axes, respectively.
We see that the trajectory crosses \(\delta(\vec{x},S)=\delta_c = \delta_c(a=1)\) at \(S \approx 5\). Therefore, at \(z=0\), the mass element is part of a halo of mass \(M \approx 2\times 10^{12}\,M_\odot\) in our \(\Lambda\)CDM cosmology. Then the trajectory meanders up and down further and, for example, reaches \(\delta_c(a) \approx 7 \approx 1.686/0.24\) at \(S\approx 10\), meaning that at \(z\approx 4.3\)—when \(D_+(a)/D_+(1) \approx 0.24\)—the mass element was part of a halo with mass \(M \approx 1.4\times 10^{11}\,M_\odot\).
Like before, downward movements in \(\delta(\vec{x},S)\) once it has exceeded \(\delta_c\) are unphysical, because in our new interpretation, these would correspond to a mass element being part of a smaller halo at a lower redshift than we had determined earlier in the trajectory. For example, in Figure 19.6, we see that the trajectory first crosses \(\delta_c(a) = 4.5\) at \(S\approx 7.5\), corresponding to \(M\approx 4.4\times 10^{11}\,M_\odot\). But because the trajectory turns down at \(S\approx 10\), it crosses \(\delta_c(a) = 4.5\) again at \(S\approx 12.5\), corresponding to \(M\approx 5\times 10^{10}\,M_\odot\). Thus, we would conclude that the mass element is both in a halo of mass \(M\approx 4.4\times 10^{11}\,M_\odot\) and of mass \(M\approx 5\times 10^{10}\,M_\odot\) at the redshift \(z\approx 2.4\) that corresponds to \(\delta_c(a) = 4.5\). To resolve this inconsistency, we proceed similarly as above: we do not allow the trajectory to decline once it has exceeded \(\delta(\vec{x},S) = \delta_c = 1.686\). For our example trajectory, we do this in Figure 19.6 by converting it to the dashed curve.
We see that once \(\delta(\vec{x},S) \geq \delta_c\), the mass element’s trajectory either steadily increases when the original random-walk trajectory increases or the mass element’s trajectory remains constant where the original trajectory decreases. Because the latter corresponds to sudden increases in the mass at a given redshift—jumping from the larger value of \(S\) to the smaller value of \(S\)—it corresponds to a merger. Steady increases in the trajectory reflect accretion. Reading off the mass element’s history then proceeds by starting at the upper right of the diagram and moving to the lower left. We see that this example mass element starts out in a halo of mass \(M \approx 5\times 10^8\,M_\odot\) at \(z\approx 8\), which grows through accretion and small mergers until it has \(M\approx 2\times10^9\,M_\odot\) at \(z\approx5\). The mass element then merges with a larger halo with \(M \approx 10^{10}\,M_\odot\) and soon after at \(z\approx 4.5\) with an even larger halo with \(M\approx 10^{11}\,M_\odot\). This halo then grows through accretion and small mergers, ending up at \(M \approx 2\times 10^{12}\,M_\odot\) at \(z=0\). To get a sense of possible merger histories, it is instructive to run the code above with different random seeds to generate different trajectories.
The distinction between mergers and accretion in this picture is purely semantic, because it is set by the resolution \(\Delta S\) of our random walk: if we had used a smaller \(\Delta S\), many of the accretion episodes would actually be small mergers and in the limit \(\Delta S \rightarrow 0\) that we need to perform for a rigorous calculation of the merger history, all mass growth is through mergers. Thus, the distinction between mergers and accretion is set by the resolution \(\Delta S\), which corresponds to a relative mass resolution of \(\Delta M / M \approx \Delta S / \sqrt{S}\) in the adopted cosmological model. Below we will directly specify the mass resolution \(\Delta M\) above which we consider mass growth to be a merger and below which we label it as accretion, but the threshold has no true physical meaning.
Let us now use the picture above to mathematically determine the merger history of a galaxy. The first ingredient that we need is the probability that a mass element at \(\vec{x}\) is part of a halo with mass \(M\) that corresponds to \(S\). This is given by the derivative of Equation (19.7), which using Equation (B.20) is \begin{equation}\label{eq-mergers-eps-deltaS-eq-deltac} p[\delta(\vec{x},S) = \delta_c] = {\mathrm{d}P[\delta(\vec{x},S) \geq \delta_c] \over \mathrm{d} S} = {\delta_c \over \sqrt{2\pi\,S^3}}\,\exp\left(-{\delta_c^2\over 2S}\right)\,. \end{equation} We then want to calculate the probability that a mass element that is part of a halo with scale \(S_2\) at redshift \(z_2\) when the overdensity is \(\delta_2 = \delta_c(a_2)\) was contained in a halo with scale \(S_1 > S_2\) at redshift \(z_1 > z_2\) when the overdensity is \(\delta_1 = \delta_c(a_1)\). We can do this by recognizing that the problem of determining the probability that a trajectory wanders from \((S_2,\delta_2)\) to \((S_1,\delta_1)\) is equivalent to the original problem that we solved above of determining the probability that a trajectory wanders from \((S,\delta) = (0,0)\) to \((S,\delta_c)\), but with the origin shifted to \((S_2,\delta_2)\). Thus, for example, possible values at \(z_1=2\) when \(\delta_c(a_1) \approx 4\) of the halo scale \(S_1\) for mass elements that are part of a halo with mass set by \(S_2=6\) (\(M \approx 10^{12}\,M_\odot\)) at \(z_2=0\) are given by simulating random-walk trajectories starting from \((S,\delta) = (6,1.686)\) with an absorbing barrier at \(\delta(\vec{x},S) = \delta_c(a) = 4\), as shown in Figure 19.7.
[12]:
from numpy.random import default_rng
rng= default_rng(32323)
def random_S_walk_absorb_shifted(nS,dS=0.1,absorb=1.686,S0=0.,delta0=0.):
S= numpy.linspace(dS,(nS-1)*dS,nS-1)
out= numpy.cumsum(numpy.hstack(\
([delta0],rng.standard_normal(size=len(S))*numpy.sqrt(dS))))
# argmax stops at first match
out[numpy.argmax(numpy.hstack((out,[absorb]))>=absorb):]= absorb
return S0+numpy.hstack(([0.],S)),out
figure(figsize=(8,5))
plot(*random_S_walk_absorb_shifted(100,absorb=4.,S0=6.,delta0=1.686))
plot(*random_S_walk_absorb_shifted(100,absorb=4.,S0=6.,delta0=1.686))
plot(*random_S_walk_absorb_shifted(100,absorb=4.,S0=6.,delta0=1.686))
axhline(1.686,color='k',lw=1.5,zorder=-100)
text(0.35,1.9,r'$\delta(\boldsymbol{{x}},S)=\delta_c(z=0)$',fontsize=18.)
axhline(4.,color='k',lw=1.5,zorder=-100)
text(.35,4.25,rf'$\delta(\boldsymbol{{x}},S)=\delta_c(z={z_from_deltac(4.):.0f})$',
fontsize=18.)
xlim(0,16.4)
ylim(-5.7,5.7)
xlabel(r"$S$")
ylabel(r"$\delta_S$");
Figure 19.7: Conditional Brownian motion trajectories with an absorbing barrier.
The probability that a mass element contained in a halo of scale \(S_1\) at the redshift \(z_1\) (determined from \(\delta_c(a_1) = \delta_1\)) conditional upon the mass element being in a halo of scale \(S_2\) at the redshift \(z_2\) (set by \(\delta_c(a_2) = \delta_2\)) where \(S_1 > S_2\), \(z_1 > z_2\), and thus \(\delta_1 > \delta_2\), is then given by Equation (19.10) with the origin shifted to \((S_2,\delta_2)\) \begin{equation}\label{eq-mergers-eps-deltaS-eq-deltac-conditional} p[\delta(\vec{x},S_1) = \delta_1 | \delta(\vec{x},S_2) = \delta_2] = {(\delta_1-\delta_2) \over \sqrt{2\pi\,(S_1-S_2)^3}}\,\exp\left(-{[\delta_1-\delta_2]^2\over 2[S_1-S_2])}\right)\,. \end{equation} The fractional contribution of mass elements from halos of mass \(M_1\) to the mass \(M_2\) can be converted to a number density of halos of \(M_1\) contributing to \(M_2\) by scaling by \(M_2/M_1\) and this then gives the conditional mass function \begin{equation}\label{eq-mergers-eps-deltaS-eq-deltac-conditional-mass-function} {\mathrm{d}n(M_1\,\mathrm{at}\,z_1|M_2\,\mathrm{at}\,z_2)\over \mathrm{d}M_1} = {M_2\over M_1}\,{(\delta_1-\delta_2) \over \sqrt{2\pi\,(S_1-S_2)^3}}\,\exp\left(-{[\delta_1-\delta_2]^2\over 2[S_1-S_2])}\right)\,\left|{\mathrm{d}S_1\over \mathrm{d}M_1}\right|\,. \end{equation} We can write this conveniently using the multiplicity function \(f_{\mathrm{PS}}(\nu)\) given by Equation (17.122) used in the Press-Schechter formalism \begin{equation}\label{eq-mergers-eps-deltaS-eq-deltac-conditional-mass-function-alt} {\mathrm{d}n(M_1\,\mathrm{at}\,z_1|M_2\,\mathrm{at}\,z_2)\over \mathrm{d}\ln M_1} = {M_2\over M_1}\,f_{\mathrm{PS}}(\nu_{12})\,\left|{\mathrm{d}\ln\nu_{12}\over \mathrm{d}\ln M_1}\right|\,, \end{equation} where \(\nu_{12} = (\delta_1-\delta_2)/\sqrt{S_1-S_2}\). This expression is very similar to the Press-Shechter halo mass function in the form of Equation (17.121) (note that \(\mathrm{d}\ln\nu_{12}/ \mathrm{d}\ln M_1 = S_1/(S_1-S_2)\,\mathrm{d}\ln \sigma^{-1}_1 / \mathrm{d} \ln M_1\)).
We can, thus, easily implement it starting from the code that we used above to compute the halo mass function:
[13]:
from scipy import special
def conditional_halo_mass_function(M1,a1,M2=1e12*u.Msun,a2=1.,
S2=None,return_cumulmass=False,
window_kind='tophat'):
# Compute nu_12
Dpluszero= growth_factor(1.)
delta1= 1.686/(growth_factor(a1)/Dpluszero)
delta2= 1.686/(growth_factor(a2)/Dpluszero)
# Note, in this EPS, always evaluate S at a=1 (z=0)
S1= sigma2(M1,a=1.,window_kind=window_kind)
if S2 is None:
S2= sigma2(M2,a=1.,window_kind=window_kind)
nu12= (delta1-delta2)/numpy.sqrt(S1-S2)
# Now compute d ln nu_12 / d ln M1
# = -S1/(S1-S2) d ln sigma1 / dln M1
# and the last factor is what we used in the HMF code
dlnnu12dlnM= -S1/(S1-S2)*dlnsigmadlnM(M1,window_kind=window_kind)
# f function
def fnu(nu):
return numpy.sqrt(2./numpy.pi)*nu*numpy.exp(-nu**2./2.)
if not return_cumulmass:
return M2/M1*fnu(nu12)*dlnnu12dlnM
else:
# Also compute the cumulative mass in halos >M
return M2/M1*fnu(nu12)*dlnnu12dlnM, special.erfc(nu12/numpy.sqrt(2.))
For example, the mass function of halos at \(z=1\) that go on to form a Milky-Way-sized halo with mass \(M = 10^{12}\,M_\odot\) at \(z=0\) is displayed in Figure 19.8.
[14]:
from scipy import integrate
# Pre-calculate S2
S2= sigma2(1e12*u.Msun,a=1.,window_kind='tophat')
# Use higher resolution between 1e11 and 1e12 Msun
ms_cond_hmf= numpy.hstack((numpy.geomspace(1e5,1e11,41)[:-1]*u.Msun,
numpy.geomspace(1e11,9.5e11,21)*u.Msun))
cond_hmf_z1, cumulmass_hmf_z1= u.Quantity([conditional_halo_mass_function(\
m,0.5,S2=S2,window_kind='tophat',
return_cumulmass=True)
for m in ms_cond_hmf])\
.to(u.dimensionless_unscaled).T
cond_hmf_z4, cumulmass_hmf_z4= u.Quantity([conditional_halo_mass_function(\
m,0.2,S2=S2,window_kind='tophat',
return_cumulmass=True)
for m in ms_cond_hmf])\
.to(u.dimensionless_unscaled).T
figure(figsize=(8,5.5))
line_dndlnm_z1= loglog(ms_cond_hmf,cond_hmf_z1,
label=r'$z_1=1$')
line_dndlnm_z4= loglog(ms_cond_hmf,cond_hmf_z4,
label=r'$z_1=4$')
xlabel(r'$M/M_\odot$')
ylabel(r'$\mathrm{d} n(M_1\,\mathrm{at}\,z_1|M_2=10^{12}\,M_\odot\,\mathrm{at}\,z_2=0) / \mathrm{d}\ln M_1$')
ylim(1e-2,1e6)
legend_dndm= legend(frameon=False,fontsize=17.,loc='lower left')
gca().add_artist(legend_dndm)
ax2= gca().twinx()
line_cumul_z1= ax2.semilogx(ms_cond_hmf,cumulmass_hmf_z1,
color=line_dndlnm_z1[0].get_color(),ls='--')
line_cumul_z4= ax2.semilogx(ms_cond_hmf,cumulmass_hmf_z4,
color=line_dndlnm_z4[0].get_color(),ls='--')
ax2.set_ylabel(r'$\mathrm{Mass\ in\ halos\ with}\ M_1 > M$')
ax2.set_ylim(0.,1.)
legend_cumul= legend([Line2D([0],[0],color='k',ls='-'),
Line2D([0],[0],color='k',ls='--')],
[r'$\mathrm{d} n/\mathrm{d}\ln M_1$',r'$\mathrm{Cumulative}$'],
frameon=False,fontsize=17.,loc='upper right')
ax2.add_artist(legend_cumul);
Figure 19.8: The \(z=1\) and \(z=4\) conditional mass functions of \(z=0\) Milky-Way halos.
We see that like the halo mass function itself, the conditional mass function is dominated by small halos and the conditional mass function in that regime is almost the same as the unconditional halo mass function (essentially because \(S_1 \gg S_2\) and, thus, \(S_1-S_2 \approx S_1\) in this regime). Close to \(M_1 = M_2=10^{12}\,M_\odot\), the conditional mass function first increases starting at \(M_1 \approx 3\times10^{11}\,M_\odot\), but then quickly declines to zero as \(M_1\) approaches \(M_2\). We also show the cumulative mass in halos with \(M_1 > M\) at \(z=1\), which demonstrates that most of the mass \(M_2\) resides in massive halos \(M_1\) at \(z_1\)—the half-mass \(M_2/2\) point is at \(M_1 \approx 2\times 10^{11}\,M_\odot\). Thus, while most of the progenitors of a \(10^{12}\,M_\odot\) halo are low-mass halos that have merged together into a large halo, most of the mass \(M_2\) was already in a massive halo \(M_1\) at redshift \(z_1=1\) that grew into a slightly more massive halo at \(z_2=0\). However, if we go back to higher redshifts, a massive progenitor at that time becomes very unlikely. For example, we also display the conditional halo mass functions for an \(M_2 = 10^{12}\,M_\odot\) halo at redshift four in Figure 19.8. We see that massive \(M_1 > 10^{11}\,M_\odot\) progenitors are now highly unlikely and that many of the progenitors have \(M < 10^5\,M_\odot\).
We can similarly compute the merger rate as a function of time and halo mass or the typical formation time of a halo (defined as when it has half of its current mass) in the EPS formalism. A comparison between the predictions from the EPS formalism and the results from \(N\)-body simulations (e.g., van den Bosch 2002; Wechsler et al. 2002; Cole et al. 2008) shows that overall the agreement is quite good, but that the EPS formalism predicts a redshift evolution that is somewhat too fast, leading it to underestimate the number of massive progenitors of a given halo at high redshift as well as the formation time of halos (defined as when it has half of its current mass).
To now draw random realizations of a halo’s merger history, we start with the halo at the final time that we want to consider—often \(z=0\), but it could be a higher redshift when considering galaxies at those redshifts—and we take small steps backwards in time. At each time step, we draw a random evolutionary step backward in time for a halo of mass \(M_2\) at \(z_2\) by taking a small step \(\Delta t = t_1-t_2\) and drawing from the progenitor mass distribution given by Equation (19.12), while conserving mass (that is, the sum of the drawn progenitor masses must equal \(M_2\)). To build an entire merger history, we repeat this procedure during the next time step by applying it to all progenitors, and so on. Thus, we start from a single halo of mass \(M\) and when mergers occur at a time step, this halo splits into two or more branches. At the next step, each branch can again split into two or more branches itself. Thus, the merger history of a single halo becomes a tree of progenitor halos and their history. Because of the tree morphology of the merger history of a halo, we call it a merger tree.
However, building a merger tree using Equation (19.12) is not correct, because it does not satisfy mass conservation. This is because we should not, in fact, sample from the conditional mass function \(\mathrm{d}n(M_1\,\mathrm{at}\,z_1|M_2\,\mathrm{at}\,z_2)/ \mathrm{d}M_1\), but instead from the joint conditional mass function \(\partial n(M^1_1,M^2_1,\ldots,M^N_1\,\mathrm{at}\,z_1|M_2\,\mathrm{at}\,z_2)/ \partial M^1_1/\partial M^2_1/\ldots/\partial M^N_1\) of all possible progenitors with masses \(M^1_1,M^2_1,\ldots,M^N_1\) that satisfy \(\sum_i M^i_1 = M_2\). The EPS formalism does not allow for the calculation of this joint probability distribution. It is difficult to see how it could, given that the number of possible progenitors is unbounded as the mass of each progenitor approaches zero. For this reason, various algorithms have been developed that approximately sample from \(\partial n(M^1_1,M^2_1,\ldots,M^N_1\,\mathrm{at}\,z_1|M_2\,\mathrm{at}\,z_2)/ \partial M^1_1/\partial M^2_1/\ldots/\partial M^N_1\) by sampling from \(\mathrm{d}n(M_1\,\mathrm{at}\,z_1|M_2\,\mathrm{at}\,z_2)/ \mathrm{d}M_1\) in a way that conserves mass and approximately matches the conditional mass function of progenitors at different redshifts (e.g., Lacey & Cole 1993; Cole et al. 2000; Somerville & Kolatt 1999; Nusser & Sheth 1999). The simplest of these only consider binary mergers in a given time step. That is, they allow a halo in the tree to either remain as a single entity or to split into two progenitor halos (e.g., Lacey & Cole 1993; Cole et al. 2000). But while simple, such algorithms are not successful at matching the conditional progenitor mass function. Instead, it is necessary to allow for the possibility that a halo splits into \(N>2\) halos.
As an example of a Monte Carlo merger-tree algorithm, we implement that of Somerville & Kolatt (1999). In this algorithm, one samples masses from the conditional mass function at time \(z_1\) at each time step until the sum of the sampled masses equals \(M_2\). Sampled masses that cause the sum to exceed \(M_2\) are rejected. Like all merger-tree algorithms, this algorithm has a minimum halo mass \(M_\mathrm{res}\) that it considers, so in practice sampling proceeds until the total mass is within \(M_\mathrm{res}\) of \(M_2\); the remaining mass is considered to be part of smooth accretion (as is all mass sampled in halos with mass \(M < M_\mathrm{res}\) during the mass sampling). After sampling, we are left with a set of progenitor masses \(\{M_1^i\} > M_\mathrm{res}\). Each of these becomes the base \(M_2\) of the sampling in the subsequent step. It is typical to identify the most massive of the \(\{M_1^i\}\) as the main progenitor of the base halo, but note that this choice does not mean that the main progenitor is always the most massive progenitor at earlier times (see an example below). When taking small time steps, a typical step will involve no merger and instead lead to a single progenitor \(M_1\) that differs from \(M_2\) due to accretion of mass below \(M_\mathrm{res}\). But some steps have mergers with one or more other halos.
The sampling procedure is easiest to implement when we use \(S\) as the mass variable and \(\delta_c(a)\) as the time variable. We can then sample halos from the conditional mass function by sampling from the mass-weighted conditional mass function of Equation (19.11) written in the form \begin{equation}\label{eq-mergers-eps-deltaS-eq-deltac-conditional-in-mc-mergertree} p(S_1,\delta_1|S_2,\delta_2) = {\Delta \delta_c \over \sqrt{2\pi\,\Delta S^3}}\,\exp\left(-{[\Delta \delta_c]^2\over 2\Delta S}\right)\,, \end{equation} where \(\Delta \delta_c = \delta_1-\delta_2 = \delta_c(a_1)-\delta_c(a_2)\) and \(\delta S = S_1-S_2\). This distribution can easily be sampled by noting that under the transformation \(x = \Delta \delta_c / \sqrt{S}\), the variable \(x>0\) has a half-normal distribution \(f_\mathrm{hn}(0;1)\) (see Equation (B.21) in Appendix B.3.1), which can be sampled by sampling a normal distribution and taking the absolute value. In our case, we need to square \(x\) to obtain \(\Delta S\) and the absolute value, therefore, does not matter; we can, thus, obtain samples \(\Delta S = (\Delta \delta_c/x)^2\) where \(x\) is drawn from the standard normal distribution \(\varphi(x) = \mathcal{N}(x;0,1)\). To perform the mass-conservation accounting, we convert all \(\Delta S\) to mass using \(M(S_2+\Delta S)\). To be able to translate easily between \(S\) and \(M\) and between \(\delta_c(a)\) and \(z\), it is expedient to implement interpolated conversion functions.
[15]:
from scipy import interpolate
# M <-> S
ms_sigma= numpy.geomspace(1e5,1e16,21)*u.Msun
sigma_tophat= numpy.array([numpy.sqrt(sigma2(m,a=1.,window_kind='tophat'))
for m in ms_sigma])
_M_from_S_tophat= interpolate.InterpolatedUnivariateSpline(\
sigma_sharpk[::-1]**2.,
numpy.log(ms_sigma[::-1].to_value(u.Msun)),
k=3,ext=0)
M_from_S_tophat= lambda S: numpy.exp(_M_from_S_tophat(S))*u.Msun
_S_from_M_tophat= interpolate.InterpolatedUnivariateSpline(\
numpy.log(ms_sigma.to_value(u.Msun)),sigma_sharpk**2.,
k=3,ext=0)
S_from_M_tophat= lambda M: _S_from_M_tophat(numpy.log(M.to_value(u.Msun)))
# z <-> delta_c(a)
Dpluszero= growth_factor(1.)
as_for_growthfactor_interp= numpy.geomspace(1e-2,1.,101)[::-1]
z_from_deltac= interpolate.InterpolatedUnivariateSpline(\
[1.686/(growth_factor(a)/Dpluszero)
for a in as_for_growthfactor_interp],
1./as_for_growthfactor_interp-1.,
k=3,ext=0)
deltac_from_z= interpolate.InterpolatedUnivariateSpline(\
1./as_for_growthfactor_interp-1.,
[1.686/(growth_factor(a)/Dpluszero)
for a in as_for_growthfactor_interp],
k=3,ext=0)
We then implement the merger-tree algorithm to evolve a halo with mass \(M_0\) backwards in time from redshift \(z_0\) to redshift \(z_\mathrm{end}\). Because internally we use \(\delta_c(a)\) as the time variable, we generate output on a grid in \(\delta_c(a)\) starting with \(\delta_c(a_0)\) and using a user-defined step \(\Delta \delta_c\); if \(\delta_c(a_\mathrm{end})\) is not an exact integer factor of \(\Delta \delta_c\) larger than \(\delta_c(a_0)\), the end point is not exactly \(z_\mathrm{end}\) (but we can force this to be the case by setting \(\Delta \delta_c\) judiciously). While we work in our standard \(\Lambda\)CDM cosmology, where the growth factor is given by Equation (17.67) and the dependence of \(\delta_c(a)\) on \(z\) is somewhat complex, in fact the growth factor in our \(\Lambda\)CDM cosmology is not that different from the \(D_+(a) = a\) behavior for a flat, matter-only (Einstein–de-Sitter) Universe. Therefore, we have that \(\delta_c(a) \approx \delta_c\,(1+z)\) and \(\Delta \delta_c \approx 1.686\,\Delta z\).
We implement the merger-tree algorithm in the class below. Initialization of the tree only stores the basic parameters of the tree (starting mass and redshift \(M_0\) and \(z_0\), final redshift \(z_{\mathrm{end}}\), mass resolution \(M_{\mathrm{res}}\), and the step size \(\Delta \delta_c\)) and some derived quantities. Random tree realizations are then generated using the construct method, which implements the merger-tree algorithm. While tree algorithms are naturally implemented using recursion, languages such as Python do not handle recursion efficiently, so instead we implement the algorithm by keeping a stack of halos to be processed and then working through the stack. For each halo in the stack, we follow the main progenitor’s mass accretion history and we push any mergers with lower-mass halos onto the stack to be processed similarly. The tree is fully constructed when the stack is empty. To record all of the information about the merger tree in a way that allows us to easily investigate the mass history of individual halos in the tree and the overall distribution of progenitor masses at earlier times, we can store the mass history of all halos in the tree in a single large matrix with size \(N_\mathrm{halo} \times N_{\delta_c}\) (setting all values to a negative number at all times after a halo has merged with another halo). We also keep track of the entire merger hierarchy with a simple string-based scheme (using, e.g., "1-0" to indicate the first descendant of the main progenitor, "3-2-0" to denote the third descendant of the main progenitor's second descendant, etc.). Because the mass history of the main progenitor is often useful on its own (without caring about how it obtained its mass through either mergers or accretion), we can also implement an option to only follow the main progenitor by not pushing any mergers onto the stack in this case.
[16]:
class mergertree:
def __init__(self,M0,z0=0.,zend=3.,Mres=1e9*u.Msun,ddeltac=0.05):
self._M0= M0
self._S0= S_from_M_tophat(self._M0)
# 'time' grid
self._z0= z0
self._deltac0= deltac_from_z(self._z0)
self._zend= zend
self._deltacend= deltac_from_z(self._zend)
self._ddeltac= ddeltac
self.deltacs= numpy.arange(self._deltac0,self._deltacend,self._ddeltac)
self._ndeltacs= len(self.deltacs)
self._Mres= Mres
def construct(self,main=False,random=numpy.random.RandomState(seed=None)):
# Stack contents: M,S,index in deltacs,branch indicator
stack= [(self._M0,self._S0,0,"0")]
out= []
out_branch= []
while len(stack) > 0:
M,S,deltac_ii,branch= stack[0]
# Prepare output
mass_history= numpy.zeros(self._ndeltacs)-1.
mass_history[deltac_ii]= (M/self._M0).to(u.dimensionless_unscaled)
out_branch.append(branch)
ndesc= 0 # To keep track of how many descendants we are creating
for ii in range(deltac_ii+1,self._ndeltacs):
# Sample from the progenitor distribution
newM, newS= [], []
cumul_mass= 0.
while True:
this_newS= S+(self._ddeltac/random.normal(size=1))**2.
this_newM= M_from_S_tophat(this_newS)
if cumul_mass+this_newM > M:
continue
cumul_mass+= this_newM
if cumul_mass > M-self._Mres \
or (main and this_newM/M > 0.5): # can stop here
newM.append(this_newM)
newS.append(this_newS)
break
if this_newM < self._Mres:
continue
newM.append(this_newM)
newS.append(this_newS)
newM= u.Quantity(newM).flatten()
newS= numpy.array(newS).flatten()
# Stop following halos that fall below the mass limit
if numpy.all(newM < self._Mres):
mass_history[ii:]= \
(newM[0]/self._M0).to(u.dimensionless_unscaled)
break
# Sort by mass
sindx= numpy.argsort(newM)[::-1]
newM, newS= newM[sindx], newS[sindx]
# Most massive one is the one we follow here
mass_history[ii]= (newM[0]/self._M0).to(u.dimensionless_unscaled)
M= newM[0]
S= newS[0]
# Push the rest onto the stack if we're not just
# following the main progenitor
if not main and len(newM) > 1:
stack.extend([(tM,tS,ii,f"{ndesc+tl}-"+branch)
for tM,tS,tl in zip(newM[1:],newS[1:],
range(1,len(newM)))])
ndesc+= len(newM)-1
out.append(mass_history)
# Remove the current one from the stack
stack.pop(0)
self.mass_history= out*self._M0
self.branch= out_branch
As a first example of using the merger trees, let us look at the typical dark-matter mass accretion history of a Milky-Way-sized halo with \(M_0 = 10^{12}\,M_\odot\) in Figure 19.9. We generate a large number of merger trees to look at the mean accretion history while also showing a few examples of individual trajectories. To avoid resolution effects, we only consider masses at least ten times larger than the resolution (which we set to \(M_\mathrm{res} = 10^8\,M_\odot\)).
[17]:
from matplotlib.ticker import FuncFormatter
from numpy.random import default_rng
rng= default_rng(11111) # For reproducibility
Mres= 1e8*u.Msun
mtree= mergertree(1e12*u.Msun,Mres=Mres,ddeltac=0.05,zend=12)
ntrees= 150
figure(figsize=(6,4))
for ii in range(ntrees):
mtree.construct(main=True,random=rng)
if ii == 0:
histories= mtree.mass_history[0].copy()/ntrees
else:
histories+= mtree.mass_history[0].copy()/ntrees
if ii < 4:
loglog(1+z_from_deltac(mtree.deltacs),mtree.mass_history[0])
loglog(1+z_from_deltac(mtree.deltacs),histories,
color='k',lw=2.,zorder=100)
xlabel(r'$1+z$')
ylabel(r'$M/M_\odot$')
ylim(10*Mres.to_value(u.Msun),2e12)
xticks(numpy.arange(1,11))
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)));
Figure 19.9: The mass accretion history of Milky-Way halos.
We see that there is a well-defined overall trend, with Milky-Way-sized halos typically having half their present-day mass at \(z\approx 1\) and a tenth of their mass at \(z\approx 3.5\). But the individual trajectories demonstrate that there is a lot of variance around this mean trend, with some halos assembling most of their mass much earlier and some much later than average. When normalized by the mass \(M_0\) of the target halo, these mass accretion histories are quite universal, with only minor differences as a function of \(M_0\) and the cosmological parameters (van den Bosch 2002).
Next, we can check whether the trees’ progenitor mass distributions at earlier times are consistent with that from Equation (19.12). For a single time step, the progenitor mass distribution of the tree matches Equation (19.12) by construction, but going back multiple time steps, this does not have to be the case. Thus, we can look at the progenitor mass distribution of a \(M_0 = 10^{12}\,M_\odot\) halo at \(z=1\) in Figure 19.10 and compare it to what we found in Figure 19.8 above; to obtain sufficient numbers of progenitor halos, we again use a large number of trees.
[18]:
Mres= 1e8*u.Msun
# Set ddeltac such that we end up at zend=1 exactly
mtree= mergertree(1e12*u.Msun,Mres=Mres,ddeltac=0.0494,zend=1)
ntrees= 30
figure(figsize=(8,5.5))
bins= numpy.geomspace(Mres.to_value(u.Msun),1e12,31)
for ii in range(ntrees):
mtree.construct()
if ii == 0:
dndlog10m,_= numpy.histogram(mtree.mass_history[:,-1].to_value(u.Msun),
bins=bins)
else:
dndlog10m+= numpy.histogram(mtree.mass_history[:,-1].to_value(u.Msun),
bins=bins)[0]
dndlog10m= dndlog10m.astype('float')
dndlog10m/= numpy.sum(dndlog10m*(numpy.roll(bins,-1)-bins)[:-1])
stairs(dndlog10m,bins,lw=1.5)
# Also plot the EPS prediction that we computed above
line_dndlnm_indx= ms_cond_hmf > Mres
line_dndlnm= loglog(ms_cond_hmf[line_dndlnm_indx],
cond_hmf_z1[line_dndlnm_indx]\
/numpy.sum(cond_hmf_z1[line_dndlnm_indx]\
*(ms_cond_hmf-numpy.roll(ms_cond_hmf,1))[line_dndlnm_indx]))
xlabel(r'$M/M_\odot$')
ylabel(r'$\mathrm{d} n(M_1\,\mathrm{at}\,z_1|M_2=10^{12}\,M_\odot\,\mathrm{at}\,z_2=0) / \mathrm{d}\ln M_1$');
Figure 19.10: The \(z=1\) conditional mass function of \(z=0\) Milky-Way halos from merger trees.
We see that they match quite well. This is not the case in many alternative merger-tree algorithms, especially those that employ a binary structure (e.g., Lacey & Cole 1993; Cole et al. 2000).
Finally, we plot a representation of the entire merger tree, generating an image consisting of bands with widths proportional to the mass of each halo in the tree. This would again be naturally implemented using recursion (because each new halo is the base of its own tree, so its merger-tree image is a sub-image of the overall one), so we can employ a stack-based algorithm like for the merger-tree itself.
[19]:
def plot_mergertree(mtree,resolution_factor=10,
random=numpy.random.RandomState(seed=None),
**plot_kwargs):
nres_elements= int(numpy.ceil(resolution_factor\
*(mtree._M0/mtree._Mres)\
.to(u.dimensionless_unscaled)))
img= numpy.zeros((nres_elements,mtree._ndeltacs))+numpy.nan
# give each halo a different color, randomize to increase contrast
base_colors= numpy.hstack(([0],
random.permutation(numpy.arange(1,len(mtree.branch)))))
stack= [(mtree.mass_history,mtree.branch,0,nres_elements,0,base_colors)]
while len(stack) > 0:
mass_history,branch,start_indx,end_indx,deltac_indx,colors= stack[0]
# Plot the current halo's mass history
for ii in range(deltac_indx,mtree._ndeltacs):
img[start_indx:start_indx+1+int(numpy.floor(resolution_factor\
*mass_history[0,ii]/mtree._Mres)),ii]\
= colors[0]
# split into the next level and add all of those to the stack
# current one is the base so has zero -
levels= numpy.array([b.count('-') for b in branch],dtype='int')
todo_indx= levels != 0
new_mass_history= mass_history[todo_indx]
new_colors= colors[todo_indx]
#remove previous level
new_branch= [b[:-2] for jj,b in enumerate(branch) if todo_indx[jj]]
levels= levels[todo_indx]
# find new trunks
next_indices= [jj for jj,level in enumerate(levels) if level == 1]
# New halos are put on the right of the present base, if there are
# multiple, we need to shift them a bit to the left by
# the mass of the previous
mass_already_accounted_for= \
u.Quantity(numpy.zeros(mtree._ndeltacs)*u.Msun)
for next_indx in next_indices:
next_branch_trunk= new_branch[next_indx]
new_indices= [next_indx]
new_indices.extend([jj for jj in range(len(new_branch))
if new_branch[jj][-2:] == f"-{next_branch_trunk}"])
# Figure out where this merger occurred
merger_indx= list(new_mass_history[new_indices][0] > 0).index(True)
new_end_indx= start_indx\
+int(numpy.floor(resolution_factor\
*(mass_history[0,merger_indx-1]
-mass_already_accounted_for[merger_indx-1])/mtree._Mres))
new_start_indx= new_end_indx\
-int(numpy.floor(resolution_factor\
*new_mass_history[next_indx,merger_indx]/mtree._Mres))
stack.append((new_mass_history[new_indices],
[new_branch[jj] for jj in new_indices],
new_start_indx,new_end_indx,
merger_indx,
new_colors[new_indices]))
mass_already_accounted_for[merger_indx-1]+= \
new_mass_history[next_indx,merger_indx]
# Remove the current one from the stack
stack.pop(0)
# Now plot
plot_kwargs['interpolation']= plot_kwargs.get('interpolation','nearest')
plot_kwargs['cmap']= plot_kwargs.get('cmap','cividis')
galpy_plot.dens2d(img.T,origin='upper',
yrange=[mtree.deltacs[0],mtree.deltacs[-1]],
xrange=[0.,1.],
**plot_kwargs)
yticks([mtree.deltacs[-1]+mtree.deltacs[0]-deltac_from_z(pz)
for pz in numpy.arange(mtree._z0+1e-3,mtree._zend+1)])#1e-3 to avoid -0
def z_label(x,pos):
return rf"${z_from_deltac(mtree.deltacs[-1]+mtree.deltacs[0]-x):.1f}$"
gca().yaxis.set_major_formatter(FuncFormatter(z_label))
xlabel(r"$M/M(z=0)$")
ylabel(r"$\mathrm{redshift}$")
return img
Two example merger trees for a \(M_0 = 10^{12}\,M_\odot\) halo out to \(z=5\) going down to \(M_\mathrm{res} = 3\times 10^{10}\,M_\odot\) are shown in Figure 19.11.
[20]:
mtree= mergertree(1e12*u.Msun,Mres=3e10*u.Msun,
ddeltac=0.0494,zend=5)
from numpy.random import default_rng
figure(figsize=(8.5,4))
for ii,seed in enumerate([11111,15151]):
mtree.construct(main=False,random=default_rng(seed))
subplot(1,2,ii+1)
img= plot_mergertree(mtree,random=default_rng(12121),gcf=True)
tight_layout(w_pad=-2.3);
Figure 19.11: Example merger trees for a \(M=10^{12}\,M_\odot\) halo.
These histories illustrate a few important aspects of merger trees. We see that multiple mergers in a given time step do occur quite commonly, which is why binary-tree algorithms are not as successful in reproducing the analytical results from the EPS formalism as algorithms that allow for \(N>2\) progenitors at each time step. We also notice that the main progenitor (the leftmost band) is not always the most massive progenitor at each time. These merger trees give a good sense of the busy merging histories of the dark-matter halos that host large galaxies. Generate some more to explore the variety of merger histories of galaxies!
Merger trees allow for the quick and easy calculation of many properties of the merger and accretion history of dark matter halos. Some of these can be calculated analytically in the EPS formalism (e.g., the merger rate as a function of mass and redshift), but a random sample of merger trees allows us to compute such quantities as the probability of a major merger for a halo of mass \(M\) during a given time interval (say, since redshift one), the ratio of the incidence of minor and major mergers vs. redshift, and many others.
.
.
.