Spherical Harmonic Methods¶

Spherical harmonic methods are the natural spectral approach for problems defined on the unit sphere $S^2$. They are the foundation of global weather and climate models, where the governing equations live on the spherical Earth.

1. Spherical Coordinates and the Laplace-Beltrami Operator¶

Points on the unit sphere are parameterised by colatitude $\theta \in [0, \pi]$ (or latitude $\phi = \pi/2 - \theta$) and longitude $\lambda \in [0, 2\pi)$:

\[\mathbf{x} = (\sin\theta \cos\lambda, \, \sin\theta \sin\lambda, \, \cos\theta)\]

The surface gradient of a scalar $f$ on $S^2$ is:

\[\nabla_S f = \frac{1}{\sin\theta} \frac{\partial f}{\partial \lambda} \hat{\boldsymbol{\lambda}} + \frac{\partial f}{\partial \theta} \hat{\boldsymbol{\theta}}\]

The Laplace-Beltrami operator (surface Laplacian) is:

\[\nabla_S^2 f = \frac{1}{\sin\theta} \frac{\partial}{\partial \theta}\!\left(\sin\theta \frac{\partial f}{\partial \theta}\right) + \frac{1}{\sin^2\theta} \frac{\partial^2 f}{\partial \lambda^2}\]

This operator is the generalisation of the Laplacian to the curved spherical geometry.

2. Spherical Harmonics¶

The spherical harmonics $Y_\ell^m(\theta, \lambda)$ are the eigenfunctions of the Laplace-Beltrami operator:

\[\nabla_S^2 Y_\ell^m = -\ell(\ell+1) Y_\ell^m\]

They are defined as:

\[Y_\ell^m(\theta, \lambda) = N_\ell^m \, P_\ell^m(\cos\theta) \, e^{im\lambda}\]

where: - $\ell = 0, 1, 2, \ldots$ is the total wavenumber (degree) - $m = -\ell, \ldots, \ell$ is the zonal wavenumber (order) - $P_\ell^m$ are the associated Legendre polynomials - $N_\ell^m$ is a normalisation constant

Normalisation¶

Using the fully normalised convention (standard in geodesy and atmospheric science):

\[N_\ell^m = \sqrt{\frac{(2\ell+1)}{4\pi} \frac{(\ell-|m|)!}{(\ell+|m|)!}}\]

Orthogonality¶

Spherical harmonics form a complete orthonormal basis on $L^2(S^2)$:

\[\int_{S^2} Y_\ell^m(\hat{\mathbf{x}}) \overline{Y_{\ell'}^{m'}(\hat{\mathbf{x}})} \, d\Omega = \delta_{\ell\ell'} \delta_{mm'}\]

where $d\Omega = \sin\theta \, d\theta \, d\lambda$ is the solid angle element.

Properties¶

Parity: $Y_\ell^m(-\hat{\mathbf{x}}) = (-1)^\ell Y_\ell^m(\hat{\mathbf{x}})$
Complex conjugate: $\overline{Y_\ell^m} = (-1)^m Y_\ell^{-m}$
Eigenvalue: $\nabla_S^2 Y_\ell^m = -\ell(\ell+1) Y_\ell^m$
Truncation: Truncation at degree $L$ gives $(L+1)^2$ basis functions

3. Associated Legendre Polynomials¶

The associated Legendre polynomials $P_\ell^m : [-1,1] \to \mathbb{R}$ satisfy the differential equation:

\[\frac{d}{d\mu}\!\left[(1-\mu^2) \frac{dP_\ell^m}{d\mu}\right] + \left[\ell(\ell+1) - \frac{m^2}{1-\mu^2}\right] P_\ell^m = 0, \quad \mu = \cos\theta\]

They are generated by the recurrence:

\[P_0^0(\mu) = 1, \quad P_m^m(\mu) = -(2m-1)\sqrt{1-\mu^2} P_{m-1}^{m-1}(\mu)\]

\[P_{\ell+1}^m(\mu) = \frac{(2\ell+1)\mu P_\ell^m(\mu) - (\ell+m) P_{\ell-1}^m(\mu)}{\ell - m + 1}\]

The Legendre recurrence is numerically stable when evaluated at the Gauss-Legendre quadrature nodes.

4. The Spherical Harmonic Transform (SHT)¶

A scalar field $f(\theta, \lambda)$ on the sphere is expanded as:

\[f(\theta, \lambda) = \sum_{\ell=0}^{L} \sum_{m=-\ell}^{\ell} \hat{f}_\ell^m \, Y_\ell^m(\theta, \lambda)\]

The spectral coefficients are:

\[\hat{f}_\ell^m = \int_0^{2\pi}\!\int_0^\pi f(\theta, \lambda) \overline{Y_\ell^m(\theta, \lambda)} \sin\theta \, d\theta \, d\lambda\]

Discrete SHT Algorithm¶

The SHT is split into two steps:

Zonal Fourier Transform (in $\lambda$): For each colatitude $\theta_j$, compute: $$F_m(\theta_j) = \frac{1}{2\pi} \int_0^{2\pi} f(\theta_j, \lambda) e^{-im\lambda} \, d\lambda \approx \text{FFT}$$
Legendre Transform (in $\theta$): For each $m$, compute: $$\hat{f}_\ell^m = \int_{-1}^{1} F_m(\mu) P_\ell^m(\mu) \, d\mu \approx \sum_j w_j F_m(\mu_j) P_\ell^m(\mu_j)$$

where $\mu_j$ are Gauss-Legendre nodes and $w_j$ are Gauss-Legendre weights.

The total cost is $O(L^3)$ for the Legendre transform step (dominant) and $O(L^2 \log L)$ for the zonal FFT step.

5. Gauss-Legendre Quadrature¶

The meridional integrals in the SHT are evaluated using Gauss-Legendre (GL) quadrature. For $J$ quadrature points $\mu_j = \cos\theta_j$ (zeros of $P_J(\mu)$) with weights $w_j$:

\[\int_{-1}^{1} g(\mu) \, d\mu \approx \sum_{j=1}^{J} w_j g(\mu_j)\]

This quadrature is exact for polynomials of degree $\leq 2J-1$. For the SHT with truncation $L$, we need $J \geq \lceil (3L+1)/2 \rceil$ quadrature points to exactly integrate the products $P_\ell^m P_{\ell'}^{m'}$ without aliasing.

The GL nodes and weights are computed from the eigenvalues and eigenvectors of the symmetric tridiagonal Jacobi matrix.

6. Spectral Differentiation on the Sphere¶

Zonal Derivative ($\partial/\partial\lambda$)¶

\[\frac{\partial f}{\partial \lambda} = \sum_{\ell, m} im \hat{f}_\ell^m \, Y_\ell^m\]

This is a simple multiplication by $im$ in spectral space, analogous to Fourier differentiation.

Meridional Derivative ($\partial/\partial\theta$)¶

The meridional derivative is more involved. In spectral space:

\[\frac{\partial Y_\ell^m}{\partial \theta} = -\sin\theta \frac{dP_\ell^m}{d(\cos\theta)} e^{im\lambda}\]

Using the recurrence for the derivative of $P_\ell^m$:

\[\frac{dP_\ell^m}{d\mu} = \frac{1}{1-\mu^2}\left[(\ell+1)\mu P_\ell^m - (\ell-m+1)P_{\ell+1}^m\right]\]

The Gradient and Divergence¶

The surface gradient and divergence can be computed using vector spherical harmonics $\mathbf{Y}_\ell^m$. For the vorticity-streamfunction formulation in atmospheric dynamics, the key relation is:

\[\nabla_S^2 \psi = \omega \quad \Longleftrightarrow \quad \hat{\psi}_\ell^m = \frac{-\hat{\omega}_\ell^m}{\ell(\ell+1)}\]

(with $\hat{\psi}_0^0 = 0$ for uniqueness).

7. Helmholtz/Poisson Equations on the Sphere¶

The Poisson equation on the sphere is:

\[\nabla_S^2 u = f\]

Since $Y_\ell^m$ are eigenfunctions with eigenvalue $-\ell(\ell+1)$, the solution is:

\[\hat{u}_\ell^m = \frac{-\hat{f}_\ell^m}{\ell(\ell+1)}, \quad \ell \geq 1\]

The $\ell = 0$ mode (global mean) is set to zero.

The Helmholtz equation:

\[(\nabla_S^2 - \alpha) u = f\]

\[\hat{u}_\ell^m = \frac{\hat{f}_\ell^m}{-\ell(\ell+1) - \alpha}\]

For $\alpha \neq 0$, the system is non-singular for all $\ell$.

8. Applications in Atmospheric and Geophysical Modelling¶

Spherical harmonic methods are the foundation of global spectral models in weather and climate prediction.

Barotropic Vorticity Equation¶

The simplest atmospheric model evolves the vorticity $\omega$ on the sphere:

\[\frac{\partial \omega}{\partial t} + J(\psi, \omega + f) = \nu \nabla_S^{2n} \omega\]

where $f = 2\Omega\cos\theta$ is the Coriolis parameter ($f$-plane/beta-plane models) and $J$ is the spherical Jacobian.

Quasigeostrophic Model on the Sphere¶

The QG potential vorticity equation:

\[q = \nabla_S^2\psi - \frac{\psi}{R^2} + f\]

\[\frac{\partial q}{\partial t} + J(\psi, q) = \text{dissipation} + \text{forcing}\]

In spectral space, PV inversion reduces to:

\[\hat{\psi}_\ell^m = \frac{-\hat{q'}_\ell^m}{\ell(\ell+1) + R^{-2}}\]

Triangular Truncation¶

Weather models typically use triangular truncation $T_L$: retain all modes with $\ell \leq L$ and $|m| \leq \ell$, giving $(L+1)^2$ complex coefficients.

The associated horizontal resolution is approximately $\Delta x \approx \pi a / L$ where $a$ is the Earth's radius. For example, T85 gives $\sim 150$ km resolution.

References¶

Durran, D.R. (2010). Numerical Methods for Fluid Dynamics. Springer.
Holton, J.R. & Hakim, G.J. (2012). An Introduction to Dynamic Meteorology. Academic Press.
Jakob-Chien, R. et al. (1995). Spectral transform solutions to the shallow water test suite. Monthly Weather Review, 123(6).
Shen, J., Tang, T. & Wang, L.L. (2011). Spectral Methods: Algorithms, Analysis and Applications. Springer.