Some time ago I stumbled on an alleged identity $$\int\limits_0^\infty \frac{dx}{x} \int\limits_0^x \frac{dy}{y} \int\limits_0^y \frac{dz}{z} [\sin{x}+\sin{(xy)}\sin{(xz)}\sin{(xy+z)}]= \frac{\pi^3}{12}.$$ How this identity can be proved? The context under which this integral has emerged is described in Multiple Integral (American Mathematical Monthly problem 11621 and its generalization) .
I consider the function $$f(t):=\int_0^t\frac{dx}{x}\int_0^x\frac{dy}{y}\int_0^y\frac{dz}{z}\bigl\{ \sin x+\sin(xy)\sin(xz)\sin(xy+z)\bigr\}.$$ It has an asymptotic expansion with main terms $$f(t)= \frac{\cos t}{2}\frac{\log^2t}{t}+O\Bigl(\frac{\log t}{t}\Bigr).$$ Therefore, the limit $\lim_{t\to+\infty}f(t)=0\ne\pi^3/12$.
My proofs are too long to post them here, and I will only briefly explain the main points: The function $f(t)$ extends to an entire function with power series $$f(t)=\sum_{n=1}^\infty(1)^n \Bigl(\sum_{k=1}^{2n+1}\frac{H_{k1}}{k}\Bigr)\frac{t^{2n+1}}{(2n+1)!(2n+1)}$$ where $H_k$ are the harmonic numbers.
The power series can be transformed into an alternative integral representation $$f(t)=\int_0^tdu\int_0^1\int_0^1\Bigl(\frac{x\sin u}{u(1x)(1x y)}\frac{\sin(ux)}{u(1x)(1y)}+\frac{\sin(uxy)}{u(1y)(1xy)}\Bigr)\,dx\,dy.$$ This can be used to write $f(t)$ as a Fourier transform $$f(t)=\int_0^1\Bigl(\frac12\log^2(1x)+\sum_{n=1}^\infty\frac{(1x)^n1}{n^2}\Bigr) \frac{\sin(tx)}{x}\,dx.$$ Since the integrand is in $L^1(0,1)$, the RiemannLebesgue Lemma implies that $\lim_{t\to\infty}f(t)=0$, and the asymptotic expansion is obtained from this representation by standard methods.
The formulas above and the asymptotic expansion allow us to compute the function $f(t)$ easily.
The details can be found in http://arxiv.org/abs/1505.00440
We shall denote the integral with the letter $I$, i.e. :
$$I:=\int_0^\infty \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \int_0^y \frac{\mathrm{d}z}{z} [\sin{x}+\sin{(xy)}\sin{(xz)}\sin{(xy+z)}]$$
The key point is to insert a damping exponential $e^{\alpha z}$ (where $\alpha \geq 0$) into the integrand, this makes the integral directly dependent on the $\alpha$ factor, then apply the Leibnitz formula for differentatiaon under the integral sign, finally integrate again recalling $I(\infty)=0$ and $I(0)=I$. Albeit one can view on this procedure, respectively it can be done, using replacement $$\frac{1}{z}=\int_0^\infty \! e^{\alpha z} \; {\mathrm{d}\alpha}$$
Then, the change of the order of integration would imply this reformulation to appear :
$$I=\int_0^\infty \! {\mathrm{d}\alpha} \int_0^\infty \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \int_0^y {\mathrm{d}z} \, [\sin{x}+\sin{(xy)}\sin{(xz)}\sin{(xy+z)}] \, e^{\alpha z}$$
The inner integral is now directrly evaluable, holding the line :
$$I=\int_0^\infty \! {\mathrm{d}\alpha} \! \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \! \frac{1e^{\alpha y}}{\alpha(1+\alpha^2)} (\sin{x}+\sin{(xy)}) +\frac{1+e^{\alpha y}}{1+\alpha^2}(\cos{x}\cos{(xy)}) $$
As you can fine in this article (http://arxiv.org/pdf/1201.1975v1.pdf), integral similar to one contained in $I$ is computed there, namely :
$$\zeta{(2)} = \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} (\cos{(xy)}\cos{x}) $$
Althought, from representation of $I$ we are able to derive this result (denoted $\tilde{I}$ in the article) regardless of using the final one as a fact or from a different source. The actual approach is however very similar as the one described in the article : Just consider a parametric integral $\tilde{I}(\alpha)$ defined as:
$$\tilde{I}(\alpha) := \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{1+e^{\alpha y}}{y} (\cos{(xy)}\cos{x}) \, \mathrm{d}y $$
Evidently : $\tilde{I}(0)=2\tilde{I} \wedge \tilde{I}(\infty)=\tilde{I}$, using Fundamental Theorem of Calculus we can write $$\tilde{I}=\left( \tilde{I}2\tilde{I} \right)=\left( \tilde{I}(\infty)\tilde{I}(0) \right)=\int_0^\infty \tilde{I'}(\alpha) \, \mathrm{d}\alpha$$
Then differentiation under the integral sign and some tricky integration will lead us to the value $\zeta{(2)}$ as well.
Nevertheless, we will now continue in evaluation of $I$, the effort done on the discussion of the $\zeta{(2)}$ integral is not meaningless at all, because now we will separate this result from our integral representation of $I$, because of the relation : $$1+e^{\alpha y}=2(1e^{\alpha y})$$
we now split the $I$integral into two pieces :
$$I=J2\int_0^\infty \mathrm{d}\alpha \frac{\zeta{(2)}}{1+\alpha^2} =J\frac{\pi^3}{6}$$
where
$$J:=\int_0^\infty \! {\mathrm{d}\alpha} \! \int_0^\infty \! \frac{\mathrm{d}x}{x} \int_0^x \frac{\mathrm{d}y}{y} \! \frac{1e^{\alpha y}}{\alpha(1+\alpha^2)} (\sin{x}+\sin{(xy)})\frac{1e^{\alpha y}}{1+\alpha^2}(\cos{x}\cos{(xy)})$$
We are doing this for preparation to make an unified replacement (another), namely :
$$\frac{1e^{\alpha y}}{y}=\int_0^\alpha \! e^{\beta y} \; {\mathrm{d}\beta}$$
After change of the order of integration in $\beta$ :
$$J=\int_0^\infty \! \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \! \int_0^\alpha \! {\mathrm{d}\beta} \! \int_0^\infty \! \frac{\mathrm{d}x}{x} {\int_0^x} \mathrm{d}y \,[\sin{x}+\sin{(xy)}\alpha (\cos{x}\cos{(xy)})] e^{\beta y}$$
Again, by direct integration of inner integral :
$$J\!=\!\int_0^\infty \!\!\!\! \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \!\! \int_0^\alpha \!\! {\mathrm{d}\beta} \! \int_0^\infty \!\! \mathrm{d}x \frac{\sin{x}\!\!\alpha\cos{x}}{\beta}\frac{1\!\!e^{\beta y}}{x}\!+\!\frac{\alpha \beta\!1}{1\!+\!\beta^2}\frac{\cos{x}\!\!e^{\beta y}}{x}\!+\!\frac{\alpha\!+\!\beta}{1\!+\!\beta^2}\frac{\sin{x}}{x}$$
Resulting integrals are somewhat elementary (i.e. Dirichlet integral etc.), however we can use that dumb replacement rule again using
$$\frac{1}{x}=\int_0^\infty \! e^{\gamma x} \; {\mathrm{d}\gamma}$$
The other replacement however could also be done, i.e.
$$\frac{1e^{\beta x}}{x}=\int_0^\beta \! e^{\gamma x} \; {\mathrm{d}\gamma}$$
Define $$ \begin{split} S(\beta) & :=\int_0^\infty \mathrm{d}x \frac{1e^{\beta y}}{x}\sin{x} \\ C(\beta) & :=\int_0^\infty \mathrm{d}x \frac{1e^{\beta y}}{x}\cos{x} \\ L(\beta) & :=\int_0^\infty \mathrm{d}x \frac{\cos{x}e^{\beta y}}{x} \\ D & :=\int_0^\infty \mathrm{d}x \frac{\sin{x}}{x} \end{split}$$
By definition then
$$ J = \int_0^\infty \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \int_0^\alpha {\mathrm{d}\beta} \; \frac{1}{\beta}S(\beta)\frac{\alpha}{\beta} C(\beta)+\frac{\alpha \beta1}{1+\beta^2} L(\beta)+\!\frac{\alpha+\beta}{1+\beta^2} D $$
Now on the evaluation  for the first two integrals $S$ and $C$ the quickest way to compute them is via second type of replacement, The third an fourth one are done (due to the convergence issues) by the "dumb" replacement. For each integral we have: $$S(\beta)=\int_0^\beta {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x \sin{x} e^{\gamma x} = \int_0^\beta {\mathrm{d}\gamma} \frac{1}{1+\gamma^2}=\arctan{\beta} \\ \\ C(\beta)=\int_0^\beta {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x \cos{x} e^{\gamma x} = \int_0^\beta {\mathrm{d}\gamma} \frac{\gamma}{1+\gamma^2}=\frac{1}{2} \ln{(1+\beta^2)} \\ L(\beta)=\int_0^\infty {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x (\cos{x}\!\!e^{\beta x}) e^{\gamma x} = \int_0^\infty {\mathrm{d}\gamma} \frac{\gamma}{1\!+\!\gamma^2}\!\!\frac{1}{\beta\!+\!\gamma} = \ln{\frac{\sqrt{1 \!+\!\gamma^2}}{\beta \! + \! \gamma}} {\Big}_0^\infty = \ln{\beta} \\ D=\int_0^\infty {\mathrm{d}\gamma} \int_0^\infty \mathrm{d}x \sin{x} e^{\gamma x} = \int_0^\infty \frac{\mathrm{d}\gamma}{1+\gamma^2} = \frac{\pi}{2} $$
Substituting these results to the $J$ integral :
$$ J = \int_0^\infty \frac{\mathrm{d}\alpha}{\alpha(1+\alpha^2)} \int_0^\alpha {\mathrm{d}\beta} \; \frac{\arctan{\beta}}{\beta}\frac{\alpha}{2\beta} \ln{(1+\beta^2)}+\frac{\alpha \beta1}{1+\beta^2} \ln{\beta}+\!\frac{\alpha+\beta}{1+\beta^2} \frac{\pi}{2} $$
Next step will turn out to be really important. Just consider how far we are, we ended up with a double integral with its domain of integration  infinite triangle, i.e. half of the first quadrant. Let us perform a variable change in such manner it transforms to rectangle. There are many candidates (e.g. polar one), however we chose this particular :
$$\alpha=t \\ \beta=tk$$
With Jacobian of transformation equal to $J=t$, domain itself now consists of a semiinfinite strip $$(t,k)\in\{(0,\infty)\times(0,1)\}$$
After this transformation :
$$ J \!=\! \int_0^1 \!\mathrm{d}k \int_0^\infty \!\mathrm{d}t \; \frac{\arctan{kt}}{kt(1+t^2)}\frac{\ln{(1\!+\!k^2t^2)}}{2k(1\!+\!t^2)}\!\!\frac{1kt^2}{(1\!+\!t^2)(1\!+\!k^2t^2)} \ln{(kt)}\!+\!\!\frac{t(1+k)}{(1\!+\!t^2)(1\!+\!k^2t^2)} \frac{\pi}{2} $$
Defining ($0<s<4$) :
$$ \begin{split} J_1(k) & :=\int_0^\infty \frac{\arctan{kt}}{t(1+t^2)} \mathrm{d}t \\ J_2(k) & :=\int_0^\infty \frac{\ln{(1\!+\!k^2t^2)}}{2(1+t^2)} \mathrm{d}t \\ P_s(k) & :=\int_0^\infty \frac{t^{s1}\,\mathrm{d}t}{(1\!+\!t^2)(1\!+\!k^2t^2)} \\ L_s(k) & :=\int_0^\infty \frac{t^{s1}\ln{t} \, \mathrm{d}t}{(1+t^2)(1+k^2t^2)} \end{split} $$
Then by direct definition :
$$ J \!=\! \int_0^1 \!\mathrm{d}k \; \frac{J_1(k)}{k}\frac{J_2(k)}{k}\!P_1(k) \ln{k}+ P_3(k) k \ln{k}  L_1(k) + L_3(k) k + \frac{\pi}{2} (1+k)P_2(k) $$
Let us compute the $P_s(k)$ value only (for arbitrary $s$), note that :
$$\frac{1}{(1+t^2)(1+k^2t^2)}=\frac{1}{1k^2}\frac{1}{1+t^2}\frac{k^2}{1k^2}\frac{1}{1+k^2t^2}$$
Substituting that decomposition to the definition of $P_s(k)$ and by inserting the substitution $tk \rightarrow t$ into the second integral we can get the $s$dependency out of both integrals, recollecting :
$$P_s(k) =\frac{1k^{2s}}{1k^2} \int_0^\infty \frac{t^{s1}}{1+t^2} \mathrm{d}t$$
Unfortunately this representation restricts the value of $s$ greatly, now $0<s<2$
Recall the definition of the beta function, using a simple tranformation $1+t^2 \rightarrow 1/t$ and applying the reflection formula one gets:
$$ P_s(k) =\frac{1k^{2s}}{1k^2} \frac{\pi}{2 \sin{(\frac{\pi s}{2})}} $$
For $s=1$ we have :
$$P_1(k) = \frac{\pi}{2} \frac{1}{1+k}$$
Using the L'Hopital's rule also extraction of $P_2(k)$ is possible :
$$P_2(k) =  \frac{\ln{k}}{1k^2}$$
From the definition is clear that $P_1(k)+k^2P_3(k)=\int_0^\infty \! \frac{\mathrm{d}t}{1+t^2} = \pi/2$, therefore :
$$P_3(k) = \frac{\pi}{2k} \frac{1}{1+k}$$
Amazingly this result could be obtained also by evaluating the $P_s$ fromula at $s=3$
Evidently $L_s(k)=\frac{\mathrm{d}}{\mathrm{d}s}P_s(k)$, performing that differentiation on the precomputed $P_s(k)$ one gets :
$$ L_s(k) =\frac{k^{2s} \ln{k}}{1k^2} \frac{\pi}{2 \sin{(\frac{\pi s}{2})}}  \frac{\pi^2}{4} \frac{1k^{2s}}{1k^2} \frac{\cos{(\frac{\pi s}{2})}}{\sin^2{(\frac{\pi s}{2})}} $$
Evaluation at $s=1$ gives:
$$L_1(k) = \frac{\pi}{2} \frac{k \ln{k}}{1k^2}$$
Using the same trick as in $P_3$ one can also evalute $L_3$ (However direct substitution also could have been done)
$$L_3(k) =  \frac{\pi}{2k} \frac{\ln{k}}{1k^2}$$
Let us now make a closer look at $J_1$ and $J_2$ integrals, as we can see, both at $k=0$ vanish, i.e. $J_1(0)=J_2(0)=0$
It seems natural to apply differentiation under the integral sign, performing that one can obtain :
$$J'_1(k) = \int_0^\infty \frac{\mathrm{d}t}{(1+t^2)(1+k^2t^2)} = P_1(k) = \frac{\pi}{2} \frac{1}{1+k} \\ J'_2(k) = k\int_0^\infty \frac{t^2 \mathrm{d}t}{(1+t^2)(1+k^2t^2)} = k P_3(k) = \frac{\pi}{2} \frac{1}{1+k} $$
Using Fundamental theorem of calculus $f(x)=f(0)+\int_0^x f'(y) \mathrm{d}y$, resulting in :
$$J_1(k)=J_2(k)= \frac{\pi}{2}\ln{(1+k)}$$
Putting all these informations ($J_s,P_s,L_s$) togehter, in $J$ we have now, due to "huge cancelerin" :
$$ J = \pi \int_0^1 \! \frac{\ln{k}}{1k} \, \mathrm{d}k = \pi \mathrm{Li}_2(1k) {\Big}_0^1 = \frac{\pi^3}{6} $$
Tracing back to the original definition of $J$ gives the desired result :) ...

$\begingroup$ Quite impressive "physicists'" style computation! Juan's "mathematicians'" style computation is also beautiful. I wonder where I made errors in my calculation which can be found at wwwsnd.inp.nsk.su/~silagadz/I4.pdf Will you get the same result zero if the integration limits were $$\int\limits_0^\infty dx \int\limits_x^\infty dy \int\limits_y^\infty dz?$$ $\endgroup$ May 19 '15 at 15:37

$\begingroup$ I have to apologize for adding and answer so late, I was in the meantime studying at uni so I forgot where this page was until now, the answer is I think according to same approach that the integral diverges to infinity. The terms will not cancel out. $\endgroup$ Feb 16 '18 at 22:01

$\begingroup$ Thank! However I don't understand now what was my mistake in wwwsnd.inp.nsk.su/~silagadz/I4.pdf In fact I'm interesting to evaluate (1) from this manuscript in direct way without physics input. Maybe you can tackle this integral too? $\endgroup$ Feb 19 '18 at 6:57