Table of contents

Back to the main page

1 Background Cosmology 2 Thermodynamics and statistical mechanics 3 Cosmological perturbation theory 4 The full Einstein-Boltzmann system

This is a collection of equations and relations and quantities computed in this course. Below $x=\log a$ and a prime $f^\prime$ denotes differentiation wrt $x$. Likewise a dot $\dot{f}$ represents differentiation wrt $t$. In most of the equations we use units for which $\hbar = c = k_c = 1$, but they are restored in some key equations relating to the numerical project. We always work with scalar perturbations in the Newtonian gauge.

Background cosmology

Here is an overview of all the relevant equations for background cosmology.

Conformal time: $$a d\eta \equiv dt \implies \eta = \int_0^t \frac{dt}{a} = \int_0^a \frac{da}{a^2H} = \int_z^\infty \frac{dz}{H} = \int_{-\infty}^x \frac{dx}{aH}$$ Comoving distance: $$\chi = \eta_0 - \eta$$ where $\eta_0$ is the value today.

Angular diameter distance: $$d_A = aS_k(\chi)$$ where $$ S_k(\chi) = \begin{cases} \sin \left( \sqrt{|\Omega_{k 0}}| H_0 \chi /c \right)/\left(H_0\sqrt{|\Omega_{k 0}|} / c\right) & \Omega_{k 0} < 0\\ \chi & \Omega_{k 0} =0 \\ \sinh \left( \sqrt{|\Omega_{k 0}|} H_0 \chi / c\right)/\left(H_0\sqrt{|\Omega_{k 0}|} / c\right) & \Omega_{k 0} > 0 \end{cases} $$

Luminosity distance: $$d_L = \frac{d_A}{a^2}$$

Redshift: $$1+z = \frac{1}{a}$$ Cosmic time: $$t = \int_0^a \frac{da}{aH}$$ Very rough approximation for the temperature of the primordial plasma as function of time in the early Universe: $$\frac{T}{1 \rm MeV} \sim \left(\frac{1 {\rm sec}}{t}\right)^{1/2}$$ Photon temperature: $$T_\gamma = \frac{T_{\rm CMB 0}}{a}$$ where $T_{\rm CMB 0} = 2.7255$ K.

Neutrino temperatur (after decoupling): $$T_\nu = T_\gamma \left(\frac{4}{11}\right)^{1/3}$$ A more accurate calculation gives that its effectivelty a factor $1.0038$ higher than this (this is why $N_{\rm eff} = 3\cdot (1.0038)^4 = 3.046$ instead of just $3$).

The Hubble function: $$H = \frac{1}{a} \frac{da}{dt} = \frac{d\log a}{dt}$$ and the conformal Hubble function is: $$\mathcal{H} = \frac{1}{a}\frac{da}{d\eta} = \frac{da}{dt} = aH$$ The metric components: For a flat FRLW metric in cartesian coordinates $$ds^2 = -dt^2 + a^2(t)(dx^2 + dy^2 + dz^2)$$ $$g_{00} = -1$$ $$g_{ij} = a^2 \delta_{ij}$$ $$g^{00} = -1$$ $$g^{ij} = a^{-2} \delta^{ij}$$ The non-zero Christoffel symbols: $$\Gamma^0_{ij} = a\frac{da}{dt}\delta_{ij} = a^2 H \delta_{ij}$$ $$\Gamma^i_{j0} = \frac{1}{a}\frac{da}{dt}\delta^i_j = H\delta^i_j$$ The Ricci tensor components: $$R_{00} = -3\frac{\ddot{a}}{a}$$ $$R_{i0} = 0$$ $$R_{ij} = (a\ddot{a} + 2\dot{a}^2)\delta_{ij}$$ The Einstein tensor components: $$G_{00} = 3\frac{\dot{a}^2}{a^2}$$ $$G_{i0} = 0$$ $$G_{ij} = \left(-2\ddot{a}a - \dot{a}^2\right)\delta_{ij}$$ The energy momentum tensor for a perfect fluid: $$T^{\mu\nu} = U^\mu U^\nu (\rho + p) + p g^{\mu\nu} $$ where $U^\mu = (1,0,0,0)$ when co-moving with the expansion. The (upper) components are: $$T^{00} = \rho$$ $$T^{0j} = 0$$ $$T^{ij} = a^{-2}p \delta_{ij} $$ The (lower) components are: $$T_{00} = \rho$$ $$T_{0j} = 0$$ $$T_{ij} = a^{2}p \delta_{ij} $$ The (mixed) components are: $$T^0_0 = -\rho$$ $$T^0_j = T^j_0 = 0$$ $$T^i_j = p \delta^i_j $$ The equation of state: $$w \equiv \frac{p}{\rho}$$ For dust and cold dark matter $w=0$, for radiation $w = 1/3$, for curvature $w = -1/3$ and for dark energy $w=-1$.

The continuity equation: $$\nabla_\mu T^{\mu\nu} = 0 \implies \frac{d\rho}{dt} + 3H(\rho + p) = 0$$ For a constant equation of state $w$ it evoloved as $\rho \propto a^{-3(1+w)}$. The different energy components (baryons, cold dark matter photons, massless neutrinos, cosmological constant / dark energy) evolve as: $$\rho_b \propto a^{-3}$$ $$\rho_{\rm CDM} \propto a^{-3}$$ $$\rho_\gamma \propto a^{-4}$$ $$\rho_\nu \propto a^{-4}$$ $$\rho_\Lambda = {\rm constant}$$ The Friedmann equations: $$H^2 + \frac{k}{a^2} = \frac{8\pi G}{3} \rho$$ $$\frac{\ddot{a}}{a} = -\frac{4\pi G}{3}(\rho + 3p)$$ where $\rho$ ($p$) is the sum of the energy density (pressure) over all the different components. For a flat Universe $k=0$.

Critical density: $$\rho_c = \frac{3H^2}{8\pi G}$$ Density parameters: $$\Omega_i \equiv \frac{\rho_i}{\rho_c} = \frac{8\pi G \rho_i}{3H^2} = \frac{\Omega_{i 0}}{a^{3(1+w)} (H/H_0)^2}$$ and $$\Omega_k = -\frac{3k}{8\pi G \rho_c}$$ Total matter density is sum of all matter $\Omega_M = \Omega_b + \Omega_{\rm CDM}$ (plus additional massive neutrinos if that is included). Total radiation density is sum of all radiation $\Omega_R = \Omega_\gamma + \Omega_\nu$. Physical density parameters is the combination $\Omega_i h^2$. The Hubble function in terms of density parameters today: $$H = H_0\sqrt{\frac{\Omega_{k0}}{a^2} + \frac{\Omega_{b0}}{a^3} + \frac{\Omega_{\rm CDM 0}}{a^3} + \frac{\Omega_{\gamma 0}}{a^4} + \frac{\Omega_{\nu 0}}{a^4} + \Omega_{\Lambda 0}}$$ The sum of all density parameters is unity $\sum \Omega_i = 1$ (so $\Omega_k = 1 - \sum_{i\not= k} \Omega_i = 1 - \frac{\rho_{\rm total}}{\rho_c}$). The scale-factor evolves as $a \propto t^{1/2}$ in the radition era, $a \propto t^{2/3}$ in the matter era and $a\propto e^{Ht}$ in the dark energy dominated era. When dominated by a fluid of equation of state $w$ we have $a\propto t^{\frac{2}{3(1+w)}}$.

Matter-radiation equality: $$\Omega_R = \Omega_M \implies a = \frac{\Omega_{R0}}{\Omega_{M0}}$$ Matter-Dark energy equality: $$\Omega_M = \Omega_\Lambda \implies a = \left(\frac{\Omega_{M0}}{\Omega_{\Lambda 0}}\right)^{1/3}$$ Onset of accelerated expansion: $$\ddot{a} = 0 \implies a \simeq \left(\frac{\Omega_{M0}}{2\Omega_{\Lambda 0}}\right)^{1/3}$$

The 4-momentum in the background:
The 4-momentum $P^\mu = \frac{dx^\mu}{d\lambda}$ satisfy $P^\mu P_\mu = -m^2$ and $p^2 \equiv g_{ij}P^iP^j$. In the background $P^0 = E = \sqrt{p^2+m^2}$, $P^i = \frac{p}{a}\hat{p}^i$ where $\hat{p}^i\hat{p}^j\delta_{ij} = 1$ ($\hat{p}$ denotes the direction).

The geodesic equation: $$\frac{dP^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta}P^\alpha P^\beta = 0$$ For $\mu=0$ this reduces to $$p \propto \frac{1}{a}$$ For $\mu=i$ this reduces to $$\frac{d\hat{p}^i}{dt} = 0$$

Thermodynamics and statistical mechanics

Distribution function in thermal equilibrium: $$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 (Fermi-Dirac distribtion).

Fluid quantities from the distribution function: $$n = \frac{g}{(2\pi)^3}\int f d^3p$$ $$\rho = \frac{g}{(2\pi)^3}\int Ef d^3p$$ $$P = \frac{g}{(2\pi)^3}\int \frac{p^2}{3E}f d^3p$$ and in general $$ 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. For a flat FRLW metric we have $-\det g = a^6$ and with perturbations (Newtonian gauge metric with scalar perturbations) $$-\det g = a^6(1+2\Psi)(1+2\Phi)^3 \simeq a^6(1+2\Psi+6\Phi)$$ so $$ T^\mu_\nu = \frac{g}{(2\pi)^3} \int d^3p \frac{P^\mu P_\nu}{E} f $$

For relativistic ($T \gg m,\mu$) bosons: $$n = \frac{\zeta(3)}{\pi^2}gT^3$$ $$\rho = \frac{\pi^2}{30} gT^4$$ $$P = \frac{1}{3}\rho$$ For relativistic ($T \gg m,\mu$) fermions: $$n = \frac{3}{4}\frac{\zeta(3)}{\pi^2}gT^3$$ $$\rho = \frac{7}{8}\frac{\pi^2}{30} gT^4$$ $$P = \frac{1}{3}\rho$$ Non-relativistic limit ($T \ll m$) for massive particles of mass $m$ (bosons and fermions): $$n \approx g\left(\frac{mT}{2\pi}\right)^{3/2}e^{-(m-\mu)/T}$$ $$\rho \approx m n + \frac{3}{2}nT \approx m n$$ $$P \approx nT \implies w = \frac{P}{\rho} \ll 1$$

The collision term in the Boltzmann equation, $\frac{df_1}{d\lambda} = C(f_1,f_2,f_3,f_4)$, for a $1+2\leftrightharpoons 3+4$ process is: $$ \begin{align} C &= \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 $f_i = f_i(p_i)$. We can usually approximate $1 \pm f \approx 1$.

The master equation: (zeroth moment of the Boltzmann equation above) for a $1+2\leftrightharpoons 3+4$ process is $$\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)$$ where $\left<\sigma v\right>$ is the thermally averaged cross-section.

The Saha approximation:
$$n_1n_2 - n_3 n_4 \left(\frac{n_1n_2}{n_3n_4}\right)_{\rm eq} \approx 0$$

Saha equation for hydrogen recombination: $$\frac{X_e^2}{(1 - X_e)} = \frac{1}{n_b}\left(\frac{k_bT m_e}{2\pi c^2\hbar^2}\right)^{3/2} e^{-\frac{\epsilon_0}{k_bT}}$$ where $X_e = \frac{n_e}{n_b}$ is the free electron fraction, $n_b$ is the total baryon number density. This assumes no helium $Y_p = 0$.

Saha equation for helium and hydrogen recombination: $$x_{\rm He+} = \frac{n_{\rm He+}}{n_{\rm He}},\,\,\, x_{\rm He++} = \frac{n_{\rm He++}}{n_{\rm He}},\,\,\, x_{\rm H+} = \frac{n_{\rm H+}}{n_{\rm H}},\,\,\,$$ $$n_e\frac{x_{\rm He+}}{1-x_{\rm He+}-x_{\rm He++}} = 2\left(\frac{m_eT_b}{2\pi}\right)^{3/2}e^{-\chi_0/T_b}$$ $$n_e\frac{x_{\rm He++}}{x_{\rm He+}} = 4\left(\frac{m_eT_b}{2\pi}\right)^{3/2}e^{-\chi_1/T_b}$$ $$n_e\frac{x_{\rm H+}}{1-x_{\rm H+}} = \left(\frac{m_eT_b}{2\pi}\right)^{3/2}e^{-\epsilon_0/T_b}$$ where $m_e$ is the electron mass, $\chi_0 = 24.5874$eV is the ionization energy of neutral helium and $\chi_1 = 4\epsilon_0 = 54.42279$ eV the ionization energy of $He^+$. The free electron number density is $n_e = 2n_{\rm He++} + n_{\rm He+} + n_{\rm H+}$ or $$\frac{n_e}{n_b} = (2x_{\rm He++} + x_{\rm He+})\frac{Y_p}{4} + x_{\rm H+}(1-Y_p) \equiv f_e \implies X_e = \frac{f_e}{1-Y_p}$$

Peebles equation: $$\frac{dX_e}{dx} = \frac{C_r(T_b)}{H} \left[\beta(T_b)(1-X_e) - n_H \alpha^{(2)}(T_b)X_e^2\right],$$ where $T_b$ is the baryon temperature. In this course we simply use the approximation that this is always equal to the photon temperature $T_b \approx T_{\rm CMB 0} / a$ (see appendix of Milestone 2 for more info regarding this). After restoring the constants $c,k_b,\hbar$ this becomes $$ \begin{align} C_r(T_b) &= \frac{\Lambda_{2s\rightarrow1s} + \Lambda_{\alpha}}{\Lambda_{2s\rightarrow1s} + \Lambda_{\alpha} + \beta^{(2)}(T_b)}, \,\,{\rm (dimensionless)}\\ \Lambda_{2s\rightarrow1s} &= 8.227 \textrm{s}^{-1},\\ \Lambda_{\alpha} &= H\frac{(3\epsilon_0)^3}{(8\pi)^2 c^3\hbar^3 n_{1s}}, \,\,{\rm (dimension~1/s)}\\ n_{1s} &= (1-X_e)n_H, \,\,{\rm (dimension~1/m^3)}\\ n_H &= (1-Y_p)n_b, \,\,{\rm (dimension~1/m^3)}\\ n_b &= (1-Y_p)\frac{3H_0^2 \Omega_{b0}}{8\pi G m_H a^3}, \,\,{\rm (dimension~1/m^3)}\\ \beta^{(2)}(T_b) &= \beta(T_b) e^{\frac{3\epsilon_0}{4k_bT_b}}, \,\,{\rm (dimension~1/s)}\\ \beta(T_b) &= \alpha^{(2)}(T_b) \left(\frac{m_e k_b T_b}{2\pi \hbar^2}\right)^{3/2} e^{-\frac{\epsilon_0}{k_bT_b}}, \,\,{\rm (dimension~1/s)}\\ \alpha^{(2)}(T_b) &= \frac{8}{\sqrt{3\pi}} c\sigma_T\sqrt{\frac{\epsilon_0}{k_bT_b}}\phi_2(T_b), \,\,{\rm (dimension~m^3/s)}\\ \phi_2(T_b) &= 0.448\ln\left(\frac{\epsilon_0}{k_bT_b}\right), \,\,{\rm (dimensionless)} \end{align} $$

Reionization model: $$X_e = X_e^{\rm Peebles} + \frac{1+f_{\rm He}}{2}\left(1+\tanh\frac{y_{\rm reion}-y}{\Delta y_{\rm reion}}\right)$$ where $y = (1+z)^{3/2} = e^{-3x/2}$, $f_{\rm He} = \frac{Y_p}{4(1-Y_p)}$, $y_{\rm reion} = (1+z_{\rm reion})^{3/2}$, $\Delta y_{\rm reion} = \frac{3}{2}\sqrt{1+z_{\rm reion}}\Delta z_{\rm reion}$. In addition to this Helium probably gets doubly reionized at late reshifts which we can model by adding a term (again just using CAMB's model) $$\frac{f_{\rm He}}{2}\left(1 + \tanh \frac{z_{\rm He reion}-z}{\Delta z_{\rm He reion}}\right)$$ to the equation above. The best-fit values from Planck 2018: $z_{\rm reion} = 8$, $\Delta z_{\rm reion} = 0.5$, $z_{\rm He reion} = 3.5$ and $\Delta z_{\rm He reion} = 0.5$.

Cosmological perturbation theory

The 4-momentum in perturbation theory:
$P^\mu$ always satisfy $P^\mu P_\mu = -m^2$ and the magnitude of the 3-momentum is as usual defined as $p^2 \equiv g_{ij}P^iP^j$. This gives $P^0 = \frac{E}{\sqrt{-g_{00}}}$ where $E = \sqrt{p^2+m^2}$. For a diagonal metric $P^i = \frac{p}{\sqrt{g_{ii}}}\hat{p}^i$ where $ \hat{p}^i$ is a unit vector satisfying $\hat{p}^i \hat{p}^j \delta_{ij} \equiv 1$ and with $g_{00} = -(1+2\Psi)$, $g_{ii} = a^2(1+2\Phi)$ (Newtonian gauge with scalar perturbations) we have $$P^0 = E(1-\Psi)$$ $$P^i = \frac{p\hat{p}^i}{a}(1-\Phi)$$ $$P_0 = E(1+\Psi)$$ $$P_i = ap\hat{p}_i(1+\Phi)$$ Where indiced on the vectors here are lowered and raised with the Kroncker symbol so $\hat{p}^i = \hat{p}_i$.

Newtonian gauge metric (scalar perturbations): $$\mathrm{d} s^{2}= -dt^2(1+2\Psi) + a^2(1+2\Phi)(dx^2 + dy^2 + dz^2)$$ Metric components: $$g_{00} = -(1+2\Psi)$$ $$g_{i0} = 0$$ $$g_{ij} = a^2(1+2\Phi)\delta_{ij}$$ $$g^{00} = - \frac{1}{1+2\Psi} \simeq -1 + 2\Psi$$ $$g^{i0} = 0$$ $$g^{ij} = \frac{1}{a^2(1+2\Phi)} \delta^{ij} \simeq a^{-2}(1-2\Phi)\delta_{ij}$$ Christoffel symbols: $$\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}$$ 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}$$ 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} - 6H(\dot{\Psi} - 4\dot{\Phi}) - 4\frac{\nabla^2\Phi}{a^2}$$ Einstein tensor: $$G^0_0 = -3H^2 -6H\dot{\Phi} + 6\Psi H^2 + \frac{2\nabla^2\Phi}{a^2}$$ Defining the projection tensor $P^i_j = \frac{k^ik_j}{k^2} - \frac{1}{3}\delta^i_j$ then $$P^i_j G^j_i = -\frac{2}{3a^2}\nabla^2(\Phi + \Psi)$$ Einstein equations:
From $G^0_0 = 8\pi G T^0_0$ we get the Poisson equation $$-\frac{1}{a^2}\nabla^2\Phi + 3H\frac{\partial \Phi}{\partial t} - 3H^2\Psi = 4\pi G \sum \overline{\rho}_i \delta_i$$ where $\delta_\gamma = 4\Theta_0$ and $\delta_\nu = 4\mathcal{N}_0$. From $P^i_j G^j_i = 8\pi GP^i_j T^j_i$ we get an equation for the anisotropic stress $$\frac{1}{a^2}\nabla^2(\Phi + \Psi) = 32\pi G[ \overline{\rho}_\gamma \Theta_2 + \overline{\rho}_\nu \mathcal{N}_2]$$

Perturbations:
For photons and neutrinos we perturb the temperature in the distribution function $$T_\gamma \equiv \overline{T}_\gamma(1+ \Theta)$$ $$T_\nu \equiv \overline{T}_\nu(1+ \mathcal{N})$$ The density perturbations $\delta_i$ are defined via $$\rho_i \equiv \overline{\rho}_i(1+\delta_i)$$ For photons $\delta_\gamma = 4\Theta_0$ and for neutrinos $\delta_\nu = 4\mathcal{N}_0$ (the monopole of the distribution).

Gauge invariant density contrast: $$\Delta = \delta - 3(1+w)\frac{\mathcal{H}}{ck}v$$

Multipole expansion:
For the photon and neutrino temperature perturbations we can expand it in multipoles $$\Theta(t,k,\mu) = \sum \frac{2\ell + 1}{i^\ell}\Theta_\ell(t,k) P_\ell(\mu)$$ $$\Theta_\ell = \frac{i^\ell}{2}\int_{-1}^1 \Theta(t,k,\mu) P_\ell(\mu)d\mu$$ where $\mu = \frac{\vec{k}\cdot\vec{p}}{kp}$ is the angle between the Fourier mode and the propagation direction and $P_\ell$ is the Legendre polynomials.

Transforming linear equations to Fourier space:
$$\nabla \to i\vec{k}\,\,\,, \frac{\partial}{\partial x^j} \to ik_j \,\,\,, \nabla^2 \to -k^2$$

The energy momentum tensor from the distribution function:
$$T^\mu_\nu = \frac{g}{(2\pi)^3} \int \frac{dP_1dP_2dP_3}{\sqrt{-\det g}} \frac{P^\mu P_\nu}{P^0} f$$ which for scalar perturbations in the Newtonian gauge becomes $$T^\mu_\nu = \frac{g}{(2\pi)^3} \int d^3p\, \frac{P^\mu P_\nu}{E} f$$ and gives us $$T^0_0 = -\frac{g}{(2\pi)^3} \int d^3p\, E f \equiv -\rho$$ $$T^i_j = \frac{g}{(2\pi)^3} \int d^3p\, \frac{p^2\hat{p}^i\hat{p}_j}{E} f$$ $$T^i_0 = -\frac{g}{(2\pi)^3} \int d^3p\, \frac{p\hat{p}^i}{a}$$ $$T^0_i = \frac{g}{(2\pi)^3} \int d^3p\, ap\hat{p}^i f$$

The energy momentum tensor for radiation:
$$T^0_0 = -\overline{\rho}(1+4\Theta_0)$$ $$P^j_i T^i_j = \left(\frac{k_i k^j}{k^2} - \frac{1}{3}\delta_i^j\right)T^i_j = -\frac{8}{3}\overline{\rho}\Theta_2$$ $$\frac{ik_i}{k}T^i_0 = -(\overline{\rho} + \overline{P})\frac{3\Theta_1}{a}$$ and the same holds for neutrinos with $\Theta\to \mathcal{N}$.

The energy momentum tensor for non-relativistic particles:
$$T^0_0 = -\overline{\rho}(1+\delta)$$ $$T^i_0 = -\overline{\rho} \frac{\vec{v}^i}{a}$$ $$T^0_i = \overline{\rho} a\vec{v}_i$$ $$T^i_j = 0$$ In Fourier space the scalar velocity $v$ is defined via $v^i = \frac{ik^i}{k} v$ (scalar perturbations are curl free so $v^i \propto k^i$) which gives $$\frac{ik_i}{k} T^i_0 = \overline{\rho} \frac{v}{a}$$

The energy momentum tensor in general:
The perturbation of the energy momentum tensor for a perfect fluid is: $$\delta T^\mu_\nu = (\delta\rho + \delta P)\overline{U}^\mu \overline{U}_\nu + (\overline{\rho}+\overline{P})(\delta U^\mu \overline{U}_\nu + \overline{U}^\mu \delta U_\nu) + \delta P\delta^\mu_\nu + \Pi^\mu_\nu$$ where $g_{\mu\nu}U^\mu U^\nu = -1$, $\overline{U}^\mu = \delta^\mu_0$, $\overline{U}_\mu = -\delta_\mu^0$, $\delta U^0 = -\Psi$, $\delta U^i \equiv \frac{v^i}{a}$ where $v^i = \frac{dx^i}{d\tau}$ is the coordinate velocity so $U^\mu = (1-\Psi,v^i/a)$. Thus $$\delta T^0_0 = -\delta\rho$$ $$\delta T^i_0 = (\overline{\rho}+\overline{P})v^i$$ $$\delta T^i_j = \delta P \delta^i_j + \Pi^i_j$$ Here $\Pi^\mu_\nu$ is the anisotropic stess satisfying $\Pi^0_0 = \Pi^i_0 = 0$ with $\Pi^i_i = 0$. For scalar perturbations this leads to a scalar potential $\sigma$ defined below.

We get $$T^0_0 = -\overline{\rho}(1+\delta)$$ $$\frac{ik_i}{k} T^i_0 = (\overline{\rho}+\overline{P}) \frac{v}{a}$$ $$P_i^jT^i_j = \left(\frac{k_i k^j}{k^2} - \frac{1}{3}\delta_i^j\right)T^i_j = -(\overline{\rho}+\overline{P})\sigma$$ where $\delta$ is the overdensity, $v$ is the scalar velocity and $\sigma$ the scalar anisotropic stress (i.e. this is how these quantities are defined). Comparing to this gives us that the photon velocity is therefore $v_\gamma = -3\Theta_1$, the photon overdensity is $\delta_\gamma = 4\Theta_0$ and the anisotropic stress $\sigma_\gamma = 2\Theta_2$. In some texts the velocity divergence $\theta = \nabla\cdot \vec{v}$ is used as a variable. This is related to our $v$ via $\theta = -kv$. In real space the second equation above is $$\frac{\partial}{\partial x^i} T^i_0 = (\overline{\rho}+\overline{P}) \frac{\nabla\cdot \vec{v}}{a}$$ where $\nabla\cdot \vec{v} = \frac{\partial}{\partial x^i} v^i = \theta$ is the velocity divergence (the curl-free part of the velocity).

The Einstein equations in general:
$$-\frac{k^2}{a^2}\Phi + 3H(-\dot{\Phi} + H\Psi) = -4\pi G \sum_i\overline{\rho}_i\delta_i$$ $$\frac{k^2}{a^2}(-\dot{\Phi} + H\Psi) = -4\pi G \sum_i (\overline{\rho}_i + \overline{P}_i) kv_i$$ $$\frac{k^2}{a^2}(\Phi + \Psi) = -12\pi G \sum_i (\overline{\rho}_i + \overline{P}_i)\sigma_i$$ where only two are independent (one follows from the others and the continuity equation). The first and the second above gives $$\frac{k^2}{a^2}\Phi = 4\pi G\sum_i \overline{\rho}_i\Delta_i$$ where $\Delta$ is the gauge invariant density perturbations (thus expressed in terms of $\Delta$ then $\Phi$ is not dynamical). The total matter density contrast is $$\overline{\rho}_M\Delta_M \equiv \sum_i \overline{\rho}_i\Delta_i \to \Delta_M = \frac{k^2\Phi}{\frac{3}{2} a^{-1}\Omega_M H_0^2}$$

The conservation equations in general:
$$\nabla^\mu T_{\mu\nu} = 0$$ gives $$\delta^\prime = (1+w)\left(\frac{ck}{\mathcal{H}}v - 3\dot{\Phi}\right) - 3\left(c_s^2 - w\right)\delta$$ $$v^\prime = -(1-3w)v - \frac{w^\prime}{1+w}v - \frac{c_s^2}{1+w}\frac{ck}{\mathcal{H}}\delta + \frac{ck}{\mathcal{H}}\sigma - \frac{ck}{\mathcal{H}}\Psi = 0$$ where $c_s^2 = \frac{\delta P}{\delta \rho}$ which for adiabatic fluctuations is just $\frac{\overline{P}^\prime}{\overline{\rho}^\prime} = w - \frac{w^\prime}{3(1+w)} = w$ since $w^\prime = 0$ for the cases of interest. For baryons and CDM its a good approximation to use $c_s^2 = 0$.

The full Einstein-Boltzmann system

Photon temperature multipoles: $$ \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 \lt \ell_{\textrm{max}} \\ \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}}\\ \end{align} $$ Photon polarization multipoles: $$ \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 \lt \ell_{\textrm{max}} \\ \Theta_{P,\ell}^\prime &= \frac{ck}{\mathcal{H}} \Theta_{\ell-1}^P-c\frac{\ell+1}{\mathcal{H}\eta(x)}\Theta_\ell^P+\tau^\prime\Theta_\ell^P, \quad\quad \ell = \ell_{\textrm{max}}\\ \end{align} $$ Neutrino multipoles: $$ \begin{align} \mathcal{N}^\prime_0 &= -\frac{ck}{\mathcal{H}} \mathcal{N}_1 - \Phi^\prime, \\ \mathcal{N}^\prime_1 &= \frac{ck}{3\mathcal{H}} \mathcal{N}_0 - \frac{2ck}{3\mathcal{H}}\mathcal{N}_2 + \frac{ck}{3\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 \lt \ell_{\textrm{max},\nu} \\ \mathcal{N}_{\ell}^\prime &= \frac{ck}{\mathcal{H}} \mathcal{N}_{\ell-1}-c\frac{\ell+1}{\mathcal{H}\eta(x)}\mathcal{N}_{\ell}, \quad\quad \ell = \ell_{\textrm{max},\nu}\\ \end{align} $$ Cold dark matter and baryons: $$ \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 \\ \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(3\Theta_1 + v_b) \\ \end{align} $$ Metric potentials: $$ \begin{align} \Phi^\prime &= \Psi - \frac{c^2k^2}{3\mathcal{H}^2} \Phi + \frac{H_0^2}{2\mathcal{H}^2} \left[\Omega_{\rm CDM 0} a^{-1} \delta_{\rm CDM} + \Omega_{b 0} a^{-1} \delta_b + 4\Omega_{\gamma 0} a^{-2}\Theta_0 + 4\Omega_{\nu 0}a^{-2}\mathcal{N}_0\right] \\ \Psi &= -\Phi - \frac{12H_0^2}{c^2k^2a^2}\left[\Omega_{\gamma 0}\Theta_2 + \Omega_{\nu 0}\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}$.

Tight coupling regime:
In the tight coupling regime, the only differences are 1) that one should only include $\ell = 0$ and 1 for $\Theta_\ell$ (and none for polarization) - all higher moments are given by those, and 2) that the expressions for $\Theta'_1$ and $v_b'$ are quite a bit more involved (see the appendix for how to derive this): $$ \begin{align} q &= \frac{-[(1-R)\tau^\prime + (1+R)\tau^{\prime\prime}](3\Theta_1+v_b) - \frac{ck}{\mathcal{H}}\Psi + (1-\frac{\mathcal{H}^\prime}{\mathcal{H}})\frac{ck}{\mathcal{H}}(-\Theta_0 + 2\Theta_2) - \frac{ck}{\mathcal{H}}\Theta_0^\prime}{(1+R)\tau^\prime + \frac{\mathcal{H}^\prime}{\mathcal{H}} - 1}\\ v_b^\prime &= \frac{1}{1+R} \left[-v_b - \frac{ck}{\mathcal{H}}\Psi + R(q + \frac{ck}{\mathcal{H}}(-\Theta_0 + 2\Theta_2) - \frac{ck}{\mathcal{H}}\Psi)\right]\\ \Theta^\prime_1 &= \frac{1}{3} (q - v_b^\prime) \end{align} $$ In the tight coupling regime, we get the same expressions for the higher-ordered photon moments as given by the initial conditions, \begin{align} \Theta_2 &= \left\{ \begin{array}{l} -\frac{8ck}{15\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(with polarization)} \\ -\frac{20ck}{45\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(without polarization)} \end{array}\right. \\ \Theta_\ell &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau'} \Theta_{\ell-1}, \quad\quad \ell \gt 2\\ \Theta_0^P &= \frac{5}{4}\Theta_2 \\ \Theta_1^P &= -\frac{ck}{4\mathcal{H}\tau'}\Theta_2 \\ \Theta_2^P &= \frac{1}{4}\Theta_2 \\ \Theta_\ell^P &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau'} \Theta_{\ell-1}^P, \quad\quad \ell \gt 2 \end{align}

Initial conditions:
$$ \begin{align} \Psi &= \frac{1}{\frac{3}{2} + \frac{2f_\nu}{5}}\\ \Phi &= -(1+\frac{2f_\nu}{5})\Psi \\ \delta_{\rm CDM} &= \delta_b = -\frac{3}{2} \Psi \\ v_{\rm CDM} &= v_b = -\frac{ck}{2\mathcal{H}} \Psi\\ &\text{Photons:}\\ \Theta_0 &= -\frac{1}{2} \Psi \\ \Theta_1 &= +\frac{ck}{6\mathcal{H}}\Psi \\ \Theta_2 &= \left\{ \begin{array}{l} -\frac{8ck}{15\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(with polarization)} \\ -\frac{20ck}{45\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(without polarization)} \end{array}\right. \\ \Theta_\ell &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau^\prime} \Theta_{\ell-1}\\ &\text{Photon Polarization:}\\ \Theta_0^P &= \frac{5}{4} \Theta_2 \\ \Theta_1^P &= -\frac{ck}{4\mathcal{H}\tau'} \Theta_2 \\ \Theta_2^P &= \frac{1}{4}\Theta_2 \\ \Theta_\ell^P &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau^\prime} \Theta_{\ell-1}^P \\ &\text{Neutrinos:}\\ \mathcal{N}_0 &= -\frac{1}{2} \Psi \\ \mathcal{N}_1 &= +\frac{ck}{6\mathcal{H}}\Psi \\ \mathcal{N}_2 &= +\frac{c^2k^2 a^2 \Psi}{30H_0^2(\Omega_{\nu 0}+\Omega_{\gamma 0})}\\ \mathcal{N}_\ell &= \frac{ck}{(2\ell+1)\mathcal{H}} \mathcal{N}_{\ell-1}, \quad\quad \ell \ge 3 \end{align} $$ where $f_{\nu} = \frac{\Omega_{\nu 0}}{\Omega_{\gamma 0} + \Omega_{\nu 0}}$. The normalization is such that the gauge invariant curvature perturbations $\mathcal{R} \equiv 1$ where the (dimensionless) power-spectrum of $\mathcal{R}$ is the usual $\mathcal{P}_\mathcal{R}(k) = A_s(k/k_{\rm pivot})^{n_s-1}$.

Line of sight integral:
The multipoles today are computed from the LOS integral $$\Theta_\ell(k) = \int_{-\infty}^{0} \tilde{S}(k,x) j_\ell[k(\eta_0-\eta)] dx,$$ where the source function is defined as $$\tilde{S}(k,x) = \tilde{g}\left[ \Theta_0 + \Psi + \frac{1}{4}\Pi\right] + e^{-\tau} \left[\Psi^\prime-\Phi^\prime\right] - \frac{1}{ck}\frac{d}{dx}(\mathcal{H}\tilde{g}v_b) + \frac{3}{4c^2k^2} \frac{d}{dx} \left[\mathcal{H}\frac{d}{dx} (\mathcal{H}\tilde{g}\Pi)\right].$$ where $\tilde{g} = -\tau^\prime e^{-\tau}$. Expanding in spherical harmonics:
$$\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})$$ $$a_{\ell m}(\vec{x},t) = \int d\Omega_{\hat{p}} Y_{\ell m}^*(\hat{p})\Theta(\vec{x},\hat{p},t)$$ The CMB power-spectrum: $$C_\ell \equiv \left\lt |a_{\ell m}|^2 \right\gt = \left\lt a_{\ell m}a_{\ell m}^* \right\gt.$$ which gives us $$C_\ell = 4\pi\int \mathcal{P}_{\mathcal{R}}(k) \Theta_\ell^2(k)\frac{dk}{k}$$ where $$\mathcal{P}_{\mathcal{R}}(k) = A_s \left(\frac{k}{k_{\rm pivot}}\right)^{n_s-1},$$ This assumes $\Theta_\ell$ is normalized such that the $\mathcal{R} \equiv 1$ initally.

Matter power-spectrum:
$$P(k,x) = |\Delta_M(k,x)|^2 \frac{2\pi^2}{k^3} \mathcal{P}_{\mathcal{R}}(k)$$ where $\Delta_M(k,x) \equiv \frac{c^2k^2\Phi(k,x)}{\frac{3}{2}\Omega_{M 0} a^{-1} H_0^2}$ and $\Omega_{M 0}$ is the total matter density parameter today. Power-spectrum of individual omponents are defined the same way with $\Delta_i = \delta_i - 3(1+w_i)\frac{\mathcal{H}}{ck}v_i$ being the gauge invariant density contrast.