Triple Euler sum result $\sum_{k\geq 1}\frac{H_k^{(2)}H_k }{k^2}=\zeta(2)\zeta(3)+\zeta(5)$ In the following thread
I arrived at the following result 
$$\sum_{k\geq 1}\frac{H_k^{(2)}H_k }{k^2}=\zeta(2)\zeta(3)+\zeta(5)$$
Defining 
$$H_k^{(p)}=\sum_{n=1}^k \frac{1}{n^p},\,\,\, H_k^{(1)}\equiv H_k $$
But, it was after long evaluations and considering many variations of product of polylogarithm integrals.
I think there is an easier approach to get the solution, any ideas ?
 A: We are going to evaluate our sum by establishing a system of two relations. 
Lets establish the first relation and using the derivative of beta function ( see here) , we have
$$-\int_0^1x^{n-1}\ln^3(1-x)\ dx=\frac{H_n^3}{n}+3\frac{H_nH_n^{(2)}}{n}+2\frac{H_n^{(3)}}{n}$$
divide both sides by $n$ then take the sum with respect to $n$, we get
\begin{align}
R_1&=\sum_{n=1}^\infty\frac{H_n^3}{n^2}+3\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}+2\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}=-\int_0^1\frac{\ln^3(1-x)}{x}\sum_{n=1}^\infty\frac{x^n}{n}\ dx\\
&=\int_0^1\frac{\ln^4(1-x)}{x}\ dx=\int_0^1\frac{\ln^4x}{1-x}\ dx=24\zeta(5)
\end{align}
Then $$\boxed{R_1=\sum_{n=1}^\infty\frac{H_n^3}{n^2}+3\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}+2\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}=24\zeta(5)}$$

To get the second relation, we need to use the sterling number formula ( check here)
$$ \frac{\ln^k(1-x)}{k!}=\sum_{n=k}^\infty(-1)^k \begin{bmatrix} n \\  k \end{bmatrix}\frac{x^n}{n!}$$
letting $k=4$ and using $\displaystyle\begin{bmatrix} n \\  4 \end{bmatrix}=\frac{1}{3!}(n-1)!\left[\left(H_{n-1}\right)^3-3H_{n-1}H_{n-1}^{(2)}+2H_{n-1}^{(3)}\right],$ we get $$\frac14\ln^4(1-x)=\sum_{n=1}^\infty\frac{x^{n+1}}{n+1}\left(H_n^3-3H_nH_n^{(2)}+2H_n^{(3)}\right)$$
differentiate both sides with respect to $x$, we get
$$-\frac{\ln^3(1-x)}{1-x}=\sum_{n=1}^\infty x^n\left(H_n^3-3H_nH_n^{(2)}+2H_n^{(3)}\right)$$
multiply both sides by $\ln x/x$ then integrate with respect to $x$, we get
\begin{align}
R_2&=\sum_{n=1}^\infty\frac{H_n^3}{n^2}-3\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}+2\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}\\
&=\int_0^1\frac{\ln^3(1-x)\ln x}{x(1-x)}\ dx=\int_0^1\frac{\ln^3x\ln(1-x)}{x(1-x)}\ dx\\
&=-\sum_{n=1}^\infty H_n\int_0^1x^{n-1}\ln^3x\ dx=6\sum_{n=1}^\infty\frac{H_n}{n^4}=18\zeta(5)-6\zeta(2)\zeta(3)
\end{align}
Then $$\boxed{R_2=\sum_{n=1}^\infty\frac{H_n^3}{n^2}-3\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}+2\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}=18\zeta(5)-6\zeta(2)\zeta(3)}$$
now we are ready to calculate our sum:
$$R_1-R_2=6\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}=6\zeta(5)+6\zeta(2)\zeta(3)$$
or 

$$\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}=\zeta(5)+\zeta(2)\zeta(3)$$

And as a bonus:
$$R_1+R_2=2\sum_{n=1}^\infty\frac{H_n^3}{n^2}=42\zeta(5)-6\zeta(2)\zeta(3)-4\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}$$
using the definition of $H_n^{(3)}$ and the partial fraction decomposition, its easy to prove $$\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^2}=\frac{11}2\zeta(5)-2\zeta(2)\zeta(3)$$
which follows 

$$\sum_{n=1}^\infty\frac{H_n^3}{n^2}=10\zeta(5)+\zeta(2)\zeta(3)$$

A: Consider the integral $$I= - \int_0^1 \frac{\ln(1-x^2)}{\sqrt{1-x^2}} (\sin^{-1} x)^4 \,dx.$$
Since $$(\sin^{-1} x)^4 = \frac32 \sum_{n=1}^{\infty} \cfrac{2^{2n} H_{n-1}^{(2)}}{n^2 \binom{2n}{n}} \,x^{2 n} \tag{1}$$
and 
$$-\int_0^1 \frac{\ln(1-x^2)}{\sqrt{1-x^2}} x^{2n}\,dx= \frac{\pi}{2} \binom{2n}{n} \frac{(H_n + 2\ln2)}{2^{2n}}, \tag{2}$$
we have 
$$\begin{align*}  
&I= - \frac32 \sum_{n=1}^{\infty} \cfrac{2^{2n} H_{n-1}^{(2)}}{n^2 \binom{2n}{n}} \int_0^1 \frac{\ln(1-x^2)}{\sqrt{1-x^2}} x^{2n}\,dx
\\&= \frac{3 \pi}{4} \sum_{n=1}^{\infty} \frac{H_{n-1}^{(2)}}{n^2} ( H_n +2 \ln2 )
\\& = \frac{\pi^5}{80} \ln2 + \frac{3 \pi}{4} \sum_{n=1}^{\infty} \frac{H_{n-1}^{(2)}\,H_n}{n^2}.
\end{align*}$$
However, substituting $x\mapsto \sin x$ and employing the fourier expansion of $\ln \cos x$:
$$\begin{align*}
& I= -2 \int_0^{\pi/2} x^4 \, \ln\cos x\, dx
\\&= 2 \int_0^{\pi/2} x^4 \left(\ln2 + \sum_{n=1}^{\infty} \frac{(-1)^n \cos(2 x n)}{n} \right)dx
\\&=  \frac{\pi^5}{80}\ln2  + 2 \sum_{n=1}^{\infty} \frac{(-1)^n}{n} \int_0^{\pi/2} x^4 \cos(2 x n) dx
\\&= \frac{\pi^5}{80}\ln2 + 2 \sum_{n=1}^{\infty} \frac{(-1)^n}{n} \frac{(-1)^n}{n^2}\left(\frac{\pi^3}{8}-\frac{3 \pi}{4 n^2}\right)
\\&= \frac{\pi^5}{80}\ln2 + \frac{\pi^3}{4}\zeta(3) - \frac{3 \pi}{2} \zeta(5). \end{align*}$$
Therefore,

$$\sum_{n=1}^{\infty} \frac{H_{n-1}^{(2)}\,H_n}{n^2} = 2\zeta(2)\,\zeta(3)-2\zeta(5).$$

Finish off using Euler's formula for $\sum H_n/n^q $.

Notes.
You may find a proof of $(1)$ here, and $(2)$ is just the derivative of a beta function. The swap of the sum and the integral should be justified.
I found this proof while exploring series involving $H_n^{(2)}$. Using the same method, I also obtained the following
related results:
$$\sum_{n=1}^{\infty} \frac{H_{n-1}^{(2)}\,H_{2n}}{n^2} =\frac{11}{4}\zeta(2)\,\zeta(3)-\frac{47}{16}\zeta(5) \tag{3}$$
$$\sum_{n=1}^{\infty} \frac{H_{n-1}^{(2)}\,H_{n}^2}{n^2} = 4 \zeta(3)^2 - \frac{5}{8} \zeta(6) \tag{4}$$
$$\sum_{n=1}^{\infty} \frac{H_n \left(H_{n-1}^{(2)2}-H_{n-1}^{(4)}\right)}{n^2} = 3\,\zeta(3)\,\zeta(4)-4\,\zeta(2)\,\zeta(5)+4\,\zeta(7) \tag{5}$$
and others.
A: A neat way of dealing with the series by series manipulations (no use of integrals at all) and that circumvents all single harmonic series except the ones of the type $\displaystyle \sum_{n=1}^{\infty}\frac{H_n}{n^m}$ may be found in the book  (Almost) Impossible Integrals, Sums, and Series pages $398-401$ (definitely one of the best solutions in the book that is worth seeing).
A: Lets start with the following double sum
$$\sum_{k=1}^\infty\sum_{n=1}^\infty\frac{H_{n-1}}{n^2(n+k)^2}=\sum_{k=1}^\infty\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}\left(-\int_0^1 x^{n+k-1}\ln xdx\right)$$
$$=-\int_0^1 \ln x\left(\sum_{k=1}^\infty x^{k-1}\right)\left(\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}x^n\right)dx$$
$$=-\int_0^1 \ln x\left(\frac{1}{1-x}\right)\left(\frac12\int_0^x \frac{\ln^2(1-y)}{y}dy\right)dx$$
$$=-\frac12\int_0^1 \frac{\ln^2(1-y)}{y}\left(\int_y^1\frac{\ln x}{1-x}dx\right)dy$$
$$=\frac12\int_0^1 \frac{\ln^2(1-y)\text{Li}_2(1-y)}{y}dy\overset{1-y=x}{=}\frac12\int_0^1 \frac{\ln^2x\text{Li}_2(x)}{1-x}dx$$
$$=\frac12\sum_{n=1}^\infty H_n^{(2)}\int_0^1 x^n \ln^2xdx=\sum_{n=1}^\infty\frac{H_n^{(2)}}{(n+1)^3}=\sum_{n=1}^\infty\frac{H_{n-1}^{(2)}}{n^3}=\sum_{n=1}^\infty\frac{H_{n}^{(2)}}{n^3}-\zeta(5)\tag1$$
On the other hand and by writing $\sum_{k=1}^\infty\frac{1}{(n+k)^2}=\zeta(2)-H_n^{(2)}$ we have
$$\sum_{k=1}^\infty\sum_{n=1}^\infty\frac{H_{n-1}}{n^2(n+k)^2}=\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}\left(\zeta(2)-H_n^{(2)}\right)$$
$$=\zeta(2)\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}-\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}+\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^3}\tag2$$
$(1)-(2)$ gives us
$$\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}=\zeta(5)+\zeta(2)\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}=\zeta(5)+\zeta(2)\zeta(3)$$
where the last sum follows from dividing both sides of $\sum_{n=1}^\infty\frac{H_{n-1}}{n}x^n=\frac12\ln^2(1-x)$ by $x$ then $\int_0^1$ , i.e
$$\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}=\frac12\int_0^1\frac{\ln^2(1-x)}{x}dx=\frac12\int_0^1\frac{\ln^2x}{1-x}dx=\zeta(3)$$
A: I think it is reasonable to start with:
$$\sum_{k=1}^{+\infty}\frac{H_k^{(2)}H_k}{k^2}=\sum_{k=1}^{+\infty}\frac{H_k}{k^4}+\sum_{k=1}^{+\infty}\frac{H_k}{k^2}\sum_{1\leq j< k}\frac{1}{j^2},\tag{1}$$
that leads to:
$$\sum_{k=1}^{+\infty}\frac{H_k^{(2)}H_k}{k^2}=\left(\sum_{k=1}^{+\infty}\frac{H_k}{k^2}\right)\left(\sum_{j=1}^{+\infty}\frac{1}{j^2}\right)-\sum_{k=1}^{+\infty}\frac{1}{k^2}\sum_{1\leq j< k}\frac{H_j}{j^2},\tag{2}$$
Now since:
$$\operatorname{Li}_2(x)+\frac{\log^2(1-x)}{2}=\sum_{k=1}^{+\infty}\frac{H_k}{k}x^k,\tag{3}$$
$$\frac{\log^2(1-x)}{2}=\sum_{k=1}^{+\infty}\frac{H_{k-1}}{k}x^k,\tag{4}$$
follows. By dividing by $x$ and integrating between $0$ and $1$ we get:
$$\sum_{k=1}^{+\infty}\frac{H_{k-1}}{k^2}=\frac{1}{2}\int_{0}^{1}\frac{\log^2(x)}{1-x}dx=\frac{1}{2}\int_{0}^{+\infty}\frac{u^2}{e^u-1}du=\zeta(3),\tag{5}$$
so:
$$\sum_{k=1}^{+\infty}\frac{H_k^{(2)}H_k}{k^2}=2\zeta(2)\zeta(3)-\sum_{k=1}^{+\infty}\frac{1}{k^2}\sum_{1\leq j< k}\frac{H_j}{j^2}.\tag{6}$$
For the last term consider:
$$-\frac{\log(1-xy)}{y(1-xy)}=\sum_{k=1}^{+\infty}H_k x^k y^{k-1}, \tag{7}$$
multiply both terms by $-\log(y)$ and integrate between $0$ and $1$ with respect to $y$:
$$\int_{0}^{1}\frac{\log(y)\log(1-xy)}{y(1-xy)}dy = \sum_{k=1}^{+\infty}\frac{H_k}{k^2}x^k.\tag{8}$$ 
Multiplying both sides by $-\frac{\log x}{1-x}$ and integrating between $0$ and $1$ with respect to $x$ should do the trick. For the last part it is only required to find an appropriate birational diffeomorphism of the unity square that puts the integral in a nicer form - a sort of "reverse Viola-Rhin method".
A: I think this is a shorter solution
By Cauchy product we have 
$$\ln(1-x)\operatorname{Li}_2(x)=-\sum_{n=1}^\infty\left(2\frac{H_n}{n^2}+\frac{H_n^{(2)}}{n}-\frac3{n^3}\right)x^n$$
multiply both sides by $\frac{\ln(1-x)}{x}$ then integrate from $x=0$ to $x=1$ and use the fact that $-\int_0^1x^{n-1}\ln(1-x)\ dx=\frac{H_n}{n}$ we get
$$\sum_{n=1}^\infty\left(2\frac{H_n}{n^2}+\frac{H_n^{(2)}}{n}-\frac3{n^3}\right)\left(\frac{H_n}{n}\right)=\int_0^1\frac{\operatorname{Li}_2(x)\ln^2(1-x)}{x}dx=I$$
or $$\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}=3\sum_{n=1}^\infty\frac{H_n}{n^4}-2\sum_{n=1}^\infty\frac{H_n^2}{n^3}+I\tag1$$
I proved in this solution
$$\int_0^1\frac{\zeta(2)-\operatorname{Li}_2(x)}{x}\ln^2(1-x)\ dx=2\sum_{n=1}^\infty\frac{H_n^2}{n^3}-2\sum_{n=1}^\infty\frac{H_n}{n^4}$$
or 
$$I=2\zeta(2)\zeta(3)-2\sum_{n=1}^\infty\frac{H_n^2}{n^3}+2\sum_{n=1}^\infty\frac{H_n}{n^4}\tag2$$
plugging $(2)$ in $(1)$ we get
$$ \sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}=5\sum_{n=1}^\infty\frac{H_n}{n^4}-4\sum_{n=1}^\infty\frac{H_n^2}{n^3}+2\zeta(2)\zeta(3)$$
Substituting 
$$\sum_{n=1}^\infty\frac{H_n}{n^4}=3\zeta(5)-\zeta(2)\zeta(3)$$
$$\sum_{n=1}^\infty\frac{H_n^2}{n^3}=\frac72\zeta(5)-\zeta(2)\zeta(3)$$
gives
$$\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{n^2}=\zeta(5)+\zeta(2)\zeta(3)$$
where $\sum_{n=1}^\infty\frac{H_n}{n^4}$ can be obtained from using Euler identity and $\sum_{n=1}^\infty\frac{H_n^2}{n^3}$ can be found here.
