I think the most accessible argument for completeness of spherical harmonics is to start with the Weierstrass approximation theorem, which asserts that on compact subsets $E$ of $R^n$ polynomials are dense in _sup_norm_, meaning the metric on functions $d(f,g)=\sup_{x\in E} |f(x)-g(x)|$. This applies to the sphere $E=S^{n-1}$. At some point in the development (as in the corresponding chapter in Stein-Weiss, for example), it is shown that every homogeneous polynomial $f$ can be written as $f=f_d + r^2f_{d-2} + r^4f_{d-4}+\ldots$ where $f_i$ is harmonic, and $r$ is radius. Thus, restricted to the sphere, every polynomial is pointwise-equal to a harmonic polynomial. Combining these two points, harmonic polynomials are dense in continuous functions on the sphere, with respect to sup-norm.
Depending what definition of "integral" one uses, one shows that continuous functions on the sphere are dense in $L^2$, so harmonic polynomials are dense in $L^2$. Then it is a small exercise to infer completeness.
A somewhat more sophisticated viewpoint has us note that the sphere is a homogeneous space $S^{n-1}=SO(n)/O(n-1)$ of orthogonal groups, so (thinking in terms of Frobenius reciprocity) the representations of $SO(n)$ appear which restrict to the trivial repn on $O(n-1)$, with multiplicity equal to the dimension of $O(n-1)$-fixed vectors, which is found to be $1$ (or $0$). The completeness in this version of the story is part of the general decomposition of $L^2(SO(n))$ (under the compact (Hilbert-Schmidt) operators coming from compactly-supported continuous functions on the group), and, instead, the issue becomes identifying irreducibles as corresponding to spherical harmonics. This is probably best done by using the fact that every irreducible has a highest weight, and, conversely, isomorphism classes are uniquely determined by highest weights.
This fancier viewpoint applies more broadly.
