Spherical harmonics in Python¶
Author: Eduardo Martín Calleja
My interest on the spherical harmonics is motivated, as I told in my previous post because these are essential mathematical objects to solve many problems, especially for the analysis of the anisotropy of the cosmic microwave background. In this post I will study what spherical harmonics are, and I will show examples of how to calculate them with the Python SciPy library, as well as a series of graphic representations of the same.
This article is written entirely using the IPython Notebook.
Imports and references¶
%matplotlib inline
from __future__ import division
import scipy as sci
import scipy.special as sp
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm, colors
# This IPython magic generates a table with version information
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information scipy, numpy, matplotlib
References:
Chapter XI. "Great Balls of PDEs" from the Evans Harrell y James Herod WWW textbook : Linear Methods of Applied Mathematics
Chapter XII from the Peter J. Olver's book Introduction to Partial Differential Equations
A quick summary of the definition and properties of spherical harmonics, along with some graphics in Clem Pryke's Logbook
Theoretical derivation of the definition of spherical harmonics¶
Separation of variables¶
We intend to seek solutions of the Laplace equation:
$$\nabla^2 u = 0$$when it is expressed in spherical coordinates. We already saw in this previous post that the problem can be formulated as follows:
Being $(\rho, \theta, \phi)$ the radius, colatitude (from $0$ to $\pi$), and longitude (from $0$ to $2 \pi$) respectively, It is about finding functions $u(\rho, \theta, \phi)$ satisfying:
$$ \nabla^2 u = \frac{\partial^2 u}{\partial \rho^2} + \frac{2}{\rho} \frac{\partial u}{\partial \rho} +\frac{1}{\rho^2} \frac{\partial^2 u}{\partial \theta^2} + \frac{\cos(\theta)}{\sin(\theta) \rho^2} \frac{\partial u}{\partial \theta} + \frac{1}{\rho^2 \sin^2 (\theta)} \frac{\partial^2 u}{\partial \phi^2} = 0 $$To find the solution, we apply the separation of variables method, looking for now solutions that are the product of two functions, one that depends only on $ \rho $ and one that depends only on the angular variables $ (\theta, \phi ) $:
$$u(\rho, \theta, \Phi) = R(\rho) Q(\theta, \phi)$$Substituting calculating derivatives, and dividing by$u = RQ$ we come to:
$$ 0 = \frac{\nabla^2 u}{u} = \frac{R''}{R} + \frac{2}{\rho} \frac{R'}{R} + \frac{1}{\rho^2 Q} \underset{\nabla^2_S Q}{\underbrace{\left[ \frac{\partial^2 Q}{\partial \theta^2} + \frac{\cos(\theta)}{\sin(\theta)} \frac{\partial Q}{\partial \theta} + \frac{1}{\sin^2(\theta)} \frac{\partial^2 Q}{\partial \phi^2} \right]}} $$Where we recognize the expression in brackets as the spherical laplacian $\nabla^2_S Q$, in accordance with this post, that is: the Laplacian operator on a spherical surface of unit radius, which therefore only depends on the coordinates $(\theta, \phi)$. Then multiplying by $ \rho^2 $ and grouping variables on each side of the equation, we arrive at:
$$ \frac{1}{R(\rho)} \left[ \rho^2 R''(\rho) + 2 \rho R'(\rho) \right] = - \frac{\nabla^2_S Q(\theta, \phi)}{Q(\theta, \phi)} $$As the two members depend on different variables, both have to be equal to a constant $\mu$:
$$ \rho^2 R''(\rho) + 2 \rho R'(\rho) - \mu R(\rho) = 0 $$$$ \nabla^2_S Q(\theta, \phi) + \mu Q(\theta, \phi) = 0 $$Focusing for now on the second equation, it is an eigenvalue equation $- \nabla^2_S Q = \mu Q$, where $\mu$ is an eigenvalue of the operator $-\nabla^2_S$. This equation is known as the spherical Helmholtz equation. The solutions we are looking for are then eigenvectors of the spherical laplacian operator. We call this $Q(\theta, \phi)$ eigenvectors the spherical harmonics.
Solving the Helmholtz equation¶
Expanding the Helmholtz equation, we can write:
$$\frac{\partial^2 Q}{\partial \theta^2} + \frac{\cos(\theta)}{\sin(\theta)} \frac{\partial Q}{\partial \theta} + \frac{1}{\sin^2(\theta)} \frac{\partial^2 Q}{\partial \phi^2} + \mu Q = 0 $$For its study, we will apply again the separation of variables, looking for solutions of the form:
$$Q(\theta, \phi) = \Theta(\theta) \Phi(\phi) $$Substituting, dividing by $Q = \Theta \Phi$ and multiplying by $\sin^2(\theta)$ we obtain:
$$ \frac{1}{\Theta} \left[ \sin^2(\theta) \Theta'' + \sin(\theta) \cos(\theta) \Theta' \right] + \mu \sin^2 \theta = - \frac{\Phi''}{\Phi}$$Both terms will equal to a constant, let's call it $\lambda$, so that the spherical Helmholtz equation can be divided into two equations, such as:
$$\sin^2(\theta) \Theta'' + \sin(\theta) \cos(\theta) \Theta' + (\mu \sin^2 \theta - \lambda) \Theta = 0$$$$\Phi'' = - \lambda \Phi $$Let's look first to the second equation. If we look for solutions in the form of complex functions (which will greatly simplify trigonometric expressions) we get solutions of the following kind:
$$\Phi(\phi) = A e^{i \sqrt\lambda \phi} = A \left[ \cos(\sqrt \lambda \: \phi) + i \sin(\sqrt \lambda \: \phi) \right] $$Being $A$ a complex constant.
But, $\Phi(\phi)$ must be a periodic function of period $2 \pi$, that is, $\Phi(\phi) = \Phi(\phi + 2 \pi)$, and we can see that, by increasing $\phi$ in $2 \pi$ the argument of the sine and cosine functions increases by $2 \pi \sqrt \lambda$. so the only possibility that the periodicity is fulfilled is that $\sqrt \lambda$ is an integer value $m$, positive or negative. $\lambda = m^2$, and the solutions of the azimuthal component of the Helmholtz equation will be the product of a complex constant by:
$$\Phi^m(\phi) = e^{i m \phi} $$Now we will replace the value of $ \lambda $ in the first equation:
$$\sin^2(\theta) \Theta'' + \sin(\theta) \cos(\theta) \Theta' + (\mu \sin^2 \theta - m^2) \Theta = 0$$And we make the change of variables:
$$z = \cos(\theta) \quad , \quad \Theta(\theta) = P(\cos(\theta)) = P(z)$$And so the trigonometric functions of the equation vanish, leaving:
$$(1-z^2) \frac{d^2 P}{d z^2} - 2z \frac{dP}{dz} + \left( \mu - \frac{m^2}{1-z^2} \right) P = 0 $$Let us now consider some properties of the solutions of this equation: in the first place, $z = \cos(\theta)$, and so $-1 \leq z \leq 1$, and the function $P(z)$ only need to be defined in the interval $[-1, 1]$. On the other side, $\Theta(\theta)$ must be defined in $\theta = 0$ and in $\theta = \pi$ so $P(z)$ must be bounded in $z=-1$ and in $z=1$:
$$| P(-1) | < \infty \quad , \quad | P(1) | < \infty$$If we now assume $m$ non-negative integer: $m = 0, 1, 2, \dots$, it holds that (we will admit it without proof):
- The only values of $ \mu $ for which the above equation admits nontrivial solutions (ie, eigenvalues) are of the form $\mu = l(l+1)$ with $l$ integer and $ l \geq |m|$, resulting in the following equation, known as associated Legendre differential equation:
- The solutions compatible with the above conditions are called associated Legendre functions of the first kind, with degree $l$ and order $m$, which are usually designated by $P_l^m(z)$ (they are also called Ferrers functions). There are also the associated functions of the second kind $Q_l^m(z)$ that need not be considered now.
The simplest case is the one corresponding to $ m = $ 0 in which case the equation is called the Legendre differential equation:
$$(1-z^2) \frac{d^2 P}{d z^2} - 2z \frac{dP}{dz} + l(l+1) P = 0 $$And in this case the corresponding Legendre functions of order $ m = $ 0 coincide with the Legendre polynomials: $P_l^0(z) = P_l(z) \:, \: l=0, 1, \dots$.
Returning to the functions $\Theta(\theta)$, as $\Theta(\theta) = P(\cos \theta)$ we obtain the following functions as solutions of the component of the solution of the Helmholtz equation that depends only on the angle $\theta$:
$$ \mu_l = l(l+1), \quad \Theta_l^m (\theta) = P_l^m ( \cos \theta) \quad \text{for} \: -l \leq m \leq l$$Spherical harmonics¶
We are now ready to find the solutions of the spherical Helmholtz equation:
$$ \nabla^2_S Q(\theta, \phi) + \mu Q(\theta, \phi) = 0 $$That is, the eigenvalues and eigenvectors of the Lagrangian operator on the spherical surface: $-\nabla^2_S Q$. We looked for solutions such as:
$$Q(\theta, \phi) = \Theta(\theta) \Phi(\phi) $$And combining the above results, we obtain the spherical harmonics $Y_l^m(\theta, \phi)$ (degree $l$, order $m$), which are usually multiplied by a normalization constant in order to get functions with norm 1 with respect to the scalar product that we will define later.
$$ \text{Eigenvalues: } \mu_l = l(l+1)$$$$ \text{Eigenfunctions: } Y_l^m(\theta, \phi) = \sqrt{\frac{2l+1}{4 \pi} \frac{(l-m)!}{(l+m)!}} \; \Theta_l^m(\theta) e^{i m \phi} = P_l^m ( \cos \theta) e^{i m \phi} \quad 0 \leq l < \infty, \; -l \leq m \leq l$$Orthogonality of the spherical harmonics¶
This is the most interesting property of the spherical harmonics, that is going to allow us to perform with them developments in series similar to the Fourier series. In order to be able to speak of orthogonality we need to define an inner product in the space $\mathsf L^2(\mathbb S^2)$ of functions, in general with complex values, that are square integrable on the unit sphere. The definition will use the solid angle differential element in spherical coordinates (or if you prefer, the differential element of area on a spherical surface of unit radius), which is given by $d \Omega = \sin (\theta) \; d \theta \; d \phi$. So, the inner product of two functions in $\mathsf L^2(\mathbb S^2)$ is defined as:
$$ \langle f,g \rangle = \int_{\phi=0}^{2 \pi} \int_{\theta=0}^\pi f(\theta, \phi) \overline{g(\theta, \phi)} \sin (\theta) \; d \theta \; d \phi $$In other words, it is the integral of the product of the first function by the complex conjugate of the second, extended to all the surface of the sphere.
Being eigenfunctions of the spherical laplacian differential operator, which is self-adjoint, the spherical harmonics corresponding to distinct eigenvalues (ie, different values of the degree $l$, are orthogonal, ie, their inner product is zero. But the orthogonality property goes beyond that, in the sense that spherical harmonics with the same degree $l$ but different order $n$ are also orthogonal. Furthermore, due to the normalization constant that we have added in their definition, we can talk about orthonormality:
$$ \langle Y_l^m, Y_{l'}^{m'} \rangle = \delta_{l,l'} \: \delta_{m, m'}$$Spherical multipole expansion¶
The spherical harmonics system constitutes a complete orthonormal set in $\mathsf L^2(\mathbb S^2)$, so any scalar square integrable function on the sphere $f(\theta, \phi) : \mathbb S^2 \longrightarrow \mathbb C $ will admit a series expansion of the form::
$$f(\theta, \phi) = \sum_{l=0}^\infty \sum_{m=-l}^l a_{l,m} Y^m_ l(\theta, \phi)$$The series converges to the function according to the norm defined in $\mathsf L^2(\mathbb S^2)$ by the inner product. The coefficients can be obtained by inner-multiplying the series by each function $Y_l^m$ and using the fact that they form an orthonormal set:
$$ a_{l,m} = \langle f, Y_l^m \rangle $$These coefficients are generally complex, and the following relationships are met::
$$ Y_l^{-m} = (-1)^m \overline{Y_l^m} \quad \text{para} \: m>0$$$$ a_{l,-m} = (-1)^m \overline{a_{l, m}} \quad \text{para} \: m>0, \; f \; \text{real}$$This type of development is called multipole expansion, and the value of $l$ indicates the degree of the multipole. Specifically, the term $ a_ {0,0} $ corresponds to the monopole, the terms on $a_{1,-1}, a_{1,0}, a_{1,1}$ correspond to the dipole mode, the terms for $l=2$ would be the quadrupole mode, for $l = 3$ we would have the octupole mode, and so on.
The spherical harmonics in Python SciPy¶
The SciPy function that calculates the spherical harmonics is: scipy.special.sph_harm(m, l, phi, theta) . Note that I changed the angles notation used in the documentation of the function to accommodate it to the one that has been used so far in these notes, replacing $n, \theta, \phi$ by $l, \phi, \theta$.
Both phi and theta can be Numpy arrays.
#Example to calculate Y_4^2
l = 4
m = 2
theta, phi = 0.6, 0.75 # Some arbitrary values of angles in radians
Y42 = sp.sph_harm(m, l, phi, theta)
Y42
Let's see how the above value matches the definition, including the normalization constant. On one side the associated Legendre functions of first kind $P_l^m(z)$ can be calculated by scipy.special.lpmv(m, l, z)
z = np.cos(theta)
P42 = sp.lpmv(m,l,z)
P42
Then we calculate the normalization constant:
f = sci.math.factorial
K_norm = np.sqrt((2*l+1)/(4 * np.pi) * f(l-m)/f(l+m))
K_norm
If we now compare the value obtained for Y42 with the one obtained by directly applying the definition, we can see that they match:
K_norm * P42* np.exp(m*phi*1j) == Y42
Orthonormality of the spherical harmonics¶
We will check numerically that the spherical harmonics are an orthonormal family. Let's define first a function that calculates the inner product $\langle f, g \rangle$ of two functions of $\theta, \phi$, taking complex values, and defined over $\mathbb S^2$
def dotprod(f,g):
#Scipy does not directly integrates complex functions.
#You have to break them down into two integrals of the real and imaginary part
integrand_r = lambda theta, phi: np.real(f(theta, phi) * np.conj(g(theta, phi)) * np.sin(theta))
integrand_i = lambda theta, phi: np.imag(f(theta, phi) * np.conj(g(theta, phi)) * np.sin(theta))
rr = sci.integrate.dblquad(integrand_r, 0, 2 * np.pi,lambda theta: 0, lambda theta: np.pi)[0]
ri = sci.integrate.dblquad(integrand_i, 0, 2 * np.pi,lambda theta: 0, lambda theta: np.pi)[0]
if np.allclose(rr,0):
rr = 0
if np.allclose(ri,0):
ri=0
return rr + ri*1j
# We check the orthogonality of the spherical harmonics:
# Si (l,m) =! (l',m') the inner product must be zero
Y = lambda l, m, theta, phi: sp.sph_harm(m, l, phi, theta)
f = lambda theta, phi: Y(4,3,theta, phi)
g = lambda theta, phi: Y(4,2,theta, phi)
dotprod(f,g)
# And, if (l,m) = (l',m') the inner product is one.
f = lambda theta, phi: Y(4,3,theta, phi)
g = lambda theta, phi: Y(4,3,theta, phi)
dotprod(f,g)
Graphical representation of spherical harmonics¶
Now we will build below several graphical representations of the same spherical harmonic, in order to grasp some idea of its geometry. We'll begin with some usual graphic representations in three dimensions. Since these are complex valued functions, we will represent first the module, and then the real part.
l = 4 #degree
m = 2 # order
PHI, THETA = np.mgrid[0:2*np.pi:200j, 0:np.pi:100j] #arrays of angular variables
R = np.abs(sp.sph_harm(m, l, PHI, THETA)) #Array with the absolute values of Ylm
#Now we convert to cartesian coordinates
# for the 3D representation
X = R * np.sin(THETA) * np.cos(PHI)
Y = R * np.sin(THETA) * np.sin(PHI)
Z = R * np.cos(THETA)
We wish to color the surface depending on the value of the radial distance from the origin, that is, of $| Y^m_ l(\theta, \phi)|$. So, the color is going to be defined by the R array, for which an object ScalarMappable with the RGB values will be created following the instructions of this answer in stackoverflow
N = R/R.max() # Normalize R for the plot colors to cover the entire range of colormap.
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(12,10))
im = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=cm.jet(N))
ax.set_title(r'$|Y^2_ 4|$', fontsize=20)
m = cm.ScalarMappable(cmap=cm.jet)
m.set_array(R) # Assign the unnormalized data array to the mappable
#so that the scale corresponds to the values of R
fig.colorbar(m, shrink=0.8);
And then we will represent only the real part of $Y^m_ l$
l = 4 # degree
m = 2 # order
PHI, THETA = np.mgrid[0:2*np.pi:200j, 0:np.pi:100j]
R = sp.sph_harm(m, l, PHI, THETA).real
X = R * np.sin(THETA) * np.cos(PHI)
Y = R * np.sin(THETA) * np.sin(PHI)
Z = R * np.cos(THETA)
#As R has negative values, we'll use an instance of Normalize
#see http://stackoverflow.com/questions/25023075/normalizing-colormap-used-by-facecolors-in-matplotlib
norm = colors.Normalize()
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(14,10))
m = cm.ScalarMappable(cmap=cm.jet)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=cm.jet(norm(R)))
ax.set_title('real$(Y^2_ 4)$', fontsize=20)
m.set_array(R)
fig.colorbar(m, shrink=0.8);
In the figure above, the negative and positive values of $\text{Real}( Y_l^m)$ are shown clearly differentiated, with negative values represented in blues gamut, and positive values in green to red, depending the values on the orientation in space. Nevertheless I prefer to imagine the spherical harmonics as scalar fields (pressure, temperature, etc.) defined on the unit sphere, where they take positive or negative values that can be displayed as peaks or valleys on the surface. The following plot goes in this direction, by just adding the radius of the spherical surface ($R = 1$) to the value of the real part of the spherical harmonic. Otherwise the method of obtaining the graph is the same as before:
l = 4 # degree
m = 2 # order
PHI, THETA = np.mgrid[0:2*np.pi:300j, 0:np.pi:150j]
R = sp.sph_harm(m, l, PHI, THETA).real
s = 1
X = (s*R+1) * np.sin(THETA) * np.cos(PHI)
Y = (s*R+1) * np.sin(THETA) * np.sin(PHI)
Z = (s*R+1) * np.cos(THETA)
norm = colors.Normalize()
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(14,10))
m = cm.ScalarMappable(cmap=cm.jet)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=cm.jet(norm(R)))
ax.set_title('1 + real$(Y^2_ 4)$', fontsize=20)
m.set_array(R)
fig.colorbar(m, shrink=0.8);
However, my preferred plot is a two dimensional representation with Mollweide projection, which has the property of preserving the ratio of the areas. I also believe that it can be an interesting example on how to build a "heatmap" of a scalar field over the sphere by using a Mollweide projection
# Coordinate arrays for the graphical representation
x = np.linspace(-np.pi, np.pi, 100)
y = np.linspace(-np.pi/2, np.pi/2, 50)
X, Y = np.meshgrid(x, y)
# Spherical coordinate arrays derived from x, y
# Necessary conversions to get Mollweide right
phi = x.copy() # physical copy
phi[x < 0] = 2 * np.pi + x[x<0]
theta = np.pi/2 - y
PHI, THETA = np.meshgrid(phi, theta)
l = 4
m = 2
SH_SP = sp.sph_harm(m, l, PHI, THETA).real # Plot just the real part
What follows is a representation of the spherical harmonic as a heatmap of a field with a Mollweide projection. The value $l-m$ gives the number of horizontal "bands" or lines of constant latitude where the spherical harmonic is null. The order $ m $ is the number of zero circles passing through the poles (beware with the count, in the flat Mollweide projection you'll see only half circles). To check this you can see a very nice interactive 3D visualization here.
#This is to enable bold Latex symbols in the matplotlib titles, according to:
#http://stackoverflow.com/questions/14324477/bold-font-weight-for-latex-axes-label-in-matplotlib
matplotlib.rc('text', usetex=True)
matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
xlabels = ['$210^\circ$', '$240^\circ$','$270^\circ$','$300^\circ$','$330^\circ$',
'$0^\circ$', '$30^\circ$', '$60^\circ$', '$90^\circ$','$120^\circ$', '$150^\circ$']
ylabels = ['$165^\circ$', '$150^\circ$', '$135^\circ$', '$120^\circ$',
'$105^\circ$', '$90^\circ$', '$75^\circ$', '$60^\circ$',
'$45^\circ$','$30^\circ$','$15^\circ$']
fig, ax = plt.subplots(subplot_kw=dict(projection='mollweide'), figsize=(10,8))
im = ax.pcolormesh(X, Y , SH_SP)
ax.set_xticklabels(xlabels, fontsize=14)
ax.set_yticklabels(ylabels, fontsize=14)
ax.set_title('real$(Y^2_ 4)$', fontsize=20)
ax.set_xlabel(r'$\boldsymbol \phi$', fontsize=20)
ax.set_ylabel(r'$\boldsymbol{\theta}$', fontsize=20)
ax.grid()
fig.colorbar(im, orientation='horizontal');
Here ends this post, which has already been too long. Until next time!.
