This question (similar other questions) have troubled me for last 5 years. Two months ago I thought I would simply solve the Hydrogen Atom problem and see these Associated Legendre Polynomials come to life by themselves. I haven't solved it yet.
I began by writing the Schrodinger equation for Hydrogen atom.
$$
-\frac{\hbar^2}{2m} \nabla^2 \psi + V \psi = E\psi
$$
I followed a well know method of separation of variables.
$\psi(r,\theta,\phi) = R(r)P(\theta)Q(\phi)$. I plugged it in the above Schrodinger equation. The result was,
$$
\underset{\text{function of r}}{\underbrace{\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) + \frac{2mr^2}{\hbar^2}(E-V)}}
+ \underset{\text{function of x}}{\underbrace{\frac{1}{P}\frac{d}{dx}\left((1-x^2)\frac{dP}{dx}\right)
+ \frac{1}{(1-x^2)}\underset{\text{function of }\phi}{\underbrace{\frac{1}{Q}\frac{d^2Q}{d\phi^2}}}}}
= 0
$$
Here $x = \cos\theta$. Now we can proceed step by step to solve it.
Step 1. $\frac{1}{Q}\frac{d^2Q}{d\phi^2} = -m^2$ where m is an integer (due to the periodic nature of the $Q(\phi)$).
Step 2. This gives us the two following equations,
$$
\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) + \frac{2mr^2}{\hbar^2}(E-V) = A
$$
$$
\frac{1}{P}\frac{d}{dx}\left((1-x^2)\frac{dP}{dx}\right)
– \frac{m^2}{(1-x^2)} = -A
$$
where A is some constant. We have not got a clue about this constant at all this point.
I simply am interested in associated Legendre so I went for the later equation.
Step 3.
$$
(1-x^2)^2P'' -2x(1-x^2)P' +\left(A(1-x^2) – m^2\right)P = 0
$$
Now this equation can be solved by assuming power series solution of the form,
$P = \sum_{j=0}^{\infty} a_j x^j$.
Step 4.
When this expansion is plugged into the the Associated Legendre DE I got,
$(j+2)(j+1)a_{j+2}+[A-m^2-2j^2]a_j+[(j-2)(j-3)-A]a_{j-2} = 0$
Remembering that $a_j = 0\ \forall\ j<0$ I found few terms,
$$a_2 = -\frac{(A-m^2)}{2\cdot1}a_0$$
$$a_3 = -\frac{(A-m^2-2)}{3\cdot2}a_1$$
$$a_4 = \left[A + \frac{(A-m^2)(A-m^2-8)}{2\cdot1}\right]\frac{a_0}{4\cdot3}$$
$$a_5 = \left[A + \frac{(A-m^2-2)(A-m^2-18)}{3\cdot2}\right]\frac{a_1}{5\cdot4}$$
Now I don't see any pattern emerging through which I could place condition on the constant A. Kindly someone help me this question is driving me crazy.
Best Answer
The Legendre polynomials $P_n(x)$ are the starting point for most developments. These polynomials were defined in the late 1700's by Legendre when performing potential expansions in $\mathbb{R}^3$: $$ \frac{1}{|x-x_1|}=P_0(\cos\theta)\frac{1}{|x_1|}+P_1(\cos\theta)\frac{|x|}{|x_1|^2}+P_2(\cos\theta)\frac{|x|^2}{|x_1|^3}+\cdots. $$ In Legendre's expansion, $\theta$ is the angle between the vectors at the origin. Legendre was able to show that such an expansion would involve $P_n$, which would be an $n$-th order polynomial. Everything can be reduced to the case where $x_1$ is a unit vector along the $z$-axis. Then, using spherical coordinates, with $\theta$ being the angle from the $z$ axis, $$ \frac{1}{|x-\hat{z}|}=\frac{1}{\sqrt{1-2r\cos\theta+r^2}} = \sum_{n=0}^{\infty}r^{n}P_n(\cos\theta). $$ The Laplacian in spherical coordinates is $$ \nabla^2= \frac{1}{r^2}\left[\frac{\partial}{\partial r}r^2\frac{\partial}{\partial r}+\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}\right] $$ Plugging the gravitational potential into $\nabla^2 \frac{1}{|x-\hat{z}|}$ then led to differential equations for the Legendre polynomials $$ \sum_{n=0}^{\infty}n(n+1)r^nP_n(\cos\theta)+\frac{1}{\sin\theta}\frac{d}{d\theta}\sin\theta\frac{d}{d\theta}P_{n}(\cos\theta) = 0. $$ The substitution $x=\cos\theta$ then leads to $$ \frac{d}{d\theta}=-\sin\theta\frac{d}{dx} \\ \sin\theta\frac{d}{d\theta}=-\sin^2\theta\frac{d}{dx}=-(1-x^2)\frac{d}{dx} \\ \frac{1}{\sin\theta}\frac{d}{d\theta}\sin\theta\frac{d}{d\theta}=\frac{d}{dx}(1-x^2)\frac{d}{dx} $$ Therefore, the Legendre polynomials $P_n(x)$ satisfy $$ \frac{d}{dx}(1-x^2)\frac{d}{dx}P_n(x)+n(n+1)P_n(x)=0,\;\;\; n=0,1,2,3,\cdots. $$ In 1816 Olinde Rodrigues discovered an elegant way to write these polynomials: $$ P_n(x) = \frac{1}{2^nn!}\frac{d^n}{dx^n}(x^2-1)^n. $$ So you have the eigenvalues $n(n+1)$ and eigenfunctions that form an orthogonal system. Because the polynomials are dense in $L^2(-1,1)$, these polynomials form a complete orthogonal basis of $L^2(-1,1)$. The normalization constants required for the $P_n$ are not difficult to compute using the Rodrigues formula and integration by parts.
By playing with the generating function $$ \frac{1}{\sqrt{1-2rx+r^2}} = \sum_{n=0}^{\infty}r^nP_n(x), $$ the higher order Legendre functions were discovered. For example, differentiating $m$ times in $x$ gives $$ \frac{(1)(3)(5)(\cdots)(2m+1)r^{m}}{(1-2rx+r^2)^{m/2}}=\sum_{n=m}^{\infty}r^{n}P_n^{(m)}(x). $$ The polynomials $P_n^{(m)}$ are the associated Legendre polynomials. If you start with the Legendre equation and differentiate it $m$ times, you end up with a new differential equation for the functions $y=P_n^{(m)}$ of the form $$ (1-x^2)y''-2(m+1)xy'+(n(n+1)-m(m+1))y = 0. $$ It is this form of the associated Legendre equation where you can get a series solution. And this equation is equivalent to the equation you want to study. This equation is no longer in Sturm-Liouville form. However, the substitution $$ y = (1-x^2)^{-m/2}g $$ puts the equation into standard form with respect to $g$: $$ \frac{d}{dx}\left((1-x^2)\frac{dg}{dx}\right)-\frac{m^2}{1-x^2}g+n(n+1)g=0. $$ In this way, the standard solutions of the associated Legendre equation were found to be $$ P_n^m(x) = (-1)^m\frac{1}{2^nn!}(1-x^2)^{m/2}\frac{d^{n+m}}{dx^{n+m}}(x^2-1)^{n}. $$ So, quite remarkably, your constant $A$ is found to not depend on $m$; A depends only on $n$ because $$ \left[-\frac{d}{dx}(1-x^2)\frac{d}{dx}+\frac{m^2}{1-x^2}\right]P_{n}^{m}(x)=n(n+1)P_{n}^{m}(x),\;\;\; n=m,m+1,m+2,\cdots. $$ (Note that for $n < m$, the functions $P_n^m$ are $0$.) So the values of your separation constant $A$ are $n(n+1)$ for $n=0,1,2,3,\cdots$, but the only valid values for a given $m$ are $n(n+1)$ for $n=m,m+1,m+2,\cdots$. For a fixed $n$, there are non-zero solutions for $m=0,1,2,3\ldots,n-1$. And is this structure that imposes itself on the orbits of the electron. So, finally, your radial equation is known. $A=n(n+1)$ is required, which comes out of $\frac{d}{dr}\left(r^2\frac{d}{dr}\right)r^n=n(n+1)r^n$.
This analysis was firmly established by the mid to late part of the 19th century, well before Quantum Mechanics. Many attempts have been made to replace this analysis, and to make it more intuitive or abstract, but it seems to me that nothing beats the original story or the closed-form solutions. Quantum ladder operators are really nice, though, I must admit. Of course the radial equation for the Hydrogen isotope was not studied in that time because no standard equations would have led to it. That's why the radial equation has no name attached to it, other than the radial Hydrogen equation.