Harvard Astro215hf

Harvard Astro215hf - Lecture 1

Spherical Harmonic Multipoles

Any scalar pattern on a sphere \(T(\theta,\phi)\) can be expressed as a series of spherical harmonic multipoles, \begin{equation} T(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{+l} a_{l,m} Y_{l,m}(\theta,\phi) \end{equation} where the \(a_{l,m}\) are complex as are the \(Y_{l,m}(\theta,\phi)\). The latter are the spherical harmonics which are defined as, \begin{equation} Y_{l,m}(\theta,\phi) = \sqrt{\frac{1}{2\pi}} P_{l,m}(\cos\theta) \, e^{im\phi} \end{equation} where the \(P_{l,m}\) are the (real) associated Legendre polynomials (in Matlab "fully normalized" form).

(Note that the use of \(\cos(\theta)\) here dictates that \(\theta=0\) is perp to the plane of \(\phi\) - the usual physics convention for spherical polar coords - different to the astronomy/geography convention where Dec=0 is in the plane of Right Ascension.)

Generate single \(Y_{l,m}(\theta,\phi)\)

Going to Matlab let's generate the \(Y_{l,m}(\theta,\phi)\) patterns on a rectangular grid of spherical polar theta/phi coordinates:

% number of phi values and angular increment
nphi=90; da=(2*pi)/nphi;

% calc theta/phi values
theta=da/2:da:pi;
phi=da/2:da:2*pi;

% get normalized associated Legendre polynomials
l=4; P_lm=legendre(l,cos(theta),'norm');

% look at the output
plot(theta,P_lm); legend({'0','1','2','3','4'});

% select desired poly and multiply by complex exponential in phi
% note use of the single quote symbol to "flip" to column vector and
% that column times row = matrix
m=3; Y_lm=sqrt(1/(2*pi))*P_lm(m+1,:)'*exp(i*m*phi);

% plot on a rectangular grid
imagesc(phi,theta,real(Y_lm));

% plot on Mollweid projection
clf; axesm('MapProjection','mollweid'); pcolorm([-90,90],[-180,180],real(Y_lm))

% plot on sphere
clf; axesm('MapProjection','globe'); pcolorm([-90,90],[-180,180],real(Y_lm));
view(30,30)
axis off

One of the best features of Matlab is that one can extend the language very easily by writing new functions. If we make a file called multipole.m and place the following code in it:

function [Y_lm,theta,phi]=multipole(l,m,nphi)
%
% Calculate spherical harmonic multipole

% angular increment
da=(2*pi)/nphi;

% calc theta/phi values
theta=da/2:da:pi;
phi=da/2:da:2*pi;

% normalized associated Legendre polynomials
P_lm=legendre(l,cos(theta),'norm');

% select desired poly and multiply by complex exponential in phi
Y_lm=sqrt(1/(2*pi))*P_lm(abs(m)+1,:)'*exp(i*m*phi);

return

then magically a new command line function "multipole" is available:

% plot the l=4, m=+3 pattern
[Y_lm,theta,phi]=multipole(4,3,90);
subplot(2,1,1); imagesc(phi,theta,real(Y_lm)); subplot(2,1,2); imagesc(phi,theta,imag(Y_lm));

% the l=4, m=-3 pattern just flips the sign of imag part
[Y_lm,theta,phi]=multipole(4,-3,90);
subplot(2,1,1); imagesc(phi,theta,real(Y_lm)); subplot(2,1,2); imagesc(phi,theta,imag(Y_lm));

Note that the \(Y_{l,m}(\theta,\phi)\) is really two patterns in one container - the real and imaginary parts.

Find single \(a_{l,m}\)'s

So for a given set of \(a_{l,m}\)'s we can calculate \(T(\theta,\phi)\). But what about the reverse? Because the \(Y_{l,m}\)'s are orthogonal to one another the \(a_{l,m}\)'s are just \begin{equation} a_{l,m} = \int_0^{2\pi} \int_0^\pi T(\theta,\phi) \, Y_{l,m}^\ast(\theta,\phi) \sin(\theta) \, \mathrm{d}\theta \, \mathrm{d}\phi \end{equation} So we can evaluate the \(a_{l,m}\)'s by just generating the \(Y_{l,m}(\theta,\phi)\)'s one at a time and taking the sum of the product against some map \(T(\theta,\phi)\). For the purposes of illustration we will use a topographic map of the Earth which ships with matlab:

% get the data
load topo

% plot on a rectangular grid
imagesc(topo); axis xy

% look at the distribution - why is this a bit wrong?
hist(topo(:),100)

% see how big the grid is
whos topo

% make the multipole on the same size grid
[Y_lm,theta,phi]=multipole(4,3,size(topo,2));

% find the pixel area
[phig,thetag]=meshgrid(phi,theta);
da=theta(2)-theta(1);
pixarea=da^2*sin(thetag);

% integrate product of map and multipole over theta and phi -
% this gives the value of a_{4,3}
x=topo.*conj(Y_lm).*pixarea; a_lm=sum(x(:))

% look at the weighted map
imagesc(real(x)); axis xy

% integrate product of multipole and its conjugate -
% this checks the normalization - should be unity
x=Y_lm.*conj(Y_lm).*pixarea; sum(x(:))

Find all the \(a_{l,m}\)'s

Now let's automate and plot the \(a_{l,m}\)'s. First make a function decompose_direct.m:

function [a_lm,l,m]=decompose_direct(map,lmax)
%
% Decompose a pattern on the sphere to yield the a_lm
% coefficients by generating each Y_lm pattern, multiplying
% and integrating

% make non-existant a_lm's be NaN
a_lm=complex(NaN,NaN)*zeros(lmax+1,2*lmax+1);

for l=0:lmax
  for m=-l:l
    [l m]
    
    % get the multipole pattern
    [Y_lm,theta,phi]=multipole(l,m,size(map,2));
    
    % find the pixel area
    [phig,thetag]=meshgrid(phi,theta);
    da=theta(2)-theta(1);
    pixarea=da^2*sin(thetag);
    
    % take the integral of the product - note conjugation
    x=map.*conj(Y_lm).*pixarea; a_lm(l+1,lmax+m+1)=sum(x(:));
  end
end

% for plotting convenience
l=0:lmax; m=-lmax:lmax;

return

Using this look at the \(a_{l,m}\)'s.

load topo

% get the a_lm's
[a_lm,l,m]=decompose_direct(topo,4);

% look at a_lm's as number
% note imaginary part is zero for m=0 and +/-m are conjugates
% of one another
a_lm

% visualize
subplot(2,1,1); imagesc(m,l,real(a_lm)); colorbar
subplot(2,1,2); imagesc(m,l,imag(a_lm)); colorbar

Reconstruct the map from the \(a_{l,m}\)'s

Now let's rebuild the map one multipole at a time...


% get the a_lm's up to a higher l
[a_lm,l,m]=decompose_direct(topo,30);

% watch it reconstruct
topor=zeros(size(topo));
lmax=size(a_lm,1)-1;
for l=0:lmax
  for m=-l:l
    % get the multipole
    [Y_lm,theta,phi]=multipole(l,m,size(topo,2));
    topor=topor+Y_lm*a_lm(l+1,lmax+m+1);
  end
  clf; axesm('MapProjection','globe'); pcolorm([-90,90],[-180,180],real(topor));
  view(-160,30); colormap jet; gridm; axis off
  title(l)
  drawnow
end

% check normalization
subplot(2,1,1); imagesc(topo); axis xy; colorbar
subplot(2,1,2); imagesc(real(topor)); axis xy; colorbar

% check imaginary part small
subplot(2,1,2); imagesc(imag(topor)); axis xy; colorbar

A smarter pixelization - Healpix

The rectangular grid of spherical polar coordinates used above is clumsy - there are lots of very small pixels around the poles. A better option developed specifically for CMB data analysis is Healpix - also see the Healpix paper. Three main features are: i) equal area pixels ii) pixels arranged along equal latitude rings iii) resolution scales by splitting each quadrilateral cell into 4 starting from basic 12 cells.

Healpix was originally Fortran, and was then recoded as C/C++. There is a library of functions and then a few top-level programs built from these - pricipally anafast which goes from \(a_{l,m}\) to map and synfast which goes the other way. There is an incomplete IDL library. There is now a complete library called healpy which can be used interactively with ipython - this is probably the best bet for new users. Unfortunately there is no proper Matlab port. But let's play with a port of a simple low-level function pix2ang_ring.m which allows us to investigate the pixel layout on the sphere:

% resolution specified by "nside" which is power of 2
[theta,phi]=pix2ang_ring(1);

% plot flat
plot(phi,theta,'.');

% plot on sphere
clf; axesm('MapProjection','mollweid','angleunits','radians'); gridm
[theta,phi]=pix2ang_ring(1); plotm(pi/2-theta,phi,'b.');
[theta,phi]=pix2ang_ring(2); plotm(pi/2-theta,phi,'r.');
[theta,phi]=pix2ang_ring(4); plotm(pi/2-theta,phi,'g.');

% values of theta along rings
unique(theta)

% number of rings = nside*4-1
length(unique(theta))

% theta spacing of rings - clearly can't be equally spaced...
diff(unique(theta))

% number of phi points around each ring
u=unique(theta);
for i=1:length(u)
  n(i)=sum(theta==u(i));
end
clf; plot(n,'.-')
% so we have polar caps with varying number of pixels around rings
% and equatorial block with constant 4*nside around each ring