How to derive the power series solution for Associated Legendre Differential Equation

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.