# 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

• I really don't understand what this is all about. I only post questions I don't know the answers and which has a interesting background in physics Aug 8, 2019 at 16:28
• @stocha On the one hand, it's wrong to use the site as a blog for interesting proofs we've found but I'm not vehemently opposed to the question since someone future person might genuinely want the answer.
– Jam
Aug 12, 2019 at 11:57
• @Jam since when and who decided that it's wrong for one to answer his own question with interesting proofs? On the contrary, this should be encouraged. Aug 12, 2019 at 12:17
• @Ali_Shather I've rolled back the edit since "Very advanced sum" is not a useful or descriptive title for the question. "Advanced" is relative; questions on math.stackexchange is less advanced compared to questions on mathoverflow.
– Jam
Aug 12, 2019 at 12:19
• @Jam I tried to express my impression - "shared with other users?", that the appearance of the site - without rating - has changed over the last years tremendously. I tried to understand, why interest and "upvotes" for contributions in my opinion decreased in the last years that much. I would like to make my contribution to change this. Aug 12, 2019 at 17:01

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.

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$$.

• Nice answer, interesting reference. (+1) Aug 8, 2019 at 20:14

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$$.