Double harmonic series $\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\frac{H_{n+m}^{(p)}}{(n+1)^{q}(m+1)^{r}}$

Question

Do these sums exist in the literature and have been investigated before? The same question for the odd variant, that is $$ \sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\frac{O_{n+m}^{(p)}}{(2n+1)^{q}(2m+1)^{r}}. $$ Here $H_{n}^{(s)}=\sum_{k=1}^{n}\frac{1}{k^{s}}$ and $O_{n}^{(s)}=\sum_{k=1}^{n}\frac{1}{(2k-1)^{s}}$, $H_{0}^{(s)}=O_{0}^{(s)}=0.$

Answer 1

This is not a complete solution (lacking closed expressions) but shows possible first steps towards it.

We shall calculate the sum

$$s(p,q,r) = \sum_{n=1}^{\infty}\sum_{m=1}^{\infty}\frac{H_{n+m}^{(p)}}{m^{q} n^{r}}\tag{1}$$

The sum is related to the sum in the OP considering the relations

$$H_{m+n-2}^{(p)}=H_{m+n}^{(p)}- \frac{1}{k^p} - \frac{1}{(k-1)^p}$$

Notice this basic reference for the calculation of Euler sums: http://algo.inria.fr/flajolet/Publications/FlSa98.pdf (Euler Sums and Contour Integral Representations, Philippe Flajolet and Bruno Salvy)

Integral representation

As a first step I have derived the following integral representation of the sum

$$s_i(p,q,r) = \frac{1}{\Gamma (p)} \int_0^1 \log ^{p-1}\left(\frac{1}{x}\right) \frac{Li_q(1) Li_r(1)-Li_q(x) Li_r(x)}{1-x} \, dx\tag{2}$$

Here $Li_q(x)=\sum_{k=1}^{\infty} \frac{x^k}{k^q}$ is the polylog function.

The derivation uses the representation of the generalized harmonic number

$$H_{m+n}^{(p)}=\sum _{k=1}^{\infty } \left(\frac{1}{k^p}-\frac{1}{(k+m+n)^p}\right)\tag{3}$$

replaces denominators by integrals like

$$k^{-s} = \frac{1}{\Gamma (s)}\int_0^{\infty } t^{s-1} \exp (-t k) \, dt\tag{4}$$

and swaps integration and (double) summation.

The double sum factorizes under the integral, and we have to do sums like

$$\sum_{n=1}^{\infty} \frac{e^{-n t}}{n^p} = Li_p(e^{-t})$$

giving the polylog function, as mentioned.

The convergence of the integral in $(2)$ depends on the behaviour of the Integrand close to $x=1$.

We have for $q=r=2$

$$\frac{Li_2(1) Li_2(1)-Li_2(x) Li_2(x)}{1-x} \underset{x \to 1} \simeq -\frac{1}{3} \pi ^2 (12 x+\log (1-x)-13)\tag{5}$$

and this is integrable at $x=1$. For greater $q$ and $r$ convergence is similar (integrable logarithmic divergence)..

Numerical results

I have found that for numerical purposes the integral is much better suited than the double sum.

For example for $q=2, r=2$ I find for $p=1..5$ the following numericial values in the format $\{p,s_i(p)\}$:

{{1, 5.83536}, {2, 3.75249}, {3, 3.14183}, {4, 2.90331}, {5, 2.79902}}

The odd veriant

Since

$$O_{n}^{(p)}=H_{2n}^{(p)}-\frac{1}{2^p}H_{n}^{(p)}\tag{6}$$

the odd variant is easily obtained from the results shown here.