Table of contents

Back to the main page 1. Background cosmology

2. Thermodynamics and statistical physics 3. CMB angular power-spectrum

Background Cosmology

Here are some notes on the mathematics that goes into background cosmology. For more details see Chapter 1 and Chapter 3 in Baumann.

Derivation of the Friedman equations

The brute-force algorithm for getting to the Einstein equations are as follows. Start with the metric $ds^2 = g_{\mu\nu}dx^\mu dx^\nu$ (whose form follows from the assumption about homogenity and isotropy) in a given coordinate system. We assume the curvature is zero in this derivation for which doing the derivation in Cartesian coordinates $(t,x,y,z)$ is simplest. We have $g_{00} = -1$ and $g_{ij} = a^2(t)\delta_{ij}$. First compute the inverse metric $g^{\mu\nu}$. Since the metric is diagonal this is easy $g^{00} = -1$ and $g^{ij} = a^{-2}(t) \delta^{ij}$. Next compute the Christoffel symbols $$\Gamma^\mu_{\alpha\beta} = \frac{g^{\mu\delta}}{2}(g_{\delta\beta,\alpha} + g_{\alpha\delta,\beta} - g_{\alpha\beta,\delta})$$ To do this systematically divide it into cases where $\mu = 0$ (time index) and $\mu = i$ (a space index). Then notice that since $g$ is diagonal only the term with $\mu = \delta$ will be non-zero. This simplifies the expression to $$\Gamma^0_{\alpha\beta} = \frac{g^{00}}{2}(g_{0\beta,\alpha} + g_{\alpha 0,\beta} - g_{\alpha\beta,0})$$ $$\Gamma^i_{\alpha\beta} = \frac{g^{ii}}{2}(g_{i\beta,\alpha} + g_{\alpha i,\beta} - g_{\alpha\beta,i})$$ and note that there is no implicit summation over $i$ in the last expression. Since the Christoffel symbol is symmetric in the two lower indices we have three cases for each of the two terms: $\alpha\beta=00$ (time-time), $\alpha\beta = 0j$ (one time-one space) and $\alpha\beta = jk$ (space-space). What further simplifies the derivation is that the metric only depends on $t$ so any space derivative ($,i$) vanishes. This gives you six different symbols to compute and only two of these are non-zero namely $$\Gamma^0_{ij} = a\dot{a}\delta_{ij}$$ $$\Gamma^i_{0j} = \Gamma^i_{j0} = \frac{\dot{a}}{a}\delta^i_j$$ Having computed these we can go on to evaluate the Ricci tensor which is given by $$R_{\mu\nu} = \Gamma^\alpha_{\mu\nu,\alpha} - \Gamma^\alpha_{\mu\alpha,\nu} + \Gamma^\alpha_{\mu\nu}\Gamma^\beta_{\alpha\beta} - \Gamma^\beta_{\mu\alpha}\Gamma^\alpha_{\beta\nu}$$ To evaluate this we again split into two cases $R_{00}$ and $R_{ij}$ (there is also $R_{0j}$, but since the metric is diagonal this is forced to be zero). This gives us $$R_{00} = -3\frac{\ddot{a}}{a}$$ $$R_{ij} = (a\ddot{a} + 2\dot{a}^2)\delta_{ij}$$ From this we can compute the Ricci scalar $$R = g^{\mu\nu}R_{\mu\nu} = g^{00}R_{00} + g^{ij}R_{ij} = 6\left(\frac{\ddot{a}}{a} + \frac{\dot{a}^2}{a^2}\right)$$ which finally gives us the Einstein tensor $G_{\mu \nu} = R_{\mu\nu} - \frac{1}{2}g_{\mu\nu}R$. We can go ahead and compute all the different component of this. The $G_{0i}$ component will be zero, the $00$ component gives $$G_{00} = 3\frac{\dot{a}^2}{a^2}$$ and the $ij$ component becomes $$G_{ij} = \left(-2\ddot{a}a - \dot{a}^2\right)\delta_{ij}$$ We are now in the positon to evaluate the Einstein equation $G_{\mu\nu} = 8\pi G T_{\mu\nu}$. Taking the energy content of the Universe to consist of the components $(n)$ (denoting baryons, cold dark matter, radiation, etc.) then $T_{\mu\nu} = \sum_n T_{\mu\nu}^{(n)}$ where each component $(n)$ is assumed to be a perfect fluid $$T_{\mu\nu}^{(n)} = (\rho_n + p_n)u_\mu u_\nu + p_ng_{\mu\nu}$$ where $u_\mu$ is the 4-velocity of the fluid which in our co-moving frame is simply $u^\mu = (1,0,0,0)$ thus $T_{00}^{(n)} = \rho_n$ and $T_{jk}^{(n)} = p_n a^2\delta_{ij}$. The $00$ component of the Einstein equations then gives us the Hubble equation $$H^2 = \frac{8\pi G}{3} \sum_n \rho_n$$ where $H \equiv \frac{\dot{a}}{a}$. The $11$, $22$ and $33$ components gives us 3 more equations (because of isotropy all of these are equal) which reads $$\frac{\ddot{a}}{a} = -\frac{4\pi G}{3}\sum_n (\rho_n + 3p_n)$$ In addition to this we have the conservation equation $\nabla_\mu T^{\mu\nu} = 0$, which in the absence of any interactions between the different species simply reduces to one equation for each component $\nabla_\mu T^{\mu\nu}_{(n)} = 0$ whose $0$ component gives us $$\nabla_\mu T^{\mu 0}_{(n)} = \partial_\mu T^{\mu 0}_{(n)} + \Gamma^\mu_{\mu\alpha} T^{\alpha 0}_{(n)} + \Gamma^0_{\mu\alpha} T^{\mu \alpha}_{(n)} \implies \dot{\rho}_n + 3H(\rho_n + p_n) = 0$$

Thermodynamics and statistical physics

Here is a brief summary of the thermodynamics / statistical physics we need in this course. For more details see Baumann (for equilibrium and close to equilibirum thermodynamics in an expanding Universe) and Dodelson (for the out of equilibirum thermodynamics and the treatment of the Boltzmann equation).

Distribution function

The distribution function tells us how many particles we have in a small phase space volume $$\Delta N = \frac{g}{(2\pi)^3} (\Delta \vec{x})^3 (\Delta \vec{p})^3$$ where $g$ is the number of internal degrees of freedom (spin/polarization). We can obtain the usual macroscopic quantities (number-density, energy density, pressure) as weigthed integrals over the distribution function: $$n = \frac{g_i}{(2\pi)^3}\int f d^3p$$ $$\rho = \frac{g_i}{(2\pi)^3}\int Ef d^3p$$ $$P = \frac{g_i}{(2\pi)^3}\int \frac{p^2}{3E}f d^3p$$ In thermal equilibrium the distribution function takes a simple form $$f = \left(e^{\frac{E-\mu}{T}} \pm 1\right)^{-1}$$ where $E =\sqrt{p^2 + m^2}$ is the relativistic energy relation, $\mu$ is the chemical potential, $T$ is the temperature of the system and $+$ is for bosons (Bose-Einstein distribution) and $-$ is for fermions (Pauli-Dirac distribtion). The exact form of the chemical potential (in general temperature dependent) will not play a big role in this course, but its important to know that if we have an interaction $1+2\leftrightarrow 3+4$ then in chemical equilibrium the chemical potential of each of the species involved has to satisfy $\mu_1 + \mu_2 = \mu_3 + \mu_4$ and a consequence of this is that $\mu = 0$ for photons (since photon numbers are not conserved). In the low temperature limit both these distributions approaches the Maxwell distribition $$f \approx e^{-\frac{E-\mu}{T}}$$ We can from this compute what $n,P,T$ are as functions of temperature. For photons ($m=0$) the integrals can be evaluated directly while for massive particles we can only compute it analytically (in terms of elementary functions) in the limits $T\gg m$ and $T\ll m$. Going beyond thermal equilibrium we need to solve the Boltzmann equation.

This was a very brief summary: see Baumann for a more complete discussion.

Boltzmann equation in general relativity

The Boltzmann equation tells us how the distribution function of some particle species evolves. This is a beast: it depends on time $t$, position $\vec{x}$ and momentum $\vec{p}$ so a seven dimensional function. Not always easy to work with, but its a very general framework that allows us to construct a macroscopic description from knowledge of the microscopic interactions that takes place.

The general relativistic Boltzmann equation reads (note that Dodelson simplifies the discussion of this by phrasing it dirrctly in term of $\frac{df}{dt}$. The difference to this is a factor $\frac{dt}{d\lambda} = P^0 \sim p$ which he has included in his definition of the collision term which has a prefactor $1/p$ compared to what we have here) $$\frac{df}{d\lambda} = C(f)$$ where $\lambda$ is an affine parameter. 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 that satisfy $g_{\mu\nu}P^\mu P^\nu = -m^2$ and follows the geodesic equation $$\frac{dP^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta}P^\alpha P^\beta = 0$$ 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$. We also have the relativistic energy relation $E = \sqrt{p^2 + m^2}$. We will choose to work with $E$ and $\hat{p}$ instead of $\vec{p}$. Both ways are equivalent, and the reason for this split is convenience. Note that for photons we simply have $E = p$.

The next thing we need to know is how particles move in general relativity. This is given by the constraint $g_{\mu\nu}P^\mu P^\nu = -m^2$ together with the geodesic equation. For a given metric this first relation can be used to compute $P^0$ and $P^i$ in terms of $E,p,\hat{p}^i$ and the metric components (for example for a diagonal metric we have $P^0 = E(-g_{00})^{-1/2}$ and $P^i = p\hat{p}^i(-g_{ii})^{-1/2}$). The distribution function depends on $(t,\vec{x},E,\hat{p})$ so we can expand the left hand side of the Boltzmann equation as $$\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}$$ which gives $$\frac{df}{d\lambda} = \frac{\partial f}{\partial t} P^0 + \frac{\partial f}{\partial x^i} P^i + \frac{\partial f}{\partial E}\frac{dE}{d\lambda} + \frac{\partial f}{\partial \hat{p}^i}\frac{d\hat{p}^i}{d\lambda}$$ and finally $$\frac{df}{d\lambda} = P^0\left[\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} + \frac{\partial f}{\partial \hat{p}^i}\frac{1}{P^0}\frac{d\hat{p}^i}{d\lambda}\right]$$ Theterms $\frac{dE}{d\lambda}$ and $\frac{d\hat{p}^i}{d\lambda}$ follows from the geodesic equation. Now we can see the advantage of the choosen set of coordinates: once we look at this equation to first order in perturbation theory the last term won't contibute as it is a second order term (photons don't change direction if there are no perturbations so this is a first order term and the distribution function in the smooth background does not depend on direction so this is also a first order term).

The geodesic equation for $\mu = 0$ $$\frac{dP^0}{d\lambda} + \Gamma^0_{\alpha\beta} P^\alpha P^\beta = 0$$ used with the expression for $P^0$ in terms of $E$ gives us an equation for $\frac{dE}{d\lambda}$ and (not needed to first order in perturbation theory) the $i$'th component $$\frac{dP^i}{d\lambda} + \Gamma^i_{\alpha\beta} P^\alpha P^\beta = 0$$ gives us an equation for $\frac{d\hat{p}^i}{d\lambda}$. For metric perturbations one needs to use that $\frac{dA(t,x)}{d\lambda} = \frac{\partial A(t,x)}{\partial t}P^0 + \frac{\partial A(t,x)}{\partial x^i}P^i$ which is just the chain rule together with the definition of $P^\mu$.

The last ingredient we need to understand is the collision term. This is where most of the complications of the Boltzmann equation enters. The most common process we will enocunter is a $1+2\leftrightarrow 3+4$ process (e.g. $e^- + \gamma \leftrightarrow e^- + \gamma$ or $p^+ + e^- \leftrightarrow H + \gamma$) and for this it reads $$\begin{align} C(f_1) &= \int \frac{d^3p_2}{(2\pi)^32E_2}\int \frac{d^3p_3}{(2\pi)^32E_3}\int \frac{d^3p_4}{(2\pi)^32E_4} \times\\ & \times |\mathcal{M}|^2 (2\pi)\delta(E_1 + E_2 - E_3 - E_4)(2\pi)^3 \delta^{(3)}(\vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4) \times \\ & \times[f_3f_4(1\pm f_1)(1\pm f_2) - f_1f_2 (1\pm f_3)(1\pm f_4)] \end{align} $$ where $E_i = \sqrt{m_i^2 + p_i^2}$, $f_i = f_i(t,\vec{x},E_i,\hat{p}_i)$ and the $\pm$ is a $+$ for bosons and $-$ for fermions (these factors represents quantum effects like stimulated emission and Pauli blocking respectively - the latter being the fact that two fermions cannot be in the same quantum state). To first order in perturbation theory these quantum effects will not be important and we can approximate $1\pm f_i \approx 1$. Note the form of the equation above: an integral over all the momenta involved in the process (apart from the momenta $p_1$ that we are studying), delta functions enforcing energy-momentum conservation, the Feynman amplitude $|\mathcal{M}|^2$ describing the probabillity/efficientcy of the interaction (this follows from quantum field theory calculations) and the two terms $f_3f_4$ and $-f_1f_2$ describing the interaction going to the left and right respectively (creation of particles $1$ and $2$ and destruction of particles $1$ and $2$). It's a beast of an equation, but its quite intuitiven and from this form its not hard to see how this should be generalized to other types of processes (or to different dimensions for that matter).

This was a very brief summary: see Dodelson for a more complete discussion.

Moments of the Boltzmann equation

The distribution function provides a complete statistical description of the system in question, but can seem a bit abstract. The thermodynamical quantities we have much better intuition about like pressure, velocity, energy- and number-density can be computed as weigthed integrals of the distribution function. For example number-density, energy-density, pressure and velocity is given as follow $$n = \frac{g_i}{(2\pi)^3}\int f d^3p$$ $$\rho = \frac{g_i}{(2\pi)^3}\int Ef d^3p$$ $$P = \frac{g_i}{(2\pi)^3}\int \frac{p^2}{3E}f d^3p$$ $$n\vec{v} = \frac{g_i}{(2\pi)^3}\int \frac{p\hat{p}}{E} f d^3p$$ where $g_i$ is the number of internal degrees of freedom (spin/polarization states) each particle species has (e.g. $2$ for both photons, protons and electrons and $2\times 2 = 4$ for hydrogen). Multiplying the Boltzmann equation by certain combinations of $p$ and $\hat{p}^i$ and then integrating is called taking moments of the Boltzmann equation (this is similar to when we take moments of a random distribution in statistics). The zeroth moment - i.e. we multiply by $1$ and integrate - gives us the continuity equation. In the absence of interactions this simply reads $$\frac{1}{a^3}\frac{d(n a^3)}{dt} = 0$$ and describes conservation of mass (number of particles). The first order moment is to multiply by the momentum and integrate. For a smooth Universe this gives us nothing (there are no peculiar velocities), but in general this leads to an Euler equation which describes momentum conservation.

The Boltzmann equation has the general property that the equation for one moment will depends on the next moment leading to an infinite set of equations. This requires us to find a way of truncating this infinite hierarchy. Luckily for many cases the higher order moments are small and can be ignored giving us a consistent closed set of equations we can solve. We shall see examples of this below.

Boltzmann equation in a smooth Universe

In a smooth Universe - as long as the particles remain in thermal equilibrium with the thermal bath - the collision term is zero. The distribution function only depends on $p$ and $t$ (isotropy and homogenity = no directional or positional dependence). For photons moving in the FRLW metric show that $$P^0 = E$$ $$P^i = p\hat{p}^i a^{-1}$$ where $E=p$. The geodesic equation reads $$\frac{dE}{d\lambda} = -\Gamma^0_{\alpha\beta} P^\alpha P^\beta$$ Evaluate the right hand side. Show that the Boltzmann equation becomes $$\frac{df}{d\lambda} = P^0\left[\frac{\partial f}{\partial t} + \frac{\partial f}{\partial E}\frac{1}{E}\frac{dE}{d\lambda}\right] = 0$$ and that this reduces to $$\frac{1}{T}\frac{dT}{dt} = -\frac{1}{a}\frac{da}{dt}$$ Solve this for $T(a)$ and explain the result physically.

This was a very brief summary: see Baumann for a more complete discussion.

Boltzmann equation for the number density in a smooth Universe

For a $1+2\leftrightarrow 3+4$ process the Boltzmann equation in a smooth Universe can be written on the form $$\frac{1}{a^3}\frac{d(n_1 a^3)}{dt} = -\alpha n_1n_2 + \beta n_3n_4$$ where $\alpha = \left<\sigma v\right>$ ($\alpha$ is basically the $|\mathcal{M}|^2$ term in the collision term) is the thermally averaged cross-section. A full derivation gives you $\alpha$ and $\beta$, but we can see what $\beta$ has to be just from studying this equation. If we are in equilibrium then the left hand side vanishes and so must the right hand side. $n_i$ takes their equilibrium values so $$\beta = \alpha \left(\frac{n_1n_2}{n_3n_4}\right)_{\rm eq}$$ and we get $$\frac{1}{a^3}\frac{d(n_1 a^3)}{dt} = -\left<\sigma v\right> \left(n_1n_2 - n_3 n_4 \left(\frac{n_1n_2}{n_3n_4}\right)_{\rm eq}\right)$$

If the interactions are efficient (large $\left<\sigma v\right> n_1n_2$ relative to the characteristic time-scale of changes in $n_i$) then we will be close to equilibrium. This leads us to the Saha approximation $$\frac{n_1n_2}{n_3n_4} \approx \left(\frac{n_1n_2}{n_3n_4}\right)_{\rm eq}$$ The right hand side we know how to compute in terms of the temperature and the particle masses. This equation tells us that each of the $n_i$'s may change (due to the interactions), but this particular combination must stay roughly constant. The most relevant interaction we will encounter in this course is $e^- + p^+ \leftrightarrow H + \gamma$ (relevant for recombination) and this equation will tell us how the number density of free electrons will evolve in time. When the Saha approximation is no longer valid then we have no option other than solving the full Boltzmann equation above. This is not so hard: its a simple ODE, however the complications comes about from the different forms of interactions we can have. The Universe don't just have hydrogen: there is plenty of Helium and some Litium also. On top of this there is not just one hydrogen: the electron can be in many different states with different energy and there will be transitions between the different states etc. This makes this whole study messy.

CMB angular power-spectrum

Spherical harmonics

Spherical harmonics $Y_{\ell m}$ (the indices here are just indices - its not a tensor and it does not matter if they are up or down) are orthogonal basis functions for functions defined on a sphere. They are eigenfunctions of the square of the angular momentum operator $\hat{L} = {\bf x}\times \nabla$ (i.e. $\hat{L}^2 Y_{\ell m} = \ell(\ell+1)Y_{\ell m}$) and pop up when trying to solve Laplace equation in spherical coordinates using separation of variables. These properties makes them the natural choice when we are going to expand the CMB temperature field. Let's review some properties (this depends on the particular convention one uses - there are different ways to define $Y$'s that lead to slightly different formulas):

An explicit expression for $Y_{\ell m}$ (which we will not need very much): $$Y_{\ell m}(\theta,\phi) = (-1)^m\sqrt{\frac{(2\ell+1)(\ell-m)!}{4\pi(\ell+m)!}}P_{\ell m}(\cos(\theta))e^{im\varphi}$$ The orthogonality conditions: $$\int d\Omega_{\hat{p}} Y_{\ell m} Y^*_{\ell' m'} = \delta_{\ell\ell'}\delta_{mm'}$$ The expansion of Legendre polynomials in spherical harmonics: $$\mathcal{P}_\ell(\hat{a}\cdot \hat{b}) = \frac{4\pi}{2\ell+1}\sum_{m=-\ell}^\ell Y_{\ell m}(\hat{a}) Y^*_{\ell m}(\hat{b})$$ where $P_{\ell m}$ are associated Legendre polynomials. When $m=0$ we have $P_{\ell m} = \mathcal{P}_{\ell}$ the usual Legendre polynomials (the ones we used when expanding the temperature perturbations in the Boltzmann equation in multipoles).

Another useful relation that we will need (which follows from the above is) $$\int d\Omega_{\hat{p}} \mathcal{P}_\ell(\hat{k}\cdot\hat{p}) Y_{\ell'm}^*(\hat{p}) = \delta_{\ell\ell'} \frac{4\pi}{2\ell+1} Y_{\ell m}^*(\hat{k}) $$ A word on notation: we will sometimes write $Y_{\ell m}(\theta,\phi)$ and sometimes $Y_{\ell m}(\hat{n})$ where $\hat{n}$ is some direction (unit vector). These are equivalent and with the latter we simply mean the first where the angles $(\theta,\phi)$ are those associated with $\hat{n}$ (i.e. the angles we get when $\hat{n}$ is written in spherical coordinates).

From perturbations to the angular power-spectrum

For the CMB we measure the temperature perturbations $\Theta = \frac{\delta T}{T}$ as function of direction on the sky (which is the same as saying we have it as a function on the unit sphere). We can expand this in spherical harmonics as $$\Theta(\vec{x},\hat{p},t) = \sum_{\ell=1}^\infty\sum_{m=-\ell}^{\ell} a_{\ell m}(\vec{x},t) Y_{\ell m}(\hat{p})$$ which means that $$a_{\ell m}(\vec{x},t) = \int d\Omega_{\hat{p}} Y_{\ell m}^*(\hat{p})\Theta(\vec{x},\hat{p},t)$$ Expressing $\Theta$ in terms of its Fourier components this becomes $$a_{\ell m}(\vec{x},t) = \int\frac{d^3k}{(2\pi)^3}e^{i\vec{k}\cdot\vec{x}}\int d\Omega_{\hat{p}} Y_{\ell m}^*(\hat{p})\Theta(\vec{k},\hat{p},t)$$ Inflation predicts that the $a_{\ell m}$'s is a (near) gaussian random field with mean $0$ and variance $C_\ell$, i.e. $$\left<a_{\ell m}a^*_{\ell' m'}\right>= \delta_{\ell\ell'}\delta_{mm'}C_\ell$$ where the brackets denotes an ensamble average (average over many different realizations of our Universes). And since its gaussian all the statistical information about $\Theta$ is contained in $C_\ell$ (in practice there will be some deviations from pure gaussianity depending on the model of inflation so higher order moments does potentially contain interesting information). For a fixed $\ell$, each $a_{\ell m}$ have the same variance and for each value of $\ell$ there are $2\ell+1$ possible values of $m$. This means that for the large $\ell$'s we have enough statistical power to accurately estimate $C_\ell$, while for very low $\ell$ we will run into cosmic variance: we have very few $m$'s to estimatee $C_\ell$ from so there is a variance - that is due to us only having one single Universe to measure this in - that is given by $\frac{\text{Var}(C_\ell)}{C_\ell^2} = \frac{2}{2\ell+1}$ (we have $C_\ell = \frac{1}{2\ell +1}\sum_m |a_{\ell m}|^2$ which follows a $\chi^2$ distribution).

To derive an expression for $C_\ell$ we expand $\Theta$ in multipoles $$a_{\ell m}(\vec{x},t) = \sum_\ell (2\ell + 1)(-i)^\ell\int\frac{d^3k}{(2\pi)^3}e^{i\vec{k}\cdot\vec{x}}\int d\Omega_{\hat{p}} Y_{\ell m}^*(\hat{p})\mathcal{P}_\ell(\mu)\Theta_\ell(\vec{k},t)$$ where $\mu = \hat{p} \cdot \hat{k}$ which gives us $$a_{\ell m}a_{\ell' m'}^* = \sum_{\ell_1,\ell_2} (2\ell_1 + 1)(2\ell_2+1)(-i)^{\ell_1-\ell_2}\int\frac{d^3kd^3k'}{(2\pi)^6}e^{i(\vec{k} - \vec{k}')\cdot\vec{x}}\int d\Omega_{\hat{p}} d\Omega_{\hat{p}'} Y_{\ell m}^*(\hat{p}) Y_{\ell'm'}(\hat{p}')\mathcal{P}_{\ell_1}(\mu)\mathcal{P}_{\ell_2}(\mu')\Theta_{\ell_1}(\vec{k},t) \Theta_{\ell_2}^*(\vec{k}',t)$$

Now remember that when we solve this numerically we put $\Theta_\ell(\vec{k},t) = \Theta^{\rm our}_\ell(k,t) \mathcal{R}_{\rm ini}(\vec{k})$ where the initial value of the curvature perturbation $\mathcal{R}_{\rm ini}\propto \Psi$ was set to unity, i.e. we factored out the common initial condition from inflation before solving (which is perfectly fine since the equation system is linear). We must now put it back in. When we take the ensamble average of the above this will lead to a term $\left<\mathcal{R}_{\rm ini}(\vec{k})\mathcal{R}^*_{\rm ini}(\vec{k}')\right> = (2\pi)^3\delta(\vec{k} - \vec{k}')P(k)$ with $P(k) = \frac{2\pi^2}{k^3}\mathcal{P}_\mathcal{R}(k)$ being the primordial power-spectrum (with $\mathcal{P}_\mathcal{R}(k) = A_s(k/k_{\rm pivot})^{n_s-1}$) and the delta-function will enforce $\vec{k} = \vec{k}'$ taking care of one of the integrals. $$\left<a_{\ell m}a_{\ell' m'}^*\right> = \sum_{\ell_1,\ell_2} (2\ell_1 + 1)(2\ell_2+1)(-i)^{\ell_1-\ell_2}\int\frac{d^3k}{(2\pi)^3}\int d\Omega_{\hat{p}} d\Omega_{\hat{p}'} Y_{\ell m}^*(\hat{p}) Y_{\ell'm'}(\hat{p}')\mathcal{P}_{\ell_1}(\mu)\mathcal{P}_{\ell_2}(\mu')\Theta^{\rm our}_{\ell_1}(k,t)\Theta^{\rm our}_{\ell_2}(k,t) \frac{2\pi^2\mathcal{P}_\mathcal{R}(k)}{k^3}$$

For the next step we need $\int d\Omega_{\hat{p}} \mathcal{P}_\ell(\mu) Y_{\ell'm}^*(\hat{p}) = \delta_{\ell\ell'} \frac{4\pi}{2\ell+1} Y^*_{\ell m}(\hat{k})$ which takes care of two of the integrals (and the two Kroneker delta kills both of the sums) leaving us with $$\left<a_{\ell m}a_{\ell' m'}^*\right> = \delta_{\ell\ell'}(4\pi)^2\int\frac{d^3k}{(2\pi)^3}Y_{\ell m}(\hat{k})Y_{\ell m'}^*(\hat{k})|\Theta^{\rm our}_\ell(k,t)|^2 \frac{2\pi^2\mathcal{P}_\mathcal{R}(k)}{k^3}$$

Finally splitting the $k$-integral into an integral over angles ($\Omega_{\hat{k}}$ now means over the angles of $k$) using that $\int d\Omega_{\hat{k}} Y_{\ell m}(\hat{k})Y_{\ell m'}^*(\hat{k}) = \delta_{mm'}$ we get $$\left<a_{\ell m}a_{\ell' m'}^*\right> = \delta_{\ell\ell'}\delta_{mm'}4\pi\int\frac{dk}{k}\mathcal{P}_\mathcal{R}(k)|\Theta^{\rm our}_\ell(k,t)|^2$$ so $$C_\ell = 4\pi \int \mathcal{P}_\mathcal{R}(k)|\Theta^{\rm our}_\ell(k)|^2\frac{dk}{k}$$ where the time is taken to be today (when is when we observe the CMB, so $\Theta^{\rm our}_\ell(k) = \Theta^{\rm our}_\ell(k,t=t_{\rm today})$). The difference wrt what is in Dodelson is that he expresses this in terms of $P(k)$ which gives us the equivalent formulation $$C_\ell = \frac{2}{\pi}\int k^2 P(k)|\Theta^{\rm our}_\ell(k)|^2 dk$$