# 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.

Figure 1: The CMB power-spectrum and the physical effects that leaves its imprints on different scales.

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, conformal time and distance measures as function of scale factor and the variable $x = \log(a)$ which will be our main time variable in this course.

The deliverables are the following (for the old 2021 version see this PDF). There is 55 points availiable for the actual deliverables and 45 points for how this is presented/discussed in the report:

• The source code used for the evaluation. This should be possible for us to easily compile and run. A common way students have chosen to do this in previous years is to make a GitHub account, upload your code there and just provide the link. This is no requirement, you can just as easily include this together with your report in Devilry.
• A copy of your paper. This should have a section about this milestone which should contain: introduction (what we do and why this is relevant for the project), a short theory section with all relevant equations, tests showing the code gives expected results, the actual plots/results asked for and any other things you deem relevant and finally a proper description and discussion of the results tying it together with the physics.
• Implement all the required functions in the template: $H(x)$, $\mathcal{H}(x)$, $\frac{d\mathcal{H}(x)}{dx}$, $\frac{d^2\mathcal{H}(x)}{dx^2}$, $\eta(x)$, $\Omega_i(x)$ (for b,CDM,$\gamma$,$\nu$,$\Lambda$) and distance functions $d_L$ (luminosity distrance), $d_A$ (angular diameter distance) and $\chi$ (conformal distance). Solve ODEs and make splines for $\eta(x)$ and $t(x)$. When you are done with everything run the MCMC fits and produce plots.
• Plots that demonstrate that your code works properly: $\frac{1}{\mathcal{H}}\frac{d\mathcal{H}(x)}{dx}$, $\frac{1}{\mathcal{H}}\frac{d^2\mathcal{H}(x)}{dx^2}$ and $\frac{\eta(x)\mathcal{H}(x)}{c}$ compared with analytical expectations in different regimes. (10p)
• Plot of the conformal Hubble factor $\mathcal{H}(x)$. (5p)
• Plot of cosmic time $t(x)$ and conformal time $\eta(x)/c$. (10p)
• Plot of density parameters $\Omega_b(x) + \Omega_{\rm CDM}(x)$, $\Omega_\gamma(x) + \Omega_\nu(x)$ and $\Omega_{\Lambda}(x)$ together. (5p)
• 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). (5p)
• Plots from MCMC fits to supernova data. 1) the $1\sigma$ constraints from MCMC fits to supernova data in the $\Omega_\Lambda$ vs $\Omega_M$ plane 2) the posterior PDF of the Hubble parameter $H_0$. (10p)
• Remember to use sensible units in the plots (e.g. 100km/s/Mpc for $H$, Mpc/Gpc for distances and Gyr for times).
• Make a table showing times (both $x$, redshift $z$ and $t$) for when we have 1) radiation-matter equality 2) matter-dark energy equality and 3) when the Universe starts to accelerate and 4) the age of the Universe today and 5) the conformal time today ($\eta(0)/c$). Note 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 or solve numerically if you prefer this. Mark these as vertical lines in the plots where relevant (if it's relevant for discussing the plot). These values will be useful also in future milestones to understand/explain results. (10p)

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 neutrino perturbations, helium and the effects of reionization - but the rest of you should still keep neutrinos in the background). 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}$$

NB: Even if we are going to work with $\Omega_k = 0$ then you still need to implement curvature to work as this is required when doing the supernova fits.

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.

Page: /

# 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), $$$$\boxed{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}}},$$$$ 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$.) We also introduce a scaled Hubble parameter, $\mathcal{H}\equiv aH$ ("H prime" or "Hp" for short).

In this project we don't include curvature, but you still have to implement it in general. Master students don't have to deal with neutrinos later on, but do include it in the Hubble equation so use the correct value of $N_{\rm eff}$.

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 for the fluids we consider in this course). 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 \boxed{ \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 \boxed{ \begin{align} \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}, \end{align} } 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 \boxed{ \begin{align} \frac{d\eta}{dx} = \frac{da}{dx}\frac{d\eta}{da} = \frac{c}{\mathcal{H}} \end{align} } 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 and use the analytical approximation $\eta(x_{\rm start}) = \frac{c}{\mathcal{H}(x_{\rm start})}$ (we can solve the ODE 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 \boxed{ \begin{align} \chi = \eta_0 - \eta \end{align} } Evaluating the integral on the right the full equation above can then be written \boxed{ \begin{align} r &= \begin{cases} \chi \cdot \frac{\sin \left( \sqrt{|\Omega_{k 0}}| H_0 \chi /c \right)}{\left(\sqrt{|\Omega_{k 0}|}H_0\chi / c\right)} & \Omega_{k 0} < 0\\ \chi & \Omega_{k 0} =0 \\ \chi \cdot \frac{\sinh \left( \sqrt{|\Omega_{k 0}|} H_0 \chi / c\right)}{\left(\sqrt{|\Omega_{k 0}|}H_0\chi / c\right) } & \Omega_{k 0} > 0 \end{cases} \end{align} } 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$ so \boxed{ \begin{align} d_A &= \frac{D}{\theta} = ar \end{align} } 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 \boxed{ \begin{align} d_L = \frac{r}{a} = \frac{d_A}{a^2} \end{align} } 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.

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 \boxed{ \begin{align} \frac{dt}{dx} = \frac{1}{H} \end{align} } 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.

After everything is implemented we would like to use our code to get constraints on our cosmological parameters by comparing to data.

The data we will use us a set of supernova with associated redshift $z_i$, luminosity distance $d_L^{\rm obs}(z_i)$ and associated measurement errors $\sigma_i$. Under the assumption that the measurements are Gaussian distributed and uncorrelated between different redshifts the Likelihood function, telling us how well the data fits the theory, is given by $L \propto e^{-\chi^2/2}$ where the chi-squared function is $$\chi^2(h,\Omega_{m0}, \Omega_{k0}) = \sum_{i=1}^N \frac{[(d_L(z_i,h,\Omega_{m0},\Omega_{k0}) - d_L^{\rm obs}(z_i)]^2 }{\sigma_i^2}$$ A low value of $\chi^2$ means a good fit (high likelihood for the choice of parameters). What we want to do is to basicsally check all possible values of the parameters to find the best-fitting model and the range of parameters around it which are in agreement with the observed values. One can do this brute-force, but this step is in practice most commonly done by performing a so-called Monte Carlo Markov Chain to randomly sample from the likelihood (e.g. using the Metropolis algorithm as in the code template). The result from performing this step is a chain of values: parameter-points and likelihood values. The set of parameters with the lowest likelihood is our best-fit model. We can check if this is really a good fit by comparing the $\chi^2$ to the number of datapoints (here $N=31$). A good fit has $\chi^2/N \sim 1$ (a much higher value denotes a bad fit). When you have found the best-fit then you can find the $1\sigma$ ($68.4\%$) confidence region by looking at all values that satisfies $\chi^2 - \chi^2_{\rm min} \lt 3.53$ (see e.g. here). We can also plot the histogram of all the accepted samples to get the probabillity distribution (PDF) for a given parameter.

# 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 - and stores them. Make functions that are able to give you the cosmological parameters, the Hubble function and $\mathcal{H} = aH$ "Hp" plus the first two derivatives (these can be computed analytically) together with the different distance measures (comoving, luminosity and angular diameter distance).

Once this is done compute $\eta(x)$, spline the result and make a function that returns this function. 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.

Then finally compute the times asked for: equality times and acceleration time (compute this analytically) plus 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.

Once the distance measures are working you are ready to do the parameter fits to supernova data. First plot the results of the luminosity distance, using the fiducial cosmological parameters we use in the course, together with the data points from supernova observations (you might want to divide $d_L$ by $z$ when making the plot to better see the comparison; if so then you need to do the same with the data/errorbars!). To do the fitting you simply need to call the fitting routine provided in SupernovaFitting.h. If you don't use the template then you need to implement this yourself (it's a simple Metropolis MCMC sampler which you might have implemented in other courses before). See the header of this file if you have any problems making it work and for how to analyze the resulting chains to produce the plots asked for.

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. See the TODO's in BackgroundCosmology.cpp for hints on what to do in each routine.

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. You can use the info() method for this (just add whatever you want to print here). 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 evaluating the function you have made and check that they work and give reasonable results.

For testing $\eta$: you can 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 \simeq \frac{c}{\mathcal{H}}$. Plot $\eta \mathcal{H}/c$ and show that it converges to $1$ the further back in time you go. 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} \equiv 1$.

For testing $\mathcal{H}$ and it derivatives: you can derive an expression for what $\frac{\mathcal{H}^\prime}{\mathcal{H}}$ should be 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}}$.

Figure 2: Here are some plots to compare your results to. We show $\mathcal{H}(x),\eta(x)$ and $\frac{\eta(x)\mathcal{H}(x)}{c}$ and the density parameters for a toy-cosmology with parameters $\Omega_{m 0} = \Omega_{b 0} + \Omega_{\rm CDM 0} = 0.5$, $\Omega_{\Lambda 0} = 0.5$, $\Omega_{\nu 0} = 0$ and $h=0.7$. You can also run my code online here (down for the moment) to get results to compare to (NB: cannot guarantee there are no bugs in that code). Below we show some example results of fits to supernovadata.

# 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)$ (opposed 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$ (the non-relativistic limit) then $$B \approx y\int_0^\infty \frac{x^2}{e^{x}+1} dx$$ so $\rho_\nu \propto 1/a^3$. The amount of massive non-relativistic neutrinos (assuming all of them are non-relativitstic today) is $$\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.