# Table of contents

Back to the main page 1. Overview: Milestone I 2. Theoretical background 3. What you have to do 4. Testing your code 5. Appendix

# Overview: Milestone I

The goal of this whole project is to be able to predict the CMB (and matter) fluctuations - described by the so-called power-spectrum - from first principles and learn about all the different physical processes that goes on to be able to explain the results.

We will do this step by step and the topic of the first project is the evolution of the uniform background in the Universe. The goal here is to make a class/module that takes in the cosmological parameters and has functions for getting the Hubble parameter and the conformal time as function of scale factor and the variable $x = \log(a)$ which will be our main time variable in this course to compute the Hubble parameter as $x$, the logarithm of the scale factor $a$, and the conformal time $\eta$ as a function of $x$.

The deliverables are the following (edited for 2022, for the 2021 version see this PDF):

- In the paper add a short description of the algorithms used and provide plots of 1) $H(x)$ and $\mathcal{H}(x)$, 2) $\frac{1}{\mathcal{H}}\frac{d\mathcal{H}(x)}{dx}$ and $\frac{1}{\mathcal{H}}\frac{d^2\mathcal{H}(x)}{dx^2}$, 3) $\eta(x)$, 4) $\frac{\eta(x)\mathcal{H}(x)}{c}$, 5) $t(x)$ and 6) $\Omega_b(x) + \Omega_{\rm CDM}(x)$, $\Omega_\gamma(x) + \Omega_\nu(x)$ and $\Omega_{\Lambda}(x)$ together. Use sensible units (e.g. 100km/s/Mpc for $H$, Mpc/Gpc for distances and Gyr for times). Describe the plots using the theory we have been through.
- Implement the distance functions and add one plot of the luminosity distance versus redshift with data from supernova observations (taken from Betoule et al. 2014) overplotted (see header line for file-format).
- Compute times (both $x$, redshift $z$ and $t$) for when we have radiation-matter equality, matter-dark energy equality and when the Universe starts to accelerate (matter means baryons+CDM and radiation means photons+massless neutrinos). You can find the expressions for these times analytically and then use this to compute the numerical values. Give these numbers in a table in the report and also mark these as vertical lines in the plots. These values will be useful in future milestones to understand the evolution of perturbations.
- Compute the age of the Universe today and the conformal time today ($\eta(0)/c$) and explain the difference.
- The light we observe as the CMB was emitted around $z\sim 1100$ (next milestone) and it has a near uniform temperature in all directions. Compute the angular size of the horizon $\theta_{\rm CMB}$ at this redshift and discuss the implications of this result.
- A transcript of the source code used for the evaluation.

The fiducial cosmology we are going to use in this course - which is the one you are going to use for all the results in the paper - is the best-fit cosmology found from fits to Planck 2018 data: $$ \begin{align} h &= 0.67, \\ T_{\rm CMB 0} &= 2.7255\,K, \\ N_{\rm eff} &= 3.046, \\ \Omega_{\rm b 0} &= 0.05, \\ \Omega_{\rm CDM 0} &= 0.267,\\ \Omega_{k 0} &= 0, \\ \Omega_{\nu 0} &= N_{\rm eff}\cdot \frac{7}{8}\left(\frac{4}{11}\right)^{4/3}\Omega_{\gamma 0}, \\ \Omega_{\Lambda 0} &= 1 - (\Omega_{k 0}+\Omega_{b 0}+\Omega_{\rm CDM 0}+\Omega_{\gamma 0}+\Omega_{\nu 0}),\\ n_s &= 0.965, \\ A_s &= 2.1\cdot 10^{-9}, \\ Y_p &= 0.245, \\ z_{\rm reion} &= 8, \\ \Delta z_{\rm reion} &= 0.5, \\ z_{\rm He reion} &= 3.5, \\ \Delta z_{\rm He reion} &= 0.5. \end{align} $$ where $n_s,A_s$ will not be relevant until the last milestone and the last few parameters are only relevant for the next milestone (and only for PhD students that have to include neutrinos, helium and the effects of reionization). For master students you can take $N_{\rm eff} = 0$ and $Y_p = 0$. The radiation and neutrino density parameters follows from the above as $$\Omega_{\gamma 0} = 2\cdot \frac{\pi^2}{30} \frac{(k_bT_{\rm CMB 0})^4}{\hbar^3 c^5} \cdot \frac{8\pi G}{3H_0^2},\,\,\,\Omega_{\nu 0} = N_{\rm eff}\cdot \frac{7}{8}\cdot \left(\frac{4}{11}\right)^{4/3}\Omega_{\gamma 0}$$

One note about the $\Omega$'s: often in the litterature people use $\Omega_X$ where they really mean $\Omega_{X0}$. I have tried to be consistent with this, but there might be some omissions. To make it more clear I have also tried to use the convention to always include the argument (i.e. write $\Omega_X(a)$) if I mean the function and not the value today, but always double-check the equations and ask if something is not clear (I mention this as its been a common source of confusion).

If you haven't done so already read Callin (2006). It goes through all the things we are going to do in this course, it has all the equations and it has some notes on possible algorithms to use for the different tasks. It also have many plots that are useful to compare to (just remember that it uses a slightly different cosmology than our fiducial one). Another very good paper (though its written on a much higher level) is this one. A short introduction to the project can be found in this PDF (keynote) which is displayed below.

# Theoretical background

See the lecture notes for the theory you should know for this milestone.

Lets review the theory basics and you can find more details in the lecture notes mentioned above. The task in this project is to compute the expansion history of the
universe, and look at the uniform background densities of the various
matter and energy components. Let us first define the
Friedmann-Robertson-Walker metric (here for a flat space where $k=0$),
$$
\begin{align}
ds^2 &= -c^2 dt^2 + a^2(t) \left(dr^2 + r^2 (d\theta^2 + \sin^2\theta
d\phi^2)\right) \\
&= a^2(t) (-d\eta^2 + dr^2 + r^2 (d\theta^2 + \sin^2\theta
d\phi^2)
\end{align}
$$
or in Cartesian coordinates
$$
\begin{align}
ds^2 &= -c^2dt^2 + a^2(t)(dx^2 + dy^2 + dz^2)\\
&= a^2(t)(-d\eta^2 + dx^2 + dy^2 + dz^2)
\end{align}
$$
where $a(t)$ is the scale factor, which measures the size of the
universe relative to today ($a_0 = a(\textrm{today}) = 1$), and $\eta$
is called conformal time. One thing to note: it's called conformal *time*, but its usually given in units of length for it to have the same dimension as the spatial coordinates. The conversion factor for this is the speed of light $c$. In this course the conformal time is a distance (and the corresponding time is this distance divided by the speed of light). As we will be looking at phenomena that
varies strongly over a wide range of time scales, we will mostly be
using the logarithm of the scale factor, $x \equiv \ln a $, as our
time variable. A fifth time variable is the redshift, $z$, which is
defined as $1+z = a_0 / a(t)$.

Einstein's General Relativity describes how the metric evolves with time, given some matter and density components. The relevant equation for our purposes here is the Friedmann equation, which may be written - when we don't assume $k=0$ - on the following form (see the theory page for some notes on how to do this derivation), $$ \begin{equation} H = H_0 \sqrt{(\Omega_{b0}+\Omega_{\rm CDM 0})a^{-3} + (\Omega_{\gamma 0} + \Omega_{\nu 0}) a^{-4} + \Omega_{k 0} a^{-2} + \Omega_{\Lambda 0}}, \end{equation} $$ where $H\equiv \dot{a}/{a}$ is the Hubble parameter (and dot denotes derivatives with respect to physical time, $\dot{} = d/dt$), and $\Omega_{b 0}$, $\Omega_{\rm CDM 0}$, $\Omega_{\gamma 0}$, $\Omega_{\nu 0}$, and $\Omega_{\Lambda 0}$ are the present day relative densities of baryonic matter, dark matter, radiation, neutrinos and dark energy, respectively. A subscript $0$ denotes the value at the present time. The term $\Omega_{k 0} = -\frac{kc^2}{H_0^2}$ denotes curvatuve and acts in the Friedmann equation as if it were a normal matter fluid with equation of state $\omega = -1/3$. This term follows from the other density parameters which can be seen from taking $a=1$ to get $\Omega_{k 0} = 1 - (\Omega_{b 0}+\Omega_{\rm CDM 0}) - (\Omega_{\gamma 0} + \Omega_{\nu 0}) - \Omega_{\Lambda 0}$ (Recall that $\Omega_x = \rho_x / \rho_c$, where $\rho_c = 3H^2/8\pi G$.)

In this project we don't include curvature, but its good to just implement it in general and put $\Omega_{k 0} = 0$ when running the code (in that way your code will be general). Master students don't have to deal with neutrinos, but do include it in the Hubble equation and just set $N_{\rm eff} = 0$ so that you have it. We also introduce a scaled Hubble parameter, $\mathcal{H}\equiv aH$ ("H prime" or "Hp" for short).

The Friedmann equations also describe how each component evolve with
time
$$\dot{\rho} + 3H(\rho + P) = 0$$
where $P$ is the pressure. It's useful to define the *equation of state* $\omega \equiv P/\rho$ (which is constant). In terms of this the solution reads $\rho \propto a^{-3(1+\omega)}$. For cold dark matter and baryons we have $\omega = 0$, for relativistic matter (radiation and massless neutrinos) we have $\omega = 1/3$ and for a cosmological constant we have $\omega = -1$.

This gives us: $$ \begin{align} \rho_{\rm CDM} &= \rho_{{\rm CDM},0} a^{-3} \\ \rho_b &= \rho_{b,0} a^{-3} \\ \rho_\gamma &= \rho_{\gamma,0} a^{-4} \\ \rho_{\nu} &= \rho_{\nu,0} a^{-4} \\ \rho_{\Lambda} &= \rho_{\Lambda,0}. \end{align} $$ Here, quantities with subscripts 0 indicate today's values. The density parameters $\Omega_X(a) = \rho_X / \rho_c$ can be written $$ \begin{align} \Omega_{k}(a) &= \frac{\Omega_{k0}}{a^2H(a)^2/H_0^2}\\ \Omega_{\rm CDM}(a) &= \frac{\Omega_{\rm CDM 0}}{a^3H(a)^2/H_0^2} \\ \Omega_b(a) &= \frac{\Omega_{b 0}}{a^3H(a)^2/H_0^2} \\ \Omega_\gamma(a) &= \frac{\Omega_{\gamma 0}}{a^4H(a)^2/H_0^2} \\ \Omega_{\nu}(a) &= \frac{\Omega_{\nu 0}}{a^4H(a)^2/H_0^2} \\ \Omega_{\Lambda}(a) &= \frac{\Omega_{\Lambda 0}}{H(a)^2/H_0^2} \end{align} $$

Two of the density parameters above follows from the observed temperature of the CMB. We have that $\Omega_{\gamma 0}$ and $\Omega_{\nu 0}$ are given by $$\Omega_{\gamma 0} = 2\cdot \frac{\pi^2}{30} \frac{(k_bT_{\rm CMB 0})^4}{\hbar^3 c^5} \cdot \frac{8\pi G}{3H_0^2}$$ $$\Omega_{\nu 0} = N_{\rm eff}\cdot \frac{7}{8}\cdot \left(\frac{4}{11}\right)^{4/3}\Omega_{\gamma 0}$$ where $T_{\rm CMB 0} = 2.7255$K is the temperature of the CMB today and $N_{\rm eff} = 3.046$ is the effective number of massless neutrinos (slightly larger than $3$ due to the fact that neutrinos had not completely decoupled when electrons and positrons annhialate and the $0.046$ accounts for the extra energy pumped into the neutrinos). We define matter-radiation equality as the time when $\Omega_b + \Omega_{\rm CDM} = \Omega_\gamma + \Omega_{\nu}$ and the matter-dark energy transition as the time when $\Omega_b + \Omega_{\rm CDM} = \Omega_{\Lambda}$. The onset of acceleration is the time when $\ddot{a} = 0$.

Another crucial concept for CMB computations is that of the
*horizon*. This is simply the distance that light may have travelled
since the Big Bang, $t=0$. If the universe was static, this would
simply have been $ct$, but since the universe also expands, it will be
somewhat larger. Note that the horizon is a strictly increasing
quantity with time, and we can therefore use it as a time
variable. This is often called *conformal time*, and is denoted
$\eta$.

To find a computable expression for $\eta$, we note that $$\frac{d\eta}{dt} = \frac{c}{a}.$$ The left-hand side of this equation may be rewritten into $$\frac{d\eta}{dt} = \frac{d\eta}{da}\frac{da}{dt} = \frac{d\eta}{da} aH,$$ such that $$\frac{d\eta}{da} = \frac{c}{a^2H} = \frac{c}{a\mathcal{H}}$$ and $$\frac{d\eta}{dx} = \frac{da}{dx}\frac{d\eta}{da} = \frac{c}{\mathcal{H}}$$ This is a differential equation for $\eta$, that can either be solved numerically by direct integration, i.e. $\eta(x) = \int_{-\infty}^x \frac{c dx'}{\mathcal{H}(x')}$ , or by plugging the expression into a ordinary differential equation solver. The initial condition is $\eta(-\infty) = 0$. We can't integrate from $-\infty$ so in practice we pick some very early time $x_{\rm start}$ and set either set $\eta(x_{\rm start}) = 0$ or, better, use the analytical approximation $\eta(x_{\rm start}) = \frac{c}{\mathcal{H}(x_{\rm start})}$ (we can solve it analytically in the radiation domianted era).

We also need some distance measures. Consider the line-element in spherical coordinates $$ds^2 = -c^2dt^2 + a^2(\frac{dr^2}{1-kr^2} + r^2d\theta^2 + r^2\sin\theta^2d\phi^2).$$ Photons move on $0$-geodesics $ds^2 = 0$ so if we consider a radially moving photon ($d\theta=d\phi=0$), traveling from $(t,r)$ to us today at $(t=t_{\rm today},r=0)$ we get $$cdt = \frac{adr}{\sqrt{1-kr}}\implies \int_t^{t_{\rm today}}\frac{cdt}{a} = \int_{0}^{r}\frac{dr'}{\sqrt{1-kr'}} $$ The left hand side is what we call the co-moving distance and is closely related to the conformal time $$\chi = \eta_0 - \eta$$ Evaluating the integral on the right the full equation above can then be written $$ r = \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} $$ For a flat Universe we simply have $r = \chi = \eta_0 - \eta$. With this in hand we can compute all the standard distance measures we have in cosmology. If we know an object's physical size, $D$, and its angular size, $\theta$, as viewed from earth then the angular diameter distance is defined as $d_A = \frac{D}{\theta}$. From the line-element we see that the angular distance when we move in the $\theta$-direction is $$dD = a r d\theta \to d_A = \frac{D}{\theta} = ar$$ Note that for a flat Universe the angular diameter distance reduces to $d_A = a\chi = a(\eta_0 - \eta)$. The last distance we need is the luminosity distance. If we know the intrinsic luminosity of a source and measure its flux then we can define the distance to the source via $F = \frac{L}{4\pi d_L^2}$ which is simply $$d_L = \frac{r}{a} = \frac{d_A}{a^2}$$ We see that all distance measures are given directly from the conformal time, the scale-factor and the curvature so you only need to implement the functions $r(\chi)$, $\chi(x)$, $d_A(x)$, $d_L(x)$ as given above.

From the knowledge of the distances we can compute the angular size of the horizon $$\theta = \frac{\eta}{d_A}.$$ This represents the angular size of a patch of the sky which could have been in causal contact at some point. Said differenty: if two objects at the same redshift are separated by an angle greater than $\theta$ then no information could have been transmitted from one to the other (in the tradiational Big Bang picture without inflation).

The final thing to compute is the relation between cosmic time $t$ and our time-coordinate $x$. From $H = \frac{1}{a}\frac{d a}{dt} \to dt = \frac{da}{aH}$ we get $$t(x) = \int_0^a \frac{da}{aH} = \int_{-\infty}^x \frac{dx}{H(x)}$$ We can solve this just as we did for $\eta$ by evolving the ODE $$\frac{dt}{dx} = \frac{1}{H}$$ and splining the result. In the radiation dominated era we have $t(x) = \frac{1}{2H(x)}$ so the initial condition is $t(x_{\rm start}) = \frac{1}{2H(x_{\rm start})}$. Evaluating this at $x = 0$ (today) gives us the age of the Universe.

# What you have to do

Implement a class that takes in all the cosmological parameters $h,\Omega_{b 0},\Omega_{\rm CDM 0},\Omega_{\Lambda 0},T_{\rm CMB 0},N_{\rm eff}$ (and compute $\Omega_{\gamma 0}$, $\Omega_{\nu 0}$ and $\Omega_{k 0}$ from these. If you are a master student you can put $\Omega_{k 0} = 0$ and $\Omega_{\nu 0} = 0$. Make functions that are able to give you back the cosmological parameters, the Hubble function and $\mathcal{H} = aH$ "Hp" plus the first two derivatives of each (these can be computed analytically). This you will need a lot later on.

Once this is done compute $\eta(x)$, spline the result and make a function that returns this function. Also make a function that prints out the cosmological parameters you have taken in and computed. NB: the spline eta_of_x_spline is already defined in the header BackgroundCosmology.h so to create it all you have to do is to call eta_of_x_spline.create(x_array,eta_array,"eta_of_x") with your x_array and eta_array you got from the ODE solver.

When solving for $\eta$ its a good idea to first write the ODE as an ODE in $x$ - our main time-variable throught this course - and solve this.

Then compute the age of the Universe. You can do this in the same way as you did for the conformal time by solving an ODE and making a spline of $t(x)$ (useful for computing the time at any $x$ later on if you need that even though we only ask for the time today here). Evaluating the spline at $x=0$ gives you the time of the Universe in seconds (if you use SI units), but convert it to a more sensible unit like gigayears ($10^9\cdot 365 \cdot 24 \cdot 60\cdot 60$ seconds) when presenting the results.

For the distance measured we ask for you will need to implement the functions you need yourself (not in the template).

For some of the times we ask for you will need to derive analytically the expressions to implement.

If you choose to work with the C++ template take a look at BackgroundCosmology class. The file BackgroundCosmology.h contains the definition of this class, the variables (and splines) and functions it contains. Your task is to fill in all the TODO in BackgroundCosmology.cpp. This means computing derived quantities from the input ($\Omega_{\gamma 0}, \Omega_{\nu 0}$ and $H_0$ from $h$, $T_{\rm CMB 0}$ and $N_{\rm eff}$), solving the $\eta$ ODE and splining the result.

# Testing your code

The most important thing about coding is to test your code. There are always bugs and more bugs. The more tests you have the better. Here are some tips (which does not cover everything) about things you can do to ensure that your code is working correctly:

Print out the variables you entered (and the ones you computed) and check that they are set correctly inside the class. A more concrete test is to compute the sum of all the density parameters today. Check that this sum is $1$ (or very close to $1$). This tests fails if you failed to initialize the parameter correctly. Try all the function you have made and check that they work and give reasonable results (for example check that $\mathcal{H}(x) / (e^xH(x))$ gives you $1$).

Show analytically that in the radiation dominated era (when you can ignore energy in matter and dark energy so that $H^2 \propto 1/a^4$) then we have $\eta = \frac{c}{\mathcal{H}}$. Plot $\eta \mathcal{H}/c$ and show that it converges to $1$ the further back in time you go. You can derive similar approximations in the matter and dark energy dominated era (these will be on the form $\eta = \eta(a_i) + f(a)$ where $a_i$ is some scale-factor in the matter era for some function $f$ related to the Hubble factor). Matching these together gives you an approximation you can compare to your numerical results. The result will not match perfectly (unless you set all density parameters apart from radiation to zero), but it will give you an indication if the result is reasonable or way off in the general case. A more direct test is to compute and spline $\eta$ and use the spline-routines to get $\eta'$ and check that $\frac{\eta'\mathcal{H}}{c} = 1$.

To test the derivatives of $\mathcal{H}$ derive an expression for $\frac{\mathcal{H}^\prime}{\mathcal{H}}$ when the Universe is dominated by a fluid with equation of state $w$. Check that you get results that agrees with the expectations in the radiation, matter and dark energy dominated regimes (or run the code with only one energy component having non-zero $\Omega$ at the time and check that you get the expected result). Do a similar thing with $\frac{\mathcal{H}^{\prime\prime}}{\mathcal{H}}$.

To test $t(x)$ you can compare to the analytical result for when the Universe is dominated by a fluid with equation of state $w$ (see appendix). You can also compare to the age of the Universe when the CMB is released ($z\sim 1100$) - this should give $\mathcal{O}(400.000)$ years.

# Appendix

### Appendix: More things one can compute...

These things are optional, but here are some extra things one can compute.

**Scalefactor as function of time: ** Once you have $t(x)$ (that we compute above) you can spline $(t(x), x)$ (apposed to $(x, t(x))$) to get the inverse function $x(t) = \log(a(t))$. This be compared to the analytical approximation $a \propto t^{\frac{2}{3(1+w)}}$ for when the Universe is dominated by a fluid with equation of state $w$.

**Massive neutrinos:**
We treat neutrinos as massless, however oscillation experiments have revealed that neutrinos have a small mass $\lesssim 0.1$ eV (experiments are only sensitive to mass differences so we cannot tell exactly their masses - but the growth of structures is very sensitive to the sum of the masses). They are relativistic in the early Universe and non-relativistic close to today ($T_0 \sim 10^{-4}$ eV). Neutrinos decouple when they are still relativistic and consequently they maintain their distribution function $f\propto \frac{1}{e^{\frac{p}{T_\nu}}+1}$ with the temperature decreasing as $T_\nu = T_{\nu 0}/a$ where $T_{\nu 0} = (3.046/3)^{1/4}\cdot (4/11)^{1/3}T_{\rm CMB 0}$ (see problem set). The energy density of a single massive neutrino can be written
$$\frac{\rho_\nu}{\rho_{c0}} = \frac{\Omega_{\nu 0}}{a^4} \frac{B(y)}{B(0)}$$
where $\Omega_{\nu 0} = \Omega_{\gamma 0}\left(T_{\nu 0}/T_{\rm CMB 0}\right)^4$ and
$$B(y) \equiv \int_0^\infty \frac{x^2\sqrt{x^2 + y^2}}{e^{x}+1} dx$$
with $y(a) = \frac{m}{T_\nu(a)}$. If we have $N_\nu = 3$ neutrinos with mass $m_i$ the total energy density is likewise
$$\frac{\rho_\nu}{\rho_{c0}} = \frac{\Omega_{\nu 0}}{a^4} \frac{\sum_i\frac{B(y_i)}{B(0)}}{N_\nu}$$
When we have massless neutrinos $y\equiv 0$ and we get the expected result $\frac{\rho_\nu}{\rho_{c0}} = \frac{\Omega_{\nu 0}}{a^4}$ that we use in this project. On the other hand if $y\gg 1$ then
$$B \approx y\int_0^\infty \frac{x^2}{e^{x}+1} dx$$
which gives us that the amount of massive non-relativistic neutrinos (assuming all of them are non-relativitstic today)
$$\Omega_{m\nu 0} \approx \Omega_{\nu 0}\frac{\sum m_i}{T_{\nu 0}} \frac{\int_0^\infty \frac{x^2}{e^{x}+1}dx}{\int_0^\infty \frac{x^3}{e^{x}+1} dx}$$
and evaluating the integrals ($\frac{3\zeta(3)}{2}$ and $\frac{7\pi^4}{120}$) we get
$$\Omega_{m\nu 0}h^2 = \frac{16\zeta(3)}{11\pi}\left(\frac{(k_bT_{\rm CMB 0})^3G({1 \rm eV})}{\hbar^3 c^5 (H_0/h)^2}\right)\frac{\sum m_i}{{1 \rm eV}}\left(\frac{N_{\rm eff}}{3}\right)^{3/4} = \frac{\sum m_i}{93.05 eV} \left(\frac{T_{\rm CMB 0}}{2.725 K}\right)^3 \left(\frac{N_{\rm eff}}{3.046}\right)^{3/4}$$
so for $T_{\rm CMB 0} = 2.725$ K, $N_{\rm eff} = 3.046$ and $h=0.674$ we get $\Omega_{m\nu 0} \approx \frac{\sum m_i}{42 {\rm eV}}$. This estimate assumes the distribution function is exactly the Fermi-Dirac one, however there are some tiny distrortions produced when electrons and positrons anhiallate after neutrinos have almost (but not fully) decoupled and a much more detailed analysis gives $93.14$ eV instead of $93.05$ eV (and $94.12$ eV if we assumed instantanous decoupling with $N_{\rm eff} = 3$) so the results above is very accurate (to $0.2\%)$.

Note that this lets us constrain the mass of the neutrinos directly from what we know about the background. We know that atmost $\sim 0.3$ of the energy budget of our Universe is mass so it follows directly that $\sum m_i \lt 12$ eV. We can of course do much better that this with a more detailed analysis (the best constraints from cosmology today gives us $\sum m_i \lesssim 0.1-0.2$ eV), but its remarkable that such a simple theoretical calculation like we did above is enough to get a constrain that is of the same order of magnitude as particle physics experiments can currently provide! A common way of including neutrinos in a CMB code is to take $\Omega_{m\nu 0}$ as a free parameter, use this to compute $\sum m_i$ and compute the energy density from the above using $N_\nu$ neutrinos with mass $\sum m_i / N_\nu$. Alternatively take in the three masses directly. Also note that with massive neutrinos the total matter density parameter $\Omega_{M 0} = \Omega_{b 0} + \Omega_{\rm CDM 0} + \Omega_{m\nu 0}$. Why we don't bother to include the mass of the neutrinos will become more clear when we do perturbations - its significantly harder with massive neutrinos as we are forced to sample the full $p$-dependence of the distribution function and then integrate it to just get things like the number density or energy density perturbations.