# Triple Euler sum result $\sum_{k\geq 1}\frac{H_k^{(2)}H_k }{k^2}=\zeta(2)\zeta(3)+\zeta(5)$

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 ?

• Related: math.stackexchange.com/q/606070 Dec 20, 2013 at 22:47
• @Grigory M. That was my topic and the link of the paper I provided there gave a citation to a paper by Borwein and R. Girgensohn, 'Evaluation of triple Euler sums paper'. The question in hand is a special case of a general formula but I was wondering whether we can solve it without that. The author used combintorics which seemed complicated for me. Dec 20, 2013 at 23:53
• This is a really hard nut to crack. I only noticed this question yesterday, and I have made several attempts, but I still haven't found the proper approach.
– robjohn
Dec 28, 2013 at 20:48
• Please post the symbolic input you entered. Dec 29, 2013 at 3:39
• $\sum_{n=1}^\infty\frac{H_n^{(2)}H_n}{n^2}$ and $\sum_{n=1}^\infty\frac{H_n^3}{n^2}$ can be evluated together using two nice identities and in a very short way. May 9, 2019 at 16:49

Here's a derivation that, while fairly long, is self-contained and uses only basic series manipulation techniques, like partial fractions decomposition, telescoping, swapping the order of summation, etc. It leans heavily on ideas from Borwein and Girgensohn's paper "Evaluation of Triple Euler Sums" (Electronic Journal of Combinatorics 3(1) 1996).

First, some notation. Define the multiple zeta functions by \begin{align} \zeta_N(a) &= \sum_{x=1}^N \frac{1}{x^a}, \:\:\: \zeta_N(a,b) = \sum_{x=1}^N \sum_{y=1}^{x-1} \frac{1}{x^a y^b}, \:\:\: \zeta_N(a,b,c) = \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^a y^b z^c}, \\ \zeta(a,b) &= \sum_{x=1}^{\infty} \sum_{y=1}^{x-1} \frac{1}{x^a y^b}, \:\:\: \zeta(a,b,c) = \sum_{x=1}^{\infty} \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^a y^b z^c}. \end{align}

We will need the following symmetry relation, as well as expressions for $\zeta(4,1)$ and $\zeta(2,2,1) + \zeta(2,1,2)$. Proofs for all of these are given at the end of the post. \begin{align} \zeta_N(a,b) + \zeta_N(b,a) &= \zeta_N(a) \zeta_N(b) - \zeta_N(a+b) \tag{1}\\ \zeta(4,1) &= \zeta(5) - \zeta(3,2) - \zeta(2,3) \tag{2}\\ \zeta(2,2,1) + \zeta(2,1,2) &= \zeta(2,3) + \zeta(3,2) \tag{3} \end{align}

Given these, we have

The Main Proof: \begin{align} \sum_{k=1}^{\infty} \frac{H^{(2)}_k H_k}{k^2} &= \sum_{k=1}^{\infty} \frac{H^{(2)}_{k-1} H_{k-1}}{k^2} + \sum_{k=1}^{\infty} \frac{H^{(2)}_{k-1}}{k^3} + \sum_{k=1}^{\infty} \frac{H_{k-1}}{k^4} + \sum_{k=1}^{\infty} \frac{1}{k^5} \\ &= \sum_{k=1}^{\infty} \frac{H^{(2)}_{k-1} H_{k-1}}{k^2} + \zeta(3,2) + \zeta(4,1) + \zeta(5). \end{align} The most complicated sum is the first, so let's look at that more closely. \begin{align} \sum_{k=1}^{\infty} \frac{H^{(2)}_{k-1} H_{k-1}}{k^2} &= \sum_{k=1}^{\infty} \frac{1}{k^2} \zeta_{k-1}(2) \zeta_{k-1}(1) \\ &= \sum_{k=1}^{\infty} \frac{1}{k^2} (\zeta_{k-1}(2,1) + \zeta_{k-1}(1,2) + \zeta_{k-1}(3)), \text{ by (1)} \\ &= \zeta(2,2,1) + \zeta(2,1,2) + \zeta(2,3), \text{ by definition of the multiple zeta functions} \\ &= 2\zeta(2,3) + \zeta(3,2), \text{ by (3)}. \end{align} Thus \begin{align} \sum_{k=1}^{\infty} \frac{H^{(2)}_k H_k}{k^2} &= 2 \zeta(2,3) + \zeta(3,2) + \zeta(3,2) + \zeta(5) - \zeta(3,2) - \zeta(2,3) + \zeta(5), \text{ by (2)} \\ &= \zeta(2,3) + \zeta(3,2) + 2 \zeta(5) \\ &= \zeta(2) \zeta(3) - \zeta(5) + 2 \zeta(5), \text{ by (1)} \\ &= \zeta(2) \zeta(3) + \zeta(5). \end{align}

Proof of (1): \begin{align} \zeta_N(a,b) + \zeta_N(b,a) &= \sum_{x=1}^N \sum_{y=1}^{x-1} \frac{1}{x^a y^b} + \sum_{x=1}^N \sum_{y=1}^{x-1} \frac{1}{x^b y^a} \\ &= \sum_{y=1}^N \sum_{x=y+1}^N \frac{1}{x^a y^b} + \sum_{x=1}^N \sum_{y=1}^{x-1} \frac{1}{x^b y^a}, \\ & \:\:\:\:\:\: \text{ swapping the order of summation on the first sum} \\ &= \sum_{x=1}^N \sum_{y=x+1}^N \frac{1}{y^a x^b} + \sum_{x=1}^N \sum_{y=1}^{x-1} \frac{1}{x^b y^a}, \text{ relabeling variables on the first sum} \\ &= \sum_{x=1}^N \sum_{y=1}^N \frac{1}{y^a x^b} - \sum_{x=1}^N \frac{1}{x^{a+b}}, \text{ combining sums} \\ &= \zeta_N(a) \zeta_N(b) - \zeta_N(a+b). \square \end{align}

Proof of (2): \begin{align} \zeta(4,1) &= \sum_{x=1}^{\infty} \sum_{y=1}^{x-1} \frac{1}{x^4 y} \\ &= \sum_{x=1}^{\infty} \sum_{y=1}^{x-1} \frac{1}{x^4 (x-y)}, \text{ reindexing the second sum} \\ &= \sum_{x=1}^{\infty} \sum_{y=1}^{x-1} \left(-\frac{1}{x^4 y} - \frac{1}{x^3 y^2} - \frac{1}{x^2y^3} - \frac{1}{x y^4} + \frac{1}{(x-y)y^4}\right), \\ &\:\:\:\:\: \text{ by partial fractions decomposition}\\ &= - \zeta(4,1) - \zeta(3,2) - \zeta(2,3) + \sum_{x=1}^{\infty} \sum_{y=1}^{x-1} \left(\frac{1}{(x-y)y^4} - \frac{1}{x y^4} \right) \\ &= - \zeta(4,1) - \zeta(3,2) - \zeta(2,3) + \sum_{x=1}^{\infty} \sum_{y=1}^{x-1} \frac{1}{y^4} \left(\frac{1}{x-y} - \frac{1}{x} \right) \\ &= - \zeta(4,1) - \zeta(3,2) - \zeta(2,3) + \sum_{y=1}^{\infty} \frac{1}{y^4} \sum_{x=y+1}^{\infty} \left(\frac{1}{x-y} - \frac{1}{x} \right), \\ & \:\:\:\:\: \text{ swapping the order of summation} \\ &= - \zeta(4,1) - \zeta(3,2) - \zeta(2,3) + \sum_{y=1}^{\infty} \frac{1}{y^4} \sum_{x=1}^y \frac{1}{x}, \text{ as the sum telescopes} \\ &= - \zeta(4,1) - \zeta(3,2) - \zeta(2,3) + \zeta(4,1) + \zeta(5) \\ &= \zeta(5) - \zeta(3,2) - \zeta(2,3). \square \end{align}

For the proof of (3), we need the following additional symmetry result: $$\zeta_N(a,b,c) + \zeta_N(a,c,b) + \zeta_N(c,a,b) = \zeta_N(c) \zeta_N(a,b) - \zeta_N(a,b+c) - \zeta_N(a+c,b) \tag{4}$$

Proof of (4): \begin{align} &\zeta_N(a,b,c) + \zeta_N(a,c,b) + \zeta_N(c,a,b) \\ &=\sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^a y^b z^c} + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^a y^c z^b} + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^c y^a z^b} \\ &= \sum_{x=1}^N \sum_{z=1}^{x-1} \sum_{y=z+1}^{x-1} \frac{1}{x^a y^b z^c} + \sum_{y=1}^N \sum_{x=y+1}^N \sum_{z=1}^{y-1}\frac{1}{x^a y^c z^b} + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^c y^a z^b}, \\ &\:\:\:\:\:\text{ swapping order of summation on the first two sums} \\ &= \sum_{z=1}^N \sum_{x=z+1}^N \sum_{y=z+1}^{x-1} \frac{1}{x^a y^b z^c} + \sum_{y=1}^N \sum_{x=y+1}^N \sum_{z=1}^{y-1}\frac{1}{x^a y^c z^b} + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^c y^a z^b}, \\ &\:\:\:\:\:\text{ swapping order of summation on the first sum} \\ &= \sum_{x=1}^N \sum_{y=x+1}^N \sum_{z=x+1}^{y-1} \frac{1}{x^c y^a z^b} + \sum_{x=1}^N \sum_{y=x+1}^N \sum_{z=1}^{x-1}\frac{1}{x^c y^a z^b} + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^c y^a z^b}, \\ &\:\:\:\:\: \text{ relabeling variables on the first two sums} \\ &= \sum_{x=1}^N \sum_{y=x+1}^N \sum_{z=1}^{y-1} \frac{1}{x^c y^a z^b} - \sum_{x=1}^N \sum_{y=x+1}^N \frac{1}{x^{b+c} y^a} + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1}\frac{1}{x^c y^a z^b}, \\ &\:\:\:\:\: \text{ combining the first two sums} \\ &= \sum_{x=1}^N \sum_{y=1}^N \sum_{z=1}^{y-1} \frac{1}{x^c y^a z^b} - \sum_{x=1}^N \sum_{z=1}^{y-1} \frac{1}{x^{a+c} z^b} - \sum_{y=1}^N \sum_{x=1}^{y-1} \frac{1}{x^{b+c} y^a}, \\ &\:\:\:\:\:\text{ combining the first and third sums and swapping the order of summation on the second} \\ &= \zeta_N(c) \zeta_N(a,b) - \zeta_N(a+c,b) - \zeta_N(a,b+c). \square \end{align}

Proof of (3): \begin{align} \zeta_N(2,2,1) &= \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \frac{1}{x^2 y^2 z} \\ &= \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \frac{1}{x^2 y^2 (y-z)}, \text{ reindexing on the third sum} \\ &= \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \left( -\frac{1}{x^2 y z^2} - \frac{1}{x^2 y^2 z} + \frac{1}{x^2(y-z)z^2} \right), \\ &\:\:\:\:\: \text{ by partial fractions decomposition} \\ &= - \zeta_N(2,1,2) - \zeta_N(2,2,1) + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \frac{1}{x^2(y-z)z^2} \tag{5}. \\ \end{align} Now, let's look at the third expression in (5). \begin{align} &\sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \frac{1}{x^2(y-z)z^2} \\ &= \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{x-y-1} \frac{1}{x^2(x-y-z)z^2}, \text{ reindexing the second sum} \\ &= \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=y+1}^{x-1} \frac{1}{x^2(x-z)(z-y)^2}, \text{ reindexing the third sum} \\ &= \sum_{x=1}^N \sum_{z=1}^{x-1} \sum_{y=1}^{z-1} \frac{1}{x^2(x-z)(z-y)^2}, \text{ swapping the order of summation} \\ &= \sum_{x=1}^N \sum_{z=1}^{x-1} \sum_{y=1}^{z-1} \frac{1}{x^2(x-z)y^2}, \text{ reindexing the third sum} \\ &= \sum_{x=1}^N \sum_{z=1}^{x-1} \sum_{y=1}^{z-1} \left(-\frac{1}{x y^2 z^2} - \frac{1}{x^2 y^2 z} + \frac{1}{(x-z)y^2 z^2} \right), \text{ by partial fractions decomposition} \\ &= - \zeta_N(1,2,2) - \zeta_N(2,1,2) + \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \frac{1}{(x-y)y^2 z^2} \tag{6}, \text{ relabeling variables}. \end{align} Let's look at the third expression in (6). \begin{align} &\sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \frac{1}{(x-y)y^2 z^2} \\ &= \sum_{x=1}^N \sum_{y=1}^{x-1} \sum_{z=1}^{y-1} \frac{1}{(x-y)y^2 z^2} + \sum_{x=1}^N \sum_{y=N+1-x}^N \sum_{z=1}^{y-1} \frac{1}{x y^2 z^2} - \sum_{x=1}^N \sum_{y=N+1-x}^N \sum_{z=1}^{y-1} \frac{1}{x y^2 z^2} \\ &= \left(\sum_{x=1}^N \frac{1}{x}\right) \left(\sum_{y=1}^N \sum_{z=1}^{y-1} \frac{1}{y^2 z^2} \right) - \sum_{x=1}^N \sum_{y=N+1-x}^N \sum_{z=1}^{y-1} \frac{1}{x y^2 z^2}, \\ &\:\:\:\:\: \text{ via the finite sum version of the Cauchy product} \\ &= \zeta_N(1) \zeta_N(2,2) - e_N(1,2,2), \tag{7} \\ \end{align} where $$e_N(1,2,2) = \sum_{x=1}^N \sum_{y=N+1-x}^N \sum_{z=1}^{y-1} \frac{1}{x y^2 z^2}.$$ Putting (5), (6), and (7) together, we have \begin{align} \zeta_N(2,2,1) =& - \zeta_N(2,1,2) - \zeta_N(2,2,1) - \zeta_N(1,2,2) - \zeta_N(2,1,2) + \zeta_N(1) \zeta_N(2,2) \\ &- e_N(1,2,2), \\ \zeta_N(2,2,1) + \zeta_N(2,1,2) &= - \zeta_N(1) \zeta_N(2,2) + \zeta_N(2,3) + \zeta_N(3,2) + \zeta_N(1) \zeta_N(2,2) \\ &- e_N(1,2,2), \text{ by (4)} \\ =& \zeta_N(2,3) + \zeta_N(3,2) - e_N(1,2,2). \\ \end{align} All that remains to complete the proof of (3) is to show that $e_N(1,2,2) \to 0$ as $N \to \infty$. We have \begin{align} e_N(1,2,2) &= \sum_{x=1}^N \sum_{y=N+1-x}^N \sum_{z=1}^{y-1} \frac{1}{x y^2 z^2} \\ &\leq \sum_{x=1}^N \sum_{y=N+1-x}^N \sum_{z=1}^N \frac{1}{x y^2 z^2} \\ &= \zeta_N(2) \sum_{x=1}^N \sum_{y=N+1-x}^N \frac{1}{x y^2} \\ &= \zeta_N(2) \sum_{y=1}^N \sum_{x=N+1-y}^N \frac{1}{x y^2}, \text{ swapping the order of summation} \\ &\leq \zeta_N(2) \sum_{y=1}^N \frac{1}{y^2} \sum_{x=N+1-y}^N \frac{1}{N+1-y} \\ &= \zeta_N(2) \sum_{y=1}^N \frac{1}{y^2} \frac{y}{N+1-y} \\ &= \zeta_N(2) \sum_{y=1}^N \frac{1}{y (N+1-y)}\\ &= \zeta_N(2) \frac{1}{N+1}\sum_{y=1}^N \left(\frac{1}{y} + \frac{1}{N+1-y} \right), \text{ by partial fractions decomposition} \\ &= \zeta_N(2) \frac{2}{N+1} \zeta_N(1), \end{align} which goes to $0$ as $N \to \infty$, since $\zeta_N(1) = O(\log N)$ and $\zeta_N(2) = O(1)$. $\square$

• I just spent two hours typesetting this, and I'm tired of looking at it. If you spot any typos, let me know, and I'll fix them when I get the chance. Jan 2, 2014 at 5:27
• Hey thanks. I really appreciate that time you spent on that . Just as a note : It can be proved that the complexity of evaluating that sum is equivalent to that of $$\sum_{k\geq 1}\frac{H_k^2}{k^3}$$ Jan 2, 2014 at 6:11
• (+1) I was trying to solve this, but I am taking a break for a while. I did come up with some nice relations for double Euler Sums. I can now compute $\sum\limits_{k=1}^\infty\frac{H_n^{(a)}}{n^b}$, so I consider the work I've done on this problem productive.
– robjohn
Jan 4, 2014 at 2:13
• @robjohn have you posted your solution for $\sum\limits_{k=1}^\infty\frac{H_n^{(a)}}{n^b}$ ? I found one proof in some article ( i forgot the name) where they used matrices. It was really boring and complicated. Their work is appreciated though. Jun 8, 2022 at 10:03

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.

• (+1) The Taylor series for $\left(\arcsin x\right)^4$ is a little gem and I have just realized how horribly efficient this approach is for tackling $$\sum_{k\geq 1}\left(\zeta(2)-H_k^{(2)}\right)^3.$$ Congrats. Jan 4, 2018 at 2:30
• awesome solution. rich of info Jul 8, 2019 at 16:21

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

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

• (+1) Good way to do it. Jul 8, 2019 at 12:33
• Thank you glad you like it. Jul 8, 2019 at 18:46

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

• Cornel always has the most unique elegant solutions and I learned a lot from him. He is a good friend of mine. Jul 10, 2019 at 5:15
• @AliShather He told me you are a great talent in the work with these series. I work hard to enter the same league. :) Jul 10, 2019 at 6:05
• thank you. its fun man you will like it but it requires a lot of patience Jul 10, 2019 at 6:26

$$\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)$$

• (+1) for the work with integrals. The double sum you start with (the version with $H_n$ in the numerator) is a critical sum that appears in the calculations by series manipulations presented in (Almost) Impossible Integrals, Sums, and Series (see the top of page $399$). Jul 12, 2020 at 8:19
• @user97357329 thank you :) . its good to see you back I think you've been away from this site for quite awhile. Jul 12, 2020 at 8:40
• very rarely around. Looks like you lately have learned a lot and become a king in the realm of harmonic series on MSE. ;) Jul 12, 2020 at 8:48
• @user97357329 I would like to see you participate in solving problems here and also in posting nice questions I learned nice tricks from your postings/solutions. Yes, thank you I've been practicing a lot lately and when i get stuck with something I ask Cornel he is the real king :) Jul 12, 2020 at 8:57

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.