Different ways to evaluate $\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{(n+1)(n+2)(n+3)}$

Question

The following question:

How to compute the harmonic series $$\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{(n+1)(n+2)(n+3)}$$

where $H_n=\sum_{k=1}^n\frac{1}{k}$ and $H_n^{(2)}=\sum_{k=1}^n\frac{1}{k^2}$, was proposed by @Tolaso (The original question had been closed and deleted).

My friend Cornel and I managed to compute it in different ways and would like to see other solutions if possible.

---

Cornel's solution:

We start with the following key identity presented in the book (Almost) Impossible Integrals, Sums, and Series) , Sect. $4.19$, pages $290-291$, $$\varphi(n)=\sum_{k=1}^{\infty} \frac{H_k H_k^{(2)}}{(k+1)(k+n+1)}$$ $$\small =2\zeta(3)\frac{H_n}{n}+\frac{\zeta(2)}{2}\frac{H_n^2}{n}-\frac{\zeta(2)}{2}\frac{H_n^{(2)}}{n}-\frac{H_n^{(4)}}{4n}-\frac{(H_n^{(2)})^2}{4n}-\frac{H_n}{n}\sum_{i=1}^n \frac{H_i}{i^2}+\frac{1}{2n}\sum_{i=1}^{n} \frac{H_i^2}{i^2}.\tag1$$ Then, based on $(1)$ we obtain that $$\sum_{k=1}^{\infty} \frac{H_k H_k^{(2)}}{(k+1)(k+2)(k+3)}=\varphi(1)-\varphi(2)=\frac{1}{2}\zeta(3)-\frac{1}{4}\zeta(2)-\frac{1}{32}.$$

End of story

A note: the whole process is completed by using series manipulations only, with no use of integrals.

---

My solution:

Replace $x$ by $xyz$ in the generating function: $$\sum_{n=1}^\infty H_nH_n^{(2)}x^n= \frac{\operatorname{Li}_3(x)+\operatorname{Li}_3(1-x)+\frac12\ln x\ln^2(1-x)-\zeta(2)\ln(1-x)-\zeta(3)}{1-x}=f(x),$$

we get

$$\sum_{n=1}^\infty H_nH_n^{(2)}(xyz)^n=f(xyz).$$

Multiply both sides by $yz^2$ then integrate w.r.t $x$, $y$ and $z$ using the fact$$\int_0^1\int_0^1\int_0^1 x^n y^{n+1}z^{n+2}dxdydz=\frac{1}{(n+1)(n+2)(n+3)}$$

we obtain

$$\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{(n+1)(n+2)(n+3)}=\int_0^1\int_0^1\int_0^1 yz^2f(xyz)dxdydz$$

$$\overset{x=t/y}{=}\int_0^1\left[\int_0^1\int_0^y z^2f(zt)dtdy\right]dz=\int_0^1\left[\int_0^1\int_t^1 z^2f(zt)dydt\right]dz$$

$$=\int_0^1\left[\int_0^1 z^2f(zt)(1-t)dt\right]dz\overset{t=u/z}{=}\int_0^1\int_0^z zf(u)\left(1-\frac{u}{z}\right)dudz$$

$$=\int_0^1\int_u^1 f(u)(z-u)dzdu=\frac12\int_0^1(1-u)^2f(u)du=\frac12\int_0^1 u^2f(1-u)du$$

$$=\frac12\underbrace{\int_0^1 u\operatorname{Li}_3(1-u)du}_{1-u\to u}+\frac12\int_0^1 u\operatorname{Li}_3(u)du$$ $$+\frac12\int_0^1 u\left[\frac12\ln(1-u)\ln^2(u)-\zeta(2)\ln(u)-\zeta(3)\right]du$$

$$=\frac12\int_0^1 \operatorname{Li}_3(u)du+\frac12\int_0^1 u\left[\frac12\ln(1-u)\ln^2(u)-\zeta(2)\ln(u)-\zeta(3)\right]du$$

$$=\frac12\left(\zeta(3)-\zeta(2)+1\right)+\frac12\left(\frac12\zeta(2)-\frac{17}{16}\right)$$

$$=\frac12\zeta(3)-\frac14\zeta(2)-\frac1{32}.$$

Answer 1

Here is another solution as requested.

§1. The sum

$$s_3=\sum_{n=1}^\infty\frac{H_nH_n^{(2)}}{(n+1)(n+2)(n+3)}$$

has two aspects of interest (a) the structure of the demoniminator as a (rising) Pochhammer symbol and (b) the quadratic expression in harmonic numbers in the summand.

As far as I know this is an original problem which was proposed and solved for the first time in the context of the previously deleted question. Please correct me if I'm wrong.

§2. My approach reduces the sum in question to a linear combination of sums over summands of the type $1/p(n)$ and $H/p(n)$ where $H$ is a (generalized) harmonic number and $p(n)$ is a polynomial.

To begin with, let us look the simpler more general sum

$$a_2=\sum_{n=1}^\infty\frac{a(n)}{(n+1)(n+2)}$$

We transform the summand as follows: partial fraction decomposition gives

$$\frac{a(n)}{(n+1)(n+2)}=\frac{a(n)}{1+n}-\frac{a(n)}{2+n}$$

Now we add and subtract a term which is the first term with the index shifted upwards by 1, i.e.

$$\frac{a(n)}{1+n}|_{n\to n+1} = \frac{a(n+1)}{2+n}$$

with the intention of creating terms that telescope under the sum.

Hence

$$\frac{a(n)}{1+n}-\frac{a(n)}{2+n}=\left(\frac{a(n)}{1+n}-\frac{a(n+1)}{2+n}\right) +\left(\frac{a(n+1)}{2+n}-\frac{a(n)}{2+n}\right)$$

we see that the first bracket telescopes term under a sum. The remaining term contains the sequence of the differences of the original seqence.

Taking the finite sum leads to

$$\sum_{n=1}^{m}\frac{a(n)}{(n+1)(n+2)}\\ = \sum_{n=1}^{m}\left(\frac{a(n)}{1+n}-\frac{a(n+1)}{2+n}\right)+\sum_{n=1}^{m}\left(\frac{a(n+1)}{2+n}-\frac{a(n)}{2+n}\right)\\ =\left(\frac{a(1)}{2}-\frac{a(n+1)}{2+n}\right)+\sum_{n=1}^{m}\left(\frac{a(n+1)-a(n)}{2+n}\right)$$

Under the assumption that $\lim_{n\to \infty} {(a(n)/n)}=0$ the corresponding infinite sum reads

$$\sum_{n=1}^{\infty}\frac{a(n)}{(n+1)(n+2)}=\frac{a(1)}{2}+\sum_{n=1}^{\infty}\left(\frac{a(n+1)-a(n)}{2+n}\right)$$

§3. The case $s_3$ can be traced back to $s_2$ setting $a(n)\to \frac{a(n)}{n+3}$.

Applying the same procedure as before we arrive at

$$\sum_{n=1}^{\infty} \frac{a(n)}{(1+n)(2+n)(3+n)} = S_0+S_1+S_2$$

where

$$S_0 = -\frac{5}{25}a(1)+\frac{1}{6}a(2)+\frac{1}{8}a(3)$$ $$S_1 = \frac{1}{2}\sum_{n=1}^{\infty}\frac{a(n+3)-a(n+1)}{n+4}$$ $$S_2 = \sum_{n=1}^{\infty}\frac{a(n+1)-a(n)}{n+3}$$

Specifically for $a(n) = H_n H_n^{(2)}$ and observing that

$$H_{k+1} = H_{k}+\frac{1}{k+1}$$ $$H_{k+1}^{(2)}=H_{k}^{(2)}+\frac{1}{(k+1)^2}$$

we get the following combination of linearized sums, together with their exact sums (based on known results):

$$A_0=\frac{719}{1728}$$

$$\begin{align} A_1& = \sum _{n=1}^{\infty } \left(\frac{H_{n+1}}{2 (n+2)^2 (n+4)}+\frac{H_{n+1}}{2 (n+3)^2 (n+4)}\\ +\frac{H_{n+1}^{(2)}}{2 (n+2) (n+4)}+\frac{H_{n+1}^{(2)}}{2 (n+3) (n+4)}\\ +\frac{1}{2 (n+2)^3 (n+4)}+\frac{1}{2 (n+2)^2 (n+3) (n+4)}\\ +\frac{1}{2 (n+2) (n+3)^2 (n+4)}+\frac{1}{2 (n+3)^3 (n+4)}\right)\\ & =\frac{3 \zeta (3)}{2}+\frac{\pi ^2}{24}-\frac{2969}{1728} \end {align}$$

$$\begin{align} A_2& =\sum _{n=1}^{\infty } \left(-\frac{H_n}{(n+1)^2 (n+3)}-\frac{H_n^{(2)}}{(n+1) (n+3)}-\frac{1}{(n+1)^3 (n+3)}\right)\\ &=-\zeta (3)-\frac{\pi ^2}{12}+\frac{61}{48} \end{align}$$

Adding the three expressions we obtain

$$s_3 = \frac{\zeta (3)}{2}-\frac{\pi ^2}{24}-\frac{1}{32}$$

which is equivalent to the expression obtained in the OP since $\zeta(2) = \frac{\pi^2}{6}$.

$4. Generalization

The obvious generalization asks for the value of

$$s_q(a,b) = \sum_{n=1}^{\infty} \frac{H_n^{(a)} H_n^{(b)}}{(n+1)_q}$$

where $a$ and $b$ can assume the values $0$ or $1$ and

$$H_n^{(a)}=\sum_{k=1}^{n}\frac{1}{k^a}$$

is the generalized harmonic number and

$$(n+1)_q = (n+1)(n+2)...(n+q)$$

is the ascending Pochhammer symbol.

By the procedure described in this answer we can reduce the quadratic sum to a combination of linear sums which can be evaluated.