In a spherically symmetric central potential why do we look for eigenfunctions of the angular momentum operator? In finding the solutions to the wave equation for a spherically symmetric potential $V(r)$, we look for the eigenfunctions of $\hat L^2$ and $\hat L_z$ operators. However, what is the reasoning behind this? The time-independent Schrodinger equation: 
$$\hat H \Psi = E \Psi$$
This is an eigenvalue equation and any eigenfunction of $\hat H$ this equation is obviously a solution for this equation. Now when there is a spherically symmetric potential $V(r)$, I understand that because of separation of variables we can look for solutions of the form:
$$\Psi (\vec r,\theta, \phi) = R(\vec r) Y(\theta, \phi)$$
Now to find the angular part why do we look for eigenstates of $\hat L^2$ and $\hat L_z$ operators (I know that they commute)?
 A: Well, to put it shortly, it's just because $L^2$ and $L_z$ are two observables that have no $r$ dependence. Since the kinetic terms of $H$ can be written as functions of $r$ and $L^2$, and the potential depends only on $r$, it is clear that the Hamiltonian commutes with the angular momenta. As such, it makes sense to write down $Y$ in terms of these two. 
But if we'd like to go a little deeper, it is a bit related to representation theory. The fact that $L^2$ and $L_z$ form a complete description of vector rotations (they represent the Lie algebra $\mathfrak{so}(3)$ of the Lie group $\mathrm{SO}(3)$) and therefore are the most natural way to express the angular dependence of a rotation-invariant space in terms of orthonormalizable functions. 
After reading @ZeroTheHero 's comment, i figured it would be instructive to provide some intuition as to why angular momentum does not depend on radius. Remember the classical definition
$$\mathbf L=\mathbf r\times\mathbf p$$
Let's then give and estimate of $L$. Suppose an atom or something has some characteristic radius $R$, and the electron moves with some characteristic momentum $p$. Then we estimate
$$L\sim Rp$$
But the de Broigle relation gives us an estimate of momentum in terms of the characteristic length
$$p\sim\frac{h}{R}$$
Where $h$ is Planck's constant. As such, out estimate for angular momentum amounts to
$$L\sim h,$$
Which has no explicit dependence on characteristic lengths of the system. As such it is reasonable to say it has no $r$ dependence, even though I've not given a rigorous proof. 
A: One objective of quantum mechanics is to fully labelled states of the system.  In some cases energy is not enough: for instance in the hydrogen atom some energy levels occur more than once; the same happens in the 3d harmonic oscillator.  In such cases, simply giving the energy is not enough to completely identify the state so we look for other ``quantum numbers'' to refine the identification.   
Of course one would prefer to label states with constants of motion, since they are constant so what we call state $n\ell m$ at the start we can use the same name at the end.  Hence the idea of using a complete set of commuting operators - the Hamiltonian and others that commute with it - so as to uniquely identify the quantum states with constants of motion.  
Of course when the potential is spherical $V(r)$, what enters in the radial part of the Schrödinger equation is actually the effective potential 
$$
V_{eff }=V(r)+\frac{\hbar^2}{2m}\frac{\ell(\ell+1)}{r^2} \tag{1}
$$
and so you see that the radial Schrödinger equation actually explicitly contains an angular momentum part so it is quite natural to look for functions $\psi(r,\theta,\phi)$ which are also eigenfunctions $Y^\ell_m(\theta,\phi)$ of $\hat L^2$ as these will provide the correct extra factor in Eq.(1) when separating the variables.   Since there are many functions $Y^\ell_m(\theta,\phi)$ with the same $\hbar^2\ell(\ell+1)$, one uses the eigenvalue $\hbar m$ of $\hat L_z$ to distinguish them.
