Computing $\int_{\Bbb{B}^3\times\Bbb{B}^3}\frac{\rm{d}^3x\,\rm{d}^3y}{\|x-y\|}=\frac{32}{15}\pi^2$, where $\Bbb{B}^3$ is the 3-dimensional unit ball

Question

Let's call $\mathbb{B}^3$ the three-dimensional unit ball, then how to compute the six fold $$ \int_{\mathbb{B}^3\times \mathbb{B}^3}\frac{\text{d}^3x\,\text{d}^3y}{\|x-y\|}=\frac{32}{15}\pi^2 $$ I saw this question on a forum, but no one seemed willing to answer it so I put it here for help. Any useful help would be greatly appreciated.

Answer 1

Let's consider a bit more general case of two balls with radii $a$ and $b$; for the convenience let's take $b>a$. Then the desired integral can be written in the form $$I(a,b)=4\pi^2\int_0^\pi\sin \theta_1d\theta_1\int_0^b y^2dy\int_0^\pi\sin \theta d\theta\int_0^a \frac{x^2dx}{\sqrt{x^2+y^2-2xy\cos\theta}}$$ where $\theta$ is the angle between $\vec x$ and $\vec y$, and $\theta_1$ is the angle between $\vec y$ and axis Z (it means that we perform integration with respect to $\vec x$ first, and direct the Z' axis in polar system of coordinates along $\vec y$, and then perform integration with respect to $\vec y$).

Splitting the interval of integration ($x<y; x>y$), denoting $z=\cos \theta$, and integrating with respect to $\theta_1$, $$=8\pi^2\int_0^by^2dy\int_{-1}^1dz\int_0^y\frac{x^2dx}{\sqrt{x^2+y^2-2xyz}}+8\pi^2\int_0^ay^2dy\int_{-1}^1dz\int_y^a\frac{x^2dx}{\sqrt{x^2+y^2-2xyz}}$$ $$=8\pi^2\int_0^by\,dy\int_0^yx^2dx\int_{-1}^1\frac{dz}{\sqrt{1-2\frac{x}{y}z+\frac{x^2}{y^2}}}+8\pi^2\int_0^ay^2dy\int_y^ax\,dx\int_{-1}^1\frac{dz}{\sqrt{1-2\frac{y}{x}z+\frac{y^2}{x^2}}}$$ Then we can perform integration with respect to $z$, or just notice that $\frac{1}{\sqrt{1-2rz+r^2}}$ is the generating function of Legendre' polynomials (and due to integration all polynomials die out, except for $P_0(z)=1;\,\int_{-1}^1P_0dz=2$). We get $$I(a,b)=16\pi^2\int_0^b\frac{y^4}{3}dy+16\pi^2\int_0^ay^2\frac{a^2-y^2}{2}dy=\frac{16\pi^2}{15}(a^5+b^5)$$

Answer 2

We can focus on one integral first - in this case consider only

$$\int_{B(0,1)}\frac{d^3x}{|x-y|} = \int_{B(y,1)}\frac{d^3x}{|x|}$$

by translating the origin to the other variable. We can always rotate to a coordinate system for $x$ where $y$ always lies on the $x_3$ axis, giving us the equation of the ball

$$x^2-2|y|x_3 = 1-y^2$$

Converting to spherical coordinates, we obtain

$$r^2-2|y|r\cos\theta = 1-y^2 \implies \begin{cases}\theta = \cos^{-1}\left(\frac{r^2+y^2-1}{2|y|r}\right)\\ r = |y|\cos\theta+\sqrt{1-y^2\sin^2\theta}\end{cases}$$

Of course this $\theta$ upper bound is only valid when $r>1-|y|$, so we need to split the integral into two

$$\int_{B(y,1)}\frac{r^2\sin\theta \: d(r,\theta,\varphi)}{r}$$

$$ = \int_0^{2\pi}\int_0^{1-|y|}\int_0^{\pi}r \sin\theta \:d\theta\:dr\:d\varphi + \int_0^{2\pi}\int_{1-|y|}^{1+|y|}\int_0^{\cos^{-1}\left(\frac{r^2+y^2-1}{2|y|r}\right)}r \sin\theta \:d\theta\:dr\:d\varphi$$

$$= 2\pi(1-|y|)^2 + 2\pi\int_{1-|y|}^{1+|y|}\frac{1-y^2}{2|y|}+r-\frac{r^2}{2|y|}\:dr$$

$$= 2\pi - 4\pi|y| + 2\pi y^2 +4\pi |y| - \frac{8\pi}{3} y^2 = 2\pi -\frac{2\pi}{3}y^2$$

We can plug that polynomial into the integral over a unit ball and get our final answer

$$2\pi\int_{B(0,1)}d^3y - \frac{2\pi}{3} \int_{\Bbb{S}^2}d\Omega \int_0^1 r^4\:dr $$

$$= 2\pi\cdot\left(\frac{4\pi}{3}\right)-\frac{2\pi}{3}\cdot(4\pi)\cdot\left(\frac{1}{5}\right) = \boxed{\frac{32\pi^2}{15}}$$