The inhomogeneous Universe

Dark matter: the driver for structure formation

See slides 20, 21 and 22 here

Dark matter was first proposed by astronomer Fritz Zwicky who observed the velocity dispersion of galaxies in the nearby Coma and Virgo cluster of galaxies. He calculated an estimate of the dynamical mass of the galaxy cluster by using the virial theorem

\begin{eqnarray} \sigma^2 \approx \frac{GM}{R} \end{eqnarray}

He found out that the velocity dispersion implied a dynamical mass for the cluster which was much bigger than the mass contained within the stars of the constituents galaxies. Even after accounting for the gas component (which can only be seen in X-rays and was discovered much later), the mass shortfall meant that something is missing in our understanding of the Universe. Zwicky termed this excess unseen mass, dark matter.

The presence of Dark matter was later confirmed by astronomer Vera Rubin who looked at the rotation velocities of gas in nearby spiral galaxies. If ordinary matter was the only component of matter, then the gas rotation velocities around a galaxy should have shown a Keplerian fall-off with \(v \propto R^{-1/2}\). Instead she found that the velocity of gas flattens out, implying that the mass distribution should grow as a function of the radius \(R\).

The existence of dark matter is also evident in the CMB temperature fluctuations. Without the presence of dark matter, the fluctuations seen in the temperature power spectrum will have a different characteristic. The height of the peaks in the CMB power spectrum can tell us the ratios of the baryonic matter to that of radiation as well as baryonic matter to matter which does not couple with radiation.

As you will see next the density fluctuations grow as the scale factor in the matter dominated regime. Given that we see density fluctuations in ordinary matter in the CMB of the order \(\delta ~ 10^{-5}\), these fluctuations would grow at most by a factor 1000 till today. Without the presence of dark matter there is no way they could have formed the large scale structure in the Universe which is highly non-linear with fluctuations of the order \(\delta ~ 100\). Dark matter could be a thermal relic which decouled early on in the history of the Universe, and because it decoupled from the thermal bath, its fluctuations were able to grow until the decoupling of protons, electrons and radiation (which happens at recombination).

The origin of fluctuations

Before we study the growth of fluctuations, let us talk about the origin of the density flcutations in the Universe.

At the very early stages the Universe also underwent a phase of an exponential expansion called inflation. Inflation is able to explain why the temperature of the CMB is so uniform. Regions which were in causal contact were expanded such that they now do not remain in causal contact. For the observable Universe to have the same temperature of the CMB throughout, requires about 56 efoldings of inflation. To give an order of magnitude estimate as to how big this is. Consider a physical region which is a cube of side 3 cm each. This small region will expand to a cube which has side length greater than the distance to the Andromeda galaxy.

Tiny quantum fluctuations during the inflationary era get stretched to astronomical scales, and these fluctuations act as the seeds for the formation of structure in the Universe. The fluctuations are Gaussian in nature. Let us try to understand the statistical properties of these fluctuations.

The overdensity, \(\delta(\vec{x})\) is given by

\begin{equation} \delta(\vec{x}) = \frac{\rho(\vec{x})}{\bar\rho} - 1 \end{equation}

Given that the fluctuations are Gaussian one can characterize them fully, by the mean (zero by construction), and the covariance \(\langle\delta(\vec{x})\delta(\vec{x}+\vec{r}) \rangle\). Invoking translational and rotational invariance, this quantity should only depend upon the modulus of the distance \(|\vec{r}|\). This quantity is called the two point correlation function, \(\xi(r)\).

It is often useful to work in Fourier space. The overdensity field in Fourier space is given by

\begin{equation} \delta(\vec{k}) = \int d^3{\vec x} \,\delta(\vec{x}) \,e^{-i\vec{k}.\vec{x}} \end{equation}

The correlation in Fourier space, \(\langle \delta(\vec{k})\delta(\vec{k'}) \rangle\) is then given by

\begin{eqnarray} \langle \delta(\vec{k})\delta(\vec{k'}) \rangle &=& \int d^3 \vec{x'} \int d^{3} \vec{x} \langle \delta(\vec{x}) \delta(\vec{x}') \rangle e^{-i(\vec{k}.\vec{x}+\vec{k}'.\vec{x}')} \\ &=& \int d^3 \vec{x'} e^{-i(\vec{k}+\vec{k'}).\vec{x'}} \left[ \int d^{3} \vec{y} \,\xi(\vec{y}) e^{-i\vec{k}.\vec{y}}\right] \\ &=& (2\pi)^3 P(\vec{k}) \delta_{\rm D} (\vec{k}+\vec{k'}) \end{eqnarray}

where we have changed variable to \(\vec{y}=\vec{x}-\vec{x'}\) in the inner integral, the power spectrum \(P(\vec{k})\) is the fourier transform of the correlation function \(\xi(\vec{r})\). Rotational invariance further demands that the power spectrum should only be a function of the modulus of \(\vec{k}\). The Power spectrum has dimensions of volume, so the dimensionless quantity \(\Delta(k)\) can be constructed which represents the power per logarithmic interval in \(k\),

\begin{equation} \Delta(k) = \frac{k^3}{2\pi^2} P(k) \end{equation}

Similar to the density fluctuations one can also consider potential fluctuations and the power spectrum of potential fluctuations. The potential fluctuations are related to the density fluctuations with the Poisson equation,

\begin{eqnarray} \nabla^2 \delta_{\rm \phi} \propto \delta \\ k^2 \delta_{\rm \phi}(k) \propto \delta(k) \end{eqnarray}

The power in the initial potential fluctuations per logarithmic interval of \(k\) set by inflation is constant due to scale invariance, i.e., $:nbsphinx-math:Delta{:nbsphinx-math:`phi`} = k^3 P{\phi}(k) = $ constant. This implies that the initial density power spectrum

\begin{equation} P(k) \propto k \end{equation}

The index in the power spectrum has been measured by data, and is quite close to unity although not exactly unity. Since inflation ends at a certain epoch, inflation in fact predicts a slope close to unity, but with slight deviations. This deviation has now been measured at more than 5\(\sigma\) in the data.

These fluctuations are set at the end of inflation. The Universe then enters an era of radiation domination. The density fluctuations set at the end of inflation will grow. There are two regimes to consider, one for fluctuations whose wavelengths are beyond the horizon, and the others where the wavelength is within the horizon. The growth of fluctuations which are larger than the horizon can be understood by understanding that such fluctuations behave as if they are in a Universe of their own, but with a different curvature. If the background cosmology is flat, then the background Universe obeys

\begin{equation} \left(\frac{\dot{a}}{a}\right)^2 = \frac{8\pi G \rho}{3} \end{equation}

while the Universe with the fluctuation will have

\begin{equation} \left(\frac{\dot{a_1}}{a_1}\right)^2 = \frac{8\pi G \rho(1+\delta)}{3} + \frac{C'}{a_1^2} \end{equation}

As the fluctuations are small, the differences between \(a_1\) and \(a\) are bound to be small. This implies that the two scale factors closely follow each other when the fluctuations are small. Then we can write,

\begin{eqnarray} \frac{8\pi G \rho(1+\delta)}{3} + \frac{C'}{a_1^2} = \frac{8\pi G \rho}{3} \end{eqnarray}

which allows us to show that

\begin{eqnarray} \delta = \frac{-3C'}{a^2\rho} \end{eqnarray}

In the radiation dominated era, \(\rho \propto a^{-4}\), thus the fluctuations beyond the horizon grow with the scale factor \(a^2\), while in the matter dominated era where \(\rho \propto a^{-3}\), the fluctuation grow as the scale factor.

For the growth of fluctuations with wavelengths within the horizon, one needs to consider the timescales for collapse of these perturbations to the time scale of expansion. The time scale for collapse is given by \((G\rho_{\rm m})^{-1/2}\) while the time scale for expansion is given by \((G\rho_{\rm tot})^{-1/2}\). In the radiation dominated era, the modes within the horizon cannot grow as the expansion time scale is shorter than the collapse time scale. In the matter dominated era the time scales are comparable and the fluctuations are able to grow at the same rate as you see for scales beyond the horizon. This results in a characteristic power spectrum which has a turnover imprinted at the scale of matter radiation equality. This scale has also been seen in observations now. These inhomogeneities then grow as a function of time in order to form structures that we see today.

Borrowed from Tegmark et al 2002.

Beyond this, the growth rates of all modes are proportional to the scale factor if \(\Omega_{\rm m}=1\). In reality there are correction factors that need to account for the presence of any vacuum energy density which reduces the variation of the growth factor with the scale factor (but we will not go into the details of this here).

Thus the recipe to understand the large scale structure of the Universe is simple. We set up tiny density perturbations at the epoch of matter radiation equality. Then we evolve them under the rules of gravity. Dark matter only interacts gravitationally, so it is easiest to simulate the behavior of dark matter with collisionless particles in an n-body simulation.

In the animation below, we can see a simulated patch of the Universe about 80 Mpc on each side, which is filled with a large number of particles. These particles are initially kept on a grid, but then slightly perturbed such that the resultant density perturbations follow the power spectrum that we derived above. We can see the result of the evolution of these density perturbations due to gravitational collapse of these density perturbations. We would like to understand this behaviour.

from IPython.display import Video, HTML, YouTubeVideo
width=600
height=400
YouTubeVideo("eDGtFRj4xXc", width=width, height=height)

Spherical collapse of a top hat perturbation

The collapse of these perturbations looks quite complicated. Therefore let us first study the collapse of a simple spherical tophat perturbation. Let us assume that we are in a flat matter only cosmological model. Let there be a spherical top hat perturbation with radius \(R_0\) in the density field centered at the origin. The density within this perturbation is given by \(\rho\) such that \(\rho/\rho_{\rm crit} = \rho/\bar\rho = \Omega_{\rm m}>1\). We would like to understand the behaviour of this density perturbation with time.

Because of the spherical symmetry, the density perturbation inside \(R_0\) will not be affected by any mass which is outside \(R_0\) and thus will evolve independently as a separate Universe. The radius \(R_0\) of this perturbation will change proportional to the scale factor in this part of the Universe. The scale factor dependence on time is given by the Friedmann equation,

\begin{eqnarray} \frac{1}{a}\frac{da}{dt} = H_0 \left[ \Omega_{\rm m} a^{-3} + (1 - \Omega_{\rm m}) a^{-2} \right]^{1/2} \end{eqnarray}

We have solved this equation before using a code. We saw that the Universe like this reaches a finite maximum scale factor and then the scale factor starts to decrease once again. It turns out that there is actually a parametric solution for this case.

If you change variable such that \(a=A(1-\cos \theta)\) where \(2 A=\Omega_{\rm m}/(\Omega_{\rm m}-1)\). Then substitute in the left hand side of the equation, we will get that

\begin{eqnarray} \sin^{2} \left(\frac{\theta}{2}\right) \frac{d\theta}{dt} = H_0\frac{\left(\Omega_{\rm m}-1\right)^{3/2}}{\Omega_{\rm m}} \end{eqnarray}

Exercise:

• Show that the above equation can be indeed obtained after substitution of \(a=A(1-\cos \theta)\) in the Friedmann equation given above.

We can get an expression for \(t\) in terms of \(\theta\), by integrating the above equation, and expressing \(\sin^2 \theta/2\) in terms of \(\cos \theta\). We obtain \(t = B (\theta - \sin \theta)\) with \(2 B= H_{0}^{-1}\Omega_{\rm m}/(\Omega_{\rm m}-1)^{3/2}\). Thus the parametric solution for the scale factor is given in terms of \(\theta\) as

\begin{eqnarray} a = A (1-\cos\theta)\\ t = B (\theta - \sin\theta) \end{eqnarray}

The maximum is reached when \(\theta=\pi\), when \(a=2A\), \(t=B\pi\), \(R=2 R_0 A\), while \(\theta=2\pi\) corresponds to \(a=0\) and \(t=2B\pi\) when the perturbation collapses and reaches a singularity.

One can use perturbation theory in the very initial stages by perturbing around \(\theta \rightarrow 0\). This allows to work out how the density within the perturbation changes as a function of time. This linear perturbation theory estimate shows that \(\delta \propto a\) in the early part of the evolution (we had worked this out before). We can then ask what would linear theory predict for the overdensity of the perturbation when it collapses in reality. This value of the overdensity predicted from linear theory is called the critical density threshold for collapse. This value of the critical density threshold for collapse based on linear theory turns out to be equal to 1.686.

Let us consider the leading order expansion of the above formulae in terms of \(\theta\). Then

\begin{eqnarray} a \approx A \left( \frac{\theta^2}{2} - \frac{\theta^4}{24} \right) \\ t \approx B \left( \frac{\theta^3}{6} - \frac{\theta^5}{120} \right) \end{eqnarray}

If we consider just the first order then we obtain

\begin{eqnarray} a &\approx& A \frac{\theta^2}{2}\\ t &\approx& B \frac{\theta^3}{6} \\ a &\approx& \frac{A}{2} \left( \frac{6t}{B} \right)^{2/3} \end{eqnarray}

We can iteratively solve for the next order by considering the \(t(\theta)\) expression and noting \begin{eqnarray} t &\approx& B \frac{\theta^3}{6} \left( 1 - \frac{\theta^2}{20} \right) \\ &\approx& B \frac{\theta^3}{6} \left( 1 - \frac{1}{20}\left(\frac{6t}{B}\right)^{2/3} \right)\\ \theta^3 &\approx& \frac{6t}{B}\left( 1 + \frac{1}{20}\left(\frac{6t}{B}\right)^{2/3} \right)\\ \theta &\approx& \left( \frac{6t}{B} \right)^{1/3} \left( 1 + \frac{1}{60}\left(\frac{6t}{B}\right)^{2/3} \right) \end{eqnarray}

Similarly for \(a(t)\), we obtain \begin{eqnarray} a &\approx& \frac{A}{2} \left( \theta^2 - \frac{\theta^4}{12} \right)\\ &\approx& \frac{A}{2} \left(\frac{6t}{B}\right)^{2/3} \left( 1 - \frac{1}{20} \left(\frac{6t}{B}\right)^{2/3} \right) \end{eqnarray}

The density is given by

\begin{eqnarray} \rho &=& \frac{M}{\frac{4}{3}\pi R(t)^3} \\ &=& \frac{M}{\frac{4}{3}\pi R_0^3 a(t)^3} \\ &=& \frac{M}{\frac{4}{3}\pi R_0^3}\frac{8}{A^3}\left(\frac{B}{6t}\right)^2 \left[ 1 + \frac{3}{20} \left(\frac{6t}{B}\right)^{2/3} \right]\\ &=& \bar \rho \left[ 1 + \frac{3}{20} \left(\frac{6t}{B}\right)^{2/3} \right] \end{eqnarray}

Therefore the overdensity contrast evolves with time in the initial stages as:

\begin{eqnarray} \delta &=& \left[\frac{3}{20} \left(\frac{6t}{B}\right)^{2/3} \right] \end{eqnarray}

If we just naively continue this extrapolation of how the overdensity grows, we will of course not be able to predict the overdensity at late times (because we are working in the \(\theta\rightarrow 0\) limit). After all we know that the overdensity collapses at \(t=2B\pi\), when formally the value of the overdensity becomes infinite. But if we naively extrapolate the linear theory prediction, we obtain that the density perturbation collapses when the linear theory prediction gives a value for the evolved perturbation equal to \(\delta_{\rm sc} \approx 1.686\).

The perturbation we considered was a spherical top hat perturbation. Perturbations in the Universe are clearly more complicated, since they do not obey spherical symmetry. Thus in reality the collapse does not proceed to a singularity. The various shells of masses start crossing each other and lead to a virialization process for the halo. Using the virial theorem we can show that the virial radius is proportional to \(R_{\rm max}/2\). At this stage we can show that the true overdensity with respect to the background is close to \(\Delta_{\rm vir} = 18 \pi^2 \approx 178\).

We did all these calculations based on the assumption of a background model with only matter and a flat cosmology. These results can be generalized for more complicated models with curvature and dark energy density. In that case, we will get a different growth rate for perturbations in the linear regime (this factor is called the growth rate factor \(D(z)\)), a slightly different collapse threshold and a different value for the virialization overdensity \(\Delta_{\rm vir}\).

Formation and statistics of dark matter halos

The fluctuations in the dark matter component at the end of radiation domination have a characteristic power spectrum \(P(k)\). We have discussed the statistics of these fluctuations before. These flcutuations grow with time leading to the growth of the power spectrum. For scales with overdensities \(\delta << 1\), the growth factor can be calculated using linear perturbation theory as we did in the last section. Scales with \(\delta \sim 1\), their fluctuation amplitudes grow non-linearly and these scales can collapse in to bound virialized regions as we saw in the simple spherical collapse model.

Now consider a region in a simulation box, evolve its densities with the linear growth factor and compute the average of the overdensity field around it within a sphere of radius \(R\). We can obtain a distribution of such overdensities if you consider different locations of the simulation box. Given that the initial density field is a Gaussian, this smoothed density field will also be a Gaussian with a certain variance. Let us try to compute this variance.

\begin{eqnarray} \delta_{R}(\vec{x}) = \int d^3 \vec{x}' \delta(\vec{x}') W_{R}(\vec{x}-\vec{x}') \end{eqnarray} where the window function \(W_{R}(\vec{x}\) can be defined as

\begin{equation} W_R(\vec{r}) = 0 \,\, \forall |\vec{r}|\gt R \end{equation} \begin{equation} W_R(\vec{r}) = \frac{3}{4\pi R^3} \,\, \forall |\vec{r}| \leq R \end{equation}

The variance of this field is then given by

\begin{eqnarray} \langle \delta_R(\vec{x})^2 \rangle &=& \int d^3 \vec{x'} \int d^3 \vec{x''} \langle \delta(\vec{x}')\delta(\vec{x}'')\rangle W_R(\vec{x}-\vec{x}') W_R(\vec{x}-\vec{x}'') \\ &=& \int d^3 \vec{x'} \int d^3 \vec{x''} \xi(\vec{x}'-\vec{x}'') W_R(\vec{x}-\vec{x}') W_R(\vec{x}-\vec{x}'') \end{eqnarray}

This double convolution can be easily evaluated in Fourier space and then transformed back to obtain

\begin{eqnarray} \sigma_{R}^2 = \langle \delta_R(\vec{x})^2 \rangle &=& {\rm FT} \left[ P(k) W_R(k)^2 \right] \end{eqnarray}

As \(R\rightarrow \infty\), the variance goes to zero, but its behaviour with \(R\) will depend upon the power spectrum. For the power spectrum in the concordance cosmological model, the power spectrum just decreases with \(R\). Now the initial conditions could be evolved with the linear growth rate, and could be smoothed with a large radius \(R\) and ask if there is any point which has a density above the critical density at a given redshift. The radius can be shrunk until there are points which have \(\delta_R(\vec{x})\) greater than the critical density threshold. Those locations will correspond to the particles in the initial density field which lead to the formation of a halo of mass \(M\). Continuing the march downwards would allow us to identify the halos starting from initial conditions.

The fraction of mass bound in halos of mass \(M\) and above can then be calculated as

\begin{equation} F(>M) = \int_{\delta_c}^{\infty} P(\delta|\sigma_{\rm R}) d\delta \end{equation} \begin{equation} F(>M) = \frac{1}{2}{\rm erfc} \left(\frac{\delta_c}{\sqrt{2}\sigma_{\rm R}}\right) \end{equation}

The number density of halos with mass in the range \([M, M+dM]\) is given by

\begin{equation} \frac{d}{dM} n(>M) = \frac{\bar\rho}{M}\frac{d}{dM}F(>M) \end{equation} \begin{equation} \frac{d}{dM} n(>M) \propto \frac{\bar\rho}{M^2} \exp\left(-\frac{\delta_{\rm c}^2}{2\sigma^2}\right)\frac{\delta_c}{\sigma}\frac{d\ln\sigma}{d\ln M} \end{equation} This particular way to derive the halo mass function was first presented by Press and Schechter in 1974. The mass function they found showed that only half of the mass in the Universe would be bound in halos. So they had added a fudge factor of 2 to the mass function they obtained with these arguments. The factor of 2 is already missing in the first equation of the above argument. This was pointed out by Bond et al. 1992. The solution for the statistical problem corresponding to this problem was first presented by Chandrasekhar some time before (he was not solving the mass function problem). It has to do with regions which are underdense than the critical threshold but are part of bigger overdensities (it is called the cloud-in-cloud problem).

The Press Schechter mass function qualitatively reproduces the halo mass function observed in simulations, you see a power law at the low mass end, and an exponential drop off at the massive end. There have been many advances in understanding the halo mass function theoretically with more quantitative success: in particular the use of ellipsoidal collapse rather than spherical collapse to set the critical threshold for collapse (e.g., Sheth, Mo, Tormen 1999, 2001), stochasticity in the critical threshold (e.g., Corasaniti et al 2011), statistics of peaks (Paranjape et al. 2013). However, these approaches often lead to better functional forms whose free parameters are tuned in order to fit the halo mass function in arbitrary cosmologies.

The mass function appears fairly universal, but there are small deviations and these have been now well characterized with the help of numerical simulations. Fitting functions or interpolation routines that can compute the mass function for a given mass and redshift are available. The following figure from Tinker et al. (2008) shows a comparison between the halo mass function as measured from numerical simulations and the Press Schechter mass function.

from IPython.display import Image
width = 600
Image("images/mass_function.png", width=width)