# Electrostatic Potential Energy integral in spherical coordinates

I'm having trouble with evaluating an integral that arises from attempting to find the total energy of an electrostatic system consisting of two point charges, which involves an integral over all space. The problem is physical, but my bone of contention is purely mathematical. There's a bit on it in the background section if context is needed. In spherical coordinates, I found the integral to be:

$$U = \frac{Q_1 Q_2}{8\pi\varepsilon_0}\int_0^\infty \hspace{-5pt}\int_0^{\pi} \frac{r - R\cos(\theta)}{(r^2-2Rr\cos(\theta)+R^2)^{\frac{3}{2}}}\sin(\theta) \space d\theta \space dr$$

I hit a brick wall upon trying to evaluate the integral - ordinarily I would use a substitution in the single integral case, and let $$u = r^2-2Rrcos(\theta)+R^2$$ but am unsure of how to do so for a double integral when the variables are all mixed up. Is this integral possible to evaluate using elementary techniques? If anybody could point me in the right direction or determine that the integral doesn't exist (which might mean I have made a grievous mistake in getting to the integral), that would be much appreciated.

I anticipate that the integral should reduce to the simple form given by the standard formula used to evaluate the energy of a system of point charges, which is $$U = \frac{1}{4\pi\varepsilon_0}\frac{Q_1Q_2}{R}$$ if that will help.

### Background

I'm trying to calculate the total energy of a simple two charge system through the integral for electrostatic energy of a system given in Griffiths' book:

$$U = \frac{\epsilon_0}{2}\int_V E^2 dV$$

Where the volume is integrated across all space so the boundary term not shown here decays to zero. I think that this should yield the same answer as the standard formula given for point charges:

$$U = \frac{1}{4\pi\varepsilon_0}\frac{Q_1Q_2}{R}$$

But I'm having trouble evaluating the integral itself. I placed $$Q_1$$ on the origin of the coordinate axes and $$Q_2$$ on the $$z$$-axis a distance $$R$$ away from the first charge, and expanded the $$E^2$$ term:

$$E = E_1 + E_2 \quad \text{ so } \quad E^2 = E_1^2 + 2E_1 \centerdot E_2 + E_2^2$$

I found that the integral of the self terms diverges when evaluated, and, after reading through Griffiths, decided to discard the self-energy terms and only retain the energy due to the exchange term.

Letting $$r = \sqrt{x^2+y^2+z^2}$$ and $$r'= \sqrt{x^2+y^2+(z-R)^2}$$, I found the integral of the interaction term to be:

$$E_1 = \frac{1}{4\pi\varepsilon_0}\frac{Q_1}{r^3}\vec{r} \quad \text{ and } \quad E_2 \frac{1}{4\pi\varepsilon_0}\frac{Q_2}{r'^3}\vec{r}'$$

$$U = \epsilon_0\int_V E_1\centerdot E_2 \space dV = \frac{Q_1 Q_2}{16\pi^2\varepsilon_0}\int_V \frac{x^2 + y^2 + z^2-zR}{(x^2 + y^2 + z^2)^{\frac{3}{2}} \space (x^2+y^2+(z-R)^2)^{\frac{3}{2}}}\space dV$$

Converting to spherical coordinates, with $$r=\sqrt{x^2+y^2+z^2}$$, $$\theta$$ the angle from the z-axis and $$\varphi$$ the azimutal angle, where I have evaluated the azimuthal integral, I get the integral in spherical coordinates above.

• Is the upper limit correct? I think it should be $\pi$ instead of $2\pi$. The integral is actually easy but it's too lengthy. Try $u=\cos\theta$ to evaluate the inner integral and separate the integrand into two terms. $$\int\left(\frac{Ru}{\sqrt{(r^2+R^2-2rRu)^3}} - \frac{r}{\sqrt{(r^2+R^2-2rRu)^3}}\right)\,du$$ The first integral can be evaluated by letting $u=\frac{r^2+R^2}{2rR}\sin^2 t$ and the second integral, let $v=r^2+R^2-2rRu$. – Anastasiya-Romanova 秀 Sep 27 '14 at 12:42
• @Anastasiya-Romanova You are correct. The integral would be zero if the upper limit was $2\pi$. – David H Sep 27 '14 at 12:45
• I got the desired result by using a substitution $u = \frac{r^2+R^2 - t}{2rR}$ for the first integral above inspired by the suggested substitution, I found that substituting $u = \frac{r^2+R^2}{2rR} sin^2(t)$ caused the terminals to be imaginary, as r and R are strictly positive in this problem. – Eweler Sep 27 '14 at 19:05

Electric fields due to $Q_1$ and $Q_2$:

$$\vec{E}_1=\frac{Q_1}{4\pi\epsilon_0}\frac{x\hat{x}+y\hat{y}+z\hat{z}}{\left(x^2+y^2+z^2\right)^{3/2}},$$

$$\vec{E}_2=\frac{Q_2}{4\pi\epsilon_0}\frac{x\hat{x}+y\hat{y}+(z-R)\hat{z}}{\left(x^2+y^2+(z-R)^2\right)^{3/2}}.$$

Dot product of electric fields in spherical coordinates:

\begin{align} \vec{E}_1\cdot\vec{E}_2 &=\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{x^2+y^2+z^2-Rz}{\left(x^2+y^2+z^2\right)^{3/2}\left(x^2+y^2+(z-R)^2\right)^{3/2}}\\ &=\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{r^2-Rr\cos{\theta}}{r^3\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\\ &=\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{r-R\cos{\theta}}{r^2\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}.\\ \end{align}

Setting up the integral:

\begin{align} U &=\epsilon_0\int_{V}\vec{E}_1\cdot\vec{E}_2\,\mathrm{d}V\\ &=\epsilon_0\int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\vec{E}_1\cdot\vec{E}_2\,r^2\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\varphi\\ &=\epsilon_0\int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{r-R\cos{\theta}}{r^2\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,r^2\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\varphi\\ &=\frac{Q_1Q_2}{16\pi^2\epsilon_0}\int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\frac{r-R\cos{\theta}}{\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\varphi\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0}\int_{0}^{\pi}\int_{0}^{\infty}\frac{r-R\cos{\theta}}{\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\\ \end{align}

Substituting $r=Rx$ and $t=\cos{\theta}$,

\begin{align} U &=\frac{Q_1Q_2}{8\pi\epsilon_0}\int_{0}^{\pi}\int_{0}^{\infty}\frac{r-R\cos{\theta}}{\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0}\int_{0}^{\pi}\int_{0}^{\infty}\frac{Rx-R\cos{\theta}}{\left(R^2x^2-2R^2x\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,R\,\mathrm{d}x\,\mathrm{d}\theta\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\pi}\int_{0}^{\infty}\frac{x-\cos{\theta}}{\left(x^2-2x\cos{\theta}+1\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}x\,\mathrm{d}\theta\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\int_{0}^{\pi}\frac{x-\cos{\theta}}{\left(x^2-2x\cos{\theta}+1\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}\theta\,\mathrm{d}x\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\int_{-1}^{1}\frac{x-t}{\left(x^2-2xt+1\right)^{3/2}}\,\mathrm{d}t\,\mathrm{d}x\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\left[\frac{x-1}{x^2\sqrt{x^2-2x+1}}+\frac{x+1}{x^2\sqrt{x^2+2x+1}}\right]\,\mathrm{d}x\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\left[\frac{x-1}{x^2|x-1|}+\frac{x+1}{x^2|x+1|}\right]\,\mathrm{d}x\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\left[\frac{\operatorname{sgn}{(x-1)}}{x^2}+\frac{1}{x^2}\right]\,\mathrm{d}x\\ &=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{1}\left[\frac{-1}{x^2}+\frac{1}{x^2}\right]\,\mathrm{d}x+\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{1}^{\infty}\left[\frac{1}{x^2}+\frac{1}{x^2}\right]\,\mathrm{d}x\\ &=\frac{Q_1Q_2}{4\pi\epsilon_0\,R}\int_{1}^{\infty}\frac{\mathrm{d}x}{x^2}\\ &=\frac{Q_1Q_2}{4\pi\epsilon_0\,R}, \end{align}

as expected.