Table of contents

Back to the main page

1. Thermodynamics and statistical mechanics

2. Thermal history of the Universe

Thermodynamics and statistical mechanics

Here is a brief summary of the thermodynamics / statistical mechanics 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). If you want a more basic introduction to statistical mechanics see for example Statistical mechanics lectures by Leonard Susskind on YouTube.

The Universe contains a hell of a lot of particles and photons that interact with each other and do all kinds of crazy stuff. Such interactions (like a proton catching an electron and becoming a hydrogen atom) are very important to take into account be able to understand how our Universe evolves. Obviously we cannot follow each and every one of the particles so to be able to describe the Universe we need a statistical framework that allows us to describe a bigger collection of particles that is also able to include the effect of interactions. This is where statistical mechanics and thermodynamics comes into play in cosmology and in particular the Boltzmann equation which describes the statistical behaviour of a thermodynamic system out of thermal equilibrium. The key quantity in the Boltzmann formalism is the distribution function.

Distribution function

The distribution function $f(t, \vec{x}, \vec{p})$ tells us how many particles (of a particular type - we have one distribution function for each species) we have in a small phase space volume (Video: Introduction to the distribution function in kinetic theory)) $$\Delta N = \frac{g}{(2\pi)^3} f(t, \vec{x}, \vec{p}) (\Delta \vec{x})^3 (\Delta \vec{p})^3$$ where $g$ is the number of internal degrees of freedom (spin/polarization). Why is this the "right" quantity to work with? We can obtain the usual macroscopic quantities (number-density, energy density, pressure) as weigthed integrals (expectation values) over 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 the average of some quantity $Q$ is given as $$\left\lt Q \right\gt = \frac{1}{n}\frac{g}{(2\pi)^3}\int Q f d^3 p = \frac{\int Q f d^3 p}{\int f d^3 p}$$ In thermal equilibrium it turns out the distribution function has a simple form that allows us to compute what $n,P,T$ are as functions of temperature and going beyond thermal equilibrium we can find how the distribution function evolves by solving the Boltzmann equation $$ \frac{df_{\rm particle~A}}{dt} = C_{\rm interactions}(f_{\rm particle~A},f_{\rm particle~B},\ldots) $$ where the right hand side contains information about interactions (inputs from quantum field theory). This is what makes it perfect for the job!

As mentioned above, in thermal equilibrium the distribution function is known and takes the 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 like photons (Bose-Einstein distribution) and $+$ is for fermions like electrons (Fermi-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 a few key facts:

  • The chemical potential of a species is energy that can be absorbed or released due to a change of the particle number of the given species, e.g. in a chemical reaction or phase transition.
  • 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$
  • $\mu = 0$ for photons. This is a consequence of the definition and the fact that photon numbers are not conserved
  • The chemical potential of an antiparticle is minus that of a particle. This is a conseqeunce of the former two points: the interaction $a + \overline{a} \leftrightarrow \gamma + \gamma$ gives $\mu + \overline{\mu} = 0 + 0$ so $\overline{\mu} = -\mu$.
If you want a more detailed overview of this quantity see e.g. Video: Lecture Chemical Potential

In the low temperature limit both these distributions approaches the classical Maxwell-Botzmann distribition $$f \approx e^{-\frac{E-\mu}{T}}$$ For our study of recombination in this course this approximation is excellent.

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 (stochastic) microscopic interactions that takes place.

The general relativistic Boltzmann equation reads (note that Dodelson simplifies the discussion of this by phrasing it directly 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 and $C$ is the collision term which contains information about interactions.

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}$ as our variables. Both ways are equivalent and the reason for this split is convenience. Note that for photons we simply have $E = p$.

The first thing we need is to express $P^0$ and $P^i$ in terms of $E,p,\hat{p}^i$ and the metric components. For this we can use the relations $p^2 \equiv g_{ij}P^i P^j$ and $g_{\mu\nu}P^\mu P^\nu = -m^2$. For example for a diagonal metric we find $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]$$ Why did we do this? It's because when we know the distribution function we can easily compute the partial derivatives and the prefactors are just terms that tell us how the particles move in the Universe, i.e. these terms ($\frac{dE}{d\lambda}$ and $\frac{d\hat{p}^i}{d\lambda}$) are determined by the geodesic equation. The advantage of the choosen set of coordinates will be seen once we look at this equation to first order in perturbation theory we will see that 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) so this makes the math a bit easier.

So how do we compute $\frac{dE}{d\lambda}$? We take the geodesic equation for $\mu = 0$ $$\frac{dP^0}{d\lambda} + \Gamma^0_{\alpha\beta} P^\alpha P^\beta = 0$$ which used with the expression for $P^0$ in terms of $E$ gives us an equation for $\frac{dE}{d\lambda}$ and 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}$ (which turns out we won't need in this course, but its useful to know how to derive this). One thing to note is that the derivative in the first term is a total derivative so we need 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$ when evaluation this in calculations (this 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 particle type $1$ this then 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 $1\pm f$ 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 considering), delta functions enforcing energy-momentum conservation, the Feynman amplitude $|\mathcal{M}|^2$ describing the probabillity/efficiency 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 form is quite intuitive and its not hard to see how this should be generalized to other types of processes.


Figure: we won't have time to cover the calculations of the different collision terms that we will encounter in this course and will just take the results as input and focus on the physical understanding of it, but see Dodelson for how to deal with them - its great practice to have gone through atleast one of the calculations.

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 weighted integrals of 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$$ Now just a the Boltzmann equation tells us how the distribution function evolves in time, integrals of the Boltzmann equation, i.e. integrals on the form $\int Q\frac{df}{d\lambda} d^3p = \int Q C d^3p $ some some function $Q$ of momenta, will tell us how the number density, velocity, velocity dispersion and so on evolve. This is as you can imagine super useful and its also helps with understanding the physics better as the distribution function can be a bit abstract while conservation equations (normal fluid equations) are something we have a lot more experience with.

Multiplying the Boltzmann equation by certain combinations of momenta and then integrating it over all momenta is called taking moments of the Boltzmann equation (this is similar to when we take moments of a random distribution in statistics). This will give us useful equations that generally describes conservation of certain quantities like mass, momentum, energy and so on.

To see how this works lets simplify and forget about general relativity for a moment and focus on a fluid of non-relativistic particles moving under Newtons laws. We have a set of particles with mass $m$ that are affected by some forces $\vec{F}$. The equation of motion is good old $m\vec{a} = \vec{F}$ or in terms of the momentum $\frac{d\vec{p}}{dt} = \vec{F}$. Since velocity is just momentum over mass the fluid velocity is therefore the average of $\vec{v} = \frac{\vec{p}}{m}$ over the distribution function: $$\vec{V} = \frac{1}{n}\frac{g}{(2\pi)^3}\int \frac{\vec{p}}{m} f d^3 p$$ Having introduced this lets look at the Boltzmann equation. The Boltzmann equation (ignoring collisions) for the distribution function of our particles $f = f(t,\vec{x},\vec{p})$ can be written $$\frac{df}{dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial \vec{x}} \cdot \frac{d\vec{x}}{dt}+ \frac{\partial f}{\partial \vec{p}} \cdot \frac{d\vec{p}}{dt} = 0$$ We can simplify this a bit by using the equation of motion for the particles making up the fluid: $\frac{d\vec{x}}{dt} = \frac{\vec{p}}{m}$ and $\frac{d\vec{p}}{dt} = \vec{F}$. We will also simplify the notation a bit by introducing the gradient $\nabla_x = \frac{d}{d\vec{x}}$ and similar for $\vec{p}$ which gives us $$\frac{\partial f}{\partial t} + \nabla_x f \cdot \frac{\vec{p}}{m}+ \nabla_p f \cdot \vec{F} = 0$$ Now we multiply this equation by some quantity $Q = Q(\vec{p})$ depending on momenta and integrate it over all momenta. We will start by doing it for a general $Q$ and then specify to some concrete examples in the end. When we integrate the first term in the Boltzmann equation is simply $$\frac{g}{(2\pi)^3}\int \frac{\partial f}{\partial t} Q d^3p = \frac{\partial}{\partial t}(n\left\lt Q\right\gt)$$ since we can pull the derivative outside the integral as $Q$ has no explicit time-dependence + we used the definition of the average. For the second term the same holds true for $x$ so $$\frac{g}{(2\pi)^3}\int Q \nabla_x f \cdot \frac{\vec{p}}{m} d^3p = \nabla_x\cdot (n\left\lt \frac{Q\vec{p}}{m}\right\gt)$$ For the third term we note that the force only depends on position so can be pulled out of the integral and it becomes $$\frac{g}{(2\pi)^3} \int Q\nabla_p f \cdot \vec{F} d^3p = -\frac{g}{(2\pi)^3} \int f \nabla_p Q d^3p \cdot \vec{F} = -(n\left\lt \nabla_p Q \right\gt) \cdot \vec{F}$$ where we performed an integration by parts to move the $\nabla_p$-derivative over to the $Q$-term. Putting it all together we get $$\frac{\partial}{\partial t}(n\left\lt Q\right\gt) + \nabla_x\cdot (n\left\lt \frac{Q\vec{p}}{m}\right\gt) - (n\left\lt \nabla_p Q \right\gt) \cdot \vec{F} = 0$$ Now we can start to look at some simple cases and extract some useful equations. If we take $Q=1$ then we get $$\frac{\partial n}{\partial t} + \nabla_x\cdot (n\vec{V}) = 0$$ Finally something we recognize! This is the well-known continuity equation that describes conservation of particle number (conservation of mass). If we take $Q = \vec{p}_i$ (the $i$'th component of the momentum) then the equation becomes $$\frac{\partial}{\partial t}(\rho \vec{V}_i) + \frac{\partial}{\partial x_j}(\rho \vec{V}_i\vec{V}_j) = - \frac{\partial}{\partial x_j}P_{ij} + n\vec{F}$$ where $\rho = mn$ and $P_{ij} \equiv \rho\left\lt (\vec{V}-\vec{v})_i(\vec{V}-\vec{v})_j \right\gt$ is the so-called pressure tensor. To make it a bit more familiar this can further be massaged into (after using the continuity equation to simplify) $$\frac{\partial \vec{V}}{\partial t} + (\vec{V} \cdot \nabla)\vec{V} = -\frac{\nabla_j P_{ij}}{\rho} + \frac{\vec{F}}{m}$$ This is the well-known Euler equation and describes conservation of momentum! We could go further and take $Q = \frac{p_ip_j}{2m}$ and get an equation for conservation of energy and there is no stop to the number of conservation equations we can extract. These equations have one annoying property: the equations for the $n$'th moment (taking $Q$ to be $n$ products of momenta $p_ip_j\cdots$) generally depends on the $n+1$'th moment! The mass equation tells us how the number density changes if we know the fluid velocity. The momentum equation tells us how the fluid velocity evolves however this requires knowledge of the velocity dispersion (in the pressure tensor) and so on and so on. The system usually never closes. However a good things is that in most systems of interest these higher order moments usually becomes smaller and smaller and we can often truncate this infinite hierachy and get a closed system. We will encounter this later in the course.

The two equations we mention above, the continuity equation and the Euler equation, will play a big role in this course so remember them and remember what they represent physially as being conservation laws for mass and momentum. A lot of what we will do with the Boltzmann equation in this course will be the same as is done above just for the case of particles and photons in an expanding Universe described by General Relativity. This compliates the algebra a bit, but the physics is largely the same (though we will also get some GR specific effects) so keep these examples in mind.

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 we have $$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$$ Evaluating the right hand side 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 this reduces to $$\frac{1}{T}\frac{dT}{dt} = -\frac{1}{a}\frac{da}{dt}$$ The solution is a fact we aready know: $T \propto 1/a$.

What about the integrated Boltzmann equation (taking the 0'th moment)? In the absence of interactions this simply becomes (see exercise) $$\frac{1}{a^3}\frac{d(n a^3)}{dt} = 0 \implies n \propto \frac{1}{a^3}$$ and as we would have expected from the example in the previous section this describes conservation of mass (number of particles), but here in an expanding Universe. In the next section we'll look at what this becomes when we finally include some interactions. What about higher moments? The first order moment is to multiply by the momentum and integrate, but for a smooth Universe this gives us nothing interesting (there are no peculiar velocities by definition so we just gets $0=0$).

See Baumann for a more complete discussion.

Boltzmann equation for the number density in a smooth Universe

The Boltzmann equation can describe any type of interaction. However in cosmology the main type of interaction that we will encounter is a so-called $1+2\leftrightarrow 3+4$ process. Some simple examples are electron-position annhialation / pair-production $e^- + e^+ \leftrightarrow \gamma + \gamma$ or hydrogen recombination / ionization $e^- + p^+ \leftrightarrow H + \gamma$. For such a 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>$ is the thermally averaged cross-section coming from the collison term. A full derivation from the collision term gives you what $\beta$ is, but we can see what it has to be in terms of $\alpha$ 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)$$ This equation is very important and will play a big role in this course. Some examples of interactions of this type that we will enounter are $$e^- + p^+ \rightleftharpoons e^- + p^+$$ $$e^- + \gamma \rightleftharpoons e^- + \gamma$$ $$e^- + p^+ \rightleftharpoons H + \gamma$$ The first one, Coulomb scattering, is key to keeping electrons and protons in equilibrium with eachother (this is this interactions main "job" in cosmology). The second one, Thompson scattering, is key for understanding the CMB as its the main form of scattering of CMB photons (we also have the same interaction with protons, but the interactions efficicency is $\propto 1/m^2$ and electrons are much less massive than protons) and the last one described how electrons and protons can go together and from a hydrogen atom (and the reverse: how light hitting a hydrogen atom can ionize it) and will be crucial for understanding how the Universe became transparent via a process we call recombination. There are many more, but most these ones are for sure the most important ones for our purposes. For each of these interactions there is a corresponding cross-section $\left<\sigma v\right>$ given by QFT that we (luckily; though its a very nice calculation that you learn how to do if you take a QFT course) will just take the result from and use for our purposes of understanding the CMB.

Decoupling, freeze-out and the Saha approximation

The Boltzmann equation we introduced in the previous section $$\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)$$ can be rewritten on a more suggestive form $$\frac{1}{n_1a^3}\frac{d(n_1 a^3)}{d x} = - \frac{\Gamma_1}{H}\left(1 - \frac{n_3 n_4}{n_1n_2} \left(\frac{n_1n_2}{n_3n_4}\right)_{\rm eq}\right)$$ where $\Gamma_1 \equiv n_2\left<\sigma v\right>$ is the interaction rate, $H$ the expansion rate and the second factor just tells us how close we are to equilibrium. Lets try to understand the solution of this general equation a bit better which will reveal a few concepts like decouping and freeze-out that will be very useful later on. If $\Gamma_1 \gg H$ then the prefactor on the right hand side is huge and the system is efficiently driven towards chemical equlibrium. This is not hard to see from the form of the equation and is also easy to understand physically: if the interaction rate is large relative to the expansion rate of the Universe then we have many interactions in the time it takes the Universe to double in size and we are easily able to maintain equilibrium $n_i \approx n_i^{\rm eq}$ even though the number density of each species gets smaller and smaller. However this cannot last forever. Remember that the equilibrium distribution for massive partices $n_i^{\rm eq} \propto e^{-\frac{m}{T}}$ gets exponensialy suppressed when the temperature goes down (so if equilibrium had lasted we would have no more massive particles left in the Universe today). Thus the equilibrium we have must eventually fail. Once $\Gamma_1$ drops below $H$ then the interactions are no longer efficient enough to maintain equilibrium and we say that it decouples. After it decouples the interaction is no longer relevant so knowing when ($\Gamma_1 \sim H$) this happens already gives us important information. If a particle species has no more efficient interactions that ties it to the rest of the primordial plasma then it can no longer stay in equilibrium with it and will evolve on its own like all the other stuff wasn't there (and if anything happens in the plasma that changes, say, its temperature then that change will not be reflected in stuff that has decoupled - see the problem of computing the neutrino temperature today for an example). As we reach $\Gamma_1 \ll H$ the right hand side is practically zero and in this regime the equation reduces to $\frac{d(n_1 a^3)}{d x} \approx 0 \implies n_1 \propto 1/a^3$ which is just normal volume dilution in an expanding Universe (i.e. the number of particles in a co-moving volume stays constant, its the physical volume that gets bigger $\propto a^3$) which is what we would expect if there was no interactions at all. When this happens for a massive particle species we say that it freezes out. The evolution is exponential decay of $n_1 a^3$ until decoupling, $\Gamma_1 \sim H$, and then it flattens out us and leaves us with a freeze-out abundance $n_1a^3 \approx $ constant. What this freeze-out abundance really is typucally requires numerical solutions of the Boltzmann equation (though we shall see an example of dark matter production where we can do some analytical estimates). We will encounter freeze-out in two main situations 1) when we look at recombination we will see that not all of the free electrons and protons will end up in atoms 2) freeze-out is one of the most common production mechanisms for producing dark matter (for example WIMP's).

The discussion above shows how we can qualitatively understand the solution to the Boltzmann equation, but we can also get a bit more quantitative understanding by extracting a very useful equation. Recall that if the interactions are efficient then we will be close to equilibrium and $\left(1 - \frac{n_3 n_4}{n_1n_2} \left(\frac{n_1n_2}{n_3n_4}\right)_{\rm eq}\right) \approx 0$. This gives us to the Saha approximation $$\frac{n_1n_2}{n_3n_4} \approx \left(\frac{n_1n_2}{n_3n_4}\right)_{\rm eq}$$ and 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, but this particular combination must stay roughly constant. One of the most important interaction we will encounter in this course is $e^- + p^+ \leftrightarrow H + \gamma$ and we will use this equation to tell us how the number density of free electrons will evolve in time, i.e. how fast the Universe goes from a plasma to a transparent Universe with atoms. The Saha approximation is great in that it's simple, doesn't require detailed knowledge of the interactions (QFT calculations) and allows us to get a better understanding for how the number density evolves in the presence of interactions. However it can only do so much, it will eventually break down and when the Saha approximation is no longer valid then we have no other option than solving the full Boltzmann equation above to be able to accurately pin down the evolution (and for our purposes of using this to study recombination to understand the CMB, for which observables are now known to subpercent accuracy, we can no longer get away with simple approximation when confronting theories with data so we will do it properly).

Thermal history of the Universe

Look through the slides below (or download the PDF if it doesn't show) to get an overview over the thermal history of the Universe (the new stuff starts at page 14) and some key events like:

For this course we won't do through and do the calculations for BBN (even though we have all the mathematical tools to do so - the Boltzmann language is great in that its super general, see Baumann Chapter 3.3.4 for this) so its enough to know the events mentioned above and some key predictions from it. This last point is the next topic and we will go through this in much more detail.

Page: /

Recombination

When the temperature of the primordial plasma starts to get lower than the binding energy of hydrogen then we can finally start making hydrogen atoms without them being broken up again imediately by energetic photons. This process in which the Universe goes from being a tightly coupled plasma of electrons and protons to having most of these elements tied up in hydrogen is called recombination. After this is completed then photons are able to travel freely without constantly scattering off electrons and the cosmic microwave background is released. For the numerical code we are making the most important thing is to learn how to compute the free electron fraction which is the basis for computing the optical depth and the visibility function.

The optical depth and the visibility function

Light that travels through a medium can be absorbed by the medium. If we have a source emiting an intensity $I_0$ then an observer a distance $r$ away from the source will observe an intensity $$I(r) = I_0 e^{-\tau(r)}.$$ The (dimensionless) quantity $\tau$ is called the optical depth. If $\tau \ll 1$ then the medium does nothing (we say its optically thin) and if $\tau \gg 1$ then we will see nothing (the medium is optically thick). The transition between these two regimes is when $\tau \sim 1$. In cosmology the main 'absorption' is Thompson scattering of photons of free electrons. The optical depth is defined as $$\tau(\eta) = \int_{\eta}^{\eta_0} n_e \sigma_T a d\eta'$$ and quantifies the probability for scattering a photon between some previous time and today (the number of scatterings of photons by electrons per unit time is $cn_e\sigma_T$). The components involved in this expression are: $n_e=n_e(\eta)$; the electron density (i.e. the number of free electrons per cubic meter) at time $\eta$, $\sigma_T = \frac{8\pi}{3}\frac{\alpha^2\hbar^2}{m_e^2c^2}$; the Thompson cross-section (dimension ${\rm m}^2$) and the scale factor $a$. The transition between the Universe being optically thick and optically thin happens around recombination when most of the free electrons are captured by free protons to form neutral hydrogen.

The expression for $\tau$ may also be written on an differential form (useful for when we are going to numerically solve it) which after restoring $c$ reads $$\frac{d\tau}{dx} = -\frac{c n_e \sigma_T}{H}$$ where $x = \log a$.

Another important quantity is the so-called visibility function, $$ \tilde{g}(x) = \frac{d}{dx}e^{-\tau} = -\frac{d\tau}{dx}e^{-\tau},\,\,\,\,\,\,\int_{-\infty}^{0} \tilde{g}(x)dx = 1. $$ This function is a true probability distribution, and describes the probability density for a given photon to last time having scattered at time $x$. As you will see, this function is sharply peaked around $x=-7$, corresponding to redshifts of $z\sim 1100$. The fact that this is sharply peaked is the reason why we often refer to the recombination period as the surface of last scattering: this process happened during a thin shell all around us, at a redshift of $z \sim 1100$.

To be able to compute the optical depth we need to know what the electron fraction is at all times which we will look at next.

Hydrogen Recombination: Saha approximation

The relevant interaction keeping eletrons and protons in equilibrium with photons is $$e^- + p^+ \rightleftharpoons H + \gamma$$ Let $n_i$ be the number-density of species $i$ then the Saha approximation we derived above tells us that as long as we are close to equilibrium then $$\frac{n_en_p}{n_H n_\gamma} = \left(\frac{n_en_p}{n_H n_\gamma}\right)_{\rm eq}$$ Lets massage this into a simpler form starting with the left hand side. The Universe is electrically neutral so we must have $n_e = n_p$. Photons are still in equilibrium with the primodial plasma at these times so $n_\gamma = n_\gamma^{\rm eq}$. If we assume that the only baryons is hydrogen (i.e. no helium and heavier elements) then $n_b \approx \frac{\rho_b}{m_H}$ where $m_H$ is the mass of hydrogen. The baryon density is given as $\rho_b = \rho_{c0}\frac{\Omega_{b0}}{a^3}$ where $\rho_{c0} = \frac{3H_0^2}{8\pi G}$ so given the cosmological parameters we can compute this. Finally we define the free electron fraction as $$X_e = \frac{n_e}{n_b}$$ where $n_b \approx n_p + n_H = n_e + n_H$ is the total baryon number density. When $X_e = 1$ then there is no hydrogen and when $X_e = 0$ all of the free electrons have been tied up in hydrogen atoms. Using all of this the Saha equation gives us $$n_b\frac{X_e^2}{(1 - X_e)} = \left(\frac{n_en_p}{n_H}\right)_{\rm eq}$$ Now lets deal with the right hand side. At the relevant temperatures the electrons, protons and hydrogen are non-relativistic so that $$n_i^{\rm eq} = \left(\frac{T m_i}{2\pi}\right)^{3/2}e^{-\frac{m_i}{T}}$$ which gives us $$\left(\frac{n_en_p}{n_H}\right)_{\rm eq} = \left(\frac{T m_em_p}{2\pi m_H}\right)^{3/2} e^{-\frac{m_e + m_p - m_H}{T}} $$ The expression in the exponential is nothing but the binding energy of hydrogen $\epsilon_0 = m_e + m_p - m_H \approx 13.6$ eV and $\frac{m_p}{m_H} \approx 1$ since the electron is much less massive which gives us the final result (with the constants $c,\hbar,k_b$ restored) $$\frac{X_e^2}{(1 - X_e)} = \frac{1}{n_b}\left(\frac{k_bT m_e}{2\pi \hbar^2}\right)^{3/2} e^{-\frac{\epsilon_0}{k_bT}}$$ This is a quadratic equation for $X_e$ that is easily solved (as we shall do numerically). We leave it as an exercise to deal with the case where we also have helium in the Universe (which we do). This is just slightly more complicated as we have three interactions to deal with instead of just the one above: $$e^- + {\rm He}^{++} \rightleftharpoons {\rm He}^{+} + \gamma$$ $$e^- + {\rm He}^{+} \rightleftharpoons {\rm He}^{0} + \gamma$$ $$e^- + p^{+} \rightleftharpoons H + \gamma$$ See the exercises for doing this calculation.

Hydrogen Recombination: Peebles equation

When $X_e$ moves away from $1$ then the Saha approximation fails to be accurate and we must resort to the full Boltzmann equation. This leads to what is known as the Peebles equation. To understand how this equation comes about lets start with the general form for thee Boltzmann equation for the number density for a $e^- + p^+ \rightleftharpoons H + \gamma$ process $$\frac{1}{a^3}\frac{d(n_e a^3)}{dt} = -\left<\sigma v\right> \left(n_e n_p - n_H n_\gamma \left(\frac{n_en_p}{n_Hn_\gamma}\right)_{\rm eq}\right)$$ Recall what we did when deriving the Saha approximation: $n_e = n_p$, $n_\gamma = n_\gamma^{\rm eq}$, introducing $X_e = \frac{n_e}{n_b}$ also noting that $n_ba^3 \approx$ constant. Doing this to the equation above gives us $$\frac{dX_e}{dx} = \frac{1}{H} \left((1 - X_e) \beta - n_b \alpha^{(2)} X_e^2\right)$$ where $\alpha^{(2)} = \left<\sigma v\right>$ and $\beta = \left<\sigma v\right>\left(\frac{k_bT m_em_p}{2\pi m_H}\right)^{3/2} e^{-\frac{m_e + m_p - m_H}{k_bT}}$. As a check of the calculation we can set the right hand side to $0$ and check that we obtain the Saha approximation (which we do). This wasn't too hard, but its a bit misleading as all the complications (the atomic physics - which we don't derive) is hidden inside the $\left<\sigma v\right>$ factor. Recombination in the Saha approximation is very easy; going beyond it suddenly gets much more complicated.

The first thing that complicates things is that there is no such things the hydrogen atom. The electron in the atom can be in many different energy and spin states (but will over time decay down to the ground state) and depending on the energy of the photon we can end up in any one of these states. For recombination we cannot recombine directly into the ground state of hydrogen. This process will emit a photon that has energy larger than the ground state (quantum number $n=1$) energy of hydrogen, $13.6$ eV, so even though we have created a new hydrogen atom this photon will likely ionize another atom so the net effect will be close to zero. This process is therefore too inefficient to recombine the Universe. Electrons therefore only efficiently recombine to the excited states of hydrogen, from which they decay quickly down to the first excited state with quantum number $n=2$. From this first excited state, electrons can reach the ground state in two ways: 1) decay from the $2p$ state by emitting a photon. This photon will likely hit another hydrogen atom in its ground state. However, cosmological redshifting systematically decreases the photon frequency, so there is a small chance that it escapes reabsorption if it gets redshifted far enough before encountering another hydrogen atom. 2) decay from the $2s$ state by emitting two photons. This process is very slow (for atomic physics), with a rate of $\Lambda_{2s\to 1s} = 8.22 s^{-1}$. What complicates things further is that atoms in the first excited state may also be reionized before they reach the ground state. When this is the case, it is as if the recombination to the excited state did not happen in the first place. To account for this possibility, Peebles defines a factor $C_r$ as the probability that an atom in the first excited state reaches the ground state through either of the two pathways described above before being photoionized. This is the basic elements of the Peebles model of recombination. It was first proposed in 1968 by Jim Peebles (Nobel prize in 2019) and Zel'dovich et al. in the USSR who independently computed the non-equilibrium recombination history of hydrogen for the first time. Putting it all together and doing the atomic physics calculations the resulting equation, the Peebles equation, reads $$\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} $$ where we have also taken into account that not all baryons are in fact hydrogen by introducing a new parameter $Y_p$: the primoridal helium (mass) fraction $Y_p = \frac{4n_{\rm He}}{n_H}$ which observations tell us is roughly $0.24$ and is a key prediction of BBN (a calculation we won't go through in this course - its a straight forward application of the theory we have already covered - but can be found in Baumann Chapter 3.3.4). This might look complicated, but the various constants are simply physical constants describing simple atomic physics. For even higher-accuracy many more states should be included, and also other atoms, most notably the helium atom. But here we will be satisfied with the Peebles' equation which is accurate at the $\sim {\rm a few}\%$ level which is more than good enough for us. The Peebles' equation is yet another linear first-order differential equation which is easy to solve numerically.

In the numerical code we first solve the Saha equation until its no longer valid ($X_e \sim 0.99$) and then we switch to the Peebles equation (we don't solve the Peebles equation very early on because the ODE is extremely stiff so its numerically unstable so we have to start with Saha). From this we obtain the free electron fraction and can use this to solve for $\tau$ and finally the visibility function.

Now that we have derived the relevant equations lets try to understand the solution physically. Lets start with Saha: $$\frac{X_e^2}{(1 - X_e)} = \frac{1}{n_b}\left(\frac{k_bT m_e}{2\pi \hbar^2}\right)^{3/2} e^{-\frac{\epsilon_0}{k_bT}}$$ When the temperature is very large the exponential is $\approx 1$, however the prefactor $\propto T^{3/2}$ is large so the right hand side is very large. Thus we must have $1-X_e \approx 0$ for the left hand side to be equally large (the free electron fraction $X_e$ is atmost $1$ - slightly larger if we include helium as its the number free electrons to protons). Thus $X_e \approx 1$ when the temperature is very large (very early in the Universe). This also makes sense: if the temperature is large and the density high then we have a lot of energetic photons around so whenever we form a hydrogen atom it will immediately get ionized. In the opposite limit when the temperature is very small then $\frac{\epsilon_0}{T}$ is very large and consequently the exponential $e^{-\frac{\epsilon_0}{k_bT}} \ll 1$ giving us $X_e \ll 1$. If the temperature is very low it means there are very few photons around with an energy large enough to inonize a hydrogen atom so the interaction will only go in the direction $e^- + p^+ \to H + \gamma$ and there will be less and less free electrons (and more and more neutral hydrogen). So what do we mean by high and low temperature? The critical temperature is the one that corresponds to the binding energy of hydrogen $T \sim \epsilon_0 \approx 13.6$ eV which again corresponds to a redshift $z \sim 60000$. However there is one thing we have ignored for this and that is the large size of the prefactor. The right hand side of the Saha equation can be written (see problem set) $$\frac{2\cdot 10^6}{\eta} \left(\frac{\epsilon_0}{T}\right)^{3/2} e^{-\frac{\epsilon_0}{T}}$$ where $\frac{1}{\eta} = \frac{n_\gamma}{n_b} \sim 10^9$ (there is a billion photons per baryons). For $T \sim 13.6$ eV the prefactor above is $\sim 10^{15}$ so recombination is delayed until the temperature is around $0.3$ eV corresponding to a redshift of about $z \sim 1200$. What is going on physically is that even though the mean energy of the photons are lower than the ionization energy of hydrogen when $T \lt 13.6$ eV there is still (remember that photons have a black-body distribution of energies) plenty of photons with energy larger than $13.6$ eV to ionize hydrogen atoms. Another feature of the Saha result that is worth noticing is that as the temperature gets lower and lower then the number of free electrons is exponentially suppressed so it predicts that all the electrons ends up in hydrogen atoms. However as we mention at some point the Saha approximation breaks down. When the abundance of free electrons and protons gets really low then the interaction is no longer efficient (its hard for a free electron to find a free proton to bond with). The interaction gets out of equilibrium and electrons eventually decouple from the primordial plasma and freezes-out, i.e. $X_e$ does not go to $0$ but stabilizes at something like $X_e \sim 10^{-3}-10^{-4}$. You can see a sketch of the evolution of $X_e$ below. To accurately capture this evolution we need to solve the Peebles equation numerically which is the task for Milestone 2.


Figure taken from Baumann. The evolution of the free electron from the early Universe till today.

Summary

See the slides below (PDF) for a brief summary of what we have gone through above.

Page: /