A multiple integral problem in $6$D: $I=\iint_{B_r\times B_r}\frac{1}{|x-y|}dxdy$ Let $B_r$ to be a ball centered at origin with radius $r$ in dimension $3$. My problem is the following multiple integral:
$$I=\iint_{B_r\times B_r}\frac{1}{|x-y|}dxdy$$
The result says that $I$ is exactly equivalent to $r^5$. The lower bound for $I$ is trivial but I have trouble in the upper bound. Moreover, can this result be extended to dimension $n$?
 A: Solution for 3-dimension space
$I(A,B)=\int\limits_{V_{x}}\int\limits_{V_y} d^3\textbf x\,d^3\textbf y \frac{1}{|\textbf x-\textbf y|}=\int\limits_{V_{y}}(\int\limits_{V_x} d^3\textbf x\, \frac{1}{|\textbf x-\textbf y|})d^3\textbf y=$$=\int\limits_{0}^A\int\limits_{0}^{2\pi}\int\limits_{0}^{\pi}r^2dr\sin\theta_1d\theta_1d\phi_1\left(\int\limits_{0}^B\int\limits_{0}^{2\pi}\int\limits_{0}^{\pi}R^2dR\sin\theta_2d\theta_2d\phi_2\frac{1}{\sqrt{R^2+r^2-2rR\cos\theta_2}}\right)$
$$I(A,B)=\int\limits_{0}^A\int\limits_{0}^{2\pi}\int\limits_{0}^{\pi}r^2dr\sin\theta_1d\theta_1d\phi_1F(B,r)$$
$A$ and $B$ are the radii of balls; for the convenience $B\geqslant{A}$ and  $z$ axis of the polar system of coordinates ($X$ space) is oriented along the vector $\textbf y$.
Integrating over $\theta_2$ and $\phi_2$
$$F(B,r)=2\pi\int\limits_{0}^B\int\limits_{0}^{\pi}R^2dR\sin\theta_2d\theta_2\frac{1}{\sqrt{R^2+r^2-2rR\cos\theta_2}}=$$$$\frac{2\pi}{\sqrt{2r}}\int\limits_{0}^B\int\limits_{-1}^1\frac{R^2}{\sqrt{R}\sqrt{\frac{r^2+R^2}{2Rr}-x}}dxdR=$$
$$=\frac{2\pi}{r}\int\limits_{0}^BR\Bigl(R+r-|R-r|\Bigr)dR=\frac{2\pi}{r}\Bigl(rB^2-\frac{r^3}{3}\Bigr)$$
$$I(A,B)=\int\limits_{0}^A\int\limits_{0}^{2\pi}\int\limits_{0}^{\pi}\frac{2\pi}{r}\Bigl(rB^2-\frac{r^3}{3}\Bigr)r^2dr\sin\theta_1d\theta_1d\phi_1=\frac{8\pi^2}{3}A^3\Bigl(B^2-\frac{1}{5}A^2)\quad (B>A)$$
For $B=A$ $$I=\frac{32\pi^2}{15}A^5$$
A: By setting $(ra,rb )=(x ,y)$ we see that $|x-y|^{-1}dxdy = r^{-1}|a-b| (r^3da)(r^3db) = r^5 |a-b|^{-1}dadb$ and the integration domain becomes the product of unit balls. So if the integral is finite when $r=1$, then for general $r$ it scales like $r^5$. In higher dimensions the answer is of course $r^{2n-1}. $
As for this unit radius case we can change variables to $u=x-y$ and $v=x+y$. Being a linear orthogonal change of coordinates it introduces a constant scale factor which I ignore. Then the integral is like
$$\int_{B-B}\int_{B+B}\frac{ dudv }{|u|}$$
Which is finite because $1/|u|$ is locally integrable in dimension 3. Written in this form it could be computed explicitly but I’ll leave that to you.
