Double Integral $\int\limits_0^1\!\!\int\limits_0^1\frac{(xy)^s}{\sqrt{-\log(xy)}}\,dx\,dy$
Solution 1:
Here is an approach.
We have, for $s>0$, $\sigma>0$,
$$ \int_0^1\!\int_0^1\frac{x^{s-1}y^{\sigma-1}}{\sqrt{-\log (xy)}}\,dx\,dy=\frac{\sqrt{\pi}}{s\sqrt{\sigma}+\sigma\sqrt{s}}.\tag1 $$
Proof. One may observe that $$ \frac1{\sqrt{a}}=\frac2{\sqrt{\pi}}\int_0^\infty e^{\large-at^2}dt, \qquad a>0. \tag2 $$ Assume $s>0$, $\sigma>0$ and $\sigma\neq s$. Then we obtain $$ \begin{align} \iint_{[0,1]\times[0,1]} \frac{x^{s-1}y^{\sigma-1}}{\sqrt{-\log (xy)}}\,dx\,dy&=\iint_{[0,1]\times[0,1]} x^{s-1}y^{\sigma-1}\left(\frac1{\sqrt{-\log (xy)}}\right)dx\,dy\\\\ &=\iint_{[0,1]\times[0,1]} x^{s-1}y^{\sigma-1}\left(\frac2{\sqrt{\pi}}\int_0^\infty e^{\large t^2\log(xy)}dt\right)dx\,dy\\\\ &=\frac2{\sqrt{\pi}}\int_0^\infty\iint_{[0,1]\times[0,1]} x^{s-1}y^{\sigma-1}\times x^{t^2}y^{t^2}dx\,dy\,dt\\\\ &=\frac2{\sqrt{\pi}}\int_0^\infty\left(\int_0^1x^{t^2+s-1}dx\right)\times \left(\int_0^1y^{t^2+\sigma-1}dy\right)\,dt\\\\ &=\frac2{\sqrt{\pi}}\int_0^\infty\frac1{(t^2+s)(t^2+\sigma)}\,dt\\\\ &=\frac2{\sqrt{\pi}(\sigma-s)}\left(\int_0^\infty\frac{1}{\left(t^2+s\right)}\,dt-\int_0^\infty\frac{1}{\left(t^2+\sigma\right)}\,dt\right)\\\\ &=\frac2{\sqrt{\pi}(\sigma-s)}\left(\frac{\pi}{2\sqrt{s}}-\frac{\pi}{2\sqrt{\sigma}}\right)\\\\ &=\frac{\sqrt{\pi}}{s\sqrt{\sigma}+\sigma\sqrt{s}}. \end{align} $$ We obtain the desired integral by setting $s:=s+1$ in $(1)$ and $\sigma :=s+1$ giving
$$ \int_0^1\!\int_0^1\frac{(xy)^s}{\sqrt{-\log (xy)}}\,dx\,dy=\frac{\sqrt{\pi}}{2(1+s)^{3/2}}.\tag3 $$
Solution 2:
\begin{align} u & = xy \\[10pt] v & = \frac y x \\[10pt] du\,dv = \frac{\partial(u,v)}{\partial(x,y)} \, dx \, dy & = \frac {2y} x \, dx\,dy = 2v\,dx\,dy \\[10pt] \frac{du\,dv}{2v} & = dx\,dy \end{align}
\begin{align} & \int_0^1 \int_u^{1/u} \frac{u^s}{\sqrt{-\log u}} \, \frac{dv \, du}{2v} \\[10pt] = {} & \int_0^1 \frac{u^s}{\sqrt{-\log u}} \cdot \frac 1 2 \left( \log \frac 1 u - \log u \right) \, du \\[10pt] = {} & \int _0^1 u^s \sqrt{-\log u} \, du \tag 1 \end{align}
\begin{align} w & = \sqrt{-\log u} \\[10pt] e^{-w^2} & = u \\[10pt] -2w e^{-w^2}\,dw & = du \end{align}
So $(1)$ becomes $\displaystyle \int_\infty^0 e^{-sw^2} w ( -2w e^{-w^2} \, dw) = 2\int_0^\infty w^2 e^{-(s+1)w^2} \, dw$.
We can reduce this to the expected value of a function of a random variable with a standard normal distribution:
\begin{align} (s+1) w^2 & = \frac{z^2} 2 \\[10pt] dw & = \frac{dz}{\sqrt2\sqrt{s+1}} \end{align}
Our integral becomes $$ 2\int_0^\infty \frac{z^2}{2(s+1)} e^{-z^2/2} \frac{dz}{\sqrt2\sqrt{s+1}} = \frac 1 {(s + 1)^{3/2}\sqrt 2} \int_0^\infty z^2 e^{-z^2/2}\,dz. \tag 1 $$ The last integral is $$ \int_0^\infty z^2 \left( \frac{e^{-z^2/2}}{\sqrt{2\pi}} \right) \, dz \cdot \sqrt{2\pi}. $$
This integral is half the expected value of a $\chi^2_1$ random variable; thus it is equal to $1/2$, so we get $(1/2)\sqrt{2\pi},$ so the whole expression $(1)$ is $$ \frac 1 {(s + 1)^{3/2}\sqrt 2} \cdot \frac 1 2 \cdot \sqrt{2\pi} = \frac{\sqrt \pi}{2(s + 1)^{3/2}}. $$
Solution 3:
Our aim is to show that
$$ \int_0^1\int_0^1\frac{(xy)^s}{\sqrt{-\log xy}}\,dx\,dy =2\int_0^1\int_y^1\frac{(xy)^s}{\sqrt{-\log xy}}\,dx\,dy =\frac{\sqrt{\pi}}{2(1+s)^{3/2}}. $$
We will need the error function, $$ \def\erf{\,\text{erf}\,} \erf(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt. $$ Let us start to calculate a primitive with respect to $x$. With $$ u=\sqrt{1+s}\sqrt{-\log(xy)} $$ we find that $$ \int\frac{(xy)^s}{\sqrt{-\log xy}}\,dx=-\frac{\sqrt{\pi}}{y\sqrt{1+s}}\erf\bigl(\sqrt{1+s}\sqrt{-\log(xy)}\bigr). $$ Inserting limits, we find that $$ \int_y^1\frac{(xy)^s}{\sqrt{-\log xy}}\,dx =-\frac{\sqrt{\pi}}{y\sqrt{1+s}}\erf\bigl(\sqrt{1+s}\sqrt{-\log(y)}\bigr) +\frac{\sqrt{\pi}}{y\sqrt{1+s}}\erf\bigl(\sqrt{1+s}\sqrt{-2\log(y)}\bigr) $$ Next, we want to integrate with respect to $y$. The integrals are very similar. We find that ($u=\sqrt{-a\log y}$ and integrating by parts twice) $$ \begin{split} \int \frac{1}{y}\erf\bigl(\sqrt{1+s}\sqrt{-a\log y}\bigr)\,dy &=-\frac{2}{a}\int u\erf(\sqrt{1+s}u)\,du\\ &=-\frac{1}{a}u^2\erf(\sqrt{1+s}u)+\frac{\sqrt{1+s}}{a\sqrt{\pi}}\int u^2 e^{-(1+s)u^2}\,du\\ &=-\frac{1}{a}u^2\erf(\sqrt{1+s}u) -\frac{1}{2a\sqrt{\pi}\sqrt{1+s}}ue^{-(1+s)u^2}\\ &\qquad+\frac{1}{4a(1+s)}\erf(\sqrt{1+s}u). \end{split} $$ I leave it to you to insert $a=1$ and $a=2$, and to work with the limits.
Solution 4:
$\sqrt{-\log(\cdot)}$ strikes me as the inverse function of a Gaussian function.
So it may be worth to set $x=e^{-u^2}$ and $y=e^{-v^2}$ to get:
$$ I =\iint_{(0,1)^2}\frac{(xy)^s}{\sqrt{-\log(xy)}}\,dx\,dy = \iint_{(0,+\infty)^2}\frac{4uv\,e^{-(s+1)(u^2+v^2)}}{\sqrt{u^2+v^2}}\,du\,dv $$ and that looks good. If we switch to polar coordinates we have:
$$ I = \left(\int_{0}^{\pi/2}4\sin\theta\cos\theta\,d\theta\right)\cdot\left(\int_{0}^{+\infty}\rho^2 e^{-(s+1)\rho^2}d\rho\right)$$ hence:
$$ I = \frac{\sqrt{\pi}}{2}\cdot (s+1)^{-3/2}.$$