A triple integral dancing in the unit cube

Straight integration seems pretty tedious and difficult, and I guess the symmetry may open some new ways of which I'm not aware. What would your idea be?

Calculate $$\int_0^1 \int_0^1 \int_0^1 \frac{x^2}{\sqrt{x^2+1} \left(x^2-y^2\right) \left(x^2-z^2\right)}+\frac{y^2}{\sqrt{y^2+1} \left(y^2-x^2\right) \left(y^2-z^2\right)}+\frac{z^2}{\sqrt{z^2+1} \left(z^2-x^2\right) \left(z^2-y^2\right)} \, dx \ dy \ dz.$$

A 300 points bounty moment: After 2 years and 10 months since the problem was posed no full solution was provided yet. Is it possible to find a slick solution (like a bolt of lightning)? Good luck!

Let us set $A=x^2,B=y^2,C=z^2$ and: $$ S_n(A,B,C) = \frac{A^n}{(A-B)(A-C)}+\frac{B^n}{(B-A)(B-C)}+\frac{C^n}{(C-A)(C-B)} .\tag{1}$$ By partial fraction decomposition, it is straightforward to check that $S_0=S_1=0$ and $S_2=1$.

For every $n>2$ induction gives: $$ S_n(A,B,C) = h_{n-2}(A,B,C)\tag{2} $$ where $h_k$ is a complete homogeneous symmetric polynomial, satisfying: $$ \sum_{k\geq 0}h_k(A,B,C)t^k=\frac{1}{(1-A t)(1-B t)(1- Ct)}.\tag{3} $$ Assume now that the Taylor series of $g(u)=\frac{u}{\sqrt{1+u}}$ around $u=0$ is given by $\sum_{n\geq 0}g_n u^n$.

Our integral is given by: $$ \iiint_{(0,1)^3}\sum_{n\geq 0}g_n\,S_n(x^2,y^2,z^2)\,d\mu=g_2+\iiint_{(0,1)^3}\sum_{n\geq 3}g_n\cdot h_{n-2}(x^2,y^2,z^2)\,d\mu\tag{4} $$

Now there are two possible approaches. The first approach is to look for a linear operator that maps $t^k$ to $g_k$, then exploit $(3)$ in order to write the RHS of $(4)$ as a simpler integral. The second approach is to apply a double-counting argument on the monomials appearing in the RHS of $(4)$. Since $\iiint_{(0,1)^3}x^{2a}y^{2b}z^{2c}\,d\mu=\frac{1}{(2a+1)(2b+1)(2c+1)}$, we get:

$$ \iiint_{(0,1)^3}\sum_{n\geq 3}g_n\cdot h_{n-2}(x^2,y^2,z^2)\,d\mu =\\= \sum_{\substack{a,b,c\geq 0\\(a,b,c)\neq (0,0,0)}}\frac{g_{a+b+c+2}}{(2a+1)(2b+1)(2c+1)}\tag{5} $$ and the last series can be further simplified by computing: $$ l_n=\sum_{\substack{a,b,c\geq 0\\a+b+c=n}}\frac{1}{(2a+1)(2b+1)(2c+1)}\tag{6} $$ then computing $\sum_{n\geq 1}l_n g_n$, that is obviously related with the integral: $$\iiint_{(0,1)^3}g(x^2 y^2 z^2)\,d\mu = \iiint_{(0,1)^3}\frac{x^2 y^2 z^2}{\sqrt{1+x^2 y^2 z^2}}\,d\mu.\tag{7}$$

2018 update. Following the first approach, the given integral equals the opposite of $$\iiint_{(0,1)^3}\frac{2}{\pi}\int_{0}^{\pi/2}\frac{\sin^2(\theta)\,d\theta}{(1+x^2\sin^2\theta)(1+y^2\sin^2\theta)(1+z^2\sin^2\theta)}\,dx\,dy\,dz$$ which by Fubini's theorem boils down to $$ \frac{2}{\pi}\int_{0}^{\pi/2}\frac{\arctan^3(\sin\theta)}{\sin\theta}\,d\theta = \frac{2}{\pi}\int_{0}^{\pi/2}\frac{\arctan^3(\cos\theta)}{\cos\theta}\,d\theta=\frac{4}{\pi}\int_{0}^{1}\frac{\arctan^3\left(\frac{1-t^2}{1+t^2}\right)}{1-t^2}\,dt $$ then to $\frac{2}{\pi}\int_{0}^{1}\frac{\arctan^3\left(\frac{1-t}{1+t}\right)}{\sqrt{t}(1-t)}\,dt = \frac{2}{\pi}\int_{0}^{1}\frac{\left(\frac{\pi}{4}-\arctan t\right)^3}{\sqrt{t}(1-t)}\,dt = \frac{2}{\pi}\int_{0}^{+\infty}\left(\frac{\pi}{4}-\arctan\tanh^2\frac{u}{2}\right)^3\,du$ which is related to the Gudermannian function. Probably it can be expressed in a closed form through high order polylogarithms / Euler sums with high weight. Accurate numerical approximation are pretty simple since the last integrand function essentially behaves like $\frac{\pi^3}{64}\exp\left(-\frac{3}{\pi}u^2\right)$ on $\mathbb{R}^+$.

We may also notice that the two approaches are pretty obviously equivalent, since $l_n$ appearing in $(6)$ is clearly given by a coefficient of the MacLaurin series of $\text{arctanh}^3(z)$.

The value of your triple integral is

$$\color{royalblue}{-\frac{3}{16}\text{Li}_3\left(\frac{1}{2}\left(1-\sqrt{2}\right)\right)-\frac{3\text{Li}_3\left(2-\sqrt{2}\right)}{16}+\frac{9}{16}\text{Li}_3\left(2\left(-1+\sqrt{2}\right)\right)-\frac{3}{16}\text{Li}_3\left(\frac{1}{2+\sqrt{2}}\right)+\frac{3}{16}\text{Li}_3\left(\frac{1}{4}\left(2+\sqrt{2}\right)\right)+\frac{3}{16}\text{Li}_3\left(\frac{1}{2\left(2+\sqrt{2}\right)}\right)-\frac{3}{64}\text{Li}_3\left(3-2\sqrt{2}\right)+\frac{3}{64}\text{Li}_3\left(17-12\sqrt{2}\right)-\frac{3}{16}\text{Li}_2\left(\frac{1}{2}\left(1-\sqrt{2}\right)\right)\ln2-\frac{9}{32}\text{Li}_2\left(-1+\sqrt{2}\right)\ln2-\frac{9}{16}\text{Li}_2\left(2\left(-1+\sqrt{2}\right)\right)\ln2+\frac{3}{8}\text{Li}_2\left(\frac{1}{4}\left(2+\sqrt{2}\right)\right)\ln2+\frac{3}{8}\text{Li}_2\left(\frac{1}{\sqrt{2}}\right)\ln\left(1+\sqrt{2}\right)+\frac{3}{16}\text{Li}_2\left(-1+\sqrt{2}\right)\ln\left(1+\sqrt{2}\right)-\frac{3}{16}\text{Li}_2\left(\frac{1}{2+\sqrt{2}}\right)\ln\left(2\left(2+\sqrt{2}\right)\right)-\frac{3}{16}\text{Li}_2\left(\frac{1}{4}\left(2+\sqrt{2}\right)\right)\ln\left(2+\sqrt{2}\right)+\frac{3}{16}\text{Li}_2\left(\frac{1}{2\left(2+\sqrt{2}\right)}\right)\ln\left(2\left(2+\sqrt{2}\right)\right)+\frac{3}{8}\text{Li}_2\left(-3+2\sqrt{2}\right)\ln\left(1+\sqrt{2}\right)-\frac{3}{16}\text{Li}_2\left(3-2\sqrt{2}\right)\ln\left(1+\sqrt{2}\right)-\frac{21\zeta(3)}{512}+\frac{73\ln^32}{128}-\frac{7}{16}\ln^3\left(1+\sqrt{2}\right)-\frac{1}{16}\ln^3\left(2+\sqrt{2}\right)+\frac{39}{32}\ln\left(1+\sqrt{2}\right)\ln^22-\frac{45}{64}\ln\left(2+\sqrt{2}\right)\ln^22-\frac{21}{32}\ln^2\left(1+\sqrt{2}\right)\ln2+\frac{15}{32}\ln^2\left(2+\sqrt{2}\right)\ln2+\frac{3}{16}\ln\left(1+\sqrt{2}\right)\ln^2\left(2+\sqrt{2}\right)+\frac{3}{16}\ln^2\left(1+\sqrt{2}\right)\ln\left(2+\sqrt{2}\right)+\frac{47}{256}\pi^2\ln2-\frac{9}{16}\ln\left(1+\sqrt{2}\right)\ln\left(2+\sqrt{2}\right)\ln2+\frac{3}{16}\pi^2\ln\left(1+\sqrt{2}\right)-\frac{37}{128}\pi^2\ln\left(2+\sqrt{2}\right)}$$

Here is the Mathematica code of the result, I hope someone can simplify it further.

47/256\[Pi]^2Log[2]+(73Log[2]^3)/128+3/16\[Pi]^2Log[1+Sqrt[2]]+39/32Log[2]^2Log[1+Sqrt[2]]-21/32Log[2]Log[1+Sqrt[2]]^2-7/16Log[1+Sqrt[2]]^3-37/128\[Pi]^2Log[2+Sqrt[2]]-45/64Log[2]^2Log[2+Sqrt[2]]-9/16Log[2]Log[1+Sqrt[2]]Log[2+Sqrt[2]]+3/16Log[1+Sqrt[2]]^2Log[2+Sqrt[2]]+15/32Log[2]Log[2+Sqrt[2]]^2+3/16Log[1+Sqrt[2]]Log[2+Sqrt[2]]^2-1/16Log[2+Sqrt[2]]^3+3/8Log[1+Sqrt[2]]PolyLog[2,1/Sqrt[2]]-3/16Log[1+Sqrt[2]]PolyLog[2,3-2Sqrt[2]]-3/16Log[2]PolyLog[2,1/2(1-Sqrt[2])]-9/32Log[2]PolyLog[2,-1+Sqrt[2]]+3/16Log[1+Sqrt[2]]PolyLog[2,-1+Sqrt[2]]-9/16Log[2]PolyLog[2,2(-1+Sqrt[2])]+3/16Log[2(2+Sqrt[2])]PolyLog[2,1/(2(2+Sqrt[2]))]-3/16Log[2(2+Sqrt[2])]PolyLog[2,1/(2+Sqrt[2])]+3/8Log[2]PolyLog[2,1/4(2+Sqrt[2])]-3/16Log[2+Sqrt[2]]PolyLog[2,1/4(2+Sqrt[2])]+3/8Log[1+Sqrt[2]]PolyLog[2,-3+2Sqrt[2]]+3/64PolyLog[3,17-12Sqrt[2]]-3/64PolyLog[3,3-2Sqrt[2]]-3/16PolyLog[3,1/2(1-Sqrt[2])]-3/16PolyLog[3,2-Sqrt[2]]+9/16PolyLog[3,2(-1+Sqrt[2])]+3/16PolyLog[3,1/(2(2+Sqrt[2]))]-3/16PolyLog[3,1/(2+Sqrt[2])]+3/16PolyLog[3,1/4(2+Sqrt[2])]-(21Zeta[3])/512

I will continue from an integral representation given by Jack D'Aurizio $$\int_0^1 {\frac{{{{\arctan }^3}x}}{{x\sqrt {1 - {x^2}} }}dx} = \int_0^{\pi/4} {\frac{{{x^3}}}{{\sin x\sqrt {\cos 2x} }}dx} = \frac{1}{8}\int_{ - \pi /2}^{\pi /2} {\frac{{{x^3}\cos \frac{x}{2}}}{{\sin x\sqrt {\cos x} }}dx} $$

Next assume $a,b,c \in \mathbb{R}$, not necessarily integer with $a-b-c>0$.

$$\begin{aligned} &\int_{ - \pi /2}^{\pi /2} {{e^{iax}}{{( - 2i\sin x)}^b}{{(2\cos x)}^c}dx} = - i\int_{ - \pi /2}^{\pi /2} {{e^{i(a - 1)x}}{{({e^{ - ix}} - {e^{ix}})}^b}{{({e^{ix}} + {e^{ - ix}})}^c}d({e^{ix}})} \\ &= - i\int_C {{z^{a - 1}}{{\left( {\frac{1}{z} - z} \right)}^b}{{\left( {z + \frac{1}{z}} \right)}^c}dz} = - i\int_{ - i}^i {{z^{a - 1}}{{\left( {\frac{1}{z} - z} \right)}^b}{{\left( {z + \frac{1}{z}} \right)}^c}dz} \\&= {e^{ - \frac{\pi }{2}i(a - b - c - 1)}}\int_0^1 {{u^{a - 1}}{{\left( {u + \frac{1}{u}} \right)}^b}{{\left( {\frac{1}{u} - u} \right)}^c}du} + {e^{\frac{\pi }{2}i(a - b - c - 1)}}\int_0^1 {{u^{a - 1}}{{\left( {u + \frac{1}{u}} \right)}^b}{{\left( {\frac{1}{u} - u} \right)}^c}du} \\ &=\sin \left( {\frac{\pi }{2}(a - b - c)} \right)\int_0^1 {{x^{(a - b - c)/2 - 1}}{{(1 - x)}^c}{{(1 + x)}^b}dx} \end{aligned}$$ where $C$ is the semicircle with ends at $\pm i$ and lying in the right half plane. You may want to convince yourself that all details of the contour integration is justified.

Take $a=1/2,b=-1,c=-1/2$, $$\int_{ - \pi /2}^{\pi /2} {{e^{iax}}\frac{1}{{( - 2i\sin x)\sqrt {2\cos x} }}dx} = \sin \left( {\frac{\pi }{2}(a + \frac{3}{2})} \right)\int_0^1 {{x^{\frac{a}{2} - \frac{1}{4}}}{{(1 - x)}^{ - \frac{1}{2}}}{{(1 + x)}^{ - 1}}dx} $$ note that integral does not converge, but we can remedy it by first differentiating with respect to $a$, then let $b \to 1$. In any case, we have $$\frac{1}{{\sqrt 2 }}\int_0^{\pi /2} {\frac{{{x^3}\cos \frac{x}{2}}}{{\sin x\sqrt {\cos x} }}dx} = - \frac{{3\pi }}{8}\color{red}{\int_0^1 {\frac{{{{\ln }^2}x}}{{\sqrt {1 - x} (1 + x)}}dx}} + \frac{\pi^3}{8}\underbrace{\int_0^1 {\frac{1}{{\sqrt {1 - x} (1 + x)}}dx}}_J $$ Thanks to the fact that $a-b-c = 2$ in our situation, so we don't need to handle integral involving $\ln^3 x$: the red integral already poses major hindrance.

Note that $J = \sqrt{2} \ln(1+\sqrt{2})$, for the red one: $$I = \int_0^1 {\frac{{{{\ln }^2}(1 - x)}}{{\sqrt x (2 - x)}}dx} = 2\int_0^1 {\frac{{{{\ln }^2}(1 - {x^2})}}{{2 - {x^2}}}dx} $$ so it boils down to the four integrals: $$\int_0^1 {\frac{{{{\ln }^2}(1 - x)}}{{2 - {x^2}}}dx} \quad \int_0^1 {\frac{{{{\ln }^2}(1 + x)}}{{2 - {x^2}}}dx} \quad \int_0^1 {\frac{{\ln (1 + x)\ln (1 - x)}}{{\sqrt 2 - x}}dx} \quad \int_0^1 {\frac{{\ln (1 + x)\ln (1 - x)}}{{\sqrt 2 + x}}dx} $$ each one is straightforward, albeit tedious evaluation of trilogarithm/dilogarithm. Combining them gives the result.

A note: The integral is only apparently improper. Substituting $$x={1\over2}\left(u-{1\over u}\right)\quad(1\leq u\leq1+\sqrt{2}),\quad\sqrt{1+x^2}={1\over2}\left(u+{1\over u}\right),\quad dx={1\over2}\left(1+{1\over u^2}\right)\>du\ ,$$ and similarly for the other variables leads to an integrand with a polynomial in the numerator and $$(u+v)(1+uv)(v+w)(1+vw)(w+u)(1+wu)$$ in the denominator. Mathematica computes the following numerical value of the integral (in $(u,v,w)$-terms as well as in the original form): $$-0.287162469\ .$$