Elementary proof of completeness of the spherical harmonics for a fixed $L$ This question asked how to prove the completeness of the spherical harmonics in the sense that $\{Y_{Lm}\}$ spans the set of square-integrable functions on the sphere. I'm looking for an elementary proof of what I think should be a much simpler theorem. Suppose that we're interested in the set of all functions on the sphere that have the same $L$, i.e., the set of all complex-valued functions that are eigenfunctions of the Laplacian with eigenvalue $-L(L+1)$. I'm looking for the most elementary possible proof that for this fixed $L$, the set $\{Y_{Lm}\}$ for varying $m$ is a basis. That is, I want an elementary proof that the multiplicity is $2L+1$. Ideally I would like this to be at a freshman physics level, intelligible even to people who haven't had linear algebra. (So the language would actually be less sophisticated that in my statement of the problem above.) I would be perfectly happy with a proof restricted to $L=1$.
The closest I've been able to come is an approach that is kind of complicated, and I'm not sure it's even right. First we reduce the problem to real-valued functions by forming linear combinations of $Y_{Lm}$ and $Y_{L,-m}$. Next we argue that if the hypothesis were false, we would be able to witness its failure with a state $W$ of definite $m$ (i.e., one whose azimuthal variation is $e^{im}$) that is independent of $Y_{Lm}$. Now form linear combinations of $W$ with $Y_{Lm}$ such that the linear combination $X$ vanishes at some $\theta$. This can be done at almost all $\theta$. Finally, tune the choice of $\theta$ so that for this particular value of $\theta$, not just $X$ but $\partial X/\partial \theta$ vanishes as well. Now argue that by uniqueness of solutions to 2nd order differential equations, $X$ vanishes identically, which is a contradiction. I'm not sure that the last part even works, since the standard uniqueness theorems may not work for a space with this topology.
Is there a simpler way to go about this?
I also thought about physical arguments involving state counting. If the $L=1$ state has multiplicity 3, then coupling two spin-1 states gives 9 states, and this makes sense because you get spin-0, spin-1, and spin-2 couplings, for a total of 1+3+5=9 states. But all this really seems to do is prove that if the multiplicity of spin $L$ is $2L+1$, we get self-consistency --- it doesn't seem to show that $2L+1$ is necessary.
Maybe there is some sort of nice argument saying that if we can couple two states of spin 1 to make spin 0, then this spin 0 state must be unique?
 A: As you know, the spherical harmonics$$Y_\ell^m(\theta, \phi) = \sqrt{{2\ell + 1(\ell + |m|)!}\over{4\pi (\ell + |m|)!}} \cdot P_\ell^m(\cos \theta)e^{im\phi}$$(here $\ell = 0, 1, 2, \ldots$ and $m = 0, \pm1, \pm2, \ldots, \pm\ell$) are orthonormal on the sphere$$\int_0^{2\pi} \int_0^\pi (Y_{\ell'}^{m'} (\theta, \phi))^* (Y_{\ell^{\prime\prime}}^{m^{\prime\prime}}(\theta, \phi)) \sin\theta\,d\theta\,d\phi = \delta^{m'm^{\prime\prime}}\delta_{\ell' \ell^{\prime\prime}}$$and are alleged to be complete.
https://dspace.mit.edu/bitstream/handle/1721.1/90373/8-07-fall-2005/contents/readings/completeness.pdf
$$\sum_{\ell = 0}^\infty \sum_{m = -\ell}^\ell (Y_\ell^m(\theta, \phi))^* (Y_\ell^m(\theta', \phi')) = {1\over{\sin \theta}} \delta(\theta - \theta') \delta(\phi - \phi')$$Set $m = 0$ and get$$Y_\ell^0(\theta, \phi) = \sqrt{{2\ell + 1}\over{4\pi}} \cdot P_\ell(\cos\theta),$$$$2\pi \int_0^\pi (Y_{\ell'}^0(\theta, \phi))^* (Y_{\ell^{\prime\prime}}^0(\theta, \phi)) \sin\theta\,d\theta = \delta_{\ell'\ell^{\prime\prime}},$$or$${{2\ell + 1}\over2} \int_0^\pi P_{\ell'}(\cos\theta)P_{\ell^{\prime\prime}}(\cos\theta)\sin\theta\,d\theta = \delta_{\ell' \ell^{\prime\prime}},$$which by $x = \cos \theta$ becomes the orthogonality statement$${{2\ell + 1}\over2} \int_0^\pi P_{\ell'}(x)P_{\ell^{\prime\prime}}(x)\,dx = \delta_{\ell' \ell^{\prime\prime}}.$$The Wikipedia article on Legendre polynomials reports the completeness relation$$\sum_{\ell = 0}^\infty {{2\ell + 1}\over2} P_\ell(x) P_\ell(y) = \delta(x - y),$$which I have offhand no idea how to prove.
https://en.wikipedia.org/wiki/Legendre_polynomials
I gather from what you have written that you are well aware of the theory of spherical harmonics is central to the quantum theory of angular momentum. From the $3$-dimensional analogues of the fundamental commutation relation$$[\textbf{x}, \textbf{p}] \equiv \textbf{xp} - \textbf{px} = i\hbar \textbf{I},$$it follows that that the operators$$\textbf{L}_x = \textbf{yp}_z - \textbf{zp}_y, \quad\textbf{L}_y = \textbf{zp}_x - \textbf{xp}_z, \quad\textbf{L}_z = \textbf{xp}_y - \textbf{yp}_x$$(inherited from the classical construction $\textbf{L}= \textbf{r} \times \textbf{p}$) fail to commute:$$[\textbf{L}_x, \textbf{L}_y] = i\hbar \textbf{L}_z, \quad [\textbf{L}_y, \textbf{L}_z] = i\hbar\textbf{L}_x, \quad [\textbf{L}_z, \textbf{L}_x] = i\hbar\textbf{L}_y.$$But from this it follows that $\textbf{L}^2 \equiv \textbf{L}_x^2 + \textbf{L}_y^2 + \textbf{L}_z^2$ does commute with each of those operators, so $\{\textbf{L}^2, \textbf{L}_z\}$ (one could replace $\textbf{L}_z$ with any real linear combination of $\{\textbf{L}_x, \textbf{L}_y, \textbf{L}_z\}$; it is by standard convention that we select $\textbf{L}_z$) comprise a commuting set of operators, and it becomes possible to ask for simultaneous eigenvectors/functions of these operators, objects $f$ such that$$\textbf{L}^2f = \lambda f \text{ and }\textbf{L}_zf = \mu f.$$An algebraic argument—that makes use of the operators$$\textbf{J}_+ = \textbf{L}_x + i\textbf{L}_y, \quad \textbf{J}_- = \textbf{L}_x - i\textbf{L}_y$$and that is as pretty as it is elementary—leads quickly to the conclusion that necessarily$$\begin{align}
\lambda = \hbar^2\ell(\ell + 1) & : \ell = 0, {1\over2}, 1, {3\over2}, 2, \ldots \\
\mu = \hbar m & : m = -\ell, -\ell + 1, \ldots, \ell - 1, \ell
\end{align}$$and accounts for the $(2\ell + 1)$-fold degeneracy of the states with a given $\ell$. For details, see i.e. §4.3.1 in David Griffiths's "Introduction to Quantum Mechanics". Elaboration of the argument leads in the cases $\ell = 0, 1, 2,\ldots$ explicit constructions of the eigenfunctions $f_\ell^m(r, \theta, \phi)$, whence to the spherical harmonics. But those arguments fail when $\ell = {1\over2}, {3\over2}, {5\over2}, \ldots$; one is led to $(2\ell + 1)$-dimensional vectors; the associated operators are then represented not by differential operators, but by matrices; and one is led to spinors and to spin representations of the rotation group.
Vector representations are available also for every fixed integral value of $\ell$. "Completeness" acquires then the meaning standard to the theory of finite-dimensional vector spaces, as when says of a basis $\{\textbf{e}_1, \textbf{e}_2, \ldots, \textbf{e}_{2\ell + 1}\}$ that it is complete (spans the vector space).
Moving now closer to your problem, as I understand it: one says of functions $\{u_1(x), u_2(x), \ldots\}$ that are orthogonal on the interval $[a, b]$$$\int_a^b u_i(x)u_j(x)\,dx = \delta_{ij}$$that they are complete if every well-behaved function $f(x)$ on $[a, b]$ can be developed$$f(x) = \sum_i f_i u_i(x) \text{ with }f_i = \int_a^b f(y)u_i(y)\,dy.$$Then$$\begin{align}
f(x) & = \sum_i f_i u_i(x) \\
& = \sum_i \int_a^b f(y)u_i(y)u_i(x)\,dy \\ 
& = \int_a^b \left\{\sum_i u_i(x)u_i(y)\right\}f(y)\,dy : \text{all nice }f(x),
\end{align}$$whence formally $\sum_i u_i(x)u_i(y) = \delta(x - y)$. If one or several of the elements of the set $\{u_i(x)\}$ are deleted then the remaining functions are complete only with respect to the space of functions $g(x)$ that are orthogonal to the deleted basis functions:$$\int_a^b g(y)u_{\text{deleted}}(y)\,dy = 0.$$You propose to delete all spherical harmonics except those with some prescribed value of $\ell$. You look, therefore, to the restricted class $\mathcal{G}$ of functions of the form$$g(\theta, \phi) = P_\ell(\cos\theta) \sum_{m = -\ell}^{m = +\ell} g_m e^{im\phi}.$$That the functions $e^{im\phi} : m = -\ell, \ldots, +\ell$ are orthogonal and complete in $\mathcal{G}$ is, I take it, evident. So your problem, as I understand it, seems to have evaporated. I have for this reason not attempted to work through the details of the arguments sketched in your final paragraphs.

The spherical harmonics are $Y_l^m(\theta,\varphi) = C_{lm} P_l^m(\cos\theta) e^{im\varphi},$ where $C_{lm}$ are some constants for normalization and $l = 0, 1, 2, \ldots$ and $m = -l, \ldots, +l.$
The completeness of $\mathcal B = \{ Y_l^m(\theta,\varphi)  \}$ depends on the completenesses of $\mathcal B_\theta = \{ P_l^m(\cos\theta) \}$ and $\mathcal B_\varphi = \{ e^{im\varphi} \}.$ If $H_1, H_2$ are two Hilbert spaces with bases $\mathcal B_1, \mathcal B_2$ then $H_1 \otimes H_2$ has basis $\mathcal B_1 \otimes \mathcal B_2 = \{ b1 \otimes b_2 \mid b_1 \in \mathcal B_1, \, b_2 \in \mathcal B_2 \}$.
The completenesses of $\mathcal B_\theta$ and $\mathcal B_\varphi$ follow from Sturm-Liouville theory.
Here we however have a small complication. The physicality of the solutions require that $|m| \leq l.$ Therefore we don't use the full $\mathcal B_\varphi$ for a given $l$ so we actually don't use a complete basis in $\varphi$. But if we took the whole $\mathcal B_\varphi$ we would have to set the coefficients for $P_l^m$ to $0$ when $|m| > l$ to get physical (bounded) solutions.

Imagine that you have the space of all $L^2$-integrable (with the right area-proportional measure) functions on the two-sphere, i.e. functions of $\theta$, $\phi$ with the right conditions, okay?
Now, on this Hilbert space of complex functions, you may define the operators $\textbf{L}^2$ and $\textbf{L}_z$. The operator $\textbf{L}_z$ is basically just $-i\hbar$ times the derivative with respect to $\phi$. The operator $\textbf{L}^2$, here it means the square of the angular momentum, is just some proper Laplacian on the two-sphere.
You may prove that $\textbf{L}^2$ and $\textbf{L}_z$ are Hermitian. So you may simultaneously diagonalize them.
The eigenvalues of $\textbf{L}^2$ may be seen to be $\ell(\ell+1)$ times $\hbar^2$ where $\ell=0,1,2,3,4\ldots.$ It is because you may take the simultaneous eigenstate of $\textbf{L}^2$ and $\textbf{L}_z$. The eigenvalues of $\textbf{L}_z$ are integers, $m$, okay? You may also notice that if you define the operators $\textbf{L}_+$ and $\textbf{L}_-$, i.e. $\textbf{L}_x$ plus minus $i\textbf{L}_y$, then the action of $\textbf{L}_+$ and $\textbf{L}_-$ does not change a vector's eigenvalue of $\textbf{L}^2$. But it changes the eigenvalue of $\textbf{L}_z$ by plus minus one, so they are raising and lowering operators.
Because the eigenvalue of $\textbf{L}_z^2$ cannot be greater than that of $\textbf{L}^2$, because their difference is positively or negatively definite, depending which is subtracted from which, it follows that something must prevent you from creating new eigenstates by the action of $\textbf{L}_+$. What prevents you is that $\textbf{L}_+$ acting on the highest-$\textbf{L}_z$ eigenstate in the multiplet with a given $\textbf{L}^2$ eigenvalue is zero.
But if it is zero, you may prove that $\textbf{L}^2$ acting on that state is the same as $\textbf{L}_z(\textbf{L}_z+1)$ acting on it, so the eigenvalue of $\textbf{L}^2$ must be $\ell(\ell+1)$—times $\hbar^2$, as always.
It is always guaranteed that any element of the Hilbert space—the $L^2$-integrable complex function on the sphere—may be decomposed to eigenstates of a Hermitian operator. The completeness relation holds for any Hermitian operator on any Hilbert space etc.
So the sum of $|\ell,m\rangle\langle \ell,m|$ over all allowed quantum numbers $\ell$, $m$ must give you the identity operator.
A: The spherical harmonics are $Y_l^m(\theta,\varphi) = C_{lm} P_l^m(\cos\theta) e^{im\varphi},$ where $C_{lm}$ are some constants for normalization and $l = 0, 1, 2, \ldots$ and $m = -l, \ldots, +l.$
The completeness of $\mathcal B = \{ Y_l^m(\theta,\varphi)  \}$ depends on the completenesses of $\mathcal B_\theta = \{ P_l^m(\cos\theta) \}$ and $\mathcal B_\varphi = \{ e^{im\varphi} \}.$ If $H_1, H_2$ are two Hilbert spaces with bases $\mathcal B_1, \mathcal B_2$ then $H_1 \otimes H_2$ has basis $\mathcal B_1 \otimes \mathcal B_2 = \{ b1 \otimes b_2 \mid b_1 \in \mathcal B_1, \, b_2 \in \mathcal B_2 \}$.
The completenesses of $\mathcal B_\theta$ and $\mathcal B_\varphi$ follow from Sturm-Liouville theory.
Here we however have a small complication. The physicality of the solutions require that $|m| \leq l.$ Therefore we don't use the full $\mathcal B_\varphi$ for a given $l$ so we actually don't use a complete basis in $\varphi$. But if we took the whole $\mathcal B_\varphi$ we would have to set the coefficients for $P_l^m$ to $0$ when $|m| > l$ to get physical (bounded) solutions.
