Computing the energy of the electrostatic field for a Gaussian distribution

Question

I want to compute, as a function of $\sigma$, the electrostatic energy of a Gaussian distribution of charges (omitting the constants): $${1\over 2}\int \rho\ \Phi\ d{\bf y} = {1\over 2}\int_{\mathbb R^3} e^{-{\bf y^2}\over \sigma^2} d{\bf y}\int_{\mathbb R^3} e^{-{\bf x}^2\over \sigma^2} {1\over ||\bf x - y||} d{\bf x} ={1\over 2}\int_{\mathbb R^3\times \mathbb R^3} e^{-{{\bf x^2 + y^2}\over \sigma^2}} {1\over ||\bf x - y||} d{\bf x}d{\bf y} . $$ Also, according to the energy formulae, this amounts to compute $$ \int_{\mathbb R^3} E^2 d{\bf y}= \int_{\mathbb R^3} d{\bf y} \left (\int_{\mathbb R^3} e^{-{\bf x}^2\over \sigma^2} {1\over ({\bf x - y})^2} d{\bf x}\right)^2 $$ (not obvious if this helps). Any idea? I would be satisfied also, by a numerical simulation that would give this function of $\sigma$ for sigma in a range $(10^{-12}, 10^{-3})$, but I don't know how to do such a simulation.

Answer 1

Denoting $x=||\bf x||, \,y=||\bf y||$ and choosing the polar system of coordinates ($\bf y$ is directed along the axis $Z$) $$I(\sigma)=\int_{\mathbb R^3\times \mathbb R^3} e^{-{{\bf x^2 + y^2}\over \sigma^2}} {1\over ||\bf x - y||} d{\bf x}d{\bf y}$$ $$=8\pi^2\int_0^\infty y^2dy\int_0^\infty x^2e^{-{{x^2 + y^2}\over \sigma^2}}dx\int_0^\pi\frac{\sin\theta \,d \theta}{\sqrt{x^2+y^2-2xy\cos\theta}}$$ Taking $\sin\theta d\theta=-d(\cos\theta)$ and integrating, $$I(\sigma)=8\pi^2\int_0^\infty\int_0^\infty e^{-{{x^2 + y^2}\over \sigma^2}}x\,y\,\big(|x+y|-|x-y|\big)dxdy$$ Using the polar system of coordinates ($x=r\cos\phi;\,y=r\sin\phi$) $$=8\pi^2\int_0^\infty e^{-\frac{r^2}{\sigma^2}}r^4dr\int_0^\frac{\pi}2\sin\phi\cos\phi\big(|\sin\phi+\cos\phi|-|\sin\phi-\cos\phi|\big)d\phi$$ $$=8\pi^2\frac{\sqrt 2}3\int_0^\infty e^{-\frac{r^2}{\sigma^2}}r^4dr=\sqrt 2\,\pi^\frac52\sigma^5$$