I have a problem where I am supposed to calculate the volume charge density of a neutral hydrogen atom. The potential is given to be
$$
\Phi = k \frac{e^{-ar}}{r} \left(1 + \frac{ar}{2}\right)
$$
Now I tried to use the Poisson-Equation stating
$$
\Delta \Phi = \frac{\rho}{\varepsilon_0}
$$
which leads me to
$$
\rho = \Delta \left( \underbrace{\frac{q}{4 \pi \epsilon_0}}_{= k} \frac{e^{-\alpha r}}{r} \left( 1 + \frac{\alpha r}{2} \right) \right)
= k \Delta \Big( \frac{e^{-\alpha r}}{r} + \frac{\alpha e^{-\alpha r}}{2} \Big)
= k \left( \Delta \left( \frac{e^{-\alpha r}}{r} \right) + \frac{\alpha}{2} \Delta \left( e^{-\alpha r} \right) \right)
$$
Now I define $f = e^{-\alpha r}$ and $g = \frac{1}{r}$. The Laplacian of the product $fg$ is then
$$
\Delta(fg) = g \Delta(f) + f \Delta(g) + \nabla (f) \cdot \nabla (g)
$$
and the derivatives are
$$
\nabla(f) = – \alpha e^{-\alpha r} \hat{ \mathbf r} \qquad \qquad \Delta(f) = \alpha^2 e^{-\alpha r}
$$
$$
\nabla(g) = – \frac{1}{r^2} \hat{ \mathbf r} \qquad \qquad \Delta(g) = – 4 \pi \delta(r)
$$
$$
\implies \nabla(f) \cdot \nabla(g) = \frac{\alpha e^{-\alpha r}}{r^2}
$$
Inserting this back into the original equation yields
$$
\rho = k e^{-\alpha r}\Big( \frac{\alpha^2}{r} – 4 \pi \delta(r) + \frac{\alpha}{r^2} + \frac{\alpha^3}{2} \Big)
$$
However, this seems somewhat wrong to me since I would have expected the expression to be increasing from the origin and then decreasing after some $r=R$ since the potential of the electron hull should take over.
Can anybody either confirm that this is correct or show me where I made the mistake?
Apart from taking the derivatives like in cartesian coordinates, I have tried calculating the Laplacian by calculation in spherical coordinates as well using the spherical Laplacian
$$
\Delta \Phi = \frac{1}{r^2} \frac{\partial}{\partial r} \left(r^2 \frac{\partial \Phi}{\partial r}\right)
$$
but still got the same result.
Best Answer
It's not clear to me exactly why you're unhappy with the answer you get. I would suggest phrasing your expectations in terms of the total charge contained in a sphere of radius $r$, $$4\pi\int_{0^-}^r \rho(r')r'^2\,\text dr'.$$ This should give the positive charge of the nucleus at $r\rightarrow 0^+$ (because the proton is point-sized in this model!) and then drop monotonically to zero as $r$ increases through $1/\alpha$ and beyond to $r\rightarrow\infty$, and your sphere includes more and more of the electron cloud that surrounds (and neutralizes) the nucleus. If this fails then you definitely need to check your calculations.