How do we obtain radial part and angular part with spherical harmonics? In a website (sorry, I can't find it anymore) I've read this when talking about spherical harmonics:

if we want to determine simultaneous eigenfuntions of $L^2$ and $L_z$
we have to solve this system:
$$L^2|l,m\rangle=\hbar^2l(l+1)|l,m\rangle$$ $$L_z|l,m\rangle=\hbar m|l,m\rangle$$
in Schrodinger's representation:
$$-(\frac{1}{sin(\theta)}\frac{\partial}{\partial{\theta}}(sin(\theta)\frac{\partial}{\partial{\theta}})+\frac{1}{sin^2(\theta)}\frac{\partial^2}{\partial{\phi^2}})\psi_{lm}(r,\theta,\phi)=l(l+1)\psi_{lm}(r,\theta,\phi)$$
$$-i\frac{\partial}{\partial{\phi}}\psi_{lm}(r,\theta,\phi)=m\psi_{lm}(r,\theta,\phi)$$
we can note that operators $L^2$ and $L_z$ don't contain variable $r$
so the $\psi_{lm}(r,\theta,\phi)$ must be of the form
$\psi_{lm}(r,\theta,\phi)=f(r)Y_{lm}(\theta,\phi)$.

Question: Why if operators $L^2$ and $L_z$ don't contain variable $r$ the $\psi_{lm}(r,\theta,\phi)$ must be of the form
$\psi_{lm}(r,\theta,\phi)=f(r)Y_{lm}(\theta,\phi)$?
I think it is related to the fact that $L^2$ and $L_z$ have a complete set of simultaneous eigenfunctions but I don't understand why this leads to the separated form of the $\psi_{lm}(r,\theta,\phi)$.
 A: There are several possible more or less elaborated answers, I will discuss here a version which may be of interest, I hope, because it can be extended to other cases.  It happens because this approach relies upon general (though elementary) abstract results of Hilbert space theory.
Let $d\Omega$ be the unique rotationally invariant (Borel) measure on the two-dimensional sphere of unit radius $S^2$ such that $\int_{S^2} 1 d\Omega = 4\pi$.
We have an Hilbert space isomorphism,
$$L^2(\mathbb{R}^3, dx^3) \equiv L^2(S^2, d\Omega)\otimes L^2([0,+\infty), r^2dr)\tag{1}$$ where the isomorphism is the unique linear and bounded  extension of the map
$$L^2(S^2, d\Omega)\otimes L^2([0,+\infty), r^2dr) \ni Y\otimes f \mapsto Y \cdot f \in L^2(\mathbb{R}^3, d^3x)$$
where
$$(Y \cdot f)(\theta, \phi, r) := Y(\theta,\phi) f(r)\:.$$
That is a special case of the general result
$$ L^2(X, dx) \otimes L^2(Y, dy)\equiv L^2(X\times Y, dx\otimes dy) \tag{2}$$
where the isomorphism
is again the unique linear and bounded extension of the map
$$L^2(X,dx)\otimes L^2(Y,dy) \ni f\otimes g \mapsto f \cdot g \in L^2(X\times Y, dx\otimes dy)$$
In the case (1), $\mathbb{R}^3 = S^2 \times [0,+\infty)$ and $d^3x= d\Omega \otimes r^2 dr$ as is well known passing from Cartesian to polar coordinates.
Let us focus attention on (1). If one writes down the explicit expression of the three selfadjoint angular momentum operators $L_x,L_y,L_z$ in $L^2(\mathbb{R}^3, dx^3)\equiv  L^2(S^2, d\Omega)\otimes L^2([0,+\infty), r^2dr)$ he/she discovers the following interesting fact.
$$L_k = {\cal L}_k \otimes I_{L^2([0,+\infty), r^2dr)}\:,\tag{3}$$
for some operators ${\cal L}_k$ which are differential operators on the dense subspace of smooth (complex valued) functions in $L^2(S^2, d\Omega)$.
For instance
$${\cal L}_z = - i \hbar \frac{\partial}{\partial \phi}$$ and similar identities are valid for the other components, where only derivatives with respect to the angles $\phi$ and $\theta$ take place: these are  the variables used to define wavefunctions in  $L^2(S^2, d\Omega)$.
From the general theory of strongly continuous unitary representations of compact topological groups, in this case $SU(2)$, as the three $-i{\cal L}_k$ define a representation of its Lie algebra, one can construct a Hilbert basis of
simultaneous eigenvectors of ${\cal L}^2$ and ${\cal L}_z$ in $  L^2(S^2, d\Omega)$. This  basis is the very well known family of spherical harmonics $$Y^\ell_m(\theta,\phi).$$
Notice that the coordinate $r$ does not appear just because we are working in a Hilbert space where it does not exist!
Let us conclude. As a general result related with (2), if $L^2(X\times Y, dx\otimes dy)$ is a Hilbert space and $\{Y_a\}_{a\in A} \subset L^2(X, dx)$ is a Hilbert basis
whereas $\{R_b\}_{b\in B} \subset L^2(Y, dy)$ is another Hilbert basis, then
$$\{Y_a \otimes R_b\}_{(a,b) \in A\times B} \subset L^2(X, dx) \otimes L^2(Y, dy)\equiv L^2(X\times Y, dx\otimes dy)$$
is a Hilbert basis as well for $ L^2(X\times Y, dx\otimes dy)$.
Specializing the result to (1), we see that a Hilbert basis of $L^2(\mathbb{R}^3, d^3x)$ is made of the products $$Y^\ell_m(\theta, \phi) f_n(r)$$ where $\{f_n\}_{n \in N}$ is any given Hilbert basis of $L^2([0,+\infty). r^2dr)$.
Notice that, in view of (3),
$$L^2 (Y^\ell_m f_n) = ({\cal L}^2 \otimes I) Y^\ell_m \otimes f_n = \hbar^2 \ell(\ell+1) (Y^\ell_m f_n)$$
and
$$L_z (Y^\ell_m f_n) = ({\cal L}_z \otimes I) Y^\ell_m \otimes f_n = \hbar m (Y^\ell_m f_n)$$
More generally, with the same argument,
$$L^2 (Y^\ell_m f)  = \hbar^2 \ell(\ell+1) (Y^\ell_m f)$$
and
$$L_z (Y^\ell_m f) =  \hbar m (Y^\ell_m f)$$
for every $f\in L^2([0,+\infty), r^2dr)$.
In summary,  the functions $Y^\ell_m(\theta,\phi)f_n(r)$ form a symultaneous Hilbert basis of eigenvectors of $L^2$ and $L_z$ nomatter the choice of the Hilbert basis $\{f_n\}_{n\in N}$ in $L^2([0,+\infty), r^2 dr)$.
With a slightly more elaborated argument, it is possible to choice the functions $f_n$ also depending on $\ell$ and possibly $m$, changing the original functions $f_n$ in every eigenspace of $L^2$ and $L_z$.
