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)
['Solarize_Light2', '_classic_test_patch', 'bmh', 'classic', 'dark_background', 'fast', 'fivethirtyeight', 'ggplot', 'grayscale', 'seaborn', 'seaborn-bright', 'seaborn-colorblind', 'seaborn-dark', 'seaborn-dark-palette', 'seaborn-darkgrid', 'seaborn-deep', 'seaborn-muted', 'seaborn-notebook', 'seaborn-paper', 'seaborn-pastel', 'seaborn-poster', 'seaborn-talk', 'seaborn-ticks', 'seaborn-white', 'seaborn-whitegrid', 'tableau-colorblind10']
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
--2020-10-06 17:27:01--  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
Resolving irsa.ipac.caltech.edu (irsa.ipac.caltech.edu)... 134.4.54.87
Connecting to irsa.ipac.caltech.edu (irsa.ipac.caltech.edu)|134.4.54.87|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 205647 (201K) [text/plain]
Saving to: ‘COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt.13’

COM_PowerSpect_CMB- 100%[===================>] 200.83K  --.-KB/s    in 0.1s    

2020-10-06 17:27:01 (1.92 MB/s) - ‘COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt.13’ saved [205647/205647]

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
TT TE EE BB PP
L
2 1016.7300 2.61753 0.030883 0.000002 5.013520e-08
3 963.7270 2.93806 0.039690 0.000004 6.099430e-08
4 912.6080 2.75866 0.034496 0.000006 7.025920e-08
5 874.4770 2.35185 0.023094 0.000009 7.829210e-08
6 848.5090 1.89605 0.012951 0.000013 8.530200e-08
... ... ... ... ... ...
2504 77.6127 -2.96537 2.910040 0.000000 0.000000e+00
2505 77.3908 -2.97153 2.913030 0.000000 0.000000e+00
2506 77.1691 -2.97827 2.916090 0.000000 0.000000e+00
2507 76.9485 -2.98537 2.919000 0.000000 0.000000e+00
2508 76.7311 -2.99293 2.921890 0.000000 0.000000e+00

2507 rows × 5 columns

input_cl.head()
TT TE EE BB PP
L
2 1016.730 2.61753 0.030883 0.000002 5.013520e-08
3 963.727 2.93806 0.039690 0.000004 6.099430e-08
4 912.608 2.75866 0.034496 0.000006 7.025920e-08
5 874.477 2.35185 0.023094 0.000009 7.829210e-08
6 848.509 1.89605 0.012951 0.000013 8.530200e-08
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()
TT TE EE BB PP
L
2 1016.730 2.61753 0.030883 0.000002 5.013520e-08
3 963.727 2.93806 0.039690 0.000004 6.099430e-08
4 912.608 2.75866 0.034496 0.000006 7.025920e-08
5 874.477 2.35185 0.023094 0.000009 7.829210e-08
6 848.509 1.89605 0.012951 0.000013 8.530200e-08
len(input_cl)
2507
lmax = input_cl.index[-1]
print(lmax)
2508
cl = input_cl.divide(input_cl.index * (input_cl.index+1) / (np.pi*2), axis="index")
cl /= 1e12
cl.head()
TT TE EE BB PP
L
2 1.064717e-09 2.741071e-12 3.234029e-14 1.904297e-18 5.250146e-20
3 5.046063e-10 1.538365e-12 2.078179e-14 1.904554e-18 3.193654e-20
4 2.867043e-10 8.666586e-13 1.083730e-14 1.904889e-18 2.207258e-20
5 1.831500e-10 4.925703e-13 4.836817e-15 1.905307e-18 1.639746e-20
6 1.269366e-10 2.836484e-13 1.937495e-15 1.905810e-18 1.276115e-20

$\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()
TT TE EE BB PP
L
0 NaN NaN NaN NaN NaN
1 NaN NaN NaN NaN NaN
2 1.064717e-09 2.741071e-12 3.234029e-14 1.904297e-18 5.250146e-20
3 5.046063e-10 1.538365e-12 2.078179e-14 1.904554e-18 3.193654e-20
4 2.867043e-10 8.666586e-13 1.083730e-14 1.904889e-18 2.207258e-20
cl = cl.fillna(0)
cl.head()
TT TE EE BB PP
L
0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
2 1.064717e-09 2.741071e-12 3.234029e-14 1.904297e-18 5.250146e-20
3 5.046063e-10 1.538365e-12 2.078179e-14 1.904554e-18 3.193654e-20
4 2.867043e-10 8.666586e-13 1.083730e-14 1.904889e-18 2.207258e-20

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)
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin

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