Unit square, Riemann zeta function: $\int_0^1\int_0^1\frac{(-\log(xy))^s}{1-xy}\,dx\,dy = \Gamma(s+2) \zeta(s+2), \operatorname{Re}(s)>1$ I want to prove following integral: For $\operatorname{Re}(s)>1,$
$$
\int_0^1 \int_0^1 \frac{(-\log(xy))^s}{1-xy}\,dx\,dy = \Gamma(s+2)\,\zeta(s+2). 
$$
I am having trouble with proving this equation.
What can be a good parametrization for this integral?  My trial was treating $t = -\log(xy),$ but this does not work well...
 A: Consider it as an iterated integral, substitute $x=z/y$ and integrate by parts:
$$\int_0^1\int_0^1\frac{(-\log xy)^s}{1-xy}\,dx\,dy=\int_0^1\int_0^y\frac{(-\log z)^s}{1-z}\,dz\,\frac{dy}{y}\\=\underbrace{\log y\left.\int_0^y\frac{(-\log z)^s}{1-z}\,dz\,\right|_{y=0}^{y=1}}_{=0}+\int_0^1\frac{(-\log y)^{s+1}}{1-y}\,dy\\=\int_0^1\frac{(-\log y)^{s+1}}{1-y}\,dy\color{gray}{\left[=\int_0^\infty\frac{x^{s+1}\,dx}{e^x-1}\right]}=\Gamma(s+2)\zeta(s+2).$$
A: Since it's highly probable that one will use this formula often for $s\in N$, here's an approach for this case (underneath is a more general one).
Consider the following integral:
$$I(n)=\int_0^1 \int_0^1 (xy)^{n-1} dxdy=\frac{1}{n^2}$$
By taking $s$ derivatives with respect to $n$ and multiplying both sides by $(-1)^s$ we get:
$$\int_0^1\int_0^1 (xy)^{n-1}(-\ln(xy))^sdxdy=\frac{\Gamma(s+2)}{n^{s+2}}$$
Finally, notice that:
$$\sum_{n=1}^\infty \int_0^1\int_0^1 (xy)^{n-1}(-\ln(xy))^sdxdy=\int_0^1\int_0^1 \frac{(-\ln(xy))^s}{1-xy}dxdy$$
$$=\sum_{n=1}^\infty \frac{\Gamma(s+2)}{n^{s+2}}=\Gamma(s+2) \zeta(s+2) 
$$

$$\begin{align}
\int_0^1\int_0^1 \frac{(-\ln(xy))^{s}}{1-xy}\,\mathrm{d}x\,\mathrm{d}y
&\overset{xy=t}=\int_0^1\int_0^y \frac{(-\ln t)^s}{y(1-t)}\,\mathrm{d}t\,\mathrm{d}y\\
&=\int_0^1 \int_t^1\frac{(-\ln t)^s}{y(1-t)}\,\mathrm{d}y\,\mathrm{d}t\\
&=\int_0^1 \frac{(-\ln t)^{s+1}}{1-t}\,\mathrm{d}t \\
&\overset{-\ln t= x}=\int_0^\infty \frac{x^{n+1}e^{-x}}{1-e^{-x}}\,\mathrm{d}x\\
&=\sum_{n=0}^\infty \int_0^\infty x^{n+1} e^{-(n+1)x}\,\mathrm{d}x\\
&=\sum_{n=0}^\infty \frac{1}{(n+1)^{s+2}} \int_0^\infty x^{n+1} e^{-x}dx \\
&=\sum_{n=1}^\infty \frac{\Gamma(s+2)}{n^{s+2}}\\
&=\Gamma(s+2)\zeta(s+2)\\
\end{align}$$
