Visualizing the real forms of the spherical harmonics

The spherical harmonics are a set of special functions defined on the surface of a sphere that originate in the solution to Laplace's equation, $\nabla^2f=0$. Because they are basis functions for irreducible representations of SO(3), the group of rotations in three dimensions, they appear in many scientific domains, in particular as the angular part of the wavefunctions of atoms (atomic orbitals). They are described in terms of an integer degree $l=0,1,2,\ldots$ and order $m=-l,-l+1,\ldots,l$.

In quantum mechanics, $l$ is identified with the (orbital) angular momentum quantum number, and $m$ with the azimuthal quantum number. In this domain, they are usually defined including a factor of $(-1)^m$ (the Condon–Shortley phase convention): $$ Y_l^m( \theta , \varphi ) = (-1)^m \sqrt{{(2l+1)\over 4\pi}{(l-m)!\over (l+m)!}} \, P_{l m} ( \cos{\theta} ) \, e^{i m \varphi } $$ where $P_{l m} ( \cos{\theta} )$ is an associated Legendre polynomial (without the factor of $(-1)^m$.)

Since $Y_l^m( \theta , \varphi )$ are complex functions of angle, it is often considered more convenient to use their real forms for their depiction in figures and in some calculations. A suitable real basis of spherical harmonics may be defined as:

$$ Y_{lm} = \left\{ \begin{array}{ll} \displaystyle \sqrt{2} \, (-1)^m \, \operatorname{Im}[{Y_l^{|m|}}] & \text{if}\ m<0\\ \displaystyle Y_l^0 & \text{if}\ m=0\\ \displaystyle \sqrt{2} \, (-1)^m \, \operatorname{Re}[{Y_l^m}] & \text{if}\ m>0. \end{array} \right. $$

The code below uses SciPy's special.sph_harm routine to calculate the spherical harmonics, which are then cast into these real functions and displayed in a three-dimensional Matplotlib plot.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
# The following import configures Matplotlib for 3D plotting.
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import sph_harm
plt.rc('text', usetex=True)

# Grids of polar and azimuthal angles
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
# Create a 2-D meshgrid of (theta, phi) angles.
theta, phi = np.meshgrid(theta, phi)
# Calculate the Cartesian coordinates of each point in the mesh.
xyz = np.array([np.sin(theta) * np.sin(phi),
                np.sin(theta) * np.cos(phi),
                np.cos(theta)])

def plot_Y(ax, el, m):
    """Plot the spherical harmonic of degree el and order m on Axes ax."""

    # NB In SciPy's sph_harm function the azimuthal coordinate, theta,
    # comes before the polar coordinate, phi.
    Y = sph_harm(abs(m), el, phi, theta)

    # Linear combination of Y_l,m and Y_l,-m to create the real form.
    if m < 0:
        Y = np.sqrt(2) * (-1)**m * Y.imag
    elif m > 0:
        Y = np.sqrt(2) * (-1)**m * Y.real
    Yx, Yy, Yz = np.abs(Y) * xyz

    # Colour the plotted surface according to the sign of Y.
    cmap = plt.cm.ScalarMappable(cmap=plt.get_cmap('PRGn'))
    cmap.set_clim(-0.5, 0.5)

    ax.plot_surface(Yx, Yy, Yz,
                    facecolors=cmap.to_rgba(Y.real),
                    rstride=2, cstride=2)

    # Draw a set of x, y, z axes for reference.
    ax_lim = 0.5
    ax.plot([-ax_lim, ax_lim], [0,0], [0,0], c='0.5', lw=1, zorder=10)
    ax.plot([0,0], [-ax_lim, ax_lim], [0,0], c='0.5', lw=1, zorder=10)
    ax.plot([0,0], [0,0], [-ax_lim, ax_lim], c='0.5', lw=1, zorder=10)
    # Set the Axes limits and title, turn off the Axes frame.
    ax.set_title(r'$Y_{{{},{}}}$'.format(el, m))
    ax_lim = 0.5
    ax.set_xlim(-ax_lim, ax_lim)
    ax.set_ylim(-ax_lim, ax_lim)
    ax.set_zlim(-ax_lim, ax_lim)
    ax.axis('off')

fig = plt.figure(figsize=plt.figaspect(1.))
ax = fig.add_subplot(projection='3d')
l, m = 3, 0
plot_Y(ax, l, m)
plt.savefig('Y{}_{}.png'.format(l, m))
plt.show()

el_max = 3
figsize_px, DPI = 800, 100
figsize_in = figsize_px / DPI
fig = plt.figure(figsize=(figsize_in, figsize_in), dpi=DPI)
spec = gridspec.GridSpec(ncols=2*el_max+1, nrows=el_max+1, figure=fig)
for el in range(el_max+1):
    for m_el in range(-el, el+1):
        print(el, m_el)
        ax = fig.add_subplot(spec[el, m_el+el_max], projection='3d')
        plot_Y(ax, el, m_el)
plt.tight_layout()
plt.savefig('sph_harm.png')
plt.show()

Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

Arend Lammertink 3 years, 6 months ago

Seems like there is an error in the first script, I get a None object error.

Should the call to fig.add be:

ax = fig.add_subplot(111, projection='3d')

?

Link | Reply

christian 3 years, 6 months ago

Hmm – this works for me. Which version of Matplotlib are you using? They may have made the subplot position argument optional (ie the default, 111) in the version I use.

lewis 2 years, 10 months ago

I've consulted this many times - thanks a lot!

Lai Hoigin 2 years, 4 months ago

I have been searching for a method to visualize the spherical harmonic functions with python for a long time! Thanks a lot!

Mike Reid 2 years, 3 months ago

Thanks for this. Your blog and book are proving very useful to my migration to python, and I appreciate the physics/chemistry -related examples.

A small convention comment is that in this code phi=0 is along the y axis, whereas the usual convention in quantum mechanics is for phi=0 to be along x.

christian 2 years, 3 months ago

I'm glad you find it useful! You're probably right about the phi convention – I'm not sure it's too important for visualization, though, since I think one can just rotate 90 deg anticlockwise around the z-axis (or swap, m = -m)?