$\frac{1}{|x-y|}=\frac{1}{\pi^3} \int_{R^3}\frac{1}{|x-z|^2}\frac{1}{|y-z|^2}dz$ Prove the following identity for  $x,y \in R^3:$ $$\dfrac{1}{|x-y|}=\dfrac{1}{\pi^3} \int_{R^3}\dfrac{1}{|x-z|^2}\dfrac{1}{|y-z|^2}dz$$
I tried multiple ways for example use green function,
put $x=(x_1,x_2,x_3),y=(y_1,y_2,y_3),z=(z_1,z_2,z_3)$ and calculate above formula but I failed
How can I solved this problem? thanks very much
 A: @Svyatoslav did admirably, but my intention with the hint was a simpler calculation. In spherical coordinates the integral after the hint and rotating to the appropriate coordinates system is
$$\frac{2\pi}{\pi^3}\int_0^\pi\int_0^\infty\frac{\sin\theta}{s^2-2as\cos\theta+a^2}\:ds\:d\theta$$
But now integrating $s$ first gives
$$\frac{2}{\pi^2}\int_0^\pi\int_0^\infty \frac{\sin\theta}{(s-a\cos\theta)^2+a^2\sin^2\theta}\:ds\:d\theta=\frac{2}{a\pi^2}\int_0^\pi\left[\tan^{-1}\left(\frac{s-a\cos\theta}{a\sin\theta}\right)\right]_0^\infty\:d\theta$$
$$= \frac{2}{a\pi^2}\int_0^\pi\frac{\pi}{2}+\tan^{-1}(\cot\theta)\:d\theta = \frac{2}{a\pi^2}\int_0^\pi\frac{\pi}{2}\:d\theta + 0 = \frac{1}{a}$$
as desired. The second term vanishes by symmetry because $\cot\theta$ is an odd function about $\theta=\frac{\pi}{2}$.
A: One of the options, probably, not the most rational one, is to evaluate the integral directly. If we choose this option, we have to define the system of coordinate - to choose the center point and axis direction. We can choose $y$ as a center point, and direct the axis $Z$ along the vector $\vec a=\vec x-\vec y$. Denoting $|\vec a|=a=|\vec x-\vec y|$, our integral takes the form
$$I=\dfrac{1}{\pi^3} \int_{R^3}\dfrac{1}{|x-z|^2}\dfrac{1}{|y-z|^2}dz=\dfrac{1}{\pi^3} \int_{R^3}\dfrac{1}{|s-(x-y)|^2}\dfrac{1}{|s|^2}ds$$
In the polar system of coordinates, and in accordance with our choice of the polar axis $Z$,
$$=\dfrac{1}{\pi^3}\int_0^{2\pi}d\phi\int_0^\infty \frac{s^2ds}{s^2}\int_0^\pi\frac{\sin\theta \,d\theta}{s^2-2sa\cos \theta+a^2}=\frac{1}{a\pi^2}\int_0^\infty \frac{ds}{s}\int_{-1}^1\frac{dx}{\frac{s^2+a^2}{2sa}-x}$$
where we made the substitution $\cos \theta =x$. Integrating with respect to $x$
$$=\frac{1}{a\pi^2}\int_0^\infty \frac{ds}{s}\ln\frac{s^2+2sa+a^2}{s^2-2sa+a^2}=\frac{2}{a\pi^2}\int_0^\infty \frac{ds}{s}\ln\frac{1+s}{|1-s|}$$
Splitting the interval by $[0;1]$ and $[1;\infty)$, and making the change $t=\frac{1}{s}$ in the second integral,
$$I=\frac{4}{a\pi^2}\int_0^1 \frac{ds}{s}\ln\frac{1+s}{1-s}=-\frac{4}{a\pi^2}\int_0^1 \frac{ds}{s}\ln\frac{1-s}{1+s}$$
Making the substitution $t=\frac{1-s}{1+s}$
$$I=-\frac{8}{a\pi^2}\int_0^1\frac{\ln t}{1-t^2}dt=-\frac{8}{a\pi^2}\int_0^1\big(1-t^2-t^4-...\big)\ln t\,dt$$
Integrating by part every term
$$=-\frac{8}{a\pi^2}\Big(-1-\frac{1}{3^2}-\frac{1}{5^2}-...\Big)=\frac{8}{a\pi^2}\Big(1+\frac{1}{2^2}+\frac{1}{3^2}+..\Big)-\frac{8}{a\pi^2}\Big(\frac{1}{2^2}+\frac{1}{4^2}+..\Big)$$
$$=\frac{8}{a\pi^2}\Big(\zeta(2)-\frac{1}{4}\zeta(2)\Big)=\frac{8}{a\pi^2}\frac{3}{4}\frac{\pi^2}{6}=\frac{1}{a}$$
