Skip to main content
Physics LibreTexts

31: Structure Formation

  • Page ID
    58989
  • \( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

    Shortly before the epoch of the CMB, matter became the dominant constituent of the universe. Throughout the matter-dominated epoch---i.e., from a few hundred thousand years after the Big Bang until dark energy--matter equality nearly 13 billion years later---dark matter and baryonic structure formed and grew. This Chapter describes how, under the gravitational influence of dark matter, small fluctuations in the matter density field evolve. Particularly dense regions collapse into nonlinear, self-gravitating systems called dark matter halos, which form the nodes of the cosmic that cosmologists observe today.

    Linear Structure Formation

    A key quantity used to describe dark matter structure is the density field \(  \rho(\vec{r},a) \), which specifies the matter density \( \rho \) (i.e., the mass contained in some volume \( V \)) as a function of three-dimensional position \( \vec{r} \) and scale factor \( a \). Specifically, consider the density contrast

    \[ \delta(\vec{r},a) = \frac{\rho(\vec{r},a)-\bar{\rho}(a)}{\bar{\rho}(a)}, \]

    where \( \bar{\rho}(a) \) is the mean matter density at a given epoch. The density contrast describes how matter is distributed relative to the mean density; at a given location and time, \( \delta > 0 \) (\( \delta < 0 \)) corresponds to an overdense (underdense) region, and \( \delta = 0 \) corresponds to a region of average density. We refer to structure on a given scale as linear when \( \lvert \delta \rvert \ll 1 \) and nonlinear when \( \lvert \delta \rvert \approx 1 \). Linear structure formation can largely be modeled analytically, while nonlinear structure formation requires numerical simulations to model accurately.

    How does \( \delta \) evolve over time, and what are its typical values at early times and today? Recall the Friedmann equation from Chapter 11:

    \[ H^2 = \left( \frac{\dot{a}}{a} \right)^2 = \frac{8\pi G \bar{\rho}}{3}, \]

    where we assume that the universe is flat (\( k = 0 \)) and use \(\bar{\rho}\) to denote the density is averaged over a large region of the universe. Technically, the density that enters the Friedmann equation is the sum of matter, radiation, and dark energy, but this sum is dominated by matter during the matter-dominated epoch.

    Consider a region of the universe with some local matter density \( \rho \). That region will either be overdense and behave locally like a closed universe, or underdense and behave locally like an open universe. Thus, this region obeys

    \[ H^2 = \left( \frac{\dot{a}}{a} \right)^2 = \frac{8\pi G \rho}{3} - \frac{k}{a^2} \label{eqn:friedmann_1} \]

    \[ \implies \rho - \bar{\rho} \propto \frac{k}{a^2} \label{eqn:friedmann_2} \]

    \[ \implies \delta = \frac{\rho-\bar{\rho}}{\bar{\rho}} \propto \frac{1}{\bar{\rho} a^2} \propto \frac{a^3}{a^2} \propto a, \label{eqn:friedmann_3} \]

    where Equation \ref{eqn:friedmann_3} follows by subtracting Equations \ref{eqn:friedmann_2} and \ref{eqn:friedmann_1}, and the final line uses \(\bar{\rho} \propto a^{-3}\) as appropriate for non-relativistic matter. For a more rigorous derivation of this scaling, see Problem 31.1.

    As the universe expands, overdense regions with \( \delta > 0 \) become denser relative to the background according to Equation \ref{eqn:friedmann_3}. Gravity attracts matter surrounding a given overdensity towards its center, slowing down the expansion of overdense regions relative to the background and causing their density contrast to grow.

     

    Box \(\PageIndex{1}\)

    Exercise 31.1.1: Qualitatively describe how underdense regions evolve in the linear regime. Does their density contrast increase or decrease as the universe expands? Justify your answer mathematically.

     

    At the time of the CMB, baryon density fluctuations have an average size of

    \[ \delta_{\mathrm{baryon}}(a_{\mathrm{CMB}}) \sim \frac{ \Delta \rho_{\mathrm{baryon}} }{\bar{\rho}_{\mathrm{baryon}}} \sim \frac{ \Delta T_{\mathrm{baryon}} }{\bar{T}_{\mathrm{baryon}}} \approx 10^{-5}. \]

    According to Equation \ref{eqn:friedmann_3}, the average size of these fluctuations today is therefore

    \[ \delta_{\mathrm{baryon}}(a_{\mathrm{today}}) \approx \frac{a_{\mathrm{today}}}{a_{\mathrm{CMB}}} \delta_{\mathrm{baryon}}(a_{\mathrm{CMB}}) \approx 10^{3} \delta_{\mathrm{baryon}}(a_{\mathrm{CMB}}) \approx 10^{-2}. \]

    This is not what we observe! Instead, today's universe is filled with galaxies, which have densities larger than the background by orders of magnitude. In other words, large, non-baryonic density fluctuations must be present at the time of the CMB in order for nonlinear structure to form by today. These density fluctuations are sourced by dark matter. Problem 31.2 explores why baryonic perturbations are much smaller than dark matter perturbations at the time of the CMB.

     

            
    Dark matter drives the growth of small baryonic overdensities observed in the CMB (left) into nonlinear structure traced by galaxies today (right).

     

    Nonlinear Structure Formation

    What happens to overdensities as they near \( \delta \approx 1 \)? Consider a toy model in which a region of the universe has a uniformly higher matter density than the rest,

    \[ \rho(r,t_i) = \bar{\rho}(t_i) \begin{cases} 1+\delta_i, r < R_i \\ 1, r > R_i \end{cases} \label{eqn:tophat} \]

    where \(r\) is the radial distance from the center of the overdensity, \(R_i\) is the initial size of the overdensity, \(t_i\) is the initial time, \(\bar{\rho}(t_i) \) is the background density, and \(\delta_i\) is the size of the overdensity. This is referred to as a top-hat overdensity.

    Problem 31.3 explores the evolution of this system. The overdensity initially expands with the Hubble flow; if \(\delta_i > \Omega_{m}(t_i)^{-1} -1\), gravitational attraction overwhelms the expansion, and the region collapses. This picture schematically describes how self-gravitating dark matter halos form.

     

    Box \(\PageIndex{2}\)

    Exercise 31.2.1: Realistic cosmological density fluctuations are more complex than the top-hat model for several reasons:

    • Overdensities are not spatially uniform; instead, the density contrast is higher near the centers of overdensities and decreases until joining onto the mean density;
    • Fluctuations are not spherically symmetric; instead, collapse occurs successively along the major, intermediate, and minor axes of the initial perturbation;
    • Collapse is not purely radial; particles orbiting an overdensity have angular momentum.

    Qualitatively describe how each effect influences how overdensities evolve. Assuming that the overdensity will collapse, do you expect each effect to speed up or slow down the process?

     

    Gravitational collapse naturally leads to bottom-up (or hierarchical) structure formation. Consider a shell of matter at distance \(r\) from the center of an overdensity as it collapses after turnaround. The time taken for this shell to collapse is given by

    \[ t_{\mathrm{ff}} \sim \sqrt{r/g} \sim \sqrt{r^3/(GM)} \sim \sqrt{1/(G\rho)}, \label{eqn:tff} \]

    where \(M\) is the enclosed mass and \(\rho\) is the enclosed density. This timescale, referred to as the free-fall time, is relevant in many other astrophysical settings where gravity plays an important role.

    Matter fluctuations have higher average densities at early times because the background matter density is higher (even though overdensities are smaller at early times). Thus, according to Equation \ref{eqn:tff}, small overdensities at early times collapse more quickly than large overdensities at later times. These small overdensities collapse, forming small dark matter halos, which merge together to form ever larger structures. The typical masses of these first-generation halos depend on the properties of dark matter and the primordial power spectrum. They are typically at least 10 orders of magnitude smaller than the most massive (\(\approx 10^{15}\ M_{\mathrm{\odot}}\)) clusters today; in many models, including weakly interacting massive particle (WIMP) dark matter, the first halos are even smaller (\(\approx 10^{-6}\ M_{\mathrm{\odot}}\)). In contrast, theories of neutrino dark matter considered in the 1970s predict top-down structure formation, in which "superclusters" form and subsequently fragment into smaller pieces. This scenario is ruled out by observations, which show that small galaxies and halos form first and subsequently merge to form massive objects.

    Numerical simulations are used to predict the distribution and properties of halos over a wide range of mass scales and redshifts. In particular, cosmological N-body simulations solve for the gravitational evolution of discretized particles, which represent the dark matter density field, in an expanding FLRW universe. A key prediction of these simulations is the mass distribution of collapsed halos, or the halo mass function, which is predicted to increase as halo mass decreases in a roughly scale-invariant fashion:

    \[ \frac{dn}{dM} \propto M_{\mathrm{halo}}^{-\alpha}, \]

    where \( n \) is the comoving number density of dark matter halos, \( M \) is halo mass, and \( \alpha \approx 1.9 \).

    Of course, it's difficult to compare predictions for the distribution of dark matter halos directly to observations. Chapter 32 explores how galaxies form in dark matter halos, and shows that the hierarchical structure formation paradigm describes observations of the nonlinear galaxy distribution today extremely well.

     

            
    Growth of dark matter structure in a cosmological N-body simulation. Pink (blue) regions correspond to overdensities (underdensities). Dark matter halos form, grow, and merge within the overdense regions, resulting in the cosmic web of structure cosmologists observe today. 

     

    HOMEWORK Problems

    Problem \(\PageIndex{1}\): The growth of linear dark matter overdensities.

    This problem walks through a more rigorous derivation of the result \(\delta \propto a\). The dynamics of the dark matter density field can be described using the following equations:

    \[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\vec{v}) = 0 \label{eqn:continuity} \]

    \[ \frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla)\vec{v} = -\frac{\nabla p}{\rho} - \nabla \Phi \label{eqn:euler} \]

    \[ \nabla^2 \Phi = 4\pi G \rho, \label{eqn:poisson} \]

    where \(\rho\) is the dark matter density, \(\vec{v}\) is its velocity, \(p\) is its pressure, and \(\Phi\) is the gravitational potential. Equations \ref{eqn:continuity} (the continuity equation) enforces conservation of mass, Equation \ref{eqn:euler} (the Euler equation) enforces conservation of momentum, and Equation \ref{eqn:poisson} (the Poisson equation) describes how matter sources the gravitational potential. Let

    \[ \rho(\vec{r},t) = \bar{\rho}(t) + \Delta\rho(\vec{r},t) = \bar{\rho}(t)[1+\delta(\vec{r},t)] \]

    \[ \vec{v}(\vec{r},t) = \vec{v}_0(\vec{r},t) + \Delta\vec{v}(\vec{r},t) \]

    \[ \Phi(\vec{r},t) = \Phi_0(\vec{r},t)  + \Delta\Phi(\vec{r},t). \]

    By plugging these into the fluid equations and keeping terms up to first order in the fluctuations \(\delta,\ \Delta\vec{v},\) and \(\Delta\Phi\), combining the equations yields a differential equation for the evolution of the density contrast:

    \[ \ddot{\hat{\delta}} + 2H\dot{\hat{\delta}} = \left( 4\pi G \bar{\rho} - \frac{c_s^2 k^2}{a^2} \right)\hat{\delta}, \label{eqn:delta_diffeq_cs} \]

    where \(\hat{\delta}\) is the Fourier transform of delta, \(k \) is the wavenumber associated with the Fourier transform, and \( c_s^2 =\partial p/\partial\rho \) is the speed of sound. On large scales (small \(k\)), dark matter behaves as a pressureless fluid with \(c_s=0\). Thus,

    \[ \ddot{\hat{\delta}} + 2H\dot{\hat{\delta}} =\frac{3}{2}H^2\hat{\delta}. \label{eqn:delta_diffeq} \]

    a) Let \(\hat{\delta} \propto t^n\) and solve for \(n\) by plugging this into Equation \ref{eqn:delta_diffeq}; you'll obtain two solutions for \(n\) because this is a second-order differential equation. Assume that the universe only contains matter to plug in \(H(t)\).

    b) Interpret both solutions obtained in part a) physically. Show that one of the solutions yields the expected behavior \(\delta \propto a\).

    Problem \(\PageIndex{2}\): The growth of linear baryonic overdensities.

    We argued that dark matter overdensities at the time of the CMB are large enough to facilitate nonlinear structure formation by today. But why are baryonic fluctuations smaller than dark matter fluctuations at early times? The key difference is that baryons can exchange energy and momentum by interacting with each other, and are therefore not a pressureless fluid. Thus, for baryons, we can't drop the \(c_s^2\) term in Equation \ref{eqn:delta_diffeq_cs}.

    The kinds of solutions to Equation \ref{eqn:delta_diffeq_cs} are determined by the sign of \(4\pi G \bar{\rho} - \frac{c_s^2 k^2}{a^2}\). At a given time, this depends on whether \(k\) is smaller or larger than a critical wavenumber \(k_J\); the corresponding scale \(\lambda_J = 2\pi/k_J \) is referred to as the Jeans length.

    a) Derive the Jeans length in terms of \(c_s\), \(a\), and \(\bar{\rho}\).

    b) Qualitatively describe how baryonic overdensities evolve on scales larger than and smaller than the Jeans length.

    c) Before recombination, baryons feel pressure exerted by photons; the corresponding speed of sound is \(c_s \approx c/\sqrt{3}\). Calculate the corresponding Jeans length assuming recombination occurs at \(z = 1100\). Relate this scale to the wavenumber of the peak in the linear matter power spectrum shown in the figure below (adapted from Chabanier et al. 2019).

    Problem \(\PageIndex{3}\): Top-hat collapse.

    Consider the top-hat perturbation defined in Equation \ref{eqn:tophat} A shell at radius \(r_i\) initially expands with the Hubble Flow, \(v_i = H(t_i)r_i\), and thus has an initial energy per unit mass (or specific energy) of

    \[E_i = K_i + U_i = \frac{1}{2}v_i^2 - \frac{GM(<r_i)}{r_i},\]

    where \(M(<r_i)\) is the enclosed mass. The integrated equation of motion for this shell is

    \[ \frac{1}{2}\left(\frac{\mathrm{d}r}{\mathrm{d}t} \right)^2 - \frac{GM}{r} = E.\]

    a) Show that the system collapses (corresponding to \(E_i < 0)\) if \(\delta_i > \Omega_m(t_i)^{-1} - 1\).

    b) Show that \(r(\theta) = A(1-\cos\theta)\), \(t(\theta) = B(\theta-\sin\theta)\) solves the integrated equation of motion, where \(\theta \in [0,2\pi]\), \(A=GM/2\lvert E\rvert \), and \(B = GM/(2\lvert E \rvert)^{3/2}\).

    c) Run the code snippet and qualitatively describe the plotted evolution of the system. Is this qualitatively consistent with the description near Equation \ref{eqn:tophat}?

    import numpy as np
    import matplotlib.pyplot as plt
    %matplotlib inline
    %config InlineBackend.figure_format='retina'
    
    ###
    
    #parametric solution in GM = 1 units
    theta_array = np.linspace(1e-1,2*np.pi-1e-1,100)
    
    def r_tophat(theta,E):
        A = 1./(2.*np.abs(E))
        return A*(1.-np.cos(theta))
    
    def t_tophat(theta,E):
        B = 1./(2.*np.abs(E))**(1.5)
        return B*(theta-np.sin(theta))
    
    ###
    
    #Plot evolution of system
    plt.figure(figsize=(8,6))
    E_array = np.linspace(-0.5,-0.1,5)
    
    for E in E_array:
        plt.plot(t_tophat(theta_array,E),r_tophat(theta_array,E),label=r'$E=${:0.1f}'.format(E))
        
    plt.xlabel(r'$t$',fontsize=16)
    plt.ylabel(r'$r$',fontsize=16)
    plt.legend(loc=1,fontsize=15)
        
    plt.show()

    This page titled 31: Structure Formation is shared under a CC BY 4.0 license and was authored, remixed, and/or curated by Ethan Nadler.

    • Was this article helpful?