Wednesday, February 1, 2017

Spherical Harmonics - learn to roll your own

It doesn't take much reading about the cosmic microwave background before you butt up against the concept of spherical harmonics (aka those funny looking $Y_l^m(\theta,\phi)$ guys) and if you're like me they threw you for a loop at first glance. Typically introduced in the context of studying anisotropies of the CMB they are essentially an alternate way of representing functions defined on the surface of a sphere. In this case the sphere we're interested in is the sky around the Earth and the function to be represented is the intensity of CMB radiation finally arriving at Earth after billions of years flying through space between photon decoupling during recombination and today. I won't go into too much detail about the physics of the CMB in this post but suffice it to say that these anisotropies are the definitive evidence for the Big Bang and cosmological inflation.

Spherical harmonics can be used to represent any arbitrary function as a (potentially infinite) list of coefficients that weight how much each harmonic contributes. If you're coming from an engineering background this may make more sense when thought of like an extension to the Fourier series where an arbitrary function can be decomposed into a (potentially infinite) weighted sum of sinusoids, each representing a specific frequency component of the original signal. Spherical harmonics are much the same but rather than using one-dimensional sinusoids defined along a circle the decomposition uses two-dimensional functions defined along the surface of a sphere. What is important to grok is that any signal can be totally defined solely in terms of the values of the weights for each harmonic.

To get a better understanding of how exactly this works in practice I wrote a small program to generate image maps to visually see how each unweighted harmonic varies across the surface of the sphere. The first step in doing so was getting a basic handle on the math involved. We'll start with the formal definition of our signal of interest $T(\theta,\phi)$ where $Y_l^m$ is a spherical harmonic (specific to degree $l$ and order $m$) and $a_l^m$ is the weighting factor for that specific harmonic.

$$T(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{+l} a_l^m Y_l^m(\theta,\phi)$$

We're interested in displaying each harmonic individually so we'll need the definition of the spherical harmonic $Y_l^m$ with normalization factor $n_l^m$.

$$Y_l^m(\theta,\phi) = n_l^m P_l^m(cos(\theta)) e^{i m \phi}$$
$$n_l^m = \sqrt{\frac{2 l + 1}{4 \pi} \frac{(l - m)!}{(l + m)!}}$$

Keep in mind $a_l^m$ and $e^{i m \phi}$ are both complex quantities even if the end result we're displaying will just be the real value.

$a_l^m Y_l^m(\theta,\phi) = a_l^m n_l^m P_l^m(cos(\theta)) e^{i m \phi}$
$a_l^m Y_l^m(\theta,\phi) = n_l^m P_l^m(cos(\theta)) a_l^m (cos(m \phi) + i\ sin(m \phi))$
$a_l^m Y_l^m(\theta,\phi) = n_l^m P_l^m(cos(\theta)) (Re(a_l^m) + i\ Imag(a_l^m)) (cos(m \phi) + i\ sin(m \phi))$
$Re(a_l^m Y_l^m(\theta,\phi)) = n_l^m P_l^m(cos(\theta)) (Re(a_l^m) cos(m \phi) - Imag(a_l^m) sin(m \phi))$

For the purposes of the image generation I'm simply using $Re(a_l^m) = Imag(a_l^m) = 1$.

$P_l^m$ is the associated Legendre polynomial. For this the easiest way to deal with the associated Legendre polynomials was to simply precompute them as I wasn't able to find a free C# library for computing them dynamically. Since I'm planning on generating 6000x3000 pixel images that are rectilinear projections of the surface of the sphere I'll need 3000 equally spaced samples. This is relatively simple to generate with numerical computing software such as GNU Octave or MATLAB. The script I used in Octave to generate the values for the first 50 associated Legendre polynomials $P_l^m$ of $cos(\theta)$ and export them to CSVs is below.

for i = 0:49 
csvwrite(strcat(strcat("legendres-",num2str(i)),".csv"), legendre(i, cos(0:pi/3000:pi*2999/3000)));

Note that this only gets you results for $m >= 0$ but we can extend this with the simple relation below.

$$P_l^{-m} = (-1)^m \frac{(l - m)!}{(l + m)!} P_l^m$$

That's all the math we're using! I wrote up a quick script that reads in the CSV output from the Octave script and uses it to calculate $Re(a_l^m Y_l^m(\theta,\phi))$ for each $l$ where $l \leq 50$ and each $m$ where $-l \leq m \leq l$ linearly sampling $\theta$ between $0$ and $\pi$ 3000 times and $\phi$ between $0$ and $2 \pi$ 6000 times. It then takes the resulting data and outputs an image where the minimum value is blue and the maximum value is red. The entire project is available here:

The whole purpose of doing this was to get a better handle on the concept of spherical harmonics by visual inspection so let's take a look at what we end up with!

$l = 0$, $m = 0$

$l = 1$, $m = 0$

$l = 1$, $m = 1$

$l = 2$, $m = 0$

$l = 2$, $m = 1$

$l = 2$, $m = 2$

$l = 3$, $m = 0$

$l = 3$, $m = 1$

$l = 3$, $m = 2$

$l = 3$, $m = 3$

A quick look at the results shows us that $l$ is related to the number of repetitions on the surface of the sphere and $m$ is related to the "orientation" of the repetitions. I've left out the images for negative $m$ as they aren't particularly illuminating but a quick look will confirm that they are simply the same as the positive $m$ results but with a rotation of $\frac{90°}{m}$ applied.

One last look at these, this time from the perspective of an observer inside the sphere looking outward as we do at the cosmos!