Figure 1: The spherical coordinates.
In the spherical coordinates, the Laplace equation reads: \[\begin{eqnarray} \nabla ^2\psi = \left[ \frac{1}{r^2}\frac{\partial }{\partial r}\left( {r^2\frac{\partial }{\partial r}} \right)+ \frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left( {\sin \theta \frac{\partial }{\partial \theta }} \right) +{ \frac{1}{r^2\sin ^2\theta }\frac{\partial ^2 }{\partial \phi ^2}} \right]\psi \tag{1}. \end{eqnarray}\]
We can separate the variables as \(\psi(r,\theta,\phi)=R(r)\Theta(\theta)\Phi(\phi)\): \[\begin{eqnarray} \frac{\Theta\Phi}{r^2}\frac{d}{d r}\left(r^2\frac{d R }{d r}\right) +\frac{R\Phi}{r^2\sin\theta}\frac{d }{d\theta}\left(\sin\theta \frac{ d \Theta}{d\theta}\right)+ \frac{ R\Theta}{r^2\sin ^2\theta }\frac{\partial ^2 \Phi}{\partial \phi ^2}=0 \tag{2}, \end{eqnarray}\] or equivalently \[\begin{eqnarray} \frac{\sin^2\theta}{R}\frac{d}{d r}\left(r^2\frac{d R }{d r}\right) +\frac{\sin\theta}{\Theta}\frac{d }{d\theta}\left(\sin\theta \frac{ d \Theta}{d\theta}\right)=- \frac{ 1}{\Phi }\frac{\partial ^2 \Phi}{\partial \phi ^2}. \tag{3} \end{eqnarray}\] Since the left-hand side of Eq. (3) depends on \(r\) and \(\theta\) only, and the right one depends on \(\phi\) only, overall they can only be equal to a constant, which we will call \(m^2\). This separates out the \(\Phi\) function. Furthermore, since \(\phi\) is the angle, the solutions have to be \(2\pi\) periodic, which gives: \[\begin{eqnarray} \Phi(\phi)=e^{i m\phi} \tag{4}. \end{eqnarray}\] Putting this back in Eq. (3) and dividing the it by \(\sin^2\theta\) we get: \[\begin{eqnarray} \frac{1}{R}\frac{d}{d r}\left(r^2\frac{d R }{d r}\right)= -\frac{1}{\sin\theta \Theta}\frac{d }{d\theta}\left(\sin\theta \frac{ d \Theta}{d\theta}\right) -\frac{m^2}{\sin^2\theta} \tag{5}. \end{eqnarray}\] Similarly, since the left-hand side of Eq. (5) depends on \(r\) only, and the right one depends on \(\theta\) only, overall they can only be equal to a constant, which we will call \(c\).
Radial dependence
The form of the solution for \(R(r)\) is easy to guess since the derivatives are balanced by the powers of \(r\), and therefore, a function of the form \(r^l\) will preserve its form up to a coefficient.
\[\begin{eqnarray} \frac{d}{d r}\left(r^2\frac{d r^l }{d r}\right)-c r^l=\left[ l(l+1)-c\right]r^l=0 \implies c=l(l+1). \tag{6} \end{eqnarray}\] However, notice the unexpected symmetry of \(l(l+1)\) under \(l\to-l-1\). This means, if \(r^l\) is a solution, so is \(r^{-l-1}\). This suggests the following form of solution for \(R\): \[\begin{eqnarray} R(r)=a_l r^l +\frac{b_l}{r^{l+1}}. \tag{7} \end{eqnarray}\]
Angular part
Putting this back in Eq. (5) yields \[\begin{eqnarray} \frac{1}{\sin\theta}\frac{d }{d\theta}\left(\sin\theta \frac{ d \Theta}{d\theta}\right) +\left[l(l+1)-\frac{m^2}{\sin^2\theta} \right]\Theta =0. \tag{8} \end{eqnarray}\] Now define \(\cos\theta=x\), which gives \(\frac{d}{d\theta}=\frac{dx}{d\theta}\frac{d}{dx}=-\sin\theta \frac{d}{dx}\) and insert this back in Eq.(8): \[\begin{eqnarray} (1-x^2)\frac{d^2 \Theta}{dx^2}-2x\frac{d\Theta}{dx} +\left[l(l+1)-\frac{m^2}{1-x^2} \right]\Theta =0. \tag{9} \end{eqnarray}\] Let’s first attempt to solve this for \(m=0\) using power series expansion[1]:
\[\begin{eqnarray} \Theta(x)=\sum_{k=0}^\infty c_k x^k. \tag{10} \end{eqnarray}\] Inserting this back to Eq. (9) we get:
\[\begin{eqnarray} 0&=&\sum_{k=0}^\infty c_k x^{k-2}k(k-1)-\sum_{k=0}^\infty c_k x^{k}k(k-1)+\sum_{k=0}^\infty c_k x^k \left[l(l+1)-2k\right]\nonumber\\ &=&\sum_{k=2}^\infty c_k x^{k-2}k(k-1)-\sum_{k=0}^\infty c_k x^{k}k(k-1)+\sum_{k=0}^\infty c_k x^k \left[l(l+1)-2k\right]\nonumber\\ &=& \sum_{k=0}^\infty c_{k+2} x^{k}(k+2)(k+1)+\sum_{k=0}^\infty c_k x^k \left[l(l+1)-k(k+1)\right]\nonumber\\ &=& \sum_{k=0}^\infty \left\{ c_{k+2} (k+2)(k+1)- c_k \left[k(k+1)-l(l+1)\right]\right\} x^{k} . \tag{11} \end{eqnarray}\] This implies \[\begin{eqnarray} c_{k+2}=\frac{k(k+1)-l(l+1)}{ (k+2)(k+1)} c_k =\frac{(k-l)(l+k+1)}{ (k+2)(k+1)} c_k , \tag{12} \end{eqnarray}\] which is the recurrence equation for the expansion coefficients.
This is a remarkable equation since it has profound consequences. Earlier in this blog, we looked at the Quantum Harmonic Oscillator and showed that for a similar series expansion to converge, we had to have the energy quantized. In this particular problem, until this point, we have no indication of \(l\) being an integer. But now, we see that it has to be an integer so that the series truncates for \(k>l\) (for every other \(k\)). That’s the first observation.
The second observation is associated with the parity symmetry of the original differential Eq. (9), which is invariant under \(x\to-x\) upto the overall sign. This shows that the solutions will also be eigenstates of the parity operator, i.e., odd and even \(k\) terms should not mix.
We have a couple of ways of terminating the series. The first one is what we have discussed above, i.e., settin \(l\) to an integer \(k^*\), which will zero out every other \(c_k\). The \(c_k\)’s, with \(k>l\), not addressed by this truncation need to be eliminated directly by their root coefficient, \(c_0\) or \(c_1\). To be more specific, take an example \(l=1\). The \(c_k\)’s with odd \(k\) quickly terminate after \(k=1\): \(c_1,0,\cdots\). The even ones will keep growing: \(c_0, \alpha c_0,\beta\alpha c_0,\cdots\). The only way to tame this series is by killing it at its root, i.e., by setting it \(a_0=0\) so that all the even terms drop out. This shows that even odd powers of \(x\) will not mix preserving respecting the parity symmetry of the original equation.
From Eq. (12), we can explicitly write the fist few Legendre polynomials: \[\begin{eqnarray} P_0(x)&=&1,\nonumber\\ P_1(x)&=&x,\nonumber\\ P_2(x)&=&\frac{1}{2}(3x^2-1),\nonumber\\ P_3(x)&=&\frac{1}{2}(5x^3-3x),\nonumber\\ P_4(x)&=&\frac{1}{8}(35x^4-30x^2+3). \tag{13} \end{eqnarray}\]
Having shown that Legendre polynomials solve the differential equation (with \(m=0\)), \[\begin{eqnarray} (1-x^2)\frac{d^2 P_l}{dx^2}-2x\frac{dP_l}{dx} +l(l+1)P_l =0. \tag{14} \end{eqnarray}\] We now need to address the full equation with \(m\neq 0\). The idea would be to differentiate \(m\) times to create the \(m^2\) term. For this, we will need the Leibniz’s formula: \[\begin{eqnarray} \frac{d^m}{dx^m}\left[f(x)g(x)\right]=\sum_{k=0}^n {n \choose k} \frac{d^k f}{dx^k}\frac{d^{n-k}g}{dx^{n-k}} \tag{15} \end{eqnarray}\] Let’s dive into the differentiation: \[\begin{eqnarray} 0&=&\frac{d^m}{dx^m}\left[(1-x^2)\frac{d^2 P_l}{dx^2}-2x\frac{dP_l}{dx} +l(l+1) P_l\right]\nonumber\\ &=&(1-x^2)\frac{d^m}{dx^m}\frac{d^2 P_l}{dx^2}+m\frac{d}{dx}(1-x^2)\frac{d^{m-1}}{dx^{m-1 }}\frac{d^2 P_l}{dx^2} +\frac{m(m-1)}{2}\frac{d^2}{dx^2}(1-x^2)\frac{d^{m-2}}{dx^{m-2 }}\frac{d^2 P_l}{dx^2} \nonumber\\ &&-2x\frac{d^m}{dx^m}\frac{dP_l}{dx}-2m\frac{d}{dx}( x) \frac{d^{m-1} }{dx^{m-1} }\frac{dP_l}{dx} +l(l+1) \frac{d^m}{dx^m} P_l\nonumber\\ &=&(1-x^2)u''-2mx u' -m(m-1)u-2xu' -2mu +l(l+1)u\nonumber\\ &=&(1-x^2)u''-2(m+1)x u' - m(m+1)u +l(l+1)u\nonumber\\ &=&(1-x^2)u''-2(m+1)x u' - (l-m)(l+m+1)u \tag{16}, \end{eqnarray}\] where \(u\equiv\frac{d^m P_l}{dx^m}\). We still need to modify the equation further so that it matches Eq.(9). First of all, we note that the equation we want to get at was self-adjoint, and we kind of destroyed it as we acted with \(\frac{d^m}{dx^m}\). Let’s restore it and see where it takes us.
Sturm–Liouville theory
We are going to use some machinery from Sturm–Liouville theory on second order differential equations. Consider the second order differential operator \(\mathcal{L}\)[1]:
\[\begin{eqnarray} \mathcal{L} u=\left(p_0(x)\frac{d^2 }{dx^2}+p_1(x) \frac{d }{dx} +p_2(x)\right)u \tag{17}. \end{eqnarray}\] We are going to define the inner product in the function space as an integral in a range \([a,b]\). \[\begin{eqnarray} \langle u| \mathcal{L}| u\rangle\equiv \int_a^b dx u(x) \mathcal{L} u(x)= \int_a^b dx u(x)\left(p_0\frac{d^2 }{dx^2}+p_1 \frac{d }{dx} +p_2\right)u= \int_a^b dx \left( p_0 u u'' + p_1 u u' +u^2 p_2\right) \tag{18}. \end{eqnarray}\] We can integrate Eq.(18) by parts. Let’s look at each term one by one:
\[\begin{eqnarray} u p_0 u'' &=& u \frac{d^2 }{dx^2}(p_0 u)-\frac{d}{dx}(up'_0 u)+2 u p'_0u' \nonumber\\ u p_1 u' &=& - u \frac{d }{dx}(p_1 u)+\frac{d}{dx}(u p_1 u) \tag{19}. \end{eqnarray}\] Putting this back in Eq.(18) gives:
\[\begin{eqnarray} \langle u| \mathcal{L}| u\rangle &=& \left[u(p_1-p_0')u \right]_a^b+ \int_a^b dx u\left(\frac{d^2 }{dx^2}(p_0 u) -\frac{d }{dx}(p_1 u) +p_2 u\right)\nonumber\\ &=&\left[u(p_1-p_0')u \right]_a^b+ \int_a^b dx u \bar{\mathcal{L}}u, \tag{20} \end{eqnarray}\] where the adjoint operator \(\bar{\mathcal{L}}\) is defined as: \[\begin{eqnarray} \bar{\mathcal{L}}u=\frac{d^2 }{dx^2}(p_0 u) -\frac{d }{dx}(p_1 u) +p_2 u. \tag{21} \end{eqnarray}\] Although \(\bar{\mathcal{L}}u\) looks pretty different from \({\mathcal{L}}u\) in Eq. (17), they can actually be the same if \(p_1=p'_0\). Such \(\mathcal{L}\) operators are self-adjoint. Furthermore, note that the boundary term also drops out for self-adjoint operators.
The good news is that if an equation is not self adjoint, it can be converted into that form if it gets multiplied by the following factor: \[\begin{eqnarray} \frac{1}{p_0(x)}\text{exp}\left\{\int^x dt \frac{p_1(t)}{p_0(t)}\right\}. \tag{22} \end{eqnarray}\]
Let’s revisit Eq. (16) to find the factor that will make the equation self-adjoint: \[\begin{eqnarray} \frac{1}{1-x^2}\text{exp}\left\{-\int^x dt \frac{2(m+1)t}{1-t^2}\right\}= \frac{1}{1-x^2}\text{exp}\left\{(m+1)\int^x \frac{d(1-t^2)}{1-t^2}\right\} =(1-x^2)^{m}. \tag{23} \end{eqnarray}\] We will take this factor, multiply Eq. (16) with it to get:
\[\begin{eqnarray} \frac{d}{dx}\left((1-x^2)^{m+1}u'\right) - (l-m)(l+m+1)u=0 \tag{24}. \end{eqnarray}\] Finally, we will want to absord half power of that coefficient into \(u\) by defining \[\begin{eqnarray} v(x)=(1-x^2)^{\frac{m}{2}} u(x) \iff u(x)=(1-x^2)^{-\frac{m}{2}} v(x) \tag{25}, \end{eqnarray}\] to get \[\begin{eqnarray} u'&=&\left[ v' +\frac{m v x}{1-x^2}\right](1-x^2)^{-\frac{m}{2}},\nonumber\\ u''&=&\left[ v'' +\frac{2m v' x}{1-x^2}+\frac{m x}{1-x^2}+\frac{m (m+2) x^2 v}{(1-x^2)^2} \right](1-x^2)^{-\frac{m}{2}} \tag{26}. \end{eqnarray}\] The new function \(v\) satisfies the following equation:
\[\begin{eqnarray} (1-x^2)\frac{d^2 v}{dx^2}-2x\frac{d v}{dx} +\left(l(l+1)-\frac{m^2}{1-x^2} \right) v =0 \tag{27}, \end{eqnarray}\] which is identical to Eq. (9). In conclusion, the angular part of the solution is given by the associated Legendre polynomials as below: \[\begin{eqnarray} P_l^m(x)=(1-x^2)^{\frac{m}{2}}\frac{d^m }{dx^m} P_l(x) \tag{28}. \end{eqnarray}\]
Full solution
Note that the highest power in \(P_l\) is \(l\), and for \(m>l\), we run out of \(x\)’s to differentiate. This automatically limits \(|m|\) to \(l\). Putting all pieces together, the full solution to the Laplace equation reads:
\[\begin{eqnarray} \psi(r,\theta,\phi)=R(r)\Theta(\theta)\Phi(\phi)= \sum_{l=0}^\infty \sum_{m=-l}^l \left(a_{lm} r^l +\frac{b_{lm}}{r^{l+1}}\right) P_l^m(\cos\theta)e^{im \phi}. \tag{29} \end{eqnarray}\]