Prove the following identity for the Apéry constant Perhaps this kind of integral is well knonw, or can be easily deduced from other. I don't know it but I would like to see the computation of this to refresh the computation of iterated integrals. I was inspired to change the factors of those integrals in the Wikipedia Page for Apéry 's constant. Firht I tried modify those integrals writting some factors $\sqrt{x}$. After with factors $1-x^2$ and computing with the online calculator of  Wolfram Alpha one find with the code

$\displaystyle{\int_0^1\int_0^1\int_0^1\int_0^1 (1-y^2)(1-t^2)/(1-xyz)
\,dx\,dy\,dz\,dt}$

an integral involving the cited constant.
I am assuming that these families were well known , because was easy to find those identities by similarity, when one modify such integrands as I've said. See also this MathWorld Page for the Apéry's constant, truly those formulas and identites are more complicated than my question, thus the kind of integral that 

Question. Show us how to prove 
  $$\int_0^1\int_0^1\int_0^1\int_0^1\frac{(1-y^2)(1-t^2)}{1-xyz}dxdydzdt=\frac{1}{36}(24\zeta(3)+9-2\pi^2).$$
   Thanks in advance.

Since the idea was easy, see comments for other example, it is easy that it was in the literature. 
Also you are welcome to answer the following optional question,

Question (Optional) Can you find a integral now involving at least a factor $\sqrt{\text{variable}}$ of previous variables $x,y,z,\ldots$, to get an identity, as the previous, for $\zeta(3)$?

 A: Because the $t$ integral is trivial, we will not bother with it and just calculate
$$
I=\int_{[0,1]^3} dx dy dz  \frac{1-y^2}{1-xyz}
$$
the first integration is standard an yields (let's take $z$)
$$
I=\int_{[0,1]^2} dx dy (1-y^2)\frac{\log(1-xy)}{xy}
$$
now, using one of the standard properties of Polylogarithms, namely $z \text{Li}_{\nu}'(z)=\text{Li}_{\nu-1}$, this becomes
$$
I=-\int_{[0,1]} dy (1-y^2)\frac{\text{Li}_2(y)}{y}
$$
Using the above mentioned property again this yields
$$
I=\text{Li}_3(1)+\int_{[0,1]} dy y\text{Li}_2(y)
$$
the last integral is easily tackeld using integration by parts with $y=u'$ and $v=\text{Li}_2(y)$. After putting back in the limits of integration this becomes
$$
I=\text{Li}_3(1)+\frac{1}{2}\text{Li}_2(1)+\frac{3}{8}
$$
which is equal to 

$$
I=\zeta(3)-\frac{\pi^2}{12}+\frac{3}{8}
$$

which gives us the announced result if we take a factor $2/3$ into account which stems from the integral over $t$
A: $\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{{#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{\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}
&\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}
{\pars{1 - y^{2}}\pars{1 -t^{2}} \over 1 - xyz}\,\dd x\,\dd y\,\dd z\,\dd t =
{2 \over 3}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}
{1 - y^{2} \over 1 - xyz}\,\dd x\,\dd y\,\dd z
\\[5mm] = &\
{2 \over 3}\sum_{n = 0}^{\infty}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}
\pars{x^{n}y^{n}z^{n} - x^{n}y^{n + 2}z^{n}}\,\dd x\,\dd y\,\dd z =
{2 \over 3}\sum_{n = 0}^{\infty}\bracks{{1 \over \pars{n + 1}^{3}} -
{1 \over \pars{n + 1}^{2}\pars{n + 3}}}
\\[5mm] = &\
{2 \over 3}\,\zeta\pars{3} +
\left.{2 \over 3}\,\totald{}{\mu}
\sum_{n = 0}^{\infty}{1 \over \pars{n + \mu}\pars{n + 3}}
\,\right\vert_{\ \mu\ =\ 1}
\\[5mm] = &\
{2 \over 3}\,\zeta\pars{3} +
{2 \over 3}\,\ \underbrace{%
\totald{}{\mu}\bracks{\Psi\pars{\mu} - \Psi\pars{3} \over \mu - 3}
_{\ \mu\ =\ 1}}_{\ds{\ {1 \over 24}\pars{9 - 2\pi^{2}}}}\qquad
\pars{~\Psi:\ Digamma Function~}
\\[5px]= &\ \ \bbox[#ffe,15px,border:1px dotted navy]
{\ds{{2 \over 3}\,\zeta\pars{3} + {1 \over 4} - {1 \over 18}\,\pi^{2}}}
\end{align}
A: EDIT:  I know this isn't as slick as the answer using polylogs but I always liked this approach to any integration with $\frac{1}{1-xyz}$ from [0,1].  Also I should note I got this approach from the book "Proofs from the Book"  in the section "Three Times $\frac{\pi^2}{6}$", I think is the name.  They use this approach to show $\zeta(2)=\frac{\pi^2}{6}$.
The t portion of the integral can be done separately using Fubini's theorem so noting that $\int_0^1(1-t^2)dt=t-\frac{t^3}{3}|_0^1=\frac{2}{3} $ this will be a prefactor that I will insert at the end. 
So the remaining integral is:
$\int_0^1 \int_0^1 \int_0^1 \frac{1-y^2}{1-xyz} dxdydz$.  We can rewrite this into two integrals:  $\int_0^1 \int_0^1 \int_0^1 \frac{1}{1-xyz} dxdydz-\int_0^1 \int_0^1 \int_0^1 \frac{y^2}{1-xyz} dxdydz$.  Now since we're working on $[0,1]$ we can use the following geometric series $\frac{1}{1-xyz}=\sum_{n=0}^\infty(xyz)^n$.  So the first integral becomes $\int_0^1 \int_0^1 \int_0^1 \sum_{n=0}^\infty (xyz)^n dxdydz$, exchanging the integrals and the sum and breaking up into three integrals we have, $\sum_{n=0}^\infty \int_0^1 x^ndx \int_0^1 y^n dy \int_0^1 z^n dz$ this is simply, $\sum_{n=0}^\infty (\frac{x^{n+1}}{n+1}|_0^1)(\frac{y^{n+1}}{n+1}|_0^1)(\frac{z^{n+1}}{n+1}|_0^1)$ = $\sum_{n=0}^\infty \frac{1}{(n+1)^3}$ or reindexing we have $\zeta(3)$.  The second integral can be tackled in the same way as the first changing everything into sums just as before we have now -$\sum_{n=0}^\infty \int_0^1 x^ndx \int_0^1 y^{n+2} dy \int_0^1 z^n dz$ or $-\sum_{n=0}^\infty (\frac{x^{n+1}}{n+1}|_0^1)(\frac{y^{n+3}}{n+3}|_0^1)(\frac{z^{n+1}}{n+1}|_0^1)=-\sum_{n=0}^\infty \frac{1}{(n+3)(n+1)^2}=-\sum_{n=1}^\infty \frac{1}{(n+2)n^2}$.  Using partial fractions this becomes: $\sum_{n=1}^\infty -\frac{1}{2n^2}-\frac{1}{4(n+2)}+\frac{1}{4n}=-\frac{1}{2} \zeta(2) - \frac{1}{4}(\sum_{n=1}^\infty\frac{1}{(n+2)}-\frac{1}{n})=-\frac{\pi^2}{12}+\frac{3}{8}$.  Where I used that I know what the zeta function is at 2 (it's $\frac{\pi^2}{6}$) and also that the last set of sums telescopes to $\frac{3}{2}$.  So grouping all this together we get $\zeta(3)-\frac{\pi^2}{12}+\frac{3}{8}$.  Multiplying this by $\frac{2}{3}$ from the t integral gives the desired result.  
