Integrating this expression over two spheres I'd like to analytically perform the following integral:
$$ I=\int\limits_V d^3\textbf x\,d^3\textbf y \frac{1}{(\textbf x-\textbf y)^2}, $$
where the volume $V$ of the integration is such that: $|\textbf x|\in [0,A]$ and $|\textbf y|\in [0,B]$, i.e. spherical for $x$ and $y$ separately.
I tried to get a solution from Mathematica, but the output confuses me a bit:
$$ I = \begin{cases} \frac{\pi^2}{3}\left[ 10A^3B+18AB^3+3\text{i} B^4\pi+6(A^4+2A^2B^2-3B^4)\text{ arctanh}(\frac{B}{A}) \right] & A>B \\ 4A^4\pi^2 & A=B \\ 2\pi^2\left[ AB(A^2+B^2)-(A^2-B^2)^2\text{ arctanh}(\frac{A}{B}) \right] & B>A \end{cases} $$
I'd appreciate any help. Thanks in advance!

I think I've managed to get a better expression out of Mathematica:
$$ I = 2\pi² \left[ AB(A^2+B^2) -(A^2 -B^2)^2 \text{ arctanh}\left(\frac{<}{>}\right)  \right] $$
...where $>$ and $<$ are placeholders for $A$ and $B$, depending on which is larger or smaller.
Still, I'd love to know how to do this analytically!

Let's start with $I=\int\limits_{V_{x}}\int\limits_{V_y} d^3\textbf x\,d^3\textbf y \frac{1}{(\textbf x-\textbf y)^2}$$=\int\limits_{V_{y}}(\int\limits_{V_x} d^3\textbf x\, \frac{1}{(\textbf x-\textbf y)^2})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}{R^2+r^2-2rR\cos\theta_2}\right)$
$$I=\int\limits_{0}^A\int\limits_{0}^{2\pi}\int\limits_{0}^{\pi}r^2dr\sin\theta_1d\theta_1d\phi_1F(B,r)$$
For the convenience we define $B>A$ and orient the $z$ axis of the polar system of $X$ coordinates along the vector $\textbf y$. We can do this, because, as soon as we intergate over $R^3$ it is our choice of how to orient the coordinates system in the space. Integrating over $\theta_2$ and $\phi_2$ we get:
$F(B,r)=2\pi\int\limits_{0}^B\int\limits_{0}^{\pi}R^2dR\sin\theta_2d\theta_2\frac{1}{R^2+r^2-2rR\cos\theta_2}$$=2\pi\int\limits_{0}^B\frac{R^2}{2Rr}\log\frac{(R+r)^2}{(R-r)^2}dR$
Evaluating the last integral we have to be careful with the singular point $R=r$ (in fact we evaluate the integral principal value). To show our actions step by step let's regularize the integral by means of small $\epsilon$ near $R=r$, and in the end we will set $\epsilon$ to zero.
$F(B,r)=\pi\int\limits_{0}^B\frac{R}{r}\log\frac{(R+r)^2}{(R-r)^2}dR$$=\frac{2\pi}{r}\int\limits_{0}^BR\log(R+r)dR-\frac{2\pi}{r}\int\limits_{0}^{r-\epsilon}R\log(r-R)dR-\frac{2\pi}{r}\int\limits_{r+\epsilon}^BR\log(R-r)dR$.
Integrating and taking  limit at $\epsilon\to0$ we get $$F(b,r)=\frac{\pi}{r}(B^2-r^2)\log\frac{B+r}{B-r}+2\pi{B}$$ And, integrating over $\theta_1$ and $\phi_1$we finally get the desired integral in the form $$I=4{\pi}^2\int_0^A(B^2-r^2)\log(\frac{B+r}{B-r})rdr+8{\pi}B\int_0^Ar^2dr$$
Integrating by part (to get rid of logarithm) we get$$I=4{\pi}^2(\frac{B^2r^2}{2}-\frac{r^4}{4})\log(\frac{B+r}{B-r})|_0^A-4{\pi}^2\int_0^A(\frac{B^2r^2}{2}-\frac{r^4}{4})(\frac{1}{B+r}+\frac{1}{B-r})rdr+\frac{8{\pi}^2BA^3}{3}=$$ $$={\pi}^2(2B^2A^2-A^4)\log(\frac{B+A}{B-A})-4{\pi}^2B\int_0^A\frac{\frac{B^2}{2}-\frac{r^2}{2}+\frac{B^2}{2}}{B^2-r^2}r^2dr+\frac{8{\pi}^2BA^3}{3}$$ Finally, $$I={\pi}^2(2B^2A^2-A^4)\log(\frac{B+A}{B-A})-{\pi}^2\frac{2BA^3}{3}-{\pi}^2B^3\left(-2A+B\log(\frac{B+A}{B-A})\right)+\frac{8{\pi}^2BA^3}{3}$$ $$I=2{\pi}^2BA(B^2+A^2)-{\pi}^2(B^2-A^2)^2\log\frac{B+A}{B-A};   B>A$$
