Table of contents


Back to the main page

1 Background Cosmology 2 Thermodynamics and statistical mechanics 3 Cosmological perturbation theory 4 The 4-Momentum and the geodesic equation 5 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.


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 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}$$


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)$$

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


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} - 6(\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:
$$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{\dot{w}}{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 us $=0$ for baryons and CDM and $=1/3$ for radiation).


The 4-Momentum and the geodesic equation


The 4-momentum:
$P^\mu$ satisfy $P^\mu P_\mu = -m^2$. The magnitude of the 3-momentum is defined as $p^2 \equiv g_{ij}P^iP^j$. $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 $g_{00} = -1$, $g_{ii} = a^2$ at the background level and $g_{00} = -(1+2\Psi)$, $g_{ii} = a^2(1+2\Phi)$ in perturbation theory (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$.

The geodesic equation: $$\frac{dP^0}{d\lambda} + \Gamma^0_{\alpha\beta}P^\alpha P^\beta = 0 \implies p \propto \frac{1}{a}$$ $$\frac{dP^i}{d\lambda} + \Gamma^i_{\alpha\beta}P^\alpha P^\beta = 0 \implies \frac{d\hat{p}^i}{dt} = 0$$ The implications are for the background.


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].$$ 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.