Harvard Astro215hf

Harvard Astro215hf - Lecture 4

Get some CMB power spectra

Now let's simulate a sky with a CMB power spectrum. We'll switch to flat sky to make life simples - i.e. Fourier Transforms instead of Spherical Harmonic Transforms.

First we need a spectrum from CAMB:

x=load('camb_spectrum.txt');

% first col is ell, then TT, TE, EE, BB spectra
% (throw out monopole and dipole to prevent touble later)
l=x(3:end,1); Cs_l=x(3:end,2:5);
plot(l,Cs_l(:,1));

% we want raw C_l - remove the conventional prefactor
pf=2*pi./(l.*(l+1));
C_l=Cs_l.*repmat(pf,[1,4]);
semilogy(l,C_l(:,1));

Make Gaussian random numbers with real FT

We want to setup a grid of Gaussian distributed random numbers, apply an appropriate scaling as a function of spatial frequency, and then Fourier transform them to get a simulated sky pattern. First part is easy:

imagesc(randn(8));

However in order to have a real Fourier transform the random numbers need to be complex and conjugate symmetric about zero frequency and some of them need to be purely real - Cf the \(a_{l,m}\)'s above. Here is a function randn_ft_of_real which delivers the necessary symmetries:

f=randn_ft_of_real(8)
imagesc(imag(f)==0)

The Fourier transform of these is simply real Gaussian random numbers:

i=fftshift(fft2(fftshift(f)))
imagesc(i)

The fft function assumes that the zero frequency is in the upper-left corner, with frequency increasing to the middle, and then wrapping around to maximally negative frequency and increasing again back to zero. I find it more convenient to keep zero frequency in the middle - hence the fftshift functions.

Define "axis data"

Let's make a convenience function which sets up "axis data" for both image and Fourier planes - calc_ad.m. Remember: when going from image to Fourier space (and vice versa) reciprocal of span is resolution, and vice versa. The Fourier space analog of angle \(\theta\) is conventionally denoted \(u\).

% a 20x20 chunk of sky with 256x256 pixels
ad=calc_ad(20,256);

% look at the radial distance from the origin in F plane
imagesc(ad.u_val,ad.u_val,ad.u_r)

Envelope the random numbers with sqrt power spectrum

Now take some random F-plane values and "envelope" them with TT power spectrum. It turns out that in the small sky approximation \(\ell = 2\pi u\).

% complex random modes with real FT
f=randn_ft_of_real(ad.npix);

% 2d grid of ell values
% - if we increase the image space resolution we increase the ell space span
imagesc(ad.l_val,ad.l_val,ad.u_r*2*pi); colorbar

% interpolate to these floating point ell values from integer ell values
% - this is just 1d power spectrum in 2d radial form
plot(l,sqrt(Cs_l(:,1)));
imagesc(ad.l_val,ad.l_val,interp1(l,sqrt(Cs_l(:,1)),ad.u_r*2*pi));

% envelope the random numbers - we really want C_l rather than Cs_l
fp=f.*interp1(l,sqrt(C_l(:,1)),ad.u_r*2*pi,[],0);
imagesc(ad.l_val,ad.l_val,real(fp))

Render a simulated CMB T sky

Now take the Fourier transform as we did earlier:

im=fftshift(fft2(fftshift(fp)));
imagesc(ad.t_val_deg,ad.t_val_deg,real(im)); colorbar; axis image

We can compare this to the maps shown in SPT-3G paper fig1 - although that is much smaller chunk of sky.

Now for polarization - E/B to Q/U

What the experiments measure are the x/y components of the polarization (pseudo)vectors over the sky - we call these Q/U (Stokes parameters). However the things which the theory predicts are the E/B gradient and curl of the polarization pattern. The transformation from Q/U to E/B is of course non-local in map space, but it turns out to be really simple in Fourier space. We simply rotate the Q/U F-modes into one another by the angle of the location of each mode in the \(u/v\) plane (analog of the \(x/y\) plane).

ad=calc_ad(10,256)

% make some random E-modes
f=randn_ft_of_real(ad.npix);

% envelope with the EE spectrum
Ef=f.*interp1(l,sqrt(C_l(:,3)),ad.u_r*2*pi,[],0);

% take the B-modes as zero
Bf=zeros(size(Ef));

% take the angle of each F-mode in the u/v plane
[u,v]=meshgrid(ad.u_val,ad.u_val);
chi=-atan2(v,u)+pi/2;
imagesc(ad.u_val,ad.u_val,chi); axis xy; colorbar

% rotate from E/B to Q/U
Qf=+Ef.*cos(2*chi)+Bf.*sin(2*chi);
Uf=+Ef.*sin(2*chi)-Bf.*cos(2*chi);

% render Q/U maps
Qi=fftshift(fft2(fftshift(Qf)));
Ui=fftshift(fft2(fftshift(Uf)));

% look at the result
imagesc(ad.t_val_deg,ad.t_val_deg,real(Qi)); colorbar; axis image
imagesc(ad.t_val_deg,ad.t_val_deg,real(Ui)); colorbar; axis image

Perhaps counter-intuitively Q/U maps for an all-E, zero-B sky pattern do have a preferred direction - Q is vert/horiz stripy, and U diagonal stripy. But this is just due to the direction that we select to measure the polarization pseudovector angle from (which in this case is the direction of the north celestial pole).

Before we move on make T map to match Q/U maps:

f=randn_ft_of_real(ad.npix);
Tf=f.*interp1(l,sqrt(C_l(:,1)),ad.u_r*2*pi,[],0);
Ti=fftshift(fft2(fftshift(Tf)));

Make scan trajectory on the sky

Our telescopes can only measure temperature differences - so we "scan" them across the sky - as we see in this speeded up movie. Let's make a fake raster map scanning trajectory:

nx=200; ny=33;

% make the out-back scan in azimuth (x) direction
% (imagine the telescope can turn around instantly)
x=linspace(-4,4,nx/2);
x=[x,fliplr(x)];
plot(x);

% make the steps in elevation (y) direction
y=linspace(-4,4,ny)';
plot(y);

% expand out to the full grid
x=repmat(x,[ny,1]);
y=repmat(y,[1,nx]);

% flip them and then vectorize
% (seems to be no way to vectorize along rows...)
x=x'; y=y'; x=x(:); y=y(:);

% overlay trajectory on T map
imagesc(ad.t_val_deg,ad.t_val_deg,real(Qi)); colorbar; axis image
hold on
plot(x,y,'r'); xlim([-5,5]); ylim([-5,5]);
hold off

Apply beam smoothing to input maps

If we were to sample simulated detector timestream off these maps there would be a lot of aliasing since there is structure in the map at scales smaller than the y step size. In reality this doesn't happen because of the instrument beam smoothing. Let's apply this smoothing in map space. First we make a little function called gauss and generate a little "thumbnail" beam:

% first in 1d
n=ad.npix/2+1;
xv=ad.t_val_deg(n-15:n+15);
plot(xv,gauss([1,0,0.5/2.35],xv));

% now in 2d
[xg,yg]=meshgrid(xv,xv);
r=sqrt(xg.^2+yg.^2);
imagesc(xv,xv,r)
beam=gauss([1,0,0.5/2.35],r);
imagesc(xv,xv,beam)

Now convolve the maps with beam:

% ensure convolution does not change normalization
beam=beam/sum(beam(:));

% same switch means output map same size as input
Ts=conv2(real(Ti),beam,'same');
Qs=conv2(real(Qi),beam,'same');
Us=conv2(real(Ui),beam,'same');
imagesc(ad.t_val_deg,ad.t_val_deg,Ts); colorbar; axis image

Make simulated timestream data along this trajectory

A linearly polarized detector scanning across the sky sees a mix of the T+Q/U maps depending on its orientation angle: \begin{equation} a(t) = T(\mathbf{r}(t)) + Q(\mathbf{r}(t))\cos(2\psi) + U(\mathbf{r}(t))\sin(2\psi) \end{equation} (assuming \(\psi\) constant with time). We can interpolate this off the maps using Matlab's nifty interp2 function:

% interpolate from maps
psiA=pi/4;

% make the map as seen by a detector at this angle
As=Ts+Qs*cos(2*psiA)+Us*sin(2*psiA);
imagesc(ad.t_val_deg,ad.t_val_deg,As); colorbar; axis image

at=interp2(ad.t_val_deg,ad.t_val_deg,As,x,y);

% now do the same for an orthogonal detector
psiB=psiA+pi/2;
Bs=Ts+Qs*cos(2*psiB)+Us*sin(2*psiB);
bt=interp2(ad.t_val_deg,ad.t_val_deg,Bs,x,y);

% plot A/B timestreams for first two out-back moves
plot(at)
hold on
plot(bt,'r')
hold off
xlim([0,400])

Add noise to timestream

We can easily add white noise.

at=at+3*randn(size(at));
bt=bt+3*randn(size(bt));

plot(at)
hold on
plot(bt,'r')
hold off
xlim([0,400])

Assignment

Reverse the process! Take the pair-difference timestream, bin it into grid of pixels and figure out how to invert to get reconstructed Q/U maps. Then go to Fourier space, go from Q/U to E/B, and take the power spectrum. Hints:

Make a bigger area input map such that you can have 128 scan rows spaced by 0.25 deg in your reconstructed maps.
To go from binned pair diff to Q/U each pixel needs to have been observed at multiple detector angles - ideally separated by \(\pi/4\). (This can be achieved by rotating the instrument or by having more than one pixel pair in the focal plane.) A 2x2 matrix inversion is then required per pixel to go from the binned observed quantities to Q/U.
Instead of doing the beam convolve in map space think about doing it in Fourier space by enveloping the random modes with an additional multiplier. Then on the way back one can divide out this beam convolution to undo its effect. (This of course causes things to blow-up at high \(\ell\) - especially if noise has been added.)

At a minimum I am looking for a well commented and neatly arranged file containing the Matlab code which starts from the power spectrum, goes to smoothed T/Q/U maps, samples the sim A/B timestreams, takes the pair sum/difference, bins this into a grid of pixels, inverts to recover smoothed Q/U maps, goes back to Fourier space, rotates from Q/U to E/B, undoes the beam smoothing, and finally takes the variance in annular bins to measure the power spectrum.

Better would be a web page or similar which intersperses the Matlab code with figures showing the state at appropriate points through the process.

Assignment is due 10am Monday March 7 by email to pryke@physics.umn.edu