(This question may belong on Physics. I put in on Mathematics because I saw more similar questions here.)
I have a seperable density function that is expressed in spherical polar coordinates as: $$n_{lm}(\mathbf{r})=Y_{lm}(\Omega)\rho_{lm}(r)$$ where $Y_{lm}(\Omega)$ is a spherical harmonic, $\Omega$ is the angular components of $\mathbf{r}$, $(\theta,\phi)$ and $n_{lm}(r)$ is the radial component of the density function.
I can use any one of the many ways to solve the Poisson equation
$$ \nabla^2 V_{lm}(\mathbf{r})=-\frac{1}{\epsilon_0}\rho_{lm}(\mathbf{r}) \tag{1} $$
to obtain the potential, $U(\mathbf{r})$ from the density, and indeed I do this in my code, successfully. My question is not about how to solve the Poisson equation, but about the form of the radial Poisson equation. The radial Poisson equation is given in several papers [e.g [Goddard, 1995] (https://doi.org/10.1103/PhysRevB.52.2348)] as:
$$\left[\frac{d^2}{dr^2}-\frac{l(l+1)}{r^2}\right]U_{lm}(r) =-4\pi r\rho_{lm}(r)$$
I want to know what I am missing to derive this form of the radial Poisson equation. Here is my attempt.
I will start with some statements. $\epsilon_0=0$ and I will write the potentials in the form $$ V_{lm}(\mathbf{r})=Y_{lm}(\Omega)U_{lm}(r)\tag{2} $$.
Substituting the definition of $\nabla^2$ in spherical polar coordinates, $$\nabla^2 = \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 \varphi^2}$$
into the Poisson equation \eqref{1}, together with \eqref{2} gives $$ Y_{lm}(\Omega)\frac{1}{r^2} \frac{\partial}{\partial r} \left(r^2 \frac{\partial}{\partial r} \right)U_{lm}(r) + U_{lm}(r)\left[\frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta} \left(\sin \theta \frac{\partial}{\partial \theta} Y_{lm}(\Omega)\right) + \frac{1}{r^2 \sin^2 \theta} \frac{\partial^2}{\partial \varphi^2}Y_{lm}(\Omega)\right] = -Y_{lm}(\Omega)\rho_{lm}(r)$$
Then comparing the angular terms with the spherical harmonic differential equation: $$\left[\frac{1}{\sin \theta} \frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial}{\partial \theta}\right)+\frac{1}{\sin ^2 \theta} \frac{\partial^2}{\partial \phi^2}+l(l+1)\right] Y_{lm}(\Omega)=0$$
I get
$$ \left[\frac{1}{r^2} \frac{d}{d r} \left(r^2 \frac{d}{d r} \right) - \frac{l(l+1)}{r^2}\right] U_{lm}(r) = -\rho_{lm}(r)$$
Expanding the derivative gives
$$ \left[\frac{d^2}{d r^2} +\frac{2}{r}\frac{d}{dr} - \frac{l(l+1)}{r^2}\right] U_{lm}(r) = -\rho_{lm}(r)$$
So, can someone show me how the extra $\frac{2}{r}\frac{d}{dr}$ disappears, and where the extra $r$ on the right hand side comes from?
Thanks.
P.s As a point of interest, while looking for similar questions, I found this question, where the left-hand side is the same as mine (although I would argue that theirs is a Laplace equation, rather than a Poisson equation).