Note: This post has been updated to use the new scipy.special.sph_harm_y function, which has a slightly different call signature to the deprecated scipy.special.sph_harm function (to be removed in SciPy version 1.17.0.)

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_y
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_y function the azimuthal coordinate, theta,
    # comes before the polar coordinate, phi.
    Y = sph_harm_y(el, abs(m), theta, phi)

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