Unit square integral $\int_{0}^{1} \int_{0}^{1} \left( - \frac{\ln(xy)}{1-xy} \right)^{m} dx dy $ In an article by Guillera and Sondow, one of the unit square integral identities that is proved (on p. 9) is: $$\int_{0}^{1} \int_{0}^{1}  \frac{\left(-\ln(xy) \right)^{s}}{1-xy}  dx dy = \Gamma(s+2) \zeta(s+2),$$ which holds for $\mathfrak{R}(s) > -1 $.
Let's define $$I_{m} := \int_{0}^{1} \int_{0}^{1} \left(  \frac{-\ln(xy) }{1-xy} \right)^{m}  dx dy  $$ for $m \in \mathbb{Z}_{\geq 1} $. Then WA finds $$I_{1} = -2 \zeta(3), \qquad I_{2} = 6 \zeta(3). $$ However, so far I haven't been able to find any closed-form evaluations for $I_{m}$ when $m \geq 3$, neither in the literature nor with CAS software.

Questions:

*

*What is $I_{m}$ for $m>2$ ?

*Does this family of definite integrals appear in the literature somewhere?


 A: We can start by considering the following integral:
$$\int_0^1\int_0^1 \frac{(xy)^a}{z-xy}dxdy=\frac{1}{z}\sum_{n=0}^\infty \int_0^1\int_0^1(xy)^a\left(\frac{xy}{z}\right)^ndxdy$$
$$=\sum_{n=0}^\infty \frac{1}{z^{n+1}}\int_0^1x^{n+a}dx\int_0^1y^{n+a}dy=\sum_{n=1}^\infty \frac{1}{z^{n+2}}\frac{1}{(n+a)^2}$$
Taking $m$ derivatives w.r.t. $a$ and setting $a$ to $0$ gives:
$$\int_0^1\int_0^1 \frac{(-\ln(xy))^m}{z-xy}dxdy=(m+1)!\operatorname{Li}_{m+2}\left(\frac{1}{z}\right)$$
Where $\operatorname{Li}_m(x)$ is the polylogarithm function.
Furthermore, we can take $m-1$ derivatives w.r.t. $z$ (and set it to $1$) in order to find:
$$\int_0^1\int_0^1 \frac{(-\ln(xy))^m}{(1-xy)^m}dxdy=(-1)^{m-1}m(m+1)\frac{d^{m-1}}{dz^{m-1}}\operatorname{Li}_{m+2}\left(\frac{1}{z}\right)\bigg|_{z=1}$$
And with the help of OEIS we can obtain the following closed form:
$$\boxed{\int_0^1\int_0^1 \left(\frac{-\ln(xy)}{1-xy}\right)^mdxdy=m(m+1)\sum_{k=1}^{m-1}|s(m-1,m-k)|\zeta(k+2)}$$
Where $s(n,m)$ is the Stirling number of the first kind.
In particular, we have:
\begin{array}{|c|c|}
\hline
I_2 & 6\zeta(3) \\ \hline
I_3 & 12(\zeta(3)+\zeta(4)) \\ \hline
I_4 & 20(\zeta(3)+3\zeta(4)+2\zeta(5)) \\ \hline
I_5 & 30(\zeta(3)+6\zeta(4)+11\zeta(5)+6\zeta(5)) \\ \hline
I_6 & 42(\zeta(3)+10\zeta(4)+34\zeta(5)+50\zeta(6)+24\zeta(7)) \\ \hline
I_7 & 56(\zeta(3)+15\zeta(4)+85\zeta(5)+225\zeta(6)+274\zeta(7)+120\zeta(8)) \\ \hline
\end{array}
A: \begin{align}
&\int_{0}^{1} \int_{0}^{1} \left( - \frac{\ln(xy)}{1-xy} \right)^{m}\  \overset{xy=t}{dy} \ dx \\
=& \int_0^1 \int_0^x \left( - \frac{\ln t}{1-t} \right)^{m}\frac1x \ dt\ dx
= \int_0^1 \int_t^1 \left( - \frac{\ln t}{1-t} \right)^{m}\frac1x \ dx\ dt\\
= &\int_0^1  \frac{(-\ln t)^{m+1}}{(1-t)^m} \ dt
\overset{t=e^{-u}}= \int_0^{\infty} \frac{u^{m+1} e^{-t}}{(1-e^{-u})^m}\ du\\
 =& \int_0^{\infty} {u^{m+1} e^{-t}} \sum_{k=0}^{\infty} \binom{m+k-1}{k} e^{-uk} du\\
=& \sum_{k=0}^{\infty} \binom{m+k-1}{k}\int_0^1 u^{m+1} e^{-(k+1)u}\ du
\\=&\sum_{k=0}^{\infty} \binom{m+k-1}{k} \frac{(m+1)!}{(k+1)^{m+2}}
\end{align}
