Spherical harmonics are a powerful tool used in mathematics and physics to describe a wide range of phenomena. They can be used to model physical phenomena such as sound waves, heat transfer, and even quantum systems, and they appear whenever you are working with partial differential equations in spherical coordinates. In this article, we will explore what spherical harmonics are, their most important properties, how they work, and why they’re so useful. Finally, we will see how we can use them in Python, both numerically and symbolically. And we will create our own visualizations of these functions so that you can play with them yourselves.

Spherical harmonics are eigenfunctions

The spherical harmonics are defined in terms of the angular coordinates (latitude and longitude) on the sphere and have the property that they are eigenfunctions of the angular part of the Laplace operator. What does that mean?

In spherical coordinates, the Laplace operator when applied to a function 𝑓 reads

where θ is the polar angle and φ is the azimuthal angle. So the Laplacian naturally splits into a radial and an angular part ( L² ).

Spherical harmonics are the eigenfunctions of the angular part of the Laplace operator. This means that when the angular part of the Laplace operator is applied to a spherical harmonic function, the result is a scalar multiple of the original function:

where L² is the angular part of the Laplace operator, Yₗₘ(θ,φ) is the spherical harmonic of degree l and order m, and -l(l+1) is the eigenvalue of the function.

There are different conventions about the ordering of the indices. Sometimes they are 𝑙𝑚, sometimes reversed, sometimes both indices downstairs, sometimes one upstairs one downstairs. The notation used here agrees with the common notation used in physics and in particular quantum mechanics. That’s also where the notation L² for the angular part of the Laplacian comes from. In quantum mechanics, L² is the square of the angular momentum operator and this gives an intuitive meaning to the indices 𝑙 and 𝑚. 𝑙 (or actually 𝑙(𝑙+1)) stands for the magnitude of the angular momentum of a particle and 𝑚 stands for the z-component of the angular momentum. The important thing to remember is the following: 𝑙 can be any non-negative integer 0, 1, 2, … etc., whereas 𝑚 can only have the values −𝑙,−𝑙+1,−𝑙+1,…,𝑙−2,𝑙−1,𝑙 for given 𝑙.

Spherical harmonics are orthogonal

The orthogonality property of spherical harmonics states that the integral of the product of two spherical harmonics over the surface of the sphere is equal to zero if the degree and order of the two harmonics are not equal. This is the same property that sines and cosines have, too. Think of Fourier transforms!

Formally, the orthogonality property of spherical harmonics can be expressed as:

where δₘₙ is the Kronecker delta, which is 1 if m=n and 0 otherwise, and 𝑑Ω=sin𝜃𝑑𝜃𝑑𝜑 is the integration measure for integrating over the surface of a sphere. Note that many of the spherical harmonics are complex, that’s why there is a conjugation star in the formula.

This property means that the spherical harmonics of different degree and order are independent of each other. This independence allows us to represent any function defined on the surface of the sphere as a sum of spherical harmonics of different degree and order, without any cross-term interference.

This orthogonality property also allows to simplify the calculation of the integral of a product of two functions defined on the sphere, by decomposing each function into a series of spherical harmonics and then using the orthogonality property to calculate the integral term by term.

Spherical harmonics are normalized

Spherical harmonics are also normalized. That means that the integral of the square of a spherical harmonic over the surface of the sphere is equal to one:

This property ensures that the spherical harmonics are properly scaled. This also means that the sum of the squares of the coefficients of the spherical harmonics expansion of a function is equal to one. This normalization property also allows to simplify the calculation of the integral of the square of a function defined on the sphere, by decomposing the function into a series of spherical harmonics and then using the normalization property to calculate the integral term by term.

Spherical harmonics form a basis

The completeness property of spherical harmonics means that any function defined on the surface of the sphere can be represented as a sum of spherical harmonics of different degrees and orders. This property is also known as the “spherical harmonics expansion”. This can be expressed as:

Where f(θ,φ) is any function defined on the surface of the sphere and fₗₘ are the expansion coefficients. The expansion coefficients can be calculated by integrating the product of the function and each spherical harmonics over the surface of the sphere using the orthogonality property:

This is absolutely analogous to how you would compute a Fourier transform. But in this case, you don’t expand in terms of sines and cosines but in terms of spherical harmonics.

The completeness property allows to represent any function defined on the surface of the sphere in terms of a sum of simpler, computationally manageable functions. It’s important to note that the completeness of the set holds under certain conditions, such as the function being square integrable.

Spherical harmonics are symmetric

The parity property of spherical harmonics states that the harmonics with even degree have even parity and the harmonics with an odd degree have odd parity. Formally, the parity property of spherical harmonics can be expressed as:

and also

This means that the spherical harmonics of even degree are unchanged when the polar angle is reflected about the equator (i.e. when θ is replaced by π-θ ), whereas the spherical harmonics of odd degree change sign.

How to use spherical harmonics in Python

When you are using spherical harmonics in any software library, you have to watch out carefully, for how they are defined. Which index has the meaning of 𝑙, and which corresponds to 𝑚? What is the order of the angles, is it (𝜃,𝜑) or (𝜑,𝜃)? In the world of scientific Python, there are mainly two packages which are commonly used for spherical harmonics. The one is SymPy, which is a computer algebra package that allows for symbolic computation. The other package is SciPy, which is for numerical computation.

SymPy

Let’s start with SymPy. In a Jupyter notebook, first import the symbols that we need:

from sympy import *

θ, φ = symbols(r"\theta, \phi")
l, m = symbols("l, m", integer=True)

The spherical harmonics can be accessed by the function Ynm . In the list of arguments, 𝑙 comes first and 𝑚 second. The order of the angles is (𝜃,𝜑):

Y = Ynm(l, m, θ, φ)

The spherical harmonics can be expressed in terms of trigonometric functions. We can easily list them using SymPy, in this case, up to 𝑙=3:

for l_ in range(3):
    for m_ in range(0, l_+1):
        Y = Ynm(l_, m_, θ, φ)
        display(
            Eq(
                Y, Y.rewrite(cos)
            )
        )

Note that we didn’t list the spherical harmonics for negative 𝑚 because they are closely related to the corresponding function with positive 𝑚 (see the section on “parity”).

SciPy

When it comes to numeric computation in Python, SciPy is an excellent choice. The spherical harmonics are part of the special subpackage for special functions. Be careful, because, in SciPy, the arguments are given in a different order than in SymPy. 𝑚 comes first and l second, and also the angle are given in reversed order. When you only want to compute the function for a single point, you pass single numbers for 𝜑 and 𝜃:

from scipy.special import sph_harm

l_ = 1
m_ = 1

theta = np.pi/2
phi = 0

# Watch out for order of arguments, it's different than in SymPy!
sph_harm(m_, l_, phi, theta)

But usually, in numeric computation, you want to give arrays of angular values, which can easily be done in NumPy:

import numpy as np

# Define the 1-dimensional coordinate axes:
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)

# Make a 2D-meshed grid out of the axes:
Theta, Phi = np.meshgrid(theta, phi, indexing="ij")


sph_harm(m_, l_, Phi, Theta)

Here is a little function that creates a plot of any desired spherical harmonic using the Python package plotly:

import plotly.graph_objects as go


def create_figure(l, m):
    """Create a figure of Y_lm using Plotly."""

    thetas = np.linspace(0, np.pi, 100)
    phis = np.linspace(0, 2*np.pi, 100)
    
    Thetas, Phis = np.meshgrid(thetas, phis)
    
    ylm = sph_harm(m, l, Phi, Theta)

    fcolors = ylm.real
    fmax, fmin = fcolors.max(), fcolors.min()
    fcolors = (fcolors - fmin)/(fmax - fmin)

    R = abs(ylm)
    X = R * np.sin(Theta) * np.cos(Phi)
    Y = R * np.sin(Theta) * np.sin(Phi)
    Z = R * np.cos(Theta)

    
    fig = go.Figure(
                data=[
                    go.Surface(
                        x=X, y=Y, z=Z, 
                        surfacecolor=fcolors,
                        colorscale='balance', showscale=False, 
                        opacity=1.0, hoverinfo='none',
                    )
                ]
    )

    fig.update_layout(title='$Y_{%d%d}$' % (l, m), autosize=False,
                      width=700, height=700,
                      margin=dict(l=65, r=50, b=65, t=90)
    )
    
    return fig

Now we can create figures in Jupyter notebooks by calling that function:

create_figure(1, 0)

create_figure(3, 1)

Conclusion

We have explored the most important properties of the spherical harmonics. Like for sines and cosines, you should know a handful of properties like addition theorems that are sufficient as working knowledge. Of course, like for trigonometric functions, there are a lot of additional properties that are needed in rare cases and can be looked up if needed. But the handful you know by heart usually keeps you going. It’s the same with spherical harmonics. If you know the properties we discussed above, this should cover 95% of all you need in computations.