An interesting investigation started here and it showed that $$ \sum_{k\geq 1}\left(\zeta(m)-H_{k}^{(m)}\right)^2 $$ has a closed form in terms of values of the Riemann $\zeta$ function for any integer $m\geq 2$.
I was starting to study the cubic analogue $ \sum_{k\geq 1}\left(\zeta(m)-H_{k}^{(m)}\right)^3 $ and I managed to prove through summation by parts that $$ \sum_{k\geq 1}\frac{\left(H_k^{(2)}\right)^2}{k^2} =\frac{1}{3}\zeta(2)^3-\frac{2}{3}\zeta(6)+\zeta(3)^2 $$ where the LHS, according to Flajolet and Salvy's notation, is $S_{22,2}$. An explicit evaluation of $\sum_{k\geq 1}\left(\zeta(2)-H_{k}^{(2)}\right)^3 $ is completed by the computation of $$\boxed{ S_{12,2} = \sum_{k\geq 1}\frac{H_k H_k^{(2)}}{k^2} }$$ which has an odd weight, hence it is not covered by Thm 4.2 of Flajolet and Salvy. On the other hand they suggest that by the kernel $(\psi(-s)+\gamma)^4$ the previous series and the cubic Euler sum $\sum_{n\geq 1}\frac{H_n^3}{(n+1)^2}=\frac{15}{2}\zeta(5)+\zeta(2)\,\zeta(3)$ are strictly related.

Question: can you help me completing this sketch, in order to get an explicit value for $S_{12,2}$ and for $\sum_{k\geq 1}\left(\zeta(2)-H_{k}^{(2)}\right)^3 $? Alternative techniques to summation by parts and residues are equally welcome.

Update: I have just realized this is solved by Mike Spivey's answer to Zaid's question here. On the other hand, Mike Spivey's approach is extremely lengthy, and I would be happy to see a more efficient derivation of $S_{12,2}=\zeta(2)\,\zeta(3)+\zeta(5)$.


I am skipping the derivation of Euler Sums $\displaystyle S(1,1;3) = \sum\limits_{n=1}^{\infty} \frac{H_n^2}{n^3}$ and $\displaystyle S(2;3) = \sum\limits_{n=1}^{\infty} \frac{H_n^{(2)}}{n^3}$ (for a derivation see solutions to $\textbf{problem 5406}$ proposed by C.I. Valean, Romania in SSMA).

Now, consider the partial fraction decomposition, \begin{align*}\sum\limits_{k=1}^{n-1}\left(\frac{1}{k(n-k)}\right)^2 = \frac{2}{n^2}\left(H_n^{(2)} + \frac{2H_n}{n}-\frac{3}{n^2}\right) \qquad \cdots (\star)\end{align*}

Multiplying both sides of $(\star)$ with $H_n$ and summing over $n \ge 1$ (and making the change of varible $m = n+k$):

\begin{align} \sum\limits_{n=1}^{\infty}\frac{H_n}{n^2}\left(H_n^{(2)}+2\frac{H_n}{n} - \frac{3}{n^2}\right) &= \sum\limits_{n=1}^{\infty}\sum\limits_{k=1}^{n-1} \frac{H_n}{k^2(n-k)^2} \\&= \sum\limits_{m,k=1}^{\infty} \frac{H_{m+k}}{m^2k^2} \tag{0}\\&= \sum\limits_{m,k,j=1}^{\infty} \frac{mj+kj}{m^2k^2j^2(m+k+j)} \tag{1} \\&= 2\sum\limits_{m,k,j=1}^{\infty} \frac{jk}{m^2k^2j^2(m+k+j)} \tag{2} \\&= 2\sum\limits_{j,k=1}^{\infty} \frac{1}{jk}\sum\limits_{m=1}^{\infty} \frac{1}{m^2(m+j+k)} \tag{3} \\&= 2\sum\limits_{k,j=1}^{\infty} \frac{1}{kj(j+k)}\left(\zeta(2) - \frac{H_{k+j}}{k+j}\right) \tag{4} \\&= 4\sum\limits_{n=2}^{\infty}\frac{H_{n-1}}{n^2}\left(\zeta(2) - \frac{H_n}{n}\right)\\&= 4\zeta(2)\sum\limits_{n=1}^{\infty}\frac{H_{n-1}}{n^2} + 4\sum\limits_{n=1}^{\infty} \frac{H_n}{n^4}-4\sum\limits_{n=1}^{\infty} \frac{H_n^2}{n^3}\end{align}

Where, in steps $(0)$ and $(3)$ we used the identity, $\displaystyle \frac{H_q}{q} = \sum\limits_{\ell = 1}^{\infty} \frac{1}{\ell (\ell + q)}$. In step $(1)$ we used the symmetry w.r.t. the variables, in step $(2)$ interchanged order of summation and in step $(4)$ made the change of variables $n = j+k$.

Thus, $$ \sum\limits_{n=1}^{\infty} \frac{H_nH_n^{(2)}}{n^2} = 4\zeta(2)\zeta(3) + 7\sum\limits_{n=1}^{\infty} \frac{H_n}{n^4} - 6\sum\limits_{n=1}^{\infty} \frac{H_n^2}{n^3}$$


Definitions

$$\label{eq:1} \displaystyle L^m_n(p,q)=\int^1_0 \frac{\mathrm{Li}_p(x)^m\, \mathrm{Li}_q(x)^n}{x} \, dx \tag{1}$$

$$S_{p,q} = \sum_{n=1}^\infty \frac{H_n^{(p)}}{k^q} \tag{2}$$

$$\mathscr{H}(p,q) =\sum_{k\geq 1}\frac{1}{k^q}\sum_{n\geq 1}\frac{1}{n^{p}(n+k)} \tag{3}$$

$$\displaystyle \mathscr{C}(p, k) =\sum_{n\geq 1}\frac{1}{n^{p}(n+k)}\,\,\, ; \,\,\,\,\mathscr{C}(1, k)=\frac{H_k}{k} \tag{4}$$

Evaluations

Then it is easy to conclude \begin{align} \mathscr{C}(p, k) &=\sum_{n\geq 1}\frac{1}{k\, n^{p-1}}\left( \frac{1}{n}-\frac{1}{n+k}\right)\\ &= \frac{1}{k}\zeta(p)-\frac{1}{k}\mathscr{C}(p-1 , k)\\ &= \frac{1}{k}\zeta(p)-\frac{1}{k^2}\zeta(p-1)+\frac{1}{k^2}\mathscr{C}(p-2 , k)\\ &= \sum_{n=1}^{p-1}(-1)^{n-1}\frac{\zeta(p-n+1)}{k^n}+(-1)^{p-1}\frac{H_k}{k^p} \tag{6} \end{align}

Divide by $k^q$ and sum with respect to $k$ to conclude

\begin{align} \mathscr{H}(p,q)&= \sum_{n=1}^{p-1}(-1)^{n-1}\zeta(p-n+1)\zeta(q+n) -\frac{1}{2}\sum_{n=1}^{{p+q}-2}(-1)^{p-1}\zeta(n+1)\zeta({p+q}-n)\\ &+(-1)^{p-1}\left(1+\frac{{p+q}}{2} \right)\zeta({p+q}+1) \tag{7}\end{align}

Note that by the infinite sums of polylogs we have

$$L^1_1(p,q) = \mathscr{H}(p,q)$$

Starting by

$$\sum_k H_k \, x^{k-1} \mathrm{Li}_q(x)= \frac{\mathrm{Li}_p(x) \mathrm{Li}_q(x)}{x}+\frac{\mathrm{Li}_p(x) \mathrm{Li}_q(x)}{1-x}$$

We can show that \begin{align} \int^1_0\frac{\mathrm{Li}_{p-1}(x) \mathrm{Li}_q(x) \log(1-x)}{x}\, dx+\int^1_0\frac{\mathrm{Li}_p(x) \mathrm{Li}_{q-1}(x) \log(1-x)}{x}\, dx &= \sum_{m=2}^{q-1}(-1)^{m-1} \zeta(q-m+1) S_{p,m} \\&+ (-1)^{q-1} \sum_{k\geq 1} \frac{H_k^{(p)} H_k}{k^q}-\mathscr{H}(p,q)\\& +\zeta(q) \zeta(p+1)-\zeta(q)\mathscr{H}(p-1,1) \tag{8}\end{align}

Let $p=q$ to conclude

\begin{align}2\int^1_0\frac{\mathrm{Li}_q(x) \mathrm{Li}_{q-1}(x) \log(1-x)}{x}\, dx &= \sum_{m=2}^{q-1}(-1)^{m-1} \zeta(q-m+1) S_{q,m} + (-1)^{q-1} \sum_{k\geq 1} \frac{H_k^{(q)} H_k}{k^q}-\mathscr{H}(q,q)\\& +\zeta(q) \zeta(q+1)-\zeta(q)\mathscr{H}(q-1,1) \tag{9} \end{align}

By letting $q=2$

\begin{align}-2\int^1_0\frac{\mathrm{Li}_2(x) \log^2(1-x)}{x}\, dx &= - \sum_{k\geq 1} \frac{H_k^{(2)} H_k}{k^2}-\mathscr{H}(2,2)+\zeta(2) \zeta(3)-\zeta(2)\mathscr{H}(1,1) \end{align}

Note that

$$\int^1_0\frac{\mathrm{Li}_2(x) \log^2(1-x)}{x}\, dx = L^1_2(1,2) $$

Then

$$\displaystyle \sum_{k\geq 1}\frac{H_k^{(2)}H_k }{k^2}=\zeta(2)\zeta(3)+\zeta(5)$$

We could use the same approach to show

$$\displaystyle \sum_{k\geq 1} \frac{H_k^{(2)} H_k}{k^3} =- \frac{97}{12} \zeta(6)+\frac{7}{4}\zeta(4)\zeta(2) + \frac{5}{2}\zeta(3)^2+\frac{2}{3}\zeta(2)^3$$