A small thread on spherical harmonics & why there are always 2n+1 harmonics of degree n. However, let's start a simpler case - Fourier series ("circular harmonics"), of which spherical harmonics are a direct generalization to a sphere.

(1/n)

Fourier series are defined for periodic functions f(t), which are the same as functions on a circle. Let's take the usual unit circle x² + y² = 1 and work with functions f(x,y) instead; x = cos(t), y = sin(t) performs the conversion from f(t) to f(x,y).

(2/n)

As a consequence of the celebrated Stone–Weierstrass theorem, any function on a circle can be approximated by polynomials to any degree. Thus, we could try to use the standard basis of polynomials {1, x, y, x², xy, y², x³, x²y, xy², y³, ...} for this.

(3/n)

However, this basis is highly redundant on a circle, because of x² + y² = 1: we can express y² = 1 - x², then xy² = x(1-x²) = x - x³, then y³ = y(1-x²) = y - x²y, and so on. Whenever we have xⁿyᵐ with m ≥ 2, we can reduce it to (xⁿ-xⁿ⁺²)yᵐ⁻².

(4/n)

Thus, any polynomial on a circle can be uniquely represented as p(x) + q(x)y, where p and q are univariate polynomials. This implies that instead of taking {xⁿ, xⁿ⁻¹y, xⁿ⁻²y², ..., xyⁿ⁻¹, yⁿ} as a basis for n-th degree polynomials, we can take just ...

(5/n)

... {xⁿ, xⁿ⁻¹y}, since all others are redundant. So, we have one basis function {1} for the zero-degree polynomials, two more functions {x, y} to include degree 1, two more {x², xy} to include degree 2, and so on.

(6/n)

Now, compare this to Fourier series: a constant for degree zero, and two functions {cos(nt), sin(nt)} for any positive degree. This is no coincidence: this simply corresponds to taking a different basis of the same vector space (polynomials with degree ≤ n).

(7/n)

Remember that cos(nt) can be expressed as a polynomial in cos(t) and sin(t) via trigonometric formulas, or in x and y when we change to cartesian coordinates. On a circle, sines/cosines of the angular coordinate are the same as polynomials in cartesian coordinates.

(8/n)

Again: {1, x, y, x², xy} and {1, cos(t), sin(t), cos(2t), sin(2t)} are two different bases for the vector space of polynomials of degree ≤ 2. The latter is better since it is an orthonormal basis, but it is still polynomial, and can be written as {1, x, y, x² - y², 2xy}.

(9/n)

The idea of all this is that the cos/sin basis used in Fourier is better, but we don't need it to just calculate the dimension of the vector space (i.e. to calculate the number of parameters required to approximate a function).

(10/n)

What does this have to do with spherical harmonics? They are exactly the same, except that angular coordinates are are more weird, but we can still use polynomials all along.

(11/n)

The defining equation now is x² + y² + z² = 1, and by the very same logic as before we can represent any polynomial uniquely as p(x,y) + q(x,y)z. Thus, to form a basis for polynomials up to degree n, we can ignore monomials that include zⁿ with n ≥ 2.

(12/n)

Explicitly, we could take (n+1) functions {xⁿ, xⁿ⁻¹y, xⁿ⁻²y², ..., xyⁿ⁻¹, yⁿ} together with n functions {xⁿ⁻¹z, xⁿ⁻²yz, ..., xyⁿ⁻²z, yⁿ⁻¹z} (2n+1 total) for degree n. Now we have {1} for degree 0, then add {x,y,z} for degree 1, then ...

(13/n)

... add {x², xy, y², xz, yz} for degree 2, and so on. Compare this to the number of spherical harmonics in each dimension: 1, 3, 5, and so on. They are the same!

(14/n)

Again, this is no coincidence: spherical harmonics of degree up to n are a basis of spherical polynomials in (x,y,z) of degree up to n. To find explicit cartesian expressions for the harmonics, see https://en.wikipedia.org/wiki/Spherical_harmonics#Spherical_harmonics_in_Cartesian_form or https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics.

(15/n)

A few notes. Firstly, I'm always talking about "polynomials of degree ≤ n", since "polynomials of degree n" don't form a vector space. Secondly, the bases we constructed are not orthonormal, as opposed to Fourier basis or spherical harmonics.

(16/n)

However, you can orthonormalize them easily using the Gram-Schmidt procedure, and get a basis which is no worse than spherical harmonics for computational purposes.

(17/n)

You can even evaluate all the integrals (the L² dot products needed for orthogonalization) without actually doing integration! E.g. the integral of x² over the sphere is the same as the integral of y² or z² (by a symmetry argument).

(18/n)

Thus integral(x²) = integral(x²+y²+z²)/3 = integral(1)/3 = (area of sphere)/3 = 4π/3. Other integrals are harder, but I suspect they still can be evaluated using symmetry & x² + y² + z² = 1.

(19/19)