Integral of the form $\int_\Omega |\mathbf{r}_1 - \mathbf{r}_2| \exp\left(-\frac{1}{2}(r_1^2+r_2^2)\right)\mathrm{d}\mathbf{r}_2$ I am trying to derive the density of Hooke's atom.
I know this is a physics model, but I think my problem is more in the realm of mathematics than physics. Let me give a quick description of the model and problem.
For this model the wave function is given as:
$$ \Psi\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)=A\left(1+\frac{1}{2}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|\right)\exp\left(-\frac{1}{4}\left(r_{1}^{2}+r_{2}^{2}\right)\right) $$
Or, expanded out:
$$ \Psi\left(\mathbf{r}_{1},\mathbf{r}_{2}\right) = A\exp\left(-\frac{1}{4}\left(r_{1}^{2}+r_{2}^{2}\right)\right)+\frac{A}{2}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|\exp\left(-\frac{1}{4}\left(r_{1}^{2}+r_{2}^{2}\right)\right) $$
With $\mathbf{r}=\left<x,y,z\right>$, and $A$ being a constant that normalizes the wave function.
Also note that everything is real, i.e. no imaginary part of either the constant or the coordinates.
The density is now defined as:
$$ \rho\left(\mathbf{r}_{1}\right) = \int_\Omega \Psi\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)\Psi\left(\mathbf{r}_{1},\mathbf{r}_{2}\right) \mathrm{d}\mathbf{r}_{2} $$
Here $\Omega$ is over all space.
Inserting the wave function in the above it is easily found that:
$$  \rho\left(\mathbf{r}_{1}\right) =A^{2}\exp\left(-\frac{1}{2}r_{1}^{2}\right)\underset{I_{1}}{\underbrace{\int_{\Omega}\exp\left(-\frac{1}{2}r_{2}^{2}\right)\mathrm{d}\mathbf{r}_{2}}} \\
 +\frac{A^{2}}{4}\exp\left(-\frac{1}{2}r_{1}^{2}\right)\underset{I_{2}}{\underbrace{\int_{\Omega}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2}\exp\left(-\frac{1}{2}r_{2}^{2}\right)\mathrm{d}\mathbf{r}_{2}}} \\
 +A^{2}\exp\left(-\frac{1}{2}r_{1}^{2}\right)\underset{I_{3}}{\underbrace{\int_{\Omega}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|\exp\left(-\frac{1}{2}r_{2}^{2}\right)\mathrm{d}\mathbf{r}_{2}}} $$
Given that:
$$ \int_{\Omega}f\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f\left(x,y,z\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z $$
and,
$$ r=\sqrt{x^{2}+y^{2}+z^{2}} $$
The two first integrals $I_1$ and $I_2$ can easily be evaluated by separation of the cartesian coordinates to give:
$$ I_1 = \int_{\Omega}\exp\left(-\frac{1}{2}r_{2}^{2}\right)\mathrm{d}\mathbf{r}_{2} = \sqrt{2\pi}^{3} $$
and,
$$ I_2 =\int_{\Omega}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2}\exp\left(-\frac{1}{2}r_{2}^{2}\right)\mathrm{d}\mathbf{r}_{2}= r_{1}^{2}\sqrt{2\pi}^{3}+3\sqrt{2\pi}^{3} $$
Now my problem arises when trying to evaluate the third integral:
$$ I_3 = \int_{\Omega}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|\exp\left(-\frac{1}{2}r_{2}^{2}\right)\mathrm{d}\mathbf{r}_{2} $$
Writing it out in cartesian coordinates I get:
$$ I_3 = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\sqrt{\left(x_{1}-x_{2}\right)^{2}+\left(y_{1}-y_{2}\right)^{2}+\left(z_{1}-z_{2}\right)^{2}} \\
\exp\left(-\frac{1}{2}x_{2}^{2}\right)\exp\left(-\frac{1}{2}y_{2}^{2}\right)\exp\left(-\frac{1}{2}z_{2}^{2}\right)\mathrm{d}x_{2}\mathrm{d}y_{2}\mathrm{d}z_{2} $$
But from here on I am stuck.
I cannot see any way to separate the variable in $\sqrt{\left(x_{1}-x_{2}\right)^{2}+\left(y_{1}-y_{2}\right)^{2}+\left(z_{1}-z_{2}\right)^{2}}$.
From the solution on Wikipedia https://en.wikipedia.org/wiki/Hooke%27s_atom, the solution is going to be proportional to:
$$ I_3 \propto \left(r + \frac{1}{r}\right)\mathrm{erf}(r) $$
But I have been unable to find any relations with the error-function that could point me towards the solution.
What would be my next step to figure out how to solve $I_3$?
 A: Let $\mathbf{u}=\mathbf{r}_1-\mathbf{r}_2$, then
$$
I_3=\int d^3u \ |\mathbf{u}| \exp \left(- \frac{1}{2} (\mathbf{r}_1-\mathbf{u})^2\right)
$$
In spherical $u$ co-ordinates, we may align the z axis with the constant $\mathbf{r}_1$ to find
$$
I_3=2\pi \int\limits_0^\pi d\theta \ \sin\theta \int\limits_0^\infty du \ u^3 \exp \left( -\frac{1}{2}(u^2+r_1^2-2ur_1 \cos\theta)\right) 
$$
The $\theta$ integral may be performed using the substitution $z=\cos\theta$, resulting in
$$
I_3=\frac{2 \pi}{r_1} \int\limits_0^\infty du \ u^2 (e^{2ur_1}-1)e^{-(u+r_1)^2/2}
$$
And finally
$$
I_3=4 \pi e^{-r_1^2/2}+(2\pi)^{3/2}\left(\frac{1}{r_1}+r_1 \right)\operatorname{Erf}(r_1/\sqrt{2})
$$
A: Not a full solution, but too long for a comment.
$$I_3 = \int_{\Omega}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|\exp\left(-\frac{1}{2}r_{2}^{2}\right)\mathrm{d}\mathbf{r}_{2}$$
Choose a spherical system of coordinates with $Z$ axis directed along $\vec r_1$ (this is your right to identify the coordinate system; the result does not depend on your choice).
In this system the integral
$$I_3=\int_0^{2\pi}\int_0^\infty\int_0^\pi\exp\left(-\frac{1}{2}r_{2}^{2}\right)\sqrt{r^2_1+r_2^2-2r_1r_2\cos\theta}\,\,\sin\theta \,d\theta\,r_2^2\,dr_2\, d\phi$$
You can integrate first over $\phi$ and $\theta$:
$$I_3=2\pi\int_0^\infty\int_{-1}^1\exp\left(-\frac{1}{2}r_{2}^{2}\right)\sqrt{r^2_1+r_2^2-2r_1r_2x}\,\,dx\,r_2^2\,dr_2$$
$$=\frac{2\pi}{3\,r_1}\int_0^\infty\exp\left(-\frac{1}{2}r_{2}^{2}\right)\biggl(\bigl((r_1+r_2)^2\bigr)^{\frac{3}{3}}-\bigl((r_1-r_2)^2\bigr)^{\frac{3}{3}}\biggr)r_2 dr_2$$
When integrating $\bigl((r_1-r_2)^2\bigr)^{\frac{3}{3}}$ you should be careful with the sign and limits: the expression is positive at any $r_2$. You should split the integration: from $0$ to $r_1$ and from $r_1$ to $\infty$ (and here the function $\mathrm{erf}(r)$ appears).
A: $$\newcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\abs}[1]{\left\lvert #1 \right\rvert}
\newcommand{\d}{\mathrm{d}}
\newcommand{\e}{\mathrm{e}}$$
Let $b$ be a half-integer. If $\Re a < b$, then
$$\int_{\mathbb{R}^{2b}} \abs{\vec{r}-\vec{x}}^{-2a}\e^{-r^2}\d^{2b} r
=\frac{\pi^{b}\Gamma(b-a)}{\Gamma(b)}M(a,b,-x^2)$$
where $M$ is a Kummer function:
$$\begin{aligned}
  \int_{\mathbb{R}^{2b}} \abs{\vec{r}-\vec{x}}^{-2a}\e^{-r^2}\d^{2b} r
  &=\frac{1}{\Gamma(a)}\int_{\mathbb{R}^{2b}}\int_0^{\infty}
    t^{a-1}\e^{-\abs{\vec{r}-\vec{x}}^2t-r^2}\d t\,\d^{2b}r \\
  &=\frac{1}{\Gamma(a)}\int_{\mathbb{R}^{2b}}\int_0^{\infty}
    t^{a-1}\e^{-(1+t)\abs{\vec{r}-\tfrac{t}{1+t}\vec{x}}^2
    -\tfrac{t}{1+t}x^2}\d t\,\d^{2b}r \\
  &=\frac{\pi^b}{\Gamma(a)}\int_0^{\infty}
    t^{a-1}(1+t)^{-b}\e^{-\tfrac{t}{1+t}x^2}\d t \\
  &=\frac{\pi^b}{\Gamma(a)}\int_0^1 t^{a-1}(1-t)^{b-a-1}\e^{-tx^2}\d t \\
  &=\frac{\pi^b\Gamma(b-a)}{\Gamma(b)}M(a,b,-x^2)
\end{aligned}$$
where we use
$$\begin{aligned}
  \frac{1}{k^a}
    &=\frac{1}{\Gamma(a)}\int_0^{\infty}
      t^{a-1}\e^{-kt}\d t \\
  {\textstyle\sum_it_i\abs{\vec{r}-\vec{x}_i}^2}
    &= ({\textstyle\sum_it_i})
      \abs{\vec{r}-\tfrac{\sum_it_i\vec{x}_i}{\sum_it_i}}^2
      +{\textstyle\sum_it_ix_i^2}
      -\tfrac{\left(\sum_it_i\vec{x}_i\right)^2}{\sum_it_i} \\
  \int_{\mathbb{R}^{2b}}\e^{-k r^2}\d^{2b}r
    &= \frac{\pi^b}{k^b} \\
  \int_0^{\infty}f(t)\d t
    &= \int_0^1 f(\tfrac{t}{1-t})\tfrac{\d t}{(1-t)^2} \\
  \int_0^1t^{a-1}(1-t)^{b-a-1}\e^{zt}\d t
    &=\frac{\Gamma(a)\Gamma(b-a)}{\Gamma(b)}M(a,b,z)\text{.}
\end{aligned}$$

$$\DeclareMathOperator{\erf}{erf}$$
When $2a$ and $2b$ are both odd integers, the Kummer function $M(a,b,z)$ is
expressible in terms of the error function and exponential functions. For
example,
$$\begin{aligned}
  M(\tfrac{1}{2},\tfrac{3}{2},-x^2) &= \tfrac{\pi^{1/2}}{2x}\erf x \\
  M(-\tfrac{1}{2},\tfrac{3}{2},-x^2)
    &=\tfrac{\pi^{1/2}}{2x}(\tfrac{1}{2}+x^2)\erf x +\tfrac{1}{2}\e^{-x^2}
  \text{.}
\end{aligned}$$
When $a$ is a nonpositive integer and $2b$ is an odd integer, $M(a,b,z)$ is
expressible in terms of Hermite polynomials. For example,
$$\begin{aligned}
  M(0,\tfrac{3}{2},x^2) &= \tfrac{1}{2x}H_1(x) &&= 1 \\
  M(-1,\tfrac{3}{2},x^2) &= -\tfrac{1}{12x}H_3(x) &&= 1 -\tfrac{2}{3} x^2
  \text{.}\end{aligned}$$
Therefore
$$\begin{aligned}
  \int_{\mathbb{R}^{3}} \abs{\vec{r}-\vec{x}}^{-1}\e^{-r^2}\d^{3} r
    &=\pi^{3/2}\frac{\erf x}{x} \\
  \int_{\mathbb{R}^{3}} \e^{-r^2}\d^{3} r
    &=\pi^{3/2} \\
  \int_{\mathbb{R}^{3}} \abs{\vec{r}-\vec{x}}\e^{-r^2}\d^{3} r
    &=\pi^{3/2}(\tfrac{1}{2x}+x)\erf x +\pi \e^{-x^2} \\
  \int_{\mathbb{R}^{3}} \abs{\vec{r}-\vec{x}}^2\e^{-r^2}\d^{3} r
    &=\pi^{3/2}(\tfrac{3}{2}+x^2)\text{.}
\end{aligned}$$
