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.
Best Answer
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.