Ok. Let me "correct" (or specify) some things in your answer. Noting that the system has this type of symmetry (spherical symmetry). Decompose the wave function in products of wave functions, each of they related to some variable of the coordinate system is already a big step. But I assume that you are aware of this subtitles.
Let's see your second equation. The fact that $\Phi(\phi)=\Phi(\phi+2\pi)$ tell us that $m_{\ell}^2$ are a natural number. So is better (here comes the "correction") specify $\Phi$ as a $\Phi_{\ell}$, after all we have a family of solutions not only one (diverses eigenstate of the hamiltonian).
In the bulk of your text you said that mathematical and physical game are in some way different ones. This is not true here. Depend on what you call mathematics. And here, the fact that "...On mathematical ground, this magnetic quantum number can be any constant by the theory of partial differential equation..." is in the same game of physics.
In the physical game you are searching eigenstates of the energy observable and the correspondent spectra. You have a hamiltonian like that:
$$
H=-\frac{\hbar^2}{2m}\nabla^2+V(r)
$$
and then the eigenvector are some wavefunction that obeys:
$$
H\psi=-\frac{\hbar^2}{2m}\nabla^2\psi+V(r)\psi=E\psi
$$
Knowing that the laplacian in spherical coordinates are:
$$
\nabla^2
= {1 \over r^2} {\partial \over \partial r}
\left(r^2 {\partial \over \partial r} \right)
+ {1 \over r^2 \sin \theta} {\partial \over \partial \theta}
\left(\sin \theta {\partial \over \partial \theta} \right)
+ {1 \over r^2 \sin^2 \theta} {\partial^2 \over \partial \varphi^2}
$$
$$
= D_r
+ {1 \over r^2 \sin \theta} D_{\theta}
+ {1 \over r^2 \sin^2 \theta} {\partial^2 \over \partial \varphi^2}
$$
With $[D_r,\,D_{\theta}]=[D_r,\,{\partial^2 \over \partial \varphi^2}]=[D_{\theta},\,{\partial^2 \over \partial \varphi^2}]$=0. So, we can determine simultaneously the eigenvectors of this three operators. Physically, $i\hbar{\partial^2 \over \partial \varphi^2}$, is the $z$-component of the angular momentum and
$$
-\frac{\hbar^2}{2m}{1 \over r^2 \sin \theta} {\partial \over \partial \theta}
\left(\sin \theta {\partial \over \partial \theta} \right)
-\frac{\hbar^2}{2m}{1 \over r^2 \sin^2 \theta} {\partial^2 \over \partial \varphi^2}
$$
is the term $\frac{L^2}{2mr^2}$ on energy conservation, where $L$ is the angular momentum.
Ok. You can play with this representation of the angular operator and derive a lot of relations. But one thing is very clear here. The eigenvector of the angular part of the laplacian can be determined no matter what is contained in the potential energy (if the potential energy depend on $r$ only). This is because $V(r)$ commute with the angular part of the laplacian. We have something as:
$$
r^2\nabla^2_{\theta \phi}Y(\theta,\phi)=c Y(\theta,\phi)
$$
This is the equation of the spherical harmonics. They coninstrain the $m_{\ell}$ with $\ell$. Physically is because $\ell$ determine the $L^2$ and $m_\ell$ one component of the vector $L$.