Potential energy of a spherically symmetric charge density in a spherically symmetric electrostatic potential

Question

I'm interested in calculating the potential energy of a spherically symmetric charge density in a spherically symmetric electrostatic potential. More specific, I'm currently trying to calculate the potential energy of a Debye-distributed charge distribution (but displaced!)

$$ \rho(\vec{r})\propto \frac{1}{|\vec{r}-\vec{r}'|}\exp^{-\frac{|\vec{r}-\vec{r}'|}{\lambda_{d1}}} $$

in either a potential of a point charge or again a Debye potential

$$ \Phi(r)\propto 1/r \;\;\;\text{or}\;\;\;\Phi(r)\propto \frac{1}{r}\exp^{-r/\lambda_{d2}} \;. $$ Note that the charge distribution is allowed to be displaced. Physically I'm looking at two simplified atomic shells from different atoms (projectile-target) interacting.

In principal this can be done with the integral $$ W = \frac{1}{2}\int\rho\Phi dV \; , $$

I'm basically at a loss how to simplify the integral. More or less the best I got was something like (exact position of $\rho$ doesn't matter, only distance, so put on z-axis $\vec{r}'=z'\vec{e}_z$; do $\varphi$ integration)

$$ W = \frac{2\pi}{2}\int \rho(|r\vec{e}_r-z'\vec{e}_z|)\Phi(r)r^2\sin{\theta}\, drd\theta \; . $$

I can't remember of figure out how to simplify the $r\vec{e}_r-z'\vec{e}_z$-stuff more than just a few simple steps or useless substitution. I guess there are some typical "tricks" for those kind of (Green-function-like) integrals which I'm obviously missing.

Answer 1

I'm basically at a loss how to simplify the integral.

I am not sure this needs to be done. In spherically symmetric cases, the volume element, $dV$ is quite simple to deal with. In your case, you have: $$ \begin{align} W & = \frac{1}{2} \int dV \ \rho\left(r\right) \ \Phi\left(r\right) \\ & = \frac{4 \pi \ C}{2} \int dr \ r^{2} \left( r^{-2} e^{-r/\lambda_{D}} \right) \left( r^{-1} e^{-r/\lambda_{D}} \right) \\ & = \frac{4 \pi \ C}{2} \int dr \ \left( r^{-1} e^{-2 r/\lambda_{D}} \right) \end{align} $$ where $C$ is some constant. The direct solution to the integral is the exponential integral, $Ei(z)$, given by: $$ Ei\left(z\right) = - \int_{-z}^{\infty} dt \ \frac{e^{-t}}{t} $$ which can be solved numerically quite easily. If you want an analytical approach, you can use the Laplace transform after integration by parts. To do this, start by getting rid of the annoying extra terms by substituting $r \rightarrow \lambda_{D} x/2$, then integrate by parts to get: $$ C \int dx \ x^{-1} \ e^{-x} = C\left( e^{-x} \ln{x} + \int dx \ e^{-x} \ln{x} \right) $$

Now we can use the Laplace transform on the 2nd term on the right-hand side to get: $$ \mathcal{L} \int dx \ e^{-x} \ln{x} = -\frac{ \ln\lvert{1 + s}\rvert }{ s } + \frac{ \gamma + \ln\lvert{1 + s}\rvert }{ 1 + s } $$ where $s$ is the transform variable and $\gamma$ is the Euler-Mascheroni constant.

Note the numerical integral above is valid for arbitrary values of $z$ but the analytical result below is probably limited in its range of validity. I am sure some of my more mathematically inclined peers could explain in more detail.

Answer 2

So I actually integrated this a while ago, I thought I should share the result. It turned out rather simple with the right substitution.

By the way I'm talking about both potential and charge disrtibution as a Debye-like. Point charge is straight forward...

First we put the center of the charge density on the z-axis as mentioned in my question. I get

$$ \vec{r}'=z'\vec{e}_z $$ $$ |r\vec{e}_r-z'\vec{e}_z|^2 = (r\cos\varphi\sin\theta)^2+(r\sin\varphi\sin\theta)^2+(r\cos\theta-z')^2=r^2-2rz'\cos\theta+z'^2\,. $$

Then I substitute $$ u^2=r^2-2rz'\cos\theta+z'^2 $$ $$ \frac{d}{d\theta}u=\frac{rz'\sin\theta}{u} $$

$$ I=\int \rho(|\vec{r}-\vec{r}'|)\Phi(r)dV = 2\pi\int_0^\infty\left(\int_{u(0)}^{u(\pi)}\frac{ru}{z'}\rho(u)\Phi(r)\,du\right)dr $$ $$ u(0)^2=r^2-z'^2\;\;\;\;\;u(\pi)^2=r^2+z'^2 $$

One has to watch out for the correct signs because of the squares after that. Integration $du$ and $dr$ and some simplification and replacing back $z'=r'$ because of symmetry yields

$$ I=\frac{4\pi}{r'}\frac{(e^{-r'/\lambda_{d1}}-e^{-r'/\lambda_{d2}})\lambda_{d1}^2\lambda_{d2}^2}{(\lambda_{d1}^2-\lambda_{d2}^2)} $$