# Quadruple Integral $\int\limits_0^1\!\!\int\limits_0^1\!\!\int\limits_0^1\!\!\int\limits_0^1\frac1{1-xyzw}\,dw\,dz\,dy\,dx$

In page 122 of a book by William J. LeVeque, namely Topics in Number Theory (1956), there is an exercise for evaluating the following integral in two ways.

$$\int_0^1\!\!\!\int_0^1\frac1{1-xy}\,dy\,dx$$

First way is to write the integrand as a geometric series,

$$\int_0^1\!\!\!\int_0^1\frac1{1-xy}\,dy\,dx=\int_0^1\!\!\!\int_0^1\left(\sum_{n=1}^\infty(xy)^{n-1}\right)\,dy\,dx=\sum_{n=1}^\infty\frac1{n^2}$$

and the second way by use of a suitable change of variables ($y:=u-v,x:=u+v$) which is also published by Tom M. Apostol in this paper.

Hence, the second way together with the result of the first way is a proof for the famous Basel problem, in fact to show that $\sum_{n=1}^\infty\frac1{n^2}=\frac{\pi^2}6$.

Now, the main question is, if there is a suitable change of variables for the following integral

$$\int_0^1\!\!\!\int_0^1\!\!\!\int_0^1\!\!\!\int_0^1\frac1{1-xyzw}\,dw\,dz\,dy\,dx~?$$

Unfortunately, I think a similar change of variables like ($w:=p\pm q\pm r\pm s,\cdots$) doesn't work here, while I'm not really sure!

• Artin, I have made a very significant edit to this question, which I hopes makes your intentions much clearer. Feel free to roll it back if you do not like it. – Potato Aug 30 '13 at 17:20
• The Maple command $$with(VectorCalculus): printlevel := 15: int(1/(1-x*y*z*w), [x, y, z, w] = Parallelepiped(0 .. 1, 0 .. 1, 0 .. 1, 0 .. 1))$$ produces it, finding in sequence the antiderivatives $$-\frac {\ln \left( -wyz+1 \right) }{wyz},$$ $$dilog(-w*z+1)/(w*z) ,$$ and $$\frac {polylog(3, w)} w .$$ – user64494 Aug 30 '13 at 17:37
• what a cute integral – Saketh Malyala Jun 28 '17 at 6:21
• You may want to take a look at some substitutions to calculate $\int_0^1\int_0^1\frac{dxdy}{1-xy}$ elementarily. See e.g. Different methods to compute $\sum_{k=1}^\infty \frac{1}{k^2}$, the cited article by R. Chapman or this answer – punctured dusk Jun 28 '17 at 11:41
• I guess the reason the 2D proof is that the 2D version of the cube and the 2D version of the octahedron are equal. This is not true for other dimensions. Perhaps there is some useful rotation of the 4D cube that is relevant. – Akiva Weinberger Jun 28 '17 at 18:32

Regarding clarification of the question, Tom Apostol's Method of proving $\zeta(2)=\pi^2/6$ involves evaluating the double integral \begin{align*}\int_{0}^{1} \int_{0}^{1} \frac{1}{1-xy} \ dy \ dx,\end{align*} in two ways: converting the integrand into a geometric series and exchanging summation and integration to obtain $\zeta(2)$ and letting $x=u+v,y=u-v$ to get $\pi^2/6.$

We use Apostol's Double Integral Method to prove the result $\zeta(4)=\pi^4/90,$ using the quadruple integral \begin{align} \label{Apostol} \tag{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \frac{1}{1-xyzw} \ dw \ dz \ dy \ dx, \end{align} a similar geometric series argument shows \eqref{Apostol} is $\zeta(4).$ We now hope to evaluate \eqref{Apostol}.

Apostol's Method on $\zeta(4)$

We let $x=u+v,y=u-v, z=t/w.$ Then \eqref{Apostol} becomes \begin{align} \label{COV} \tag{2} \int_{\substack{0<u+v,u-v<1} } \int_{0}^{1} \int_{t}^{1} \frac{\frac{2}{w}}{1+(v^2-u^2)t} \ dw \ dt \ dv \ du. \end{align}

We split \eqref{COV} and consider the integral over the region $-u<v<u,0<u<1/2,$ that is,

\begin{align} \label{Part} \tag{3} \int_{0}^{\frac{1}{2}} \int_{-u}^{u} \int_{0}^{1} \int_{t}^{1} \frac{\frac{2}{w}}{1+(v^2-u^2)t} \ dw \ dt \ dv \ du= \int_{0}^{\frac{1}{2}} \int_{0}^{1} \int_{t}^{1} \int_{-u}^{u} \frac{\frac{2}{w}}{1+(v^2-u^2)t} \ dv \ dw \ dt \ du. \end{align} Integrating this, we find \eqref{Part} becomes \begin{align*} \int_{0}^{\frac{1}{2}} \int_{0}^{1} \int_{t}^{1} \frac{4\tan^{-1} \left(\frac{u\sqrt{t}}{\sqrt{1-tu^2}} \right)}{w\sqrt{t}\sqrt{1-tu^2}} \ dw \ dt \ du &=-\int_{0}^{\frac{1}{2}} \int_{0}^{1} \frac{4\log(t)\tan^{-1} \left(\frac{u\sqrt{t}}{\sqrt{1-tu^2}} \right)}{\sqrt{t}\sqrt{1-tu^2}} \ dt \ du \\ &=- \int_{0}^{1} \frac{2 \log(t) \left(\tan^{-1} \left(\sqrt{\frac{t}{4-t}}\right) \right)^2}{t} \ dt, \end{align*} after reversing the order of integration. Integrating by parts with $u=\left(\tan^{-1} \left(\sqrt{\frac{t}{4-t}}\right) \right)^2$ and $dv= \frac{-2 \log(t)}{t} \ dt$ and subsequently doing a change of a variables $t \mapsto x^2$ yields \eqref{Part} is equal to \begin{align} \label{Arctan} \tag{4} \int_{0}^{1} \frac{8\log^2(x)\tan^{-1} \left(\frac{x}{\sqrt{4-x^2}} \right)}{\sqrt{4-x^2}} \ dx. \end{align} It can be shown that \eqref{Arctan} upon expanding the integrand into a binomial series, is equal to \begin{align*} C(4)=\sum_{n=1}^{\infty} \frac{1}{n^4 \binom{2n}n}.\end{align*}

We wish to show $C(4)=17 \pi^4/3240.$ Expanding \eqref{Arctan} into the double integral \begin{align} \label{Double Integral} \tag{5} \int_{0}^{1} \int_{0}^{x} \frac{8\log^2(x)}{4-x^2+y^2} \ dy \ dx, \end{align} and applying the change of variables $x=u+v, y=u-v$ and $u \mapsto x, v \mapsto y,$ we get that \eqref{Double Integral} is \begin{align} \label{DI} \tag{6} \iint_{0<x-y<x+y<1} \frac{4\log^2(x+y)}{1-xy} \ dy \ dx. \end{align}

Integrating \eqref{DI} in the order presented, we split the region of integration into two triangular regions: $0<x<1/2, 0<y<x$ and $1/2<x<1,0<y<1-x$ and then use Mathematica and the integral additivity rule to get

\begin{align*}-\int_{0}^{1} \frac{4 \log^2(x)\log(x^2+1)}{x} \ dx + \int_{0}^{1} \frac{8 \log(x) Li_2 \left(\frac{x^2}{x^2+1} \right)}{x} \ dx - \int_{0}^{1} \frac{8 Li_3 \left(\frac{x^2}{x^2+1} \right)}{x} \ dx + \int_{\frac{1}{2}}^{1} \frac{8 Li_3 \left(\frac{x}{x^2+1} \right)}{x} \ dx + J, \end{align*} where $J$ is the sum of miscellaneous integrals. On the other hand, reversing the order of integration gives \eqref{DI} is \begin{align*}-J- \int_{0}^{\frac{1}{2}} \frac{8 Li_3 \left(\frac{y}{y^2+1} \right)}{y} \ dy \end{align*} Thus, equating these formulas, we see $C(4)$ is the sum of four integrals: \begin{align*} I_1 &=\int_{0}^{1} \frac{-2 \log^2(x)\log(x^2+1)}{x} \ dx \\ I_2 & = \int_{0}^{1} \frac{4 \log(x) Li_2 \left(\frac{x^2}{x^2+1} \right)}{x} \ dx \\ I_3 & = \int_{0}^{1} \frac{-4 Li_3 \left(\frac{x^2}{x^2+1} \right)}{x} \ dx \\ I_4 & =\int_{0}^{1} \frac{4 Li_3 \left(\frac{x}{x^2+1} \right)}{x} \ dx. \end{align*} Now we evaluate each integral.

To evaluate $I_1,$ we recall the series \begin{align*}\log(x^2+1)=\sum_{n=1}^{\infty} \frac{(-1)^{n-1} x^{2n}}{n}, \quad |x|<1.\end{align*} Interchanging summation and integration, integrating term by term, and finally recalling the Eta function \begin{align*}\eta(4)=\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n^4}=\frac{7}{8} \zeta(4),\end{align*} we see that \begin{align*}I_1=-\frac{7}{16} \zeta(4).\end{align*}

Using the substitution $u \mapsto \frac{x^2}{x^2+1}$, we get \begin{align} I_2 &=\int_{0}^{\frac{1}{2}} \frac{\log(\frac{u}{1-u}) Li_2 \left(u\right)}{u-u^2} \ du \nonumber \\ &= \int_{0}^{\frac{1}{2}} \frac{\log(1-u) \log^2\left( \frac{u}{1-u} \right)}{2u} \ du \label{I2IBP} \tag{7} \\ &= \int_{0}^{\frac{1}{2}} \frac{\log(1-u) \log^2(u)}{2u} \ du -\int_{0}^{\frac{1}{2}} \frac{\log(u) \log^2(1-u)}{u} \ du +\int_{0}^{\frac{1}{2}} \frac{ \log^3(1-u)}{u} \ du \label{splittedintegrals} \tag{8} \end{align} in which \eqref{I2IBP} follows from integrating by parts (with the function we differentiate being $Li_2(u)$) and \eqref{splittedintegrals} results from expanding the $log^2$ term. We will not evaluate the three integral in \eqref{splittedintegrals}; instead we will collect only the specific constants $Li_4(1), Li_4(-1),$ and constants that are rational multiples of $\pi^4$ The other constants are not necessary since they either will cancel out with one another or will cancel out with the constants resulting from the evaluation of $I_3.$ We recall
\begin{align*} Li_2(1) &= \zeta(2)= \frac{\pi^2}{6}, \\ Li_4(1) &=\zeta(4), \\ Li_4(-1) &=-\eta(4)=-\frac{7}{8}\zeta(4),\\ \left(Li_2 \left(\frac{1}{2} \right) \right)^2 &= \left(\frac{\zeta(2)}{2} - \frac{\log^2(2)}{2} \right)^2 = \frac{\pi^4}{144} - \frac{\pi^2}{12}\log^2(2) + \frac{\log^4(2)}{4}. \end{align*} Using Mathematica, the first integral in \eqref{splittedintegrals} does not consist of any constants we look for. The second integral, upon examining the antiderivative of the integrand yields the terms: \begin{align*}-2Li_4(1-u)+2Li_4(u)+2Li_4 \left( \frac{u}{-1+u} \right).\end{align*} Evaluating these terms at the endpoints yields the desired constant from the second integral in \eqref{splittedintegrals}, which is
$-\frac{7}{4}\zeta(4)+2\zeta(4)=\frac{1}{4}\zeta(4).$ The antiderivative of the integrand from the third integral in \eqref{splittedintegrals}, consisting of the term, \begin{align*}3 Li_4(1-u),\end{align*} yields the desired constant $-3\zeta(4)$ upon evaluation of the latter endpoint $u=0.$ Hence, for $I_2$ the desired constant is \begin{align*}-3\zeta(4)+\frac{1}{4}\zeta(4)=-\frac{11}{4}\zeta(4).\end{align*}

Using the substitution $u \mapsto \frac{x^2}{x^2+1}$, we get \begin{align} I_3 &= - \int_{0}^{\frac{1}{2}} \frac{2 Li_3(u)}{u-u^2} \ du \nonumber \\ &= - \int_{0}^{\frac{1}{2}} \frac{2Li_3(u)}{u} \ du - \int_{0}^{\frac{1}{2}} \frac{2Li_3(u)}{1-u} \ du \label{PF} \tag{9}, \end{align} in which \eqref{PF} is the result of using partial fractions. The first integral in \eqref{PF} does not consist of the aforementioned constants we seek. However, the integrand from the second integral in \eqref{PF} consists of the term \begin{align*} \left(Li_2(u) \right)^2,\end{align*} and upon evaluation of the endpoints, we get the constant $\pi^4/144.$

Lastly, for the fourth term, we apply $u \mapsto \frac{x}{x^2+1}$ to see \begin{align} I_4 & = \int_{0}^{\frac{1}{2}} \frac{4 Li_3(u)}{u\sqrt{1-4u^2}} \ du \nonumber \\ & = \sum_{n=1}^{\infty} \frac{\beta(n/2,n/2)}{n^3} \label{BS Conversion} \tag{10} \\ & = \sum_{n=1}^{\infty} \frac{\beta \left(\frac{2n}{2},\frac{2n}{2} \right)}{(2n)^3}+\sum_{n=1}^{\infty} \frac{\beta \left(\frac{2n-1}{2},\frac{2n-1}{2} \right)}{(2n-1)^3} \label{evenodd} \tag{11} \\ & = \frac{C(4)}{4} + \int_{0}^{1} \frac{\pi \log^2(u)}{\sqrt{4-u^2}} \ du \label{logsine} \tag{12} \\ & = \frac{C(4)}{4} + \frac{7 \pi^4}{216} \label{7pi^4/216} \tag{13}. \end{align} in which \eqref{BS Conversion} follows from recalling the definition of $Li_3(u)$, converting the integrand into a binomial series, and integrating term by term. Then the first term in \eqref{logsine} follows from simplifying the summand in the first term in \eqref{evenodd} with the properties of the beta function. The second term in \eqref{logsine} follows from converting the integrand into a binomial series, integrating term by term, and recognizing the expression corresponding to $\beta \left(\frac{2n-1}{2},\frac{2n-1}{2} \right).$ Lastly, the second term in \eqref{7pi^4/216} is the value of the integral in \eqref{logsine}, which we will not prove here (there are proofs available on StackExchange).

We already established $C(4)=I_1+I_2+I_3+I_4$ and asserted that only $\zeta(4)$ constants and constants that are rational multiples of $\pi^4$ remain from summing $I_2$ and $I_3.$ Hence, Rearranging terms yields the relation \begin{align} \label{sys1} \tag{14} \frac{3}{4} C(4)=-\frac{51}{16} \zeta(4) + \frac{17 \pi^4}{432} \end{align}

Now we obtain another system of equations relating $C(4)$ and $\zeta(4).$ This time we will reconsider $I_3$ and analyze the integral another way. Recall \begin{align*}I_3 = \int_{0}^{1} \frac{-4 Li_3 \left(\frac{x^2}{x^2+1} \right)}{x} \ dx.\end{align*} Integration by parts (with $dv=-4/x,$) yields

\begin{align} I_3 &= \int_{0}^{1} \frac{8 \log(x) Li_2 \left( \frac{x^2}{x^2+1} \right)}{x^3+x} \ dx \nonumber \\ & = \int_{0}^{1} \frac{8 \log(x) Li_2 \left( \frac{x^2}{x^2+1} \right)}{x} \ dx - \int_{0}^{1} \frac{8x \log(x) Li_2 \left( \frac{x^2}{x^2+1} \right)}{x^2+1} \ dx \label{PF2} \tag{15}\\ & = 2I_2- \int_{0}^{1} \frac{8x \log(x) Li_2 \left( \frac{x^2}{x^2+1} \right)}{x^2+1} \ dx \label{simplified} \tag{16}, \end{align}
in which \eqref{PF2} follows from partial fractions and \eqref{simplified} follows from the definition of $I_2.$ Now for the second term in \eqref{simplified}, we make the substitution $u \mapsto \frac{x}{x^2+1},$ to see it becomes \begin{align} \label{2nd term} \tag{17} - \int_{0}^{1} \frac{8x \log(x) Li_2 \left( \frac{x^2}{x^2+1} \right)}{x^2+1} \ dx & = - \int_{0}^{\frac{1}{2}} \frac{2 \log \left(\frac{u}{1-u} \right) Li_2(u)}{1-u} \ du \end{align} Integrating by parts on the right hand side of \eqref{2nd term} (with the differentiating function being $Li_2(u)$) yields \eqref{2nd term} is equal to \begin{multline} \label{split2} \tag{18} -\frac{7 \pi^4}{72} + \frac{7}{12} \pi^2 \log^2(2)+ \int_{0}^{\frac{1}{2}} \frac{\log(1-u) \log^2(u-1)}{u} \ du -\int_{0}^{\frac{1}{2}} \frac{2\log(u-1) \log^2(1-u)}{u} \ du \\ -\int_{0}^{\frac{1}{2}} \frac{2\log(1-u) Li_2(1-u)}{u} \ du \end{multline} Repeating the analysis of process of seeking the desired constants as done with \eqref{splittedintegrals}, we find the real part of the first integral in \eqref{split2} has a $\frac{3}{2} \zeta(4)$ constant, the second integral has a $12\zeta(4)$ constant, while the third has a $\frac{11}{8} \zeta(4)$ constant. Thus, \eqref{split2} has the desired constant \begin{align*}\frac{119}{8} \zeta(4)-\frac{7 \pi^4}{72},\end{align*} and from \eqref{simplified}, we see $I_3$ has the desired constant \begin{align*}\frac{119}{8} \zeta(4)- \frac{11}{2} \zeta(4)-\frac{7 \pi^4}{72}= \frac{75}{8} \zeta(4) -\frac{7 \pi^4}{72}. \end{align*} Hence, combining $I_1,\dots, I_4,$ we see that another relation arises: \begin{align} \label{sys2} \tag{19} \frac{3}{4}C(4)= \frac{99}{16} \zeta(4) - \frac{7 \pi^4}{108}. \end{align} Hence we have the system of equations \begin{align} \label{sys} \tag{20} \frac{3}{4} C(4)=-\frac{51}{16} \zeta(4) + \frac{17 \pi^4}{432}, \quad \frac{3}{4}C(4)= \frac{99}{16} \zeta(4) - \frac{7 \pi^4}{108}, \end{align} and solving this system yields $C(4)=17\pi^4/3240$ and $\zeta(4)=\pi^4/90.$

Remark To my knowledge, this method does not generalize any further to $\zeta(2n)$ for $n\geq 3.$ This is because the general version of the central binomial coefficient is not a rational multiple of $\pi^{2n}.$ The subtle point behind Apostol's method for $\zeta(2)$ is being able to evaluate the sum \begin{align*} \sum_{n=1}^{\infty} \frac{1}{\binom{2n}n n^2}= \frac{\pi^2}{18}, \end{align*} which is the result of evaluating $$\int_{0}^{1} \int_{0}^{1} \frac{1}{1-xy} \ dy \ dx$$ over the triangle $0<u<1/2, -u<v<u$ when making the change of variables $x=u+v, y=u-v.$

Note

$$\frac{1}{1-wxyz} = 1 + wxyz + (wxyz)^2+(wxyz)^3+\cdots$$

when $|wxyz|<1$, which is the case here. Integrating the series term by term gives the answer.

• @Artin Why do you need/want to do this with a change of variable? – Potato Aug 30 '13 at 16:52
• @Artin Could you please edit your question so that it describes what you found for $\zeta(2)$? I don't know what you're talking about. – Potato Aug 30 '13 at 16:58
• @Potato I think the question is how to get $\pi^4/90$ from this integral (and thus to show that $\zeta(4)=\pi^4/90$). – user8268 Aug 30 '13 at 17:03
• @user8268 Oh, I see. I've seen this kind of thing many times before, but I've never seen the change of variable for the $n=2$ case that Artin describes. I'd really like to see it. It might help find one for the $n=4$ case. – Potato Aug 30 '13 at 17:06
• This is not really an answer, since I don't see any change of variables! – user95733 Feb 11 '16 at 11:18

• This is not an answer, since this paper contains a suitable change of variables for some other integral representation of $\zeta(2n)$! – user95733 Feb 11 '16 at 11:22