Table of contents

Back to the main page

1. Cosmological perturbation theory

Cosmological perturbation theory

So far we have worked with the flat FRLW metric $$\mathrm{d} s^{2}=-dt^2 + a^2(dx^2 + dy^2 + dz^2)$$ This metric describes a smooth Universe. However our Universe is obiously not perfectly smooth - even on fairly large scales. To study the evolution of structures we need to go beyond the smooth Universe. The main tool we have for this is cosmological perturbation theory. We imagine that the matter fluids (baryons, photons etc.) initially was not perfectly smooth, but had some small perturbations in them (we'll get back to how these initial perturbations were set up later when we go through inflation), i.e. $$\rho(t,\vec{x}) = \overline{\rho}(t) + \delta \rho(t,\vec{x})$$ for the density $$P(t,\vec{x}) = \overline{P}(t) + \delta P(t,\vec{x})$$ for the pressure $$\vec{v}(t,\vec{x}) = \overline{\vec{v}}(t) + \delta \vec{v}(t,\vec{x})$$ for the velocity field and so on. Here an overbar denotes a background quantity (e.g. $\overline{\rho}(t) \propto 1/a^3$ which is the density of matter in a perfectly smooth Universe) and $\delta \rho(t,\vec{x})$ the perturbation. In perturbation theory we assume the perturbations are small so for example a term like $(\delta \vec{v}(t,\vec{x}))^2$ would be ignored. Same if we have a product of different perturbations $\delta \rho_A(t,\vec{x}) \delta \rho_B(t,\vec{x})$. This means the resulting equations are going to be linear in the perturbations quantities which will be easy to solve. A downside of this is that the results are only valid when the perturbations are small $\delta A \ll A$. Luickly for us this is true for the CMB and also for clustering of matter on very large scales (on small scales it eventually breaks down and we would have to use numerical simulations to make progress, but that is a different course). Having introduced perturbations what we then are interested in is finding out how these perturbations $\delta\rho,\delta P,\ldots$ evolve. For example if we have a big overdensity then this will attract more matter gravitationally and grow larger unless the pressure is big enough so this will mean we finally get to some interesting dynamics and interplay between different physical effects. To do this we could (for the most parts) just go ahead and perturb these fluid quantities and derive evolution equations for them, however as we saw earlier the more fundamental quantity is the distribution function for each species $f_i$ so we will study perturbations to this quantity, i.e. we take $$f_i(t,\vec{x},\vec{p}) = \overline{f_i}(t,\vec{p}) + \delta f_i(t,\vec{x},\vec{p})$$ From this the perturbations in the energy density and pressure follows by using the Boltzmann formalism, but it has the advantage of allowing us to take into account interactions between different species. Now that we have perturbations in these quantities it means the energy momentum tensor, the right hand side of the Einstein equation, also gets perturbed which again means the metric will get perturbed. This again will change how matter and photons move and need to be taken into account when we look at the Boltzmann equation. We can already now expect to see the outline of an interesting interplay between different matter species and geometry taking place at the same time.


Figure: Interplay between different species and geometry that we need to take into account. Figure taken from Daniel Baumann's lecture notes.

To do all this consistently what we will do is to introduce perturbations to the distribution functions and the metric and derive the resulting coupled system of Boltzmann and Einstein equations to get a full set of evolution equations that we will later solve numerically (and will eventually lead us the CMB power-spectrum). Before we start with this we need to go through some important things regarding perturbations in General Relativity.

Metric perturbations

With perturbations in fluid quantities like energy density and pressure implies the energy momentum tensor, $T_{\mu\nu} = \overline{T}_{\mu\nu} + \delta T_{\mu\nu}$, is pertubed and the Einstein equations then tell us that we will also generate perturbations in the metric (and of course the other way around also: if we have perturbations in the metric we will over time generate perturbations in the fluids).

To perturb the metric we can write $$g_{\mu\nu} = \overline{g}_{\mu\nu} + \delta g_{\mu\nu}$$ where $\overline{g}_{\mu\nu}$ is the smooth FRLW metric. We could now go ahead and add a perturbation to each of the ($10$!) components. This quickly gets messy. Smart people in the past (Bardeen and friends) have found a better way of doing it. First one can show that the most general form for the perturbed metric can be written (the discussion below closely follows Baumann) $$ \mathrm{d} s^{2} = -(1+2 A) \mathrm{d} t^{2} + 2 aB_{i} \mathrm{~d} x^{i} \mathrm{~d} t + a^2\left(\delta_{i j}+h_{i j}\right) \mathrm{d} x^{i} \mathrm{~d} x^{j} $$ where $A, B_{i}$ and $h_{i j}$ are functions of space and time. We shall adopt the useful convention that Latin indices on spatial vectors and tensors are raised and lowered with $\delta_{i j},$ e.g. $h_{i}^{i}=\delta^{i j} h_{i j}$. What we can do next is to perform a scalar-vector-tensor decomposition of the perturbations. For 3-vectors, this should be familiar. It simply means that we can split any 3-vector into the gradient of a scalar and a divergenceless vector (Helmholtz decomposition; just as we do in electromagnetism when we decompose the electric field as $\vec{E} = -\nabla \phi + \nabla \times \vec{A}$) $$ B_{i}=\underbrace{\partial_{i} B}_{\text {scalar }}+\underbrace{\hat{B}_{i}}_{\text {vector }} $$ with $\partial^{i} \hat{B}_{i}=0$ (so this has only $2$ degrees of freedom). Similarly, any rank-2 symmetric tensor can be written $$ h_{i j}=\underbrace{2 C \delta_{i j}+2 \partial_{\langle i} \partial_{j\rangle} E}_{\text {scalar }}+\underbrace{2 \partial_{(i} \hat{E}_{j)}}_{\text {vector }}+\underbrace{2 \hat{E}_{i j}}_{\text {tensor }} $$ where $$ \begin{aligned} \partial_{\langle i} \partial_{j\rangle} E & \equiv\left(\partial_{i} \partial_{j}-\frac{1}{3} \delta_{i j} \nabla^{2}\right) E \\ \partial_{(i} \hat{E}_{j)} & \equiv \frac{1}{2}\left(\partial_{i} \hat{E}_{j}+\partial_{j} \hat{E}_{i}\right) \end{aligned} $$ As before, the hatted quantities are divergenceless, i.e. $\partial^{i} \hat{E}_{i}=0$ (so this has only $2$ degrees of freedom) and $\partial^{i} \hat{E}_{i j}=0$. The tensor perturbation is also traceless, $\hat{E}_{i}^{i}=0$ (so this contains only $2$ degrees of freedom corresponding to the two polarization states of the graviton). The 10 degrees of freedom of the metric have thus been decomposed into $4$ scalar, $4$ vector and $2$ tensor degrees of freedom:

  • scalars: $A, B, C, E$
  • vectors: $\hat{B}_{i}, \hat{E}_{i}$
  • tensors: $\hat{E}_{i j}$

What makes this decomposition so powerful is the fact that the Einstein equations for scalars, vectors and tensors don't mix at linear order and can therefore be treated separately! In this course, we will mostly be interested in scalar fluctuations and the associated density perturbations. Vector perturbations aren't produced by inflation and even if they were, they would decay quickly with the expansion of the universe (exam problem 2020). Tensor perturbations - gravitational waves - are an important prediction of inflation and we will get back to them in the end of this course.

The Newtonian gauge

The scalar-vector-tensor decomposition of the perturbations gives us a great simplification: we only have $4$ free functions to work with instead of $10$. However the metric perturbations as we wrote them above are not uniquely defined, but depend on our choice of coordinates, the so-called gauge choice. When we wrote down the perturbed metric, we implicitly chose a specific timeslicing of the spacetime and defined specific spatial coordinates on these time slices. Making a different choice of coordinates, can change the values of the perturbation variables. It may even introduce fictitious perturbations. These are fake perturbations that can arise by an inconvenient choice of coordinates even if the background is perfectly homogeneous. We won't go through this in great detail in this course, but you can read Baumann for more info about this. We have two options to deal with this: either choose a gauge or find combinations of variables that do not change from gauge to gauge (gauge invariant quantities). This is a similar situation to what you might have encountered before in electromagnetism: there is no unique choice for the electromagnetic four potential $A^\mu$, so you are free to make a choice for this when doing calculations. In the end observables always turn out to be gauge invariant quantities so its of no importance what gauge we pick and we pick the one that makes our life easy. Its the same story in General Relativity. The important part for us is that this coordinate freedom implies that we are free to set a few of the perturbation to zero so the end result is that we only have two independent free functions. The full form for the perturbed metric that we will use is given by $$\mathrm{d} s^{2}= -dt^2(1+2\Psi) + a^2(1+2\Phi)(dx^2 + dy^2 + dz^2)$$ where $\Psi,\Phi$ are functions of both time and space (and when they are zero we recover the smooth FRLW metric; NB: Dodelson uses a different sign convention than in Baumann with $\Phi\to -\Phi$ so be aware of this if you use both of these texts). This choice is called the Newtonian gauge and the potentials $\Psi,\Phi$ are closely related to the Newtonian gravitational potential (see the discussion on the Newtonian limit of GR). This choice is convenient as the physics is (in my opinion) much more transparent than in some other common gauge choices so we can use our Newtonian intuition when interpreting the equations we end up with.

Perturbations in real and fourier space

The perturbation equations we end up with are going to be coupled linear partial differential equations. Solving these in real space are not trivial, but luckily there is a solution to this by using Fourier transforms (now is the time to go and read the introduction to the Fourier transform). To motivate this consider the wave equation $$\frac{\partial^2 u(x,t)}{\partial t^2} - c^2 \frac{\partial^2 u(x,t)}{\partial x^2} = 0$$ whose solution depends on position $x$ and time $t$. If we take the Fourier transform of this equation, i.e. we multiply it by $e^{-ikx}$ and integrate over all $x$ using $\hat{u}(k,t) = \int_{\mathbb{R}} dx u(x,t)e^{-ikx}$, then we get $$\frac{\partial^2 \hat{u}(k,t)}{\partial t^2} + c^2k^2 \hat{u}(k,t) = 0$$ For a fixed wavenumber $k$ this is now an ordinary differential equation (and we can in fact easily write down the analytical solution if we want). This shows that the Fourier transform is able to take a more complicated PDE and turn it into a simple ODE (that we already know how to solve). The downside is that we have to solve one ODE for each $k$ to be able to find the full solution, but there are no free lunches so this we have to live with. For any linear PDE with constant coefficients (here in $1+1$ dimension) we can transform it to Fourier space by simply making the replacement $$u(x,t)\to \hat{u}(k,t)$$ and for terms that contain spatial derivatives $$\frac{d^n u(x,t)}{dx^n} \to (ik)^n \hat{u}(k,t)$$ Note that this only works if the equation is linear as the Fourier transform of a product $u(x,t)v(x,t)$ does not equal $\hat{u}(k,t)\hat{v}(k,t)$, but rather is a convolution $\hat{u}(k,t) * \hat{v}(k,t) \equiv \int_{-\infty}^\infty \hat{u}(k-k',t)\hat{v}(k',t)dk'$. For functions of more than one spatial variable we have the equivalent rule $\nabla u \to i\vec{k} \hat{u}$. For example the linear continuity equation $$\frac{d \delta(x,t)}{dt} + \nabla \cdot \vec{v} = 0 \to \frac{d \hat{\delta}(k,t)}{dt} + i\vec{k} \cdot \vec{\hat{v}}(k,t) = 0$$ This simple mapping allows us to easy go back and forth between a real and fourier description. People often do the notational abuse of using the same symbol (i.e. no hat) for both the real quantity and the fourier component and let the context determine what you are talking about (i.e. $\frac{\partial^2 u}{\partial t^2} + c^2k^2 u = 0$ for the equation above, the "$k^2$" shows that this is a fourier equation). In the lectures I'll try to be careful with the notation, but in these notes I will for sure abuse the notation like this as its too much work to add a hat everywhere so get used to this.

Our plan

We are now ready to start. The goal is to derive a full set of equations that determines how structures form in our Universe. We will perturb the distribution functions for each of the quantities baryons (electrons and protons) $f_e,f_p$, photons $f_\gamma$, massless neutrinos $f_\nu$ and cold dark matter $f_{\rm CDM}$. We will for each of these evaluate the Boltzmann equation $\frac{df_i}{dt} = C_i$ where the right hand side deals with relevant interactions between the different species. The main interactions we have to deal (which requires some results from quantum field theory) is Compton scattering $e^- + \gamma \rightleftharpoons e^- + \gamma$ and Coloumb sattering $e^- + p^+ \rightleftharpoons e^- + p^+$. We will also need to evalate the geodesic equation that tells us how things move in a perturbed Universe so the equations we find will depend on the metric perturbations. Having found equations for the different species this will give us a perturbed energy momentum tensor $\delta T_{\mu\nu}$. We will then take the perturbed metric and derive the same quantities as we did for the background metric: Christoffel symbols, Ricci tensor and finally the Einstein tensor. This will give us perturbed Einstein equations $\delta G_{\mu\nu} = 8\pi G \delta T_{\mu\nu}$ that tells us what the metric perturbations are given the perturbation in the fluids. In the end we will end up with a closed system of coupled differential equations for how the metric perturbation, the baryon density, baryon velocity, CDM density, CDM velocity and so on evolves that we will solve in Fourier space. This system of equations is also what you will implement and solve numerically and once this is done we will just be a couple of integrals away from our final goal: key observables such as the CMB angular power-spectrum and the matter power-spectrum.

Perturbations of the Einstein-Boltzmann equations

The Boltzmann equation tells us how the distribution function of some particle species evolves in time $$\frac{df}{d\lambda} = C(f)$$ where $\lambda$ is an affine parameter along a trajectory and $C$ is the collision term which contains information about interactions. In the background this was not so bad, it only depends on momentum $p$ and time $t$, but introducing perturbations it gets more messy. It will now depend on time $t$, position $\vec{x}$ and momentum $\vec{p}$ (magnitude and direction) so its a seven dimensional function.

Introducing phase space coordinates ($x^\mu, P^\mu$), where $x^\mu$ is the position spacetime-coordinates and $P^\mu = \frac{dx^\mu}{d\lambda}$ is the momentum-coordinates we can expand $$\frac{df}{d\lambda} = \frac{\partial f}{\partial t} \frac{dt}{d\lambda} + \frac{\partial f}{\partial x^i} \frac{dx^i}{d\lambda} + \frac{\partial f}{\partial E}\frac{dE}{d\lambda} + \frac{\partial f}{\partial \hat{p}^i}\frac{d\hat{p}^i}{d\lambda}$$ Note that we have choosen to use $(p,\hat{p})$ (magnitude and direction) as our momentum variables instead of the vector itself $\vec{p}$. This is just for convenience later on as we will see the last term in the equation above will be second order in the perturbations and can therefore be ignored.

The terms multiplying the derivatives of the distribution functions descibes how particles/photons move in spacetime. To understand these we need to look at the geodesic equation.

Geodesic equation in a perturbed Universe

To evaluate the terms in the Boltzmann equation we need to know how photons and particles move in a perturbed Universe. Let $x^\mu$ be the position spacetime-coordinates and $P^\mu = \frac{dx^\mu}{d\lambda}$ the 4-momentum. The 4-momentum satisfy $g_{\mu\nu}P^\mu P^\nu = -m^2$. We define the usual 3-momentum via $p^2 \equiv g_{ij}P^i P^j$ and the direction $\hat{p}^i \propto P^i$ normalized such that $\hat{p}^i\hat{p}^j\delta_{ij} = 1$. Its also convenient to introduce the energy $E$ via the energy-momentum relation $E^2 = p^2 + m^2$. This gives $$g_{00}(P^0)^2 + p^2 = -m^2 \to P^0 = E(-g_{00})^{-1/2} = \frac{E}{\sqrt{1+2\Psi}} \simeq E(1-\Psi)$$ The physical effect here is gravitational redshift: the wavelength gets blueshifted (redshifted) when we move in (out of) a gravitational well. We also find $$P^i = p\hat{p}^i(g_{ii})^{-1/2} = \frac{p\hat{p}^i}{a\sqrt{1+2\Phi}} \simeq \frac{p\hat{p}^i}{a}(1-\Phi)$$

To be able to evaluate the terms $\frac{dE}{d\lambda}$ and $\frac{d\hat{p}^i}{d\lambda}$ in the Boltzmann equation we need to look at the geodesic equation. This is just $$\frac{dP^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} P^\alpha P^\beta = 0$$ and for $\mu=0$ we get $$\frac{d[E(1-\Psi)]}{d\lambda} + \Gamma^0_{\alpha\beta} P^\alpha P^\beta = 0$$ The first term is $$\frac{d[E(1-\Psi)]}{d\lambda} = (1-\Psi)\frac{dE}{d\lambda} - E\frac{d\Psi}{d\lambda}$$ and the latter term can be expanded using the chain rule as $$\frac{d\Psi}{d\lambda} = \frac{\partial \Psi}{\partial x^\mu} \frac{dx^\mu}{d\lambda} = \frac{\partial \Psi}{\partial x^\mu} P^\mu$$ Finally we need the Christoffel symbols $$\Gamma^0_{\alpha \beta} P^\alpha P^\beta = \Gamma^0_{00} (P^0)^2 + 2\Gamma^0_{0i} P^0P^i + \Gamma^0_{ij} P^iP^j\\ \simeq \frac{\partial \Psi}{\partial t}E^2 + 2\frac{\partial \Psi}{\partial x^i}E \frac{p\hat{p}^i}{a} + p^2(H-2H\Psi + \frac{\partial \Phi}{\partial t})$$ Collecting terms gives us the desired expression for $\frac{dE}{d\lambda}$: $$\frac{dE}{d\lambda} \simeq -p^2 \left[H - H\Psi + \frac{\partial \Phi}{\partial t} + \frac{E}{p}\frac{\partial \Psi}{\partial x^i}\frac{\hat{p}^i}{a} \right]$$ This equation tells us how the energy of a, say, photon changes along a trajectory. In the absence of perturbations this just tells us that $\frac{1}{P^0}\frac{dp}{d\lambda} = \frac{dp}{dt} = -pH$ or $p\propto 1/a$ - the usual redshifting of photons due to expansion. The different terms on the right hand side represent redshifting due to the expansion of space (and some points in space are effectively a bit ahead or behind the mean expansion, $a_{\rm local} \sim a(1+\Phi)$, so there are corrections to how this looks in a smooth Universe) and gravitatational redshift: clocks ticks slower so the frequency of a photon would be larger deep in a gravitational well.

We could go on and look at the $\mu = i$ component of the geodesic equation to compute $\frac{d\hat{p}^i}{d\lambda}$, but as we mentioned previously we won't actually need this.

Boltzmann equation for Photons

For photons $E = p$ so the geodesic equation becomes $$\frac{dp}{dt} = \frac{1}{P^0}\frac{dp}{d\lambda} \simeq -p \left[H + \frac{\partial \Phi}{\partial t} + \frac{\partial \Psi}{\partial x^i}\frac{\hat{p}^i}{a} \right]$$ Before we go ahead lets talk a bit more about what this equation tells us. Introducing the total time-derivative $\frac{dX}{dt} = \frac{\partial X}{\partial t} + \frac{\partial X}{\partial x^i}\frac{dx^i}{d\lambda} \frac{d\lambda}{dt}$ with the last two terms being $P^i / P^0$ we can write it as $$\frac{d\log p}{dt} + \frac{d\log a}{dt} + \frac{d\Psi}{dt} \simeq \left[\frac{\partial \Psi}{\partial t} - \frac{\partial \Phi}{\partial t}\right]$$ We can formally integrate it along the trajectory of a photon from a source (recombination) to an observer (today) to get $$\log\left(\frac{(pa)_{\rm today}}{(pa)_{\rm recombination}}\right) = (\Psi_{\rm recombination} - \Psi_{\rm today}) + \int_{t_{\rm recombination}}^{t_{\rm today}} \left[\frac{\partial \Psi}{\partial t} - \frac{\partial \Phi}{\partial t}\right]dt$$ Writing the momentum of a photon as $p = \overline{p} + \delta p$ where $\overline{p}$ just redshifts as $1/a$ due to expansion the equation above to first order gives us $$\left(\frac{\delta p}{p}\right)_{\rm today} = \left(\frac{\delta p}{p}\right)_{\rm recombination} + (\Psi_{\rm recombination} - \Psi_{\rm today}) + \int_{t_{\rm recombination}}^{t_{\rm today}} \left[\frac{\partial \Psi}{\partial t} - \frac{\partial \Phi}{\partial t}\right]dt$$ This tells us that the frequency of photons (and thereby their temperature) change due to the presence of structure (and note that this fractional change does not depend on the frequency of photons - all photons of different frequencies are affected equally). This is called the Sachs-Wolfe effect and describes two important sources anisotropies in the cosmic microwave background: gravitational redshift of photons ("hot" photons are in dense environments with largeish $\Psi$ that they have to climb out of on their way to us) and the so-called integrated Sachs-Wolfe effect: change of energy due to time-varying gravitational wells (e.g. dark energy makes gravitational wells decay at late time leaving an imprint on the CMB). We will get back to this when we in the end of the course will try to understand the CMB power-spectrum and disucss the various physical effects that leaves imprints on it. Using the geodesic equation the Boltzmann equation to linear order now becomes $$\frac{df}{d\lambda} \simeq P^0\left[\frac{\partial f}{\partial t} + \frac{\partial f}{\partial x^i} \frac{\hat{p}^i}{a} - p(H + \frac{\partial \Phi}{\partial t} + \frac{\partial \Psi}{\partial x^i}\frac{\hat{p}^i}{a})\frac{\partial f}{\partial p}\right]$$ Next we introduce perturbations to the distribution function: $$T = \overline{T}(1 + \Theta(t,\vec{x},p,\hat{p}))$$ or in other words $\Theta = \frac{\delta T}{T}$ - the perturbation in the local photon temperature. For photons the relevant interaction is Compton scattering (or Thompson scattering as the low-energy limit of Compton scattering is often called) and to first order this does not change the momentum of the photons only the direction (see Dodelson) so we have a simplication that $\Theta(t,\vec{x},p,\hat{p}) = \Theta(t,\vec{x},\hat{p})$. The distribution function can be expanded as $$f = \frac{1}{e^{\frac{p}{\overline{T}(1+ \Theta)}}-1} \simeq \overline{f} - \Theta p\frac{\partial \overline{f}}{\partial p}$$ where $\overline{f} = \frac{1}{e^{\frac{p}{\overline{T}}}-1}$. At the background level the Boltzmann equation reads $$\frac{\partial \overline{f}}{\partial t} - pH\frac{\partial \overline{f}}{\partial p} = 0$$ and since $\overline{f}$ is a function of $\frac{p}{\overline{T}(t)}$ we have $$\frac{\partial \overline{f}}{\partial t} = -p\frac{1}{\overline{T}}\frac{d\overline{T}}{dt}\frac{\partial \overline{f}}{\partial p}$$ the 0th order solution is $\overline{T} \propto 1/a$ which just tells us what we already know: photons are redshifted as $1/a$ as the Universe expands. From this we find $$\frac{\partial f}{\partial t} \simeq \frac{\partial \overline{f}}{\partial t} - \frac{\partial \Theta}{\partial t} p \frac{\partial \overline{f}}{\partial p} - \Theta p \frac{\partial^2 \overline{f}}{\partial p \partial t}$$ $$\frac{\partial f}{\partial p} \simeq \frac{\partial \overline{f}}{\partial p} - \Theta\frac{\partial \overline{f}}{\partial p} - \Theta p\frac{\partial^2 \overline{f}}{\partial p^2}$$ $$\frac{\partial f}{\partial x^i} \simeq - \frac{\partial \Theta}{\partial x^i} p\frac{\partial \overline{f}}{\partial p}$$ Note that we have ignored the $p$ dependence on $\Theta$ as discussed above. Taking the $p$ derivative of the 0th order Boltzmann equation we also see that $$\frac{\partial^2 \overline{f}}{\partial p\partial t} = H\frac{\partial \overline{f}}{\partial p} + Hp\frac{\partial^2 \overline{f}}{\partial p^2}$$ and this equality in turn will ensure that the second derivative terms of the distribution function that we see in the formulas above will cancel out in the Boltzmann equation. The Boltzmann equation with the background subtracted (using the equation above to cancel terms) then gives $$\frac{df}{d\lambda} = -P^0p\frac{\partial \overline{f}}{\partial p}\left[\frac{\partial \Theta}{\partial t} + \frac{\partial \Theta}{\partial x^i} \frac{\hat{p}^i}{a} + (\frac{\partial \Phi}{\partial t} + \frac{\partial \Psi}{\partial x^i}\frac{\hat{p}^i}{a})\right]$$ To continue we need to look the at the collision term which is usually where most of the complications enter. We are not going to cover this in great detail here, take a look a page 95 in Dodelson for a full derivation. For photons the main interaction we need to consider is Compton scattering: scattering of light by free electrons (for why other interactions can be ignored see the home exam 2020): $$e^- + \gamma \rightleftharpoons e^- + \gamma$$ The collision term for the process to first order in perturbation theory (ignoring polarization and the angular dependence of the Thompson cross-section which will add an additional term - this is the expression Dodelson derives) is given by $$C = -p^2\frac{\partial \overline{f}}{\partial p}n_e \sigma_T[\Theta_0 - \Theta + \hat{p}\cdot \vec{v}_b]$$ where $\vec{v}_b$ is the baryon velocity (how fast electrons - that are tightly coupled to protons - are moving), $n_e$ the free electron density (the more free electrons the bigger the effect), $\sigma_T$ the Thompson cross-section (how "efficient" this interaction is) and $\Theta_0 \equiv \frac{1}{4\pi}\int d\Omega_{\hat{p}} \Theta = \frac{1}{2}\int_{-1}^1 \Theta d\mu$ is the monopole moment of the photon distribution (the local temperature averaged over all directions - and the first example of a "multipole" making its appearance). This gives us the final result $$\frac{\partial \Theta}{\partial t} + \frac{\partial \Theta}{\partial x^i} \frac{\hat{p}^i}{a} + (\frac{\partial \Phi}{\partial t} + \frac{\partial \Psi}{\partial x^i}\frac{\hat{p}^i}{a}) = n_e \sigma_T[\Theta_0 - \Theta + \hat{p}\cdot \vec{v}_b]$$ This equation tells us that when Compton scattering is efficient then the photon distribution will get reduces to $\Theta \approx \Theta_0$ so at any point in space the photons from different directions will have the same temperature. Lets talk a bit about other interactions. Why can we ignore $p^+ + \gamma \rightleftharpoons p^+ + \gamma$? The reason is that the Thompson cross-section is proportional to $1/m^2$ and with protons being much heavier than electrons this process is irrelevant in comparison. And what about $\gamma + \gamma \rightleftharpoons \gamma + \gamma$? This process does not transfer any energy and momentum away from photons (like Compton scattering can do) so its also irrelevant. If one wants to go beyond first order perturbation theory, which is significantly more complicated, then some more interactions needs to be included but that is way beyond the scope of this course. Our final equation might look innocent enough, but it is a function of a $\hat{p}$ (2 numbers), $\vec{x}$ (3 numbers) and $t$ so 6 free variables. The first thing we do to simplify it is to go to fourier space and its then convenient to introduce the variable $\mu \equiv \frac{\hat{p}\cdot \vec{k}}{k}$ - the cosine of the angle between the Fourier wavevector $\vec{k}$ and the photon direction $\hat{p}$. The baryons velocity is curl-free so in Fourier space $\vec{v_b} = iv_b \frac{\vec{k}}{k}$ where $v_b$ is the magnitude of $\vec{v_b}$ - so we only need to track the magnitude which is another simplification. With all this the equation becomes $$\frac{\partial \Theta}{\partial t} + \frac{ik\mu}{a} \Theta + (\frac{\partial \Phi}{\partial t} + \frac{i k\mu}{a}\Psi) = n_e \sigma_T[\Theta_0 - \Theta + i\mu v_b - \frac{3\mu^2-1}{4}\Pi]$$ where the final term we did not include earlier $\Pi = \Theta_2 = -\frac{1}{2}\int_{-1}^1\Theta \frac{3\mu^2-1}{2}d\mu$ comes from the angular dependence of the Compton scattering and also has additional contributions from how polarized the CMB is that is ignored in this formula, but we'll come back to this later. We now have $\Theta = \Theta(t,k,\mu)$ a function of just $3$ variables, a great simplification from the $7$ we started with. Some of the simplification comes from the convenient assuption we were allowed to make that $\Theta$ does not depend on the momenta $p$ and the rest comes from isotropy of the background and the fact that we only consider scalar perturbations which ensures that directional dependence (on the Fourier wavevector and the photon direction) always appears in the form $\mu = \hat{p} \cdot\vec{k} / k$ in our equations. Later on we will only expand this equation in Legendre multipoles leaving us with a coupled system of equations for the multipoles $\Theta_\ell(k,t)$ that is the final form we are going to implement and solve numerically. I should also mention that Compton scattering can create a linear polarisation of photons, an effect we now can measure in CMB experiments, but we'll get back to this later in this course (this will just modify the $\Pi$ term).

Boltzmann equation for Cold Dark Matter

Dark matter is nice in that it has no (or very weak) interactions with the normal standard model particles. Thus its collision term is zero and the Boltzmann equation simply reads $$\frac{\partial f}{\partial t} + \frac{\partial f}{\partial x^i} \frac{P^i}{P^0} + \frac{\partial f}{\partial E}\frac{1}{P_0}\frac{dE}{d\lambda} = 0$$ where we have ignored the term $\frac{\partial f}{\partial \hat{p}^i}\frac{d\hat{p}^i}{d\lambda}$ as this is a second order term for exactly the same reasons as for photons. We also know $$\frac{1}{P^0}\frac{dE}{d\lambda} \simeq -\frac{p^2}{E} \left[H + \frac{\partial \Phi}{\partial t} + \frac{E}{p}\frac{\partial \Psi}{\partial x^i}\frac{\hat{p}^i}{a} \right]$$ giving us $$\frac{\partial f}{\partial t} + \frac{\partial f}{\partial x^i} \frac{\hat{p}^i}{a} \frac{p}{E} - \frac{\partial f}{\partial E}\left[H\frac{p^2}{E} + \frac{\partial \Phi}{\partial t} \frac{p^2}{E} + \frac{\partial \Psi}{\partial x^i}\frac{p\hat{p}^i}{a} \right] = 0$$ This is exactly the same as for photons with the exception that we haven't taken $E=p$, since dark matter has mass so $E^2 = p^2 + m^2$ instead. Instead of perturbing the distribution function as we did for photons it's easier to work with the integrated Boltzmann equation. Why? because dark matter is non-relativistic at the times we are interested in and will only have a sizable monopole and dipole term, i.e. the evolution is fully characterized by the density and velocity, so the conservation equation that we get from the zeroth and first moment of the Boltzmann equation is enough to determine their evolution. So we multiply the equation above by $\frac{g d^3 p}{(2\pi)^3}$ and integrate over all momenta to get $$\frac{g}{(2\pi)^3}\int\frac{\partial f}{\partial t} d^3p+ \frac{g}{(2\pi)^3}\int \frac{\partial f}{\partial x^i} \frac{\hat{p}^i}{a} \frac{p}{E} d^3p - \frac{g}{(2\pi)^3}\int \frac{\partial f}{\partial E}\left[H\frac{p^2}{E} + \frac{\partial \Phi}{\partial t} \frac{p^2}{E} + \frac{\partial \Psi}{\partial x^i}\frac{p\hat{p}^i}{a} \right]d^3p = 0$$ Recall what we did when we talked about integrating the Boltzmann equation? If not go back and read that again because we will do exactly the same here. The first term is trivial, we just pull the $t$-derivative out of the integral and get $\frac{\partial n}{\partial t}$. For the second term we pull the $x$-derivative out of the integral and get $$\frac{\partial }{\partial x^i} \left(\frac{g}{(2\pi)^3}\int f \frac{\hat{p}^i}{a} \frac{p}{E} d^3 p\right) = \frac{1}{a}\frac{\partial }{\partial x^i}(n v^i)$$ where the fluid velocity is the average of the momentum over mass which in the non-relativistic case energy is just the mass to first order so $$v^i = \frac{1}{n}\left\lt \frac{p^i}{m}\right\gt \simeq \frac{1}{n}\left\lt \frac{p\hat{p}^i}{E}\right\gt$$ Finally we need to deal with the third term. The term with $\Psi$ is easy: since $\Psi$ is first order and any directional dependence is also first order this term is second order and can be neglected. The last terms requires some integration by parts to move the derivative from the distribution function over to what multiplies it. For this remember that $\Phi$ and $H$ does not depend on momenta, $d^3p = p^2 dp d\Omega$, $E^2 = p^2 + m^2\implies EdE = pdp$ and $\frac{\partial}{\partial E} = \frac{dp}{dE}\frac{d}{dp} = \frac{E}{p}\frac{d}{dp}$ thus $$\int \frac{\partial f}{\partial E} \frac{p^2}{E}d^3p = \int \frac{\partial f}{\partial p} p^3 dpd\Omega = -\int f \frac{\partial p^3}{\partial p} dpd\Omega = -3\int f d^3p$$ Thus the third term just becomes $$-[H + \frac{\partial \Phi}{\partial t}]\cdot(-3n)$$ and the full result is $$\frac{\partial n}{\partial t} + \frac{1}{a}\frac{\partial }{\partial x^i}(n v^i) + 3n[H + \frac{\partial \Phi}{\partial t}] = 0$$ This you should recognice as the usual continuity equation describing conservation of particle number / mass just with a few extra terms that comes from the fact that universe is expanding (without perturbations $n\propto 1/a^3$ which gives rise to the $3nH$ term) and a related term that comes from the fact that the local value of the scale-factor $a_{\rm local}\sim a(1+\Phi)$ so the volume in a perturbed universe is $V_{\rm local} = a^3(1+\Phi)^3 \sim a^3(1 + 3\Phi)$ which gives rise to the last term, i.e. the terms not depending on velocity can be combined into $\frac{1}{a_{\rm local}^3}\frac{\partial (na_{\rm local}^3)}{\partial t}$ so the equation is in fact on exactly the same form as we saw in the background.

We also need an equation that tells us how the velocity evolves. Again this follows the same recipe as we followed in the non-relativistic example previously discussed in the course by taking the first moment and derive the corresponding Euler equation which describes momentum conservation. Thus if we multiply the Boltzmann equation by $\frac{g}{(2\pi)^3}\frac{p\hat{p}^i}{E}$ and integrate over all momenta (we leave this calculation as an exericise) we get $$\frac{\partial v^i}{\partial t} + H v^i = -\frac{1}{a}\frac{\partial \Psi}{\partial x^i}$$ We recognize the right hand side as the gravitational force and the two terms on the left hand side is just $\frac{1}{a}\frac{\partial (av^i)}{\partial t}$ so expressed in terms of the proper velocity $av^i$ its exactly like the usual Euler equation. Note that we don't have the $(v\cdot\nabla)v$ term since this would be second order and there is not pressure term as dark matter is pressureless.

Boltzmann equation for Baryons

We are going to do the same for baryons as we did for dark matter above, but with the exception that we need to take into account interactions. For baryons (by which we in cosmology mean electons and protons) the main interactions are Coulomb scattering $$e^- + p^+ \leftrightharpoons e^- + p^+$$ and Compton scattering of electrons and protons $$e^- + \gamma \leftrightharpoons e^- + \gamma$$ $$p^+ + \gamma \leftrightharpoons p^+ + \gamma$$ We already mentioned why the latter interaction is irrelevant: the cross-section is proportional to $1/m^2$ and since protons are much more massive than electrons we can ignore this. Lets start with Coulomb scattering. The main "job" of this interaction is to keep electrons and protons tightly coupled to eachother. If electrons clustered away from protons then the resulting electric force (which is much stronger than gravity) would quickly wash out these differences. This ensures that the perturbations in the density and velocity of electrons and protons stay the same and this common value is what we call the baryon density and velocity: $$\delta_p = \delta_e \equiv \delta_b$$ $$\vec{v}_p = \vec{v}_e \equiv \vec{v}_b$$ We could now write down the collision term for these interactions and do the math, but we don't have to - we can figure this our by some simple arguments! Recall the thing we have stressed so many times about what the first two moments of the Boltzmann equation tell us: conservation of particle number and momentum. For both interactions the number of baryons don't change and consequently none of them contribute to the continuity equation. Thus the result for dark matter applies and we can immediately write it down for baryons $$\frac{\partial n_b}{\partial t} + \frac{1}{a}\frac{\partial }{\partial x^i}(n_b v_b^i) + 3n_b[H + \frac{\partial \Phi}{\partial t}] = 0$$ That was easy. What about the Euler equation? Without the interactions we know the result from dark matter $$\frac{\partial v_b^i}{\partial t} + H v_b^i = -\frac{1}{a}\frac{\partial \Psi}{\partial x^i} + {\rm momentum~transfer}$$ where the momentum transfer term is not present in dark matter and depends on the interactions. For Coulomb scattering there is no momentum transfer out of baryons: we just move the momentum from one baryon to another so this interaction does not contribute to the Euler equation. However for Compton scattering we have already seen that there is momentum transfer from/to the photons so this interaction will contribute. By how much then, surely we need to evaluate the collision term to find this? Well no, we know how exactly how much momentum is transfered to the photons and momentum is conserved so what is gained (lost) by photons must be lost (gained) by baryons! The Euler equation for photons is (this will be derived a bit later in Multipole expansion by taking moments of the photon equation we derived so just trust me for now): $$v_\gamma^\prime = -\frac{ck}{4\mathcal{H}} \delta_\gamma + \frac{2ck}{\mathcal{H}}\Theta_2 - \frac{ck}{\mathcal{H}}\Psi + \tau^\prime\left[v_\gamma - v_b\right]$$ where the last term is the momentum transfer due to Compton scattering. The momentum density (this is just $\delta T^i_0$) is $(\overline{\rho}+\overline{P})\vec{v}$ so $(\overline{\rho}_b+\overline{P}_b)\vec{v}_b + (\overline{\rho}_\gamma+\overline{P}_\gamma)\vec{v}_\gamma$ cannot depend on the interaction so the Euler equation for $\vec{v}_b$ must contain the opposite term as in the photon equation for $\vec{v}_\gamma$ weighted by the relative densities $\frac{\overline{\rho}_\gamma + \overline{P}_\gamma}{\overline{\rho}_b + \overline{P}_b} = \frac{4\overline{\rho}_\gamma}{3\overline{\rho}_b}$ (since $w = 1/3$ for photons and $w = 0$ for baryons). This argument shows that the momentum transfer in the baryon equation should be $${\rm momentum~transfer} = - \frac{4\overline{\rho}_\gamma}{3\overline{\rho}_b}\tau^\prime (v_\gamma - v_b)$$ and the equation therefore must be $$\frac{\partial v_b^i}{\partial t} + H v_b^i = -\frac{1}{a}\frac{\partial \Psi}{\partial x^i} - \frac{4\overline{\rho}_\gamma}{3\overline{\rho}_b}\tau^\prime (v_\gamma - v_b)$$ If you think this derivation is sketchy and prefer a hard mathematical derivation then you should do the full derivation as outlined in Dodelson Chapter 4.6 (its great practice).

Boltzmann equation for Massless neutrinos

This one is really easy. If we perturb the neutrino distribution function just as we did for photons setting $$\mathcal{N} = \frac{\delta T_\nu}{\overline{T}_\nu}$$ then exacty the same calculations that we did for photons applies for neutrinos except for the collision term which is zero after neutrino decoupling. There is one important difference and that is that neutrinos have a different temperature than photons $T_\nu \not = T_\gamma$, however that does not enter in the equation as we consider $\frac{\delta T_\nu}{\overline{T}_\nu}$. If you look back through the derivation you can see that the only fact we used is that $T \propto 1/a$ which is true also for neutrinos. Thus taking the photon equation $$\frac{\partial \Theta}{\partial t} + \frac{ik\mu}{a} \Theta + (\frac{\partial \Phi}{\partial t} + \frac{i k\mu}{a}\Psi) = n_e \sigma_T[\Theta_0 - \Theta + i\mu v_b - \frac{3\mu^2-1}{4}\Pi]$$ taking $\sigma_T \to 0$ (as this puts the collision term to zero) and just replacing the symbol $\Theta$ with $\mathcal{N}$ gives us the neutrino equation $$\frac{\partial \mathcal{N}}{\partial t} + \frac{ik\mu}{a} \mathcal{N} + (\frac{\partial \Phi}{\partial t} + \frac{i k\mu}{a}\Psi) = 0$$ The neutrino perturbation are soley affected by gravitational redshift and similar gravitational effects - they basially just free-stream from neutrino decoupling till today. The absence of source terms on the right hand side means we don't have the tight-coupling of neutrinos to baryons so the neutrino distribution is not driven to a monopole term as we saw would happen for photons.

What about massive neutrinos? We know from solar oscillation experiments that neutrinos have mass. Neutrinos are fully relativistic in the early Universe (for which the same treatment above applies) and, unless their mass is very low, they will be fully non-relativistic close to the present time (for which the cold dark matter treatment applies) however in the transition region we have do to do things properly. In this case $\mathcal{N} = \mathcal{N}(t,k,p,\hat{p})$ will have a $p$ dependence so we will get extra equations for each value of $p$ that we have to solve and from this we have to compute the Boltzmann integral over the distribution function to get the density perturbation needed in the Poisson equation and so on. This makes it more complicated to do numerically so we will simply skip it in this course and stick with massless neutrinos (for PhD students; master students can ignore neutrinos alltogether), but you can see the classical paper by Chung-Pei Ma and Edmund Bertschinger for more info on how to do this.

Multipole expansions

We have seen that the equation for the photon temperature perturbation (and also for the baryons) depends on angular integrals: $\Theta_0$ and $\Theta_1$. To get a simpler set of equations its convenient to expand the $\mu$-dependence of $\Theta$ in Legendre multipoles (a more convenient "Taylor" series - more convenient as its closely related to spherical harmonics which is the natural choice of basis functions for describing the CMB that we actually observe: as function of direction on the sky). That is we write (the factor $(2\ell+1)/i^\ell$ is just a convention) $$\Theta(t,k,\mu) = \sum \frac{2\ell + 1}{i^\ell}\Theta_\ell(t,k) P_\ell(\mu)$$ where $P_\ell$ is the Legendre polynomials and $\Theta_\ell$ are called the multipoles of the photon distribution and are given by $$\Theta_\ell = \frac{i^\ell}{2}\int_{-1}^1 \Theta(t,k,\mu) P_\ell(\mu)d\mu$$ The moments we have encountered so far are 1) the monopole $\Theta_0$ which (as we shall see) is nothing but the density perturbation $\delta_\gamma = 4\Theta_0$. 2) the dipole $\Theta_1$ which is nothing but the photon velocity $v_\gamma = - 3\Theta_1$ 3) the quadrupole $\Theta_2$ (the $\Pi$ term). Higher order moments describe additional anisotropies in the photon distribution (that are not always easy to give a clear physical intuition off). We are now ready to take moments of the Boltzmann equation for photons. For this we need the orthogonality relation for the Legendre multipoles, $P_0=1$, $P_1=\mu$, $P_2 = \frac{3\mu^2-1}{2}$ and Bonnet's recursion formula (for $\ell \gt 0$, but this also holds for $\ell = 0$ if we take $P_{-1} = 0$) $$\mu P_\ell = \frac{\ell+1}{2\ell+1}P_{\ell+1} + \frac{\ell}{2\ell+1}P_{\ell-1}$$ Multipying the photon equation by $\frac{i^\ell}{2} P_\ell$ and I have rewitten $1 = P_0$ and $\mu = P_1$ to more clearly see how the orthogonality formula then can be used on each term $$\frac{\partial}{\partial t} \left(\frac{i^\ell}{2}\Theta P_\ell\right) + \frac{ik}{a} \frac{i^\ell}{2} \Theta \mu P_\ell + (\frac{\partial \Phi}{\partial t} \frac{i^\ell}{2}P_0P_\ell + \frac{i k}{a}\Psi \frac{i^\ell}{2} P_1P_\ell) = n_e \sigma_T[\Theta_0 \frac{i^\ell}{2} P_0 P_\ell - \frac{i^\ell}{2} \Theta P_\ell + ik v_b \frac{i^\ell}{2} P_1 P_\ell - \frac{1}{2}\Pi \frac{i^\ell}{2} P_2 P_\ell]$$ and integrating over $\mu$ we obtain (note that $\Phi,\Psi,\Pi,\Theta_0$ does not depend on $\mu$ so can be taken outside the integral) $$\frac{\partial}{\partial t} \Theta_\ell + \frac{ik}{a} \frac{i^\ell}{2}\int_{-1}^1 \Theta \mu P_\ell d\mu + (\frac{\partial \Phi}{\partial t} \delta_{\ell 0} - \frac{k}{3a}\Psi \delta_{\ell 1}) = n_e \sigma_T[\Theta_0 \delta_{\ell 0} - \Theta_\ell - \frac{v_b}{3} \delta_{\ell 1} - \frac{1}{2}\Pi \delta_{\ell 2}]$$ and using Bonnet's formula on the second term we get $$\frac{\partial}{\partial t} \Theta_\ell + \frac{k}{a}[\frac{\ell+1}{2\ell+1}\Theta_{\ell+1} - \frac{\ell}{2\ell+1}\Theta_{\ell-1} ] + (\frac{\partial \Phi}{\partial t} \delta_{\ell 0} - \frac{k}{3a}\Psi \delta_{\ell 1}) = n_e \sigma_T[\Theta_0 \delta_{\ell 0} - \Theta_\ell - \frac{v_b}{3} \delta_{\ell 1} - \frac{1}{2}\Pi \delta_{\ell 2}]$$ Writing this out for $\ell=0,1,\ldots$ we see that we get a hierarchy of equations where the equation for $\Theta_\ell$ is sourced by both smaller and the next multipole above: $$ \begin{align} \Theta^\prime_0 &= -\frac{ck}{\mathcal{H}} \Theta_1 - \Phi^\prime, \\ \Theta^\prime_1 &= \frac{ck}{3\mathcal{H}} \Theta_0 - \frac{2ck}{3\mathcal{H}}\Theta_2 + \frac{ck}{3\mathcal{H}}\Psi + \tau^\prime\left[\Theta_1 + \frac{1}{3}v_b\right], \\ \Theta^\prime_\ell &= \frac{\ell ck}{(2\ell+1)\mathcal{H}}\Theta_{\ell-1} - \frac{(\ell+1)ck}{(2\ell+1)\mathcal{H}} \Theta_{\ell+1} + \tau^\prime\left[\Theta_\ell - \frac{1}{10}\Pi \delta_{\ell,2}\right], \quad\quad 2 \le \ell \end{align} $$ where derivatives are now with respect to $x=\log a$, we have reintroduced $c$ and we have used the definition of the optical depth $\tau^\prime = -\frac{cn_e\sigma_T}{H}$. The first equation above is nothing but the perturbed continuity equation and describes mass-energy conservation as can be seen by writing it on the more suggestive form (remember $\delta_\gamma = 4\Theta_0$ and $v_\gamma = -3\Theta_1$) $$\delta_\gamma^\prime = \frac{4ck}{3\mathcal{H}} v_\gamma - 4\Phi^\prime$$ which has the standard form $\delta^\prime = (1+w)(\frac{ck}{\mathcal{H}} v - 3\Phi^\prime)$ for the perturbed continuity equation. We could have gotten this equation much easier from simply perturbing the energy momentum tensor directly (see e.g. this problem where we do that calculation for dark matter) instead of going through the full Boltzmann formalism (but that is not true for the higher order moments so it was not it vein), but it serves as a nice double check that what we get makes sense. The second equation is the pertubed Euler equation which we can see by writing it on the more suggestive form $$v_\gamma^\prime = -\frac{ck}{4\mathcal{H}} \delta_\gamma + \frac{2ck}{\mathcal{H}}\Theta_2 - \frac{ck}{\mathcal{H}}\Psi + \tau^\prime\left[v_\gamma - v_b\right]$$ The term $\Theta_2$ is the quadrupole of the photon distribution and the latter term represents momentum transfer between photons and baryons. This last term ensures that early on when Compton scattering is efficient then $v_b = v_\gamma$, i.e. baryons and photons move together as on fluid (which we call tight-coupling). What makes photons different from baryons and CDM is that they have higher order moments that need to be taken into account and to get these we really had to do the calculations as we did above - there is no shortcut. But we have an infinite amount of equations, isn't this a problem? In pratice we cannot solve an infinite amount of equations so we must truncate the hierarcy and luckily the higher order moments are small. A naive cutoff (just setting high orders to zero) is not a good idea, we must modify the equation for the highest multipole to avoid this screwing up the solution for the smaller multipoles. See Callin for how to best do this. The punchline is that we only include multipoles up to some given $\ell_{\textrm{max}}$ and for this last multipole the equation is taken to be $$ \Theta_{\ell}^\prime = \frac{ck}{\mathcal{H}} \Theta_{\ell-1}-c\frac{\ell+1}{\mathcal{H}\eta(x)}\Theta_\ell+\tau^\prime\Theta_\ell, \quad\quad \ell = \ell_{\textrm{max}} $$ This now represents a closed system of equations for the $\Theta_\ell$'s given the other perturbations. In practice (once we have introduced line of sight integration) we can get away with $\ell_{\textrm{max}} \sim 6$ or for high precision (say $0.1\%$ accuracy) $\ell_{\textrm{max}} \sim 10$ is more than enough. What about neutrinos? We don't have to do any work for this. Massless neutrinos and photons evolve the same apart from the interations of photons with baryons so simply take the photon equations and remove the terms that contains the optical depth $\tau^\prime$, change $\Theta \to \mathcal{N}$ and voila you have the neutrino equations. For solving it numerically we need to include a few more multipoles (as for neutrinos we have no tight-coupling where the higher order multipoles are driven to zero). Including something like $\sim 10-12$ multipoles should be good enough.

Einstein equations

Having found the perturbed Boltzmann equations the final piece missing is equations tells us exactly what the metric potentials $\Phi$ and $\Psi$ are. This is just some brute-force calculations akin to what we did when computing the Friedmann equations. We follow the same steps starting with the metric $$\mathrm{d} s^{2}= -dt^2(1+2\Psi) + a^2(1+2\Phi)(dx^2 + dy^2 + dz^2)$$ where $g_{00} = -(1+2\Psi)$ and $g_{ij} = a^2(1+2\Phi)\delta_{ij}$. The metric is diagonal so the inverse is trival to compute: $g^{00} = - \frac{1}{1+2\Psi} \simeq -1 + 2\Psi$ and $g^{ij} = \frac{1}{a^2(1+2\Phi)} \delta^{ij} \simeq a^{-2}(1-2\Phi)\delta_{ij}$. Next we compute the Christoffel symbols and some of these we have already computed when evaluating the geodesic equation. The only difference with respect to the calculation we did for the background is that $\Phi,\Psi$ now depends on both time $t$ and position $x^i$. The non-zero coefficients are $$\Gamma^0_{00} = \dot{\Psi}$$ $$\Gamma^0_{0i} = \Psi_{,i}$$ $$\Gamma^0_{ij} = \delta_{ij}a^2(H + 2H(\Phi-\Psi) + \dot{\Phi})$$ $$\Gamma^i_{00} = \frac{1}{a^2}\Psi_{,i}$$ $$\Gamma^i_{0j} = \delta_{ij}(H + \dot{\Phi})$$ $$\Gamma^i_{jk} = \delta_{ij}\Phi_{,k} + \delta_{ik}\Phi_{,j} - \delta_{jk}\Phi_{,i}$$ From this we compute the Ricci tensor $$R_{00} = -3\frac{\ddot{a}}{a} - \frac{k^2}{a^2}\Psi - 3\ddot{\Phi} + 3H(\dot{\Psi} - 2\dot{\Phi})$$ $$R_{0i} = 0$$ $$R_{ij} = \delta_{ij}[(2a^2H^2 + a \ddot{a})(1+2\Phi - 2\Psi) + a^2(6\dot{\Phi} - \dot{\Psi}) + a^2\ddot{\Phi} - \nabla^2\Phi] - (\Phi + \Psi)_{,ij}$$ and the Ricci scalar $$R = g^{\mu\nu} R_{\mu\nu} = 6(H^2 + \frac{\ddot{a}}{a})(1-2\Psi) - 2\frac{\nabla^2\Psi}{a^2} + 6\ddot{\Phi} - 6(\dot{\Psi} - 4\dot{\Phi}) - 4\frac{\nabla^2\Phi}{a^2}$$ and the Einstein tensor $$\delta G^0_0 = -6H\dot{\Phi} + 6\Psi H^2 + \frac{2\nabla^2\Phi}{a^2}$$ The rest of the components are left as an exercise. Before we continue on to evaluating the Einstein equations (for which there are some tricks to extract the equations we need) we need to look at the right hand side of the Einstein equation.

The general relativistic expression for the energy momentum tensor in the Boltzmann formalism is $$ T^\mu_\nu = \frac{g}{(2\pi)^3} \int \frac{dP_1 dP_2 dP_3}{\sqrt{-\det g}} \frac{P^\mu P_\nu}{P^0} f $$ where $P^\mu = \frac{dx^\mu}{d\lambda}$ is the 4-momentum and $\det g$ is the determinant of the metric $g_{\mu\nu}$. Since the metric is diagonal we simply have $$-\det g = a^6(1+2\Psi)(1+2\Phi)^3 \simeq a^6(1+2\Psi+6\Phi)$$ The expressions for $P^\mu$ we know from the treatment above where we found $P^i = \frac{\vec{p}^i}{a}(1-\Phi)$ and therefore $P_i = a\vec{p}_i (1+\Phi)$ we therefore have the Jacobian $\left|\frac{\partial \vec{P}}{\partial \vec{p}}\right| = a^3(1+3\Phi)$. For the $00$ component this gives us the result we already know $$ T^0_0 = -\frac{g}{(2\pi)^3} \int d^3\vec{p} E f = - \rho$$ and with $f = \overline{f} + \delta f$ we get $$ T^0_0 = \overline{T^0_0} - \frac{g}{(2\pi)^3} \int d^3\vec{p} E \delta f $$ where $\overline{T^0_0} = -\overline{\rho}$. For photons $\delta f = -p\Theta \frac{\partial \overline{f}}{\partial p}$ so the integral above is $$ \frac{g}{(2\pi)^3} \int d^3\vec{p} E \delta f = -\frac{g}{(2\pi)^3} \int d^3\vec{p} \Theta p^2\frac{\partial \overline{f}}{\partial p}$$ Now using spherical coordinates we have $$ \frac{g}{(2\pi)^3} \int d^3\vec{p} E \delta f = -\frac{g}{(2\pi)^3} \int d\Omega \int dp \Theta p^4\frac{\partial \overline{f}}{\partial p}$$ Remembering that $\Theta$ does not depend on $p$ we can do an integration by parts to get $$ \frac{g}{(2\pi)^3} \int d^3\vec{p} E \delta f = \frac{g}{(2\pi)^3} 4\int d\Omega \Theta \int dp p^3 \overline{f}$$ The angular integral is just the definition of the monopole $\Theta_0 = \frac{1}{4\pi}\int d\Omega \Theta$ and what remains $\int dp 4\pi p^2 \cdot p \overline{f} = \int d^3 p p \overline{f}$ is nothing but the energy density in the background giving us $$ \frac{g}{(2\pi)^3} \int d^3\vec{p} E \delta f = -4\Theta_0 \overline{\rho}_\gamma$$ and it follows that $$ (T^0_0)_\gamma = -\overline{\rho}_\gamma(1+4\Theta_0)$$ This is the justification of the interpretation of the monopole that we mentioned previously without proof that $\frac{\delta \rho_\gamma}{\overline{\rho}_\gamma} = 4\Theta_0$. For massless neutrinos its exactly the same calculation so $$ (T^0_0)_\nu = -\overline{\rho}_\nu(1+4\mathcal{N}_0)$$ and for baryons and cold dark matter we don't have to do any calculations as we defined $\delta_b$ and $\delta_{\rm CDM}$ to be the density perturbation so $$(T^0_0)_b = -\overline{\rho}_b(1 + \delta_b)$$ $$(T^0_0)_{\rm CDM} = -\overline{\rho}_{\rm CDM}(1 + \delta_{\rm CDM})$$ Note that we derived this is real-space though exactly the same expression holds in Fourier space - the beauty of linear perturbation theory is that we can just jump back and forth and don't have to be super careful about this.

Summary

We have now derived the full set of Boltzmann equations for all the species. This includes a continuity equation ($\delta_b', \delta_{\rm CDM}', \delta_\gamma', \delta_\nu'$), an Euler equation ($v_b', v_{\rm CDM}', v_\gamma', v_\nu'$) and a set of equations for the higher multipoles for photons and neutrinos. Expressed in terms of the time-variable $x=\log a$ these can be written:

Photon temperature multipoles ($\delta_\gamma = 4\Theta_0$ and $v_\gamma = -3\Theta_1$): $$ \begin{align} \delta_\gamma^\prime &= \frac{4}{3}\left(\frac{ck}{\mathcal{H}} v_\gamma - 3\Phi^\prime\right), \\ v_\gamma^\prime &= -\frac{ck}{4\mathcal{H}} \delta_\gamma + \frac{2ck}{\mathcal{H}}\Theta_2 - \frac{ck}{\mathcal{H}}\Psi - \tau^\prime\left[v_b - v_\gamma\right], \\ \Theta^\prime_\ell &= \frac{\ell ck}{(2\ell+1)\mathcal{H}}\Theta_{\ell-1} - \frac{(\ell+1)ck}{(2\ell+1)\mathcal{H}} \Theta_{\ell+1} + \tau^\prime\left[\Theta_\ell - \frac{1}{10}\Pi \delta_{\ell,2}\right], \quad\quad 2 \le \ell \end{align} $$ Neutrino multipoles ($\delta_\nu = 4\mathcal{N}_0$ and $v_\nu = -3\mathcal{N}_1$): $$ \begin{align} \delta_\nu^\prime &= \frac{4}{3}\left(\frac{ck}{\mathcal{H}} v_\nu - 3\Phi^\prime\right), \\ v_\nu^\prime &= -\frac{ck}{4\mathcal{H}} \delta_\nu + \frac{2ck}{\mathcal{H}}\mathcal{N}_2 - \frac{ck}{\mathcal{H}}\Psi \\ \mathcal{N}^\prime_\ell &= \frac{\ell ck}{(2\ell+1)\mathcal{H}} \mathcal{N}_{\ell-1} - \frac{(\ell+1)ck}{(2\ell+1)\mathcal{H}}\mathcal{N}_{\ell+1},\quad\quad 2 \le \ell \end{align} $$ Cold dark matter: $$ \begin{align} \delta_{\rm CDM}^\prime &= \frac{ck}{\mathcal{H}} v_{\rm CDM} - 3\Phi^\prime \\ v_{\rm CDM}^\prime &= -v_{\rm CDM} -\frac{ck}{\mathcal{H}} \Psi \end{align} $$ Baryons: $$ \begin{align} \delta_b^\prime &= \frac{ck}{\mathcal{H}}v_b -3\Phi^\prime \\ v_b^\prime &= -v_b - \frac{ck}{\mathcal{H}}\Psi + \tau^\prime R(v_b - v_\gamma) \end{align} $$

On top of this we also have a set of equations of the polarization of photons that we did not derive (but for which PhD students will implement in the numerical project): $$ \begin{align} \Theta^\prime_{P0} &= -\frac{ck}{\mathcal{H}}\Theta_{P1} + \tau^\prime\left[\Theta_{P0} - \frac{1}{2}\Pi \right]\\ \Theta_{P\ell}^\prime &= \frac{\ell ck}{(2\ell+1)\mathcal{H}} \Theta_{\ell-1}^P - \frac{(\ell+1)ck}{(2\ell+1)\mathcal{H}} \Theta_{\ell+1}^P + \tau^\prime\left[\Theta_\ell^P - \frac{1}{10}\Pi\delta_{\ell,2}\right],\quad\quad 1 \le \ell \end{align} $$

The system above is closed by the Einstein equations which tell us what the metric potentials are the Poisson equation for $\Phi$ and an equation for the anisotropic stress $\Phi+\Psi$: $$ \begin{align} \Phi^\prime &= \Psi - \frac{c^2k^2}{3\mathcal{H}^2} \Phi + \frac{1}{2} \left[\Omega_{\rm CDM}(a) \delta_{\rm CDM} + \Omega_{b}(a) \delta_b + \Omega_{\gamma}(a)\delta_\gamma + \Omega_{\nu}(a)\delta_\nu\right] \\ \Psi &= -\Phi - \frac{12\mathcal{H}^2}{c^2k^2}\left[\Omega_{\gamma}(a)\Theta_2 + \Omega_{\nu}(a)\mathcal{N}_2\right] \\ \end{align} $$

In the equations above $\Pi = \Theta_2 + \Theta_0^P + \Theta_2^P$ and $R = \frac{4\Omega_{\gamma 0}}{3\Omega_{b 0} a}$ (note that our $R$ is $1/R$ in Dodelson).

The continuity equations above takes the general form $$\delta' = (1+w)\left(\frac{ck}{\mathcal{H}}v -3\Phi^\prime\right)$$ where $w$ is the equation of state and the Euler equations (here without the collision term and anisotropic stress) takes the general form $$v' = -(1-3c_s^2)v - \frac{wck}{(1+w)\mathcal{H}}\delta - \frac{ck}{\mathcal{H}}\Psi$$ where $c_s^2$ is the sound speed. For baryons and CDM $w=c_s^2 = 0$ and for photons and neutrinos $c_s^2=w=1/3$ which gives the equations above (without the collision term and the $\Theta_2$,$\mathcal{N}_2$ terms). If it wasn't for the photons and neutrino distribution function having higher order moments and the collision term we could have derived these directly by perturbing the energy momentum tensor and just using the conservation and Einstein equations to derive these.

We are almost ready to solve this system numerically, however these are a few things missing. How do we deal with having an infinite hierarchy of equations for photons / neutrinos / polarization and what are the initial conditions? The former question turns out to be fairly straight forward: the higher order multipoles goes to zero quite rapidly so we only need to include a finite amount (and modify the equations for the largest multipole slightly). For the latter question we need to talk about the theory of inflation which is the topic we will cover next.