Integral $\int_{0}^{1} \int_{0}^{1} \frac{1}{(1+x y) \ln (x y)} d x d y$ 
Evaluate $$\int_{0}^{1} \int_{0}^{1} \frac{1}{(1+x y) \ln (x y)} d x d y$$

I couldn't get very far on this one, so I would appreciate some help =)
My attempt so far (transcribed from the comments):
By extending it to the Dirichlet Eta Function I evaluated this integral to be $\ln 2$. I arrived at the identity below after differentiating a unit square integral expression for $\eta(2)$.
$$\eta(s)\Gamma(s)=\int_0^1\int_0^1 \frac{(-\ln(xy))^{s-2}}{1+xy}dxdy$$
But I can't for the love of it solve it any differently than that. I would love to find an "elementary approach" if possible.
 A: I think it is possible to avoid the Dirichlet eta function.
Put\begin{equation*}
 I=\int_{0}^{1}\int_{0}^{1}\dfrac{1}{(1+xy)\ln(xy)}\,dxdy.
\end{equation*}
By symmetry 
\begin{equation*}
 I = 2\int_{0}^{1}\left(\int_{0}^{x}\dfrac{1}{(1+xy)\ln(xy)}\,dy\right)\, dx.
\end{equation*}
Via the substitution $y=\dfrac{z}{x}$ we get
\begin{equation*}
 I = \int_{0}^{1}\left(\int_{0}^{x^2}\dfrac{1}{x(1+z)\ln(z)}\,dz\right)\,dx.
\end{equation*}
We have a double integral integrated over the domain $0<z<x^2, \, 0<x<1$. If we change the order of integration then we first have to integrate with respect to $x$ over $\sqrt{z}<x<1$ and then with respect to $z$ over $0<z<1.$ Thus
\begin{gather*}
 I = 2\int_{0}^{1}\left(\int_{\sqrt{z}}^{1}\dfrac{1}{x(1+z)\ln(z)}\, dx\right)\, dz = 2\int_{0}^{1}\dfrac{1}{(1+z)\ln(z)}\left[\ln(x)\right]_{\sqrt{z}}^{1}\, dz =\\[2ex]
 \int_{0}^{1}\dfrac{-1}{1+z}\, dz = -\ln(2).
\end{gather*}
A: $$I(a)=\int_0^1\int_0^1 \frac{(xy)^{a}}{(1+xy)\ln(xy)}dxdy\Rightarrow I'(a)=\int_0^1\int_0^1\frac{(xy)^a}{1+xy}dxdy$$
$$\overset{xy=t}=\int_0^1\frac1y\int_0^y\frac{t^a}{1+t}dtdy=\sum_{n=0}^\infty (-1)^n \int_0^1 \frac1{y}\int_0^y t^{a+n}dtdy$$
$$\require{cancel}=\sum_{n=0}^\infty (-1)^n \int_0^1 \cancel{\frac1{y}}\frac{y^{a+n+\cancel1}}{a+n+1}dy=\sum_{n=0}^\infty (-1)^n \frac{1}{(a+n+1)^2}=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(a+n)^2}$$
We are looking to find $I(0)$, but we also have that $\lim\limits_{a\to \infty}I(a)=0$ (see also the comments bellow why not $a\to -\infty$), thus by Newton-Leibniz formula:
$$I(0)=-(I(\infty)-I(0))=-\int^{\infty}_0 I'(a)da=-\sum_{n=1}^\infty(-1)^{n-1} \int^{\infty}_0 \frac{1}{(a+n)^2}da$$
$$=-\sum_{n=1}^\infty(-1)^{n-1}\left(-\frac{1}{a+n}\right)\bigg|^{\infty}_0
=-\sum_{n=1}^\infty\frac{(-1)^{n-1}}{n}=-\ln 2 $$
