Harvard Astro215hf - Lecture 2

Planck Maps

The maps from Planck CMB space mission are available for download. However they are super high resolution (nside=2048) 2GByte files in Healpix format.

So for the purposes of illustration I have downgraded the 100GHz maps and mapped it to RA/Dec rectangular grid - available here: 

% get the data
load planck_100_longlat_0p5deg

% plot flat
clf; imagesc(phi,theta,p100.T); axis xy; colormap jet;
% astronomy convention is RA increases to the left - why?
set(gca,'XDir','reverse');

% galactic emission along plane dominates - truncate the color range
caxis([-300,300])

% plot mollweid
clf; axesm('MapProjection','mollweid'); pcolorm([-90,90],[-180,180],p100.T); caxis([-300,300])
set(gca,'XDir','reverse');
gridm
We can check that this last looks sensible by comparing it to the I map shown on this page. (Note that CMB people use I (intensity) and T (temperature) map interchangeably.)

We can now get the \(a_{l,m}\)'s of this map using our direct decomposition function from above: 

[a_lm,l,m]=decompose_direct(p100.T,4);
a_lm

Mask the galactic plane

If we take the power spectrum of the above map it will be corrupted by the bright emission along the galactic plane - mask it: 
% make a grid of theta values
[phig,thetag]=meshgrid(phi,theta);
imagesc(phi,theta,thetag); axis xy

% mask all the pixels along the plane
p100.T(abs(thetag)<0.3)=0;

% visualize
clf; axesm('MapProjection','mollweid'); pcolorm([-90,90],[-180,180],p100.T); caxis([-300,300])
set(gca,'XDir','reverse');
gridm

Interruption: Ligo Announcement

Class was interupted to watch the LIGO announcement live - just too cool to miss!

Extract power spectrum

The power spectrum is simply the expectation value of the \(a_{l,m}\)'s squared: \begin{equation} C_\ell = \langle a_{l,m} a_{l,m}^\ast \rangle = \langle | a_{l,m} |^2 \rangle = \mathrm{Var}(a_{l,m}) \end{equation} The decompose_direct function is too slow to extract multipoles up to the 100's that we desire so here is an alternate faster function decompose_fft.m - you can look inside to see how it works if you like but we won't worry about that now. 
[a_lm,l,m]=decompose_fft(p100.T);

% visualize
subplot(2,1,1); imagesc(m,l,real(a_lm)); colorbar
subplot(2,1,2); imagesc(m,l,imag(a_lm)); colorbar
For a Gaussian random field the real and imaginary parts of the \(a_{l,m}\)'s for a given \(l\) are just Gaussian random numbers with zero mean: 
subplot(2,1,1); hist(real(a_lm(end,:)),100)
subplot(2,1,2); hist(imag(a_lm(end,:)),100)
Look at the \(C_\ell\)'s: 
% nanvar takes the variance ignoring NaN values - the 2 says take variance
% along each row of the a_lm array
plot(l,nanvar(a_lm,[],2))
It is conventional to plot: \begin{equation} D_\ell = \frac{\ell(\ell+1)}{2\pi}C_\ell \end{equation} 
plot(l,nanvar(a_lm,[],2)'.*l.*(l+1)./(2*pi))

% why is this slighly different at low l (and more correct)?
plot(l,nanmean(a_lm.*conj(a_lm),2)'.*l.*(l+1)./(2*pi))
This looks roughly like what we expect - peak position agrees with table 7 of Planck 2015 I. Theory spectra come from CAMB program. Why do we expect the height of the peak to be a little supressed?

If your computer can handle it we can see the second peak as well - here is a map with 2x the resolution: 

load planck_100_longlat_0p25deg

p100.T(abs(thetag)<0.3)=0;

% takes 50 seconds on a fast computer...
tic; [a_lm,l,m]=decompose_fft(p100.T); toc

plot(l,nanvar(a_lm,[],2)'.*l.*(l+1)./(2*pi))