Evaluating $\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \sqrt{{x^2+y^2+z^2}}\, dx \,dy \,dz$ by converting to spherical coordinates

Split up the cube into $6$ tetrahedral chunks. One such tetrahedron has vertices $\{(0,0,0),(0,0,1),(1,0,1),(1,1,1)\}$ and is enclosed by the planes $y=0$, $z=1$, $y=x$, and $z=x$:

$$\begin{align*} I_1 &= \int_0^\tfrac\pi4\int_0^{\arctan(\sec\theta)}\int_0^{\sec\varphi} \rho^3\sin\varphi \,\mathrm d\rho\,\mathrm d\varphi\,\mathrm d\theta \\ &= \frac14 \int_0^\tfrac\pi4 \int_0^{\arctan(\sec\theta)} \tan\varphi \sec^3\varphi \, d\varphi \, d\theta \\ &= \frac14 \int_0^\tfrac\pi4 \int_1^{\sqrt{1+\sec^2\theta}} f^2 \, df \, d\theta & f=\sec\varphi \\ &= \frac1{12} \int_0^\tfrac\pi4 \left(2+\tan^2\theta\right)^{3/2} \,d\theta - \frac\pi{48} \\ &= \frac13 \int_0^{\operatorname{arcsch}\sqrt2} \frac{\cosh^4u}{\cosh(2u)} \, du - \frac\pi{48} & \tan\theta = \sqrt2\,\sinh u \\ &= \frac{4\operatorname{arcsch}\sqrt2 - \pi}{48} + \frac1{12} \int_0^{\operatorname{arcsch}\sqrt2} \left(2\cosh^2u + \frac{1}{2\cosh^2u-1}\right) \, du & \text{partial fractions} \\ &= \frac{2\sqrt3 + 8\operatorname{arcsch}\sqrt2 - \pi}{48} + \frac1{12} \int_0^{\operatorname{arcsch}\sqrt2} \frac{\operatorname{sech}^2u}{1+\tanh^2u} \, du \\ &= \frac{6\sqrt3 + 24\operatorname{arcsch}\sqrt2 - \pi}{144} \end{align*}$$

By symmetry, the remaining integrals $I_2$ through $I_5$ (listed below) share the same value as $I_1$, and the closed form for OP's integral follows.

$$\begin{array}{l|l|l} i & \text{vertices} & \text{integration range} \\ \hline 2 & \{(0,0,0),(1,0,0),(1,0,1),(1,1,1)\} & \displaystyle \int_0^\tfrac\pi4 \int_{\arctan(\sec\theta)}^{\arctan(\csc\theta)} \int_0^{\sec\theta\csc\varphi} \\ 3 & \{(0,0,0),(1,0,0),(1,1,0),(1,1,1)\} & \displaystyle \int_0^\tfrac\pi4 \int_{\arctan(\csc\theta)}^\tfrac\pi2 \int_0^{\sec\theta\csc\varphi}\\ 4 & \{(0,0,0),(0,1,0),(1,1,0),(1,1,1)\} & \displaystyle\int_\tfrac\pi4^\tfrac\pi2\int_{\arctan(\sec\theta)}^\tfrac\pi2\int_0^{\csc\theta\csc\varphi} \\ 5 & \{(0,0,0),(0,1,0),(0,1,1),(1,1,1)\} & \displaystyle\int_\tfrac\pi4^\tfrac\pi2\int_{\arctan(\csc\theta)}^{\arctan(\sec\theta)}\int_0^{\csc\theta\csc\varphi} \\ 6 & \{(0,0,0),(0,0,1),(1,1,0),(1,1,1)\} & \displaystyle\int_\tfrac\pi4^\tfrac\pi2\int_0^{\arctan(\csc\theta)}\int_0^{\sec\varphi} \end{array}$$

Replacing $\theta$ by $\dfrac\pi2-\theta$ reveals $I_{1,2,3}$ are respectively duplicated by $I_{6,5,4}$. Fubini's theorem $(\star)$ can in turn be used to express $I_{1,2,3}$ in terms of a common integral.

$$\begin{align*} 4I_1 &= \int_0^\tfrac\pi4 \int_0^{\arctan(\sec\theta)} \sec^4\varphi \sin\varphi \, d\varphi \, d\theta \\ &= \int_0^\tfrac\pi4 \int_0^1 f\sec^2\theta\sqrt{1+f^2\sec^2\theta} \, df \, d\theta & f=\cos\theta\tan\varphi \\ &= \underbrace{\int_0^1 \int_0^1 f\sqrt{1+f^2\left(1+t^2\right)} \, df \, dt}_{\mathcal I} & t=\tan\theta \\[2ex] 4I_3 &= \int_0^\tfrac\pi4 \int_{\arctan(\csc\theta)}^\tfrac\pi2 \sec^4\theta \csc^3\varphi \, d\varphi \, d\theta \\ &= \int_0^\tfrac\pi4 \int_0^1 \sec^4\theta\sin\varphi\sqrt{1+f^2\sin^2\theta} \, df \, d\theta & f=\sin\theta\cot\varphi \\ &= \int_0^1 \int_0^1 t \sqrt{1+t^2\left(1+f^2\right)} \, df \, dt & t=\tan\theta \\ &= \int_0^1 \int_0^t \sqrt{1+f^2+t^2} \, df \, dt & f\to\frac ft \\ &= \int_0^1 \int_f^1 \sqrt{1+f^2+t^2} \, dt \, df & (\star) \\ &= \int_0^1 \int_1^\tfrac1f f\sqrt{1+f^2\left(1+t^2\right)} \, dt \, df & t\to ft\\ &= \int_0^1 \int_f^1 \frac f{t^3} \sqrt{t^2\left(1+f^2\right)+f^2} \, dt \, df & t\to\frac 1t \\ &= \int_0^1 \int_0^t \frac f{t^3} \sqrt{t^2\left(1+f^2\right)+f^2} \, df \, dt & (\star) \\ &= \mathcal I & f\to ft \\[2ex] 4I_2 &= \int_0^\tfrac\pi4 \int_{\arctan(\sec\theta)}^{\arctan(\csc\theta)} \sec^4\theta\csc^3\varphi \, d\varphi \, d\theta \\ &= \int_0^\tfrac\pi4 \int_{\sec\theta}^{\csc\theta} \sec^4\theta \frac{\sqrt{1+\varphi^2}}{\varphi^3} \, df\, d\theta & f=\tan\varphi \\ &= \int_0^\tfrac\pi4 \int_1^{\cot\theta}\frac{\sec^2\theta}{f^3}\sqrt{1+f^2\sec^2\theta} \, d\theta \, df \\ &= \int_0^\tfrac\pi4 \int_{\tan\theta}^1 \sec^2\theta \sqrt{f^2+\sec^2\theta} \, d\theta \, df \\ &= \int_0^1 \int_t^1 \sqrt{1+f^2+t^2} \, df\,dt & t=\tan\theta \\ &= \int_0^1 \int_0^f \sqrt{1+f^2+t^2} \, dt\,df & (\star) \\ &= \mathcal I & t\to ft \end{align*}$$

$$\implies \sum_{i=1}^6 I_i = \boxed{\frac{6\sqrt3-\pi+12\log\left(2+\sqrt3\right)}{24}}$$

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} &\bbox[5px,#ffd]{\int_{0}^{1}\int_{0}^{1}\int_{0}^{1} \root{{x^{2} + y^{2} + z^{2}}}\dd x\,\dd y\,\dd z} = \iiint_{\pars{0,1}^{\,\,3}}r\,\,\dd^{3}\vec{r} \end{align} Note that

$\ds{\nabla\cdot\pars{r\vec{r}} = \pars{1 \times{\vec{r} \over r}}\cdot\vec{r} + r\nabla\cdot\vec{r} = 4r \implies r = {1 \over 4} \nabla\cdot\pars{r\vec{r}}}$ such that \begin{align} &\bbox[5px,#ffd]{\int_{0}^{1}\int_{0}^{1}\int_{0}^{1} \root{{x^{2} + y^{2} + z^{2}}}\dd x\,\dd y\,\dd z} \\[5mm] = &\ {1 \over 4} \iiint_{\pars{0,1}^{\,\,3}}\nabla\cdot\pars{r\vec{r}}\,\,\dd^{3}\vec{r} = {1 \over 4}\iint_{\partial\pars{0,1}^{\,\,3}}\,\, r\,\vec{r}\cdot\dd\vec{S} \end{align} where I used Gauss-Ostrogradsky Divergence Theorem.

• The contribution to the integral from the faces which intersect the origen of coordinates vanishes out because $\ds{\vec{r} \perp \dd\vec{S}}$.
• Each contribution from the remaining three faces are equal between them. So, I evaluate just one integral in one of the faces ( for instance, the one which is perpendicular to $\ds{\hat{z}}$ and multiply the integral result by $\ds{3}$ ). Namely, \begin{align} &\bbox[5px,#ffd]{\int_{0}^{1}\int_{0}^{1}\int_{0}^{1} \root{{x^{2} + y^{2} + z^{2}}}\dd x\,\dd y\,\dd z} \\[5mm] = &\ 3\pars{{1 \over 4}\int_{0}^{1}\int_{0}^{1}\root{x^{2} + y^{2} + 1^{2}}\dd x\,\dd y} \\[5mm] = &\ {3 \over 4}\int_{0}^{\pi/2}\int_{0}^{\infty} \root{\rho^{2} + 1}\ \times \\[2mm] &\ \bracks{\rho\cos\pars{\phi} < 1} \bracks{\rho\sin\pars{\phi} < 1}\rho\,\dd\rho\,\dd\phi \\[5mm] = &\ {3 \over 8}\int_{0}^{\pi/2}\int_{0}^{\infty} \root{\rho + 1}\ \times \\[2mm] &\ \phantom{{3 \over 8}} \bracks{\rho < \min\braces{{1 \over \cos^{2}\pars{\phi}}, {1 \over \sin^{2}\pars{\phi}}}} \dd\rho\,\dd\phi \\[5mm] = &\ {3 \over 8}\int_{0}^{\pi/4} \int_{0}^{1/\cos^{2}\pars{\phi}} \root{\rho + 1}\dd\rho\,\dd\phi \\[2mm] + &\ {3 \over 8}\int_{\pi/4}^{\pi/2} \int_{0}^{1/\sin^{2}\pars{\phi}} \root{\rho + 1}\dd\rho\,\dd\phi \\[5mm] = &\ {1 \over 4}\int_{0}^{\pi/4} \braces{\bracks{{1 \over \cos^{2}\pars{\phi}} + 1}^{3/2} - 1}\dd\phi \\[2mm] + &\ {1 \over 4}\int_{\pi/4}^{\pi/2} \braces{\bracks{{1 \over \sin^{2}\pars{\phi}} + 1}^{3/2} - 1}\dd\phi \\[5mm] = &\ {1 \over 2}\ \underbrace{\int_{0}^{\pi/4} \bracks{\sec^{2}\pars{\phi} + 1}^{3/2}\,\dd\phi} _{\ds{{\root{3} \over 2} + {\pi \over 6} + 2\on{arccoth}\pars{\root{3}}}}\ -\ {\pi \over 8} \\[5mm] = &\ \bbx{{\root{3} \over 4} - {\pi \over 24} + \on{arccoth}\pars{\root{3}}} \approx 0.9606 \\ & \end{align}

Substitute:

$$z=\sqrt{x^2+y^2} t$$

$$I=\int_{0}^{1} \int_{0}^{1} \int_{0}^{1/\sqrt{x^2+y^2}} (x^2+y^2) \sqrt{{1+t^2}} dt dy dx$$

Integrating w.r.t. $t$:

$$I=\frac{1}{2} \int_{0}^{1} \int_{0}^{1} (x^2+y^2) \left(\frac{1}{\sqrt{x^2+y^2}}\sqrt{1+\frac{1}{x^2+y^2}}+\sinh^{-1} \frac{1}{\sqrt{x^2+y^2}} \right) dy dx$$

Using the symmetry:

$$I= \int_{0}^{1} \int_{0}^{x} (x^2+y^2) \left(\frac{1}{\sqrt{x^2+y^2}}\sqrt{1+\frac{1}{x^2+y^2}}+\sinh^{-1} \frac{1}{\sqrt{x^2+y^2}} \right) dy dx$$

Now we can use polar coordinates:

$$x=r \cos \phi$$

$$y=r \sin \phi$$

$$0<\sin \phi<\cos \phi, \qquad 0<\phi< \frac{\pi}{4}$$

$$0<r< \frac{1}{\cos \phi}$$

$$I= \int_{0}^{\frac{\pi}{4}} \int_{0}^{\frac{1}{\cos \phi}} r^3 \left(\frac{1}{r}\sqrt{1+\frac{1}{r^2}}+\sinh^{-1} \frac{1}{r} \right) dr d \phi$$

$$I= \int_{0}^{\frac{\pi}{4}} \int_{0}^{\frac{1}{\cos \phi}} \left(r\sqrt{r^2+1}+r^3 \sinh^{-1} \frac{1}{r} \right) dr d \phi$$

We have:

$$\int_{0}^{\frac{1}{\cos \phi}} r\sqrt{r^2+1} dr= \frac{1}{3} \left(\frac{1}{\cos^2 \phi}+1 \right)^{3/2}-\frac{1}{3}$$

$$\int_{0}^{\frac{1}{\cos \phi}} r^3 \sinh^{-1} \frac{1}{r} dr= \frac{1}{4 \cos^4 \phi} \sinh^{-1} \cos \phi+ \frac{1-2 \cos^2 \phi}{12 \cos^3 \phi} \sqrt{1+\cos^2 \phi}+\frac16$$

So we have a complicated expression:

$$I= \frac{\pi}{24}+ \frac13 \int_{0}^{\frac{\pi}{4}} \left(\left(\frac{1}{\cos^2 \phi}+1 \right)^{3/2}-1 \right) d \phi+ \\+ \frac{1}{4} \int_{0}^{\frac{\pi}{4}} \left( \frac{1}{\cos \phi}\sinh^{-1} \cos \phi+ \frac{1}{3} (1-2 \cos^2 \phi)\sqrt{1+\cos^2 \phi} \right) \frac{d \phi}{\cos^3 \phi}$$

These are elliptic kind of integrals, though some of them might be elementary. Substitution $\cos \phi=s$ seems prudent here.

I hope this might be helpful.

Maybe using spherical coordinates from the start is better, but I haven't figured out the correct bounds yet either.