Evaluate $\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\sqrt{x^2+y^2+z^2}\arctan{\sqrt{x^2+y^2+z^2}}dxdydz$ I am trying to evaluate this integral:
$$\displaystyle\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\sqrt{x^2+y^2+z^2}\arctan{\sqrt{x^2+y^2+z^2}}dxdydz$$
I coverted the region by using polar coordinate system for xy-plane and it became:
$$\displaystyle\int_{0}^{1}\int_{0}^{\frac{\pi}{4}}\int_{0}^{\frac{1}{\cos\theta}}r\sqrt{r^2+z^2}\arctan{\sqrt{r^2+z^2}}drd\theta dz+\displaystyle\int_{0}^{1}\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{\frac{1}{\sin\theta}}r\sqrt{r^2+z^2}\arctan{\sqrt{r^2+z^2}}drd\theta dz$$
After IBP
$$\int r\sqrt{r^2+z^2}\arctan{\sqrt{r^2+z^2}}dr=\frac{1}{3}(r^2+z^2)^{\frac{3}{2}}\arctan{\sqrt{r^2+z^2}}-\frac{1}{6}(r^2+z^2)+\frac{1}{6}\ln(r^2+z^2+1)+C $$ and replace upper and lower limits, it is scary-looking.
My experience in solving this tyle of integral is still weak, I think there is a way to get rid of one variable from the beginning integral, but i can't figure it out. Hope everyone can help me, thanks.
 A: A simpler way using symmetry.
Firstly denote $r=\sqrt{x^2+y^2+z^2}$ and split the integrand
$$
I=\int_0^1\int_0^1\int_0^1r\arctan r~dV=\frac\pi2\underbrace{\int_0^1\int_0^1\int_0^1r~dV}_{I_1}-\underbrace{\int_0^1\int_0^1\int_0^1r\operatorname{arccot} r ~dV}_{I_2}
$$
The former is the average distance from a point in a unit cube to the origin. Use its symmetry to split the cube into 6 parts, in each of which the integral is identical. Substitute $p=xz,q=yz$ in one of them and perform the easy integral over $z$
$$
I_1=6\int_0^1\int_0^z\int_0^y\sqrt{x^2+y^2+z^2}dxdydz\\
=6\int_0^1\int_0^1\int_0^1z^3\sqrt{p^2+q^2+1}dpdqdz=\frac32\int_0^1\int_0^1\sqrt{p^2+q^2+1}dpdq
$$
Now switch to polar coordinate $(p,q)\to (\rho,\theta)$
$$
\begin{align}
&I_1=\frac32\int_0^1\int_0^1\sqrt{p^2+q^2+1}dpdq
\\ =&\frac32\int_0^{\pi/4}\int_0^{\sec\theta}\sqrt{\rho^2+1}~\rho d\rho d\theta
\\ =&\frac12\int_0^{\pi/4}d\theta\cdot\int_0^{\sec\theta}\frac32\sqrt{\rho^2+1}~d(\rho^2+1)
\\ =&\frac12\int_0^{\pi/4}\frac{\sqrt{1+\cos^2\theta}^3}{\cos^3\theta}-1~d\theta
\\ =&\frac12\left(\arcsin\left(\frac{\sin \theta }{\sqrt{2}}\right)+\frac{\sin\theta~\sqrt{1+\cos^2\theta}}{2\cos^2\theta }\right.
\\&\left.\left.+\log \left(\frac{1+\sin\theta\sqrt{1+\cos^2\theta}}{\cos^2\theta}\right)\right)\right|_0^{\pi/4}-\frac\pi8
\\ =&\frac{\sqrt{3}}{4}-\frac{\pi}{24}+\frac12\log \left(2+\sqrt3\right)
\end{align}
$$
The last integral has an elementary antiderivative and can be evaluated quickly substituting $s=\sin\theta$.
By the identity $\displaystyle\operatorname{arccot}\lambda=\int_0^1\frac{\lambda ~dw}{\lambda^2+w^2}$ and the symmetry of the latter
$$
\begin{align}
&I_2=\int_0^1\int_0^1\int_0^1r\cdot\left(\int_0^1\frac{r~dw}{w^2+r^2}\right)dV
\\ =&\int_0^1\int_0^1\int_0^1\int_0^1\frac{x^2+y^2+z^2}{x^2+y^2+z^2+w^2}dxdydzdw
\\ =&3\int_0^1\int_0^1\int_0^1\int_0^1\frac{x^2}{x^2+y^2+z^2+w^2}dxdydzdw
\\ =&\frac34\int_0^1\int_0^1\int_0^1\int_0^1\frac{x^2+y^2+z^2+w^2}{x^2+y^2+z^2+w^2}dxdydzdw
\\ =&\frac34
\end{align}
$$
Combining the results above and the desired integral follows
$$
I=\frac{\pi \sqrt{3}}{8}-\frac{\pi ^2}{48}+\frac{\pi}{4}   \log \left(2+\sqrt{3}\right)-\frac34
$$
A: Very similar steps from (Almost) Impossible Integrals, Sums, and Series, pages 261-261, can be applied where there we have that
$$\int_0^1\left( \int_0^1 \sqrt{x^2+y^2}\arctan\left(\frac{1}{\sqrt{x^2+y^2}}\right)\textrm{d}x\right)\textrm{d}y$$ $$=\int_0^1 \left(\int_0^1 \left(\int_0^1 \frac{x^2+y^2}{x^2+y^2+z^2}\textrm{d}z\right)\textrm{d}x\right)\textrm{d}y.$$
Similarly, we have that
$$\int_0^1  \left( \int_0^1 \left(\int_0^1 \sqrt{x^2+y^2+z^2}\arctan\left(\frac{1}{\sqrt{x^2+y^2+z^2}}\right)\textrm{d}x\right)\textrm{d}y\right)\textrm{dz}$$
$$=\int_0^1 \left(\int_0^1 \left(\int_0^1\left(\int_0^1 \frac{x^2+y^2+z^2}{x^2+y^2+z^2+w^2}\textrm{d}w\right)\textrm{d}z\right)\textrm{d}x\right)\textrm{d}y,$$
where since the last quadruple integral is immediately done by symmetry, the rest of the work to do is straightforward in various ways.
