Expectation value $\bigg\langle \frac{1}{|\mathbf r_1-\mathbf r_2|} \bigg\rangle$ for the ground-state of He atom The approximate ground state for the two electrons of the He atom is given by $$\psi_0(r_1,r_2) = \frac{8}{\pi a^3}e^{-2(r_1+r_2)/a}$$, where $a$ is the Bohr radius.
I want to calculate the expectation value $\big\langle \frac{1}{|\mathbf r_1-\mathbf r_2|} \big\rangle$, but I really don't see how to proceed. The integrand seems to be correct.
We have 
\begin{align} \require{cancel}
\bigg\langle \frac{1}{|r_1-r_2|} \bigg\rangle &= \bigg( \frac{8}{\pi a^3} \bigg)^2
\cancel{\int_{r_1=0}^{\infty}\int_{r_2=0}^{\infty} \frac{1}{\sqrt{(r_1-r_2)^2}}
e^{-4(r_1+r_2)/a} dr_1 dr_2} \\
&= \ldots
\end{align}
Any hints would be appreciated.

The integrals $$\int xe^{-\alpha x}dx = e^{-\alpha x}\bigg( -\frac{x}{\alpha} - \frac{1}{\alpha} \bigg)$$
and 
$$
\int x^2 e^{-\alpha x} dx = e^{-\alpha x} \bigg( -\frac{x^2}{\alpha}-\frac{2x}{\alpha^2} - \frac{2}{\alpha^3} \bigg)
$$
are given as additional information.

As @pawel_winzig pointed out, we need to integrate over all space:
\begin{align}
\bigg\langle \frac{1}{|\mathbf r_1-\mathbf r_2|} \bigg\rangle &= \bigg( \frac{8}{\pi a^3} \bigg)^2
\int_{\mathbf{r}_1 \in \mathbb R^3} 
\int_{\mathbf{r}_2 \in \mathbb R^3}
 \frac{1}{\sqrt{(\mathbf r_1-\mathbf r_2)^2}}
e^{-4(r_1+ r_2)/a} d\mathbf r_1 d \mathbf r_2 
\end{align}
 A: Let us first take a look at the integrand 
$$
\frac{e^{-4(r_1+r_2)/a}}{|\mathbf r_1 - \mathbf r_2|}d\mathbf r_1 d\mathbf r_2
$$
and its Jacobian 
$$
d\mathbf r_1 d\mathbf r_2 = \big( r_1^2 \sin(\theta_1) \operatorname{dr_1 d\theta_1 d\phi_1}\big) \cdot \big(\operatorname{r_2^2 \sin(\theta_2) dr_2 d\theta_2 d\phi_2} \big).
$$ 
The reciprocal can be rewritten as 
$$
\frac{1}{|\mathbf r_1 - \mathbf r_2|}=\frac{1}{\sqrt{r_1^2+r_2^2-2r_1r_2\cos(\gamma)}}
$$
where $\gamma =\angle(\mathbf r_1, \mathbf r_2)$ is the angle between the vectors. 
Since the order of integration doesn't matter, we can assume that $\theta_1=0$, while we integrate over $\theta_2$; effectively placing the first particle on the $z$-axis, at first. Hence, we can set $\color{blue}{\gamma  := \theta_2 }$.
We can now write our integrand with a substitution 
$
\color{green}{u = r_1^2+r_2^2-2r_1 r_2 \cos(\theta_2) \ , \ du = 2r_1r_2\sin(\theta_2)}:
$
$$
\frac{e^{-4(r_1+r_2)/a}}{|\mathbf r_1 - \mathbf r_2|}d\mathbf r_1 d\mathbf r_2
= \frac{1}{2} \frac{e^{-4(r_1+r_2)/a}}{\sqrt{u}}r_1r_2 \sin(\theta_1) \operatorname{dr_1 dr_2 d \theta_1 du d\phi_1 d\phi_2} .
$$
After integrating over the azimuthal angles $\phi_1, \phi_2$, which yields $4\pi^2$; our substitution $u$, which yields $\frac{1}{2} \sqrt{u}$; our polar angle $\theta_1$, which yields $2$; and  we are left with:
$$
\bigg( \frac{8}{\pi^3}\bigg) ^2 \cdot \frac{1}{2} \cdot \frac{1}{2} \cdot 4\pi^2 \cdot 2 \int_{u=(r_1-r_2)^2}^{(r_1+r_2)^2}\sqrt{u} e^{-4(r_1+r_2)/a} \operatorname{dr_1 dr_2 du} = \frac{32}{\pi^4} \int_{r_1=0}^\infty \int_{r_2=0}^\infty r_2^2 r_1 e^{-4(r_1+r_2)/a}\operatorname{dr_1 dr_2}
$$
for which we have our standard integrals mentioned in the post.
