Code
import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import numpy as np
October 10, 2021
Spherical harmonics are powerful mathematical tools, allowing us to represent any function on a sphere as the sum of simpler basis functions (much like a Fourier series!). This post aims to explain spherical harmonics through the lens of Fourier series, by “lifting” them from a circle to a sphere, extensively relying on visual explanations.
Suppose we have a periodic function
Since the function is periodic, there is an equivalent way to view this: as a function defined on a circle with radius 1 (and in turn, a circumference of 2π). The period here is represented counter-clockwise (as is convention in trigonometry). Just as in the previous plot, we can see that the function first rises to a large “hill”, before descending into two small “valleys”, returning to the starting point.
One of the most widely used discoveries in mathematics is that any real-valued, periodic function1 can be represented as the weighted sum of sines and cosines2 of varying frequencies, called a Fourier Series. Specifically, for a function with a period of
The process of finding the exact weights belong to the study of Fourier Analysis. In our case, since the function we’ve chosen is already neatly written as a sum of sines and cosines,
Note that for a finite number of terms, this is exact only when the function can be expressed in terms of sines and cosines; otherwise (with function such as say, a square wave) this only exact when the sum uses an infinite number of terms.
Looping back5 to earlier, these basis functions can also be represented on a circle:
The takeaway here is that even though the circle plots look very different from the standard plot (with
If you’ve followed the visuals so far, you already know what a spherical harmonic is: instead of basis functions defined on a circle, they’re basis functions defined on a sphere 6. Just like how
fig = make_subplots(rows=1, cols=1, specs=[[{'is_3d': True}]])
s = np.linspace(0, 2 * np.pi, 100)
t = np.linspace(0, np.pi, 100)
tGrid, sGrid = np.meshgrid(s, t)
r = 1
x = r * np.cos(sGrid) * np.sin(tGrid)
y = r * np.sin(sGrid) * np.sin(tGrid)
z = r * np.cos(tGrid)
fig.add_trace(
go.Surface(x=x, y=y, z=z,
surfacecolor=np.sqrt(3/4*np.pi)*y,
colorscale='RdBu',
colorbar=dict(thickness=10)
)
)
fig.update_layout(
font_family="JuliaMono",
showlegend=False,
margin=go.layout.Margin(l=0, r=0, b=0, t=0),
paper_bgcolor='rgba(0,0,0,0)',
)
fig.update_traces(showscale=False)
fig.show()
On the circle, Fourier Analysis allows us to convert any function defined on a circle (i.e. a periodic function) into the weighted sum of sines and cosines of varying frequencies. Now, on the sphere, any function on a sphere can be decomposed into a weighted sum of the spherical harmonic functions
There’s a few things to note here:
Real vs Complex: Both indices are subscripts on
Coordinates: Note that we’re using cartesian coordinates (x, y, z) as inputs to the spherical harmonics, not angles. We could equivalently use harmonics of the form
Often, spherical harmonics are not depicted on the sphere directly, but in a different form. For instance, in the alternate depiction
fig = make_subplots(rows=1, cols=1,
specs=[[{'is_3d': True}]])
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
theta, phi = np.meshgrid(theta, phi)
x = np.sin(theta) * np.cos(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(theta)
r0a = np.sqrt(3/4*np.pi)*y
r0 = np.abs(r0a)
fig.add_trace(
go.Surface(
x=r0*x, y=r0*y, z=r0*z,
surfacecolor=r0a,
colorscale='RdBu'
)
)
fig.update_layout(
font_family="JuliaMono",
showlegend=False,
margin=go.layout.Margin(l=0, r=0, b=0, t=0),
paper_bgcolor='rgba(0,0,0,0)',
)
fig.update_layout(
scene = dict(
xaxis = dict(nticks=4, range=[-1.8,1.8],),
yaxis = dict(nticks=4, range=[-1.8,1.8],),
zaxis = dict(nticks=4, range=[-1.8,1.8],)
),
scene1_aspectratio=dict(x=1, y=1, z=1)
)
fig.update_traces(showscale=False)
fig.show()
Note that while it appears that
Now, we’re plotting a surface through all points:
Specifically, we allow the magnitude of the function
As we can see, on the plane
To really visually grasp this, it’s helpful to drop back down to the 2D case. Let’s do the shrink and expand procedure we just described to
Completing the visual analogy, here’s the original and alternate version diagrams for
The alternate visual depiction of the harmonics may seem needlessly complicated at first. But they serve an important purpose: as the harmonics become more complex as the degree goes up and we begin linearly combining them together, the alternate form allows a clearer depiction of how different values of the function on the sphere compare to each other. For example, here’s
fig = make_subplots(rows=1, cols=1,
specs=[[{'is_3d': True}]])
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
theta, phi = np.meshgrid(theta, phi)
x = np.sin(theta) * np.cos(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(theta)
r0a = 3/16*np.sqrt(1/np.pi)*(35*z**4-30*z**2+3) + 2*(1/2*np.sqrt(15/np.pi)*y*z)
r0 = np.abs(r0a)
fig.add_trace(
go.Surface(
x=r0*x, y=r0*y, z=r0*z,
surfacecolor=r0a,
colorscale='RdBu'
)
)
fig.update_layout(
font_family="JuliaMono",
showlegend=False,
margin=go.layout.Margin(l=0, r=0, b=0, t=0),
paper_bgcolor='rgba(0,0,0,0)',
)
fig.update_layout(
scene = dict(
xaxis = dict(nticks=4, range=[-1.8,1.8],),
yaxis = dict(nticks=4, range=[-1.8,1.8],),
zaxis = dict(nticks=4, range=[-1.8,1.8],)
),
scene1_aspectratio=dict(x=1, y=1, z=1)
)
camera = dict(eye=dict(x=1.5, y=1.5, z=1.5))
fig.layout.scene1.camera = camera
fig.update_traces(showscale=False)
fig.show()
You may have noticed I’ve been referring to the harmonics by their degree and order (such as
These are the Cartesian versions of the spherical harmonics, but it is important to know that conventions vary significantly across fields (for instance, the quantum mechanics community implictly places a factor of
To get a version of the harmonics in spherical coordinates (to get something akin to terms in the Fourier series), all one has to do is apply the standard change of basis formula9:
As an example,
fig = make_subplots(rows=1, cols=1,
specs=[[{'is_3d': True}]])
theta = np.linspace(0, 2*np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
theta, phi = np.meshgrid(theta, phi)
x = np.sin(theta) * np.cos(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(theta)
r0a = np.sqrt(3/4*np.pi) * np.sin(theta) * np.sin(phi)
# r0 = np.abs(r0a)
fig.add_trace(
go.Surface(x=theta, y=phi, z=r0a,
surfacecolor=r0a,
colorscale='RdBu'
)
)
fig.update_layout(
font_family="JuliaMono",
showlegend=False,
margin=go.layout.Margin(l=0, r=0, b=0, t=0),
paper_bgcolor='rgba(0,0,0,0)',
)
fig.update_layout(
scene = dict(
xaxis = dict(
tickvals=[0, np.pi/4, np.pi/2, 3*np.pi/4, np.pi],
ticktext=['0', 'pi/4', 'pi/2', '3pi/4', 'pi'],
range=[0,np.pi],
title_text="θ"),
yaxis = dict(
tickvals=[0, np.pi/4, np.pi/2, 3*np.pi/4, np.pi, 5*np.pi/4, 6*np.pi/4, 7*np.pi/4, 8*np.pi/4],
ticktext=['0', 'pi/4', 'pi/2', '3pi/4', 'pi', '5pi/4', '3/2pi', '7pi/4', '2pi'],
nticks=4,
range=[0,2*np.pi],
title_text="ϕ"),
zaxis = dict(
title_text="Y(θ, ϕ)",
nticks=3,
tickvals=[-1.5, 0, 1.5])),
scene1_aspectratio=dict(x=1, y=2, z=1),
margin=dict(r=20, l=10, b=10, t=10)
)
camera = dict(eye=dict(x=2, y=2, z=2))
fig.layout.scene1.camera = camera
fig.update_traces(showscale=False)
fig.show()
Just as how we were able to interpret a periodic function on the 1D line as being on a circle, we can now see certain periodic functions on the 2D plane as existing on a sphere11. This now also means we have three different ways of looking at periodic functions on circles and spheres!
Firstly, in the 1D case, we can represent
Similarly, in the 2D case, we can represent
Now that we’ve finished making the connection between the Fourier series and Spherical harmonics, one might ask: what are they actually used for? There’s plenty (they pop up rather frequently when functions on a sphere are involved), here are some:
Solving PDEs: Recall that the motivating reason behind the discovery of the Fourier series was to solve the heat equation, which is a partial differential equation. Likewise, we can now decompose any (well-behaved) function on a sphere into a sum of spherical harmonics as basis functions, and solve the PDE using these basis functions (instead of the more complicated original)
Physics: One particularly important PDE is Schrödinger’s equation. Focusing on hydrogen-like atoms, the spherical harmonics can be use to describe the angular component of the solutions, which can then be converted into the probability that the electron will be found in a specific direction. Note that the angular component cannot tell us the probability how far from the origin the electron will be found (that’d be the radial component), only the probability in a given direction.
Machine Learning: One problem with ordinary neural networks is that they do not respect input geometry. For example, if your input is a molecule (with atom coordinates) and your output is toxicity, if you rotate the input, the hidden layer outputs will change in an ill-defined way and return a different prediction, even though your input is the same molecule! Spherical harmonics are used to build rotationally equivariant neural networks, such that the change for a rotation properly is defined, and then ultimately produce an invariant prediction: predicting the same toxicity level for a molecule, regardless of input orientation.
Spherical harmonics were one of those things that felt intimidating to me at first, because of the dense math involved. But ultimately, they’re basically a sibling of the Fourier series, the only difference being that they’re for functions defined on a sphere instead of a circle. The math, once visualized, is beautiful, and I hope this has been a helpful read!
Appendix B, Spherical Harmonics: This chapter from Wojciech Jarosz’s Ph.D. dissertation is a great resource on spherical harmonics, particularly the real-valued ones (it’s written from a computer graphics standpoint; most physics-oriented ones would focus on the complex-valued ones).
e3nn library: This PyTorch library is geared toward equivariant neural networks broadly, but also contains optimized implementations to compute spherical harmonics.
SEGNNs paper: This is a recent paper on constructing equivariant neural networks; the appendix has a great section on how spherical harmonics are used to build them.
Or more broadly, even nonperiodic and/or complex-valued functions. We’re not taking a look at them here, but for a deeper dive there’s this great 3b1b video.↩︎
There are many different conventions on how to represent fourier series (such as the amplitude-phase form, the exponential form, etc.), which you can find here. We’re using the sine-cosine form for this article.↩︎
The curious reader could verify this directly! Each coefficient is defined as the integral over the period (here,
They’re called basis functions, as a linear combination of these functions can produce any periodic function; akin to how a linear combination of a basis in
More formally, these infinitely many sine and cosine functions together form a basis for “square-integrable” functions defined on the circle, that is, for
Terrible pun absolutely intended.↩︎
More precisely, they form a basis for
This also means in practice, efficient calculation of the values are non-trivial, since there are so many different functions to account for. The functions themselves are not arbitrary; they’re made of the derivatives of Legendre polynomials of varying degrees.↩︎
A more exhaustive table of them can be found on Wikipedia here. Note that we’re assuming r=1 throughout this article.↩︎
For a recap of Spherical coordinates, this may be helpful. Note that there’s multiple conventions! We’re using the physics convention here (in the math one,
Note that we plot
Note for this to hold, the values on the entire boundary must be the same (not just the endpoint as in the 1D case), as they get mapped to one of the two poles on the sphere; here all points on the boundary have the value 0.
For a more general periodic function (where the entire boundary doesn’t have the same value), where the period endpoints of any 1D slice (parallel to the x or y axes) of the 2D function are the same, it is possible to see that function existing on a torus, that is,