Compute $\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}$

How to prove

$$\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}=\frac{19}{2}\zeta(3)\zeta(4)-2\zeta(2)\zeta(5)-7\zeta(7)\ ?$$ where $H_n^{(p)}=1+\frac1{2^p}+\cdots+\frac1{n^p}$ is the $n$th generalized harmonic number of order $p$.

This series is very advanced and can be found evaluated in the book (Almost) Impossible Integrals, Sums and Series page 300 using only series manipulations, but luckily I was able to evaluate it using only integration, some harmonic identities and results of easy Euler sums.

Can we prove the equality above in different methods besides series manipulation and the idea of my solution below? All approaches are highly appreciated.

Solution is posted in the answer section.

Thanks

Solution 1:

To compute the target sum, we are going to establish two relations and solve them by elimination.

First Relation:

From here we have $$-\int_0^1x^{n-1}\ln^3(1-x)\ dx=\frac{H_n^3+3H_nH_n^{(2)}+2H_n^{(3)}}{n}$$

Multiply both sides by $\large \frac{H_n}{n^2}$ then sum both sides from $n=1$ to $\infty$ to get

\begin{align} R_1&=\sum_{n=1}^\infty\frac{H_n^4}{n^3}+3\sum_{n=1}^\infty\frac{H_n^2 H_n^{(2)}}{n^3}+2\sum_{n=1}^\infty\frac{H_nH_n^{(3)}}{n^3}=-\int_0^1\frac{\ln^3(1-x)}{x}\sum_{n=1}^\infty\frac{H_n}{n^2}x^n\ dx\\ &=\small{-\int_0^1\frac{\ln^3(1-x)}{x}\left(\operatorname{Li}_3(x)-\operatorname{Li}_3(1-x)+\ln(1-x)\operatorname{Li}_2(1-x)+\frac12\ln x\ln^2(1-x)+\zeta(3)\right)\ dx}\\ &\left\{\text{ let $1-x \mapsto x$ for all integrals but the first one and lets call it $I\ $}\right\}\\ &=\small{-I+\int_0^1\frac{\ln^3x\operatorname{Li}_3(x)}{1-x}-\int_0^1\frac{\ln^4x\operatorname{Li}_2(x)}{1-x}-\frac12\int_0^1\frac{\ln^5x\ln(1-x)}{1-x}-\zeta(3)\int_0^1\frac{\ln^3x}{1-x}\ dx}\\ &=\small{-I+\sum_{n=1}^\infty H_n^{(3)}\int_0^1 x^n \ln^3x-\sum_{n=1}^\infty H_n^{(2)}\int_0^1 x^n \ln^4x+\frac12\sum_{n=1}^\infty H_n\int_0^1 x^n \ln^5x+6\zeta(3)\zeta(4)}\\ &=-I-6\sum_{n=1}^\infty\frac{H_n^{(3)}}{(n+1)^4}-24\sum_{n=1}^\infty\frac{H_n^{(2)}}{(n+1)^5}-60\sum_{n=1}^\infty\frac{H_n}{(n+1)^6}+6\zeta(3)\zeta(4)\\ &=-I-6\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}+6\zeta(7)-24\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^5}+24\zeta(7)-60\sum_{n=1}^\infty\frac{H_n}{n^6}+60\zeta(7)+6\zeta(3)\zeta(4) \end{align}

Then

$$R_1=\sum_{n=1}^\infty\frac{H_n^4}{n^3}+3\sum_{n=1}^\infty\frac{H_n^2 H_n^{(2)}}{n^3}+2\sum_{n=1}^\infty\frac{H_nH_n^{(3)}}{n^3}\\=6\zeta(3)\zeta(4)+90\zeta(7)-I-60\sum_{n=1}^\infty\frac{H_n}{n^6}-24\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^5}-6\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}$$

Second Relation:

From here, we have

$$-\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)\tag{1}$$

Multiply both sides of $(1)$ by $\large-\frac{\ln x}{x}$ then integrate from $x=0$ to $1$ to get \begin{align} S&=\sum_{n=1}^\infty \frac1{n^2}\left(H_n^3-3H_nH_n^{(2)}+2H_n^{(3)}\right)=\int_0^1\frac{\ln^3(1-x)\ln x}{x(1-x)}\ dx\quad \text{let} 1-x\mapsto x\\ &=\int_0^1\frac{\ln^3x\ln(1-x)}{x(1-x)}\ dx=-\sum_{n=1}^\infty H_n\int_0^1 x^{n-1}\ln^3x\ dx=6\sum_{n=1}^\infty\frac{H_n}{n^4}=S\tag{2} \end{align}

Divide both sides of $(1)$ by $x$ then integrate from $x=0$ to $x=y$, we get

$$-\int_0^y\frac{\ln^3(1-x)}{x(1-x)}\ dx=\sum_{n=1}^\infty \frac{y^n}{n}\left(H_n^3-3H_nH_n^{(2)}+2H_n^{(3)}\right)\tag{3}$$

Now multiply both sides of $(3)$ by $-\frac{\operatorname{Li}_2(y)}{y}$ then integrate from $y=0$ to $y=1$ and use the fact that $-\int_0^1 y^{n-1}\operatorname{Li}_2(y)\ dy\overset{IBP}{=}\large\frac{H_n}{n^2}-\frac{\zeta(2)}{n}$, we get

$$\sum_{n=1}^\infty\left(\frac{H_n^3-3H_nH_n^{(2)}+2H_n^{(3)}}{n}\right)\left(\frac{H_n}{n^2}-\frac{\zeta(2)}{n}\right)=\int_0^1\int_0^y\frac{\ln^3(1-x)\operatorname{Li}_2(y)}{xy(1-x)}\ dx\ dy$$

$$\sum_{n=1}^\infty\frac{H_n^4}{n^3}-3\sum_{n=1}^\infty\frac{H_n^2 H_n^{(2)}}{n^3}+2\sum_{n=1}^\infty\frac{H_nH_n^{(3)}}{n^3}-\zeta(2)S=\int_0^1\frac{\ln^3(1-x)}{x(1-x)}\left(\int_x^1\frac{\operatorname{Li}_2(y)}{y}\ dy\right)\ dx$$

Rearranging the terms, we have

\begin{align} R_2&=\sum_{n=1}^\infty\frac{H_n^4}{n^3}-3\sum_{n=1}^\infty\frac{H_n^2 H_n^{(2)}}{n^3}+2\sum_{n=1}^\infty\frac{H_nH_n^{(3)}}{n^3}=\zeta(2)S+\int_0^1\frac{\ln^3(1-x)}{x(1-x)}\left(\zeta(3)-\operatorname{Li}_3(x)\right)\ dx\\ &=\zeta(2)S+\int_0^1\frac{\ln^3(1-x)}{x}\left(\zeta(3)-\operatorname{Li}_3(x)\right) dx+\underbrace{\int_0^1\frac{\ln^3(1-x)}{1-x}\left(\zeta(3)-\operatorname{Li}_3(x)\right) dx}_{IBP}\\ &=\zeta(2)S+\zeta(3)\int_0^1\frac{\ln^3(1-x)}{x}\ dx-I-\frac14\int_0^1\frac{\ln^4(1-x)\operatorname{Li}_2(x)}{x}\ dx, \quad 1-x\mapsto x\\ &=\zeta(2)S+\zeta(3)\int_0^1\frac{\ln^3x}{1-x}\ dx-I-\frac14\int_0^1\frac{\ln^4x\operatorname{Li}_2(1-x)}{1-x}\ dx\\ &=\zeta(2)S-6\zeta(3)\zeta(4)-I-\frac14\int_0^1\frac{\ln^4x}{1-x}\left(\zeta(2)-\ln x\ln(1-x)-\operatorname{Li}_2(x)\right)\ dx\\ &=\zeta(2)S-6\zeta(3)\zeta(4)-I-6\zeta(2)\zeta(5)+\frac14\int_0^1\frac{\ln^5x\ln(1-x)}{1-x}\ dx+\frac14\int_0^1\frac{\ln^4x\operatorname{Li}_2(x)}{1-x}\ dx\\ &=\zeta(2)S-6\zeta(3)\zeta(4)-I-6\zeta(2)\zeta(5)-\frac14\sum_{n=1}^\infty H_n\int_0^1 x^n \ln^5x+\frac14\sum_{n=1}^\infty H_n^{(2)}\int_0^1 x^n\ln^4x\\ &=\zeta(2)S-6\zeta(3)\zeta(4)-I-6\zeta(2)\zeta(5)+30\sum_{n=1}^\infty \frac{H_n}{(n+1)^6}+6\sum_{n=1}^\infty \frac{H_n^{(2)}}{(n+1)^5}\\ &=\zeta(2)S-6\zeta(3)\zeta(4)-I-6\zeta(2)\zeta(5)+30\sum_{n=1}^\infty \frac{H_n}{n^6}-30\zeta(7)+6\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^5}-6\zeta(7)\\ &=\zeta(2)S-6\zeta(3)\zeta(4)-I-6\zeta(2)\zeta(5)-36\zeta(7)+30\sum_{n=1}^\infty \frac{H_n}{n^6}+6\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^5}\\ \end{align}

Substitute the result of $S$ from $(2)$ to get

$$R_2=\sum_{n=1}^\infty\frac{H_n^4}{n^3}-3\sum_{n=1}^\infty\frac{H_n^2 H_n^{(2)}}{n^3}+2\sum_{n=1}^\infty\frac{H_nH_n^{(3)}}{n^3}\\ =-6\zeta(3)\zeta(4)-6\zeta(2)\zeta(5)-36\zeta(7)-I+6\zeta(2)\sum_{n=1}^\infty \frac{H_n}{n^4}+30\sum_{n=1}^\infty \frac{H_n}{n^6}+6\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^5}$$.

Therefore

$$ \sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}=\frac{R_1-R_2}{6}\\ =2\zeta(3)\zeta(4)+21\zeta(7)+\zeta(2)\zeta(5)-\zeta(2)\sum_{n=1}^\infty\frac{H_n}{n^4}-15\sum_{n=1}^\infty\frac{H_n}{n^6}-5\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^5}-\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}$$

We have

$$S_1=\sum_{n=1}^\infty\frac{H_n}{n^4}=3\zeta(5)-\zeta(2)\zeta(3)$$

$$S_2=\sum_{n=1}^\infty\frac{H_n}{n^6}=4\zeta(7)-\zeta(2)\zeta(5)-\zeta(3)\zeta(4)$$

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

$$S_4=\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}=18\zeta(7)-10\zeta(2)\zeta(5)$$

By plugging these results, we get

$$\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}=\frac{19}2\zeta(3)\zeta(4)-2\zeta(2)\zeta(5)-7\zeta(7)$$

Proofs:

The results of $S_1$ and $S_2$ can be obtained from using Euler's identity.

To compute $S_3$, I am going to start with $S_4$:

\begin{align} S_4&=\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}=\sum_{n=1}^\infty\frac1{n^4}\left(\zeta(3)-\sum_{k=1}^\infty\frac1{n+k)^3}\right)\\ &=\zeta(3)\zeta(4)-\sum_{k=1}^\infty\sum_{n=1}^\infty\frac{1}{n^4(n+k)^3}\\ &\small{=\zeta(3)\zeta(4)-\sum_{k=1}^\infty\sum_{n=1}^\infty-\frac{10}{k^6}\left(\frac{1}{n}-\frac{1}{n+k}\right)+\frac6{k^5n^2}+\frac{4}{k^5(n+k)^2}-\frac3{k^4n^3}+\frac1{k^4(n+k)^3}+\frac1{k^3n^4}}\\ &=\zeta(3)\zeta(4)-\sum_{k=1}^\infty-\frac{10H_k}{k^6}+\frac{6\zeta(2)}{k^5}+4\frac{\zeta(2)-H_k^{(2)}}{k^5}-\frac{3\zeta(3)}{k^4}+\frac{\zeta(3)-H_k^{(3)}}{k^4}+\frac{\zeta(4)}{n^3}\\ \color{red}{S_4}&\small{=\zeta(3)\zeta(4)+10\sum_{k=1}^\infty\frac{H_k}{k^6}-6\zeta(2)\zeta(5)-4\zeta(2)\zeta(5)+4\sum_{k=1}^\infty\frac{H_k^{(2)}}{k^5}+3\zeta(3)\zeta(4)-\zeta(3)\zeta(4)+\color{red}{S_4}-\zeta(4)\zeta(3)}\\ &0=2\zeta(3)\zeta(4)-10\zeta(2)\zeta(5)+10\sum_{k=1}^\infty\frac{H_k}{k^6}+4\sum_{k=1}^\infty\frac{H_k^{(2)}}{k^5}\\ \end{align}

Substituting $\displaystyle \sum_{k=1}^\infty\frac{H_k}{k^6}=4\zeta(7)-\zeta(2)\zeta(5)-\zeta(3)\zeta(4)\ $ gives

$$S_3=\sum_{n=1}^\infty\frac{H_k^{(2)}}{k^5}=5\zeta(2)\zeta(5)+2\zeta(3)\zeta(4)-10\zeta(7)$$

If we follow the same approach of evaluating $S_3$ above and start with $\sum_{n=1}^\infty\frac{H_n^{(5)}}{n^2}$, we can find $S_4$ but I am going to present a new way instead.

By the Cauchy product we have,

$$\operatorname{Li}_3^2(x)=\sum_{n=1}^\infty\left(\frac{12H_n}{n^5}+\frac{H_n^{(2)}}{n^4}+\frac{2H_n^{(3)}}{n^3}-\frac{20}{n^6}\right)x^n$$

Divide both sides by $x$ then integrate from $x=0$ to $1$ to get

\begin{align} I&=\sum_{n=1}^\infty\left(\frac{12H_n}{n^6}+\frac{6H_n^{(2)}}{n^5}+\frac{2H_n^{(3)}}{n^4}-\frac{20}{n^7}\right)=\int_0^1\frac{\operatorname{Li}_3^2(x)}{x}\ dx\\ &=\sum_{n=1}^\infty\frac{1}{n^3}\int_0^1x^{n-1}\operatorname{Li}_3(x)\ dx\quad \text{apply integration by parts}\\ &=\sum_{n=1}^\infty\frac{1}{n^3}\left(\frac{\zeta(3)}{n}-\frac{\zeta(2)}{n^2}+\frac{H_n}{n^3}\right)\\ &=\zeta(3)\zeta(4)-\zeta(2)\zeta(5)+\sum_{n=1}^\infty\frac{H_n}{n^6} \end{align}

Rearranging the terms we have

$$\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}=\frac12\zeta(3)\zeta(4)-\frac12\zeta(2)\zeta(5)+10\zeta(7)-\frac{11}{2}\sum_{n=1}^\infty\frac{H_n}{n^6}-3\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^5}$$

Plugging the results:

$$\sum_{n=1}^\infty\frac{H_n}{n^6}=4\zeta(7)-\zeta(2)\zeta(5)-\zeta(3)\zeta(4)$$

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

We get

$$S_4=\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}=18\zeta(7)-10\zeta(2)\zeta(5)$$

The interesting thing about this solution is that I did not use any result of advanced series and that the integral $I$ in $R_1$ and $R_2$ got cancelled out which requires results of wicked series of weight 7 to crack.

Solution 2:

The series $\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}$ can be written as $$\sum_{\substack{n_1\geq n_2\geq 1 \\ n_1\geq n_3\geq 1 \\ n_1\geq n_4\geq 1}}\frac{1}{n_1^3 n_2 n_3 n_4^2},$$ which can be recognized as a linear combination of multiple zeta values of weight $7$.

Multiple zeta values of weight $w$ are series of the form $$\zeta(s_1, \ldots, s_k) = \sum_{n_1 > n_2 > \cdots > n_k > 0} \ \frac{1}{n_1^{s_1} \cdots n_k^{s_k}},$$ such that $s_1,\dots,s_k$ are positive integers and $s_1>1$ such that $s_1+\dots+s_k=w$.

By breaking your sum into parts (depending whether $n_1>n_2>n_3>n_4$ or $n_1>n_2>n_3=n_4$ etc), your sum is equal to the following expression: \begin{align*} \sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}&=2\zeta(3,2,1,1)+2\zeta(3,1,2,1)+2\zeta(3,1,1,2)+2\zeta(5,1,1)+2\zeta(4,2,1)+2\zeta(4,1,2) \\&\quad +\zeta(3,3,1)+2\zeta(3,2,2)+2\zeta(3,1,3)+2\zeta(6,1)+2\zeta(5,2)+2\zeta(4,3) \\&\quad +\zeta(3,4)+\zeta(7). \end{align*}

Now due to the algebraic relations between multiple zeta values (the shuffle and stuffle relations), all multiple zeta values of weight $7$ or less can be computed as a weight preserving $\mathbb{Q}$-linear combination of products of single zeta values. This follows from writing out the relations found in theorems 3.1, 3.2, 3.3 in these lecture notes by Wadim Zudilin. (The weight of the product $\zeta(s_1)\dots\zeta(s_k)$ is the sum $s_1+\dots+s_k$.)

An advantage of this method is that it works in high generality. For example, if one has a series of the form $$\sum_{n=1}^\infty\frac{H_n^{(i_1)}H_n^{(i_2)}\ldots H_n^{(i_k)}}{n^s},$$ with $s, i_1,\dots, i_k$ positive integers and $s>1$, then it can be written as a $\mathbb{Z}$-linear combination of multiple zeta values of weight $w=s+i_1+\dots+i_k$. Therefore, if $w\leq 7$, then the series can be written as a $\mathbb{Q}$-linear combination of products of single zeta values of weight $w$.

Solution 3:

Here is another approach: Again, we are going to estabish two relations and solve for the target sum.

First Relation:

From here we have

$$\int_0^1x^{n-1}\ln^4(1-x)\ dx=\frac1n\left(H_n^4+6H_n^2H_n^{(2)}+8H_nH_n^{(3)}+3\left(H_n^{(2)}\right)^2+6H_n^{(4)}\right)$$

Divide both sides by $n^2$ then sum both sides from $n=1$ to $\infty$ to get

$$R_1=\sum_{n=1}^\infty \frac{1}{n^3}\left(H_n^4+6H_n^2H_n^{(2)}+8H_nH_n^{(3)}+3\left(H_n^{(2)}\right)^2+6H_n^{(4)}\right)\\=\int_0^1\frac{\ln^4(1-x)}{x}\sum_{n=1}^\infty\frac{x^n}{n^2}\ dx=\int_0^1\frac{\ln^4(1-x)\operatorname{Li}_2(x)}{x}\ dx$$

Second Relation:

From here we have

$$\frac{\ln^4(1-x)}{1-x}=\sum_{n=1}^\infty\left(H_n^4-6H_n^2H_n^{(2)}+8H_nH_n^{(3)}+3\left(H_n^{(2)}\right)^2-6H_n^{(4)}\right)x^n$$

Multiply both sides by $\large\frac{\ln^2x}{2x}$ then integrate both sides from $x=0$ to $1$ and using the fact that $ \int_0^1 x^{n-1}\ln^2x\ dx=\large\frac{2}{n^3}$ to get

$$R_2=\sum_{n=1}^\infty \frac{1}{n^3}\left(H_n^4-6H_n^2H_n^{(2)}+8H_nH_n^{(3)}+3\left(H_n^{(2)}\right)^2-6H_n^{(4)}\right)\\=\frac12\int_0^1\frac{\ln^4(1-x)\ln^2x}{x(1-x)}dx\overset{1-x\ \mapsto\ x}{=}\frac12\int_0^1\frac{\ln^4x\ln^2(1-x)}{x(1-x)}dx$$

Then

$$\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}=\frac{R_1-R_2}{12}-\sum_{n=1}^\infty\frac{H_n^{(4)}}{n^3}\\=\frac1{12}\underbrace{\int_0^1\frac{\ln^4(1-x)\operatorname{Li}_2(x)}{x}\ dx}_{\Large I_1}-\frac1{24}\underbrace{\int_0^1\frac{\ln^4x\ln^2(1-x)}{x(1-x)}dx}_{\Large I_2}-\underbrace{\sum_{n=1}^\infty\frac{H_n^{(4)}}{n^3}}_{\Large S}$$

Lets calculate each term and starting with the first one

\begin{align} I_1&=\int_0^1\frac{\ln^4(1-x)\operatorname{Li}_2(x)}{x}\ dx\overset{1-x\ \mapsto\ x}{=}\int_0^1\frac{\ln^4x\operatorname{Li}_2(1-x)}{1-x}\ dx\\ &=\int_0^1\frac{\ln^4x}{1-x}(\zeta(2)-\ln x\ln(1-x)-\operatorname{Li}_2(x))\ dx\\ &=24\zeta(2)\zeta(5)-\int_0^1\frac{\ln^5x\ln(1-x)}{1-x}\ dx-\int_0^1\frac{\ln^4x\operatorname{Li}_2(x)}{1-x}\ dx\\ &=24\zeta(2)\zeta(5)+\sum_{n=1}^\infty H_n\int_0^1 x^n \ln^5x\ dx-\sum_{n=1}^\infty H_n^{(2)}\int_0^1 x^n \ln^4x\ dx\\ &=24\zeta(2)\zeta(5)-120\sum_{n=1}^\infty \frac{H_n}{(n+1)^6}-24\sum_{n=1}^\infty \frac{H_n^{(2)}}{(n+1)^5}\\ &=24\zeta(2)\zeta(5)-120\sum_{n=1}^\infty \frac{H_n}{n^6}+120\zeta(7)-24\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^5}+24\zeta(7)\\ &\boxed{I_1=24\zeta(2)\zeta(5)+144\zeta(7)-120\sum_{n=1}^\infty \frac{H_n}{n^6}-24\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^5}} \end{align}

\begin{align} I_2&=\int_0^1\frac{\ln^4x\ln^2(1-x)}{x(1-x)}dx=\sum_{n=1}^\infty\left(H_n^2-H_n^{(2)}\right)\int_0^1 x^{n-1}\ln^4x\ dx\\ &\boxed{I_2=24\sum_{n=1}^\infty\frac{H_n^2}{n^5}-24\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^5}} \end{align} .

We can easily find the value of $S$ if we 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=4$ and $q=3$, we have

$$\boxed{S=\sum_{n=1}^\infty\frac{H_n^{(4)}}{n^3}=\zeta(3)\zeta(4)+\zeta(7)-\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}}$$.

Collecting the boxed results of $I_1$, $I_2$ and $S$ we get

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

We have the following results:

$$S_1=\sum_{n=1}^\infty\frac{H_n}{n^6}=4\zeta(7)-\zeta(2)\zeta(5)-\zeta(3)\zeta(4)$$

$$S_2=\sum_{n=1}^\infty\frac{H_n^{(3)}}{n^4}=18\zeta(7)-10\zeta(2)\zeta(5)$$

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

$$S_4=\sum_{n=1}^\infty\frac{H_n^2}{n^5}=6\zeta(7)-\zeta(2)\zeta(5)-\frac52\zeta(3)\zeta(4)$$

By substituting these results, we get

$$\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}=\frac{19}2\zeta(3)\zeta(4)-2\zeta(2)\zeta(5)-7\zeta(7)$$

Proofs:

The result of $S_1$ can be obtained from Euler Identity. $S_2$ and $S_3$ are already proved in my previous solution above. As for $S_4$, we compute it as follows

From here, we have

$$\int_0^1x^{n-1}\ln^2(1-x)\ dx=\frac{H_n^2+H_n^{(2)}}{n}$$

Divide both sides by $n^4$ then sum both sides from $n=1$ to $\infty$ to get

\begin{align} \sum_{n=1}^\infty\frac{H_n^2}{n^5}+\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^5}&=\int_0^1\frac{\ln^2(1-x)}{x}\sum_{n=1}^\infty\frac{x^n}{n^4}\ dx\\ &=\int_0^1\frac{\ln^2(1-x)\operatorname{Li}_4(x)}{x}\ dx\\ &=2\sum_{n=1}^\infty\left(\frac{H_n}{n}-\frac1{n^2}\right)\int_0^1x^{n-1} \operatorname{Li}_4(x)\ dx\\ &=2\sum_{n=1}^\infty\left(\frac{H_n}{n}-\frac1{n^2}\right)\left(\frac{\zeta(4)}{n}-\frac{\zeta(3)}{n^2}+\frac{\zeta(2)}{n^3}-\frac{H_n}{n^4}\right)\\ 3\sum_{n=1}^\infty\frac{H_n^2}{n^5}+\sum_{n=1}^\infty\frac{H_n^{(2)}}{n^5}&=2\sum_{n=1}^\infty\frac{H_n}{n^6}+2\zeta(4)\sum_{n=1}^\infty\frac{H_n}{n^2}-2\zeta(3)\sum_{n=1}^\infty\frac{H_n}{n^3}\\ &\quad+2\zeta(2)\sum_{n=1}^\infty\frac{H_n}{n^4}-2\zeta(2)\zeta(5) \end{align}

From Euler's identity, we can obtain the following results:

$$\sum_{n=1}^\infty\frac{H_n}{n^2}=2\zeta(3)$$ $$\sum_{n=1}^\infty\frac{H_n}{n^3}=\frac54\zeta(4)$$ $$\sum_{n=1}^\infty\frac{H_n}{n^6}=4\zeta(7)-\zeta(2)\zeta(5)-\zeta(3)\zeta(4)$$

By plugging these results along with the result of $S_3$, we get the closed form of $S_4$.