Distribution identity involving Dirac delta

Question

I was reading the following paper https://arxiv.org/pdf/gr-qc/0109005.pdf (page 6 equation 11), where it says:
"In the sense of distributions it holds \begin{equation} \Delta\left( \frac{1}{r}e^{-r/l}\right)=\frac{1}{rl^2}e^{-r/l}-4\pi \delta \;\;\;'' \end{equation} Where $\Delta$ is the laplacian and $\delta$ the three dimensional Dirac delta. I can see this identity is true via Fourier transform method, but I can't understand the following.
$\frac{1}{r}e^{-r/l}$ is a function only dependent of $r$, so we can apply the known identity for the laplacian $$\Delta(uv)=u\Delta v+2\nabla u\cdot\nabla v+v\Delta u $$ only taking into account the radial contribution, but applying this formula we don't get the supposed result above. What's going on here?

Answer 1

One can use the product rule $\nabla^2(uv)=u\nabla^2(v)+2\nabla(u)\cdot \nabla (v)+v\nabla^2(u)$. Proceeding, we have in distribution

$$\begin{align} \nabla^2 \left(\frac{e^{-r/\ell}}{r}\right)&=\color{blue}{e^{-r/\ell}\nabla^2 \left(\frac{1}{r}\right)}+2\nabla \left(\frac1r\right)\cdot \nabla \left(e^{-r/ \ell}\right)+\color{red}{\frac1r\nabla^2(e^{-r/\ell})}\\\\ &=\color{blue}{-4\pi \delta(r)}+2\frac{e^{-r/\ell}}{\ell r^2}+\color{red}{\left(\frac1{\ell^2 r}-\frac2{\ell r^2}\right)e^{-r/\ell}}\\\\ &=-4\pi \delta(r)+\frac{1}{\ell^2 r}e^{-r/\ell} \end{align}$$

as was to be shown!