How to evaluate this monsterous integral?

Let $r = \sqrt{x^2+y^2+z^2}$ and $r_1 = \sqrt{(x-1)^2+(y-1)^2 + (z-1)^2}$.

For any $R, \epsilon > 0$, let

• $B_0(R) = \{ (x,y,z) : r < R \}$ be the ball of radius $R$ centered at origin.
• $B_1(\epsilon) = \{ (x,y,z) : r_1 < \epsilon \}$ be the ball of radius $\epsilon$ centered at $(1,1,1)$.

Notice the integrand blows up as $\frac{1}{r_1^2}$ as $(x,y,z)$ approaches $(1,1,1)$ which get suppressed by a factor $r_1^2$ from the volume element. The triple integral should be interpreted as a limit of following form:

$$\iiint ( \cdots ) = \lim_{R\to\infty}\lim_{\epsilon\to 0}\int_{B_0(R)\setminus B_1(\epsilon)} (\cdots ) $$

By symmetry, the integral can be rewritten as

$$\begin{align} & -\frac23\iiint \frac{x(x-1)+y(y-1)+z(z-1)}{r_1^3} e^{-r^2} dV\\ = & -\frac23 \iiint \left(\vec{x}\cdot\frac{\vec{\nabla}r_1}{r_1^2}\right) e^{-r^2} dV\\ = & -\frac13 \iiint \vec{\nabla}\left(\frac{1}{r_1}\right)\cdot\vec{\nabla}e^{-r^2} dV\\ = & -\frac13 \iiint \left[ \vec{\nabla}\cdot\left(e^{-r^2}\vec{\nabla}\left(\frac{1}{r_1}\right)\right) - e^{-r^2}\nabla^2\left(\frac{1}{r_1}\right) \right] dV\\ \end{align} $$ Notice $\nabla^2\left(\frac{1}{r_1}\right) = 0$ outside $B_1(\epsilon)$, we can use divergence theorem to rewrite the integral as

$$-\frac13\lim_{R \to \infty}\lim_{\epsilon \to 0} \int_{\partial (B_0(R)\setminus B_1(\epsilon))} e^{-r^2}\vec{\nabla}\left(\frac{1}{r_1}\right)\cdot d\vec{S} = -\frac13\left[ \lim_{R \to \infty} \int_{\partial B_0(R)} - \lim_{\epsilon \to 0} \int_{\partial B_1(\epsilon)} \right] \left( e^{-r^2}\vec{\nabla}\left(\frac{1}{r_1}\right)\cdot d\vec{S} \right) $$ Because of the $e^{-r^2}$ factor, the integral over $\partial B_0(R)$ vanishes as $R \to \infty$.

For the remaining integral, since $e^{-r^2}$ is continuous near $(1,1,1)$, we can pull it outside the integral sign. As a result, the integral we seek is equal to

$$\frac13 \left(\lim_{(x,y,z)\to (1,1,1)} e^{-r^2}\right) \left(\lim_{\epsilon \to 0}\int_{\partial B_1(\epsilon)}\vec{\nabla}\left(\frac{1}{r_1}\right)\cdot d\vec{S}\right) = \frac13 e^{-3} (-4\pi) = -\frac{4\pi}{3e^3} $$

Another answer:

I noticed it also follows from this identity see identity

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {\partial f(x',y',z')}{\partial x'}\frac {x'-x}{4\pi r'^3}dx'dy'dz'+\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {\partial f(x',y',z')}{\partial y'}\frac {y'-y}{4\pi r'^3}dx'dy'dz'+ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {\partial f(x',y',z')}{\partial z'}\frac {z'-z}{4\pi r'^3}dx'dy'dz'=-f(x,y,z)$$

with $$f(x,y,z)=e^{-(x^2+y^2+z^2)}$$

it is easily seen when $x=y=z=1$ $\quad$ $f(1,1,1)=e^{-3}$

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {\partial f(x',y',z')}{\partial x'}\frac {x'-1}{4\pi r'^3}dx'dy'dz'=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(x'-x'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(y'-y'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(z'-z'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'$$

Therefore: $$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(x'-x'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\frac {-f(1,1,1)}{3}$$

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(x'-x'^2)e^{-(x'^2+y'^2+z'^2)}}{((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\frac {-4\pi e^{-3}}{3}$$