Evaluation of a dilogarithmic integral 
Problem. Prove that the following dilogarithmic integral has the indicated value:
  $$\int_{0}^{1}\mathrm{d}x \frac{\ln^2{(x)}\operatorname{Li}_2{(x)}}{1-x}\stackrel{?}{=}-11\zeta{(5)}+6\zeta{(3)}\zeta{(2)}.$$


My attempt:
I began by using the polylogarithmic expansion in terms of generalized harmonic numbers,
$$\frac{\operatorname{Li}_r{(x)}}{1-x}=\sum_{n=1}^{\infty}H_{n,r}\,x^n;~~r=2.$$
Then I switched the order of summation and integration and used the substitution $u=-\ln{x}$ to evaluate the integral:
$$\begin{align}
\int_{0}^{1}\mathrm{d}x \frac{\ln^2{(x)}\operatorname{Li}_2{(x)}}{1-x}
&=\int_{0}^{1}\mathrm{d}x\ln^2{(x)}\sum_{n=1}^{\infty}H_{n,2}x^n\\
&=\sum_{n=1}^{\infty}H_{n,2}\int_{0}^{1}\mathrm{d}x\,x^n\ln^2{(x)}\\
&=\sum_{n=1}^{\infty}H_{n,2}\int_{0}^{\infty}\mathrm{d}u\,u^2e^{-(n+1)u}\\
&=\sum_{n=1}^{\infty}H_{n,2}\frac{2}{(n+1)^3}\\
&=2\sum_{n=1}^{\infty}\frac{H_{n,2}}{(n+1)^3}.
\end{align}$$
So I've reduced the integral to an Euler sum, but unfortunately I've never quite got the knack for evaluating Euler sums. How to proceed from here?
 A: It is easy to see that
$$2\sum^\infty_{n=1}\frac{H_n^{(2)}}{(n+1)^3}=2\sum^\infty_{n=1}\frac{H_{n+1}^{(2)}}{(n+1)^3}-2\sum^\infty_{n=1}\frac{1}{(n+1)^5}=2\sum^\infty_{n=1}\frac{H_n^{(2)}}{n^3}-2\zeta(5)$$
Consider $\displaystyle f(z)=\frac{\pi\cot{\pi z} \ \Psi^{(1)}(-z)}{z^3}$. We know that
$$\pi\cot{\pi z}=\frac{1}{z-n}-2\zeta(2)(z-n)+O((z-n)^3)$$
and
$$\Psi^{(1)}(-z)=\frac{1}{(z-n)^2}+\left(H_n^{(2)}+\zeta(2)\right)+O(z-n)$$
At the positive integers,
\begin{align}
{\rm Res}(f,n)
&=\operatorname*{Res}_{z=n}\left[\frac{1}{z^3(z-n)^3}+\frac{H_n^{(2)}-\zeta(2)}{z^3(z-n)}\right]\\
&=\frac{H_n^{(2)}}{n^3}-\frac{\zeta(2)}{n^3}+\frac{6}{n^5}\\
\end{align}
At the negative integers,
\begin{align}
{\rm Res}(f,-n)&=-\frac{\Psi^{(1)}(n)}{n^3}\\&=\frac{H_{n-1}^{(2)}-\zeta(2)}{n^3}\\&=\frac{H_{n}^{(2)}}{n^3}-\frac{\zeta(2)}{n^3}-\frac{1}{n^5}\tag1
\end{align}
At $z=0$, 
\begin{align}
{\rm Res}(f,0)&=[z^2]\left(\frac{1}{z}-2\zeta(2)z\right)\left(\frac{1}{z^2}+\zeta(2)+2\zeta(3)z+3\zeta(4)z^2+4\zeta(5)z^3\right)\\
&=4\zeta(5)-4\zeta(2)\zeta(3)
\end{align}
Since the sum of the residues $=0$, we conclude that
\begin{align}
\color\red{\int^1_0\frac{\log^2{x} \ {\rm Li}_2(x)}{1-x}{\rm d}x}
&=2\sum^\infty_{n=1}\frac{H_n^{(2)}}{n^3}-2\zeta(5)\\
&=\zeta(2)\zeta(3)-6\zeta(5)+\zeta(2)\zeta(3)+\zeta(5)-4\zeta(5)+4\zeta(2)\zeta(3)-2\zeta(5)\\
&\large{\color\red{=6\zeta(2)\zeta(3)-11\zeta(5)}}
\end{align}
Explanation
$(1):$ Use the functional equation $\displaystyle \Psi^{(1)}(z+1)=-\frac{1}{z^2}+\Psi^{(1)}(z)$ which is derived by differentiating the functional equation of the digamma function, as well as the fact that $\displaystyle H_n^{(2)}=\frac{1}{n^2}+H_{n-1}^{(2)}$. 
As for how to obtain the laurent series, the series for $\Psi(z)$ was cleverly derived here by Random Variable. In essence,
$$\color{blue}{\gamma+\Psi(-z)=\frac{1}{z-n}+H_n+\sum^\infty_{k=1}(-1)^k\left(H_n^{(k+1)}+(-1)^{k+1}\zeta(k+1)\right)(z-n)^k}$$
Differentiating yields
$$\color{blue}{\Psi^{(1)}(-z)=\frac{1}{(z-n)^2}+\sum^\infty_{k=1}(-1)^{k+1}k\left(H_n^{(k+1)}+(-1)^{k+1}\zeta(k+1)\right)(z-n)^{k-1}}$$
For $\pi\cot{\pi z}$,
\begin{align}
\color{blue}{\pi\cot{\pi z}}
&=\Psi(1-z)-\Psi(z) \ \ \ \ \ \text{(reflection formula for digamma function)}\\
&=\int^1_0\frac{t^{z-1}-t^{-z}}{1-t}{\rm d}t \ \ \ \ \ \text{(recall that $\Psi(z)=-\gamma+H_{z-1}$)}\\
&=\sum_{k=0}^\infty\int^1_0\left(t^{z+k-1}-t^{-z+k}\right){\rm d}t\\
&=\sum_{k=0}^\infty\left(\frac{1}{z+k}+\frac{1}{z-k-1}\right)\\
&=\frac{1}{z}+\frac{1}{z-1}+\frac{1}{z+1}+\frac{1}{z-2}+\frac{1}{z+2}+\cdots\\
&=\frac{1}{z}+\sum^\infty_{k=1}\left(\frac{1}{z-k}+\frac{1}{z+k}\right)\\
&=\frac{1}{z}+\sum^\infty_{k=1}\frac{2z}{z^2-k^2}\\
&=\frac{1}{z}-2\sum^\infty_{k=1}\sum^\infty_{m=1}\frac{z^{2m-1}}{k^{2m}}\\
&=\color{blue}{\frac{1}{z}-2\sum^\infty_{m=1}\zeta(2m)z^{2m-1}}\\
&=\pi\cot(\pi (z-n)) \ \ \ \ \ \text{(since cotangent has a period of $\pi$)}\\
&=\color{blue}{\frac{1}{z-n}-2\sum^\infty_{m=1}\zeta(2m)(z-n)^{2m-1}}\\
\end{align}
A: This integral belongs to a wider class of integrals that can always be reduced to single zeta values and to bi-variate zeta values which in turn -- provided the weight of the zeta function in question is not too big-- can also be reduced to single zeta values. My standard way of doing this integrals is as follows. We need to calculate:
\begin{equation}
{\mathcal I}_{0,2}^{(2)} := \int\limits_0^1 \frac{[\log(1/x)]^2}{2!} \cdot \frac{Li_0(x) Li_2(x)}{x} dx
\end{equation}
which is just one half of your integral. Now we use the following identity:
\begin{equation}
\frac{[\log(1/x)]^2}{2!} = \int\limits_{x < \xi_1  < \xi_2 < 1} \prod\limits_{j=1}^2 \frac{d \xi_j}{\xi_j}
\end{equation}
Inserting this into the above and changing order of integration gives:
\begin{eqnarray}
{\mathcal I}_{0,2}^{(2)} = \int\limits_{0 < \xi_1 < \xi_2 < 1} \frac{1}{\xi_1} \frac{1}{\xi_2} \underbrace{\int\limits_0^{\xi_1} \frac{Li_0(x) Li_2(x)}{x} d x}_{\left[Li_1(\xi_1) Li_2(\xi_2) - \int\limits_0^{\xi_1} \frac{[Li_1(x)]^2}{x} dx\right]} \cdot d\xi_1 d\xi_2
\end{eqnarray}
where we integrated by $x$ using integration by parts. Since two minus zero is even we are left with an irreducible integral that we leave unevaluated for the time being. Now we have:
\begin{eqnarray}
{\mathcal I}_{0,2}^{(2)} &=& \int\limits_0^1 \frac{[\log(1/\xi_1)]^1}{1!} \left( \frac{1}{2} [Li_2(\xi_1)]^2 \right)^{'} d\xi_1 - \int\limits_0^1 \frac{[\log(1/x)]^2}{2!} \cdot \frac{[Li_1(x)]^2}{x} dx \\
&=& \frac{1}{2} \left(
\underbrace{\int\limits_0^1 \frac{[Li_2(\xi)]^2}{\xi} d\xi}_{J_1} - 
\underbrace{\int\limits_0^1 [\log(1/\xi)]^2 \cdot \frac{[Li_1(\xi)]^2}{\xi} d\xi}_{J_2}
\right)
\end{eqnarray}
Now from Compute an integral containing a product of powers of logarithms. we have that :
\begin{eqnarray}
J_2&=& -\frac{1}{3} \Psi^{(4)}(1) + 2 \Psi^{(2)}(1) \Psi^{(1)}(1)\\
&=& 8 \zeta(5) - 4 \zeta(3) \zeta(2)
\end{eqnarray}
where $\Psi^{(j)}(1)$ is the polygamma function at unity and $\Psi^{(j)}(1)=(-1)^{j+1} j! \zeta(j+1)$. On the other hand we have:
\begin{eqnarray}
J_1 &=& \sum\limits_{m\ge1,n\ge 1} \frac{1}{m^2} \frac{1}{n^2} \frac{1}{(m+n)}\\
& =& \sum\limits_{m\ge 1} \left(\frac{\zeta(2)}{m^3} - \frac{H_m}{m^4}\right) \\
&=& \zeta(2) \zeta(3) - {\bf H}^{(1)}_4(+1)\\
&=& 2 \zeta(2) \zeta(3) - 3 \zeta(5)
\end{eqnarray}
where in the last line above we used my answer to Calculating alternating Euler sums of odd powers . Therefore:
\begin{equation}
{\mathcal I}^{(2)}_{0,2}= \frac{1}{2} \left(J_1-J_2\right)= \frac{1}{2} \left(-11 \zeta(5)+6 \zeta(3) \zeta(2) \right)
\end{equation}
as expected. Note that exactly the same steps can be performed is we replace the power of the logarithm and the orders of the poly-logarithms by any nonnegative integers. The generic result is then given in An integral involving product of poly-logarithms and a power of a logarithm. .
A: From $$\sum_{n=1}^\infty H_n^{(2)}x^n=\frac{\operatorname{Li}_2(x)}{1-x}$$
it follows that
$$I=\int_0^1\frac{\ln^2(x)\operatorname{Li}_2(x)}{1-x}dx=\sum_{n=1}^\infty H_n^{(2)}\int_0^1 x^n \ln^2(x)dx\\=2\sum_{n=1}^\infty\frac{H_n^{(2)}}{(n+1)^3}=2\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}-2\zeta(5)\tag1$$
By Cauchy product we have
$$\operatorname{Li}_2(x)\operatorname{Li}_3(x)=\sum_{n=1}^\infty\left(\frac{6H_n}{n^4}+\frac{3H_n^{(2)}}{n^3}+\frac{H_n^{(3)}}{n^2}-\frac{10}{n^5}\right)x^n$$ 
set $x=1$ to get
$$\zeta(2)\zeta(3)=6\sum_{n=1}^\infty\frac{H_n}{n^4}+3\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}+\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}-10\zeta(5)\tag{2}$$
Now lets use the well-known identity 
$$\sum_{n=1}^\infty\frac{H_n^{(p)}}{n^q}+\sum_{n=1}^\infty\frac{H_n^{(q)}}{n^p}=\zeta(p)\zeta(q)+\zeta(p+q)$$
set $p=2$ and $q=3$ 
$$\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}=\zeta(2)\zeta(3)+\zeta(5)-\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}$$
Plugging this result in $(2)$ yields
$$\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}=3\zeta(2)\zeta(3)-\frac92\zeta(5)$$
and fianlly, plugging this result in $(1)$, the closed form of $I$ follows. Note that the value $\sum_{n=1}^\infty 
\frac{H_n}{n^4}=3\zeta(5)-\zeta(2)\zeta(3)$ was used in our calculations which can be found using Euler identity.
