Evaluating the integral $\int_1^{\infty} \int_1^{\infty} \frac{\Gamma(x+1)\Gamma(y+1)}{x y \Gamma(x+y+2)} \ dx \ dy$

To compute $$\int_{0}^{1}\operatorname{li}(x)\operatorname{li}(1-x)\,dx$$ i'd rather exploit the fact that the logarithmic integral function admits a representation in terms of the shifted Legendre polynomials, since: $$\begin{eqnarray*}\int_{0}^{1}\operatorname{li}(x)\,L_n(2x-1)\,dx &=& \frac{1}{n!}\int_{0}^{1}\operatorname{li}(x)\frac{d^n}{dx^n}(x^2-x)^n\,dx\\&=&\frac{1}{n!}\left.\frac{L_{n+1}(2x-1)-L_{n-1}(2x-1)}{4n+2}\operatorname{li}(x)\right|_{0}^{1}\\&-&\frac{1}{(4n+2)n!}\int_{0}^{1}\frac{L_{n+1}(2x-1)-L_{n-1}(2x-1)}{\log x}\,dx.\end{eqnarray*}$$ So, if we set $L_n(2x-1)=Q_n(x)$ for shortness, we have: $$\operatorname{li}(x) = -\log 2-\sum_{n=1}^{+\infty} c_n Q_n(x), $$ where: $$ c_n = \frac{1}{2n!}\int_{0}^{1}\frac{Q_{n+1}(x)-Q_{n-1}(x)}{\log x}\,dx \tag{1}$$ is a linear combination of logarithms of integers by Frullani's theorem and: $$\operatorname{li}(1-x) = -\log 2-\sum_{n=1}^{+\infty}(-1)^n c_n Q_n(x), $$ giving: $$\int_{0}^{1}\operatorname{li}(x)\operatorname{li}(1-x)\,dx=\log^2 2+\sum_{n=1}^{+\infty}\frac{(-1)^n c_n^2}{2n+1}.\tag{2}$$