Integral in relative and center of mass variables I'm trying to compute the following integral where $\vec{x}$ and $\vec{y}$ are integrated inside a sphere of radius $l$:
$$I = \int\ d^3xd^3y\frac{1}{|\vec{x} - \vec{y}|^2}$$
To do it, I perform the change of variables:
$$\vec{r} = \vec{x} - \vec{y}, \quad \vec{R} = \frac{\vec{x} + \vec{y}}{2}$$
This change is such that 
$$d^3xd^3y = d^3rd^3R$$
So now we have:
$$I = \int\ dr\sin\theta d\theta d\varphi\int\ R^2dR\sin\theta' d\theta' d\varphi'$$
Where $\vec{r}$ has spherical coordinates $(r, \theta, \varphi)$ and $\vec{R}$ has ($R, \theta', \varphi'$). My problem is that I do not know how to write the limits in this last integral. I've tried: $r \in (0, 2l), \theta \in (0, \pi/2), \varphi \in (0, 2\pi)$ and $R \in (0, l), \theta' \in (0, \pi), \varphi' \in (0, 2\pi)$ but these ones do not reproduce the original volume, i.e., $\int d^3xd^3y = (\frac{4\pi l^3}{3})^2$, so I'm really stuck. 

EDIT: the solution label by the green check gives a solution for I without making use of this change of variable. If anyone is able to solve I with $\{\vec{r}, \vec{R}\}$, feel free giving a new answer. I'll appreciate it.
 A: If you integrate over a sphere the distance $r = |x-y|$ can go from $0$ up to $2l$ but you have to consider the starting point $R$, which if it is already on the edge of the sphere then r can only be 0 so the limit for r can be: $r ∈ (0, 2(l-R)).$
So the formula becomes:
$$
I = \int_0^l R^2 (\int_0^{2(l-R)} dr) dR \int_0^\pi \sin(\theta) d\theta \int_0^{2\pi} d\phi \int_0^\pi \sin(\theta ') d\theta ' \int_0^{2\pi} d\phi ' = \frac{8}{3} \pi^2 l^4
$$
But this is different from the previous answer that gives $4\pi^2 l^2$ ($l$ in OP's notation, $R$ in previous answer's notation)
A: Here's an approach based only on the use of spherical coordinates:
For $\vec{x} \in B_R(0)$ we define
$$ J(\vec{x}) = \int \limits_{B_R(0)} \frac{\mathrm{d}^3 \vec{y}}{|\vec{x}-\vec{y}|^2} \, . $$
Then 
$$ I = \int \limits_{B_R(0)} \mathrm{d}^3 \vec{x} \,J(\vec{x}) \, .$$
We rotate the $\vec{y}$-coordinate system in such a way that $\vec{x}$ points in positive $y_3$-direction and introduce spherical coordinates to obtain
\begin{align}
J(\vec{x}) &= \int \limits_{B_R(0)} \frac{\mathrm{d}^3 \vec{y}}{|\vec{x}|^2 + |\vec{y}|^2 - 2 |\vec{x}| y_3} \\
&= \int \limits_0^R \mathrm{d} r \, r^2 \int \limits_{-1}^1 \mathrm{d} \cos(\theta) \int \limits_0^{2\pi} \mathrm{d} \phi \, \frac{1}{|\vec{x}|^2 + r^2 - 2 |\vec{x}| r \cos(\theta)} \\
&= \frac{2 \pi}{|\vec{x}|} \int \limits_0^R \mathrm{d} r \, r \ln \left(\frac{r + |\vec{x}|}{|r-|\vec{x}||}\right)
\end{align}
for $\vec{x} \neq 0$ and $J(0) = 4 \pi R$ . The remaining integral for $\vec{x} \neq 0$ can be split into two parts ($r < |\vec{x}|$ and $r > |\vec{x}|$), which are readily computed via integration by parts. The final result is
\begin{align} J(\vec{x}) &= 2 \pi \left[R + \frac{(R+|\vec{x}|)(R-|\vec{x}|)}{2 |\vec{x}|} \ln \left(\frac{R + |\vec{x}|}{R-|\vec{x}|}\right)\right] \\
&= 2 \pi R \left[1+(1-|\vec{x}|^2 /R^2) \frac{\operatorname{artanh}(|\vec{x}|/R)}{|\vec{x}|/R}\right]
\end{align}
for $\vec{x} \in B_R (0) \setminus \{0\}$ . 
Now we can evaluate $I$ (again using spherical coordinates and integration by parts):
\begin{align}
I &= 8 \pi^2 R \int \limits_0^R \mathrm{d} r \, r^2 \left[1+(1-r^2 /R^2) \frac{\operatorname{artanh}(r/R)}{r/R}\right] \\
&= 8 \pi^2 R^4 \int \limits_0^1 \mathrm{d} \rho \, \left[\rho^2 + \rho (1-\rho^2) \operatorname{artanh}(\rho)\right] \\
&= 8 \pi^2 R^4 \left[\frac{1}{3} + \int \limits_0^1  \mathrm{d} \rho \, \frac{1}{4} (1-\rho^2)^2 \frac{1}{1-\rho^2} \right] \\
&= 4 \pi^2 R^4 \, .
\end{align}
Similarly (but with less effort), we can compute
$$ \int \limits_{B_R(0) \times B_R(0)} \frac{\mathrm{d}^3 \vec{x} \mathrm{d}^3 \vec{y}}{|\vec{x}-\vec{y}|} = \frac{32}{15} \pi^2 R^5 \, . $$
