Read and process Planck CMB power spectra with `healpy`
In this tutorial we will load CMB power spectra generated by the cosmological parameters measured by Planck, we will create a realization of the Spherical Harmonics coefficients $a_{\ell m}$, save it to disk, and use them to create maps at different resolutions.
import pandas as pd
import healpy as hp
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# filter out all warnings
import warnings
warnings.filterwarnings("ignore")
print(plt.style.available)
plt.style.use('seaborn-talk')
First we go to the Planck data release page at IPAC, I prefer plain web pages instead of the funky javascript interface of the Planck Legacy Archive. The spectra are relegated to the Ancillary data section:
https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/
We can choose one of the available spectra, I would like the spectrum from theory, BaseSpectra, and download it locally:
!wget https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/cosmoparams/COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt
Now, I couldn't find any documentation about the format of those spectra, there are no units and we don't know if those spectra are just the $C_\ell$ or $\dfrac{\ell(\ell+1)}{2\pi} C_\ell$,
it is probably a standard output of the CAMB software, anyway I prefer to cross check with some public plot.
pandas has trouble with the # at the beginning of the title line, the easiest is to edit that out with a file editor.
input_cl = pd.read_csv("COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt",
delim_whitespace=True, index_col=0)
input_cl
input_cl.head()
input_cl.plot(logx=True, logy=True, grid=True);
We can compare this to one of the plots from the Planck wiki:
input_cl.TT.plot(logx=False, logy=False, grid=True)
plt.ylabel("$\dfrac{\ell(\ell+1)}{2\pi} C_\ell~[\mu K^2]$")
plt.xlabel("$\ell$")
plt.xlim([50, 2500]);
Very good, they agree, therefore the data are in $\dfrac{\ell(\ell+1)}{2\pi} C_\ell [\mu K^2] $, let's transform to $C_\ell [K^2]$
input_cl.head()
len(input_cl)
lmax = input_cl.index[-1]
print(lmax)
cl = input_cl.divide(input_cl.index * (input_cl.index+1) / (np.pi*2), axis="index")
cl /= 1e12
cl.head()
$\ell$ starts at 2, but all the healpy tools expect an array starting from zero, so let's extend that.
cl = cl.reindex(np.arange(0, lmax+1))
cl.head()
cl = cl.fillna(0)
cl.head()
Generate a CMB map realization
The power spectrum gives us the distribution of power with the angular scale, if we want to create a map we could use hp.synfast, however the realization will be different each time that we run synfast, we could always use the same seed to make it reproducible.
However, in case we want to have the same realization of a CMB map at different resolutions ($N_{side}$), it is better to generate a realization of the spectrum in Spherical Harmonics space, creating a set of $a_{\ell m}$, and then when needed transform them to a map with hp.alm2map.
We also want to make sure that the $a_{\ell m}$ have the same $\ell_{max}$ of the input $C_\ell$:
seed = 583
np.random.seed(seed)
alm = hp.synalm((cl.TT, cl.EE, cl.BB, cl.TE), lmax=lmax, new=True)
high_nside = 1024
cmb_map = hp.alm2map(alm, nside=high_nside, lmax=lmax)
We can double check that we got the order of magnitude correct, for example we can compare with one of the published Planck maps, the scale is $\pm 300 \mu K$:
hp.mollview(cmb_map[0], min=-300*1e-6, max=300*1e-6, unit="K", title="CMB Temperature")
Polarization is generally an order of magnitude lower:
for m in cmb_map[1:]:
hp.mollview(m, min=-30*1e-6, max=30*1e-6, unit="K", title="CMB Polarization")