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.
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 \, . $$