kuruman said:
So you're saying that post #30 is nonsense because we are trying to solve Poisson's equation here? (I dread the answer).
No, I attempted "humor". I may have erroneously assumed that many people would be at least vaguely aware that those are among the contents of my textbook. Hence, saying that there is no good textbook on the subject would be along the lines of self-derogative humor.
To be on topic though ...
The Spherical harmonics ##Y_{\ell m}(\theta,\varphi)## are eigenfunctions of the angular part of the negative of the Laplace operator ##-\nabla^2##, i.e., in spherical coordinates
$$
-\nabla^2 = -\partial_r^2 - \frac 2r \partial_r + \frac 1{r^2} \hat \Lambda,
$$
where ##\hat \Lambda Y_{\ell m}(\theta,\varphi) = \ell(\ell+1) Y_{\ell m}(\theta,\varphi)##. Being constructed from individual Sturm-Liouville problems in the ##\theta## and ##\varphi## directions, it follows that
any function on ##\mathbb R^3## may be written on the form
$$
f(\vec r) = \sum_{\ell = 0}^\infty \sum_{m = -\ell}^\ell f_{\ell m}(r) Y_{\ell m}(\theta, \varphi)
$$
and that the ##Y_{\ell m}## are linearly independent so the expansion is unique. (From now on, let me just write ##\sum_{\ell, m}## for the sums to save typing.)
So, armed with this knowledge, we can diagonalise any linear problem involving the Laplace operator, in particular on the form of Poisson's equation
$$
-\nabla^2 u(\vec r) = \rho(\vec r)
$$
where I have used convenient constant normalisation. As both ##u## and ##\rho## are functions of position, they can both be written in terms of spherical coordinates and therefore expanded in terms of the spherical harmonics as described above with expansion coefficients ##u_{\ell m}(r)## and ##\rho_{\ell m}(r)##, respectively. (Note that the spherical harmonics form the expansion basis for functions on a sphere. We are therefore expanding each function on the separate spheres of radius ##r##, which is why the expansion coefficients depend on ##r##.)
Inserting the expansion into Poisson's equation directly leads to
$$
\sum_{\ell, m} \left[ - u_{\ell m}''(r) - \frac 2r u_{\ell m}'(r) + \frac 1{r^2} \ell(\ell + 1) u_{\ell m}(r)\right] Y_{\ell m}(\theta, \varphi) = \sum_{\ell, m} \rho_{\ell m}(r) Y_{\ell m}(\theta, \varphi).
$$
Now, since the spherical harmonics are linearly independent, it directly follows that
$$
- u_{\ell m}''(r) - \frac 2r u_{\ell m}'(r) + \frac 1{r^2} \ell(\ell + 1) u_{\ell m}(r) - \rho_{\ell m}(r) = 0,
$$
which is an ordinary inhomogeneous differential equation of second order that can be solved. Boundary conditions are typically found from ##u|_{r = 0}## being regular (although in actuality this will follow from the differential equation itself in coordinates other than spherical) and ##\lim_{r \to \infty} u = 0##. Solving for each ##u_{\ell m}(r)## solves Poisson's equation for an arbitrary source term ##\rho(\vec r)##.
Note that the functions ##\rho_{\ell m}## may be found through the inner product (restricting to real functions for simplicity)
$$
\langle f, g \rangle = \int_{\theta = 0}^\pi \int_{\varphi = 0}^{2\pi} f(\theta,\varphi) g(\theta,\varphi) \sin(\theta) d\theta \, d\varphi
$$
as
$$
\rho_{\ell m}(r) = \frac{\langle Y_{\ell m},\rho\rangle }{\langle Y_{\ell m}, Y_{\ell m}\rangle},
$$
where I am leaving the normalisation term ##\langle Y_{\ell m}, Y_{\ell m}\rangle## as written since there are a relatively large number of different normalisation conventions for spherical harmonics and there will always be someone complaining if you pick a particular one ...
Regardless of normalisation, the numerator is on the form
$$
\langle Y_{\ell m},\rho\rangle = \int_{\theta = 0}^\pi \int_{\varphi = 0}^{2\pi} Y_{\ell m}(\theta,\varphi) \rho(r,\theta,\varphi) \sin(\theta) d\theta \, d\varphi,
$$
making it clear that ##\rho_{\ell m}## is indeed a function of ##r##.
Now, looking at a source distribution such that ##\rho(r) = 0## for all ##r > R## for some ##R##, the differential equations become homogeneous
$$
- u_{\ell m}''(r) - \frac 2r u_{\ell m}'(r) + \frac 1{r^2} \ell(\ell + 1) u_{\ell m}(r)= 0
$$
for ##r > R##. This has the general solution
$$
u_{\ell m}(r) = \frac{A_{\ell m}}{r^{\ell+1}} + B_{\ell m} r^{\ell}.
$$
With the boundary condition at infinity, this leads to ##B_{\ell m}= 0## and therefore ##u_{\ell m}(r) = A_{\ell m}/r^{\ell + 1}## outside of the source distribution. This leads to the multipole expansion
$$
u = \sum_{\ell, m} \frac{A_{\ell m}}{r^{\ell + 1}} Y_{\ell m}(\theta,\varphi)
$$
for ##r > R##.
(I had planned to go a bit further but it is getting late ...)
Edit: Minor typos.
