How to Fit Spherical Harmonics to a 3D Grid

A spherical harmonic is a function of direction on a sphere. A 3D field generally depends on both distance from the origin and direction, so there is no single spherical-harmonic fit until you decide how to handle the radius. The usual approach is to fit the angular pattern separately on one or more spherical shells.

This also clarifies “nearest”: it can mean the best single harmonic, or a truncated sum of harmonics that approximates the field. The fitting method below covers both.

Separate radial and angular dependence

Write a scalar field as an angular expansion whose coefficients may vary with radius:

$$f(r,\theta,\phi)=\sum_{\ell=0}^{L}\sum_{m=-\ell}^{\ell}a_{\ell m}(r)Y_{\ell m}(\theta,\phi).$$

Here $r$ is distance from the chosen origin, $\theta$ is the polar angle, and $\phi$ is the azimuth. This article uses the complex spherical harmonics normalized so that

$$\int_{S^2}Y_{\ell m}^{*}(\theta,\phi)Y_{\ell' m'}(\theta,\phi)\,d\Omega =\delta_{\ell\ell'}\delta_{mm'},$$

where the star denotes complex conjugation and $d\Omega=\sin\theta\,d\theta\,d\phi$. Other normalizations and real-valued bases are common; use one convention consistently.

For a continuous field at fixed $r$, the coefficient is

$$a_{\ell m}(r)=\int_{S^2}f(r,\theta,\phi)Y_{\ell m}^{*}(\theta,\phi)\,d\Omega.$$

On a 3D Cartesian grid, choose a shell around the radius of interest and use the field values in that shell. Repeat for several radii if you need the radial dependence of the coefficients.

Fit grid samples by weighted least squares

For each selected grid point $(x_i,y_i,z_i)$, calculate

$$r_i=\sqrt{x_i^2+y_i^2+z_i^2},\qquad \theta_i=\arccos(z_i/r_i),\qquad \phi_i=\operatorname{atan2}(y_i,x_i).$$

The origin has no defined direction, so do not use it in an angular fit. At each remaining point, evaluate the chosen basis functions $Y_{\ell m}(\theta_i,\phi_i)$. Collect the field values into a vector $\mathbf f$ and the basis values into a matrix $A$, with one row per sample and one column per retained harmonic. Estimate the coefficient vector $\mathbf a$ by minimizing

$$\sum_i w_i\left|f_i-\sum_{\ell,m}a_{\ell m}Y_{\ell m}(\theta_i,\phi_i)\right|^2.$$

In matrix form, solve $W^{1/2}A\mathbf a\approx W^{1/2}\mathbf f$, where $W$ is diagonal with entries $w_i$. Use a stable least-squares solver such as QR or singular-value decomposition rather than explicitly forming a matrix inverse.

The weights matter. A regular Cartesian grid does not sample directions uniformly: some parts of a shell may have more grid points than others. Give samples weights that approximate their represented solid angle. For small voxels in a thin shell, a useful approximation is proportional to voxel volume divided by $r_i^2$; the radial thickness is then effectively constant. A purpose-built spherical quadrature grid can instead use its supplied solid-angle weights. Simply giving every Cartesian point equal weight can bias the fit.

If the samples are well-distributed quadrature nodes and the weights approximate $d\Omega$, you can estimate a coefficient directly with $a_{\ell m}\approx\sum_i w_i f_iY_{\ell m}^{*}(\theta_i,\phi_i)$. For a general Cartesian subset, least squares is usually safer because the sampled basis functions may not be orthogonal.

Decide which harmonic is nearest

To find the best single harmonic, fit each candidate mode separately and compare its weighted residual sum of squares. Choose the mode with the smallest residual. When the sampling and weights make the basis orthonormal, this is equivalent to choosing the mode with the largest coefficient magnitude $|a_{\ell m}|$. For nonuniform or incomplete sampling, that shortcut need not hold.

To approximate the field with a sum of harmonics, fit all modes through some maximum degree $L$, then increase $L$ until the residual is sufficiently small. Higher degrees represent finer angular variation, but they can also fit noise and become poorly determined if there are too few well-distributed samples. A useful practical check is to compare residuals on held-out samples or across successive values of $L$.

Worked example: a constant plus a dipole

Suppose the field on a unit sphere is

$$f(\theta,\phi)=1+2\cos\theta.$$

Under the normalization above, $Y_{00}=1/\sqrt{4\pi}$ and $Y_{10}=\sqrt{3/(4\pi)}\cos\theta$. Therefore

$$f=\sqrt{4\pi}\,Y_{00}+2\sqrt{\frac{4\pi}{3}}\,Y_{10}.$$

The only nonzero coefficients are $a_{00}=\sqrt{4\pi}$ and $a_{10}=2\sqrt{4\pi/3}$. If you require one harmonic only, compare their squared magnitudes: the constant mode has $4\pi$, while the dipole mode has $16\pi/3$. Thus $Y_{10}$ is the better single-mode fit in the continuous, uniformly weighted spherical $L^2$ sense. It captures $4/7$ of the field’s total squared norm. A grid fit should reach the same result approximately when its shell sampling and weights represent the sphere well.

Common mistakes

• Fitting a volume as if it were one sphere. Keep radius fixed, fit separate shells, or explicitly model radial dependence.
• Using unweighted Cartesian samples. Uneven directional coverage changes the effective fit and can favor some angles.
• Comparing coefficients from incompatible conventions. Normalization and phase conventions affect coefficient values; state and keep the convention fixed.
• Using too many modes for the data. There are $(L+1)^2$ complex modes through degree $L$. Sparse or clustered samples can make the fit unstable.
• Ignoring the origin or coordinate orientation. The origin has no angular coordinates, and changing the chosen axis changes the mode labels and coefficients.