<s>[INST] 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 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.
 [/INST] 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}}$$
</s>