# 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.

• Considering the use of a $\sqrt{x^2+y^2+z^2}$, I feel a shift to spherical (instead of polar/cylindrical) coordinates will be the way to go. Commented Feb 11, 2022 at 7:51
• @EeveeTrainer Thank, can you give me a hint for the upper and lower limits when transform to spherical coordinates? Commented Feb 11, 2022 at 8:01
• Increase the symmetry of the problem by going $-1$ to $1$ and divide by the symmetry factor at the end. For the simpler problem with the integrand being indicator function for $r \leq R$ Commented Feb 11, 2022 at 8:35
• @EeveeTrainer Spherical coordinates are obviously the natural choice for the integrand, but a very unnatural choice for the region of integration. It’s not clear if the overall problem is any less difficult than where we started. Commented Feb 11, 2022 at 8:46
• One can get a numerical approximation by dividing the unit cube into N parts along each axis (N^3 smaller cubes) and using Rieman sums. So for N=10, 20, 30, .., the numerical approximations are 0.75752560813883478192 0.75855222922061113336 0.75874228522195874044 0.75880880046659419103 0.75883958675808790787 0.7588563099728932806 Commented Feb 11, 2022 at 10:10

Let $$\mathcal{I}$$ denote the value of the definite integral in question:

$$\mathcal{I}:=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}z\,\sqrt{x^{2}+y^{2}+z^{2}}\arctan{\left(\sqrt{x^{2}+y^{2}+z^{2}}\right)}\approx0.758894.$$

The derivation below will show that

$$\mathcal{I}=-\frac34-\frac{\pi^{2}}{48}+\frac{\pi\sqrt{3}}{8}+\frac{\pi}{4}\ln{\left(2+\sqrt{3}\right)}.$$

To begin, we can simplify this integral by taking advantage of its permutation symmetry, i.e., the integration intervals are the same for each variable and the integrand $$f(x,y,z)$$ is unchanged under any rearranging of $$(x,y,z)$$.

For any function $$f(x,y,z)$$ integrable over the unit cube and symmetric in its three variables, its integral can be split up into a sum of $$3!=6$$ parts as

\begin{align} \int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}z\,f(x,y,z) &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{0}^{y}\mathrm{d}z\,f(x,y,z)\\ &~~~~~+\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{y}^{x}\mathrm{d}z\,f(x,y,z)\\ &~~~~~+\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{x}^{1}\mathrm{d}z\,f(x,y,z)\\ &~~~~~+\int_{0}^{1}\mathrm{d}x\int_{x}^{1}\mathrm{d}y\int_{0}^{x}\mathrm{d}z\,f(x,y,z)\\ &~~~~~+\int_{0}^{1}\mathrm{d}x\int_{x}^{1}\mathrm{d}y\int_{x}^{y}\mathrm{d}z\,f(x,y,z)\\ &~~~~~+\int_{0}^{1}\mathrm{d}x\int_{x}^{1}\mathrm{d}y\int_{y}^{1}\mathrm{d}z\,f(x,y,z).\\ \end{align}

After changing all the orders of integration and relabeling variables as needed, we find that these six integrals are in fact all equal to each other, i.e.,

$$\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}z\,f(x,y,z)=6\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{0}^{y}\mathrm{d}z\,f(x,y,z).$$

Making use of the symmetry of $$\mathcal{I}$$, we can rewrite the triple integral as

\begin{align} \mathcal{I} &=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}z\,\sqrt{x^{2}+y^{2}+z^{2}}\arctan{\left(\sqrt{x^{2}+y^{2}+z^{2}}\right)}\\ &=6\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{0}^{y}\mathrm{d}z\,\sqrt{x^{2}+y^{2}+z^{2}}\arctan{\left(\sqrt{x^{2}+y^{2}+z^{2}}\right)};~~~\small{symmetry}\\ &=6\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\,x^{3}\sqrt{1+t^{2}+u^{2}}\arctan{\left(x\sqrt{1+t^{2}+u^{2}}\right)};~~~\small{\left[(y,z)=(xt,xu)\right]}\\ &=6\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\int_{0}^{1}\mathrm{d}x\,x^{3}\sqrt{1+t^{2}+u^{2}}\arctan{\left(x\sqrt{1+t^{2}+u^{2}}\right)}\\ &=6\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\int_{0}^{\sqrt{1+t^{2}+u^{2}}}\mathrm{d}v\,\frac{v^{3}\arctan{\left(v\right)}}{\left(1+t^{2}+u^{2}\right)^{3/2}};~~~\small{\left[x=\frac{v}{\sqrt{1+t^{2}+u^{2}}}\right]}\\ &=\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\,\frac{6}{\left(1+t^{2}+u^{2}\right)^{3/2}}\int_{0}^{\sqrt{1+t^{2}+u^{2}}}\mathrm{d}v\,v^{3}\arctan{\left(v\right)}.\\ \end{align}

The integration over $$v$$ in the last line above is in fact elementary, as can be seen by verifying the following derivative:

$$\frac{d}{dv}\left[v\left(3-v^{2}\right)+3\left(v^{4}-1\right)\arctan{\left(v\right)}\right]=12v^{3}\arctan{\left(v\right)};~~~\small{v\in\mathbb{R}}.$$

Transforming the outer double integral of $$\mathcal{I}$$ to polar coordinates, we can reduce the integral $$\mathcal{I}$$ to single-variable integrals as follows:

\begin{align} \mathcal{I} &=\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\,\frac{6}{\left(1+t^{2}+u^{2}\right)^{3/2}}\int_{0}^{\sqrt{1+t^{2}+u^{2}}}\mathrm{d}v\,v^{3}\arctan{\left(v\right)}\\ &=\int_{0}^{\frac{\pi}{4}}\mathrm{d}\theta\int_{0}^{\sec{\left(\theta\right)}}\mathrm{d}\rho\,\frac{6\rho}{\left(1+\rho^{2}\right)^{3/2}}\int_{0}^{\sqrt{1+\rho^{2}}}\mathrm{d}v\,v^{3}\arctan{\left(v\right)};~~~\small{\left[\left(t,u\right)=\left(\rho\cos{\left(\theta\right)},\rho\sin{\left(\theta\right)}\right)\right]}\\ &=\int_{0}^{\frac{\pi}{4}}\mathrm{d}\theta\int_{1}^{\sqrt{1+\sec^{2}{\left(\theta\right)}}}\mathrm{d}x\,\frac{6}{x^{2}}\int_{0}^{x}\mathrm{d}v\,v^{3}\arctan{\left(v\right)};~~~\small{\left[\sqrt{1+\rho^{2}}=x\right]}\\ &=\int_{0}^{\frac{\pi}{4}}\mathrm{d}\theta\int_{1}^{\sqrt{1+\sec^{2}{\left(\theta\right)}}}\mathrm{d}x\,\frac{1}{2x^{2}}\int_{0}^{x}\mathrm{d}v\,\frac{d}{dv}\left[v\left(3-v^{2}\right)+3\left(v^{4}-1\right)\arctan{\left(v\right)}\right]\\ &=\int_{0}^{\frac{\pi}{4}}\mathrm{d}\theta\int_{1}^{\sqrt{1+\sec^{2}{\left(\theta\right)}}}\mathrm{d}x\,\frac{1}{2x^{2}}\left[x\left(3-x^{2}\right)+3\left(x^{4}-1\right)\arctan{\left(x\right)}\right]\\ &=\frac12\int_{0}^{\frac{\pi}{4}}\mathrm{d}\theta\int_{1}^{\sqrt{1+\sec^{2}{\left(\theta\right)}}}\mathrm{d}x\,\left[\frac{3}{x}-x+\frac{3\left(x^{4}-1\right)}{x^{2}}\arctan{\left(x\right)}\right]\\ &=\frac12\int_{0}^{\frac{\pi}{4}}\mathrm{d}\theta\int_{1}^{\sqrt{2+\tan^{2}{\left(\theta\right)}}}\mathrm{d}x\,\frac{d}{dx}\left[-x^{2}+2\ln{\left(\frac{1+x^{2}}{2}\right)}+\frac{3+x^{4}}{x}\arctan{\left(x\right)}\right]\\ &=\frac12\int_{0}^{1}\mathrm{d}t\,\frac{1}{1+t^{2}}\int_{1}^{\sqrt{2+t^{2}}}\mathrm{d}x\,\frac{d}{dx}\bigg{[}-x^{2}+2\ln{\left(\frac{1+x^{2}}{2}\right)}\\ &~~~~~+\frac{3+x^{4}}{x}\arctan{\left(x\right)}\bigg{]};~~~\small{\left[\theta=\arctan{\left(t\right)}\right]}\\ &=\frac12\int_{0}^{1}\mathrm{d}t\,\frac{1}{1+t^{2}}\bigg{[}-\left(2+t^{2}\right)+2\ln{\left(\frac{3+t^{2}}{2}\right)}\\ &~~~~~+\frac{3+\left(2+t^{2}\right)^{2}}{\sqrt{2+t^{2}}}\arctan{\left(\sqrt{2+t^{2}}\right)}-\left(-1+\pi\right)\bigg{]}\\ &=\frac12\int_{0}^{1}\mathrm{d}t\,\frac{1}{1+t^{2}}\bigg{[}-\left(1+t^{2}\right)+2\ln{\left(\frac32\right)}-\pi\\ &~~~~~+\frac{7+4t^{2}+t^{4}}{\sqrt{2+t^{2}}}\arctan{\left(\sqrt{2+t^{2}}\right)}+2\ln{\left(1+\frac13t^{2}\right)}\bigg{]}\\ &=-\frac12\int_{0}^{1}\mathrm{d}t+\ln{\left(\frac32\right)}\int_{0}^{1}\mathrm{d}t\,\frac{1}{1+t^{2}}-\frac{\pi}{2}\int_{0}^{1}\mathrm{d}t\,\frac{1}{1+t^{2}}\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\frac{7+4t^{2}+t^{4}}{2\left(1+t^{2}\right)\sqrt{2+t^{2}}}\arctan{\left(\sqrt{2+t^{2}}\right)}+\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(1+\frac13t^{2}\right)}}{1+t^{2}}\\ &=-\frac12+\frac{\pi}{4}\ln{\left(\frac32\right)}-\frac{\pi^{2}}{8}+2\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{2+t^{2}}\right)}}{\left(1+t^{2}\right)\sqrt{2+t^{2}}}\\ &~~~~~+\frac12\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{2+t^{2}}\right)}}{\sqrt{2+t^{2}}}\\ &~~~~~+\frac12\int_{0}^{1}\mathrm{d}t\,\sqrt{2+t^{2}}\arctan{\left(\sqrt{2+t^{2}}\right)}\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(1+\frac13t^{2}\right)}}{1+t^{2}}.\\ \end{align}

Define the function $$\mathcal{J}:\mathbb{R}_{>0}^{2}\rightarrow\mathbb{R}$$ via the definite integral

$$\mathcal{J}{\left(a,b\right)}:=\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{a+t^{2}}\right)}}{\left(b+t^{2}\right)\sqrt{a+t^{2}}},$$

define the function $$\mathcal{K}:\mathbb{Z}_{\ge0}\times\mathbb{R}_{\ge0}\rightarrow\mathbb{R}$$ via the definite integral

$$\mathcal{K}{\left(n,a\right)}:=\int_{0}^{1}\mathrm{d}t\,\left(a+t^{2}\right)^{n-1/2}\arctan{\left(\sqrt{a+t^{2}}\right)},$$

and define the function $$F:\mathbb{R}\rightarrow\mathbb{R}$$ via the definite integral

$$F{\left(b\right)}:=\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(1+b^{2}t^{2}\right)}}{1+t^{2}}.$$

Then,

\begin{align} \mathcal{I} &=-\frac12+\frac{\pi}{4}\ln{\left(\frac32\right)}-\frac{\pi^{2}}{8}+2\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{2+t^{2}}\right)}}{\left(1+t^{2}\right)\sqrt{2+t^{2}}}\\ &~~~~~+\frac12\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{2+t^{2}}\right)}}{\sqrt{2+t^{2}}}\\ &~~~~~+\frac12\int_{0}^{1}\mathrm{d}t\,\sqrt{2+t^{2}}\arctan{\left(\sqrt{2+t^{2}}\right)}\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(1+\frac13t^{2}\right)}}{1+t^{2}}\\ &=-\frac12+\frac{\pi}{4}\ln{\left(\frac32\right)}-\frac{\pi^{2}}{8}+2\mathcal{J}{\left(2,1\right)}\\ &~~~~~+\frac12\mathcal{K}{\left(0,2\right)}+\frac12\mathcal{K}{\left(1,2\right)}+F{\left(\frac{1}{\sqrt{3}}\right)}.\\ \end{align}

This answer by Ali Shadhar can be used (with slight modifications) to provide a closed-form expression for $$F$$:

$$F{\left(b\right)}=\frac{\pi}{2}\ln{\left(1+b\right)}-C+\operatorname{Ti}_{2}{\left(\frac{1-b}{1+b}\right)};~~~\small{0

where the inverse tangent integral $$\operatorname{Ti}_{2}$$ is defined by

$$\operatorname{Ti}_{2}{\left(z\right)}:=\int_{0}^{z}\mathrm{d}t\,\frac{\arctan{\left(t\right)}}{t};~~~\small{z\in\mathbb{R}},$$

and the Catalan constant $$C$$ may be given be the integral representation

$$C:=\operatorname{Ti}_{2}{\left(1\right)}=\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(t\right)}}{t}.$$

It is shown in this question that $$\mathcal{J}$$ obeys the following functional relation:

\begin{align} \mathcal{J}{\left(a,b\right)} &=\frac{\pi\arctan{\left(\sqrt{\frac{a-b}{(a+1)b}}\right)}}{2\sqrt{(a-b)b}}+\frac{\pi\arctan{\left(\sqrt{\frac{b}{(a+1)(a-b)}}\right)}}{2\sqrt{(a-b)b}}\\ &~~~~~-\frac{\arctan{\left(\frac{1}{\sqrt{a-b}}\right)}\arctan{\left(\frac{1}{\sqrt{b}}\right)}}{\sqrt{(a-b)b}}-\mathcal{J}{\left(a,a-b\right)};~~~\small{0

This has the implication that

$$a\,\mathcal{J}{\left(a,\frac{a}{2}\right)}=\pi\arctan{\left(\frac{1}{\sqrt{a+1}}\right)}-\arctan^{2}{\left(\sqrt{\frac{2}{a}}\right)};~~~\small{0

and in particular

$$2\,\mathcal{J}{\left(2,1\right)}=\pi\arctan{\left(\frac{1}{\sqrt{3}}\right)}-\arctan^{2}{\left(1\right)}=\frac{\pi^{2}}{6}-\frac{\pi^{2}}{16}=\frac{5\pi^{2}}{48}.$$

The value of $$\mathcal{K}{\left(1,a\right)}$$ can be related to $$\mathcal{K}{\left(0,a\right)}$$ through the recurrence relation

$$2\mathcal{K}{\left(1,a\right)}=-1+\frac{\pi}{2}\sqrt{a+1}+a\mathcal{K}{\left(0,a\right)};~~~\small{0

Proof: Given $$a>0$$,

\begin{align} 2\mathcal{K}{\left(1,a\right)} &=2\int_{0}^{1}\mathrm{d}t\,\sqrt{a+t^{2}}\arctan{\left(\sqrt{a+t^{2}}\right)}\\ &=\int_{0}^{1}\mathrm{d}t\,\frac{2(a+t^{2})}{\sqrt{a+t^{2}}}\arctan{\left(\sqrt{a+t^{2}}\right)}\\ &=\int_{0}^{1}\mathrm{d}t\,\frac{a}{\sqrt{a+t^{2}}}\arctan{\left(\sqrt{a+t^{2}}\right)}+\int_{0}^{1}\mathrm{d}t\,\frac{a+2t^{2}}{\sqrt{a+t^{2}}}\arctan{\left(\sqrt{a+t^{2}}\right)}\\ &=a\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{a+t^{2}}\right)}}{\sqrt{a+t^{2}}}+\int_{0}^{1}\mathrm{d}t\,\arctan{\left(\sqrt{a+t^{2}}\right)}\,\frac{d}{dt}\left[t\sqrt{a+t^{2}}\right]\\ &=a\mathcal{K}{\left(0,a\right)}+\int_{0}^{1}\mathrm{d}t\,\arctan{\left(\sqrt{a+t^{2}}\right)}\,\frac{d}{dt}\left[t\sqrt{a+t^{2}}\right]\\ &=a\mathcal{K}{\left(0,a\right)}+\sqrt{a+1}\arctan{\left(\sqrt{a+1}\right)}\\ &~~~~~-\int_{0}^{1}\mathrm{d}t\,t\sqrt{a+t^{2}}\frac{d}{dt}\arctan{\left(\sqrt{a+t^{2}}\right)};~~~\small{I.B.P.}\\ &=a\mathcal{K}{\left(0,a\right)}+\sqrt{a+1}\arctan{\left(\sqrt{a+1}\right)}\\ &~~~~~-\int_{0}^{1}\mathrm{d}t\,t\sqrt{a+t^{2}}\,\frac{t}{\left(1+a+t^{2}\right)\sqrt{a+t^{2}}}\\ &=a\mathcal{K}{\left(0,a\right)}+\sqrt{a+1}\arctan{\left(\sqrt{a+1}\right)}-\int_{0}^{1}\mathrm{d}t\,\frac{t^{2}}{1+a+t^{2}}\\ &=a\mathcal{K}{\left(0,a\right)}+\sqrt{a+1}\arctan{\left(\sqrt{a+1}\right)}+\int_{0}^{1}\mathrm{d}t\,\left[-1+\frac{1+a}{1+a+t^{2}}\right]\\ &=a\mathcal{K}{\left(0,a\right)}+\sqrt{a+1}\arctan{\left(\sqrt{a+1}\right)}-1+\sqrt{a+1}\arctan{\left(\frac{1}{\sqrt{a+1}}\right)}\\ &=-1+\frac{\pi}{2}\sqrt{a+1}+a\mathcal{K}{\left(0,a\right)}.\\ \end{align}

Then, for $$a>0$$ we find

\begin{align} \mathcal{K}{\left(0,a\right)} &=\mathcal{K}{\left(0,0\right)}+\int_{0}^{a}\mathrm{d}x\,\frac{d}{dx}\mathcal{K}{\left(0,x\right)}\\ &=\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(t\right)}}{t}+\int_{0}^{a}\mathrm{d}x\,\frac{d}{dx}\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{x+t^{2}}\right)}}{\sqrt{x+t^{2}}}\\ &=\operatorname{Ti}_{2}{\left(1\right)}+\int_{0}^{a}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\,\frac{\partial}{\partial x}\left[\frac{\arctan{\left(\sqrt{x+t^{2}}\right)}}{\sqrt{x+t^{2}}}\right]\\ &=C+\int_{0}^{a}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\,\left[\frac{1}{2\left(1+x+t^{2}\right)\left(x+t^{2}\right)}-\frac{\arctan{\left(\sqrt{x+t^{2}}\right)}}{2\left(x+t^{2}\right)^{3/2}}\right]\\ &=C+\int_{0}^{a}\mathrm{d}x\,\frac{\arctan{\left(\frac{1}{\sqrt{x+1}}\right)}-\arctan{\left(\sqrt{x+1}\right)}}{2x\sqrt{x+1}}\\ &=C+\int_{0}^{a}\mathrm{d}x\,\frac{\left[\frac{\pi}{2}-2\arctan{\left(\sqrt{x+1}\right)}\right]}{2x\sqrt{x+1}}\\ &=C+\int_{0}^{a}\mathrm{d}x\,\frac{\left[\frac{\pi}{4}-\arctan{\left(\sqrt{x+1}\right)}\right]}{x\sqrt{x+1}}\\ &=C+\int_{0}^{a}\mathrm{d}x\,\frac{\arctan{\left(\frac{1-\sqrt{x+1}}{1+\sqrt{x+1}}\right)}}{x\sqrt{x+1}}\\ &=C+\int_{1}^{\sqrt{a+1}}\mathrm{d}y\,\frac{2\arctan{\left(\frac{1-y}{1+y}\right)}}{\left(y^{2}-1\right)};~~~\small{\left[\sqrt{x+1}=y\right]}\\ &=C+\int_{0}^{\frac{1-\sqrt{a+1}}{1+\sqrt{a+1}}}\mathrm{d}t\,\frac{\arctan{\left(t\right)}}{t};~~~\small{\left[y=\frac{1-t}{1+t}\right]}\\ &=C-\int_{0}^{\frac{-1+\sqrt{a+1}}{1+\sqrt{a+1}}}\mathrm{d}t\,\frac{\arctan{\left(t\right)}}{t};~~~\small{\left[t\mapsto-t\right]}\\ &=C-\operatorname{Ti}_{2}{\left(\frac{\sqrt{a+1}-1}{\sqrt{a+1}+1}\right)}.\\ \end{align}

Putting it all together, we arrive finally at an expression for $$\mathcal{I}$$ in terms of well-known constants and functions:

\begin{align} \mathcal{I} &=-\frac12+\frac{\pi}{4}\ln{\left(\frac32\right)}-\frac{\pi^{2}}{8}+2\mathcal{J}{\left(2,1\right)}+\frac12\mathcal{K}{\left(0,2\right)}+\frac12\mathcal{K}{\left(1,2\right)}+F{\left(\frac{1}{\sqrt{3}}\right)}\\ &=-\frac12+\frac{\pi}{4}\ln{\left(\frac32\right)}-\frac{\pi^{2}}{8}+\frac{5\pi^{2}}{48}-\frac14+\frac{\pi\sqrt{3}}{8}+\mathcal{K}{\left(0,2\right)}+F{\left(\frac{1}{\sqrt{3}}\right)}\\ &=-\frac34-\frac{\pi^{2}}{48}+\frac{\pi\sqrt{3}}{8}+\frac{\pi}{4}\ln{\left(2+\sqrt{3}\right)}.\blacksquare\\ \end{align}

• Oh my, such a wonderful attempt, thank you for posting these. Commented Feb 25, 2022 at 6:10
• One may also look in to Ahmed' s Integral in Google. Commented Feb 3, 2023 at 14:35

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$$

• Well that’s certainly quicker… +1 Commented Feb 2, 2023 at 13:08

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.