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.
