I have the following definite integral over 6D space, in spherical coordinates:
$$\int d^3r d^3r' \frac{e^{-2a(r+r')}}{|\mathbf{r-r'}|}$$
I am unsure how to approch this but I have a couple of ideas.
The first is to use the Laplace expansion of $\frac{1}{|\mathbf{r-r'}|}$ in spherical harmonics.
The second is to look for an 6D coordinate system but I have not been able to find anything on commonly-used coordinate systems in $D>3$.
On
Here is a PDE-centric approach. Denote
$$f(\mathbf{r}) = \int_{\Bbb{R}^3}d^3\mathbf{r}'\frac{e^{-2ar'}}{|\mathbf{r-r'}|}$$
Notice that this is the solution to the inhomogeneous Laplace equation, or in other words
$$\Delta f(\mathbf{r}) = 4\pi e^{-2ar}$$
And choosing $f$ to be purely radial, we can use
$$\Delta f(\mathbf{r}) = \Delta f(r) = \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{df}{dr}\right)$$
to compute $f(r)$:
$$\nabla f = \frac{4\pi\hat{\mathbf{r}}}{r^2}\int dr \:r^2e^{-2ar} = -\frac{\pi e^{-2ar}}{a^3}\left(2a^2+\frac{2a}{r}+\frac{1}{r^2}\right)\hat{\mathbf{r}}$$
$$\implies f = -\frac{\pi}{a^3}\int dr \:\left(2a^2e^{-2ar}\right)+\left(\frac{2ae^{-2ar}}{r}+\frac{e^{-2ar}}{r^2}\right) = \frac{\pi e^{-2ar}}{a^3}\left(a+\frac{1}{r}\right)$$
choosing all $+C$s to be zero to allow the functions to vanish as $r\to\infty$. Thus the integral we want is
$$U = \frac{\pi}{a^3}\int_{\Bbb{R}^3}d^3\mathbf{r}\:e^{-4ar}\left(a+\frac{1}{r}\right) = \frac{4\pi^2}{a^3}\int_0^\infty dr\:ar^2e^{-4ar}+re^{-4ar}$$ $$ = \frac{\pi^2}{a^5}\left(\frac{2!}{4^2}+\frac{1!}{4}\right) = \boxed{\frac{3\pi^2}{8a^5}}$$
Let's denote $$I(a)=\int d^3r_1 d^3r_2 \frac{e^{-2a(r_1+r_2)}}{|\mathbf{r_1-r_2}|}$$ Switching to the polar system of coordinate and choosing the direction of $\mathbf{r_2}$ as the polar axis in the system of coordinates of $\mathbf{r_1}$ $$I(a)=\int d^3r_2\int_0^\infty e^{-2a(r_1+r_2)}r_1^2dr_1\int_0^\pi\frac{\sin\theta_1d\theta_1}{\sqrt{r_1^2+r_2^2-2r_1r_2\cos\theta_1}}\int_0^{2\pi}d\phi_1$$ $$=8\pi^2\int_0^\infty r_2^2dr_2\int_0^\infty r_1^2dr_1\int_{-1}^1\frac {e^{-2a(r_1+r_2)}}{\sqrt{2r_1r_2}}\frac{dx}{\sqrt{\frac{r_1^2+r_2^2}{2r_1r_2}-x}}$$ Integrating with respect to $x$ $$=8\pi^2\int_0^\infty dr_2\int_0^\infty dr_1\frac{r_1^2r_2^2}{2r_1r_2}e^{-2a(r_1+r_2)}\big(r_1+r_2-|r_1-r_2|\big)$$ Making the substitution $r_1+r_2=x\,,\,\,r_1-r_2=y$ $$I(a)=\frac{\pi^2}2\int_0^\infty dx\int_{-x}^xdy(x+y)(x-y)e^{-2ax}\big(x-|y|\big)$$ $$=\pi^2\int_0^\infty e^{-2ax}dx\int_0^xdy(x^3-xy^2-x^2y+y^3)=\pi^2\frac5{12}\int_0^\infty e^{-2ax}x^4dx$$ $$\boxed{\,\,I(a)=\frac{5\pi^2}{16\,a^5}\,\,}$$