Spherical harmonics and spherical Fourier neural operators

Before this lesson: lesson 08 (the FFT, spectral derivatives, aliasing and the Fourier neural operator). Helpful: linear algebra lesson 19 (equivariance) and calculus lesson 11 (Gauss–Hermite quadrature, a close relative of the Gauss–Legendre rule used here). Time: about 100 minutes, plus the notebook.

---

Why this matters #

• Weather and climate models are on the sphere. ERA5, the standard reanalysis used to train ML weather models, has $721 \times 1440$ points (a $0.25°$ grid). Every row has 1440 points. The east–west spacing is about 28 km at the equator and about 120 m on the row next to the pole. Methods that treat the grid as a flat image waste points near the poles and see a false seam.
• Spectral models. ECMWF's IFS computes derivatives, diffusion and its semi-implicit time step in spherical-harmonic space. Its high-resolution forecast uses the truncation TCo1279, about 9 km grid spacing (ECMWF documentation, as of 2026). The dynamical core of the hybrid model NeuralGCM (Kochkov et al., 2024) is also spectral, at resolutions from 2.8° to 0.7°.
• ML weather and climate models. The SFNO (Bonev et al., 2023) replaced the flat 2D FFT of FourCastNet by the SHT. The paper reports stable autoregressive rollouts over a simulated year (1460 steps of 6 hours). SFNO is also the base of the ACE climate emulator (Watt-Meyer et al., 2023). Lesson 16 puts these models side by side.
• Geophysics and seismology. Global models of the mantle from seismic tomography (for example S40RTS, up to degree 40), Earth's gravity field (EGM2008, to degree 2159) and the main geomagnetic field (IGRF, to degree 13) are all published as spherical-harmonic coefficients. The free oscillations of the Earth after a large earthquake are labelled by their degree $l$.
• Equivariant networks. Real spherical harmonics are the features of the E(3)-equivariant networks of lesson 10.

1. Predict first #

2. The picture #

The grid problem #

Use the colatitude $\theta \in [0, \pi]$ ($\theta = 0$ at the North Pole, $\theta = \pi/2$ on the equator) and the longitude $\phi \in [0, 2\pi)$. The latitude is $90° - \theta$. The area element of the unit sphere is $d\Omega = \sin\theta\,d\theta\,d\phi$: a band near a pole has little area because $\sin\theta$ is small there.

A regular latitude–longitude grid puts the same number of points on every row. So the points crowd together near the poles.

Interactive figure: loading…

So a lat–lon field is not a nice periodic image. Latitude is not periodic at all: the row at the North Pole is not a neighbour of the row at the South Pole. The right basis must be built for the sphere.

From the circle to the sphere #

On one circle of latitude, a field is a periodic function of $\phi$. So it is a Fourier series in $\phi$, exactly as in lesson 08: modes $\cos m\phi$ and $\sin m\phi$ with an integer order $m$. The coefficient of each mode is a function of the colatitude. For each $m$ we expand this function in a second family, the associated Legendre functions $\bar P_l^m(\cos\theta)$ with $l = m, m+1, \dots$ The product of the two is a spherical harmonic:

The degree $l = 0, 1, 2, \dots$ is the total wavenumber. The order $m$, with $-l \le m \le l$, is the east–west wavenumber. Each $Y_l^m$ is a standing wave on the sphere. Its zero set (the nodal lines) has $|m|$ great circles through the poles and $l - |m|$ circles of latitude.

The first few, up to constant factors:

• $Y_0^0$: a constant (the global mean);
• $Y_1^0 \propto \cos\theta$: warm in the north, cold in the south;
• $Y_1^1 \propto \sin\theta\cos\phi$ and $Y_1^{-1} \propto \sin\theta\sin\phi$: an east–west tilt;
• $Y_2^0 \propto 3\cos^2\theta - 1$; $Y_2^2 \propto \sin^2\theta\cos 2\phi$.

Interactive figure: loading…

Why the poles cause no trouble. $\bar P_l^m(\cos\theta)$ contains the factor $\sin^{|m|}\theta$. So a mode with a large order $m$ is almost zero near the poles. Fine east–west waves cannot live on the tiny circles near a pole, and the basis does not need them there. A degree $l$ has a wavelength of about $2\pi a/l$ on a sphere of radius $a$: on Earth ($a = 6371$ km), $l = 639$ gives about 63 km.

3. Definitions and main results #

3.1 The Laplacian on the sphere #

Write $x = \cos\theta$ (so $x = 1$ at the North Pole and $x = -1$ at the South Pole). On the unit sphere the Laplace–Beltrami operator (the Laplacian of the sphere) is

For the second form use $\partial_\theta = -\sin\theta\,\partial_x$ and $\sin^2\theta = 1 - x^2$. The area element is $d\Omega = \sin\theta\,d\theta\,d\phi = dx\,d\phi$ (up to the sign of $dx$).

3.2 Associated Legendre functions and spherical harmonics #

Definition. For integers $0 \le m \le l$, with $P_l$ the Legendre polynomial of degree $l$, the normalised associated Legendre function is

It solves the associated Legendre equation

and is bounded at $x = \pm 1$. The real spherical harmonics $Y_l^m$ are the products of section 2. They are orthonormal: $\int_{S^2} Y_l^m\,Y_{l'}^{m'}\,d\Omega = \delta_{ll'}\,\delta_{mm'}$. The factor $\sqrt2$ for $m \ne 0$ is there because $\int_0^{2\pi}\cos^2 m\phi\,d\phi = \pi$, while $\int_0^{2\pi} 1\,d\phi = 2\pi$.

3.3 Theorem 1: spherical harmonics are eigenfunctions of the Laplacian #

Why this is true (intuition). In longitude, $\cos m\phi$ is an eigenfunction of $\partial_\phi^2$ with eigenvalue $-m^2$, as on a circle. In latitude, the associated Legendre equation is exactly the condition that the remaining part gives $-l(l+1)$.

Proof.

1. Take $m \ge 0$ and write $Y = \bar P(x)\,\cos m\phi$ with $\bar P = \bar P_l^m$ (the $\sin m\phi$ case is the same). Then $\partial_\phi^2 Y = -m^2\,Y$.
2. With the $x$ form of the operator, $\nabla^2 Y = \Big[\big((1 - x^2)\,\bar P'\big)' - \dfrac{m^2}{1 - x^2}\,\bar P\Big]\cos m\phi$.
3. $\bar P$ solves the associated Legendre equation, so the bracket equals $-l(l+1)\,\bar P$. A constant factor does not change this, so the normalisation plays no role.
4. Hence $\nabla^2 Y = -l(l+1)\,Y$. On a sphere of radius $a$, each of the two derivatives carries a factor $1/a$, which gives $-l(l+1)/a^2$. $\blacksquare$

Consequences.

• The eigenvalue depends on $l$ only. All $2l + 1$ harmonics of degree $l$ share it. This is the rotation symmetry of the sphere at work (section 3.8).
• Spectral calculus. If $f = \sum_{l,m} a_l^m Y_l^m$, then $\nabla^2 f$ has the coefficients $-\tfrac{l(l+1)}{a^2}\,a_l^m$.
• Poisson's equation is a division. The streamfunction $\psi$ of a vorticity field $\zeta$ solves $\nabla^2\psi = \zeta$, so $\psi_l^m = -a^2\,\zeta_l^m/(l(l+1))$ for $l \ge 1$. At $l = 0$ the eigenvalue is 0: the global mean of $\zeta$ must be zero, and the mean of $\psi$ is free (set it to 0).
• Diffusion is a damping per degree. $\partial_t u = \nu\nabla^2 u$ multiplies each coefficient by $e^{-\nu l(l+1)t/a^2}$. Spectral models apply (hyper)diffusion this way.
• The longitude derivative is also diagonal: it maps the pair $(Y_l^m, Y_l^{-m})$ to $(-m\,Y_l^{-m}, m\,Y_l^m)$. The latitude derivative is not diagonal: $\sin\theta\,\partial_\theta$ maps degree $l$ to a combination of degrees $l - 1$ and $l + 1$.

3.4 Triangular truncation #

A field is band-limited to $N$ if it has no harmonics of degree above $N$. The triangular truncation $T_N$ keeps all $0 \le l \le N$ and $-l \le m \le l$: that is $\sum_{l=0}^N (2l + 1) = (N + 1)^2$ real coefficients. "T79" or "T639" in the weather literature means this truncation.

Triangular truncation is the natural choice because its resolution is the same everywhere on the globe. The addition theorem says $\sum_{m=-l}^{l} \big(Y_l^m(\theta, \phi)\big)^2 = (2l + 1)/(4\pi)$ at every point. So no place, and no direction, is preferred. The set of all $T_N$ fields is also mapped to itself by every rotation (section 3.8).

Interactive figure: loading…

3.5 Computing $\bar P_l^m$: the stable recursion #

The closed formula of section 3.2 is useless for large $l$. The factorials overflow (171! is larger than the largest float64), and high derivatives of polynomials lose all digits to cancellation. Use a recursion for each order $m$, directly on the normalised functions:

1. Start: $\bar P_0^0 = 1/\sqrt{4\pi}$.
2. Diagonal: $\bar P_m^m = \sqrt{\dfrac{2m+1}{2m}}\,\sqrt{1 - x^2}\;\bar P_{m-1}^{m-1}$.
3. First step off the diagonal: $\bar P_{m+1}^m = \sqrt{2m + 3}\;x\,\bar P_m^m$.
4. Three-term recursion in $l$, for $l \ge m + 2$:

The numbers $\alpha_l^m$ and $\beta_l^m$ are fixed; they are not the coefficients $a_l^m$ of a field. All values stay of moderate size, so nothing overflows. The values agree with SciPy (after the phase and the normalisation are matched) to about $10^{-14}$ for $l \le 20$. This is the recursion of torch-harmonics and, in an equivalent form, of SHTns.

Very high degrees

For $N$ in the thousands, $\bar P_m^m \propto \sin^m\theta$ underflows near the poles. At the polar node of a 1280-latitude Gauss grid, $\sin\theta \approx 1.9 \times 10^{-3}$, and $\sin^{1279}\theta \approx 10^{-3487}$. These values are negligible: the harmonic is zero there for all practical purposes. Libraries skip them, or use scaled (extended-exponent) arithmetic.

3.6 Gauss–Legendre quadrature #

The SHT must compute integrals over latitude, $\int_{-1}^1 g(x)\,dx$, from samples. A quadrature rule with $K$ nodes $x_j$ and weights $w_j$ approximates $\int_{-1}^1 g(x)\,dx \approx \sum_{j=1}^K w_j\,g(x_j)$. With equispaced nodes and well-chosen weights you can integrate every polynomial of degree $K - 1$ exactly. If you may also choose the nodes, you can do twice as well.

Definition. The Gauss–Legendre rule with $K$ nodes uses the $K$ roots of the Legendre polynomial $P_K$ as nodes, and the weights $w_j = \int_{-1}^1 \ell_j(x)\,dx$, where $\ell_j$ is the Lagrange polynomial that is 1 at $x_j$ and 0 at the other nodes.

Proof

1. Let $p$ have degree $\le 2K - 1$. Divide by $P_K$: $p = q\,P_K + r$, with $\deg q \le K - 1$ and $\deg r \le K - 1$.
2. $P_K$ is orthogonal to every polynomial of lower degree, so $\int_{-1}^1 q\,P_K\,dx = 0$. Hence $\int p = \int r$.
3. At the nodes $P_K(x_j) = 0$, so $p(x_j) = r(x_j)$. Hence $\sum_j w_j\,p(x_j) = \sum_j w_j\,r(x_j)$.
4. $r$ has degree $\le K - 1$, so it equals its Lagrange interpolant $\sum_j r(x_j)\,\ell_j$. Integrate: $\int r = \sum_j r(x_j)\int\ell_j = \sum_j w_j\,r(x_j)$.
5. Steps 2–4 give $\sum_j w_j\,p(x_j) = \int p$. For degree $2K$, take $p = P_K^2$: its integral is positive, but the rule gives 0. $\blacksquare$

Interactive figure: loading…

3.7 Theorem 2: the SHT on a Gauss grid is exact #

Grid. $N_{\text{lat}}$ Gauss–Legendre nodes $x_j$ (colatitudes $\theta_j = \arccos x_j$) with weights $w_j$, and $N_{\text{lon}}$ equispaced longitudes $\phi_k = 2\pi k/N_{\text{lon}}$. This is a Gaussian grid.

Analysis (grid → coefficients). For each order $m \ge 0$:

1. FFT along each row: $F_m(x_j) = \dfrac{1}{N_{\text{lon}}}\sum_k f(x_j, \phi_k)\,e^{-im\phi_k}$.
2. Legendre step, one matrix–vector product per $m$: $c_l^m = 2\pi\sum_j w_j\,F_m(x_j)\,\bar P_l^m(x_j)$.
3. Real coefficients: $a_l^0 = c_l^0$ and, for $m > 0$, $a_l^m = \sqrt2\,\operatorname{Re} c_l^m$, $a_l^{-m} = -\sqrt2\,\operatorname{Im} c_l^m$.

Synthesis (coefficients → grid): set $c_l^0 = a_l^0$ and $c_l^m = (a_l^m - i\,a_l^{-m})/\sqrt2$; then $F_m(x_j) = \sum_{l=m}^N c_l^m\,\bar P_l^m(x_j)$, and $f(x_j, \phi_k) = F_0 + 2\operatorname{Re}\sum_{m>0}F_m(x_j)\,e^{im\phi_k}$ (an inverse real FFT).

Proof of Theorem 2

1. The continuous formula. By orthonormality $a_l^m = \int f\,Y_l^m\,d\Omega$. Do the longitude integral first. With $F_m(x) = \frac{1}{2\pi}\int_0^{2\pi} f\,e^{-im\phi}\,d\phi$ and $c_l^m = 2\pi\int_{-1}^1 F_m(x)\,\bar P_l^m(x)\,dx$, one checks $a_l^m = \sqrt2\operatorname{Re} c_l^m$ and $a_l^{-m} = -\sqrt2\operatorname{Im} c_l^m$ (use $\cos m\phi = \operatorname{Re} e^{-im\phi}$ and $\sin m\phi = -\operatorname{Im} e^{-im\phi}$).
2. Longitude. On a row, $f$ is a trigonometric polynomial with orders $|m'| \le N$. The discrete sum $\frac{1}{N_{\text{lon}}}\sum_k e^{i(m' - m)\phi_k}$ is 1 if $m' - m$ is a multiple of $N_{\text{lon}}$ and 0 otherwise (lesson 08, Theorem 1). Since $|m' - m| \le 2N < N_{\text{lon}}$, only $m' = m$ survives: step 1 gives the exact $F_m(x_j)$.
3. Latitude. $F_m(x)$ is a combination of $\bar P_{l'}^m(x)$ with $l' \le N$. Each $\bar P_l^m$ is $(1 - x^2)^{m/2}$ times a polynomial of degree $l - m$. So the product $\bar P_{l'}^m\bar P_l^m$ of two functions with the same $m$ is $(1 - x^2)^m$ times a polynomial of degree $l + l' - 2m$: a polynomial of degree $l + l' \le 2N$. (For odd $m$ one $\bar P_l^m$ alone is not a polynomial; the product is.)
4. Quadrature. $N_{\text{lat}}$ Gauss nodes are exact up to degree $2N_{\text{lat}} - 1 \ge 2N + 1$. So step 2 equals the integral of proof step 1. The synthesis is the exact evaluation of a finite sum, and the coefficients determine the field. $\blacksquare$

Two traps show that both bounds are sharp (we use them again in section 4):

• With $N_{\text{lat}} = N$ nodes, the nodes are the roots of $P_N$, and $Y_N^0 \propto P_N(\cos\theta)$ is zero at every grid point. The grid cannot see it.
• With $N_{\text{lon}} = 2N$, $\sin(N\phi_k) = \sin(\pi k) = 0$ at every grid longitude. The mode $Y_N^{-N}$ is invisible.

Cost. The FFTs cost $O(N_{\text{lat}}N_{\text{lon}}\log N_{\text{lon}}) = O(N^2\log N)$. The Legendre step costs about $(N + 1 - m)\,N_{\text{lat}}$ per order, so $O(N^3)$ in total: it dominates for large $N$. The symmetry $\bar P_l^m(-x) = (-1)^{l+m}\,\bar P_l^m(x)$ saves a factor 2, and on GPUs the step is a batch of dense matrix products.

Products, dealiasing and other grids

The product of two $T_N$ fields has degrees up to $2N$, as in 1D (lesson 08). To get the $T_N$ coefficients of a product without aliasing, the integrands have degree up to $3N$. Then you need $N_{\text{lat}} \ge (3N + 1)/2$ and $N_{\text{lon}} \ge 3N + 1$: the quadratic grid, the spherical form of the 2/3 rule (Boyd, §18.12). The minimal grid of Theorem 2 is the linear grid. ECMWF's IFS uses an even finer cubic octahedral grid, which also has fewer points per row near the poles. An equiangular grid such as ERA5's has no Gauss nodes. Exact quadrature on it needs about twice as many latitudes (Driscoll and Healy, 1994).

3.8 Rotations, spherical convolution and the SFNO layer #

Rotations keep the degree. Let $R$ be a rotation and $f_R(\hat x) = f(R^{-1}\hat x)$ the rotated field. The rotated version of a degree-$l$ harmonic is again a combination of degree-$l$ harmonics. Reason: the Laplacian commutes with rotations (it is defined by the geometry of the sphere), so it maps the eigenspace of $-l(l+1)$ to itself. A rotation also keeps the inner product $\int fg\,d\Omega$. Hence the vector $a_l = (a_l^{-l}, \dots, a_l^l)$ of the degree-$l$ coefficients changes by an orthogonal $(2l+1) \times (2l+1)$ matrix:

$D^l(R)$ is the real Wigner-D matrix; lesson 10 studies it in detail. Two consequences:

• The power spectrum $E_l = \sum_m (a_l^m)^2 = \|a_l\|^2$ does not change under rotations. It is the standard diagnostic of small-scale energy. An ML forecast trained with the MSE often loses power at high $l$: it is blurry.
• A rotation about the polar axis does not mix orders: it turns each pair $(a_l^m, a_l^{-m})$ by the angle $m\alpha$. Any other rotation mixes all orders $m$ of the same degree.

Proof. On degree $l$, $T$ is the matrix $c_l I$. Rotate first, then filter: $c_l\,D^l(R)\,a_l$. Filter first, then rotate: $D^l(R)\,(c_l\,a_l)$. A multiple of the identity commutes with every matrix, so the two are equal for every $l$. $\blacksquare$

A diagonal factor $c_{lm}$ that changes with $m$ is a diagonal matrix that is not a multiple of the identity. It does not commute with a $D^l(R)$ that mixes the orders. Schur's lemma (lesson 10) says more: the only linear maps of the degree-$l$ space that commute with all rotations are multiples of the identity.

Interactive figure: loading…

Spherical convolution. On the circle, a convolution is a multiplication of Fourier coefficients (lesson 08). On the sphere, the convolution of a field with a kernel that is symmetric about the pole (a zonal kernel) multiplies the coefficients of degree $l$ by one number per degree. That number is proportional to the $m = 0$ coefficient of the kernel at degree $l$ (Driscoll and Healy, 1994; SFNO paper, §4.1). So degree-wise filters are exactly the spherical convolutions.

The SFNO layer. With channels $v(\hat x) \in \R^{d_v}$ at each point $\hat x$ of the sphere, one layer is

where $\hat v_l^m \in \R^{d_v}$ are the SHT coefficients of the channels and $R_l \in \R^{d_v \times d_v}$ is a learned matrix for each degree, with no dependence on $m$. It is the FNO layer of lesson 08 with the FFT replaced by the SHT. $R_l$ acts on the channel index and $D^l(R)$ on the order index, so they commute: the linear part is exactly rotation-equivariant (Theorem 3 for each channel pair). The pointwise skip $W$ and the activation $\sigma$ are equivariant on the continuous sphere; on a discrete grid the activation creates high degrees that alias, so the full network is only approximately equivariant (SFNO paper, §4.2).

Answers to section 1. (1) No. Latitude is not periodic: a 2D FFT makes the North-Pole row a neighbour of the South-Pole row and creates a false seam. The rows near the poles are tiny circles with as many points as the equator, and the flat Fourier modes are not smooth at the poles. Use the FFT only in longitude, and associated Legendre functions with Gauss quadrature in latitude. (2) $N + 1$ Gauss latitudes and $2N + 1$ longitudes (Theorem 2). For exact products, use $(3N + 1)/2$ latitudes and $3N + 1$ longitudes. (3) Yes. A factor $c_l$ that depends only on the degree commutes with every rotation (Theorem 3). A factor that depends on $m$ does not, except for rotations about the polar axis.

4. Worked example #

Take $f(\theta, \phi) = 1 + 3\cos\theta + \sin^2\theta\,\cos 2\phi$. Its highest degree is $N = 2$.

Step 1: the coefficients by hand. We need three harmonics:

• $Y_0^0 = \bar P_0^0 = 1/\sqrt{4\pi}$. So $1 = \sqrt{4\pi}\,Y_0^0$ and $a_0^0 = \sqrt{4\pi} = 3.5449$.
• $\bar P_1^0 = \sqrt{3/(4\pi)}\,x$, so $Y_1^0 = \sqrt{3/(4\pi)}\cos\theta$. Then $3\cos\theta = 3\sqrt{4\pi/3}\,Y_1^0$ and $a_1^0 = 3\sqrt{4\pi/3} = 2\sqrt{3\pi} = 6.1400$.
• $\bar P_2^2 = \sqrt{\tfrac{5}{4\pi}\tfrac{0!}{4!}}\,(1 - x^2)\cdot 3 = \tfrac14\sqrt{15/(2\pi)}\,\sin^2\theta$ (because $\tfrac{d^2}{dx^2}P_2 = 3$). So $Y_2^2 = \sqrt2\,\bar P_2^2\cos 2\phi = \tfrac14\sqrt{15/\pi}\,\sin^2\theta\cos 2\phi$, and $a_2^2 = 4\sqrt{\pi/15} = 1.8306$.

All other coefficients are 0.

Step 2: the grid. Theorem 2 needs $N_{\text{lat}} \ge 3$ and $N_{\text{lon}} \ge 5$. The 3-node Gauss rule has nodes $x = -\sqrt{3/5},\ 0,\ \sqrt{3/5}$ (that is $\pm 0.7746$ and 0) and weights $5/9,\ 8/9,\ 5/9$. Check its limit:

• degree 4: $\sum_j w_j x_j^4 = 2 \cdot \tfrac59 \cdot \tfrac{9}{25} = \tfrac25 = \int_{-1}^1 x^4\,dx$ (exact);
• degree 6: $\sum_j w_j x_j^6 = 2 \cdot \tfrac59 \cdot \tfrac{27}{125} = 0.24$, but $\int_{-1}^1 x^6\,dx = \tfrac27 = 0.2857$.

Degree $2N = 4$ is all we need.

Step 3: analysis of $a_1^0$. On each row, the mean over the 5 longitudes is $F_0(x) = 1 + 3x$: the $\cos 2\phi$ term sums to zero over 5 equispaced points. Then

which equals $\int_{-1}^1 (1 + 3x)\,x\,dx = 2$. So $c_1^0 = 2\pi\sqrt{3/(4\pi)} \cdot 2 = 2\sqrt{3\pi} = 6.1400 = a_1^0$. ✓

Step 4: analysis of $a_2^2$. On each row, $F_2(x) = \tfrac15\sum_k \sin^2\theta\cos 2\phi_k\,e^{-2i\phi_k} = \tfrac12(1 - x^2)$. (Only the $e^{2i\phi}$ half of $\cos 2\phi = \tfrac12(e^{2i\phi} + e^{-2i\phi})$ survives; $e^{-4i\phi}$ sums to 0 on 5 points.) Then

The sum is $2 \cdot \tfrac59 \cdot (0.4)^2 + \tfrac89 \cdot 1 = \tfrac{16}{15}$, which is $\int_{-1}^1(1 - x^2)^2\,dx$. So $c_2^2 = 1.2944$ and $a_2^2 = \sqrt2\,\operatorname{Re} c_2^2 = 1.8306$. ✓

Step 5: the Laplacian. Multiply each coefficient by $-l(l+1)$: $a_0^0 \to 0$, $a_1^0 \to -2 \times 6.1400 = -12.280$, $a_2^2 \to -6 \times 1.8306 = -10.983$. Synthesis gives

You can check this with the operator of section 3.1. The inverse Laplacian gives back $f$ up to a constant: the $l = 0$ part is lost.

Step 6: the two traps.

• $N_{\text{lat}} = 2$: the nodes $\pm 1/\sqrt3$ are the roots of $P_2$, so $Y_2^0 \propto P_2(\cos\theta)$ is zero at every grid point. Also, the 2-node rule gives $\sum_j w_j x_j^4 = 2/9$, not $2/5$.
• $N_{\text{lon}} = 4$: $\sin 2\phi_k = \sin(k\pi) = 0$ at all four longitudes. The orders $m = 2$ and $m = -2$ alias, and the $\sin 2\phi$ part of a field is invisible.

5. In code #

The code checks the quadrature limit and computes the worked example with the recursion and the SHT of Theorem 2. Only NumPy and SciPy's roots_legendre are used.

import numpy as np
from scipy.special import roots_legendre

# (a) Gauss-Legendre with K nodes is exact up to degree 2K - 1.
x, w = roots_legendre(3)
print(x, w)                            # [-0.7746 0. 0.7746] [0.5556 0.8889 0.5556]
print(w @ x**4, w @ x**6, 2 / 7)       # 0.4 (= 2/5, exact)  0.24  0.2857

# (b) Normalised associated Legendre functions by the recursion (no (-1)^m phase).
def pbar_table(N, x):
    s, P = np.sqrt(1 - x**2), np.zeros((N + 1, N + 1, x.size))
    P[0, 0] = 1 / np.sqrt(4 * np.pi)
    for m in range(1, N + 1):
        P[m, m] = np.sqrt((2 * m + 1) / (2 * m)) * s * P[m - 1, m - 1]
    for m in range(N):
        P[m + 1, m] = np.sqrt(2 * m + 3) * x * P[m, m]
        for l in range(m + 2, N + 1):
            a = np.sqrt((4 * l * l - 1) / (l * l - m * m))
            b = np.sqrt(((l - 1) ** 2 - m * m) * (2 * l + 1) / ((l * l - m * m) * (2 * l - 3)))
            P[l, m] = a * x * P[l - 1, m] - b * P[l - 2, m]
    return P

# (c) Analysis on the Gauss grid: rfft in longitude, quadrature in latitude.
def sht(f, N):
    nlat, nlon = f.shape
    x, w = roots_legendre(nlat)
    P, F = pbar_table(N, x), np.fft.rfft(f, axis=1) / nlon
    A = np.zeros((N + 1, 2 * N + 1))                   # A[l, N + m] = a_l^m
    for m in range(N + 1):
        c = 2 * np.pi * (P[m:, m] * w) @ F[:, m]
        if m == 0:
            A[:, N] = c.real
        else:
            A[m:, N + m], A[m:, N - m] = np.sqrt(2) * c.real, -np.sqrt(2) * c.imag
    return A

x, _ = roots_legendre(3)
th, ph = np.arccos(x)[:, None], 2 * np.pi * np.arange(5)[None, :] / 5
f = 1 + 3 * np.cos(th) + np.sin(th) ** 2 * np.cos(2 * ph)
A = sht(f, 2)
print(A[0, 2], A[1, 2], A[2, 4])       # 3.5449 6.1400 1.8306
print(np.sqrt(4 * np.pi), 2 * np.sqrt(3 * np.pi), 4 * np.sqrt(np.pi / 15))   # the same
print(np.abs(A).sum() - A[0, 2] - A[1, 2] - A[2, 4])   # all other coefficients: 4e-15

# (d) The Laplacian: multiply degree l by -l(l+1).
print(A[[0, 1, 2], [2, 2, 4]] * -np.array([0, 2, 6]))   # [0. -12.2799 -10.9835]

In the notebook you also write the synthesis, check a round trip with random $T_{20}$ coefficients (error about $10^{-14}$), solve the Poisson equation and test the rotation property of Theorem 3.

6. Where it shows up #

• Spectral dynamical cores. The IFS (ECMWF) and the core of NeuralGCM compute derivatives, diffusion and the semi-implicit solver in spectral space, and the nonlinear terms on a Gaussian grid (the "transform method"). NeuralGCM's core is written in JAX, so the learned physics is trained through it (lesson 03).
• ML weather and climate models. SFNO-based models (FourCastNet with SFNO, the ACE climate emulator) use the degree-wise filters of section 3.8. The open-source library torch-harmonics provides a differentiable SHT on GPUs. Spherical power spectra are a standard check of blurriness in lesson 16.
• Positional encodings on the globe. Spherical harmonics of the latitude and longitude of a point are used as location encodings for geospatial ML (Rußwurm et al., 2024).
• Seismology and geophysics. Global tomography models of seismic velocity in the mantle, gravity models of the geoid, and the geomagnetic reference field are spherical-harmonic expansions. A seismic "spheroidal mode" $_nS_l$ of the free oscillations of the Earth has the angular shape of a degree-$l$ harmonic.
• Equivariant networks. In lesson 10, the real harmonics $Y_l$ of an edge direction are the features of E(3)-equivariant networks, and the matrices $D^l(R)$ of section 3.8 appear again.
• Practice tools. torch-harmonics (PyTorch), SHTns (fast C library), pyshtools (geodesy and geophysics).

7. Common mistakes #

• Colatitude versus latitude. Gauss nodes are values of $x = \cos\theta = \sin(\text{latitude})$. Mixing the two flips or shifts the field.
• Wrong argument order in SciPy. The new function is sph_harm_y(l, m, colatitude, longitude). The older sph_harm(m, l, longitude, colatitude), still in old code, used the opposite order in both pairs.
• Mixing conventions. Orthonormal, $4\pi$-normalised and Schmidt semi-normalised functions differ by factors, and the phase $(-1)^m$ may be on or off. Match both before comparing.
• Too few latitudes, or Gauss weights on the wrong grid. Fewer than $N + 1$ Gauss latitudes, or Gauss weights on an equiangular grid, silently drop modes (section 4).
• Dividing by $l(l+1)$ at $l = 0$, or forgetting the radius. The constant mode is the null space of the Laplacian. On Earth the eigenvalue is $-l(l+1)/a^2$.
• Spectra from $m \ge 0$ storage. If you store complex $c_l^m$ for $m \ge 0$ only, then $E_l = |c_l^0|^2 + 2\sum_{m>0}|c_l^m|^2$: count each $m > 0$ twice.
• Factorial formulas at high degree. Use the normalised recursion.
• Treating the lat–lon grid as a periodic image. A 2D FFT (as in FNO or AFNO) creates a seam between the poles and artefacts near them. This is what SFNO fixes.
• Weights that depend on $m$. They break rotation equivariance (Theorem 3 and the quiz above).

8. Check yourself #

1. What is the eigenvalue of $\nabla^2$ for degree $l = 3$ on the unit sphere, and on Earth ($a = 6371$ km)?
   Answer

   $-l(l+1) = -12$ on the unit sphere. On Earth $-12/a^2 = -12/(6.371 \times 10^6\ \text{m})^2 \approx -2.96 \times 10^{-13}\ \text{m}^{-2}$. All seven harmonics of degree 3 share it.
2. A model uses the truncation T79. Give the minimal Gauss grid for an exact transform, and the minimal grid for exact (dealiased) quadratic products.
   Answer

   Linear grid: $N_{\text{lat}} \ge 80$ and $N_{\text{lon}} \ge 159$ (160 in practice, for the FFT). Quadratic grid: $N_{\text{lat}} \ge (3 \cdot 79 + 1)/2 = 119$ and $N_{\text{lon}} \ge 238$ (often $120 \times 240$).
3. On the unit sphere the vorticity is $\zeta = 6\sin^2\theta\cos 2\phi$. Find the streamfunction $\psi$ with $\nabla^2\psi = \zeta$.
   Answer

   $\zeta$ has degree 2 only (it is a multiple of $Y_2^2$), so divide by $-l(l+1) = -6$: $\psi = -\sin^2\theta\cos 2\phi$, plus any constant.
4. Why does a weight $c_{lm}$ that depends on $m$ break rotation equivariance, while $c_l$ does not?
   Answer

   A rotation acts on the degree-$l$ coefficients by the orthogonal matrix $D^l(R)$, which mixes the orders $m$. The filter is the diagonal matrix $\operatorname{diag}(c_{l,-l}, \dots, c_{l,l})$. It commutes with every $D^l(R)$ only if it is a multiple of the identity, that is, if $c_{lm} = c_l$ (Schur's lemma, lesson 10).
5. Which step of the SHT dominates at large $N$, and how does it scale?
   Answer

   The Legendre step: about $N^2/2$ coefficients times $N_{\text{lat}} \approx N$ latitudes, so $O(N^3)$. The FFTs cost only $O(N^2\log N)$. The parity symmetry halves the Legendre step, and GPUs run it as batched matrix products.
6. You compare the power spectra $E_l$ of an ML forecast and of the analysis at a lead time of 5 days, and the forecast has much less power for $l > 100$. What does this tell you?
   Answer

   The forecast is too smooth at scales shorter than about $2\pi a/100 \approx 400$ km. This is typical of models trained with the MSE: under uncertainty the MSE-optimal forecast is the conditional mean, which averages out small features. The rotation invariance of $E_l$ makes it a clean, orientation-free diagnostic.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is the eigenvalue of the Laplacian for $Y_l^m$ on a sphere of radius $a$?

   Show answer

2. How many real coefficients does the triangular truncation $T_N$ have?

   Show answer

3. Up to which degree is a $K$-node Gauss–Legendre rule exact?

   Show answer

4. What are the two steps of the spherical harmonic transform on a Gauss grid?

   Show answer

5. Which grid makes the SHT exact for a band limit $N$?

   Show answer

6. Why compute $\bar P_l^m$ with the three-term recursion instead of the closed formula?

   Show answer

7. Which spectral filters $a_l^m \mapsto c\,a_l^m$ commute with all rotations of the sphere?

   Show answer

8. How do you solve $\nabla^2\psi = \zeta$ on the unit sphere with the SHT?

   Show answer

9. What does the SFNO change compared with the FNO?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why spherical harmonics and not a 2D FFT on the lat–lon grid?" Latitude is not periodic and the grid is singular at the poles. Spherical harmonics are the eigenfunctions of the Laplacian on the sphere, and the triangular truncation has the same resolution everywhere.
• "What does SFNO change compared with FNO?" The SHT replaces the flat FFT, and the weights depend on the degree $l$ only, so the linear part commutes with rotations. The paper reports stable rollouts over one simulated year (1460 steps of 6 hours).
• "How would you show that an ML forecast is blurry?" Compare the spherical power spectrum $E_l$ of the forecasts with that of the analysis (for example ERA5). Models trained with the MSE lose power at high $l$, more at longer lead times.
• "Why is Poisson's equation easy on the sphere, and what is the catch?" Divide each coefficient by $-l(l+1)/a^2$. The catch is $l = 0$: the constant is the null space, so the mean of the right-hand side must be zero.
• "How do you choose a grid for the truncation $N$?" $N + 1$ Gauss latitudes and $2N + 1$ longitudes for linear terms; the quadratic grid ($(3N + 1)/2$ latitudes and $3N + 1$ longitudes) to avoid aliasing in products; about twice the latitudes on an equiangular grid.

10. Practice #

The notebook for this lesson runs in Colab and in the browser. It has five exercises, each with a check cell and a folded solution:

1. Build the Gauss–Legendre rule with the Golub–Welsch eigenvalue method and test its exactness up to degree $2K - 1$.
2. Compute the normalised associated Legendre functions with the recursion; compare with scipy.special.lpmv and check orthonormality.
3. Write the forward and inverse SHT on a Gauss grid; reproduce the worked example and get a round-trip error near $10^{-14}$.
4. Solve Poisson's equation on the sphere (vorticity to streamfunction) and handle the $l = 0$ mode.
5. Test numerically that a filter $c_l$ commutes with rotations and a filter $c_{lm}$ does not.

11. Watch and read #

Spherical harmonics as a solution to Laplace's equation on a sphereJeremy Maurer · 9 min

Watch how separation of variables gives the Fourier part in longitude and the associated Legendre part in latitude.

Spherical Harmonics and the Multipole ExpansionThe Gravity Hub · 35 min

A geophysics view: the gravity field of a planet as a sum of spherical harmonics, with the same degree and order.

Fourier Neural Operator (FNO) [Physics Informed Machine Learning]Steve Brunton · 18 min

A reminder of the FNO layer; the SFNO keeps it and changes only the transform and the weight sharing.

Reading

• Bonev et al., Spherical Fourier Neural Operators: Learning Stable Dynamics on the Sphere (ICML 2023): §3–4 for the SHT, spherical convolution and equivariance. arXiv 2306.03838
• Boyd, Chebyshev and Fourier Spectral Methods (2nd ed., free PDF), chapter 18: spherical harmonics, triangular truncation, Gaussian grids and dealiasing. PDF
• Schaeffer, Efficient spherical harmonic transforms aimed at pseudo-spectral numerical simulations (G-cubed, 2013; the SHTns paper): definitions, the recursion and the exactness argument. arXiv 1202.6522
• torch-harmonics notebooks (NVIDIA): the SHT, quadrature and a Poisson solver in PyTorch. GitHub