# Laplace's equation and separation of variables in spherical coordinates done rigorously

Consider a region $$\Omega \subset \mathbb{R}^3$$. A well-known method of solving the Laplace equation $$\Delta f = 0$$ on $$\Omega$$ (we're looking for clasical, i.e. $$C^2$$ solutions) is to employ the separation of variables, i.e. to look for functions $$f: \mathbb{R}^3 \setminus \{ 0 \} \to \mathbb{R}$$ of the form $$f(r, \theta, \varphi) = R(r)\Theta(\theta)\Phi(\varphi)$$ and then claim that a general solution can be expressed by a series of such ansatz solutions. However, most references I found (Jackson, Griffiths and various resources found online) skip the details of this approach. I'd like to find a reference that treats this method rigorously

Precise description of the details I'm missing below.

Since the Laplace equation in spherical coordinates reads $$\nabla^2 f = \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 \frac{\partial f}{\partial r}\right) + \frac{1}{r^2 \sin\theta} \frac{\partial}{\partial \theta}\left(\sin\theta \frac{\partial f}{\partial \theta}\right) + \frac{1}{r^2 \sin^2\theta} \frac{\partial^2 f}{\partial \varphi^2} = 0$$

Substituting $$\psi(r, \theta, \varphi) = R(r)\Theta(\theta)\Phi(\varphi)$$ we get $$\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) + \frac{1}{\Theta}\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta \frac{\partial \Theta}{\partial\theta}\right) + \frac{1}{\Phi}\frac{1}{\sin^2\theta}\frac{\partial^2\Phi}{\partial\varphi^2} = 0$$

After multiplying by $$\sin^2\theta$$ we immediately get that for some constant $$\lambda$$ $$\frac{1}{\Phi}\frac{\partial^2\Phi}{\partial\varphi^2} = \lambda$$ Since $$\Phi$$ is continuous and $$2\pi$$-periodic, then necessarily $$\lambda < 0$$. At this point we're allowed to substitute $$\frac{1}{\Phi}\frac{\partial^2\Phi}{\partial\varphi^2} = -m^2$$ which yields solutions $$\Phi(\varphi) = \operatorname{Re} Ae^{\pm im \varphi}$$. For some reason most resources write $$\Phi(\varphi) = e^{\pm im \varphi}$$, but such a function is not a real function.

Similarly, we must have $$\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) = \alpha, \qquad \frac{1}{\Theta}\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta \frac{\partial \Theta}{\partial\theta}\right) = -\alpha$$

At this point one does a substitution $$\alpha = l(l+1)$$. However, this substitution is only if $$\alpha \geq - \frac {1}{4}$$. Moreover, the texts seem to assume that $$l$$ is an integer, and it's completely unclear why we can assume it. None of the texts I found really comments why this substitution is allowed.

Furthermore, $$Y_{lm}$$ is defined as $$Y_{lm}(\theta, \varphi) = c_{lm} P_{lm}(\cos \theta) e^{im\varphi}$$ However, we were looking for real solutions, and this solution is a complex one, the entity are mismatched.

Finally the general solution is claimed to be of the form: $$f(r, \theta, \varphi) = \sum_{l=0}^\infty \sum_{m=-l}^\ell a_{lm} r^l Y_{lm} (\theta, \varphi )$$ However, all the texts I found don't justify why functions of the form $$r^l Y_{lm}$$ form a complete set of functions. It's a priori unclear why functions $$r^l$$ form an orthogonal basis of $$L^2([0, \infty))$$.

• There are a lot of good questions here. I will post a detailed answer as soon as I can. Nov 12 '20 at 12:41
• Related Nov 12 '20 at 20:16
• Note that if you're just using this method for solving a particular BVP and you don't want to show that this method always works, you can just "guess" that $\alpha = l(l+1)$, find some solution and use the uniqueness of the solution. Nov 30 '20 at 21:14

There are a lot of great questions here, many of which I struggled with myself when I was first learning about the physical applications of PDEs. I'll go point by point and try to give the best explanations that I can. First, some housekeeping - whenever you use spherical coordinates, it is always good to explicitly state the coordinate transformation and notation you're using for reasons I mention in this answer and also discussed on Wolfram Mathworld. So, here we go:

I will use the $$\displaystyle ( r,\theta ,\phi )$$ (radial, azimuthal, polar) convention for spherical coordinates, as used by Wolfram Mathworld. Explicitly, the coordinate transformation is $$\begin{equation*} \begin{bmatrix} x\\ y\\ z \end{bmatrix} =\begin{bmatrix} r\cos \theta \sin \phi \\ r\sin \theta \sin \phi \\ r\cos \phi \end{bmatrix} \end{equation*}$$

Yes, I know the theta and phi are in the opposite order as to what you're using, but this is just the notation that I'm the most comfortable with and trying to switch it around will only make me more prone to errors. Let's get on with actually answering your questions then.

First question: Why do we have complex solutions to what should be a real valued PDE?

This is probably the easiest question to answer. Take for instance the ordinary differential equation $$y''+y=0$$ We can solve this in the complex plane using the usual methods as $$y(x)=c_1e^{-ix}+c_2e^{ix}$$ This is the most general solution and allows for complex valued $$y$$, and, critically, complex valued coefficients. Going a bit further, we can break $$c_1,c_2$$ and $$e^{-ix},e^{ix}$$ into their real and imaginary parts: $$y(x)=(a_1+ib_1)(\cos(x)-i\sin(x))+(a_2+ib_2)(\cos(x)+i\sin(x))$$ Which now allows us to split $$y$$ into its real and imaginary parts: $$y(x)=\operatorname{Re}y+\operatorname{Im}y=\bigg[(a_1+a_2)\cos(x)+(b_1-b_2)\sin(x)\bigg]+i\bigg[(a_2-a_1)\sin(x)+(b_1+b_2)\cos(x)\bigg]$$ If we make the artificial restriction that $$y$$ must be real valued, i.e $$\operatorname{Im}y=0$$, we get $$a_1=a_2~;~b_1=-b_2\implies c_1=c_2^*\implies y(x)=A\sin(x)+B\cos(x)$$. In short, if we were to write $$y''+y=0\implies y(x)=c_1 e^{-ix}+c_2e^{ix}$$ This is still valid even if we only want real valued $$y$$! It's just implied that $$c_1=c_2^*$$. Why do we do this? Sometimes, writing the solution in terms of exponentials can make expressions a little bit shorter and easier to manipulate, even if it means we have to introduce complex valued coefficients.

Second question: Why does everyone just assume $$l$$ is an integer? And what do we do if $$\alpha<-1/4$$?

The second part of this question is in some ways easier to answer than the first. the answer to the first question is, sort of, "because that's the simplest possible case", but there's more to it than that (which I'll get to after I answer the second part.)

In short, complex numbers come to save the day again. Suppose $$\alpha=-2$$. Then, $$l(l+1)=-2$$ $$l^2+l+2=0$$ $$l=\frac{-1\pm\sqrt{1-4\cdot2}}{2}=\frac{-1}{2}\pm i\frac{\sqrt{7}}{2}$$ This is totally allowed! While this will no longer be a Legendre polynomial, we do have what is known as the Legendre function (of the first kind): $$\mathscr{P}_{l,m}(z)=\frac{1}{\Gamma(1-m)}\left(\frac{1+z}{1-z}\right)^{m/2}{}_2F_1\left([-l,l];[1-m];\frac{1-z}{2}\right)$$ Which uses the (Gaussian) hypergeometric function. It can be shown (though not easily) that in the case of $$m,l\in \mathbb{N}$$ and $$m this reduces to the normal associated Legendre polynomial.

On to the first part then: why does everyone just assume $$l$$ is an integer? To answer this, let's solve the simpler problem of the spherical harmonics where we assume there is no $$\theta$$ (azimuthal) angle dependency. That is, let's have a look at the PDE $$\nabla^2 U(r,\phi)=\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial U}{\partial r}\right)+\frac{1}{r^2\sin\phi}\frac{\partial}{\partial \phi}\left(\sin\phi \frac{\partial U}{\partial \phi}\right)$$ Assuming separability, i.e $$U(r,\phi)=R(r)\Phi(\phi)$$, we get (also note that I've multiplied through by $$r^2$$ which is fine as long as we don't care about the origin) $$\frac{\partial}{\partial r}\left(r^2\frac{\partial R}{\partial r}\right)=\alpha R$$ And $$\frac{1}{\sin\phi}\frac{\partial}{\partial \phi}\left(\sin\phi \frac{\partial \Phi}{\partial \phi}\right)=-\alpha \Phi$$ A bit of an inspired substitution $$\alpha=l(l+1)$$. This substitution is particularly nice if $$\alpha$$ is a positive integer of the form $$n^2+n$$, so that $$l$$ is a positive integer, but in general we just have $$l=(-1\pm \sqrt{1-4\alpha})/2$$. Recognizing that the $$R$$ equation is a Cauchy-Euler ODE, we get $$R(r)=Ar^l+B\frac{1}{r^{l+1}}$$ Now for the $$\Phi$$ equation. We have, after some manipulation $$\Phi''+\cot(\phi)\Phi'+l(l+1)\Phi=0$$ Which, after a substitution $$\mu=\cos\phi$$, becomes (I'll let you work out the details) $$(1-\mu^2)\frac{\mathrm{d}^2\Phi}{\mathrm{d}\mu^2}-2\mu\frac{\mathrm{d}\Phi}{\mathrm{d}\mu}+l(l+1)\Phi=0$$ The reduced form of Legendre's differential equation.

So we can now resort to the Frobenius method and try to come up with a solution $$\Phi(\mu)=\sum_{n=0}^\infty a_n \mu^n$$. I'll cut to the chase (you can find the details on Wolfram Mathworld) and assert that after some work we get a recurrence $$a_{n}=\frac{\left(l+(n+1)\right)(l-n)}{(n+1)(n+2)}a_{n-2}$$ So what's the upshot here? Yes, you are correct: $$l$$ doesn't have to be an integer for the above to simply be a solution to our differential equation. BUT, If $$l$$ wasn't an integer, then the equation $$a_n=0$$ has no solutions since $$l+n+1=0$$ or $$l-n=0$$ is clearly impossible if $$l$$ is not an integer. So, as long as $$a_0\neq 0$$, none of the coefficients in the expansion of $$\Phi(\mu)$$ would zero, i.e, it would be an infinite sum with quickly growing coefficients and thus would diverge for $$\mu=1$$ which corresponds to $$\phi=0,\pi$$. Since typically we are solving for an electric potential, it doesn't make any physical sense to have a electrical potential that blows up near the poles. However if $$l$$ is an integer, we will at some point reach some value of $$n$$ such that $$a_n=0$$ which would lead to $$a_{n+2}=a_{n+4}=...=0$$. All the previous even coefficients will be nonzero and their precise form will follow from the recurrence relation. As for the odd coefficients, i.e $$a_{n+2m+1}$$, we can simply choose all of these to be $$0$$. Since the recurrence is of second order, we need only worry about every second coefficient and we can choose the odd ones to all be zero which trivially satisfies the recurrence. Which half of the coefficients we set to $$0$$ depends on the parity of $$l$$. The precise series expansions for each case are listed on Wolfram Mathworld. To illustrate the above, the expansions will look like $$P_l(\mu)=a_0\mu^0+0\mu^1 +a_2\mu^2+0\mu^3...+0\mu^{l-1}+a_l\mu^l+0\mu^{l+1}+0\mu^{l+2}+0+...\text{ for even }l$$ $$P_l(\mu)=0\mu^0+a_1\mu^1+0\mu^2+a_3\mu^3+...0\mu^{l-1}+a_l\mu^l+0\mu^{l+1}+0\mu^{l+2}+0+...\text{ for odd }l$$ As you can see, the expansions terminate which will mean that $$\Phi$$ will remain bounded at the poles.

To summarize: No, $$l$$ need not be an integer in order to get something that solves the ODE. But, in order to get a solution of the ODE that remains bounded, $$l$$ must be an integer.

Finally, onto complete function sets and orthogonal bases. It's pretty obvious that $$r^l$$ form an basis (but not an orthogonal basis) for $$L^2[0,\infty)$$ - simply consider Taylor expansions. However, the Legendre polynomials do provide an orthogonal basis for $$L^2[-1,1]$$. We can arrive at this by applying the Gram-Schmidt process to $$\{1,x,x^2,x^3,...\}$$ on the interval $$[-1,1]$$. (REFERENCE: https://en.wikipedia.org/wiki/Orthogonal_functions#Polynomials) If we instead want an orthogonal basis on $$[0,\infty)$$, we can perform a Cayley Transform which brings us to the Legendre rational functions. Finally, about the completeness of $$\{r^lY_{l,m}\}$$ - I don't think this is required. $$\{r^l Y_{l,m}\}$$ need only form a basis, not necessarily an orthogonal one.

• It's pretty obvious that $r^l$ form an basis (but not an orthogonal basis) for $L^2[0,\infty)$ - simply consider Taylor expansions. The Weierstrass approximation theorem only works for compact intervals. If the function vanishes in infinity, it can't be approximated by polynomials in either $L^2$ or $\sup$ norm. Take for instance $f(x) = \min(x, 1/x)$ Nov 12 '20 at 15:33
• @marmistrz Ok, fair point. I wasn't actually aware of this since I haven't seriously studied functional analysis. Hopefully the rest of the answer was still helpful - any questions on that? Nov 12 '20 at 15:50
• Yes, your answer is very helpful; it appears that most texts have a ton of hidden assumptions. Thanks a lot for the detailed answer. Nov 12 '20 at 19:03
• Finally, about the completeness of - this is what completeness exactly means, that the linear span of the set is dense. Of course every complete set can be orthogonalized. Nov 12 '20 at 19:07