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!
The deliverables are the following. 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 solve for the multipoles and the power-spectra.
- Plot of the transfer function $\Theta_\ell(k)$ for a few different values of $\ell$. Hint: use an "appropriate" $x$-axis (e.g. $k\eta_0$) (10p).
- Plot of the integrand of the $C_\ell$ integral $|\Theta_\ell(k)|^2/k$ for a few different values of $\ell$ (10p).
- Plot of the CMB power-spectrum. Optional (but reccomended as its useful for debugging): include the contribution to $C_\ell$ from each individual components in the source function (SW, ISW, doppler and polarization term; you simply set the other terms to zero in your source function and run the code as usual) (20p).
- Plot of the matter power-spectrum today with the equality scale $k_{\rm eq}$ marked in (10p).
- 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 (5p).
- Plot of the EE and TE polarization spectrum if you include polarization.
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 and the fourth term is a consequence of the angular dependence of Thompson scattering (photons are more likely to be scattered in certain directions than others).
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.