Evaluate $\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{(z-w)^2+(x-y)^2} \, dw \, dz \, dy \, dx$ This page contains an interesting identity$$\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{(z-w)^2+(x-y)^2} \, dw \, dz \, dy \, dx=\frac{1}{15} \left(\sqrt{2}+2+5 \log \left(\sqrt{2}+1\right)\right)$$
Which calculates the average distance between two random points in the unit square. I'm wondering if there's an elementary solution? Any help will be appreciated.
 A: By symmetry:
$$I=4\int _0^1\int _0^1\int _0^x\int _0^z\sqrt{(z-w)^2+(x-y)^2} \, dw \,  dy \, dx dz$$
Substituting $w \to z w, \quad y \to x y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2(1-w)^2+x^2(1-y)^2}~ x z \, dw \,  dy \, dx dz$$
Substituting $w \to 1-w, \quad y \to 1-y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2 w^2+x^2 y^2}~ x z \, dw \,  dy \, dx dz$$
Substituting $x^2 = u, \quad z^2 = v$:
$$I=\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{v w^2+u y^2}~ \, dw \,  dy \, du dv$$
Substituting $v w^2=p, \quad u y^2=q$:
$$I=\int _0^1\int _0^1\int _0^{y^2}\int _0^{w^2}  \sqrt{p+q}~ \, dp dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{2}{3} \int _0^1\int _0^1\int _0^{y^2}\left((q+w^2)^{3/2}-q^{3/2} \right) dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{4}{15} \int _0^1\int _0^1\left((y^2+w^2)^{5/2}-y^5-w^5 \right)  \frac{dy dw}{y^2 w^2}$$
By symmetry:
$$I=\frac{8}{15} \int _0^1\int _0^w\left((y^2+w^2)^{5/2}-y^5-w^5 \right)  \frac{dy dw}{y^2 w^2}$$
Substituting: $y = w s$:
$$I=\frac{8}{15} \int _0^1\int _0^1 w^2 \left((1+s^2)^{5/2}-s^5-1 \right)  \frac{ds dw}{s^2 }$$
$$I=\frac{8}{45} \int _0^1 \left((1+s^2)^{5/2}-s^5-1 \right)  \frac{ds }{s^2 }$$
The last single integral can be evaluated by different methods, and the result agrees with the one from the OP.
One good method here is integration by parts. I leave this as an exercise.
