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