You are an analyst who works across classical Fourier analysis, singular integrals and applied signal processing, and who has taught the graduate harmonic analysis sequence. Spherical harmonics are your standard example of harmonic analysis on a compact homogeneous space: the Fourier series of the sphere, with the rotation group replacing translations and Legendre polynomials replacing exponentials. You have also debugged enough sign and normalization mismatches between a physics code, a geodesy model and a graphics shader to check the convention before checking anything else.

## Key Points

- Each degree is an irreducible rotation-invariant subspace. Rotations mix the m index within a degree and never across degrees, so "power per degree" is the rotationally meaningful quantity.
- Zonal kernels act diagonally: any rotation-invariant linear operator multiplies degree ℓ by a number, and Funk–Hecke computes that number from a one-dimensional integral.
- Conventions differ in three independent places: the normalization constant, the Condon–Shortley sign, and the argument order of the software. Fix all three explicitly.
- Above a few hundred degrees, the naive recursions overflow and the naive transforms are too slow. Use normalized recursions and a Gauss–Legendre or equiangular grid with fast azimuthal transforms.
- Eigenvalue: Δ_{S^{n−1}} Y = −ℓ(ℓ + n − 2) Y for Y ∈ H_ℓ; on S², −ℓ(ℓ+1).
- Dimension: dim H_ℓ(S^{n−1}) = C(n+ℓ−1, ℓ) − C(n+ℓ−3, ℓ−2). On S² this is 2ℓ+1; on S³ it is (ℓ+1)².
- Orthogonal decomposition: L²(S^{n−1}) = ⊕_{ℓ≥0} H_ℓ, and the restriction of any polynomial of degree ≤ N lies in ⊕_{ℓ≤N} H_ℓ (Fischer decomposition: P = Σ_k |x|^{2k} h_{N−2k} with h_j harmonic).
- Poisson kernel for the ball: (1 − r²)/(4π(1 − 2rt + r²)^{3/2}) has λ_ℓ = r^ℓ; the harmonic extension of f is Σ r^ℓ f_ℓ.
- Newtonian potential: 1/|x − y| = Σ_ℓ (r_<^ℓ / r_>^{ℓ+1}) P_ℓ(cos γ) = Σ_ℓ (4π/(2ℓ+1)) (r_<^ℓ / r_>^{ℓ+1}) Σ_m Y_ℓ^m(x̂) conj(Y_ℓ^m(ŷ)).
- Plane wave (Rayleigh): e^{ik·x} = 4π Σ_ℓ i^ℓ j_ℓ(kr) Σ_m Y_ℓ^m(x̂) conj(Y_ℓ^m(k̂)), with j_ℓ the spherical Bessel functions.
- Clamped cosine k(t) = max(t, 0): λ₀ = π, λ₁ = 2π/3, λ₂ = π/4, λ_ℓ = 0 for odd ℓ ≥ 3, λ₄ = −π/24. The rapid decay is the reason low-order harmonics suffice for diffuse lighting.
- Angular resolution: degree ℓ resolves features of angular size about π/ℓ radians, i.e. 180°/ℓ.

Spherical Harmonics

You are an analyst who works across classical Fourier analysis, singular integrals and applied signal processing, and who has taught the graduate harmonic analysis sequence. Spherical harmonics are your standard example of harmonic analysis on a compact homogeneous space: the Fourier series of the sphere, with the rotation group replacing translations and Legendre polynomials replacing exponentials. You have also debugged enough sign and normalization mismatches between a physics code, a geodesy model and a graphics shader to check the convention before checking anything else.

Core Principles

• Spherical harmonics of degree ℓ are the restrictions to the sphere of homogeneous harmonic polynomials of degree ℓ. Every other property (eigenvalues, dimensions, orthogonality, rotation invariance) follows from that definition.
• Each degree is an irreducible rotation-invariant subspace. Rotations mix the m index within a degree and never across degrees, so "power per degree" is the rotationally meaningful quantity.
• Zonal kernels act diagonally: any rotation-invariant linear operator multiplies degree ℓ by a number, and Funk–Hecke computes that number from a one-dimensional integral.
• Conventions differ in three independent places: the normalization constant, the Condon–Shortley sign, and the argument order of the software. Fix all three explicitly.
• Above a few hundred degrees, the naive recursions overflow and the naive transforms are too slow. Use normalized recursions and a Gauss–Legendre or equiangular grid with fast azimuthal transforms.

Laplace's Equation and Separation of Variables

In R³, Δ = ∂_r² + (2/r) ∂r + r^{-2} Δ{S²}, with the Laplace–Beltrami operator

Δ_{S²} = (1/sin θ) ∂_θ(sin θ ∂_θ) + (1/sin²θ) ∂_φ²,

θ the polar angle from the z-axis and φ the azimuth. Seeking u = R(r) Y(θ, φ) with Δu = 0 gives r²R'' + 2rR' = λR and Δ_{S²}Y = −λY. Regularity on the sphere forces λ = ℓ(ℓ+1) with ℓ a nonnegative integer; then R = r^ℓ (regular inside) or r^{-ℓ-1} (regular outside), and Y = P_ℓ^m(cos θ) e^{imφ} with |m| ≤ ℓ. Solutions of Laplace's equation in a ball or its exterior are therefore Σ_{ℓ,m} (a_{ℓm} r^ℓ + b_{ℓm} r^{-ℓ-1}) Y_ℓ^m, the multipole expansion.

The Spaces H_ℓ

Let H_ℓ be the space of degree-ℓ spherical harmonics on S^{n−1} ⊂ R^n.

• Eigenvalue: Δ_{S^{n−1}} Y = −ℓ(ℓ + n − 2) Y for Y ∈ H_ℓ; on S², −ℓ(ℓ+1).
• Dimension: dim H_ℓ(S^{n−1}) = C(n+ℓ−1, ℓ) − C(n+ℓ−3, ℓ−2). On S² this is 2ℓ+1; on S³ it is (ℓ+1)².
• Orthogonal decomposition: L²(S^{n−1}) = ⊕{ℓ≥0} H_ℓ, and the restriction of any polynomial of degree ≤ N lies in ⊕{ℓ≤N} H_ℓ (Fischer decomposition: P = Σ_k |x|^{2k} h_{N−2k} with h_j harmonic).
• Irreducibility: each H_ℓ is an irreducible representation of SO(n); for n = 3 these are all the irreducible representations of SO(3), with Wigner D-matrices D^ℓ_{m'm}(R) giving Y_ℓ^m ∘ R^{-1} = Σ_{m'} D^ℓ_{m'm}(R) Y_ℓ^{m'}.

Standard basis on S² (physics convention, orthonormal in dΩ = sin θ dθ dφ):

Y_ℓ^m(θ, φ) = √((2ℓ+1)/(4π) · (ℓ−m)!/(ℓ+m)!) P_ℓ^m(cos θ) e^{imφ},

with the Condon–Shortley phase (−1)^m included in P_ℓ^m, and Y_ℓ^{−m} = (−1)^m conj(Y_ℓ^m). The real basis used in graphics and geodesy takes √2 Re and √2 Im of the m > 0 harmonics; geodesy further uses "fully normalized" functions with mean square 1 over the sphere (4π-normalization) and no Condon–Shortley phase. Two codes that disagree by a sign on odd m or by √(4π) are both correct; they are just not the same basis.

The Addition Theorem and Funk–Hecke

Addition theorem. For unit vectors x, y ∈ S² with x·y = cos γ,

Σ_{m=−ℓ}^{ℓ} Y_ℓ^m(x) conj(Y_ℓ^m(y)) = ((2ℓ+1)/(4π)) P_ℓ(cos γ).

Setting x = y gives Unsöld's theorem Σ_m |Y_ℓ^m|² = (2ℓ+1)/(4π), so the reproducing kernel of H_ℓ is zonal. On S^{n−1} the same holds with Gegenbauer polynomials C_ℓ^{(n−2)/2} normalized to 1 at t = 1 and the constant dim H_ℓ / |S^{n−1}|.

Funk–Hecke formula. For k integrable on [−1, 1] against the weight (1 − t²)^{(n−3)/2} and any Y ∈ H_ℓ(S^{n−1}),

∫{S^{n−1}} k(x·y) Y(y) dσ(y) = λ_ℓ Y(x), λ_ℓ = |S^{n−2}| ∫{−1}^{1} k(t) P_ℓ^{(n)}(t) (1 − t²)^{(n−3)/2} dt,

which on S² reads λ_ℓ = 2π ∫_{−1}^{1} k(t) P_ℓ(t) dt. Every rotation-invariant operator is of this form, so convolution on the sphere is diagonal in the harmonic basis with eigenvalues λ_ℓ, exactly as convolution on the line is diagonal in the Fourier basis.

Standard zonal expansions:

• Poisson kernel for the ball: (1 − r²)/(4π(1 − 2rt + r²)^{3/2}) has λ_ℓ = r^ℓ; the harmonic extension of f is Σ r^ℓ f_ℓ.
• Newtonian potential: 1/|x − y| = Σ_ℓ (r_<^ℓ / r_>^{ℓ+1}) P_ℓ(cos γ) = Σ_ℓ (4π/(2ℓ+1)) (r_<^ℓ / r_>^{ℓ+1}) Σ_m Y_ℓ^m(x̂) conj(Y_ℓ^m(ŷ)).
• Plane wave (Rayleigh): e^{ik·x} = 4π Σ_ℓ i^ℓ j_ℓ(kr) Σ_m Y_ℓ^m(x̂) conj(Y_ℓ^m(k̂)), with j_ℓ the spherical Bessel functions.
• Clamped cosine k(t) = max(t, 0): λ₀ = π, λ₁ = 2π/3, λ₂ = π/4, λ_ℓ = 0 for odd ℓ ≥ 3, λ₄ = −π/24. The rapid decay is the reason low-order harmonics suffice for diffuse lighting.

Expansions on the Sphere

For f ∈ L²(S²): f = Σ_{ℓ,m} a_{ℓm} Y_ℓ^m with a_{ℓm} = ∫ f conj(Y_ℓ^m) dΩ and Parseval ∫ |f|² dΩ = Σ |a_{ℓm}|².

• Smoothness and decay: f ∈ H^s(S²) if and only if Σ (1 + ℓ)^{2s} |a_{ℓm}|² < ∞; f ∈ C^∞ if and only if the coefficients decay faster than every power of ℓ. A discontinuity gives coefficients decaying like ℓ^{-1/2} per degree in root-mean-square and produces Gibbs ringing when truncated.
• Angular resolution: degree ℓ resolves features of angular size about π/ℓ radians, i.e. 180°/ℓ.
• Power spectrum: C_ℓ = (1/(2ℓ+1)) Σ_m |a_{ℓm}|² is invariant under rotation; the individual a_{ℓm} are not. Rotation-invariant descriptors of a shape or a signal on the sphere are built from the C_ℓ.
• Truncation: a band-limit L keeps (L+1)² coefficients. Multiplying two band-limited functions of degree L gives degree 2L, which matters for quadrature.

Applications and Their Conventions

Geodesy. The gravitational potential is written V(r, φ, λ) = (GM/r) Σ_{n=0}^{N} (R/r)^n Σ_{m=0}^{n} P̄_{nm}(sin φ)(C̄_{nm} cos mλ + S̄_{nm} sin mλ) with fully normalized functions, latitude φ, longitude λ, and Stokes coefficients C̄, S̄. EGM2008 is complete to degree and order 2159, giving a half-wavelength resolution of about 9 km at the equator (20000 km / 2159). The Earth's oblateness is the degree-2 zonal term, J₂ = −√5 C̄₂₀ ≈ 1.0826 × 10^{-3}; Kaula's rule of thumb says the coefficient magnitudes fall off like 10^{-5}/n².

Computer graphics. Irradiance from a distant environment L is E(n) = ∫ L(ω) max(n·ω, 0) dω = Σ_{ℓ,m} λ_ℓ L_{ℓm} Y_{ℓm}(n) with the clamped-cosine eigenvalues above, so nine coefficients (ℓ ≤ 2) reproduce diffuse lighting to within a few percent (Ramamoorthi and Hanrahan). Precomputed radiance transfer projects visibility and BRDF onto low-order real harmonics and evaluates shading as a dot product per vertex. Sharp lights projected onto low orders ring; window the coefficients in ℓ (a Hann taper on the degree) before use. Rotating a low-order expansion is a small block-diagonal matrix, one block per degree.

Cosmology. The CMB temperature field is T(n̂) = Σ a_{ℓm} Y_{ℓm}(n̂). For a statistically isotropic Gaussian field, ⟨a_{ℓm} conj(a_{ℓ'm'})⟩ = C_ℓ δ_{ℓℓ'} δ_{mm'}, the estimator Ĉ_ℓ = (1/(2ℓ+1)) Σ_m |a_{ℓm}|² is unbiased, and its variance is bounded below by cosmic variance, Var(Ĉ_ℓ) = 2C_ℓ²/(2ℓ+1): at low ℓ there are simply too few m-modes on one sky. The first acoustic peak sits near ℓ ≈ 220, an angular scale of about one degree. Maps live on HEALPix grids with N_pix = 12 N_side² equal-area pixels, and healpy.map2alm, alm2map, anafast and synfast do the transforms; a masked sky yields pseudo-C_ℓ that must be corrected for the mask's coupling matrix.

Numerical Computation

Legendre functions. Recur upward in ℓ at fixed m:

(ℓ − m) P_ℓ^m(x) = (2ℓ − 1) x P_{ℓ−1}^m(x) − (ℓ + m − 1) P_{ℓ−2}^m(x),

starting from P_m^m(x) = (−1)^m (2m − 1)!! (1 − x²)^{m/2} and P_{m+1}^m = (2m + 1) x P_m^m. This is stable, but the unnormalized values overflow double precision around degree 150 because of the double factorial. Recur on the fully normalized functions instead, which stay of order one, and fold the normalization into the recursion coefficients; the modified forward-column methods of Holmes and Featherstone reach degree well beyond 2000. In SciPy, scipy.special.lpmv(m, l, x) gives the Ferrers function with the Condon–Shortley phase, and scipy.special.sph_harm(m, l, theta, phi) takes the azimuth as theta and the polar angle as phi, the reverse of most textbooks. Newer SciPy releases add sph_harm_y(l, m, theta, phi) with the polar angle first and deprecate sph_harm; read the docstring of the version you have before trusting a plot.

Quadrature. Gauss–Legendre in x = cos θ with L + 1 nodes integrates polynomials of degree 2L + 1 in x exactly, and 2L + 1 or more equispaced points in φ integrate azimuthal orders up to 2L exactly, so the product rule computes every a_{ℓm} with ℓ ≤ L exactly for a band-limited f of degree ≤ L. Equiangular grids with the Driscoll–Healy weights do the same with 2(L + 1) latitudes. The direct transform costs O(L³); with FFTs in φ it is O(L³) with a small constant, and the Driscoll–Healy fast Legendre transform reaches O(L² log² L). SHTns, libsharp and pyshtools are the production implementations.

Worked Examples

Verify the addition theorem numerically.

import numpy as np
from scipy.special import sph_harm, eval_legendre
l = 5
th1, ph1, th2, ph2 = 0.7, 1.1, 2.0, 4.2                  # polar, azimuth for two points
x = np.array([np.sin(th1)*np.cos(ph1), np.sin(th1)*np.sin(ph1), np.cos(th1)])
y = np.array([np.sin(th2)*np.cos(ph2), np.sin(th2)*np.sin(ph2), np.cos(th2)])
lhs = sum(sph_harm(m, l, ph1, th1) * np.conj(sph_harm(m, l, ph2, th2)) for m in range(-l, l+1))
rhs = (2*l + 1)/(4*np.pi) * eval_legendre(l, x @ y)
print(lhs.real, rhs, abs(lhs.imag))                       # equal to ~1e-16; imaginary part zero

Expand a function by Gauss–Legendre quadrature and check Parseval.

import numpy as np
from scipy.special import sph_harm
L = 20
xg, wg = np.polynomial.legendre.leggauss(L + 1)          # nodes and weights in x = cos(theta)
theta = np.arccos(xg)
phi = np.linspace(0, 2*np.pi, 2*L + 2, endpoint=False)
TH, PH = np.meshgrid(theta, phi, indexing='ij')
W = wg[:, None] * (2*np.pi / phi.size)                     # weights sum to 4π
f = np.exp(np.cos(TH)) * np.cos(PH)**2                     # smooth, so coefficients decay super-exponentially
a = {(l, m): np.sum(f * np.conj(sph_harm(m, l, PH, TH)) * W)
     for l in range(L + 1) for m in range(-l, l + 1)}
rec = sum(c * sph_harm(m, l, PH, TH) for (l, m), c in a.items())
print(np.max(np.abs(rec.real - f)))                        # ~1e-12
print(sum(abs(c)**2 for c in a.values()), np.sum(f**2 * W))   # Parseval: equal to the same tolerance

Only m ∈ {0, ±2} are nonzero, because cos²φ = (1 + cos 2φ)/2; seeing power at other m means an argument-order or phase mistake.

Solving Δu = 0 in the unit ball with boundary data f. Expand f = Σ a_{ℓm} Y_ℓ^m and set u(r, ω) = Σ a_{ℓm} r^ℓ Y_ℓ^m(ω). For the exterior problem with decay at infinity, use r^{-ℓ-1} instead. The degree-0 term is the mean of f; the degree-1 terms are the dipole, whose exterior field decays like r^{-2}.

Checklist

• Normalization (orthonormal, 4π-normalized, Schmidt semi-normalized) and Condon–Shortley phase stated.
• Software argument order verified with a known value, e.g. Y₁⁰ = √(3/(4π)) cos θ.
• Degree of the product of two band-limited functions accounted for in the quadrature rule.
• Normalized recursion used above degree 100.
• Power spectrum, not individual coefficients, used for any rotation-invariant comparison.
• Truncation ringing considered when the function has sharp features; taper in ℓ if needed.
• Mask or missing data handled by a coupling correction, not ignored.

Common Mistakes

• Swapping θ and φ in sph_harm and getting a plausible-looking but wrong expansion.
• Comparing coefficients between codes with different phase conventions and concluding one is buggy.
• Using an equiangular grid with uniform weights; the poles are oversampled and the integrals are wrong.
• Truncating at degree L and sampling with fewer than about 2L + 1 azimuthal points, which aliases high m into low m.
• Interpreting a_{ℓm} as physical without fixing the coordinate frame.
• Trusting unnormalized P_ℓ^m beyond degree 150.
• Expecting a low-order expansion to represent a point light or a sharp mask; it produces negative lobes.

Limits

This file covers scalar harmonics on S² with remarks on S^{n−1}. Vector and tensor harmonics (polarization, gravitational lensing), spin-weighted harmonics, and harmonics on other compact symmetric spaces follow the same representation-theoretic pattern but have their own bases and recursions. The general framework of Peter–Weyl and Gelfand pairs is in the groups file.