On the Euler sum $\sum \limits_{n=1}^{\infty} \frac{H_n^{(4)} H_n^2}{n^6}$

Here is an Euler sum I ran into.

$$\mathcal{S} = \sum_{n=1}^{\infty} \frac{\mathcal{H}_n^{(4)} \mathcal{H}_n^2}{n^6}$$

where $\mathcal{H}_n^{(s)}$ is the generalised harmonic number of order $s$. I have no idea to what this evaluates neither do I have the appropriate techniques to begin cracking it.

We have seen quite plenty Euler sums but I don't think we have seen something like this before. Correct me if I am wrong!

Any help?

Edit: Maybe this link ( inverse sin ) is helpful for an integral approach.

Update (by Editor): By stuffle relations of Multiple Zeta Values the result is: >$$S=\small -\frac{3}{20} \pi ^4 \zeta(6,2)-\frac{5}{3} \pi ^2 \zeta(8,2)+\frac{143}{4} \zeta(10,2)-6\zeta(8,2,1,1)+\frac{\zeta (3)^4}{3}-\frac{23 \pi ^6 \zeta (3)^2}{1620}-\frac{8}{45} \pi ^4 \zeta (5) \zeta (3)-\frac{31}{3} \pi ^2 \zeta (7) \zeta (3)+\frac{2531 \zeta (9) \zeta (3)}{18}-\frac{9 \pi ^2 \zeta (5)^2}{2}+\frac{1115 \zeta (5) \zeta (7)}{8}-\frac{964213 \pi ^{12}}{8756748000}$$

A partial answer, designed to simplify further attempts. By Thm 4.2 of Flajolet and Salvy the Euler sum $S_{42,6}=\sum_{n\geq 1}\frac{H_n^{(4)}H_n^{(2)}}{n^6}$ can be written in terms of values of the $\zeta$ function. It follows that it is enough to evaluate $$ \mathcal{J} = \sum_{n\geq 1}\frac{H_n^{(4)}\left(H_n^2-H_n^{(2)}\right)}{n^6}=\sum_{n\geq 1}\frac{H_n^{(4)}}{n^5}\int_{0}^{1}\log^2(1-x)x^{n-1}\,dx. $$ and in order to deal with $H_n^{(4)}$ we may steal a very interesting idea from nospoon:

$$ \arcsin(x)^6 = \frac{45}{8}\sum_{k\geq 1}\left(\left(H_{k-1}^{(2)}\right)^2-\color{red}{H_{k-1}^{(4)}}\right)\frac{4^k x^{2k}}{\color{red}{k^2}\binom{2k}{k}} $$

see also equation $(21)$ here. It might be useful to run an implementation of the PSLQ algorithm for reaching a conjectural closed form for $S_{411,6}$; in the meanwhile I'll keep digging in the literature.

The question can also be solved by evaluating $$ \sum_{a,b,c\geq 1}\frac{(2 a+b) (b+2 c) \left(b^2+2 b c+2 c^2\right)}{a^2 b^4 (a+b)^2 c^4 (b+c)^4} $$ or $$ \iint_{(0,1)^2}\frac{\text{Li}_4(xy)\log^4(y)\log^2(1-x)}{xy(1-xy)}\,dx\,dy ,$$ good luck. Apart from the reflection formulas for the tetralogarithm, we may notice that

$$\iint_{(0,1)^2}\frac{\text{Li}_4(xy)\log^4(y)\log^2(1-x)}{xy}\,dx\,dy = 24\int_{0}^{1}\frac{\text{Li}_9(x)\log^2(1-x)}{x}\,dx $$ depends on the Euler sum $S_{2,10}$, which might be irreducible just like $S_{2,6}$ and $S_{2,8}$. So, unless we have some cancellation, we cannot be sure that the original Euler sum with weight $12$ can be expressed in terms of values of the $\zeta$ function only. Such impossibility for $\zeta(6,4,1,1)$ has been proved by Houches here.

This is not going to be a full answer however I believe that this approach is viable and can lead to an answer. By using the generating function of squares of harmonic numbers Integral involving a dilogarithm versus an Euler sum. and also by using the following integral representation of the generalized harmonic number: \begin{equation} H^{(p)}_m = (-1)^{p-1} \int\limits_0^1 \sum\limits_{l=1}^p Li_l(\xi) \cdot \frac{[\log(\xi)]^{p-l}}{(p-l)!}\cdot (-1)^{l-1} \cdot m \xi^{m-1} d\xi \end{equation} we arrive at the following integral representation: \begin{eqnarray} &&\sum\limits_{n=1}^\infty [H_n]^2 \cdot H_n^{(4)} \cdot \frac{1}{n^6} = \\ && \frac{1}{3! 4!} \sum\limits_{j=1}^4 (-1)^j \binom{3}{j-1} (j-1)! \int\limits_{(0,1)^2}[\log(t)]^{4-j} [\log(\xi)]^4 Li_j(t) \cdot \frac{Li_2(\xi t)+\log(1-\xi t)^2}{t \xi (1-t \xi)} dt d\xi =\\ &&-\frac{1}{3! 4!} \int\limits_{[0,1]^3} \log(\xi)^4 \log(t \eta)^3 \frac{(Li_2(\xi t)+\log(1-\xi t)^2)}{\xi(1-\xi t)(1-t \eta)} d\xi dt d\eta=\\ &&-\frac{1}{3! 4!} \int\limits_0^1 \int\limits_0^\xi \int\limits_0^{t/\xi} \log(\xi)^4 \log(\eta)^3 \frac{(Li_2(t)+\log(1-t)^2)}{\xi(1-t)(1-\eta) t} d\eta dt d\xi=\\ &&-\frac{1}{3! 5!} \int\limits_0^1 \int\limits_0^1 (-\log(t)^5 1_{\eta < t} + (\log(t/\eta))^5-\log(t)^5) 1_{\eta > t}) \log(\eta)^3 \frac{Li_2(t)+\log(1-t)^2}{(1-t)(1-\eta) t} d\eta dt=\\ &&-\frac{1}{3! 5!} \int\limits_0^1 \left(40320 \text{Li}_9(t)-120 \text{Li}_6(t) \log ^3(t)+2160 \text{Li}_7(t) \log ^2(t)-15120 \text{Li}_8(t) \log (t)+\right.\\ &&\left.-120 \zeta (5) \log ^4(t)-7200 \zeta (7) \log ^2(t)-\frac{80}{63} \pi ^6 \log ^3(t)-\frac{8}{3} \pi ^8 \log (t)-40320 \zeta (9) \right) \frac{Li_2(t)+\log(1-t)^2}{(1-t) t} dt \end{eqnarray} In the third line from the top we used the integral representation of $L_j(t)$ and carried out the sum over $j$. In the forth line from the top we substituted for $(\xi,t \xi, t\eta)$. In the fifth line we carried out the integral over $\xi$ and in the sixth line the integral over $\eta$. Now it is clear that many of the terms in the bottom line can be integrated in closed form. We will complete this asap.

Meanwhile we only mention that the sum in question converges slowly. In fact basing on the last integral our sum equals $1.0440866271501388232534111083784540300372112557593..$. On the other hand if we were to compute the sum from its definition we need to sum up ten thousand terms in order to reach an accuracy of eighteen decimal digits.

Update: Now, we deal with the last integral above. It splits into different terms which all fall into three groups only. Let us start with the easiest term and then end with the hardest one. We have: \begin{eqnarray} &&\int\limits_0^1 [\log(t)]^p \frac{Li_2(t)+\log(1-t)^2}{(1-t)t}dt=\\ &&(-1)^p p! \zeta(p+3) + \\ &&(-1)^p p! (2 {\bf H}^{(1)}_{p+2}(1)+ {\bf H}^{(2)}_{p+1}(1)-3 \zeta(p+3))\cdot 1_{p\ge 2} -\frac{1}{2} \zeta(2)^2 \cdot 1_{p=1}+\\ &&\left(-\frac{1}{3} \Psi^{(p+2)}(1) + \frac{1}{2} \sum\limits_{j=1}^{p-1} \binom{p}{j} (\Psi^{(j+1)}(1) \Psi^{(p-j)}(1)+\Psi^{(j)}(1) \Psi^{(p+1-j)}(1))+\sum\limits_{1 \le j < j_1 \le p-1} \frac{p!}{j! (j_1-j)!(p-j_1)!} \Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(p-j_1)}(1)\right) \end{eqnarray} for $p=1,2,3,\cdots$. Here $\Psi$ is the polygamma function.

Since in our case the biggest value of $p$ equals four all the generating functions ${\bf H}^{(.)}_{.}(1)$ are expressed through single zeta functions only as seen in Calculating alternating Euler sums of odd powers.

The remaining terms are as follows. Firstly the term of medium difficulty: \begin{eqnarray} &&\int\limits_0^1 \left(Li_{9}(t) -Li_9(1)\right) \frac{Li_2(t)+\log(1-t)^2}{(1-t)t}dt=\cdots\\ \end{eqnarray} and then the hardest terms: \begin{eqnarray} &&\int\limits_0^1 Li_{9-p}(t) [\log(t)]^p \frac{Li_2(t)+\log(1-t)^2}{(1-t)t}dt=\cdots\\ \end{eqnarray} for $p=1,2,3$.

We will deal with those terms later on. Now we can only note that these terms will involve multiple zeta functions of depth three and as such may not necessarily reduce to single zeta values.