Integral $\iint \limits_{{x,y \ \in \ [0,1]}} \frac{\log(1-x)\log(1-y)}{1-xy}dx\,dy=\frac{17\pi^4}{360}$ Hi I am trying to integrate $$
\mathcal{I}:=\iint \limits_{{x,y \ \in \ [0,1]}} \frac{\log(1-x)\log(1-y)}{1-xy}dx\,dy=\int_0^1\int_0^1 \frac{\log(1-x)\log(1-y)}{1-xy}dx \,dy
$$
A closed form does exist. I tried to write
\begin{align}
\mathcal{I} &=\int_0^1 \log(1-y)dy\int_0^1 \log(1-x)\frac{dx}{1-xy} \\
&= \int_0^1\log(1-y)dy \ \int_0^1 \sum_{n\geq0}(xy)^n\, \ln(1-x) \ dx \\
&= \sum_{n\geq 0}\frac{1}{n+1}\int_0^1 \log(1-y) y^n\, dy \\
&= \sum_{n\geq 0}\frac{1}{n+1}\int_0^1 \log(1-y)y^n\, dy = ?
\end{align}
I was able to realize that 
$$
\mathcal{I}=\sum_{n\geq 1}\left(\frac{H_n}{n}\right)^2=\frac{17\zeta_4}{4}=\frac{17\pi^4}{360},\qquad \zeta_4=\sum_{n\geq 1} n^{-4} 
$$
however this does not help me solve the problem.  How can we calculate $\mathcal{I}$? Thanks.
 A: First notice that
$$ \int_{0}^{1} x^{n} \log(1-x)  \ dx = -\int_{0}^{1} x^{n} \sum_{k=1}^{\infty} \frac{x^{k}}{k}  \ dx $$
$$ = -\sum_{k=1}^{\infty} \frac{1}{k} \int_{0}^{1} x^{n+k} \ dx = -\sum_{k=1}^{\infty} \frac{1}{k(n+k+1)}$$
$$ = - \frac{1}{n+1} \sum_{k=1}^{\infty} \left(\frac{1}{k} - \frac{1}{n+k+1} \right)$$
$$ = -\frac{H_{n+1}}{n+1}$$
Then
$$ \int_{0}^{1} \int_{0}^{1} \frac{\log(1-x) \log(1-y)}{1-xy} \ dx \ dy $$
$$ =\sum_{n=0}^{\infty} \int_{0}^{1} x^{n} \log(1-x) \ dx \int_{0}^{1} y^{n} \log(1-y)  \ dy $$
$$ = \sum_{n=0}^{\infty} \left( \frac{H_{n+1}}{n+1} \right)^{2}$$
A closed form for the sum $\sum_{n=1}^{\infty}\left(\frac{H_n}{n}\right)^2$
A: Consider the double integral 
\begin{align}
I = \int_{0}^{1} \int_{0}^{1} \frac{\ln(1-x) \ \ln(1-y)}{1-xy} \ dx dy
\end{align}
which, upon expansion into series form, becomes
\begin{align}
I &= \sum_{n=0}^{\infty} \ \int_{0}^{1} x^{n} \ln(1-x) dx \ \int_{0}^{1} y^{n} \ln(1-y) dy \\
&= \sum_{n=0}^{\infty} \left( \int_{0}^{1} x^{n} \ \ln(1-x) \ dx \right)^{2}.
\end{align}
Without much difficulty it can be shown that
\begin{align}
\int_{0}^{1} x^{n} \ \ln(1-x) \ dx &= \left. \partial_{\mu} B(n+1, \mu+1) \right|_{\mu = 1} \\
&= B(n+1,1) \left( \psi(1) - \psi(n+2) \right) \\
&= \frac{H_{n+1}}{n+1}. 
\end{align}
Now, returning to the primary integral, it is seen that
\begin{align}
I = \sum_{n=0}^{\infty} \left( \frac{H_{n+1}}{n+1}\right)^{2} = \sum_{n=1}^{\infty} \left( \frac{H_{n}}{n} \right)^{2} = \frac{17 \zeta(4)}{4}.
\end{align}
Hence
\begin{align}
\int_{0}^{1} \int_{0}^{1} \frac{\ln(1-x) \ \ln(1-y)}{1-xy} \ dx dy = \frac{17 \zeta(4)}{4}.
\end{align}
A: Wow, RV and L are pretty fast. Well, even though their ways are more efficient, I worked out yet another way so I may as well post it. 
$$\int_{0}^{1}\log(1-y)\int_{0}^{1}\frac{\log(1-x)}{1-xy}dxdy$$
Let $u=xy, \;\ \frac{du}{x}=dy$
$$\int_{0}^{x}\frac{\log(1-u/x)}{1-u}du\int_{0}^{1}\frac{\log(1-x)}{x}dx$$
The left integral w.r.t u is a classic dilog: $\int_{0}^{x}\frac{\log(1-u/x)}{1-u}du=-Li_{2}(\frac{x}{1-x})$
Thus,  $$\int_{0}^{1}\frac{Li_{2}(\frac{x}{1-x})\log(1-x)}{x}dx$$
But,  $$-Li_{2}(\frac{x}{1-x})=\sum_{n=1}^{\infty}\frac{H_{n}}{n}x^{n}$$
giving:
$$\int_{0}^{1}\log(1-x)\sum_{n=1}^{\infty}\frac{H_{n}}{n}x^{n-1}dx$$
Also, $\int_{0}^{1}x^{n-1}\log(1-x)dx=\frac{H_{n}}{n}$
So, we finally obtain:
$$\sum_{n=1}^{\infty}\frac{H_{n}^{2}}{n^{2}}=\frac{17\pi^{4}}{360}$$
