Visualizing the real forms of the spherical harmonics

(2 comments)

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.

Y(3,0) real spherical harmonic

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()

To plot a family of these functions, try:

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()

Spherical harmonics for l = 0, 1, 2, 3

Current rating: 4.6

Comments

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

Arend Lammertink 1 month, 2 weeks 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
Currently unrated

christian 1 month, 1 week 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.

Link | Reply
Currently unrated

New Comment

required

required (not published)

optional

required