How to evaluate this strange double integral (Vardi's integral type)? I found this strange double integral from some a site (Can't remember where).
I would like to know how this integral it is evalauted:
$$\int_{0}^{1}\int_{0}^{1}\ln^2(x) \frac{\ln\ln\left(\frac{1}{x}\right)-\ln\ln\left(\frac{1}{y}\right)}{\ln^2(x)+\ln^2(y)}\mathrm dx \mathrm dy=0.61685...\tag1$$
It looks like the value $0.61685...$ apporaches $\left(\frac{\pi}{4}\right)^2$
How would one go about to evaluate this above integral?
Intgeral $(1)$ looked similar to Vardi's integral,
$$\int_{0}^{1}\ln\ln\left(\frac{1}{x}\right)\frac{\mathrm dx}{1-x+x^2}=\frac{\pi}{\sqrt{3}}\ln\left[\frac{2\pi^{5/6}}{\Gamma(1/6)}\right]$$
 A: Let
$$ I = \int \limits_0^1 \int \limits_0^1 \ln^2(x) \frac{\ln(-\ln(x)) - \ln(-\ln(y))}{\ln^2(x)+\ln^2(y)} \, \mathrm{d} x \, \mathrm{d} y \, .$$
We can switch the names of $x$ and $y$ and average the result and the original version to get
$$I = \frac{1}{2} \int \limits_0^1 \int \limits_0^1 \left[\ln(-\ln(x)) - \ln(-\ln(y))\right] \frac{\ln^2(x) - \ln^2(y)}{\ln^2(x)+\ln^2(y)} \, \mathrm{d} x \, \mathrm{d} y \, .$$
Now we let $x = \mathrm{e}^{-s}$ and $y = \mathrm{e}^{-t}$ to obtain
$$I = \frac{1}{2} \int \limits_0^\infty \int \limits_0^\infty \ln\left(\frac{s}{t}\right) \frac{s^2 - t^2}{s^2+t^2} \,\mathrm{e}^{-(s+t)} \,\mathrm{d} s \, \mathrm{d} t = \int \limits_0^\infty \int \limits_0^t \ln\left(\frac{s}{t}\right) \frac{s^2 - t^2}{s^2+t^2} \,\mathrm{e}^{-(s+t)} \,\mathrm{d} s \, \mathrm{d} t \, ,$$
where the last step follows from the symmetry of the integrand. Introducing $u = \frac{s}{t}$ and $v = s+t$ simplifies this integral considerably. We have $u \in (0,1)$ and $v \in (0,\infty)$ and the Jacobian of the inverse transformation $s =  \frac{u v}{1+u}, t= \frac{v}{1+u}$ is $\frac{v}{(1+u)^2}$, so
$$ I = \int \limits_0^\infty \int \limits_0^1 - \ln(u) \frac{1-u^2}{1+u^2} \mathrm{e}^{-v} \frac{v}{(1+u)^2} \, \mathrm{d} u \, \mathrm{d} v  = \int \limits_0^1 - \ln(u) \frac{1-u}{(1+u)(1+u^2)} \, \mathrm{d} u \int \limits_0^\infty v \, \mathrm{e}^{-v} \mathrm{d} v \, .$$
The $v$-integral is equal to $\operatorname{\Gamma}(2) = 1$ and the $u$-integral can be computed using partial fraction decomposition and integration by parts:
\begin{align}
I &= \int \limits_0^1 - \ln(u) \left[\frac{1}{1+u} - \frac{u}{1+u^2}\right] \mathrm{d} u = \int \limits_0^1 \frac{\ln(1+u) - \frac{1}{2} \ln(1+u^2)}{u} \, \mathrm{d} u \\
&= \frac{3}{4} \int \limits_0^1 \frac{\ln(1+u)}{u} \, \mathrm{d} u = \frac{3}{4} \operatorname{\eta}(2) = \frac{3}{4} \frac{\pi^2}{12} = \frac{\pi^2}{16} \, .
\end{align}
The final integral is discussed here. By the way, the same method can be used to calculate a similar integral, which happens to have the same value:
$$ \int \limits_0^1 \int \limits_0^1 - \ln(x) \frac{\ln(-\ln(x)) - \ln(-\ln(y))}{\ln^2(x)+\ln^2(y)} \, \mathrm{d} x \, \mathrm{d} y = \frac{\pi^2}{16} \, . $$
