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)$, $\frac{1}{\mathcal{H}}\frac{d\mathcal{H}(x)}{dx}$, $\frac{1}{\mathcal{H}}\frac{d^2\mathcal{H}(x)}{dx^2}$, $\eta(x)$, and one plot with $\Omega_b(x) + \Omega_{\rm CDM}(x)$, $\Omega_r(x) + \Omega_\nu(x)$ and $\Omega_{\Lambda}(x)$ together. If you are a Ph. d. student, you should also include masseless neutrinos in your computations. If you don't implement neutrinos just put $\Omega_{\nu 0} = 0$ or remove it alltogether from the code (same goes with curvature).
  • 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.
  • A transcript of the source code used for the evaluation.

The fiducial cosmology we are going to use in this course is (the same as in Callin (2006) - so that you can compare to the plots there): $$h = 0.7,\,\,\,T_{\rm CMB 0} = 2.725,\,\,\,\Omega_{r0} = 5.042\cdot 10^{-5},\,\,\,\Omega_{\rm CDM 0} = 0.224,\,\,\,\Omega_{\Lambda 0} = 1 - (\Omega_{k 0}+\Omega_{b 0}+\Omega_{\rm CDM 0}+\Omega_{r 0}+\Omega_{\nu 0}) = 0.72995,\,\,\,\Omega_{k 0} = \Omega_{\nu 0} = 0$$

If you are a PhD student then you also need to include neutrinos and our standard value is $\Omega_{\nu 0} = N_{\rm eff}\cdot \frac{7}{8}\left(\frac{4}{11}\right)^{4/3}\Omega_{r 0}$ with $N_{\rm eff} = 3.046$.

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

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.

Theoretical background

For a textbook that covers the theoretical background needed for this course see Dodelson "Modern Cosmology". There are also some very good lecture notes online for example Daniel Baumann's notes from Cambridge (which is basically a book) is very good (the presentation here is more theoretical level than Dodelson so it depends on what you like). I strongly reccomend both of these.

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 &= -dt^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_{r 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_r$, $\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_{r 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 you can choose not to include neutrinos and curvature, but implement it in general and just put $\Omega_{k 0} = \Omega_{\nu 0} = 0$ when running the code (in that way your code will be general). 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_r &= \rho_{b,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_r(a) &= \frac{\Omega_{r 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_{r 0}$ and $\Omega_{\nu 0}$ are given by $$\Omega_{r 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_{r 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 use an analytical approximation for $\eta(x_{\rm start})$.

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_{r 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.

Make a plot of $(\Omega_{\rm CDM}(x)+\Omega_b(x)), \Omega_\Lambda(x)$ and $(\Omega_r(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_{r 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 faild 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.

The age of the Universe: $$t = \int_0^a\frac{da}{aH} = \int_{-\infty}^x \frac{dx}{H}$$ For $x=0$ we get the age of the Universe today ($\sim 14$ billion years).

Scalefactor as function of time: Once you have $t(x)$ (from 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 = \eta_0 - \eta$$ 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}$$