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$

Question

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$?

Answer 1

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).

Answer 2

$$\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}$$

Answer 3

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}) $$