17.1. Cosmological evolution of small perturbations¶
Dark matter halos—and the galaxies that reside within them—form from small seeds sown by whatever process sets the initial pattern of overdensities and underdensities in all of the components that make up the Universe. The overdensities start out as very small perturbations on top of the smooth density. They evolve under the influence of gravity and other physical processes and eventually grow so large that they form the bound, virialized objects known as dark matter halos. When the perturbations are small, their evolution can be studied using relatively simple tools, which we describe in this first section of this chapter.
17.1.1. Homogeneous expansion of the Universe¶
We start with a brief discussion of the homogeneous expansion of the Universe that provides the backdrop for the formation of galaxies. Our discussion will be brief and focused on the elements we need to understand the formation of galaxies; we refer to standard undergraduate and graduate cosmology texts for a more in-depth and expanded discussion. A formal derivation of the material in this sub-section using the general theory of relativity is presented in Appendix C.4, but is not required to understand the material in this or subsequent chapters.
In our current cosmological framework, known as \(\Lambda\)CDM, the Universe is assumed to be homogeneous and isotropic and its homogeneous evolution is described by the Friedmann–Lemaître–Robertson–Walker (FLRW) metric. “Homogeneous” here means that all constituents of the Universe have equal energy densities everywhere in space, such that the only changes in them are a function of time. The FLRW metric depends only on the scale factor \(a\)—the relative size of the Universe relative to today—and a parameter \(k \in \{-1,0,+1\}\) that indicates whether the spatial curvature of the Universe is closed, flat, or open. The Einstein equations from the general theory of relativity relate the metric to the energy components present in the Universe, which allows us to derive evolution equations for the scale factor; these are the Friedmann equations. Cosmological observations over the last few decades have determined the energy components present in our Universe to high precision. For our cosmological purposes, the energy density of the Universe consists of three main constituents: matter, dark energy, and radiation.
Matter contains both dark matter and the ordinary matter of the standard model that constitutes stars, gas, planets, and all of the matter of our life on Earth, which is typically referred to as “baryons”. Baryons only make up about 15% of the total matter content of the Universe. While, as we have seen in the previous parts, the behavior of baryons within galaxies is very different from that of dark matter, on larger, cosmological scales, baryons and dark matter act almost in concert—we discuss this in more detail below—and for the homogeneous evolution we can combine them into a single “matter” component.
Dark energy causes the acceleration of the expansion observed at \(z < 1\) (Riess et al. 1998; Perlmutter et al. 1999) and all current data is consistent with dark energy taking the form of a cosmological constant. For our purposes, this is a form of energy that does not dilute under the expansion of the Universe and that does not cluster. Thus, it mainly affects the homogeneous evolution of the Universe.
Radiation largely consists of photons, but also contains relativistic neutrinos. While there are currently only upper limits on the masses of neutrinos and they may therefore be relativistic at all times, standard cosmological models since the first release of the Planck data (Planck Collaboration et al. 2014) assume a sum of neutrino masses of \(0.06\,\mathrm{eV}\). With this non-zero mass, neutrinos are non-relativistic today, but relativistic at the epoch of matter-radiation equality (see below). This complicates cosmological calculations significantly, but the small neutrino masses mean that the difference with assuming that they are massless are small. We will simply assume that neutrinos are massless and can then lump photons and neutrinos together as “radiation”.
Because on the whole no matter or radiation is created or destroyed, the energy density of matter dilutes simply as one over volume or \(\rho_m \propto 1/a^3 = (1+z)^3\), where \(z\) is the cosmological redshift \(z = 1/a - 1\). The energy density of radiation has an additional contribution of \(1/a\) because photons redshift in an expanding Universe, so radiation dilutes as \(\rho_r \propto 1/a^4 = (1+z)^4\). With the dark energy density remaining constant (\(\rho_\mathrm{DE} \propto \mathrm{constant}\)), the Universe therefore goes through three main epochs.
Radiation domination, which starts soon after the Big Bang and lasts until the epoch of matter-radiation equality where, as the name suggests, the energy density of radiation matches that of matter. In our Universe, this happens at a redshift of \(\approx 3400\) or about 50,000 years after the Big Bang (Planck Collaboration et al. 2020b).
Matter domination, which lasts from matter-radiation equality until \(z \approx 0.3\), when the energy density of matter matches that of dark energy.
Dark-energy domination, which describes the present-day Universe and likely the future.
The time evolution of the scale factor is determined using the Friedmann equations (Equations C.130 and C.131 in Appendix C.4), which we can write as \begin{align}\label{eq-dmhaloform-friedmann-1} H^2 & = \left({\dot{a} \over a}\right)^2= {8\pi G\rho\over 3}-{c^2k \over a^2R_0^2}\,,\\ \dot{H} + H^2 & = {\ddot{a} \over a} = -{4\pi G\over 3}\left(\rho + 3{p\over c^2}\right)\,,\label{eq-dmhaloform-friedmann-2} \end{align} where \(H \equiv \dot{a}/a\) is the Hubble parameter, \(\rho\) and \(p\) are the sum of the energy densities and pressures of all of the constituents of the Universe, respectively, and unlike in Equations (C.130) and (C.131) we have decided to include dark energy as part of \(\rho\) and \(p\) rather than through the cosmological constant \(\Lambda\). By solving the Friedmann equations, we can determine the dependence of the scale factor on time for a given set of cosmological parameters. Because of this, we can express the time coordinate equivalently as time \(t\), scale factor \(a\), or redshift \(z\). In this chapter, and in cosmology more generally, we will use time, scale factor, and redshift interchangeably to act as the time coordinate.
The relation between pressure and density for each component is typically assumed to follow an equation of state of the form \(p = w \rho c^2\), where \(w = 0\) for matter, \(w=1/3\) for radiation, and \(w = -1\) for dark energy. Equation (17.1) makes it clear that when \(\rho = \rho_{c}\) with \begin{equation}\label{eq-dmhaloform-rhocrit} \rho_c = {3H^2\over 8\pi G} \end{equation} we have that \(k=0\) and the Universe is spatially flat. The density \(\rho_c\) is the critical density. Using the critical density, we can re-write the first Friedmann equation (17.1) in terms of density parameters \(\Omega_i = \rho_i / \rho_c\) as \begin{align}\label{eq-dmhaloform-friedmann-matter-de-radiation-curvature} {H(a)^2 \over H_0^2} & = \Omega_{0,m}\,a^{-3} + \Omega_{0,\Lambda}+\Omega_{0,r}\,a^{-4}+\Omega_{0,k}\,a^{-2}\,, \end{align} where the ‘\(0\)’ subscript indicates that quantities are evaluated at the present day, \(z=0\) (below, we will require the density parameters at earlier times and we will denote this often as, e.g., \(\Omega_m(a)\)). In this equation, we have also represented the curvature term using a density parameter \(\Omega_{0,k} = -c^2k/(H^2_0 R_0^2)\). The density parameters of the various energy components can be constrained using observations of the cosmic microwave background and of the large-scale structure of galaxies. These observations constrain the spatial curvature to be small \(\Omega_{0,k} = 0.001\pm0.002\) (Planck Collaboration et al. 2020b) and the standard model of cosmology therefore assumes that \(\Omega_{0,k} =0\) (there are good theoretical reason for \(\Omega_{0,k}\) to be zero or at least orders of magnitudes smaller than the observational bound as well). Assuming flatness, the latest measurements from Planck Collaboration et al. (2020b) are \begin{align}\label{eq-dmhaloform-planck2018-Om0} \Omega_{0,m} & = 0.3111 \pm 0.0056\,,\\ \Omega_{0,\Lambda} & = 0.6889 \pm 0.0056\,,\\ \Omega_{0,r} & = (9\pm0.4)\times 10^{-5}\,. \end{align} Note that because \(\Omega_{0,m} +\Omega_{0,\Lambda}+\Omega_{0,r} = 1\), we do not need to specify all three, but it is useful to list all three because the uncertainties on the matter and dark-energy components are larger than the radiation parameter. Combining these density parameters with a measurement of the Hubble parameter at redshift zero \(H_0\), the Hubble constant, we have fully specified the homogeneous expansion of the Universe. The Planck Collaboration et al. (2020b) measurements give \(H_0 = 67.7\,\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1}\). The Hubble constant is also often expressed using the dimensionless parameter \(h\) defined through \(H_0 = 100\,h\,\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1}\) and known affectionately as “little \(h\)”. For the Planck Collaboration et al. (2020b) value we, therefore, have that \(h = 0.677\). As discussed above, the Planck Collaboration et al. (2020b) results actually assume a small non-zero neutrino mass, which contributes \(\approx 0.014\) to \(\Omega_{0,m}\), but this is a subtlety that we ignore here. Because we will consider the effect of baryons below as well, we further note that the Planck Collaboration et al. (2020b) measurement of the baryon density parameter is \begin{equation} \Omega_{0,b} = 0.0490\pm 0.0007\,. \end{equation}
To obtain the age of the Universe as a function of scale factor, we can solve Equation (17.4) through direct integration. This is implemented in the following function and we also compute the current age of the Universe:
[6]:
Om0= 0.3111
Ode0= 0.6889
Or0= 9e-5
H0= 67.70*u.km/u.s/u.Mpc
from scipy import integrate
def age(a):
return (integrate.quad(lambda ap: 1./ap/numpy.sqrt(Om0/ap**3.+Ode0+Or0/ap**4.),
0.,a)[0]/H0).to(u.Gyr)
print(f"The current age of the Universe is {age(1.):.2f}")
The current age of the Universe is 13.78 Gyr
The age of the Universe as a function of redshift is shown in Figure 17.1.
[7]:
a= numpy.geomspace(1e-4,1.,101)
figure(figsize=(7,5))
galpy_plot.plot(1./a,[age(ap).to_value(u.Myr) for ap in a],
loglog=True,
xlabel=r'$1+z$',
ylabel=r'$\mathrm{Age\,(Myr)}$',
xrange=[40000.,0.9],yrange=[3e-3,95000],
gcf=True,
label=r'$\mathrm{Flat}\ \Lambda\mathrm{CDM}$')
from astropy.cosmology import Planck18 as cosmo
loglog(1./a,[cosmo.age(1./ap-1).to_value(u.Myr) for ap in a],
label=r'$\mathrm{Planck\ 2018}$')
from matplotlib.ticker import FuncFormatter
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
gca().yaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
# Mark important epochs
# Matter-radiation equality
zeq= 3387.
vlines(1.+zeq,ymin=3e-3,ymax=age(1./(1.+zeq)).to_value(u.Myr),
color='0.5',lw=3.)
hlines(age(1./(1.+zeq)).to_value(u.Myr),xmin=40000.,xmax=1.+zeq,
color='0.5',lw=3.)
text(33000,0.01,r'$\mathrm{Matter-}$'+'\n'+r'$\mathrm{radiation}$'
+'\n'+r'$\mathrm{equality}$',fontsize=14.)
# Decoupling (CMB)
zdec= 1090.
vlines(1.+zdec,ymin=3e-3,ymax=age(1./(1.+zdec)).to_value(u.Myr),
color='0.5',lw=3.)
hlines(age(1./(1.+zdec)).to_value(u.Myr),xmin=40000.,xmax=1.+zdec,
color='0.5',lw=3.)
text(33000,0.49,r'$\mathrm{Recombination}$',fontsize=14.)
# Reionization
zre= 7.82
vlines(1.+zre,ymin=3e-3,ymax=age(1./(1.+zre)).to_value(u.Myr),
color='0.5',lw=3.)
hlines(age(1./(1.+zre)).to_value(u.Myr),xmin=40000.,xmax=1.+zre,
color='0.5',lw=3.)
text(33000,900,r'$\mathrm{Reionization}$',fontsize=14.)
# Reionization
zeqde= 0.3054
vlines(1.+zeqde,ymin=3e-3,ymax=age(1./(1.+zeqde)).to_value(u.Myr),
color='0.5',lw=3.)
hlines(age(1./(1.+zeqde)).to_value(u.Myr),xmin=40000.,xmax=1.+zeqde,
color='0.5',lw=3.)
text(33000,13500,r'$\mathrm{Matter-dark\,energy\ equality}$',fontsize=14.)
# legend
legend(frameon=True,fontsize=18.,loc='lower right',
facecolor='w',edgecolor='w');
Figure 17.1: Important cosmological events over the age of the Universe.
In this figure, we compare our solution to the calculation using astropy’s cosmology module, which correctly accounts for the non-zero neutrino masses in the Planck Collaboration et al. (2020b) cosmology. They overlap entirely in the figure and in detail agree to better than a percent over the entire redshift range shown. We have also labeled a few important cosmological epochs: the epochs of matter-radiation and matter-dark energy equality already discussed above, the recombination era when photons and
baryons decouple and the cosmic microwave background is formed (\(z \approx 1100\) at \(\approx 372,000\) years), and the redshift at which the Universe is reionized after a period where baryons are largely neutral hydrogen (\(z \approx 10\) at \(\approx 1\,\mathrm{Gyr}\)).
For the cosmological mix that represents our Universe, we cannot analytically solve for the time dependence of the scale factor \(a(t)\). But galaxies in the Universe largely form between \(z \approx 10\) and \(z \approx 1\) and during this period the energy density of the (spatially flat) Universe is dominated by matter, \(\Omega_m = 1\). This is the Einstein–de Sitter model and for this case we can solve the cosmological evolution of the scale factor and Hubble parameter analytically as \begin{align}\label{eq-dmhaloform-eds-aH} a(t) & = a_0\,\left(\frac{t}{t_0}\right)^{\frac{2}{3}}\,;\quad H(t) = \frac{2}{3}\,t^{-1}\,, \end{align} where \(t_0\) is a reference time (e.g., the present time, in which case \(a_0 = 1\)). This solution will be useful in our discussion of the growth of structure in the matter-dominated era below.
Finally, because the overall expansion of the Universe causes distances to increase in proportion to the scale factor, it is often advantageous to factor this effect out and work in comoving coordinates, which are coordinates \(\vec{x}\) that are relative to the expanding Universe and are defined by \(\vec{r} = a(t)\,\vec{x}\), where \(a(t)\) is the scale factor at time \(t\) and \(\vec{r}\) are physical coordinates. Distances measured in comving coordinates are comoving distances and, similarly, we can define comoving versions of quantities such as radii, volumes, and number and mass densities.
Now that we understand the basics of the evolution of the homogeneous Universe, it is time to turn our attention to the evolution of inhomogeneities.
17.1.2. Evolution of small cold-dark-matter perturbations¶
The dark-matter halos that host galaxies form through the gravitational growth and collapse of initial, small density perturbations. These density perturbations grow at different rates during the epochs of radiation, matter, and dark-energy domination discussed above. The initial density perturbations are created at very early times, with the leading theory being that the Universe experienced a period of rapid accelerated expansion known as inflation (Guth 1981; Linde 1982; Albrecht & Steinhardt 1982) during which quantum fluctuations seeded an initial spectrum of density fluctuations (Mukhanov & Chibisov 1981; Hawking 1982; Starobinsky 1982; Guth & Pi 1982; Bardeen et al. 1983). There is currently no direct evidence that inflation occurred, although all cosmological data is consistent with the predictions from inflation, some of which are very specific. However, for our purposes all that matters is that some process laid down an initial pattern of small density fluctuations that then started to evolve in a Universe consisting of radiation, matter, and dark energy.
Calculating the evolution of the initial density perturbations from when they were first created until the epoch of galaxy formation is complicated by the fact that (i) all relevant scales were larger than the particle horizon (or comoving horizon)—the largest distance that light can travel since the Big Bang—at some point, because the horizon grows with time after inflation, and on such super-horizon scales one needs the full general theory of relativity to describe the evolution of perturbations; and (ii) the gravitational potential is dominated by photons during radiation domination and the evolution of photon perturbations is complex. We ignore these issues for now and simply assume that we can evolve the density perturbations until a redshift \(z < 1,000\) deep in the matter-domination era, at which point density perturbations are still very small and the relevant scales are well inside the horizon, so we can use the Newtonian approximation to the general theory of relativity. This allows us to focus on their subsequent growth and collapse into dark-matter halos. In Section 17.1.5 below, we will return to the evolution of the density perturbations in the early Universe to fully characterize the statistical properties of density perturbations during the epoch of galaxy formation and evolution.
In the standard model of cosmology, dark matter is collisionless and its evolution on scales small compared to the horizon is therefore given by the collisionless Boltzmann equation, which we derived in Chapter 5.3. The collisionless Boltzmann equation gives the evolution of the distribution function \(f(\vec{q},\vec{p},t)\) that represents the number of dark matter particles in a given volume of phase space centered on \((\vec{q},\vec{p})\) at time \(t\). In proper coordinates \((\vec{r},\vec{u})\)—actual coordinates, not those expanding with the Universe—the Lagrangian is the usual one \begin{equation} \mathcal{L}(\vec{r},\dot{\vec{r}},t) = \frac{m|\dot{\vec{r}}^2|}{2}-m\Psi(\vec{r},t)\,, \end{equation} where \(\Psi(\vec{r},t)\) is the gravitational potential, because the dark matter only interacts gravitationally (we use \(\Psi\) rather than \(\Phi\), because we will later introduce \(\Phi\) as the potential due to perturbations in the density, not the full density).
Because the Universe is expanding isotropically and homogeneously, it is expedient to express equations in terms of comoving coordinates \(\vec{r} = a(t)\,\vec{x}\), where \(a(t)\) is the scale factor at time \(t\). The Lagrangian can then be written in terms of \((\vec{x},\dot{\vec{x}})\) as \begin{equation} \mathcal{L}(\vec{x},\dot{\vec{x}},t) = \frac{m}{2}\,\left(\dot{a}\,\vec{x} + a\,\dot{\vec{x}}\right)^2-m\Psi(\vec{x},t)\,. \end{equation} Using Equation (3.26) from Chapter 3.3, we obtain the canonical momentum conjugate to \(\vec{x}' = \vec{x}\)—we denote this canonical coordinate pair with an apostrophe, because we will perform a further transformation to the final set of canonical coordinates): \begin{equation} \vec{p}' = \frac{\partial \mathcal{L}}{\partial \dot{\vec{x}'}} = m\,a\,\left(\dot{a}\,\vec{x}' + a\,\dot{\vec{x}}'\right)\,. \end{equation} The Hamiltonian in terms of \((\vec{x}',\vec{p}')\) is given by \begin{align} H(\vec{x}',\vec{p}',t) & = \frac{m\,a^2}{2}\,\dot{\vec{x}}'^2 -\frac{m}{2}\,\dot{a}^2\,\vec{x}'^2+m\Psi(\vec{x}',t) \end{align} To further simplify this, we canonically transform \((\vec{x}',\vec{p}')\) using a generating function of the third kind \(F_3(\vec{p}',\vec{x})\) that is a function of the old momentum and the new coordinate. Specifically, we use \(F_3(\vec{p}',\vec{x}) = -\vec{p}'\cdot\,\vec{x} + m\,a\,\dot{a}\,|\vec{x}|^2/2\), because then the transformation equations give \begin{align} \vec{p} & = -\frac{\partial F_3}{\partial \vec{x}} = m\,a^2\,\dot{\vec{x}}\,;\quad \vec{x}' = -\frac{\partial F_3}{\partial \vec{p}'} = \vec{x}\,. \end{align} Thus, the effect of the canonical transformation is to remove the term \(\propto \vec{x}\) in the canonical momentum, while keeping the coordinate the same. The momentum is now proportional to \(\vec{v} = a\,\dot{\vec{x}}\), the peculiar velocity with respect to the expansion of the Universe. The Hamiltonian is transformed under the canonical transformation specified above as \begin{align} H \rightarrow & H + \frac{\partial F_3}{\partial t}= \frac{1}{2\,m\,a^2}\,\vec{p}^2+m\Phi(\vec{x},t)\,, \end{align} where we have defined \begin{equation}\label{eq-dmhaloform-potential-perturbation} \Phi(\vec{x},t) = \Psi(\vec{x},t) + \frac{a\,\ddot{a}}{2}\,|\vec{x}|^2\,. \end{equation}
To determine the density that sources the potential \(\Phi\), we can compute the Laplacian \(\nabla^2 \Phi(\vec{x},t)\), where the derivative is with respect to \(\vec{x}\). Working for the moment in a matter-only Universe as appropriate during the matter-dominated era of galaxy formation, we have that \begin{align} \nabla^2 \Phi(\vec{x},t) & = \nabla^2 \Psi(\vec{x},t) + 3\,a\,\ddot{a}= 4\pi\,G\,\overline{\rho}(t)\,a(t)^2\,\delta(\vec{x},t)\,,\label{eq-poisson-comoving-perturb} \end{align} where we have used the second Friedmann equation (17.2) for a Universe consisting solely of matter (\(p = \Lambda = 0\)). The potential \(\Phi\) is therefore that due to the density perturbation \(\delta(\vec{x},t) = \rho(\vec{x},t)\,/\,\overline{\rho}(t)-1\) with respect to the mean background matter density. This continues to hold when dark energy becomes a significant constituent of the energy density (because then \(\nabla^2 \Psi(\vec{x},t) = 4\pi\,G\rho(\vec{x},t)\,a(t)^2 - \Lambda c^2\,a(t)^2\); see Equation C.69).
We can now write the collisionless Boltzmann equation in the comoving, canonical coordinates by using the Hamiltonian in the expression of the collisionless Boltzmann equation for general canonical coordinates from Equation (5.31): \begin{equation}\label{eq-cbe-comoving} \frac{\partial f}{\partial t} +\frac{1}{m\,a^2}\,\vec{p}\,\frac{\partial f}{\partial \vec{x}} -m\,\frac{\partial \Phi}{\partial \vec{x}}\,\frac{\partial f}{\partial \vec{p}}=0\,. \end{equation} We can then again derive Jeans equations like we did in Chapter 5.4. Integrating this equation over the momenta gives a continuity equation \begin{equation}\label{eq-collisionless-continuity} \frac{\partial \delta}{\partial t} + a^{-1}\,\nabla\cdot\langle\vec{v}\rangle = 0\,, \end{equation} where \(\nabla\) is the gradient with respect to \(\vec{x}\), \(\langle\vec{v}\rangle\) is the mean velocity, and \(\delta\) is the overdensity in the comoving number density \(n(\vec{x})\) (or, equivalently, the mass density), that is, \(n(\vec{x}) = \overline{\rho}\,a^3\,\left(1+\delta\right)/m\). Similarly, multiplying Equation (17.18) by \(v_i\) and integrating over all momenta gives \begin{equation} \frac{\partial}{\partial t} \left(\left[1+\delta\right]\,\langle v_i\rangle\right)+\frac{\dot{a}}{a}\,\left(1+\delta\right)\,\langle v_i\rangle = -\frac{\left(1+\delta\right)}{a}\,\frac{\partial \Phi}{\partial x_i}-a^{-1}\,\sum_j\,\frac{\partial}{\partial x_j}\left(\left[1+\delta\right]\langle v_i\,v_j\rangle\right)\,. \end{equation} If we then multiply the continuity equation (17.19) by \(\langle v_i\rangle\) and subtract it from this equation, we obtain \begin{equation}\label{eq-collisionless-euler} \frac{\partial \langle \vec{v}\rangle}{\partial t}+\frac{\dot{a}}{a}\,\langle \vec{v}\rangle +a^{-1}\,(\langle \vec{v}\rangle\cdot\nabla)\langle\vec{v}\rangle = -a^{-1}\,\nabla \Phi -\frac{1}{a\left(1+\delta\right)}\,\nabla\cdot \left(\left[1+\delta\right]\,\boldsymbol{\sigma}^2\right)\,, \end{equation} where \(\boldsymbol{\sigma}^2\) is the stress tensor \(\sigma_{ij}^2 = \langle v_i\,v_j\rangle - \langle v_i\rangle\,\langle v_j\rangle\). The stress tensor acts as a type of pressure—in the next subsection, we derive the evolution equations for baryons and will see that this equation is the same as the Euler equation when we interpret \(\sigma_{ij}^2\) as pressure. That dark matter is cold means that \(\boldsymbol{\sigma}^2 \ll \langle \vec{v}\rangle^2\) and this remains so as long as \(\delta \ll 1\) (dark matter only “heats” up in this sense once it falls into a collapsing halo and then \(\delta \gg 1\)). Therefore, Equation (17.21) for cold dark matter simplifies to \begin{equation}\label{eq-collisionless-euler-cold} \frac{\partial \langle \vec{v}\rangle}{\partial t}+\frac{\dot{a}}{a}\,\langle \vec{v}\rangle +a^{-1}\,(\langle \vec{v}\rangle\cdot\nabla)\langle\vec{v}\rangle = -a^{-1}\,\nabla \Phi\,. \end{equation}
Equations (17.19) and (17.22) are the basic equations that govern the evolution of dark matter density fluctuations. As long as perturbations are small, we can keep only terms linear in the perturbation (noting that peculiar velocities are perturbative as well: \(\langle\vec{v}\rangle \sim \delta\)) and these equations simplify to \begin{align}\label{eq-dm-evol-0} \frac{\partial \delta}{\partial t} + a^{-1}\,\nabla\cdot\langle\vec{v}\rangle & = 0\,,\\ \frac{\partial \langle \vec{v}\rangle}{\partial t}+\frac{\dot{a}}{a}\,\langle \vec{v}\rangle & = -a^{-1}\,\nabla \Phi\,.\label{eq-dm-evol-1} \end{align} Taking the time derivative of the first of these, we can use the second equation to substitute the mean velocity derivative. Then employing the Poisson equation (17.17), we find that \begin{equation}\label{eq-ddotdelta-dm} \frac{\partial^2 \delta}{\partial t^2} +2\,\frac{\dot{a}}{a}\,\frac{\partial \delta}{\partial t} = 4\pi\,G\,\overline{\rho}\,\delta\,. \end{equation} For a given cosmology—which specifies \(a(t)\) through the Friedmann equations—we can solve this second-order differential equation for the evolution of small perturbations of cold dark matter. Operating Equation (17.24) with \(\nabla \times\) we find that \begin{equation} \frac{\partial \nabla \times \langle \vec{v}\rangle}{\partial t}+\frac{\dot{a}}{a}\,\nabla \times \langle \vec{v}\rangle = 0\,, \end{equation} because \(\nabla \times (\nabla \Phi) = 0\)—the curl of the gradient of scalar function is always zero. Therefore \(\nabla \times \langle \vec{v}\rangle \propto a^{-1}\) and the curl of the mean velocity field decays as the Universe expands. Because there is no source of vorticity at first order, it is therefore a good assumption that the mean velocity is curl-free and can therefore be written as the gradient of a “vector potential”.
17.1.3. Evolution of small baryonic overdensities in the Universe¶
Cold dark matter is about five times more abundant than baryonic matter, so while the formalism given above describes the evolution of the dominant matter component, baryons cannot be entirely neglected. Before discussing specific solutions to the evolution equations for cold dark matter given above, we therefore describe the evolution of baryonic perturbations.
Baryons in the early Universe are in gaseous form and they can therefore not be described by the collisionless formalism above. Because collisions are so frequent, we can approximate the baryonic component as an ideal fluid—a valid assumption as long as the mean free path is much smaller than the (cosmological) scales of interest. As above, we only consider perturbations on scales small compared to the overall size of the Universe and we can therefore use a Newtonian formalism. The evolution of a fluid is given by the continuity \begin{equation}\label{eq-dmhaloform-fluideqs-continuity} \frac{\partial \rho}{\partial t} + \nabla \cdot \rho\,\vec{u} = 0\,, \end{equation} and Euler equations \begin{equation}\label{eq-dmhaloform-fluideqs-euler} \frac{\partial \vec{u}}{\partial t} + (\vec{u}\cdot \nabla)\,\vec{u} = -\frac{\nabla P}{\rho} - \nabla \Psi\,, \end{equation} where \(P\) is the pressure and \(\Psi\) the gravitational potential. These equations are again in terms of proper position-velocity coordinates \((\vec{r},\vec{u})\) and the gradients are taken with respect to \(\vec{r}\). The relation between the density and potential is given by the Poisson equation.
As above, we express these equations in terms of comoving coordinates \(\vec{x}\), in which \(\vec{u} =\dot{a}\,\vec{x}+ \vec{v}\) with \(\vec{v} = a\,\dot{\vec{x}}\) the peculiar velocity. We again express the density in terms of the perturbation away from the mean density: \(\rho(\vec{x},t) = \overline{\rho}(t)\,\left[1+\delta(\vec{x},t)\right]\), where \(\overline{\rho}(t) \propto a^{-3}\); and also express the potential \(\Psi\) in terms of the potential perturbation \(\Phi\) from Equation (17.16). Then the fluid equations and the Poisson equation become \begin{align}\label{eq-baryon-evol-preeqstat-0} \frac{\partial \delta}{\partial t} + a^{-1}\,\nabla\cdot \left(\left[1+\delta\right]\,\vec{v}\right) & = 0\,,\\ \frac{\partial \vec{v}}{\partial t}+ \frac{\dot{a}}{a}\,\vec{v}+a^{-1}\,(\vec{v}\cdot \nabla)\,\vec{v} & = a^{-1}\,\nabla \Phi-\frac{\nabla P}{a\,\overline{\rho}\left[1+\delta\right]}\,,\label{eq-baryon-evol-preeqstate}\\ \nabla^2 \Phi & = 4\pi\,G\,\overline{\rho}\,a^2\,\delta\,.\label{eq-baryon-evol-preeqstat-2} \end{align} In these equations, the spatial gradients are now with respect to \(\vec{x}\).
To solve the fluid–Poisson equations, we need one more equation: an equation of state \(P = P(\rho,T)\) that relates the pressure to the density and temperature of the gas \(T\), or equivalently, the density and specific entropy \(S\) as \(P = P(\rho,S)\). Because we believe that the Universe and perturbations within it evolve adiabatically (or isentropically), that is without changing the entropy, many of the useful equations of state simply relate pressure to density \(P = P(\rho)\); these are barotropic equations of state. Examples are those used when computing the expansion of the Universe: radiation is described by \(P = \rho\,c^2/3\), dark energy by \(P = -\rho\,c^2\), and adiabatic dust (the simplest cosmological model for dark matter and baryons) has \(P = \rho\,c_s^2\), where \(c_s\) is the sound speed, which is zero for dark matter. Because we are interested in describing the evolution of perturbations of gaseous baryons, we will consider the latter case in more detail.
If we approximate the baryonic fluid as an ideal gas, the equation of state is \begin{equation}\label{eq-ideal-gas} P = \frac{\rho\,k_B\,T}{\mu\,m_p}\,, \end{equation} where \(k_B\) is the Boltzmann constant, \(m_p\) is the proton mass, and \(\mu\) is the mean molecular weight in terms of the proton mass (\(\mu = 1.22\) for a primordial mixture of hydrogen and helium with a helium mass fraction of 0.24). To derive the adiabatic equation of state, we use the first law of thermodynamics relating changes in the specific internal energy \(u\) to changes in entropy and volume: \begin{equation}\label{eq-first-thermo} \mathrm{d} u = T\,\mathrm{d} S - p\,\mathrm{d} V\,, \end{equation} where \(V = 1/\rho\), the volume per unit mass. For particles with \(q\) degrees of freedom, the specific internal energy is \begin{equation}\label{eq-dmhaloform-idealgas-internalenergy} u = \frac{q}{2}\,\frac{k_B\,T}{\mu\,m_p}\,, \end{equation} and for a mono-atomic gas specifically, \(q=3\), such that \(u = 3 k_B T/[2\mu\,m_p]\). Plugging this into Equation (17.33) and using the ideal gas law in Equation (17.32) we find that \begin{equation} \mathrm{d} S= \frac{k_B}{\mu\,m_p}\left(\frac{q}{2}\,\mathrm{d} \ln T - \mathrm{d}\ln \rho\right)\,. \end{equation} Integrating this and using Equation (17.32), this gives the equation of state \begin{equation} P\propto \rho^{(q+2)/q}\,\exp\left(\frac{2}{q}\,\frac{\mu\,m_p}{k_B}\,S\right)\,. \end{equation} For the isentropic evolution of a mono-atomic gas, we therefore have that \begin{equation} P \propto \rho^{5/3}\,. \end{equation}
Going back to Equation (17.30), the term involving \(\nabla P\) requires \begin{equation} \frac{\nabla P}{\overline{\rho}} = c_s^2\,\nabla \delta + \frac{2}{3}\,\left[1+\delta\right]\,T\,\nabla S\,, \end{equation} where \(c_s\) is the sound speed \begin{equation}\label{eq-dmhaloform-sound-speed} c_s = \sqrt{\frac{q+2}{q}\,\frac{k_B\,T}{\mu\,m_p}}\,. \end{equation} Equation (17.30) can then be written as \begin{equation}\label{eq-baryon-evol-general} \frac{\partial \vec{v}}{\partial t}+ \frac{\dot{a}}{a}\,\vec{v}+a^{-1}\,(\vec{v}\cdot \nabla)\,\vec{v} = a^{-1}\,\nabla \Phi-\frac{c_s^2\,\nabla \delta}{a\,\left[1+\delta\right]}-\frac{2}{3}\,a^{-1}\,T\,\nabla S\,. \end{equation} When baryonic perturbations are small, they evolve adiabatically—non-adiabatic changes are only significant once baryons fall into dark matter halos and shock, but at this point the perturbation is \(\delta \gg 1\). The evolution of small perturbations is therefore given by the \(\delta, \vec{v} \ll 1\) limit of Equations (17.29)–(17.31) with \(\nabla S = 0\) and we drop all terms that are a power of two or more of the small quantities \begin{align}\label{eq-baryon-evol-general-0} \frac{\partial \delta}{\partial t} + a^{-1}\,\nabla\cdot \vec{v} & = 0\,,\\ \frac{\partial \vec{v}}{\partial t}+ \frac{\dot{a}}{a}\,\vec{v} & = a^{-1}\,\nabla \Phi-a^{-1}\,c_s^2\,\nabla \delta\,,\label{eq-baryon-evol-general-1}\\ \nabla^2 \Phi & = 4\pi\,G\,\overline{\rho}\,a^2\,\delta\,.\label{eq-baryon-evol-general-2} \end{align} Like we did for dark matter perturbations, we can take the time derivative of the first equation and simplify it using the second and third to obtain \begin{equation}\label{eq-ddotdelta-with-pressure} \frac{\partial^2 \delta}{\partial t^2}+2\,\frac{\dot{a}}{a}\,\frac{\partial \delta}{\partial t} = 4\pi\,G\,\overline{\rho}\,\delta +\frac{c_s^2}{a^2}\,\nabla^2\delta\,. \end{equation}
Operating Equation (17.40) with \(\nabla \times\), we again obtain that \(\nabla \times \vec{v} = a^{-1}\)—now because all source terms on the right-hand side are gradients of scalar functions—and we can again approximate the velocity field as being curl-free. Comparing the evolution equations of cold dark matter and baryons, we notice that they are almost exactly the same, except for the pressure term in the baryonic Euler equation.
The equations we derived here apply only after baryons have decoupled from photons at recombination at \(z \approx 1100\). Before recombination, baryons are coupled to photons through Thomson scattering and form a baryon-photon plasma that also behaves as an ideal fluid, but whose evolution equations have additional terms due to the coupling. After recombination, the electromagnetic interactions with photons cease and because radiation stopped being a significant component of the Universe at \(z \approx 3400\) and dark energy only becomes significant at \(z<1\), as far as baryons and dark matter care, the Universe only consists of matter. The co-evolution of cold dark matter and baryons is therefore given by the fluid-like equations for the dark matter (Equations 17.23 and 17.24 for the dark matter overdensity \(\delta_c\) and velocity \(\langle \vec{v}_c\rangle\)), the fluid equations for the baryons (Equations 17.41 and 17.42 for the baryonic overdensity \(\delta_b\) and velocity \(\vec{v}_b\)), and the Poisson equation relating the gravitational potential perturbation to the total density perturbation. If we define each component’s overdensity, \(\delta_c\) and \(\delta_b\), to be relative to its own mean density, then the Poisson equation becomes \begin{equation} \nabla^2 \Phi = 4\pi\,G\,\overline{\rho}\,a^2\,\left(\frac{\Omega_c}{\Omega_m}\,\delta_c+\frac{\Omega_b}{\Omega_m}\,\delta_b\right)\,, \end{equation} where \(\Omega_c\) is the density parameter of dark matter. The Poisson equation is the only equation coupling the dark matter to the baryons, because dark matter only interacts gravitationally in standard cosmological models. Because we often do not particularly care about the evolution of the peculiar velocities, the relevant evolution equations are the coupled version of the equations for the second time derivative of the density perturbation that we derived for cold dark matter and baryons \begin{align} \frac{\partial^2 \delta_c}{\partial t^2} +2\,\frac{\dot{a}}{a}\,\frac{\partial \delta_c}{\partial t} & = 4\pi\,G\,\overline{\rho}\,\left(\frac{\Omega_c}{\Omega_m}\,\delta_c+\frac{\Omega_b}{\Omega_m}\,\delta_b\right)\,,\\ \frac{\partial^2 \delta_b}{\partial t^2}+2\,\frac{\dot{a}}{a}\,\frac{\partial \delta_b}{\partial t} & = 4\pi\,G\,\overline{\rho}\,\left(\frac{\Omega_c}{\Omega_m}\,\delta_c+\frac{\Omega_b}{\Omega_m}\,\delta_b\right) +\frac{c_s^2}{a^2}\,\nabla^2\delta_b\,. \end{align}
17.1.4. Solutions for the evolution of small matter perturbations¶
Now that we have derived the equations governing the evolution of small perturbations in the Universe (on scales less than the horizon scale and for baryons only after decoupling from the photons at recombination), we discuss aspects of their solution. In general, the evolution equations must be solved numerically—and this is what is done when initializing cosmological \(N\)-body simulations—but we can make some analytical headway that sheds light on the solutions.
17.1.4.1. Fourier decomposition¶
Because the evolution equations are linear in the perturbation for small perturbations, they are typically solved by taking the spatial Fourier transformation and solving for the evolution of each Fourier mode. This is useful, because it removes all spatial derivatives from the equations. We decompose the density perturbation (of dark matter or baryons) as \begin{equation} \delta(\vec{x},t) = \frac{1}{V}\,\sum_{\vec{k}}\delta_{\vec{k}}(t)\,\exp\left(i\vec{k}\cdot\vec{x}\right)\,, \end{equation} where we compute the coefficients \(\delta_{\vec{k}}(t)\) as the Fourier transform (see Equation B.24) \begin{equation} \delta_{\vec{k}}(t) =\int_V\,\mathrm{d}\vec{x}\,\delta(\vec{x},t)\,\exp\left(-i\vec{k}\cdot\vec{x}\right)\,, \end{equation} where \(V=L^3\) is a large volume over which the perturbation is assumed to be periodic. In the limit where \(L \rightarrow \infty\), these become \begin{align}\label{eq-dmhaloform-fourier-deltax} \delta(\vec{x},t) & = {1 \over (2\pi)^3}\,\int\,\mathrm{d}\vec{k}\,\delta_{\vec{k}}(t)\,\exp\left(i\vec{k}\cdot\vec{x}\right)\,,\\ \delta_{\vec{k}}(t) & =\int\,\mathrm{d}\vec{x}\,\delta(\vec{x},t)\,\exp\left(-i\vec{k}\cdot\vec{x}\right)\,.\label{eq-dmhaloform-fourier-deltak} \end{align}
If we take the Fourier transform of Equation (17.44)—which by setting \(c_s = 0\) includes the cold-dark-matter case—we obtain the evolution equation for the modes \(\delta_{\vec{k}}\) \begin{equation}\label{eq-ddotdelta-with-pressure-ink} \frac{\partial^2 \delta_{\vec{k}}}{\partial t^2}+2\,\frac{\dot{a}}{a}\,\frac{\partial \delta_{\vec{k}}}{\partial t} = 4\pi\,G\,\overline{\rho}\,\delta_{\vec{k}} -\frac{k^2\,c_s^2}{a^2}\,\delta_{\vec{k}}\,, \end{equation} where we write \(k^2 \equiv |\vec{k}|^2\). For cold dark matter, this equation does not contain \(\vec{k}\) and the evolution of each Fourier mode \(\delta_{\vec{k}}\) is identical: perturbations on all scales grow or decay at the same rate. For baryons we can re-write Equation (17.52) as \begin{equation}\label{eq-ddotdelta-with-pressure-ink-baryons} \frac{\partial^2 \delta_{\vec{k}}}{\partial t^2}+2\,\frac{\dot{a}}{a}\,\frac{\partial \delta_{\vec{k}}}{\partial t} = 4\pi\,G\,\overline{\rho}\,\left(1 -\frac{k^2}{k_J^2}\right)\,\delta_{\vec{k}}\,. \end{equation} where \(k_J\) is defined as \begin{equation}\label{eq-dmhaloform-jeansk} k^2_J = \frac{4\pi\,G\,\overline{\rho}\,a^2}{c_s^2}\,, \end{equation} and corresponds to the Jeans length \(\lambda_J\) \begin{equation}\label{eq-dmhaloform-jeanslength} \lambda_J = \frac{2\pi a}{k_J} = c_s\,\sqrt{\frac{\pi}{G\,\overline{\rho}}}\,. \end{equation}
The time \((G\overline{\rho}/\pi)^{-1/2}\) is the typical time associated with the expansion of the Universe in the matter-dominated era—it is the dynamical time of the Universe, see Chapter 2.3. Therefore, the Jeans length is the length over which a sound wave can propagate over the typical expansion time scale. On large spatial scales with \(k \ll k_J\), Equation (17.53) reduces to the equation for cold dark matter and again does not contain \(k\); baryonic perturbations with \(k \ll k_J\) therefore grow or decay at the same rate regardless of scale as well. At scales smaller than the Jeans length, \(k \gtrsim k_J\), the baryonic pressure becomes important and the evolution of perturbations is scale dependent. We discuss the evolution of systems on scales smaller and larger than the Jeans length in more detail in Section 17.3.3 below.
It is also useful to understand the evolution of the peculiar velocities and of the gravitational-potential perturbation. Because the velocity dispersion is small for cold dark matter, we simply write \(\vec{v}\) instead of \(\langle \vec{v}\rangle\) for dark matter as well. The Fourier transform of the continuity equation (17.23) (or Equation [17.41]) gives \begin{equation}\label{eq-continuity-ink} \frac{\partial \delta_{\vec{k}}}{\partial t} -i\, a^{-1}\,\vec{k}\cdot \vec{v}_{\vec{k}} = 0\,. \end{equation} Because we have shown that we can approximate the peculiar velocity field as being curl-free, we can write the peculiar velocity in terms of a scalar velocity potential \(\Xi\) through \(\vec{v} = \nabla \Xi\) and the Fourier transform of this is \(\vec{v}_{\vec{k}} = -i\vec{k}\,\Xi_{\vec{k}}\). Substituting this into Equation (17.56), we can derive \begin{equation}\label{eq-velocity-perturbation-ink} \vec{v}_{\vec{k}} = i\,a\,\frac{\vec{k}}{k^2}\,\frac{\partial \delta_{\vec{k}}}{\partial t}\,. \end{equation} For the gravitational potential, we can take the Fourier transform of the Poisson equation (17.17) and find \begin{equation}\label{eq-dmhaloform-poisson-comoving-perturb-k} \Phi_{\vec{k}} = -\frac{4\pi\,G\overline{\rho}\,a^2}{k^2}\,\delta_{\vec{k}}\,. \end{equation}
Once the evolution of the density perturbation is known, therefore, the evolution of the peculiar velocities or the potential perturbation can easily be derived. The evolution of these quantities will typically be scale dependent (due to the appearance of \(\vec{k}\) in their equations).
17.1.4.2. General properties of the evolution of pressure-free perturbations¶
The evolution equation (17.25) for small, pressure-free matter perturbations—for cold dark matter or baryons on scales larger than the Jeans length—is a second-order equation and therefore has two solutions. We denote these as \(\delta_+(t)\) and \(\delta_-(t)\), because one of these will typically be a solution in which the perturbation increases with time—a growing mode—and one will typically decrease with time—the decaying mode. Because the part of the initial density fluctuation that follows the decaying mode quickly becomes vanishingly small, at late times, the growing mode is the only relevant one for the formation of structure and collapsed dark matter halos. We have suppressed the spatial dependence of the solution here, because pressure-free perturbations evolve in the same way regardless of scale as we demonstrated above.
We can show that the two solutions satisfy the following relation \begin{equation}\label{eq-deltapm-equation} \delta_+\,\frac{\partial \delta_-}{\partial t}-\frac{\partial \delta_+}{\partial t}\,\delta_- \propto a^{-2} \end{equation} by taking the time derivative of the left-hand side of this equation, using Equation (17.25) to simplify the expression, and solving for the left-hand side. Thus, if we have found one solution to Equation (17.25), we can obtain the other one by solving the first-order differential equation (17.59).
We can further find one general solution as follows. The time evolution of the Hubble function \(H = \dot{a}/a\) in the matter-dominated era satisfies the second Friedmann equation (17.2), which we can write as \begin{equation}\label{eq-Hubble-timederiv} \frac{\mathrm{d} H}{\mathrm{d} t} + H^2 = -\frac{4\pi G}{3}\,\overline{\rho}\,, \end{equation} in the case of a matter-only Universe. Taking the time derivative of this equation, we find that
\begin{equation}\label{eq-Hubble-timederiv2} \frac{\mathrm{d}^2 H}{\mathrm{d} t^2} + 2\,\frac{\dot{a}}{a}\,\frac{\mathrm{d}H}{\mathrm{d}t} = 4\pi G\,\overline{\rho}\,H\,, \end{equation}
where we have used that \(\overline{\rho}(t) \propto a^{-3}\). Because the presence of dark energy in the form of a cosmological constant adds a constant to the right-hand side of Equation (17.60), Equation (17.61) in fact holds even when dark energy becomes important (at least if dark energy is a cosmological constant). Therefore, Equation (17.61) holds for the Universe we live in at all times after radiation becomes unimportant (after \(z\approx 3400\)). In fact, because the spatial curvature does not appear in the second Friedmann equation, Equation (17.61) holds even in the presence of non-zero spatial curvature. Equation (17.61) has the same form as the perturbation evolution equation (17.25) and \(H(t)\) is therefore a solution! Because the Hubble function generally declines with time, we have that \(\delta_-(t) \propto H\) (although, as we will see below, in a dark-energy dominated Universe, \(\delta_+(t) \propto H\)).
The growing mode can be obtained by substituting \(\delta_-(t) \propto H\) into Equation (17.59), which leads to \begin{equation}\label{eq-growing-from-hubble} \delta_+(t) \propto H(t)\,\int_0^t\,\frac{\mathrm{d}t'}{\dot{a'}^{2}(t')}=H(a)\,\int_0^a\,\frac{\mathrm{d}a'}{\dot{a'}^{3}}=H(a)\,\int_0^a\,\frac{\mathrm{d}a'}{a'^3\,H^3(a')}\,. \end{equation} These integrals in general do not have a closed-form solution.
17.1.4.3. Dark matter perturbations in Einstein–de Sitter and other Universes¶
In the Einstein–de Sitter model we can solve the perturbation evolution equations analytically. The Einstein–de Sitter model is the flat, matter-only model with \(\Omega_m = 1\) at all times. Because observations are consistent with a flat Universe, the Einstein-de Sitter model is in fact a good description of the Universe between the end of radiation domination at \(z \approx 3400\) and the start of the dark-energy era at \(z \lesssim 1\). This encompasses most of the era during which galaxies form.
The background evolution in the Einstein–de Sitter Universe is given by Equations (17.9). The decaying mode is therefore \(\delta_-(t) \propto t^{-1}\). The growing mode is found from Equation (17.62) and is given by \begin{equation}\label{eq-dmhaloform-growingmode-eds} \delta_+(t) \propto t^{2/3} \propto a(t)\,. \end{equation} Density perturbations in the Einstein–de Sitter Universe therefore grow as the scale factor.
In cosmology, the growth of dark-matter perturbations after the end of the radiation era is generally described in terms of the (linear) growth factor \(D_+(a)\). The growth factor is nothing more than the growing mode \(\delta_+(t)\) normalized such that \(D_+(a) = a\) during matter domination (that is, Equation 17.63 with an equal sign rather than a proportionality sign). In the Einstein–de Sitter Universe, we therefore have that \begin{equation}\label{eq-dmhaloform-growthfactor-eds} D_+(a) = a \end{equation}
To calculate the growth factor for a more general Universe made up of matter and dark energy (which could have non-zero spatial curvature), we need to normalize Equation (17.62) such that it gives \(D_+(a) = a\) during matter domination. Because \(H^2(a) = H_0^2\,\Omega_{0,m}\,a^{-3}\) during matter domination (Equation 17.4), the correctly-normalized solution is \begin{equation}\label{eq-dmhaloform-growthfactor-general} D_+(a) = {5\Omega_{0,m} \over 2}\,{H(a)\over H_0}\,\int_0^a\,\frac{\mathrm{d}a'}{(a'\,H(a')/H_0)^3}\,. \end{equation} In general, this equation needs to be solved numerically, but we can make some headway towards a more general solution by noting that the ratio \(D_+(a)/a\) only depends on the scale factor through the density parameters \(\Omega_i(a)\) evaluated at that scale factor. We leave it as an exercise to demonstrate that \begin{equation}\label{eq-dmhaloform-growthfactor-overa-general} {D_+(a) \over a} = {5\Omega_{m}(a) \over 2}\,\int_0^1\,\frac{\mathrm{d}x}{\left(x\,\sqrt{\Omega_{m}(a)\,x^{-3} + \Omega_{\Lambda}(a)+\Omega_{r}(a)\,x^{-4}+\Omega_{k}(a)\,x^{-2}}\right)^3}\,. \end{equation} This form allows us to create simple and quite general fitting functions for \(D_+(a)\). For example, for plausible values of \(\Omega_{0,m}\) and \(\Omega_{0,\Lambda}\) in a flat Universe, a good approximation is given by \begin{equation}\label{eq-dmhaloform-growthfactor-general-approx} D_+(t) = {5\Omega_m \over 2}\,\left(\Omega_m^{4/7}-\Omega_\Lambda+\left[1+\frac{\Omega_m}{2}\right]\,\left[1+\frac{\Omega_\Lambda}{70}\right]\right)^{-1}\,a\,, \end{equation} where all of the density parameters on the right-hand side depend on the scale factor (Carroll et al. 1992). For the Planck Collaboration et al. (2020b) cosmological parameters, we can compute the growth factor as a function of redshift as shown in Figure 17.2.
[6]:
Om0= 0.3111
Ode0= 0.6889
H0= 67.70*u.km/u.s/u.Mpc
from scipy import integrate
def H(a):
# Assumes flat Universe
return H0*numpy.sqrt(Om0/a**3.+Ode0)
def growth_factor(a):
# Direct calculation
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)
def growth_factor_approx(a):
def Oma(a):
return Om0/a**3./(Om0/a**3.+Ode0)
def Odea(a):
return Ode0/(Om0/a**3.+Ode0)
return 5.*Oma(a)/2.*a/(Oma(a)**(4./7.)-Odea(a)+(1.+Oma(a)/2.)*(1.+Odea(a)/70.))
figure(figsize=(7,5))
a= numpy.geomspace(3e-3,1.,101)
galpy_plot.plot(1./a,[growth_factor(ap)/ap for ap in a],
semilogx=True,
xlabel=r'$1+z$',
ylabel=r'$D_+/a$',
xrange=[400.,0.8],yrange=[0.,1.1],
gcf=True,
label=r'$\mathrm{Numerical\ calculation}$')
plot(1./a,growth_factor_approx(a)/a,
label=r'$D_+ = {5\,\Omega_m \over 2}\,\left(\Omega_m^{4/7}-\Omega_\Lambda+\left[1+\frac{\Omega_m}{2}\right]\,\left[1+\frac{\Omega_\Lambda}{70}\right]\right)^{-1}\,a$')
from matplotlib.ticker import FuncFormatter
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
legend(frameon=True,fontsize=15.,loc='lower right',
facecolor='w',edgecolor='w');
Figure 17.2: The growth factor in our Universe.
We compare the direct numerical calculation to the approximation from Equation (17.67) and we see that they agree remarkably well (to better than \(0.1\%\) everywhere in fact). Until a redshift of \(z \approx 2\), the growth rate behaves as in the matter-only limit \(D_+(a) = a\), but as the Universe transitions to a state that is dominated by dark energy, the growth rate becomes smaller and deviates from the earlier \(D_+(a) = a\) behavior. To understand this better, we can compute the decaying and growing mode for a flat Universe consisting entirely of dark energy. In such a Universe, the Hubble function is constant, \(H(a) = \dot{a}/a =\) constant, and thus one solution is \(\delta_+(t) \propto\) constant—we label this one as the “growing” mode in this case, because the other solution decays—and the other can be obtained from Equation (17.62): \(\delta_-(a) \propto a^{-2}\). Thus, once the Universe becomes dominated by dark energy, perturbations stop growing altogether. In our Universe which has a mix of matter and dark energy at the present time, the growth of perturbations therefore slows down at low redshift (when dark energy becomes important) and in the not-so-distant future, perturbations will stop growing and structure formation will end.
From Equation (17.58), we find that the potential perturbation corresponding to the growing mode in the Einstein–de Sitter model remains constant in time—it neither grows nor decays. In general, the potential perturbation is \(\Phi_{\vec{k}} \propto \delta_+(a)/a\). As dark energy becomes more important in the evolution of the Universe, \(\delta_+(a)\) starts to grow slower than the scale factor and the potential perturbation decays. Equation (17.57) for the velocity perturbation then also demonstrates that peculiar velocities become small and therefore accretion onto existing structures slows down and eventually slows to a halt. The future of our Universe—assuming dark energy does not change in any way, e.g., by decaying—is therefore one where existing structures become isolated, surrounded by a web of small, frozen-in matter perturbations.
This evolution is nicely illustrated by the following movie by Diemer & Mansfield, which shows the evolution of perturbations in a \(\approx 100\times100\,\mathrm{Mpc}^2\) comoving box from \(a=0.1\) (\(z \approx 10\)) into the future (until \(a = 10\)).
17.1.5. The evolution of small perturbations from the early Universe to the present day¶
To fully compute and characterize the density perturbations that grow into dark-matter halos, we need to consider their evolution all the way from the initial time at which they formed until they enter the scale-independent growth after radiation ceases to be a significant component of the Universe’s energy density. As we already stated above, this is the domain of early-Universe cosmology and is covered in detail in standard cosmology textbooks (e.g., Dodelson & Schmidt 2020). We will not discuss this in detail, but we will go over the important aspects of the evolution of perturbations in the early Universe and cover the results that we need to study the formation of dark-matter halos.
As we already discussed in Section 17.1.2 above, in the most popular cosmological paradigm, density perturbations are created during an epoch of inflation. In the simplest inflation models, at the end of inflation (so-called “reheating”), a pattern of density perturbations is created that is a Gaussian random field with mean zero and a covariance function that does not depend on angle. These density perturbations are adiabatic or isentropic, meaning that their relative size is the same in all energy components (dark matter, baryons, photons, etc.). Because the horizon shrinks during inflation, all relevant length scales are larger than the horizon at the end of inflation and each currently-observable length scale needs to “enter the horizon” or “cross the horizon” during the subsequent expansion of the Universe. Because all energy components have fractionally the same perturbation, they can all be related to the initial curvature perturbation \(\mathcal{R}(\vec{k})\). This curvature perturbation is conserved while the length scale of the perturbation is larger than the horizon during both radiation and matter domination—that is, when the perturbation is a super-horizon mode. Because it only starts to evolve when a mode enters the horizon, it is convenient to express the initial perturbations in terms of \(\mathcal{R}(\vec{k})\). The covariance function of the curvature perturbation can be expressed as \begin{equation} \langle \mathcal{R}(\vec{k})\mathcal{R}^*(\vec{k}')\rangle = (2\pi)^3\,P_{\mathcal{R}}(k)\,\delta^{(3)}(\vec{k}-\vec{k}')\,, \end{equation} where \(k = |\vec{k}|\), because the covariance function does not depend on angle. The function \(\delta^{(3)}(\cdot)\) is a three-dimensional version of the Dirac delta function (see Appendix B.3.2). The one-dimensional function \(P_{\mathcal{R}}(k)\) is the curvature power spectrum. The power spectrum of perturbations in the density of matter is defined in a similar way through \begin{align} \langle \delta_\vec{k}\delta^*_{\vec{k}'}\rangle & = (2\pi)^3\,P_{m}(k)\,\delta^{(3)}(\vec{k}-\vec{k}')\,.\label{eq-dmhaloform-matterpowspec-define} \end{align} Because the perturbations in the curvature and the density are defined in real space as unitless relative quantities, their Fourier-space representations have units of volume (see Equation 17.51). Because of the presence of the three-dimensional delta function, which also has units of volume, the units of the power spectra are, therefore, those of volume as well, typically expressed as \(\mathrm{Mpc}^3\).
During the linear evolution of perturbations, the power spectrum evolves with time, but the distribution of perturbations remains Gaussian and independent of angle until perturbations become non-linear. While primordial deviations from Gaussianity are expected in some inflationary models, current limits on primordial Gaussianity from observations of the cosmic microwave background (CMB) are strong and rule out any significant non-Gaussianity (Planck Collaboration et al. 2020c).
The initial curvature power spectrum produced by inflation has a power-law form \begin{equation}\label{eq-dmhaloform-pkcurv-init} P_{\mathcal{R}}(k) = 2\pi^2\,A_s\,k^{-3}\,\left({k\over k_p}\right)^{n_s-1}\,, \end{equation} where \(k_p\) is a pivot scale typically defined on large scales (e.g., Planck Collaboration et al. 2020b uses \(k_p = 0.05\,\mathrm{Mpc}^{-1}\)). The dimensionless parameter \(A_s\) sets the amplitude of the initial power spectrum. The power-law exponent is written in terms of \(n_s\) and a scale-invariant spectrum is characterized by \(n_s = 1\). Inflationary models predict \(n_s\) to be close, but not quite equal to one and this expectation is borne out by CMB observations, which find (Planck Collaboration et al. 2020b) \begin{align} n_s & = 0.9665 \pm 0.0038\,;\quad \ln(10^{10}A_s) = 3.047 \pm 0.014\,. \end{align}
The perturbations seeded by inflation then grow, decay, or oscillate depending on whether they are (i) larger or smaller than the horizon, (ii) are going through radiation or matter domination, and (iii) are affected by collisions. Collisions do not affect the dark matter perturbations, but are important for photons and baryons, which until the Universe’s temperature drops enough for neutral hydrogen to remain stable (at recombination, \(z \approx 1100\)) are tightly coupled through Thomson scattering. Because photons dominate the gravitational potential, which does affect the dark matter, until we reach the end of the radiation era, computing the evolution of perturbations of any energy component is more complicated than it is during the matter- and dark-energy dominated eras that we discussed above. The main important aspects of the evolution are
Curvature perturbations are constant for modes with scales larger than the horizon; the corresponding potential perturbation is also constant in this case until matter-radiation equality, when they drop by a factor of \(9/10\). The length scales relevant for galaxy formation enter the horizon before matter-radiation equality, so the potential perturbation is constant during their super-horizon evolution. Because the curvature perturbations do not evolve before horizon entry, it does not matter for any calculation when exactly the primordial pattern of density fluctuations was created; this is why we do not have to specify when inflation ended.
For large scales that enter the horizon after matter-radiation equality, the potential remains constant during horizon crossing; because we are then in the matter era, the potential remains constant until dark energy becomes important and potentials start to decay.
For small scales, including those relevant for galaxies, that enter the horizon before matter-radiation equality, the potential starts oscillating and decaying during the radiation era (because of the acoustic modes in the photon-baryon plasma that also create the CMB); once they enter the matter era, the potential becomes constant until it again starts decaying during the dark-energy era.
As is clear from this overview, unlike the scale-independent evolution of matter perturbations in the matter-dominated era (or after recombination for large-scale baryon modes), this evolution is strongly dependent on scale. For the study of galaxy formation and other late time (\(z \ll 1000\)) applications, all of the early-Universe evolution can be conveniently summarized in a transfer function \(T(k)\), which describes the evolution of perturbations from the primordial power spectrum of Equation (17.70) until a time well into matter domination when the formalism from the previous sections applies. The transfer function is defined in terms of the relation of the potential perturbation \(\Phi(\vec{k},a)\) (where we explicitly indicate the time dependence through the scale factor) to the initial curvature perturbation \(\mathcal{R}(\vec{k})\) through \begin{equation}\label{eq-dmhaloform-potential-transfer} \Phi(\vec{k},a) = {3\over 5}\,c^2\,\mathcal{R}(\vec{k})\,T(k)\,{D_+(a)\over a}\,. \end{equation} The factor of \(3/5\) here is chosen such that \(T(k) = 1\) on large scales where until the dark-energy era, the only evolution of the potential is the drop by \(9/10\) during matter-radiation equality. The factor of \(D_+(a)/a\) fully describes the evolution of the transfer function from when matter starts dominating the potential until the present-day and the transfer function therefore encapsulates all earlier evolution. Because the evolution most strongly depends on whether or not a spatial scale enters the horizon before or after matter-radiation equality, the most important scale in the transfer function is the particle horizon at the time of matter-radiation equality, which in Fourier space is approximately given by (Eisenstein & Hu 1998; assuming that the temperature of the CMB is \(2.7255\,\mathrm{K}\), Fixsen 2009) \begin{equation} k_{\mathrm{eq}} = 7.32\times 10^{-2}\,\Omega_{0,m}\,h^2\,\mathrm{Mpc}^{-1}\,. \end{equation} For the Planck Collaboration et al. (2020b) parameters, we have that \(k_{\mathrm{eq}} = 0.0104\,\mathrm{Mpc}^{-1}\). If the Universe only consisted of radiation and non-interacting dark matter, the transfer function would only depend on \(k/k_{\mathrm{eq}}\), because all that matters for any mode is whether it enters the horizon before or after matter-radiation equality and exactly when this happens. Because baryons make up only a small fraction of the matter in the Universe and dark energy can be neglected during the early-Universe evolution that the transfer function describes, our Universe is actually not that far off from this and the transfer function is almost entirely a function of \(k/k_{\mathrm{eq}}\).
Nowadays, the transfer function is efficiently numerically calculated for any realistic cosmological model using variations on a numerical scheme originally developed by Seljak & Zaldarriaga (1996). The main of these so-called Boltzmann codes—so named because they solve the relevant Boltzmann equation—are CAMB (Lewis et al. 2000) and CLASS (Blas et al. 2011). They are mainly used to compute the anisotropies seen in the CMB for different models, but they can output the transfer function as well (note that their definition of the transfer function includes the growth factor and other factors and thus differs from our definition; to get their transfer function on the same scale as ours, one can just normalize theirs to equal one on large scales). These codes have easy-to-use Python wrappers and any serious research work should make use of them, but for our purposes, it is useful to work with simpler fitting formulae for the transfer function. The most famous of these is the transfer function from Appendix G of Bardeen et al. (1986), a paper so famous that it is generally known simply as “BBKS”. Their fitting function for the transfer function for a Universe without baryons is \begin{equation}\label{eq-dmhaloform-transfer-bbks} T(k) = {\ln \left( 1 + 2.34q \right)\over 2.34 q} \left[1 + 3.89 q + (16.1 q)^2 + (5.46 q)^3 + (6.71 q)^4 \right]^{-1/4} \, , \end{equation} where \begin{equation}\label{eq-dmhaloform-transfer-bbks-q} q = {k \over \Omega_{0,m}\,h^2\,\mathrm{Mpc}^{-1}} = 7.32\times 10^{-2}\,{k \over k_{\mathrm{eq}}}\,. \end{equation}
The BBKS form interpolates between the \(T(k) \approx 1\) behavior at \(k \ll k_{\mathrm{eq}}\) to an approximate \(T(k) \propto 1/k^2\) behavior at \(k \gg k_{\mathrm{eq}}\).
As discussed above, without baryons the transfer function is only a function of \(k/k_{\mathrm{eq}}\). The presence of baryons causes small-scale modes that enter the horizon before matter-radiation equality to start to oscillate and decay, leading to a suppression of the transfer function on small scales (large \(k\)) relative to the dark-matter-only case and to the appearance of an oscillatory pattern. The suppression is the main effect of importance for our purposes and its effect on the transfer function can be accurately represented by a slight re-definition of \(q\) to take account of the baryon fraction (Sugiyama 1995) \begin{equation}\label{eq-dmhaloform-transfer-bbks-sugiyama} q = {k \over \Omega_{0,m}\,h^2\,\mathrm{Mpc}^{-1}} \,{\left(T_{0,\mathrm{CMB}} / 2.7\,\mathrm{K}\right)^2 \over \exp\left( -\Omega_{0,b} - \sqrt{2h} {\Omega_{0,b}/ \Omega_{0,m}} \right)}\,, \end{equation} where \(\Omega_{0,b}\) is the baryon density parameter and \(T_{0,\mathrm{CMB}} = 2.7255\,\mathrm{K}\) is the CMB temperature at redshift zero (Fixsen 2009). Figure 17.3 displays the BBKS transfer function with and without the Sugiyama (1995) prescription for the presence of baryons for the Planck Collaboration et al. (2020b) parameters and compares it to the exact numerical calculation that we have computed using CAMB (Lewis et al. 2000).
[8]:
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
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
k= numpy.geomspace(1e-4,20.,101)/u.Mpc
figure(figsize=(10,6))
from matplotlib import gridspec
from matplotlib.ticker import FuncFormatter
gs= gridspec.GridSpec(2,1,height_ratios=[2,1],hspace=0.)
subplot(gs[0])
loglog(k,transfer_bbks(k,model='BBKS'),label=r'$\mathrm{BBKS}$')
loglog(k,transfer_bbks(k,model='Sugiyama'),label=r'$\mathrm{Sugiyama\ (1995)}$')
ylabel(r'$T$')
k_CAMB= numpy.array([7.07651270e-06, 7.81367726e-06, 8.62763227e-06, 9.52637674e-06,
1.05187446e-05, 1.16144893e-05, 1.28243755e-05, 1.41602986e-05,
1.56353835e-05, 1.72641303e-05, 1.90625433e-05, 2.10483013e-05,
2.32409129e-05, 2.56619314e-05, 2.83351510e-05, 3.12868360e-05,
3.45460030e-05, 3.81446807e-05, 4.21182303e-05, 4.65057128e-05,
5.13502346e-05, 5.66994167e-05, 6.26058318e-05, 6.91275127e-05,
7.63285643e-05, 8.42797526e-05, 9.30592141e-05, 1.02753242e-04,
1.13457107e-04, 1.25275998e-04, 1.38326068e-04, 1.52735563e-04,
1.68646118e-04, 1.86214063e-04, 2.05612101e-04, 2.27030818e-04,
2.50680750e-04, 2.76794279e-04, 3.05628113e-04, 3.37465550e-04,
3.72619485e-04, 4.11435496e-04, 4.54294932e-04, 5.01619070e-04,
5.53873018e-04, 6.11570256e-04, 6.75277843e-04, 7.45621917e-04,
8.23293813e-04, 9.09056747e-04, 1.00375363e-03, 1.10831531e-03,
1.22376904e-03, 1.35124975e-03, 1.49201008e-03, 1.64743362e-03,
1.81904773e-03, 2.00853893e-03, 2.21776962e-03, 2.44879583e-03,
2.70388834e-03, 2.98555358e-03, 3.29656038e-03, 3.63996485e-03,
4.01914213e-03, 4.43781819e-03, 4.90010809e-03, 5.41055528e-03,
5.97417587e-03, 6.59650844e-03, 7.28367036e-03, 8.04241467e-03,
8.88019707e-03, 9.80525184e-03, 1.08266696e-02, 1.19544901e-02,
1.31997960e-02, 1.45748267e-02, 1.60930939e-02, 1.77695211e-02,
1.96205806e-02, 2.16644667e-02, 2.39212681e-02, 2.64131594e-02,
2.91646309e-02, 3.20810936e-02, 3.49975564e-02, 3.79140191e-02,
4.08304818e-02, 4.37469482e-02, 4.66634110e-02, 4.95798737e-02,
5.24963364e-02, 5.54128028e-02, 5.83292618e-02, 6.12457283e-02,
6.48791790e-02, 6.85126334e-02, 7.21460804e-02, 7.57795274e-02,
7.94129819e-02, 8.30464289e-02, 8.66798833e-02, 9.03133303e-02,
9.39467847e-02, 9.75802317e-02, 1.01213686e-01, 1.04847133e-01,
1.08480588e-01, 1.12114035e-01, 1.15747489e-01, 1.19380936e-01,
1.23014383e-01, 1.26647845e-01, 1.30281284e-01, 1.33914739e-01,
1.37548193e-01, 1.41181648e-01, 1.44815087e-01, 1.48448542e-01,
1.52081996e-01, 1.55715436e-01, 1.59348890e-01, 1.62982345e-01,
1.66615799e-01, 1.70249239e-01, 1.73882693e-01, 1.77516133e-01,
1.81149602e-01, 1.84783041e-01, 1.88416496e-01, 1.92049935e-01,
1.95683405e-01, 1.99316844e-01, 2.02950299e-01, 2.06583738e-01,
2.10217193e-01, 2.13850662e-01, 2.17484102e-01, 2.21117541e-01,
2.24750996e-01, 2.28384435e-01, 2.32017905e-01, 2.35651359e-01,
2.39284799e-01, 2.42918238e-01, 2.46551707e-01, 2.50185162e-01,
2.53818601e-01, 2.57452041e-01, 2.61085480e-01, 2.64718950e-01,
2.68352419e-01, 2.71985859e-01, 2.75619298e-01, 2.79252768e-01,
2.82886207e-01, 2.86519647e-01, 2.90153116e-01, 2.93786556e-01,
2.97420025e-01, 3.01053464e-01, 3.04686904e-01, 3.08320343e-01,
3.11953813e-01, 3.15587252e-01, 3.19220722e-01, 3.22854161e-01,
3.26487601e-01, 3.30121070e-01, 3.33754510e-01, 3.37387949e-01,
3.41021419e-01, 3.44654858e-01, 3.48288298e-01, 3.51921737e-01,
3.55555236e-01, 3.59188676e-01, 3.62822115e-01, 3.66455585e-01,
3.70089024e-01, 3.73722464e-01, 3.77355903e-01, 3.80989343e-01,
3.84622812e-01, 3.88256282e-01, 3.91889721e-01, 3.95523190e-01,
3.99156630e-01, 4.02790070e-01, 4.06423509e-01, 4.10056978e-01,
4.13690418e-01, 4.17323858e-01, 4.20957297e-01, 4.24590796e-01,
4.42282051e-01, 5.47528386e-01, 6.81588948e-01, 8.48473907e-01,
1.05622005e+00, 1.31483233e+00, 1.63676500e+00, 2.03752184e+00,
2.53640270e+00, 3.15743303e+00, 3.93052101e+00, 4.89289665e+00,
6.09090757e+00, 7.58224726e+00, 9.43873692e+00, 1.17497826e+01,
1.46266804e+01, 1.82079773e+01, 2.26661453e+01],dtype=numpy.float32)
T_CAMB= numpy.array([1.00000000e+00, 9.99999702e-01, 9.99999404e-01, 9.99998868e-01,
9.99998450e-01, 9.99997854e-01, 9.99997079e-01, 9.99996126e-01,
9.99995053e-01, 9.99993622e-01, 9.99992013e-01, 9.99989986e-01,
9.99987543e-01, 9.99984384e-01, 9.99980748e-01, 9.99976277e-01,
9.99970794e-01, 9.99963999e-01, 9.99956012e-01, 9.99945939e-01,
9.99933898e-01, 9.99919116e-01, 9.99901235e-01, 9.99879360e-01,
9.99852657e-01, 9.99820113e-01, 9.99780774e-01, 9.99732733e-01,
9.99674320e-01, 9.99603271e-01, 9.99517083e-01, 9.99412239e-01,
9.99285340e-01, 9.99131620e-01, 9.98945475e-01, 9.98720706e-01,
9.98449683e-01, 9.98123705e-01, 9.97732460e-01, 9.97264266e-01,
9.96705532e-01, 9.96041417e-01, 9.95254695e-01, 9.94326651e-01,
9.93235946e-01, 9.91959214e-01, 9.90469933e-01, 9.88738298e-01,
9.86731350e-01, 9.84411418e-01, 9.81737137e-01, 9.78663087e-01,
9.75139678e-01, 9.71113205e-01, 9.66526210e-01, 9.61318374e-01,
9.55426872e-01, 9.48787630e-01, 9.41335559e-01, 9.33005512e-01,
9.23730910e-01, 9.13444996e-01, 9.02079999e-01, 8.89569044e-01,
8.75844955e-01, 8.60859334e-01, 8.44496965e-01, 8.26697707e-01,
8.07417512e-01, 7.86588132e-01, 7.64138579e-01, 7.40005493e-01,
7.14135289e-01, 6.86479986e-01, 6.57013237e-01, 6.25738502e-01,
5.92709422e-01, 5.58052480e-01, 5.21996319e-01, 4.84905273e-01,
4.47310716e-01, 4.09928203e-01, 3.73643249e-01, 3.39446396e-01,
3.08290511e-01, 2.81759530e-01, 2.60875255e-01, 2.43905976e-01,
2.29525447e-01, 2.16738030e-01, 2.04725415e-01, 1.92877948e-01,
1.81006402e-01, 1.69263780e-01, 1.57975391e-01, 1.47459209e-01,
1.35808721e-01, 1.26063988e-01, 1.18302770e-01, 1.12255208e-01,
1.07376561e-01, 1.03062309e-01, 9.88297760e-02, 9.44029242e-02,
8.97217318e-02, 8.48901570e-02, 8.01038221e-02, 7.56007284e-02,
7.16082305e-02, 6.82572201e-02, 6.55252710e-02, 6.32714257e-02,
6.13189489e-02, 5.94981536e-02, 5.76610006e-02, 5.57114929e-02,
5.36368564e-02, 5.14971241e-02, 4.93835099e-02, 4.73856740e-02,
4.55746800e-02, 4.39881831e-02, 4.26214673e-02, 4.14345935e-02,
4.03678641e-02, 3.93554606e-02, 3.83400954e-02, 3.72894518e-02,
3.62026915e-02, 3.51008885e-02, 3.40128504e-02, 3.29680219e-02,
3.19932252e-02, 3.11059225e-02, 3.03071029e-02, 2.95819305e-02,
2.89079025e-02, 2.82629542e-02, 2.76291426e-02, 2.69941352e-02,
2.63531543e-02, 2.57101599e-02, 2.50757858e-02, 2.44623311e-02,
2.38790940e-02, 2.33306699e-02, 2.28175409e-02, 2.23368797e-02,
2.18827687e-02, 2.14472003e-02, 2.10223440e-02, 2.06028875e-02,
2.01869700e-02, 1.97755974e-02, 1.93713680e-02, 1.89773794e-02,
1.85966156e-02, 1.82314999e-02, 1.78831629e-02, 1.75507702e-02,
1.72318798e-02, 1.69237182e-02, 1.66240986e-02, 1.63312927e-02,
1.60438139e-02, 1.57609228e-02, 1.54831391e-02, 1.52117983e-02,
1.49480570e-02, 1.46924220e-02, 1.44449789e-02, 1.42056625e-02,
1.39741516e-02, 1.37496693e-02, 1.35311242e-02, 1.33175096e-02,
1.31082814e-02, 1.29033178e-02, 1.27026318e-02, 1.25062289e-02,
1.23142218e-02, 1.21269086e-02, 1.19445641e-02, 1.17671583e-02,
1.15943598e-02, 1.14257904e-02, 1.12612136e-02, 1.11004412e-02,
1.09431762e-02, 1.07890777e-02, 1.06379725e-02, 1.04899295e-02,
1.03450790e-02, 1.02034127e-02, 1.00648059e-02, 9.92919784e-03,
9.31018963e-03, 6.62348932e-03, 4.64416808e-03, 3.23948334e-03,
2.24911119e-03, 1.55490870e-03, 1.07085682e-03, 7.34918052e-04,
5.02763898e-04, 3.42945394e-04, 2.33306288e-04, 1.58330862e-04,
1.07207408e-04, 7.24401980e-05, 4.88536243e-05, 3.28878923e-05,
2.21028804e-05, 1.48313329e-05, 9.93727554e-06], dtype=numpy.float32)
loglog(k_CAMB,T_CAMB,label=r'$\mathrm{CAMB}$')
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
gca().yaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
legend(frameon=False,fontsize=18.)
subplot(gs[1])
semilogx(k_CAMB,transfer_bbks(k_CAMB/u.Mpc,model='BBKS')/T_CAMB,
label=r'$\mathrm{BBKS}$')
semilogx(k_CAMB,transfer_bbks(k_CAMB/u.Mpc,model='Sugiyama')/T_CAMB,
label=r'$\mathrm{Sugiyama\ (1995)}$')
xlabel(r'$k/\mathrm{Mpc}$')
ylabel(r'$T/T_{\mathrm{CAMB}}$')
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
legend(frameon=False,fontsize=18.);
Figure 17.3: The transfer function.
As expected, the transfer function is equal to one on large scales and starts to decline at \(k \approx k_{\mathrm{eq}}\). Without the suppression due to baryons, the BBKS transfer function is about \(50\%\) higher than the Sugiyama (1995) one. The Sugiyama (1995) one in fact closely matches the exact numerical result. Considering the transfer function relative to the exact numerical result in the bottom panel, we see that the Sugiyama (1995) transfer function is overall accurate over the entire range of \(k\) considered here. The only difference with the exact result that is left is the oscillatory pattern due to oscillations in the baryon-photon plasma before recombination. These are the so-called acoustic oscillations that are also visible in the CMB and that leave an imprint on the transfer function. As the bottom panel makes clear, they account for a \(\lesssim 10\%\) deviation from the exact result when using the Sugiyama (1995) transfer function. While these baryon acoustic oscillations have been detected in the statistical distribution of galaxies at low redshift (Eisenstein et al. 2005; Cole et al. 2005), we can ignore them when considering the overall distribution of dark matter halos, because the scales of galaxies are \(k \gtrsim 1\,\mathrm{Mpc}^{-1}\) where the oscillations are negligible.
Returning to the statistics of density fluctuations, using the Poisson equation in Fourier space (Equation 17.58), we can relate the potential perturbation to the matter perturbation at times when matter or dark energy dominates the energy density and we then find that \begin{equation} \delta_\vec{k}(a) = -{2\over 3}{k^2a\over \Omega_{0,m}\,H_0^2}\,\Phi(\vec{k},a)\,, \end{equation} where we have expressed the mean matter density in terms of the matter density parameter and the Hubble constant. Combining this equation with Equation (17.72), we have that \begin{equation} \delta_\vec{k}(a) = -{2 \over 5}\,{c^2\,k^2\over \Omega_{0,m}\,H_0^2}\,\mathcal{R}(\vec{k})\,T(k)\,D_+(a)\,. \end{equation} From Equation (17.69), we then see that the linear matter power spectrum \(P_{m}(k,a)\) is given by \begin{equation}\label{eq-dmhaloform-linmatpowspec} P_{m}(k,a) = {8\pi^2 \over 25}\,{A_s\,k_p\,c^4 \over \Omega_{0,m}^2\,H_0^4}\,T^2(k)\,D_+^2(a)\,\left({k\over k_p}\right)^{n_s}\,. \end{equation} The linear matter power spectrum describes the statistics of small matter fluctuations at all redshifts \(z \lesssim\) a few hundred.
Because of the scale-independent growth of density perturbations in the matter- and dark-energy-dominated eras, the only scale dependence of the matter power spectrum comes from the initial curvature power spectrum on the one hand, which contributes the \(k^{n_s}\) factor, and the transfer function squared on the other hand. On large scales \(k \lesssim k_{\mathrm{eq}}\), \(T(k) \approx 1\) and \(P_{m}(k)\) is therefore approximately \(\propto k^{n_s} \approx k\) for our Universe. On small scales, the transfer function is approximately \(T(k) \propto 1/k^2\) and \(P_{m}(k)\), therefore, behaves as \(\propto k^{n_s-4} \approx k^{-3}\) for our Universe. In between, the matter power spectrum peaks. Using the BBKS transfer function with the Sugiyama (1995) prescription for \(q\) and numerically calculating the maximum, we find that the peak occurs at \(q \approx 0.095\) or \(k \approx 1.3\,k_{\mathrm{eq}}\,\exp\left( -\Omega_{0,b} - \sqrt{2h} {\Omega_{0,b}/ \Omega_{0,m}} \right) \approx 1.03\,k_{\mathrm{eq}}\).
[9]:
ns= 0.9665
from scipy import optimize
qopt= optimize.minimize_scalar(lambda q: -transfer_bbks(q/u.Mpc,model='Sugiyama',
Om0=1.,h=1.,Ob0=0.)**2.*q**ns,
bounds=[0.01,1.],options={'disp': False})
print(f"The linear matter power spectrum peaks at q = {qopt.x:.3f}")
The linear matter power spectrum peaks at q = 0.095
Thus, the matter power spectrum peaks approximately at the scale that enters the horizon at matter-radiation equality.
Using the BBKS transfer function with the Sugiyama (1995) prescription for \(q\), Figure 17.4 shows the linear matter power spectrum at redshift zero (\(a = 1\)). We also include various measurements of the linear matter power spectrum, coming from CMB and large-scale structure observations. While the raw measurements are not directly of the matter power spectrum, they can be expressed as measurements of the linear matter power spectrum using an approach pioneered by Tegmark & Zaldarriaga (2002). In particular, we make use of the recent compilation (and code) from Chabanier et al. (2019a) as also used in Planck Collaboration et al. (2020a). This compilation includes measurements originating from the temperature (TT), polarization (EE), and lensing (\(\Phi\Phi\)) of the CMB from Planck Collaboration et al. (2020a), cosmic shear measurements from the DES (Troxel et al. 2018), measurements coming from the large-scale distribution of massive galaxies in SDSS (Reid et al. 2010), and from the large-scale distribution of hydrogen gas traced through Lyman-\(\alpha\) absorption (Chabanier et al. 2019b).
[10]:
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
# First the code to get the growth factor
from scipy import integrate
from astropy.constants import c
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)
def matter_power_spectrum(k,a=1.,model='BBKS'):
return 8.*numpy.pi**2./25.*c**4.*As*kp/Om0**2./H0**4.*transfer_bbks(k,model=model)**2.\
*growth_factor(a)**2.*(k/kp)**ns
figure(figsize=(10,6))
k= numpy.geomspace(1e-4,5.,101)/u.Mpc
loglog(k,matter_power_spectrum(k,model='Sugiyama').to_value(u.Mpc**3),color='k',lw=1.5)
xlabel(r'$k\,(1/\mathrm{Mpc})$')
ylabel(r'$P_m\,(\mathrm{Mpc}^3)$')
# Planck 2018
ebar_kwargs= {'capsize':4,'elinewidth':2,
'markeredgewidth':1,
'marker':"o",'markersize':4}
ttk= h*numpy.array([0.0002785081327527984, 0.00047662761084573817, 0.0009202633962783727, 0.002112758472891597, 0.007986390185155652, 0.012435748941442042, 0.017655889415993244, 0.026243522221950928, 0.03613746066104165, 0.06583739836155929, 0.07400419971914153, 0.08690244233435813, 0.09791784405247703, 0.12659894842891328, 0.14365120047230187, 0.16508078280514932, 0.18350418598968113, 0.22250217256970695])
ttpk= h**-3*numpy.array([248.92872532592318, 1970.0911324845895, 3196.251475915504, 6834.8756344590665, 20028.807453540154, 23038.86291255816, 24327.231899744093, 21285.94122592114, 16564.446791692964, 10081.922895609325, 9102.315651666275, 7033.931572749942, 5559.22343810523, 4195.391692520878, 3387.1448698845447, 2572.9135922015275, 2154.0807306810784, 1561.364429537585])
ttkerr= h*numpy.array([[9.742282627334879e-5, 0.00013228371533774745, 0.00018788321790492654, 0.0006698523247851576, 0.0024100846252400948, 0.0032564584077827323, 0.004281591627809571, 0.004388991440177019, 0.005256909786031436, 0.01900578100567097, 0.016481079535570244, 0.017877174786641657, 0.009176575339279342, 0.021917548807312743, 0.016287890406049034, 0.013554364052033127, 0.02141118380990593, 0.02212132199113972],
[0.0001418182876166154, 0.0002185711263789327, 0.00038069566696771046, 0.0017296308409131, 0.006497292781372328, 0.00892931857858187, 0.011102093831901855, 0.013622631676662849, 0.030032136211824698, 0.030907027925588297, 0.030834300237336507, 0.03083871246273196, 0.03350900374352235, 0.033511033841113855, 0.03265286789394475, 0.033741362623676496, 0.033450133546994096, 0.03261172648458402]])
ttpkerr= h**-3*numpy.array([2644.302755531956, 1304.854046878472, 1001.5025986799458, 735.8578499951388, 765.0564740624235, 719.4854471689332, 367.5988679299325, 205.59450972894737, 138.44198641236002, 73.6148986893802, 54.66047482195898, 29.32562046457404, 19.95447845309284, 16.378025207430667, 15.266087660656067, 17.022041124423378, 22.01079468140679, 41.61084386598731])
errorbar(ttk,ttpk,yerr=ttpkerr,xerr=ttkerr,ls='none',label=r"$\mathrm{Planck\ 2018\ TT}$",zorder=-1,**ebar_kwargs)
eek= h*numpy.array([0.005279721325386101, 0.008061638344167057, 0.012574918345876276, 0.03119663160162892, 0.036768096062207725, 0.043402784101814695, 0.05116638569193338, 0.07384181032416533, 0.08491577651129392, 0.10968330119433234, 0.13515187437147763, 0.14641540516267973, 0.17629433086141355, 0.20984001238534594])
eepk= h**-3*numpy.array([21617.597888302338, 21931.46742128491, 23130.6678781801, 20043.555495896406, 16367.133962000753, 13554.702925729604, 12316.011652554476, 8889.165724074874, 7418.388421344758, 4694.622526244998, 3928.9727759724483, 3162.6866971840464, 2462.7355250982914, 2192.0618060222987])
eekerr= h*numpy.array([[0.0010077017569064953, 0.0010606954150645632, 0.002285695946883118, 0.012065346231163569, 0.003182797104423277, 0.0033678125436359717, 0.0030852160666537715, 0.004373445680041228, 0.005795938535648574, 0.005153754914506671, 0.01720646291003007, 0.005862558894875123, 0.0071761243506671835, 0.014949740492594059],
[0.0010180278206637229, 0.0013335401523355953, 0.0029824635808796966, 0.009331396247763978, 0.005743572247744892, 0.0035031046829741993, 0.013541047965053291, 0.006106914853492273, 0.021666614310030686, 0.005719273007211615, 0.012089552409595698, 0.010411220316772832, 0.00960445596736717, 0.01257501688171483]])
eepkerr= h**-3*numpy.array([3128.971322464853, 1555.9801625254102, 801.7815852513377, 651.3329121692252, 218.9652047485547, 139.15827451382117, 154.20075750311688, 76.94929880143962, 78.43284551676486, 59.12652382376625, 136.6953745991443, 140.55694826925372, 230.88004832834042, 475.38675980919817])
errorbar(eek,eepk,yerr=eepkerr,xerr=eekerr,ls='none',label=r"$\mathrm{Planck\ 2018\ EE}$",zorder=-2,**ebar_kwargs)
ppk= h*numpy.array([0.0047718232216222185, 0.011875602223730219, 0.023360766272873968, 0.05644140824347293, 0.10663696554841638, 0.18180326513588335, 0.3164144490629016])
pppk= h**-3*numpy.array([16225.62463121451, 24952.720129016016, 23994.57288271041, 10723.827345151138, 5119.484477447558, 2239.0710384061517, 605.842817397391])
ppkerr= h*numpy.array([[0.0023133933289078666, 0.005642300376655263, 0.011196893152268515, 0.026571679597277954, 0.04944550116347209, 0.08362166202832513, 0.14176500221506955],
[0.007858501567059424, 0.0192382828963272, 0.03769333768449379, 0.08761974595910088, 0.1577662785116775, 0.2527241969081614, 0.3962943271544918]])
pppkerr= h**-3*numpy.array([2702.8538465846955, 2300.087811856112, 1042.083657438043, 482.6130544248406, 302.3304666922485, 205.75884440199584, 203.30577866133117])
errorbar(ppk,pppk,yerr=pppkerr,xerr=ppkerr,ls='none',label=r"$\mathrm{Planck\ 2018}\ \Phi\Phi$",zorder=-3,**ebar_kwargs)
# DES
desk= h*numpy.array([2.0101030356219147, 0.9962461119832121, 0.5822323365284601, 0.38419189750535293, 0.22460652053883198])
despk= h**-3*numpy.array([13.389165911971853, 52.11729902060237, 221.7618198303865, 558.7982950057899, 1618.3109312967345])
deskerr= h*numpy.array([[0.7483077541632812, 0.28954166515069946, 0.1270297581161739, 0.11368869556823183, 0.09426428653076101],
[2.4833332470136087, 1.0116343911185608, 0.4185051045994231, 0.19229772403984108, 0.15560670393946008]])
despkerr= h**-3*numpy.array([2.448264380453913, 7.4247354088166055, 20.047085716357092, 47.069556982953685, 165.45183061426897])
errorbar(desk,despk,yerr=despkerr,xerr=deskerr,ls='none',label=r"$\mathrm{DES\ Y1\ cosmic\ shear}$",**ebar_kwargs)
#SDSS LRGs
lrgk= h*numpy.array([0.03285323826500297, 0.03882664884135473, 0.04479991087344029, 0.050773172905525854, 0.05674658348187761, 0.06271984551396316, 0.06869310754604874, 0.07466651812240048, 0.08063978015448604, 0.0866131907308378, 0.09258645276292336, 0.09855971479500893, 0.10453312537136067, 0.11050638740344623, 0.116479797979798, 0.12245306001188357, 0.12842632204396912, 0.13439973262032087, 0.14037299465240644, 0.146346256684492, 0.1523202614379085, 0.1582932263814617, 0.16426619132501488, 0.17023915626856806, 0.17621360665478314, 0.18218657159833632, 0.18815953654188952, 0.1941325014854427, 0.20010546642899588, 0.20607991681521096, 0.21205288175876413, 0.2180258467023173, 0.22399881164587052, 0.2299732620320856, 0.23594622697563877, 0.24191919191919195, 0.24789215686274513, 0.2538666072489602, 0.2598395721925134, 0.26581253713606656, 0.27178550207961977, 0.27775995246583485, 0.28373291740938805, 0.2897058823529412, 0.2956788472964944])
lrgpk= h**-3*numpy.array([20193.404417754333, 16283.939710724471, 12407.648706360425, 11631.231017349326, 11219.407559535824, 10191.234190552022, 9204.180152690085, 9992.333192404913, 8860.871283357645, 7235.528073974297, 6524.89067902336, 5661.911434275714, 5224.583242363387, 4585.059429700386, 4452.481593889421, 4348.74654914214, 4173.463857579495, 3702.811494492947, 3701.3554816773108, 3303.5325412606626, 3100.1808623855895, 2871.6670743511822, 2565.948726897393, 2462.591599425841, 2508.2001137544553, 2262.040698631979, 2179.0213959213734, 2120.244939217837, 1908.8099629006747, 1795.483675743998, 1795.6399359805473, 1684.9033956892754, 1512.968774700264, 1503.4324373164457, 1480.9891989780413, 1402.5854055879392, 1318.4985556236018, 1256.3829713467812, 1250.8681522287304, 1167.549354215755, 1124.2786855600032, 1049.6993736423503, 1029.9978488459979, 965.435392064056, 923.7337985084602])
lrgpkerr= h**-3*numpy.array([2289.279489903156, 1617.0351344115438, 1224.103378449522, 1004.9321163292992, 881.3621404149991, 798.3655252201441, 718.5489585323758, 634.9321250381231, 562.6749984430783, 490.35675518490996, 426.814682384283, 366.12487535994825, 327.5287782922301, 298.28253516820206, 280.04112719663914, 268.003518778657, 259.8993794472954, 251.0890527064996, 239.67430550048775, 225.8223769229168, 211.26973445931563, 195.18042256662633, 181.98286371684975, 171.59320748481483, 162.53137010489934, 154.49383054178662, 147.44523210046705, 140.19063511635244, 134.0825662124791, 128.51629211807034, 122.96960270950981, 117.82440593942808, 113.12871914092148, 110.5351277952558, 107.72153255486678, 106.50493533836806, 104.99538438546577, 102.51082432676569, 99.22433791616348, 94.80342020488322, 90.56096573909335, 86.40780032372149, 82.18992066939066, 78.35608389143455, 75.63061163858366])
errorbar(lrgk,lrgpk,yerr=lrgpkerr,ls='none',label=r"$\mathrm{SDSS\ DR7\ LRG}$",zorder=5,**ebar_kwargs)
#eBOSS Lya
lyak= h*numpy.array([0.30858868241817705, 0.49947009591095926, 0.6903515094037415, 0.8812329228965237, 1.0721143363893058, 1.262995749882088, 1.4538771633748704, 1.6447585768676523, 1.8356399903604348, 2.026521403853217, 2.2174028173459988, 2.408284230838781, 2.5991656443315634, 2.7900470578243457, 2.980928471317128, 3.17180988480991, 3.362691298302693, 3.5535727117954745, 3.744454125288257])
lyapk= h**-3*numpy.array([803.7027763677258, 321.9397736037925, 246.26769395906535, 119.66402999801294, 50.38059776050974, 61.617123985136644, 27.764507481799036, 24.16406664322154, 19.05679226534816, 11.804015286601775, 11.556402327395114, 7.688349706515366, 7.546176793320766, 6.088068633044457, 6.9705857322592415, 3.799450100162277, 7.140928984472065, 3.4364452937674415, 4.986717538057116])
lyapkerr= h**-3*numpy.array([263.2175666929971, 81.81021590478728, 49.818356802341526, 18.970907055499517, 9.173003270956299, 15.586124359247478, 5.89871510748373, 5.720750938675552, 4.801322218693053, 3.42812783972946, 3.6989361085060946, 3.1272198101280675, 2.4245815624938385, 3.7355071775222344, 4.572269667568962, 3.675056820677179, 5.55659386516595, 1.8899521920737834, 13.35011838679237])
errorbar(lyak,lyapk,yerr=lyapkerr,ls='none',label=r"$\mathrm{eBOSS\ DR14\ Ly-}\alpha\ \mathrm{forest}$",zorder=4,**ebar_kwargs)
legend(frameon=False,fontsize=18.)
from matplotlib.ticker import FuncFormatter
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
gca().yaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)));
Figure 17.4: The linear matter power spectrum.
It is clear from this figure that all measurements are in excellent agreement with the predicted linear matter power spectrum using the Planck Collaboration et al. (2020b) cosmological parameters.
In Section 17.4 below, we will use the linear matter power spectrum at different redshifts to compute the number density of dark matter halos of different masses at different times in the Universe. Because of the scale-independent growth of perturbations in the matter- and dark-energy-dominated eras, the only scale-factor or redshift dependence of the matter power spectrum is an overall scaling \(D^2_+(a)/D_+^2(a=1)\) when starting from the \(a=1\) (\(z=0\)) matter power spectrum.
.
.
.