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 
I would like to know how to evaluate the following triple integral with the help of spherical coordinates
$$\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \sqrt{{x^2+y^2+z^2}} \,dx \,dy\, dz$$

The relations between Cartesian coordinates and spherical ones are given by
${\displaystyle {\begin{aligned}x&=r\,\sin \theta \,\cos \varphi \\y&=r\,\sin \theta \,\sin \varphi \\z&=r\,\cos \theta \end{aligned}}}$
I know that a function is generally integrated over $\mathbb{R}^3$ by the following triple integral
$$ \ \int \limits _{\varphi =0}^{2\pi }\ \int \limits _{\theta =0}^{\pi }\ \int \limits _{r=0}^{\infty }f(r,\theta ,\varphi )r^{2}\sin \theta \,\mathrm {d} r\,\mathrm {d} \theta \,\mathrm {d} \varphi .$$
I found a numerical solution with Wolfram Alpha (0.960592), I tried to change the bounds of integration from Cartesian to spherical but I got a different numerical value.
Could someone please give a detailed solution showing how to change the limits of integration?
Thanks.
 A: You can split up the cube into $6$ tetrahedral chunks, so the integral is the total of the integrals over each tetrahedron, which are easier to set up. One such chunk involves the tetrahedron with vertices $(0,0,0),(1,0,0),(1,1,0),(1,1,1)$:
$$\int_0^{\pi/4}\int_{\tan^{-1}(\csc\theta)}^{\pi/2}\int_0^{\sec\theta\csc\varphi}\rho^3\sin\varphi\,\mathrm d\rho\,\mathrm d\varphi\,\mathrm d\theta=\frac{6\sqrt3-\pi+6\ln(7+4\sqrt3)}{144}\approx0.1600976746$$
The remaining integrals' bounds in spherical coordinates are
$$\int_0^{\pi/4}\int_{\tan^{-1}(\sec\theta)}^{\tan^{-1}(\csc\theta)}\int_0^{\sec\theta\csc\varphi}$$
(tetrahedron with vertices $(0,0,0),(1,0,0),(1,0,1),(1,1,1)$)
$$\int_0^{\pi/4}\int_0^{\tan^{-1}(\sec\theta)}\int_0^{\sec\varphi}$$
(tetrahedron with vertices $(0,0,0),(0,0,1),(1,0,1),(1,1,1)$)
$$\int_{\pi/4}^{\pi/2}\int_{\tan^{-1}(\sec\theta)}^{\pi/2}\int_0^{\csc\theta\csc\varphi}$$
(tetrahedron with vertices $(0,0,0),(0,1,0),(1,1,0),(1,1,1)$)
$$\int_{\pi/4}^{\pi/2}\int_{\tan^{-1}(\csc\theta)}^{\tan^{-1}(\sec\theta)}\int_0^{\csc\theta\csc\varphi}$$
(tetrahedron with vertices $(0,0,0),(0,1,0),(0,1,1),(1,1,1)$)
$$\int_{\pi/4}^{\pi/2}\int_0^{\tan^{-1}(\csc\theta)}\int_0^{\sec\varphi}$$
(tetrahedron with vertices $(0,0,0),(0,0,1),(1,1,0),(1,1,1)$)
By symmetry, the value of the original integral is

$$\frac{6\sqrt3-\pi+6\ln(7+4\sqrt3)}{24}\approx0.9605919565$$

A: $\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}
A: 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.
