Double series of Harmonic Numbers

Question

In a solution presented here a series involving the product of Harmonic numbers is involved. The intent of the problem is to determine a form of the series \begin{align} \sum_{n=1}^{\infty} \frac{H_{n}}{(n+1) (n+s+3)^{3}} \end{align} such that \begin{align} \sum_{k=1}^{\infty} \frac{(-1)^{k-1} H_{k}}{k+1} \left(\sum_{n=1}^{\infty} \frac{H_{n}}{(n+1) (n+s+3)^{3}} \right) \end{align} has a closed form. Alternatively, one could calculate the series \begin{align} A_{n} = \sum_{s=1}^{n} \frac{(-1)^{s-1} \, H_{n-s} H_{s} }{ (s+1)(n-s+1) } \end{align} and then find a closed form of the series \begin{align} \sum_{n=1}^{\infty} \frac{A_{n}}{(n+3)^{3}} \end{align}

Answer 1

I have obtained a closed form using Mathematica after some preparatory steps.

EDIT: simplification of problem and solution

Consider the more general sum

$$g(p,s) =\sum _{n=1}^{\infty } \frac{H_n}{(n+1) (n+s+3)^p}\tag{E1} $$

In the OP we had $p=3$.

The sum converges for $p \gt 0$ and, trivially, for an $s$ such that a pole due to summing is avoided.

Now we notice that higher values of $p$ can be generated by differentiating with respect to $s$. Namely

$$g(p+1,s) = -\frac{1}{p} \frac{\partial g(p,s)}{\partial s}\tag{E2}$$

But, luckily, for $p=1$ Mathematica yields immediately the closed expression

$$g(1,s) = \frac{1}{{2 s+4}}(\psi ^{(0)}(s+3)^2+2 \gamma \psi ^{(0)}(s+3)-\psi ^{(1)}(s+3)+\gamma ^2+\frac{\pi ^2}{6})\tag{E3}$$

Using (E3) and (E2) we can generate the closed expressions for a whole family of sums (E1).

Original post

First of all we define the sum to be calculated

$$f_s(s)=\sum _{n=1}^{\infty } \frac{H_n}{(n+1) (n+s+3)^3}\tag{1}$$

Then we generate part of the denominator writing

$$\frac{1}{(n+s+3)^3}=\int_0^{\infty } \frac{1}{2} t^2 \exp (t (-(n+s+3))) \, dt\tag{2}$$

Inserting this in (1) and observing

$$\sum _{n=1}^{\infty } \frac{H_n x^n}{n+1}= \frac{\log ^2(1-x)}{2 x}\tag{3}$$

we get

$$f_i(s)=\frac{1}{4} \int_0^{\infty } t^2 e^{-(s+2) t} \log ^2\left(1-e^{-t}\right) \, dt\tag{4}$$

This integral is evaluated by Mathematica to give the lengthy expression

$$f_\psi(s)\tag{5} = \frac{1}{24 (s+2)^4}\left(\left(6 \gamma ^2+\pi ^2\right) (s+2)+6 (s+2)^3 \psi ^{(0)}(s+3)^4+12 \gamma (s+2)^3 \psi ^{(0)}(s+3)^3-12 (s+2)^3 \psi ^{(0)}(s+2) \psi ^{(0)}(s+3)^3+6 \gamma ^2 (s+2)^3 \psi ^{(0)}(s+3)^2+\pi ^2 (s+2)^3 \psi ^{(0)}(s+3)^2+6 (s+2)^3 \psi ^{(0)}(s+2)^2 \psi ^{(0)}(s+3)^2-24 \gamma (s+2)^3 \psi ^{(0)}(s+2) \psi ^{(0)}(s+3)^2+6 (s+2)^3 \psi ^{(1)}(s+2) \psi ^{(0)}(s+3)^2-36 (s+2)^3 \psi ^{(1)}(s+3) \psi ^{(0)}(s+3)^2+12 \gamma (s+2)^3 \psi ^{(0)}(s+2)^2 \psi ^{(0)}(s+3)+36 (s+2)^3 \psi ^{(0)}(s+2) \psi ^{(1)}(s+3) \psi ^{(0)}(s+3)+6 \gamma ^2 (s+2)^3 \psi ^{(0)}(s+2)^2+\pi ^2 (s+2)^3 \psi ^{(0)}(s+2)^2-2 \left(6 \gamma ^2+\pi ^2\right) (s+2)^3 \psi ^{(0)}(s+2)^2-6 (s+2)^3 \psi ^{(1)}(s+2)^2+18 (s+2)^3 \psi ^{(1)}(s+3)^2+2 (s+2) \left(6 \gamma -\left(6 \gamma ^2+\pi ^2\right) (s+2)\right) \psi ^{(0)}(s+2)+6 (s+2) (1-4 \gamma (s+2)) \psi ^{(1)}(s+2)-6 (s+2)^3 \psi ^{(0)}(s+2)^2 \psi ^{(1)}(s+3)+12 (s+2)^2 (\gamma (s+2)+2) \psi ^{(2)}(s+3)+12 (s+2)^3 \psi ^{(0)}(s+2) \psi ^{(2)}(s+3)-6 (s+2)^3 \psi ^{(3)}(s+3)+36 \gamma \right)$$

Here $\gamma$ is Euler's constant and $\psi ^{(k)}(s)$ is the polygamma function which can be expressed in terms of the harmonic number.

As a numerical check we take values at $s=1.2$

$$N(f_s(1.2) = 0.00949182, N(f_\psi(1.2)) = 0.00949182$$

For some non positive integer values of $s$ we have in the format $\{s, f(s), N(f(s))\}$

$$ \begin{array}{c} \left\{0,\frac{\psi ^{(2)}(3)}{2}-\frac{\pi ^2}{6}+\frac{15}{8}-\frac{\pi ^4}{720},0.0177186\right\} \\ \left\{-1,\psi ^{(2)}(2)-\frac{\pi ^2}{3}-\frac{\pi ^4}{360}+4,0.0354373\right\} \\ \left\{-2,\frac{1}{12} \left(\pi ^2 \psi ^{(2)}(1)-\psi ^{(4)}(1)\right),0.0965512\right\} \\ \left\{-3,\frac{1}{72} \left(342 \psi ^{(2)}(1)+\pi ^4+12 \pi ^2-270 \psi ^{(2)}(2)+540\right),0.593724\right\} \\ \end{array} \tag{6}$$

For negative integer values $s\le -4$ the function $f(s)$ has triple poles with residues $H_k/(k+1)$ where $k=-s-3$.

The asymptotic behaviour is given by

$$f(s\to \infty) \sim \frac{1}{2 s^3}(\log(s)^2+2 \gamma \log(s)-3 \log(s) +\zeta (2)+\gamma ^2-3 \gamma +1)\tag{7}$$

In order to avoid typing errors in further use here is the complete Mathematica expression of $f(s)$

f[s_] = 1/(24 (2 + 
     s)^4) (36 EulerGamma + (6 EulerGamma^2 + \[Pi]^2) (2 + s) + 
    2 (2 + s) (6 EulerGamma - (6 EulerGamma^2 + \[Pi]^2) (2 + 
          s)) PolyGamma[0, 2 + s] + 
    6 EulerGamma^2 (2 + s)^3 PolyGamma[0, 
      2 + s]^2 + \[Pi]^2 (2 + s)^3 PolyGamma[0, 2 + s]^2 - 
    2 (6 EulerGamma^2 + \[Pi]^2) (2 + s)^3 PolyGamma[0, 2 + s]^2 + 
    12 EulerGamma (2 + s)^3 PolyGamma[0, 2 + s]^2 PolyGamma[0, 
      3 + s] + 
    6 EulerGamma^2 (2 + s)^3 PolyGamma[0, 
      3 + s]^2 + \[Pi]^2 (2 + s)^3 PolyGamma[0, 3 + s]^2 - 
    24 EulerGamma (2 + s)^3 PolyGamma[0, 2 + s] PolyGamma[0, 
      3 + s]^2 + 
    6 (2 + s)^3 PolyGamma[0, 2 + s]^2 PolyGamma[0, 3 + s]^2 + 
    12 EulerGamma (2 + s)^3 PolyGamma[0, 3 + s]^3 - 
    12 (2 + s)^3 PolyGamma[0, 2 + s] PolyGamma[0, 3 + s]^3 + 
    6 (2 + s)^3 PolyGamma[0, 3 + s]^4 + 
    6 (2 + s) (1 - 4 EulerGamma (2 + s)) PolyGamma[1, 2 + s] + 
    6 (2 + s)^3 PolyGamma[0, 3 + s]^2 PolyGamma[1, 2 + s] - 
    6 (2 + s)^3 PolyGamma[1, 2 + s]^2 - 
    6 (2 + s)^3 PolyGamma[0, 2 + s]^2 PolyGamma[1, 3 + s] + 
    36 (2 + s)^3 PolyGamma[0, 2 + s] PolyGamma[0, 3 + s] PolyGamma[1, 
      3 + s] - 
    36 (2 + s)^3 PolyGamma[0, 3 + s]^2 PolyGamma[1, 3 + s] + 
    18 (2 + s)^3 PolyGamma[1, 3 + s]^2 + 
    12 (2 + s)^2 (2 + EulerGamma (2 + s)) PolyGamma[2, 3 + s] + 
    12 (2 + s)^3 PolyGamma[0, 2 + s] PolyGamma[2, 3 + s] - 
    6 (2 + s)^3 PolyGamma[3, 3 + s]);