# Table of contents

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

# Overview: Milestone IV

Finally. After spending a lot of time on understanding the evolution of the background properties of the universe, its ionization history, and the growth of structure in the universe, we are now ready to pull it all together, and compute the key statistical observables in cosmology: the CMB and matter power spectrum!

*possible*to get for our fiducial cosmology (if everything is correct and the numerical accuracy is under full control). For PhD students its possible to get an almost perfect match while for the master students there will be some discrepancies due to not including helium, neutrinos, polarization and reionization (but you can cheat a bit and take $\ell \to \ell^{1.018}$ and $C_\ell \to C_\ell \exp({-0.05(\ell/200)^{1.5}})$ to model the effect of shifts of peaks due to neutrinos and damping due to helium and get a result that should match much better). Below we show the corresponding results for polarization that PhD students are going to compute.

The deliverables for this project are the following:

- In the paper define the quantities to be computed and add a short description of the algorithms used
- One plot of the transfer function $\Theta_\ell(k)$ for a few different values of $\ell$, one plot of the spectrum integrand, $\Theta_\ell(k)^2 / k$, for a few different values of $\ell$, and finally, the actual CMB power spectrum as function of $\ell$.
- It's also a good idea (but optional) to compute the CMB power-spectrum from each of the four terms in the source function (set the other terms to zero and run the code as usual).
- A plot of the matter power-spectrum today with the equality scale $k_{\rm eq}$ marked in.
- Data file with {$\ell$, $C_\ell$}. We will spend (at least) one lecture simply on comparing these from all of you, and talk about what these actually mean physically. Compare the results you get for the fiducial cosmology to the Planck spectrum (low $\ell$ TT data, high $\ell$ TT data, high $\ell$ EE data and high $\ell$ TE data - see headerline for fileformat, the power-spectrum $D_\ell$ in those files are $\ell(\ell+1)C_\ell/(2\pi)$ in units of $\mu K^2$). Master students only need to compare to the low $\ell$ data.
- A transcript of the module written for the evaluations

A short introduction to the milestone can be found in this PDF (keynote) shown below.

# Theoretical background

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

We now are in the position to actually compute the CMB power spectrum: we know how the ionization history of the universe, and we know how structure have grown from the very earliest epochs. Now we just have to pull it all together, in order to derive our preferred observable, namely the CMB power spectrum.

In order to understand how to get to the CMB power spectrum from our computed quantities, let us first recall the definition of the spherical harmonics transform of the CMB temperature field, $$T(\hat{n}) = \sum_{\ell m} a_{\ell m} Y_{\ell m}(\hat{n}),$$ where $\hat{n}$ is the direction on the sky, $a_{\ell m}$ are the spherical harmonics coefficients, and $Y_{\ell m}$ are the spherical harmonics themselves.

The CMB power spectrum, on the other hand, is simply defined as the expectation value of the square of the spherical harmonics coefficients, $$C_\ell \equiv \left\lt |a_{\ell m}|^2 \right\gt = \left\lt a_{\ell m}a_{\ell m}^* \right\gt.$$ Note that, in principle, this function should have two subscripts, $C_{\ell m}$, but because we assume that the universe is isotropic, it must have the same power spectrum towards both the $x$, $y$ and $z$ directions, and this implies full rotational invariance. As a result, there is no $m$ dependence in the power spectrum, and we simply average over $m$, and only call the spectrum $C_\ell$.

So, in order to get to the power spectrum, we need to know the temperature field we observe around us today, $T(\hat{n}, x=0)$. But fortunately, we already have this information (more or less), as we have already computed the evolution of the temperature field in the form of $\Theta_\ell(k,x)$. So all we have to do, is to read off the values of these functions at $x=0$, Fourier transform the resulting coefficients (to get $T(\vec{n})$ instead of $T(\vec{k})$), and read off the correct values.

But, you object, we only computed these functions up to $\ell = 6$,
and we are surely interested in smaller scales than that!. And you
are completely right: we want to know the power spectrum to at least
$\ell = 1200$. So what you have to then, is to rerun the code, but this
time with $\ell_{\rm max}=1200$ instead of $6$, and proceed. The only problem with
that is that it will take a *very* long time to complete.

This is where the method of Zaldarriaga and Seljak comes to our
rescue, through what they call the *line-of-sight integration*
approach. What we really need to know, is $\Theta(k,\mu,x=0)$. But
instead of first expanding the full temperature field in multipoles
and then solve the coupled equations, one can start by formally
integrating the original equation for $\dot{\Theta}$, and then expand
in multipoles at the end. For the details of this process, read
Section IVA in Callin (2006),
or Chapter 8 in Dodelson. The final expression is simply
$$\Theta_\ell(k, x=0) = \int_{-\infty}^{0} \tilde{S}(k,x) j_\ell[k(\eta_0-\eta)] dx,$$
where the *source function* is defined as
$$\tilde{S}(k,x) = \tilde{g}\left[ \Theta_0 + \Psi + \frac{1}{4}\Pi\right] + e^{-\tau} \left[\Psi^\prime-\Phi^\prime\right] -
\frac{1}{ck}\frac{d}{dx}(\mathcal{H}\tilde{g}v_b) + \frac{3}{4c^2k^2} \frac{d}{dx}
\left[\mathcal{H}\frac{d}{dx} (\mathcal{H}\tilde{g}\Pi)\right].$$
where you have computed all quantities in earlier milestones, and $\Pi
= \Theta_2 + \Theta_0^P + \Theta_2^P$. (If you skipped polarization in
the previous case, just set the two latter terms to zero.) For full
expressions for the last second derivative, see section IVA in Callin (2006)

So, the intuition behind this approach is that the CMB radiation we observe in a given direction on the sky, is basically the integral of the local CMB monopole (weighted by the visibility function) along the line of sight from us to infinity. That's the first $\Theta_0$ term in the source function. However, there are a number of corrections to this effect. First, the $\Psi$ term encodes the fact that the photons have to climb out of a graviational potential, and therefore loose energy on its way to us. $\Pi$ is a small quadrupolar (+ polarization) correction to the original monopole contribution.

The next main term is essentially the so-called Integrated Sachs-Wolfe contribution, which describes the fact that graviational potentials actually change while the photons are moving. The third term is a Doppler term.

So, this is a much more beautiful approach: Instead of evaluating what every single photon moment is at our position today, we can compute the monopole at all positions and times, and then do the integral through space. Much faster, and also quite intuitive.

The $j_\ell(x)$'s in the above expression are the so-called spherical Bessel functions, and take into account the projection of the 3D field (characterized by $k$) onto a 2D sphere (characterized by $\ell$).

But we're not quite done yet, even if we now know how to compute $\Theta_\ell(k)$ - we still need to go to actual $C_\ell$'s. This corresponds to 1) take the square of $\Theta_\ell(k)$ (since $C_\ell$ is the square of the Fourier coefficients), 2) multiply with the primordial power spectrum $P_{\rm primordial}(k)$ coming from inflation (recall that we originally set $\Phi \sim 1$ for all modes; this is now corrected by rescaling everything by $P(k)$ instead, which is perfectly valid, since all our equations are linear), and 3) integrate over all three spatial directions, instead of just the $z$ direction (but since we assume isotropy, we can use the same derived functions for all three directions). The CMB power spectrum then reads $$C_\ell = \frac{2}{\pi}\int k^2P_{\rm primordial}(k) \Theta_\ell^2(k)dk.$$ However, this can be massaged a bit further, by noting that most inflation models predict a so-called Harrison-Zel'dovich spectrum, for which $$\frac{k^3}{2\pi^2} P_{\rm primordial}(k) = A_s \left(\frac{k}{k_{\rm pivot}}\right)^{n_s-1},$$ where $n_s$ is the spectral index of scalar perturbations ($n_s \sim 0.96$), and expected to be close to unity, $k_{\rm pivot}$ is some scale for which the amplitude is $A_s$ (for $k_{\rm pivot} = 0.05/$Mpc we have $A_s \sim 2\cdot 10^{-9}$ for our Universe). The final expression for the spectrum is therefore $$C_\ell = 4\pi \int_0^{\infty} A_s \left(\frac{k}{k_{\rm pivot}}\right)^{n_s-1} \Theta_\ell^2(k) \frac{dk}{k}.$$ See this note for a full derivation.

Two final comments: First, the power spectrum is most often plotted in units of $\ell(\ell+1)/2\pi$ in $\mu \textrm{K}^2$, because it's overall trend is to drop as $\ell^{-2}$. It is therefore easier to see features when plotted in these units (i.e. you multiply $C_\ell$ by $\frac{\ell(\ell+1)}{2\pi}\left(10^6T_{\rm CMB 0}\right)^2$).

You should also compute the matter power-spectrum. This is really easy, you just need to multiply some numbers together: $$P(k,x) = |\Delta_M(k,x)|^2P_{\rm primordial}(k)$$ where $\Delta_M(k,x) \equiv \frac{c^2k^2\Phi(k,x)}{\frac{3}{2}\Omega_{M 0} a^{-1} H_0^2}$ and $\Omega_{M 0}$ is the total matter density parameter today. You only need to compute it for $x=0$, i.e. today. This is often plotted as $k$ in units of $h/\text{Mpc}$ versus $P(k)$ in units of $(\text{Mpc}/h)^3$. To explain this plot its useful to mark in the equality scale $k_{\rm eq} = \frac{a_{\rm eq} H(a_{\rm eq})}{c}$ in the plot where $a_{\rm eq}$ is the scale-factor for matter-radiation equality.

# What you have to do

Implement a class/module (PowerSpectrum.h if you use the C++ template) that takes in a BackgroundCosmology, RecombinationHistory and a Perturbations object - the ones you created in the previous three milestones - and use this to compute $\Theta_\ell(k,x=0)$ for $2 <\ell < 1200$ using line-of-sight integration and integrate this again to obtain the power-spectrum.

You will have to go back to what you did in milestone III and add a few more lines of code to make sure you compute and store the source-function defined above. This is what we need to do the line-of-sight integration. You will also need Bessel-functions (the function j_ell(ell, x) defined in Utils.h). If you haven't installed the ComplexBessel library you should probably do this (see the GitHub page for how to do this). The fallback if you don't have this is to use the GSL one and I found it to be quite buggy for large $\ell$.

There are several ways you can do the line of sight integration: 1) you can just evaluate the integral directly (see Callin for tips on this) or 2) you can evaluate the integral by treating it as an ODE and using the ODE solver (here you will have to be careful with the accuracy settings of the ODE solver), i.e. $$\frac{d\Theta_\ell(k, x)}{dx} = \tilde{S}(k,x) j_\ell[k(\eta_0-\eta)],\,\,\, \Theta_\ell(k, -\infty) = 0$$ Same goes for the integration of the $C_\ell$'s.

# Testing your code

# Appendix

### Appendix: More things to compute...

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

**Generate a CMB map: ** From the theoretical $C_\ell$ spectrum we can generate a CMB map (as shown at the beginning of this document). Generating a map is a fairly straight forward thing to do, but the hard part of this is to pixelate the sphere and project that down to a 2D map that we can plot. Luckily there is a great library for this namely the HEALPIX library (for this you also need the CFITSIO library). There is also a Python wrapper for this library called Healpy. You can supply that library with the $C_\ell$'s or a realization of the $a_{\ell m}$'s (which you can generate from the $C_\ell$'s) and have it provide you the map.

Healpix has routines for generating maps directly from a power-spectrum or from $a_{\ell m}$'s. If you want to generate your own realization of $T(\hat{n})$ then note that the $a_{\ell m} = x + iy$ with $x,y$ being random numbers drawn from a gaussian distribution with variance $C_\ell/\sqrt{2}$. A simple way to generate a realization is to draw (for each value of $\ell,m$) two random numbers $A,\theta$ from a uniform distribution on $[0,1)$ and set $$a_{\ell m} = \sqrt{-\log(A)C_\ell}e^{2\pi i \theta}$$ Note that $T = T^*$ is real so $$a_{\ell, -m} = (-1)^m a_{\ell m}^*$$ and we only need to generate them for $m\ge 0$. The other way around: if you want to estimate the angular power-spectrum from a set of $a_{\ell m}$'s by computing $$\hat{C}_\ell = \frac{1}{2\ell + 1}\sum_{m=-\ell}^\ell |a_{\ell m}|^2$$ This should give you back the input power-spectrum for large $\ell$ (and will suffer from cosmic variance for small $\ell$).

**Neutrino power-spectrum: ** We can't observe this, but you can still compute it. The computation is the same as for photons and polarization just with a different source function that is (ignoring any quadrupole moments) given by
$$\tilde{S}_\nu = (\mathcal{N}_0 + \Psi)\delta(\eta) + \frac{d\Psi}{dx} - \frac{d\Phi}{dx}$$
so its just a simple SW + ISW term without the supression of the optical depth since neutrinos have been free streaming since decoupling.

**Lensing potential: ** The source function for the CMB lensing potential is given by
$$\tilde{S}_\Psi = -\frac{2c\Psi}{\mathcal{H}(\eta_0 - \eta)} W(\chi,\chi_*)$$
$$W(\chi,\chi_*) = \frac{S_k(\chi - \chi_*)}{S_k(\chi_*)}\,\,\,\text{for}\,\, x \geq x_*\,\,\text{and zero otherwise}$$
where $x_*$ is the last scattering time ($\tau(x_*) = 1$), $\chi = \eta_0 - \eta$ and $S_k(\chi) = \chi$ for a flat Universe (see appendix of milestone 1 for a general expression when the curvature is non-zero). The usual normalization when plotting it is $\frac{\ell^2(\ell+1)^2}{2\pi}C_\Psi$.

**Angular correlation function: ** The angular correlation function $C(\theta)$ (correlations between two directions $\hat{n}_1,\hat{n}_2$ on the sky separated by an angle $\cos\theta = \hat{n}_1\cdot \hat{n}_2$):
$$C(\theta) = \frac{1}{4\pi}\sum_{\ell} (2\ell+1)C_\ell \mathcal{P}_\ell(\cos\theta)$$
where $\mathcal{P}_\ell$ is the Legendre polynomials.

**Correlation function: ** The Fourier transform of the matter power-spectrum
$$\xi(r) = \int_0^\infty \frac{\sin(kr)}{kr} \frac{k^3P(k)}{2\pi^2} \frac{dk}{k}$$
An aternative to computing the integral above is to write the integral as
$$r\xi(r) = \int_{-\infty}^\infty f(k)e^{ikr}dk$$
where $f(k) = \frac{kP(k)}{4i\pi^2}$ for $k\gt 0$ and $f(k) = f^*(-k) = -f(-k)$ for $k\lt 0$ which allows us to compute it by performing a 1D numerical Fourier transform (for this see FFTLog).

**Lensed CMB spectrum: ** Compute the effects of gravitational lensing on the CMB power spectrum (see Lewis & Challinor (Eq. 9.12 + 9.16-18) or Seljak).

The lensed correlation function is $$C^{\rm lensed}(\theta) \approx \frac{1}{4\pi}\sum_{\ell} (2\ell+1)C_\ell e^{-\ell(\ell+1)\sigma^2(\theta)/2} (\mathcal{P}_\ell(\cos\theta) + \frac{\ell(\ell+1)}{2}C_{\rm gl,2}(\theta) d^\ell_{1,-1}(\theta))$$ where $$\sigma^2(\theta) = C_{\rm gl}(0) - C_{\rm gl}(\theta)$$ $$C_{\rm gl}(\theta) = \sum_\ell \frac{2\ell + 1}{4\pi}\ell(\ell+1)C_\ell^\Psi d_{1,1}^\ell(\theta)$$ $$C_{\rm gl,2}(\theta) = \sum_\ell \frac{2\ell + 1}{4\pi}\ell(\ell+1)C_\ell^\Psi d_{-1,1}^\ell(\theta)$$ and $d^\ell_{mn}$ is reduced Wigner functions $$d^\ell_{mn}(\theta) = \sum_{i} (-1)^i\frac{[(\ell+m)!(\ell-m)!(\ell+n)!(\ell-n)!]^{1/2}}{(\ell+m-i)!(\ell-n-i)!i!(i+n-m)!}\cos^{2\ell+m-n-2i}(\theta/2)\sin^{2i+n-m}(\theta/2)$$ and the sum is over all terms where the factorials are non-negative and $\cos^2(\theta/2) = \frac{1+\cos(\theta)}{2}$ and $\sin^2(\theta/2) = \frac{1-\cos(\theta)}{2}$. From this the lensed $C_\ell$'s are given by $$C_\ell^{\Theta\,\rm lensed} = 2\pi \int_{-1}^1 C^{\rm lensed}(\theta) P_\ell(\cos\theta) d\cos\theta$$ Similar expressions apply for polarization, see Lewis & Callinor.