Prove $\int_{0}^{1}\int_{0}^{1}[-\ln(xy)]^s\left(\frac{1}{\ln(xy)}+\frac{1}{1-xy}\right)dxdy=\Gamma(s+2)\left[\zeta(s+2)-\frac{1}{s+1}\right]$ I got the idea from here
$(1)$ yield the same result as Hadjicostas. Is this $(1)$ same as Hadjicostas but just write in a different way?
$$\int_{0}^{1}\int_{0}^{1}[-\ln(xy)]^s\left(\frac{1}{\ln(xy)}+\frac{1}{1-xy}\right)dxdy=\Gamma(s+2)\left[\zeta(s+2)-\frac{1}{s+1}\right] \tag{1}$$
also special case, setting $s=-1$ in $(1)$ yeild a $\gamma=0.5772156...$
 A: Note that $$\int_{0}^{1}\int_{0}^{1}\left(-\log\left(xy\right)\right)^{s}\left(\frac{1}{\log\left(xy\right)}+\frac{1}{1-xy}\right)dxdy$$ $$=-\int_{0}^{1}\int_{0}^{1}\left(-\log\left(xy\right)\right)^{s-1}dxdy+\int_{0}^{1}\int_{0}^{1}\frac{\left(-\log\left(xy\right)\right)^{s}}{1-xy}dxdy.$$ Now in the first integral take $x=e^{-u}$ and $y=e^{-v}$. We get $$-\int_{0}^{1}\int_{0}^{1}\left(-\log\left(xy\right)\right)^{s-1}dxdy=-\int_{0}^{\infty}\int_{0}^{\infty}\left(u+v\right)^{s-1}e^{-u}e^{-v}dudv
 $$ and now taking $u+v=z$ we have $$=-\int_{0}^{\infty}\int_{v}^{\infty}z^{s-1}e^{-z}dzdv=-\int_{0}^{\infty}\Gamma\left(s,v\right)dv
 $$ which is, integrating by parts, $$=-\int_{0}^{\infty}v^{s}e^{-v}dv=-\Gamma\left(s+1\right)=-\frac{\Gamma\left(s+2\right)}{s+1}.
 $$ For the second integral we can use a similar argument. We can take again $x=e^{-u}$ and $y=e^{-v}$ and get $$\int_{0}^{1}\int_{0}^{1}\frac{\left(-\log\left(xy\right)\right)^{s}}{1-xy}dxdy=\int_{0}^{\infty}\int_{0}^{\infty}\frac{\left(u+v\right)^{s}}{e^{u+v}-1}dudv$$ $$=\sum_{k\geq1}\int_{0}^{\infty}\int_{0}^{\infty}\left(u+v\right)^{s}e^{-k\left(u+v\right)}dudv$$ and taking $z=\frac{u+v}{k}
 $ we have $$=\sum_{k\geq1}\frac{1}{k^{s+1}}\int_{0}^{\infty}\int_{kv}^{\infty}z^{s}e^{-z}dzdv=\sum_{k\geq1}\frac{1}{k^{s+1}}\int_{0}^{\infty}\Gamma\left(s+1,kv\right)dv
 $$ $$=\sum_{k\geq1}\int_{0}^{\infty}v^{s+1}e^{-kv}dv=\Gamma\left(s+2\right)\zeta\left(s+2\right)
 $$ as wanted.
