Formation of Large-scale Structures in the Universe

Until now, we have studied the evolution of a homogeneous Universe, at least at very large scales, beyond approximately $100\,\Mpc$ . However, today we observe that matter is agglomerated in the form of planets, stars, galaxies, galaxy clusters and superclusters of galaxies. The question that arises in this chapter is: how do these structures form in an expanding Universe? If the formation of the smallest objects involves many physical processes and is highly dependent on local initial conditions, it is possible to develop a simple linear model for the evolution of the largest structures $\sim 50\,\Mpc$ , such as galaxy clusters or superclusters. To do this, we will simply use a Newtonian theory of linear perturbations, and calculate the evolution of these perturbations in an expanding Universe.

Figure 1:Size and mass scales in the Universe.

1Range of Validity¶

However, let us state the necessary conditions for the validity of the Newtonian perturbation theory. A first obvious condition is that the perturbation of the metric $g_{\mu\nu}$ and of the energy-momentum tensor $T^{\mu\nu}$ must be small. In the case of non-relativistic matter, this means that density perturbations $\delta \rho$ must remain small compared to the average density $\rho_0$ i.e. $\delta \rho/\rho_0 \ll 1$ . Perturbations of the gravitational field $\delta \phi$ are also small compared to the average gravitational field $\phi_0$ i.e. $\delta \phi / \phi_0 \ll 1$ .

But to be able to return to Newtonian gravity in a homogeneous expanding spacetime and study only density perturbations, there is another condition. The size of the studied regions must be much smaller than the Hubble radius $D_H={c / H}$ , because in Newtonian gravitation the interaction propagates instantaneously. If the structure is of cosmological size, by the time the information from a collapsing region has traveled to its boundaries, the expansion factor changes significantly: a situation that obviously cannot be treated using Newtonian gravity. We therefore only consider structures of size much smaller than $D_H$ to be able to assume that the gravitational force propagates instantaneously in the structure with respect to the expansion rate of the Universe.

2Classical Acoustics Preamble¶

The non-relativistic matter field is described by a perfect fluid model. Looking at the evolution of perturbations in a perfect fluid means studying the field of acoustics. Let us recall some useful notions here.

In an Eulerian description of fluids, for a perfect fluid at rest, we define its velocity field $\vec v(t, \vec r)$ , its pressure field $P(t, \vec r)$ and its mass density field $\rho(t, \vec r)$ as constants:

\vec v(t, \vec r)=\vec 0, P(t, \vec r) = p_0, \rho(t, \vec r) = \rho_0

(1)

On top of this, we superimpose small amplitude perturbations:

\begin{align*} \vec v(t, \vec r) & = \vec v_0 + \delta \vec v(t, \vec r) \\ P(t, \vec r) & = P_0 + \delta P(t, \vec r) \\ \rho(t, \vec r) & = \rho_0 + \delta \rho(t, \vec r) \\ \end{align*}

(2)

The generated overpressures change the internal energy of the perturbations and thus their temperature, and this work of pressure forces allows the propagation of acoustic energy. In acoustics, we measure that compressions are adiabatic and reversible, meaning that the size of perturbations is larger than the mean free path of fluid particles: we can then neglect heat conduction. We therefore use the isentropic compressibility $\chi_S$ instead of the isothermal compressibility $\chi_T$ :

\chi_S = -\frac{1}{V} \left.\frac{\partial V}{\partial P}\right\vert_S \approx \frac{1}{\rho_0} \frac{\delta \rho}{\delta P}

(3)

We deduce the link between pressure and density perturbations via the isentropic compressibility:

\delta \rho = \rho_0 \chi_S \delta P

(4)

We also have local mass conservation:

\frac{\partial \rho}{\partial t} + \diverg \left(\rho \vec v\right) = 0 \Rightarrow \frac{\partial \delta \rho}{\partial t} + \rho_0 \diverg \delta \vec v = 0

(5)

Since the fluid is perfect, its dynamics is governed by the Euler equation (Navier-Stokes without viscosity):

\frac{\partial \vec v}{\partial t} + \left(\vec v \cdot \grad\right)\vec v = - \frac{1}{\rho} \grad P + \cancel{\vec g} \Rightarrow \frac{\partial \delta \vec v}{\partial t} = - \frac{1}{\rho_0} \grad \delta P

(6)

where we usually neglect gravitation. By combining these three equations, we obtain a d’Alembert equation:

\boxed{ \frac{\partial^2 \delta \vec \rho}{\partial t^2} - \frac{1}{\rho_0 \chi_S} \triangle \delta \rho = 0, \quad c_s = \frac{1}{\sqrt{\rho_0\chi_S}} }

(7)

with $c_s$ the speed of sound.

For a perfect gas and isentropic transformations, we use Laplace’s Law $PV^\gamma = \cst$ with $\gamma=C_P / C_V$ the ratio of heat capacities at constant pressure $C_P$ and constant volume $C_V$ , also called the adiabatic index. If the fluid particles possess $N_{dof}$ degrees of freedom (translations, rotations, vibrations) to store energy in kinetic (thermal) form, then the adiabatic index is $\gamma = 1+2/N_{dof}$ . For a diatomic gas like air at usual temperature, there are $N_{dof}=5$ degrees of freedom (3 translations and 2 rotations) so $\gamma = 7/5$ . We deduce a formula for the speed of sound as a function of temperature and particle mass^[1]:

PV^\gamma = \cst \Rightarrow \chi_S \approx \frac{1}{\gamma P_0}

(8)

\boxed{ c_s = \sqrt{\frac{\gamma P_0}{\rho_0}} = \sqrt{\frac{\gamma RT}{M}} = \sqrt{\frac{\gamma k_B T}{m_{\text{air}}}} }

(9)

with $M$ the molar mass of the fluid and $m_{\text{air}}$ the effective mass of an air particle if that is the medium considered, calculated by:

m_{\text{air}} = 78\%\,m_{N_2} + 22\%\,m_{O_2}

(10)

3Acoustics in Expanding Space with Gravitation¶

At the emission of the CMB, we observe that the Universe is a perfect gas of non-relativistic matter, slightly inhomogeneous, with density enhancements of reduced size (of the order of $\delta \rho / \rho_0 \sim 10^{-5}$ ). These fluctuations will a priori increase over time thanks to gravitational attraction, to result in a Universe where the density contrast is of the order of $\delta \rho / \rho_0 \sim 1$ ^[2] ( Figure 1,Figure 2). What are the mechanisms allowing this growth?

The formation of large-scale structures in the Universe can be understood through a linear perturbation theory in the perfect fluid constituted by non-relativistic matter after recombination. To study the evolution of these density perturbations in the matter field, it is therefore a matter of rewriting the previous acoustic equations, in the particular context of an expanding medium subject to the laws of gravity.

Figure 2:An artistic celebration of the Dark Energy Spectroscopic Instrument (DESI) Year 1 data, showing a slice of the larger 3D map that DESI is constructing during its five-year survey. DESI is mounted on the Nicholas U. Mayall 4-meter Telescope at Kitt Peak National Observatory. Credit: DESI Collaboration/KPNO/NOIRLab/NSF/AURA/P. Horálek/R. Proctor

Figure 3:This simulation shows how gravity affects the position of observed galaxies, thus modifying the way matter clumps together to form cosmic structures. As different gravity models predict different structure formations, DESI scientists can compare observations with predictions and thus test gravity at cosmological scales. Credit: Claire Lamman and Michael Rashkovetskyi / DESI collaboration

Coordinate Systems¶

The equations of classical physics are valid using physical (non-comoving) quantities: proper distance, proper time, etc... The use of physical distances as a coordinate system is however impractical in an expanding Universe because time and distance are linked. A more practical method is to write these equations in comoving space. Let us denote:

$\vec{r}$ the physical position vector (whose norm is the proper distance $D_p(t)$ ),
$\vec{\sigma}$ the comoving position vector (whose norm is the comoving distance $\chi$ ),
$\vec{v}_r$ the physical velocity,
$\vec{v}_{\sigma}$ the comoving velocity.

We have the following relations:

\begin{align} \vec{r}& =a(t) \vec{\sigma} \\ \notag \vec{v}_r& =\dot{a}\vec{\sigma}+a \vec{v}_{\sigma} \end{align}

(11)

Although the Hubble flow (first term of (11)) should not be considered as a real velocity, it appears in the relationship between $\vec{v}_r$ and $\vec{v}_{\sigma}$ . This would pose a problem if it could take values of the order of $c$ (or larger). But the condition that $\vec{r} \ll D_H$ guarantees that $\dot{a}\vec{\sigma} \ll c$ (see section Section 1).

How do we switch from the physical coordinate system to the comoving coordinate system? The conversion of spatial partial derivatives is quite simple:

\nabla_\mathbf{r} = {\partial \over \partial r^i} = {1 \over a} {\partial \over \partial \sigma^i} = \frac{1}{a} \nabla_\sigma,\qquad \nabla^2_\mathbf{r} = \frac{1}{a^2} \nabla^2_{\sigma}

(12)

The time partial derivative requires more careful examination. Let us write the differential of any function of time and space $f(t,\vec{r})$ :

\begin{align} \dd f&= \left. {\partial f\over \partial t} \right|_\mathbf{r} \dd t + \left.{\partial f\over \partial r^i}\right|_\mathbf{t} \dd r^i \\ &= \left. {\partial f\over \partial t} \right|_\mathbf{r} \dd t + \left.{\partial f\over \partial r^i}\right|_\mathbf{t} \dd(a \sigma^i) \\ &= \left. {\partial f\over \partial t} \right|_\mathbf{r} \dd t + {1 \over a}\left.{\partial f\over \partial \sigma^i}\right|_\mathbf{t} (\sigma^i \dd a + a \dd \sigma^i)\\ &= \left( \left. {\partial f\over \partial t} \right|_\mathbf{r} + {\dot{a} \over a} (\vec{\sigma}\cdot \vec \nabla_{\sigma})f \right) \dd t + \left.{\partial f\over \partial \sigma^i}\right|_\mathbf{t} \dd \sigma^i \end{align}

(13)

Therefore by identification:

\left. {\partial \over \partial t} \right|_\mathbf{r}= \left. {\partial \over \partial t} \right|_{\sigma} - H (\vec{\sigma}\cdot \vec \nabla_{\sigma})

(14)

Mass Conservation¶

In physical space, the continuity equation is written:

\left.{\partial \rho \over \partial t}\right|_\mathbf{r} + \vec \nabla_\mathbf{r} \cdot (\rho \vec{v}_\mathbf{r})=0

(15)

Using the comoving density $\rho_{\sigma}=\rho a^3$ , we can express this equation in comoving space:

\begin{align} \left.{\partial \rho_{\sigma}\, a^{-3} \over \partial t}\right|_{\sigma} - H (\vec{\sigma}.\nabla_{\sigma})(\rho_{\sigma} \, a^{-3})+{1 \over a} \nabla_{\sigma}\left[\rho_{\sigma} \, a^{-3} (\dot{a}\vec{\sigma}+a \vec{v}_{\sigma})\right] &=&0 \end{align}

(16)

By introducing the classical vector calculus formula

\vec \nabla \cdot (\phi \vec{A})=\vec{A}\cdot \vec \nabla\phi+ \phi\ \vec \nabla \cdot \vec{A}

(17)

where $\phi$ is a scalar and $\vec{A}$ a vector, the above equation reduces to:

\boxed{\left.{\partial \rho_{\sigma}\over \partial t}\right|_{\sigma} + \vec \nabla_{\sigma} \cdot (\rho_{\sigma} \vec{v}_{\sigma})=0}

(18)

The continuity equation therefore takes exactly the same form in comoving space, using comoving variables. However, this is not the case if we use the peculiar velocity $a \vec{v}_{\sigma}$ .

Euler Equation¶

In physical space, the Euler equation for a self-gravitating perfect (non-viscous) fluid takes the usual form:

\left. {\partial \vec{v_r} \over \partial t} \right|_\mathbf{r}+ \left(\vec{v}_r\cdot \vec \nabla_\mathbf{r}\right) \vec{v_r}=-\vec \nabla_\mathbf{r}\phi - \frac{1}{\rho} \vec \nabla_\mathbf{r} P

(19)

In comoving space, it transforms into:

\begin{align*} \left. {\partial (\dot{a}\vec{\sigma}+a\vec{v}_{\sigma}) \over \partial t} \right|_{\sigma} & - H \vec{\sigma}\cdot \vec \nabla_{\sigma}(\dot{a}\vec{\sigma} + a\vec{v}_{\sigma}) + (\dot{a}\vec{\sigma}+a\vec{v}_{\sigma}) \cdot {1 \over a}\vec \nabla_{\sigma} (\dot{a}\vec{\sigma}+a\vec{v}_{\sigma}) \\ & = -{1 \over a^2} \vec \nabla_{\sigma}\phi_{\sigma} - \frac{1}{a \rho} \vec \nabla_{\sigma} P \end{align*}

(20)

with $\phi_\sigma = a \phi$ the comoving potential. Ensuring that:

\left. {\partial (\dot{a}\vec{\sigma}+a\vec{v}_{\sigma}) \over \partial t} \right|_{\sigma}=\ddot{a}\vec{\sigma}+\dot{a}\vec{v}_{\sigma}+a \left.{\partial \vec{v}_{\sigma} \over \partial t}\right|_{\sigma}

(21)

the Euler equation in comoving coordinates reduces to:

\boxed{\left.{\partial \vec{v}_{\sigma} \over \partial t} \right|_{\sigma} + 2 H \vec{v}_{\sigma}+(\vec{v}_{\sigma} \cdot \vec \nabla_{\sigma})\vec{v}_{\sigma}= -{1 \over a^3} \vec \nabla_{\sigma}\phi_{\sigma} - \frac{1}{a^2 \rho} \vec\nabla_{\sigma}P - {\ddot{a} \over a} \vec{\sigma} }

(22)

We can see two additional terms compared to the version of the Euler equation written in physical coordinates. The additional term $2 H \vec{v}_{\sigma}$ is a drag force created by the expansion on comoving velocities. The last term depends only on kinematic quantities, so it is a fictitious force linked to the change of coordinate system, non-zero if and only if there is acceleration between the reference frames ( $\ddot{a}=0$ ). In an accelerating expanding universe ( $\ddot{a}>0$ ), it corresponds to a centripetal fictitious force.

Poisson Equation¶

In a universe without cosmological constant, from general relativity we deduce the usual Poisson equation in the weak field limit (38):

\triangle_\mathbf{r} \phi = 4 \pi \GN \rho

(23)

where $\phi$ is the gravitational potential. Given the relationship between the potential and the Newtonian force, the comoving potential is naturally defined as $\phi_{\sigma}=\phi a$ , and the Poisson equation is also unchanged in comoving space:

\triangle_{\sigma} \phi_{\sigma} = 4 \pi \GN \rho_{\sigma}

(24)

4Newtonian Theory of Perturbations¶

Zeroth Order Solution¶

In a homogeneous and isotropic universe, at zeroth order the mass density does not depend on space and matter is globally at rest. There is therefore no global peculiar velocity and the mass conservation equation in comoving space (18) admits the zeroth order solution:

\rho_{{\sigma},0}= \cst, \qquad \vec{v}_{{\sigma},0}= \vec 0

(27)

The comoving mass density therefore remains stationary and homogeneous. The integration of the Poisson equation (24) in spherical coordinates^[3] gives:

\phi_{\vec{\sigma},0} = \frac{4}{6} \pi \GN \rho_{\sigma,0} \sigma^2

(28)

If we use the Friedmann equations (60):

\frac{\ddot a}{a} = -\frac{4 \pi \GN}{3}(\epsilon + 3P)

(29)

for a matter-dominated Universe with $\Lambda=0$ and $\epsilon = c^2 \rho_{\sigma,0} a^{-3} \gg P$ , we obtain the force derived from the gravitational potential:

-\vec \nabla \phi_{{\sigma},0} = -{4 \pi \GN \rho_{\sigma,0} \over 3 } \vec{\sigma} \approx a^3 \frac{\ddot a}{a} \vec{\sigma}

(30)

The gradient of the gravitational field gives at zeroth order a centrifugal force reflecting the repulsive character of the universe’s expansion.

From the Euler equation (22) and the previous results, we logically arrive at finding that $\vec \nabla_\sigma P =\vec 0$ therefore:

P = \cst = P_0

(31)

consistent with a homogeneous Universe.

First Order Linear Solution¶

Let us consider small perturbations around the zeroth order solution:

\rho_{\sigma}=\rho_0(1+\delta) \qquad \vec{v}_{\sigma}=\vec{v}_0+\delta\vec{v} \qquad \phi_{\sigma}=\phi_0 + \delta \phi_{\sigma}, \qquad P = P_0 + \delta P

(32)

The quantity $\delta=\delta \rho/\rho$ is called the density contrast and is widely used in structure formation theory. Note that it has the same value in physical space and in comoving space. By injecting these expressions into the continuity, Poisson and Euler equations, the zeroth order solution disappears (as expected) and by removing second order terms, we obtain the set of linearized equations:

{\partial \delta \over \partial t} + \vec{\nabla}_{\sigma} \cdot \delta\vec{v}_r = 0

(33)

{\partial \delta\vec{v}_\sigma \over \partial t} + 2 H \delta\vec{v}_\sigma = - {1 \over a^3} \vec \nabla_{\sigma} \delta \phi_{\sigma} - \frac{1}{a^2 \rho_0} \vec \nabla_{\sigma} \delta P

(34)

\triangle_{\sigma}^2 \delta \phi_{\sigma} = 4 \pi \GN \rho_{\sigma,0} \delta

(35)

By taking the divergence of the linearized Euler equation (34), the time derivative of the linearized mass conservation equation (33) and reinjecting the Poisson equation (35), we finally obtain:

\ddot{\delta}+ 2 H \dot{\delta} - \frac{1}{a^2 \rho_0} \triangle_{\sigma} \delta P = 4 \pi \GN {\rho_0 } \delta

(36)

We introduce the speed of sound $c_s$ in the medium for adiabatic compressions:

c_s^2 = \left.\frac{\partial P}{\partial \rho}\right\vert_{S} \approx \frac{\delta P}{\delta \rho}\Rightarrow \delta P \approx c_s^2 \rho_0 \delta

(37)

Then we obtain the perturbation growth equation in comoving space:

\boxed{ \ddot{\delta}+ 2 H \dot{\delta} - \frac{c_s^2}{a^2} \triangle_{\sigma} \delta = 4 \pi \GN \rho_0 \delta }

(38)

Let us compare this equation to the d’Alembert equation obtained in classical acoustics (7). We have two additional terms:

a friction term $2 H \dot{\delta}$ which will slow down structure formation all the more as the expansion rate is high, on cosmological time scales $\sim 1/H$ ;
a source term proportional to $\delta$ which models gravitational collapse (the more marked the structure, the more it collapses). The evolution of matter perturbations is therefore dictated by the relative intensities of the gravitational force (attractive), pressure forces (repulsive), and the expansion of the Universe (damping).

Speed of Sound¶

In a radiation-dominated Universe, the properties of the relativistic gas are:

\rho_r \propto T^4, \quad P_r = \rho_r c^2 /3

(39)

hence the speed of sound in the relativistic plasma:

c_s^2 = \left.\frac{\partial P_r}{\partial \rho_r}\right\vert_{S} \approx \frac{c^2 \delta \rho_r}{3\delta \rho_r} \Rightarrow \boxed{c_s = \frac{c}{\sqrt{3}} \approx 0.58 c}

(40)

Just before recombination, the propagation of sound waves is therefore relativistic in the primordial plasma.

In a non-relativistic matter-dominated Universe, the thermodynamic properties of baryons are:

\rho_b = n_b m_b + \frac{n_b k_B T}{(\gamma_b -1)c^2},\quad P_b = n_b k_B T

(41)

with $m_b=(1-Y_p)m_\mathrm{H} + Y_p m_{\mathrm{He}} \approx 1.72 m_p$ the effective mass of baryons in a mixture of atomic hydrogen and helium gas (note the resemblance with (10)), and $\gamma_b=5/3$ the adiabatic index of baryons because hydrogen is in atomic form (and helium of course). For a reversible adiabatic transformation, with Laplace’s law we have $P_b \rho_b^{-\gamma_b}=\cst$ , the baryon temperature $T_b$ evolves as:

P_b = n_b k_B T_b \propto a^{-3\gamma_b}

(42)

But $n_b \propto a^{-3}$ so the baryon temperature evolves as:

T_b = T_{\mathrm{dec}} \left(\frac{a}{a_{\mathrm{dec}}}\right)^{-3(\gamma_b -1)}

(43)

With $\gamma=5/3$ , for baryons the pressure varies as $a^{-5}$ and the temperature as $a^{-2}$ . The speed of sound in the baryon gas is written:

c_s^2 = \left.\frac{\partial P_b}{\partial \rho_b}\right\vert_{S} \approx \gamma_b \rho^{\gamma_b-1} = \frac{\gamma_b n_b k_B T_b}{n_b m_b} \Rightarrow \boxed{c_s(T_b) = \sqrt{\frac{\gamma_b k_B T_b}{m_b}} \ll c}

(44)

Just after photon-electron decoupling, the baryon temperature is $T_{\mathrm{dec}} = 3055\,\kelvin$ so the speed of sound in the baryon gas is approximately $1.5 \times 10^{-5}c \approx 5\,\mathrm{km}/\mathrm{s}$ . Around recombination, the change is of several orders of magnitude so the evolution of perturbations changes drastically at this moment. Then the baryon temperature decreases as $a^{-2}$ , at least as long as heating by interstellar radiation does not change things (around $z\sim 10$ ). The speed of sound therefore decreases as $a^{-1}$ in our linear model.

<Figure size 640x480 with 2 Axes> — Figure 4:Speed of sound in the Universe as a function of the scale factor $a$ . The slight decrease between equivalence and decoupling is given by equation (48) (see box). The calculation is inaccurate when we enter the non-linear regime of structure growth i.e. when $\delta \sim 1$ (shaded region). Moreover, around this redshift also, the baryon gas is heated by stellar radiation so the baryon temperature no longer evolves as $a^{-2}$ .

Mixture of Baryons and Coupled Photons

As long as photons are coupled to free electrons, they are mechanically coupled to the baryon gas by the Coulomb force. During the period when the Universe is dominated by matter but photons are still coupled to electrons (Weinberg (1989) p.509,566):

\rho_{\mathrm{tot}} = \rho_b + \frac{n_b k_B T}{(\gamma_b -1)c^2} + \rho_\gamma \approx \eta m_b n_\gamma + \rho_\gamma

(45)

P_{\mathrm{tot}} = n_b k_B T + \frac{1}{3}\epsilon_\gamma \approx \frac{1}{3}\rho_\gamma c^2

(46)

\begin{align*} c_s^2 & = \left.\frac{\partial P_{\mathrm{tot}}}{\partial \rho_{\mathrm{tot}}}\right\vert_{S} \approx \frac{\delta P_{\mathrm{tot}}}{\delta \rho_{\mathrm{tot}}} \approx \frac{\frac{4}{3}\frac{\delta T}{T}\epsilon_\gamma}{3 \eta m_b n_\gamma \frac{\delta T}{T} + 4 \frac{\delta T}{T}\rho_\gamma} \\ & = \frac{c^2}{3} \frac{1}{\left(1 + \frac{3}{4} \frac{\eta n_\gamma m_b c^2}{\epsilon_\gamma} \right)} = \frac{c^2}{3} \frac{1}{\left(1 + \frac{3}{4} \frac{\rho_b}{\rho_\gamma} \right)} \end{align*}

(47)

c_s = \frac{c}{\sqrt{3}} \frac{1}{\sqrt{\left(1 + \dfrac{3}{4} \dfrac{\Omega_b^0}{\Omega_\gamma^0(1+z)} \right)} }

(48)

Jeans Instability¶

The solutions to the perturbation evolution equation (38) will be studied in the next semester, but we can already study their behavior from their dispersion relation. Let us look for plane wave solutions of the form:

\delta \propto e^{i(\omega t - \vec k_\sigma \cdot \vec \sigma)}

(49)

and study the solutions in Fourier domain. The dispersion relation is written:

-\omega^2 + 2 i H \omega - \frac{c_s^2}{a^2} k^2_\sigma = 4 \pi \GN \rho_0

(50)

By neglecting the damping term (of the order of the age of the Universe) compared to the characteristic time of perturbation evolution ( $\omega \gg H$ ), we arrive at:

\boxed{\omega^2 = c_s^2 \left( k_r^2 - \frac{4\pi \GN \rho_0}{c^2_s}\right) }

(51)

with $k_\sigma = 2 \pi / \chi = a k_r$ . This is the same dispersion relation as that of electromagnetic waves in a plasma. We define the Jeans wavenumber (The Royal Society (1902)) $k_J$ and the Jeans length $\lambda_J$ by:

\boxed{k_J = \sqrt{\frac{4\pi \GN \rho_0}{c^2_s}},\quad \lambda_J = \frac{2\pi}{k_J} = c_s \sqrt{\frac{\pi}{\GN \rho_0}} }

(52)

If $k_r > k_J \Leftrightarrow \lambda < \lambda_J$ , then $\omega^2 > 0$ so we have an oscillating solution i.e. an acoustic wave that propagates: perturbations of size small compared to the Jeans length oscillate thanks to the pressure force and do not grow.
If $k_r < k_J \Leftrightarrow \lambda > \lambda_J$ , $\omega^2 < 0$ so we have a non-oscillating exponential solution: large structures evolve only under the effect of gravity and grow (and voids decrease).

Since the scale $\lambda$ evolves with expansion ( $\lambda \propto a$ ) as well as the Jeans length, the discussion is not easy to conduct along the history of the Universe. Let us define $M=4\pi \rho_0\lambda^3/3$ the mass in a sphere of radius $\lambda$ , conserved with expansion, and compare it to the Jeans mass $M_J$ :

M_J = \frac{4}{3} \pi \rho_0 \lambda_J^3 \propto \frac{c_s^3}{\GN^{3/2} \rho_0^{1/2}}

(53)

We now have only one quantity that varies with $a$ :

If $M < M_J$ , then the structure is too light and pressure can compensate for the force of gravity: the structure oscillates.
If $M > M_J$ , then the structure is too heavy and collapses gravitationally.

The plot of the Jeans mass as a function of the scale factor $a$ allows us to predict which structures can collapse gravitationally and those whose growth is prevented Figure 5 (Weinberg (1989) p. 565).

Describing structure growth before recombination contradicts the non-relativistic hypothesis on which our study is based. However, the speed of sound scale remains a relevant scale and moreover non-relativistic matter dominates after equivalence. Knowing this (and that the relativistic study does not give fundamentally different results from what follows), we will allow ourselves to discuss the growth of large structures in the primordial and recent Universe (although it is non-linear today).

According to Figure 5, we see that before decoupling structures with the mass of a galaxy or galaxy cluster are not heavy enough to collapse (or only in the very first moments of the Universe). Pressure waves travel through them and they oscillate. However, after decoupling, the speed of sound drops by 5 orders of magnitude. The decoupling of photons freezes the pressure waves and structures with the mass of galaxies and even dwarf galaxies can begin their growth.

Statistical Description of Perturbations¶

Gaussian Field¶

Let us describe the field $\delta(\vec r)$ in Fourier space. We use the following convention for direct and inverse Fourier transforms:

\hat{\delta}(\vec{k})= \int \delta(\vec{r}) e^{i \vec k \cdot \vec r} \dd^3\vec r

(54)

\delta(\vec{r})= {1 \over (2 \pi)^3}\int \hat{\delta}(\vec{k}) e^{-i \vec k \cdot \vec r} \dd^3\vec k

(55)

The field in Fourier is a complex number described by its modulus and phase which we will denote:

\hat \delta(\vec k) = \vert \delta_k \vert e ^{i \phi_k}

(56)

At the moment of the Big Bang, we can think that the phases are randomly and independently distributed (no coupling between scales $\vec k$ , see chapter ). If perturbations follow a linear evolution, then the modes evolve independently of each other: the phases remain independent random variables. But if the phases of a field in Fourier are random and independent, then this field is called a Gaussian field. If we examine the inverse Fourier transform, we note indeed that the field in real space is the result of the sum of independent random variables so we can apply the central limit theorem: the result of the sum follows a Gaussian distribution. This field is therefore statistically entirely defined by the two parameters of the Gaussian distribution, its mean and its variance:

\left\langle \delta \right\rangle = 0, \quad \sigma^2 = \left\langle \delta(\vec r)\delta(\vec r) \right\rangle = \int_{\vec r} \delta(\vec r)\delta^*(\vec r) \dd^3 \vec r

(57)

with $\delta(\vec r) = \delta^*(\vec r)$ since this field is physical hence real ( $\delta(\vec r) \in \mathbb{R}$ ). Strictly speaking, the notation $\left\langle \ldots \right\rangle$ designates an average over a statistical ensemble of realizations. However, since we only have one realization of the universe, this will be an average over a volume with samples of $\delta$ . The mean is zero by definition of the density contrast, so let us focus on the variance.

Matter Power Spectrum¶

The conjugate of a Fourier mode is written:

\delta^*(\vec{r})= {1 \over (2 \pi)^3}\int \hat{\delta}^*(\vec{k}) e^{i \vec k \cdot \vec r} \dd^3\vec k

so the variance can be rewritten:

\begin{align*} \sigma^2 & = \frac{1}{(2\pi)^3} \int_{\vec r} \delta(\vec r) \left[\int_{\vec k} \delta^*(\vec{k}) e^{i \vec k \cdot \vec r} \dd^3\vec k\right] \dd^3\vec r \\ & = \frac{1}{(2\pi)^3} \int_{\vec k} \hat \delta^*(\vec k) \left[\int_{\vec r} \delta(\vec{r}) e^{i \vec k \cdot \vec r} \dd^3\vec r \right] \dd^3\vec k\\ & = \frac{1}{(2\pi)^3} \int_{\vec k} \hat \delta^*(\vec k) \hat\delta(\vec k) \dd^3\vec k = \frac{1}{(2\pi)^3} \int_{\vec k} \vert \hat\delta(\vec k)\vert^2 \dd^3\vec k\\ \end{align*}

(58)

For an isotropic field, only the norm of the vector $\vec k$ matters therefore:

\sigma^2 = \frac{1}{2\pi^2} \int_{0}^\infty k^3\vert \hat\delta(k)\vert^2 \frac{\dd k}{k} = \int_{0}^\infty \frac{k^3\vert \hat\delta(k)\vert^2 }{2\pi^2} {\dd \ln k}

(59)

We define the matter power spectrum $P(k)$ by:

\left\langle \hat\delta(\vec k)\hat\delta^*(\vec k)\right\rangle = P(k) (2\pi)^3 \delta_D(\vec k + \vec k')

(60)

where $\delta_D$ denotes here the Dirac delta function, homogeneous to the inverse of a volume. In the sense of distributions, its Fourier transform equals 1 which allows us to deduce a formulation in the form of an integral:

\begin{align} \hat\delta_D(\vec k) & = \int \delta_D(\vec{r}) e^{i \vec k \cdot \vec r} \dd^3\vec r = 1 \\ \delta_D(\vec r) & = {1 \over (2 \pi)^3} \int \dd^3\vec k \hat\delta_D(\vec k) e^{i \vec k \cdot \vec r}= {1 \over (2 \pi)^3} \int \dd^3 \vec k e^{i \vec k \cdot \vec r} \end{align}

(61)

The power spectrum is therefore homogeneous to a volume: $\left[P(k)\right]=\mathsf{L}^3$ . For a real field, the conjugate Fourier modes verify $\hat\delta^*(k) = \hat\delta (-\vec k)$ so we can resume the definition of the power spectrum and calculate the mean of the norm of the Fourier modes:

\begin{align} \left\langle \vert \hat\delta(\vec k)\vert^2\right\rangle & = \left\langle \hat\delta(\vec k)\hat\delta^*(-\vec k)\right\rangle \\ & = \int_{\vec k'} \dd^3 \vec k' \left\langle \hat\delta(\vec k)\hat\delta^*(\vec k') \right\rangle \delta_D(\vec k + \vec k') \\ & = \int_{\vec k'} \dd^3 \vec k' (2\pi)^3 P(k) \delta_D(\vec k + \vec k') \delta_D(\vec k + \vec k')\\ & = (2\pi)^3 P(k) \int_{\vec k'} \dd^3 \vec k' \frac{1}{(2\pi)^3} \int_{\vec r} \dd^3\vec r e^{i (\vec k + \vec k') \cdot \vec r} \delta_D(\vec k + \vec k')\\ & = P(k) \int_{\vec r} e^{i \times 0} \dd^3\vec r = V P(k) \end{align}

(62)

with $V$ an integration volume, in practice finite because measurements or simulations are discrete and realized in a finite volume. For an isotropic field, we can perform the ensemble average by averaging over modes $\vec k$ of same norm $k$ . We obtain:

\boxed{P(k) = \left\langle \vert \hat\delta(\vec k)\vert^2\right\rangle_{\Vert \vec k\Vert = k} / V = \vert \hat\delta(k)\vert^2 / V }

(63)

The matter power spectrum at a scale $k$ is therefore the average of Fourier modes of same norm. We also define the modal logarithmic variance as:

\boxed{\Delta^2(k) = \frac{k^3}{2\pi^2}P(k)}

(64)

This is the contribution to the variance of modes of norm $k$ per logarithmic interval:

\sigma^2 = \int_{0}^\infty \Delta^2(k) {\dd \ln k}

(65)

In the context of linear structure growth, Fourier modes are independent. The matter power spectrum at a redshift $z$ can then be written:

P(k, z) = T^2(k, z) P_{\rm ini}(k)

(66)

where we define the transfer function $T(k,z)$ and the initial matter power spectrum $P_{\rm ini}(k)$ . What could it be? If we assume that at the beginning of the Universe density perturbations are random and there is no particular distance scale, a power law is appropriate:

P_{\rm ini}(k) \sim k^{n_s}

(67)

For initial conditions, we will prefer to focus on scalar perturbations of the metric $\delta \phi$ . But by Poisson’s equation $\hat\delta \sim k^2 \hat\delta \phi_k$ . We deduce the normalized power spectrum of scalar fluctuations:

\Delta_s^2(k) = A_s \left(\frac{k}{k_s}\right)^{n_s-1}

(68)

In the early 1970s, cosmologists Harrison and Zel’dovich proposed that all scales should contribute equally to the variance at the beginning of the Universe, no spatial scale counts more than others (Harrison (1970), Zeldovich (1972)). If this hypothesis is true, then we must have $n_s \approx 1$ .

The theoretical prediction of the spectrum $P(k)$ by looking at the evolution of different density perturbation scales along the history of the Universe is a convex curve, presenting a form in $k^{n_s}$ at large scales (too large to have evolved) and in $k^{-3}$ at small scales (those that collapse). The switch between the two regimes is the mode $k_{eq}$ corresponding to the horizon scale at equivalence (see Figure 6).

Matter power spectrum measured by different techniques and compared with theoretical linear (solid line) and non-linear (dotted line) predictions, adapted from doi:10.1093/mnras/stz2310. — Figure 6:Matter power spectrum measured by different techniques and compared with theoretical linear (solid line) and non-linear (dotted line) predictions, adapted from Chabanier *et al.* (2019).

Correlation Function and $\sigma_8^0$ ¶

The power spectrum is related to the correlation function in physical space by the following formulas:

\begin{align} \xi(r) & = \left\langle \delta(\vec r')\delta(\vec r'+ \vec r)\right\rangle_{\Vert \vec r\Vert = r} = {1 \over 2 \pi^2 } \int_0^\infty k^2 \frac{\sin(kr)}{kr} P(k) \dd k \\ P(k) & = \int_O^\infty r^2 \frac{\sin(kr)}{kr} \xi(k) \dd k \end{align}

(69)

We note that autocorrelation gives back the variance of the field: $\xi(0) = \sigma^2$ , and that the latter corresponds roughly to the integral of the matter power spectrum hence to its normalization. The measurement of the power spectrum normalization is therefore a way to measure matter clumping: the higher the power spectrum, the higher the variance of the matter field, hence the more matter is condensed into overdensities. From an experimental point of view, the power spectrum is measured by galaxy counting in spheres of radius $R=8\,\Mpc$ :

\delta_R(\vec r) = \int \dd^3\vec r \delta(\vec r') W_{R}(\vec r - \vec r'), \quad W_R(\vec r) = \left\lbrace \begin{array}{ll} 3/(4\pi R^3) & \text{ if }\Vert \vec r \Vert < R \\ 0 & \text{ if }\Vert \vec r \Vert \geqslant R \\ \end{array}\right.

(70)

\sigma_R^2 = \left\langle \delta_R^2(r)\right\rangle = \frac{1}{(2\pi)^3} \int_0^\infty P(k) \hat{W}^2_R(k) 4\pi k^2 \dd k = \frac{9}{2\pi^2 R^6} \int_0^\infty P(k) \left[ \sin(kR) - (kR)\cos(kR) \right]^2 \frac{\dd k}{k}

(71)

Today, Tristram, M. et al. (2024) measures:

\sigma_8^0 = 0.8186 \pm 0.0221

(72)

Since this value is close to 1, this means that at scales smaller than $8\,\Mpc$ we are already in a non-linear regime of structure growth since $\sigma \sim \delta \sim 1$ . Intuitively, we also understand that since structures collapse gravitationally over time, the variance of the matter field increases hence the normalization of $P(k)$ also (Figure 7).

<Figure size 640x480 with 1 Axes> — Figure 7:Evolution of the linear matter power spectrum with redshift.

From Galaxy Counting to Measurement of

\sigma_8^0

From the counting of galaxies $N_g$ in spheres of $8\,\Mpc$ , we deduce the density contrast for this galaxy field:

\delta_g(\vec r) = \frac{N_g(\vec r) - \left\langle N_g\right\rangle}{\left\langle N_g\right\rangle}

(73)

However, counting galaxies is not directly probing the underlying, invisible dark matter field. We consider that the number of galaxies $N_g(\vec r)$ is a Poisson random realization of the underlying total matter field $\delta(\vec r)$ , with a biased mean $b \delta$ . The bias $b$ depends on the galaxy type. According to their age, their history, they are more or less concentrated where dark matter is present: the bias typically varies between 1 and 4. The matter field is deduced from the galaxy field by:

\delta = \frac{\delta_g}{b}

(74)

Acoustic Waves in the Primordial Universe¶

Before decoupling, acoustic waves travel through the primordial plasma. Let us define the comoving sound horizon as:

r_s^c = \int_0^{t_\mathrm{dec}} \frac{c_s(t) \dd t}{a(t)} = \int_{z_\mathrm{dec}}^\infty \frac{c_s(z)}{H(z)}\dd z

(75)

the maximum comoving distance traveled by an acoustic wave since the beginning of the Universe. At the moment of decoupling, it is worth $150\,\Mpc/a_\mathrm{dec}$ or $r_s = a_\mathrm{dec} r_s^c \approx 0.13\,\Mpc$ . This is therefore the distance that a wave from an overdensity present at the Big Bang could travel. This wave propagation translates into a positive correlation on the presence of matter at this fundamental spatial scale. This scale converts into an angular separation on the sky, imprinted on the CMB temperature anisotropy map and given by:

\boxed{ \theta_s = \frac{r_s}{D_A(z_{\mathrm{dec}})} }

(76)

A scale in distance space gives a sinusoidal power spectrum in reciprocal space. Measuring the position of maxima in the CMB temperature anisotropy power spectrum leads to a precise measurement of the sound horizon Tristram, M. et al., 2024:

100 \theta_s = 1.04123 \pm 0.00046\,\mathrm{rad}, \quad \mathrm{i.e.}\quad \theta_s \sim 0.6\,\degree

(77)

Measuring the amplitude of the spectrum as well as the slope at large scales (small $\ell$ of Figure 10) allows to trace back to the parameters of the initial power spectrum:

\log(10^{10} A_s) = 3.073 \pm 0.061, \quad n_s = 0.9678 \pm 0.0072

(78)

We note that $n_s$ is close to 1, confirming the initial intuition of the Harrison-Zel’dovich spectrum. The fact that $n_s$ is slightly less than 1 (at more than 5\sigma) is a signature of inflation models (see next chapter). Finally, the fact that the second peak is attenuated compared to the first and third peak is a signature of the presence of baryons, which allows to measure:

\Omega_b^0 h^2 = 0.02224 \pm 0.00025 \quad \mathrm{i.e.}\quad \Omega_b^0 \sim 4.8\%

(79)

Figure 8:Time evolution of 10 Gaussian peaks under the effect of acoustic propagation.

Figure 9:Time evolution of a Gaussian field with an initial power spectrum $\propto k^{-3}$ . Left, the evolution of the map. Center, the stacking of 1000 thumbnails centered on initial peaks. Right, the evolution of the power spectrum.

Power spectrum of CMB temperature fluctuations as a function of angular mode \ell \sim \pi / \theta . — Figure 10:Power spectrum of CMB temperature fluctuations as a function of angular mode $\ell \sim \pi / \theta$ Tristram, M. *et al.*, 2024.

:name: fig:desi_peaks :align: center :width: 60%

Illustration of BAO formation in galaxy distribution. Credit: DESI collaboration

Adiabatic or Isocurvature Initial Conditions?

Initial fluctuations can be either adiabatic or isocurvature. In the adiabatic case, we have equality of fluctuations for all components of the Universe: $\delta_b = \delta_\gamma = \delta_{\rm dm}$ , and we predict the presence of all harmonics (1,2,3,...) in the CMB power spectrum. In the isocurvature case, the total density fluctuation is zero $\delta_{\rm tot}=0$ so we have oppositions like $\delta_b = -\delta_\gamma$ . We then predict the presence of only odd peaks in the CMB spectrum (1,3,5,...). The presence of all peaks indicates that initial conditions were adiabatic (Hu & Sugiyama (1995) p. 27).

Can We See the CMB Evolve Over Time?

To be reworked: https://www.reddit.com/r/askscience/comments/100snke/is_any_movement_visible_in_the_fluctuations_of/

On what time scale should we then expect to see the CMB change? The smallest scales we can currently resolve are about $0.07\degree$ in the sky, which corresponds to about \num{50 000} light-years (15\,\kpc) at the distance of the CMB.

For a structure of \num{50 000} light-years, the light from the far end takes \num{50 000} years longer to reach us than light from the near end. Does this mean we should expect CMB temperature fluctuations at these scales to change in about \num{50 000} years? Well, not quite. The CMB is redshifted by a factor of 1100, which means its apparent evolution in time is dilated by the same factor. We therefore expect fluctuations at the (currently) smallest resolved scales to change on a time scale (1100 times longer) of about 55 million years.

Optical Depth¶

We define the optical depth $\tau$ by the ratio of the number of photons received on Earth without having undergone any Thomson scattering $N(t_0)$ to the number of photons $N(t)$ emitted at a time $t$ :

\frac{N(t_0)}{N(t)} = e^{-\tau}

(80)

with

\tau = \int_{t_{\mathrm{ls}}}^{t_0} \Gamma_\gamma(t) \mathrm{d}t

(81)

The time $t_{\mathrm{ls}}$ for which $\tau=1$ is called the time of last scattering. This is the time since which a CMB photon has not interacted with an electron. More precisely,

\tau(z) = \int_0^z \frac{\Gamma_\gamma(z)}{H(z)}\frac{\dd z}{1+z}

(82)

This is one of the six parameters of the standard $\Lambda$ CDM model. Indeed, after the emission of the cosmic microwave background, we enter the Dark Ages of the Universe, when the Universe is transparent but no stars are yet emitting light. But with the birth of the first stars and galaxies, perhaps 150 million years after the Big Bang, the neutral medium is ionized once again. Although very sparse, CMB photons again interact with electrons by Thomson scattering, which reduces the amplitude of small-scale anisotropies in the CMB power spectrum, and introduces new anisotropies in polarization anisotropies. This is the least well measured parameter of the $\Lambda$ CDM model at the moment (Figure 12) $\tau = 0.058 \pm 0.006$ Tristram, M. et al., 2024, but it informs about the appearance of the first luminous objects (Figure 13).

The optical depth to reionization \tau measured as a smearing of the CMB angular power spectrum (from https://lambda.gsfc.nasa.gov/education/graphic_history/taureionzation.html, image credit: NASA / LAMBDA Archive Team). — Figure 12:The optical depth to reionization $\tau$ measured as a smearing of the CMB angular power spectrum (from https://lambda.gsfc.nasa.gov/education/graphic_history/taureionzation.html, image credit: NASA / LAMBDA Archive Team).

Key Points

The linear growth of large-scale structures in the Universe follows the laws of acoustics, but in an expanding space with gravitational forces at work.
Pressure forces cause matter overdensities to oscillate, gravitational force causes their collapse, and the expansion of space slows their growth.
After recombination, the sound speed drops from $c/\sqrt{3}$ to a few hundred meters per second: pressure forces no longer compensate for gravitational force and overdensities with masses comparable to galaxies can collapse under their own weight.
The statistical distribution of the matter field, through the power spectrum $P(k)$ , provides information on the growth of structures and on the balance between expansion and gravitation. The normalization of $P(k)$ gives the variance of the matter field and its shape provides information on the underlying physics.
The study of the CMB power spectrum makes it possible to constrain the 6 parameters of the standard model of cosmology $\Lambda$ CDM (flat): $\Omega^0_mh^2, \Omega^0_b h^2, 100\theta_s, \log(10^{10} A_s), n_s, \tau$ thanks to the angular size of the peaks, their amplitude, the relative height of the 2nd peak compared to the 3rd, and the slopes at small and large scales Tristram, M. et al., 2024.

Footnotes¶

if we had used the isothermal hypothesis, we would instead use the equation of state $PV=\cst$ and the speed of sound would be different by a factor $\sqrt{\gamma}$ .
↩
with the orders of magnitude of Figure 1 and $\rho_0\sim 6\text{ protons/m}^3$ , for a planet $\delta \sim 10^{32}$ whereas for a supercluster $\delta \sim 6$
↩
$\triangle_\sigma = \frac{1}{\sigma^2} \frac{\partial}{\partial \sigma}\left( \sigma^2 \frac{\partial}{\partial \sigma}\right)$
(83)
.
↩

References¶

Einstein, A. (1917). Kosmologische Betrachtungen zur allgemeinen Relativitätstheorie. Sitzungsberichte Der Königlich Preussischen Akademie Der Wissenschaften, 142–152. https://adsabs.harvard.edu/pdf/1917SPAW.......142E
Weinberg, S. (1989). The cosmological constant problem. Reviews of Modern Physics, 61(1), 1–23. 10.1103/RevModPhys.61.1
(1902). Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 199(312–320), 1–53. 10.1098/rsta.1902.0012
Harrison, E. R. (1970). Fluctuations at the Threshold of Classical Cosmology. Physical Review D, 1(10), 2726–2730. 10.1103/physrevd.1.2726
Zeldovich, Y. B. (1972). A Hypothesis, Unifying the Structure and the Entropy of the Universe. Monthly Notices of the Royal Astronomical Society, 160(1), 1P-3P. 10.1093/mnras/160.1.1p
Chabanier, S., Millea, M., & Palanque-Delabrouille, N. (2019). Matter power spectrum: from Ly α forest to CMB scales. Monthly Notices of the Royal Astronomical Society, 489(2), 2247–2253. 10.1093/mnras/stz2310
Tristram, M., Banday, A. J., Douspis, M., Garrido, X., Górski, K. M., Henrot-Versillé, S., Hergt, L. T., Ilić, S., Keskitalo, R., Lagache, G., Lawrence, C. R., Partridge, B., & Scott, D. (2024). Cosmological parameters derived from the final Planck data release (PR4). A&A, 682, A37. 10.1051/0004-6361/202348015
Hu, W., & Sugiyama, N. (1995). Toward understanding CMB anisotropies and their implications. Physical Review D, 51(6), 2599–2630. 10.1103/physrevd.51.2599