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:

  • In the paper add a short description of the algorithms used and provide plots of $H(x)$, $\mathcal{H}(x)/H_0$, $\frac{1}{\mathcal{H}}\frac{d\mathcal{H}(x)}{dx}$, $\frac{1}{\mathcal{H}}\frac{d^2\mathcal{H}(x)}{dx^2}$, $\frac{\eta(x)\mathcal{H}(x)}{c}$, and one plot with $\Omega_b(x) + \Omega_{\rm CDM}(x)$, $\Omega_\gamma(x) + \Omega_\nu(x)$ and $\Omega_{\Lambda}(x)$ together.
  • Compute the time ($x$ and redshift) 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 do find the expressions for these times analytically and then use this to compute the numerical values. Give these numbers 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.
  • 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: $$h = 0.67,\,\,\,T_{\rm CMB 0} = 2.7255\,K,\,\,\,\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_{\rm eff} = 3.046,\,\,\,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$$ 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$.

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), $$ \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).

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

The final thing to compute is the age of the Universe. From $H = \frac{1}{a}\frac{d a}{dt} \to dt = \frac{da}{aH}$ we get $$t = \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}$$ In the radiation domainated 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.

Make a plot of $(\Omega_{\rm CDM}(x)+\Omega_b(x)), \Omega_\Lambda(x)$ and $(\Omega_\gamma(x)+\Omega_\nu(x))$ in the same figure and describe the figure. Also plot $H(x)$, $\mathcal{H}(x)$, $\frac{1}{\mathcal{H}}\frac{d\mathcal{H}(x)}{dx}$, $\frac{1}{\mathcal{H}}\frac{d^2\mathcal{H}(x)}{dx^2}$ and $\eta(x)$.

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

Here are some plots to compare your results to. We show $\mathcal{H}(x),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$.


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

Comoving distance: $$\chi(x) = \eta_0 - \eta(x)$$ where $\eta_0 = \eta(x=0)$.

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

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