Associated Legendre Polynomial - Part - 1

The General Legendre differential equation as we know, $$ (1-x^2)y''-2xy'+n(n+1)y=0$$ is solved using the frobenius power series method and the solutions are obtained to be the general Legendre polynomials, where suitable normalization constants and boundary conditions are used. 
With this result, we proceed further to completely exploit all the possibile solutions can be obtained from this equation. 
Before that, we will have to use the orthonormalizability property of our solution which is proved as, 
From the generating function, $$ \frac{1}{\sqrt{1-2xz+z^2}}=\sum_n P_n(x) z^n$$ or $$ \frac{1}{1-2xz+z^2} = \sum_n \sum_m z^{n+m} P_n(x) P_m(x) \\ \int\frac{1}{1-2xz+z^2} \,dx = \sum_{m,n} z^{n+m} \int{P_n(x)}P_m(x) \,dx$$ The limits can change maximum from -1 to +1 as in Legendre polynomials and summation is usually implied from 0 to infinity. 
Left hand side term is evaluated to give,$$ \frac{-1}{2z}\int_{(1+z)^2}^{(1-z)^2}\frac{1}{u} du = \frac{-1}{2z} \ln{\frac{(1-z)^2}{(1+z)^2}}\\where\,\, u=1-2xz+z^2 and \,\,du = -2z\,dx$$
and $$ \frac{-1}{2z}\ln{\left(\frac{1-z}{1+z}\right)^2} = \frac{-2}{2z} \ln{\frac{1-z}{1+z}} = \frac{1}{z}\ln{\frac{(1+z)}{(1-z)}}$$
$$\ln{\frac{(1+z)}{(1-z)}} = \ln(1+z) - \ln(1-z)$$
using taylor expansion, 
$$ f(x) = f(a) + f'(a)\, (x-a) + \frac{f''(a)}{2!}\,(x-a)^2+...$$
where the function should be infinitely differentiable at x=a.
So, $$ \ln(1+z) = \ln(1+a) + \frac{1}{1+a}(z-a) - \frac{1}{2!(1+a)^2} (z-a)^2 +...$$
where we can take a=0, $$ \ln(1+z) = z - \frac{z^2}{2} + \frac{z^3}{3} +... $$ Similarly, $$ \ln(1-z) = -z - \frac{z^2}{2}  - \frac{z^3}{3} -..$$ This gives, $$ \ln(1+z) - \ln(1-z) = 2 \left[z+\frac{z^3}{3} + \frac{z^5}{5}+...\right] = 2\sum_{n} \frac{z^{2n+1}}{2n+1}$$
Thus we have our integral equal to, $$ \sum_{m,n} z^{n+m} \int_{-1}^{1} P_n(x)\,P_m(x) \,dx = \frac{1}{z} \sum_{n} \frac{z^{2n+1}}{2n+1} = 2  \sum_n \frac{z^{2n}}{2n+1}$$
From this, the left side terms will equal the right side terms only when n=m, so we get the result,
$$ when\,\, n=m\rightarrow\,\, \sum_n z^{2n} \int_{-1}^{+1}P_n(x)\,P_m(x) \,dx = 2\sum_n \frac{z^{2n}}{2n+1}\\ \rightarrow\,\,\,\, \int_{-1}^{1} P_n(x)\,P_m(x)\, dx = \frac{2}{2n+1}$$ $$ when \,\, n\neq{m}\rightarrow\,\, \int_{-1}^{1}P_n(x)\,P_m(x) \,dx = 0 \\ or \int_{-1}^{1}P_n(x)\,P_m(x)\,dx = \delta_{nm} \frac{2}{2n+1}$$ which is our orthogonality relation. 
We can also prove the completeness relation - that any arbitrary function can be expanded in terms of the linear combination of these Legendre polynomials. 
$$ f(x) = \sum_n a_n P_n(x)\\ \int_{-1}^{1} f(x)\,P_m(x)\,dx = \sum_n a_n \frac{2}{2n+1}\delta{nm} = \sum_n a_n \int_{-1}^{1}P_n(x)\,P_m(x)\,dx = \frac{2}{2m+1}a_m $$ amd $$ a_m = \frac{2m+1}{2} \int_{-1}^{1} f(x)\,P_m(x)\,dx$$
 Now, we will get back to our differential equation where we substitute our solution and differentiate it "m" times w.r.t. x 
(m < n), $$ (1-x^2)P_n''-2xP_n'+ n(n+1)P_n=0$$
$$ \frac{d^m}{dx^m}\left((1-x^2)\frac{d^2P_n}{dx^2}\right) = (1-x^2) \frac{d^{m+2}P_n}{dx^{m+2}} - m.2x.\frac{d^{m+1}P_n}{dx^{m+1}} - \frac{m(m-1)}{2!}2 \frac{d^mP_n}{dx^m}$$ 
where we made use of the Leibniz formula, $$ \frac{d^m}{dx^m} A(x) \,B(x) = \sum_{k=0}^m \frac{m!}{k!(m-k)!} \frac{d^kA}{dx^k} \frac{d^{n-k}B}{dx^{n-k}}$$
Similarly, $$ -2\frac{d^m}{dx^m}\left[x\frac{dP_n}{dx}\right] = -2x \frac{d^{m+1}P_n}{dx^{m+1}} - 2m\frac{d^mP_n}{dx^m}$$
Which finally gives, $$ \frac{d^m}{dx^m}\left[Legendre\, eqn.\right] \\= (1-x^2)\frac{d^2}{dx^2}\left(\frac{d^mP_n}{dx^m}\right) - 2(m+1)x \frac{d}{dx}\left(\frac{d^mP_n}{dx^m}\right) + \frac{d^mP_n}{dx^m}\left[n(n+1)-m(m+1)\right] = 0 $$
If we assume, $$ \frac{d^mP_n}{dx^m} = V = \frac{W}{(1-x^2)^{\frac{m}{2}}}$$ So, $$ \frac{dV}{dx} = \frac{1}{(1-x^2)^{\frac{m}{2}}}\left[ \frac{dW}{dx} +\frac{mxW}{1-x^2}\right]$$ and $$ \frac{d^2V}{dx^2} = \frac{1}{(1-x^2)^{\frac{m}{2}}}\left[ \frac{d^2W}{dx^2} + \frac{2mx}{1-x^2}\frac{dW}{dx}+W\left(\frac{m}{1-x^2}+ \frac{mx^2(m+2)}{(1-x^2)^2}\right)\right]$$ Substituting and doing the necessary manipulation (multiply the whole equation by $(1-x^2)^{\frac{m}{2}}$, We get our associated Legendre equation, $$ (1-x^2)\frac{d^2W}{dx^2} - 2x\frac{dW}{dx} + \left[n(n+1) - \frac{m^2}{1-x^2}\right]W = 0 $$

Eigen values of $L^2$ and $L_z$ operator

The Angular momentum operator is defined by, $$ \vec{L} = \vec{r}\times\vec{p}= -i\hbar\vec{r}\times\nabla $$ In cartesian coordinates we have the usual known coordinates, $$ {L_x} = y{p_z}-z{p_y}\\ {L_y} = zp_x-xp_z\\ L_z = xp_y-yp_x$$ where in quantum mechanics, these terms are replaced by corresponding operators so we get the commutation relations as, $$ [L_x,L_y] = i\hbar{L_z} $$ with other components in cyclic form we have, $$ [L_i,L_j] = i\hbar\epsilon_{ijk}L_k$$ using levi-civita symbol. The components of the angular momentum operators are real, Hermitian since the terms on each coordinate commute with one another. And so, the angular momentum operator itself is a Hermitian operator. 
We can also express it in spherical polar coordinates as given in 
http://scientistech.blogspot.it/2015/11/angular-momentum-operator-in-spherical.html 

From, $$[r_j,p_k] = \delta_{jk}i\hbar$$ since these are operators, Identity operator is intrinsically assumed here. 
Using these we find $L^2$ operator commutes with every other component. 
$$ [L^2, L_x] =  [L_x^2 , L_x] + [L_y^2,L_x] + [L_z^2,L_x] $$
where, $$ [L_x^2,L_x] = 0 $$ $$[L_y^2,L_x] = [L_y,L_x]L_y + L_y[L_y,L_x] = - i\hbar{L_yL_z+L_zL_y}$$
$$ [ L_z^2,L_x] = [L_z,L_x]L_z + L_z[L_z,L_x] = i\hbar{L_yL_z+L_zL_y}$$ Combining all three we get, $$ [L^2, L_x] = 0$$ and similarly for other components. And so, $$ [L^2, \vec{L}]=0$$ 
So, $L^2$ and $L_z$ can be simultaneously diagonalized. It just a convention we use $L_z$. All that matters is that we can simultaneously diagonalize the $L^2$ with one of its components but not with all the three as we know other two are non-commuting with each other. 
Now, we define two new operators as. $$ L_+ = L_x + iL_y \\ L_- = L_x - iL_y$$ where one operator is the adjoint of the other. 
Once you define an operator, you can check for its commutation relation with all other operators, such as, $$[L_+, L_-] = [L_x+iL_y, L_x-iL_y] = 2\hbar{L_z}$$ and $$ [L_z,L_+] = [L_z,L_x+iL_y] = i\hbar{L_y}+ \hbar{L_x} = \hbar{L_+}$$ and $$[L_z, L_-] = -\hbar{L_-}$$ and finally, $$L_+L_- = L_x^2 + L_y^2 +iL_yL_x-iL_xL_y = L_x^2+L_y^2+\hbar{L_z} \\ \rightarrow L^2 = L_+L_- - \hbar{L_z} + L_z^2 $$ using the commutation between the two ladder operators, $$L^2 = L_-L_++ \hbar{L_z} +L_z^2$$

Remember, we have so far not used anything physical in our arguments. It is just that we define an operator as $$ \vec{L} = \vec{r}\times\vec{p} $$ and use the commutation relation $$[r_j,p_k] = i\hbar\delta_{jk}$$ for manipulation. Other than this, it can be or mean anything. It needn't even have any physical correspondence at all. 

Now, we make use of the commutation relation between $L^2,\,and\,L_z$ and seek simultaneous eigenstates[1] for both of these operator (a physical statement from the Heisenberg's principle[2]). 
We start by defining the eigenstate as, $$ L^2|X\rangle = \lambda\hbar^2|X\rangle \\ L_z|X\rangle = m\hbar|X\rangle $$
where $\hbar $ and its square is chosen for dimensional reasons. 
Note: We are just talking about the simultaneous eigenstates and we never checked whether it is the complete set of solution or not?[1].
Now, we check for the action of $L_+$ operator on any eigenstate by making use of the commutation, $$ [L_z,L_+] = L_+\hbar \\ \rightarrow L_zL_+ = L_+L_z + \hbar{L_+}$$ and it gives, $$ L_zL_+|X\rangle = \left(L_+L_z +\hbar{L_+}\right)|X\rangle = \hbar (m+1) L_+|X\rangle$$
similarly, $$L_zL_-|X\rangle = \hbar(m-1)L_-|X\rangle$$ From the above, it is understood that if $L_z$ operator acting on the general $|X\rangle$ gives the eigen value $m\hbar$. But the same operator acting on the wavefunction produced by $L_+ \,\,and\,\,L_-$ operator gives the new eigenvalue which is respectively unit value higher and lower than the old eigenvalue. i.e.These operators works as step operator which makes the eigenfunction to jump from one state to another in fixed unit values. From,$$ L^2 = L_x^2 +L_y^2 +L_z^2$$  we see that the for a fixed "$\lambda$" value there is a maximum and minimum value for "m" since the eigenvalues of Hermitian operators are real and its square is positive quantity $$ \lambda\geq{m^2}$$
Thus, there is a maximum and minimum value for "m" such that, $$ L_+|X_{max}\rangle = 0 \\ L_-|X_{min}\rangle = 0 $$ 
and the corresponding eigenvalue of $L^2$ operator is given by, $$ L^2 = \left(L_-L_++\hbar{L_z} + L_z^2\right)|X_{max}\rangle = \hbar^2m_{max}(m_{max}+1) $$
Similarly for the lowest state , $$ L^2 |X_{min}\rangle= \left(L_+L_- - \hbar{L_z} + L_z^2\right)|X_{min}\rangle = \hbar^2 m_{min}(m_{min}-1)|X_{min}\rangle$$
But $L^2$ operator should have the same eigenvalue and so $$ m_{min}(m_{min}-1) = m_{max}(m_{max}+1) \\ \rightarrow\,\,\, m_{min} = - m_{max}=-l \,\,and\,\, m_{min} = m_{max}+1$$ where the second possibility is excluded from physical reality. So, we do have, "m" varies from -l to +l in integer steps which forces the values of "l" to be the integers and the values of "l" are 0, 1/2, 1, 3/2,2..etc.(either positive or negative).   
Thus we can exploit the properties of the eigenvalues of $L^2$ and $L_z$ operators. 

Check for :
1.Simultaneous eigenstates
2.Heisenberg's Uncertainty principle - general statement