Spherical harmonics form a complete orthonormal basis for solving Laplace's equation in spherical coordinates. Expansions in Legendre polynomials and associated Legendre functions allow systematic solution of electrostatic problems with spherical symmetry, including multipole expansions.
Start with Legendre polynomials for azimuthally symmetric problems, then generalize to full angular dependence. Apply to conducting sphere and dielectric sphere problems to verify orthogonality and convergence.
Spherical harmonics are specific to electrostatics (they apply to any Laplacian problem). Assuming convergence without checking domain of validity.
Your prerequisite work on separation of variables showed that Laplace's equation ∇²V = 0 can be broken apart into independent ordinary differential equations when the geometry fits a coordinate system. In Cartesian coordinates this gave sinusoids; in spherical coordinates it gives something richer. Writing V(r,θ,φ) = R(r)·Θ(θ)·Φ(φ) and separating, you find that the radial equation gives power law solutions R ∝ rˡ or r^(−ℓ−1), the azimuthal equation gives Φ ∝ e^(imφ) (m an integer), and the polar equation gives the associated Legendre functions P_ℓᵐ(cos θ). The product Θ · Φ — normalized — is what we call a spherical harmonic Y_ℓᵐ(θ,φ). The integer ℓ ≥ 0 is the angular momentum quantum number; |m| ≤ ℓ gives the projection. For each ℓ, there are 2ℓ+1 values of m.
The key property that makes spherical harmonics so powerful is orthonormality: if you integrate the product of two different harmonics over the full sphere (all angles), you get zero; if you integrate a harmonic times its own complex conjugate, you get one. This is the same structure as Fourier series but on the surface of a sphere. Because they form a complete basis, any smooth function on the sphere — any arbitrary boundary condition you might impose on a spherical surface — can be expanded as a sum of spherical harmonics. This turns the problem of finding the electrostatic potential with a given boundary condition on a sphere into a coefficient-matching exercise.
To solve a typical problem — say, a conducting sphere in a uniform external field — you write the general solution as V = Σ (Aˡ rˡ + Bˡ r^(−ℓ−1)) Y_ℓᵐ(θ,φ). Far from the sphere, the potential must approach the uniform field −E₀z = −E₀r cos θ, which you recognize as the ℓ=1 term since P₁(cos θ) = cos θ. Near the origin, terms that blow up as r → 0 must vanish (or vice versa at r → ∞). You then apply the boundary condition on the sphere surface (V = constant for a conductor), use orthogonality to extract each coefficient, and you are done. The solution is built systematically from the expansion rather than guessed.
The importance of this technique extends far beyond electrostatics. The same Legendre polynomials and spherical harmonics appear in quantum mechanics as the angular part of atomic wave functions (the s, p, d, f orbital shapes you know are |Y_ℓᵐ|² plotted on a sphere), in gravitational potential theory for planetary shapes, and in acoustics for sound radiation patterns. The multipole expansion in electrostatics — expressing a localized charge distribution's far-field potential as a sum of monopole, dipole, quadrupole terms — is precisely an expansion in spherical harmonics: each ℓ term falls as r^(−ℓ−1) at large r. Once you recognize this structure, you see the same mathematics recurring across physics wherever a problem has spherical geometry.