Polygamma function series: $\displaystyle\sum_{k=1}^{\infty }\left(\Psi^{(1)}(k)\right)^2$

Applying the Copson’s inequality, I found:
$$S=\displaystyle\sum_{k=1}^{\infty }\left(\Psi^{(1)}(k)\right)^2\lt\dfrac{2}{3}\pi^2$$ where
$\Psi^{(1)}(k)$ is the polygamma function.
Is it known any sharper bound for the sum $S$?
Thanks.

Solutions Collecting From Web of "Polygamma function series: $\displaystyle\sum_{k=1}^{\infty }\left(\Psi^{(1)}(k)\right)^2$"

The upper bound can be improved using asymptofic series :

enter image description here

First, we have
$$
\begin{align}
\psi'(n)
&=\sum_{k=0}^\infty\frac1{(k+n)^2}\\
&=\sum_{k=n}^\infty\frac1{k^2}\tag{1}
\end{align}
$$
Then
$$
\begin{align}
\sum_{n=1}^\infty\psi'(n)^2
&=\sum_{n=1}^\infty\sum_{j=n}^\infty\frac1{j^2}\sum_{k=n}^\infty\frac1{k^2}\tag{2}\\
&=\sum_{n=1}^\infty\left(\sum_{j=n}^\infty\frac1{j^4}+2\sum_{j=n}^\infty\sum_{m=1}^\infty\frac1{j^2}\frac1{(j+m)^2}\right)\tag{3}\\
&=\sum_{j=1}^\infty\sum_{n=1}^j\frac1{j^4}+2\sum_{j=1}^\infty\sum_{m=1}^\infty\sum_{n=1}^j\frac1{j^2}\frac1{(j+m)^2}\tag{4}\\
&=\sum_{j=1}^\infty\frac1{j^3}+2\sum_{j=1}^\infty\sum_{m=1}^\infty\frac1{j(j+m)^2}\tag{5}\\
&=\zeta(3)+2\sum_{n=1}^\infty\frac{H_{n-1}}{n^2}\tag{6}\\
&=\zeta(3)-2\zeta(3)+2\sum_{n=1}^\infty\frac{H_n}{n^2}\tag{7}\\[9pt]
&=\zeta(3)-2\zeta(3)+4\zeta(3)\tag{8}\\[18pt]
&=3\zeta(3)\tag{9}
\end{align}
$$
Explanation:
$(2)$: use $(1)$
$(3)$: first sum covers $j=k$ the other $j\lt k$ and $j\gt k$
$(4)$: change order of summation
$(5)$: sum in $n$
$(6)$: $n=j+m$ and $j=n-m$ and sum in $m$
$(7)$: $H_{n-1}=H_n-\frac1n$
$(8)$: equation $(14)$ of this answer with $q=2$
$(9)$: add

You may evaluate your series in closed form.

Here are the steps.

Recall the following standard representation for the digamma function:
$$
\psi(x) = -\gamma+\int_0^1 \frac{1 – t^{x-1}}{1 – t}{\rm{d}} x, \quad x>0,
$$
giving, by differentiation,
$$
\psi'(x) = -\int_0^1 \frac{t^{x-1} \ln t}{1 – t}{\rm{d}} x, \quad x>0.
$$
One may deduce
$$
\left(\psi'(k)\right)^2 = \int_0^1\int_0^1 \frac{(uv)^{k-1} \ln u\ln v}{(1 – u)(1-v)}{\rm{d}}u \:{\rm{d}} v
$$
and
$$
\sum_{k=1}^\infty\left(\psi'(k)\right)^2 = \int_0^1 \! \!\int_0^1 \frac{\ln u\ln v}{(1-uv)(1 – u)(1-v)}{\rm{d}} u\:{\rm{d}} v
$$
By partial fraction decomposition and successive integrations involving $\text{Li}_{2}(\cdot)$, you arrive at
$$
\sum_{k=1}^\infty\left(\psi'(k)\right)^2 = 3\zeta(3) = 3.60617070947878285619921 \cdots.
$$

Continuing from Olivier Oloa’s answer,

$$ \begin{align} \sum_{k=1}^{\infty} \big(\psi^{(1)} (k) \big)^{2} &= \int_{0}^{1} \int_{0}^{1} \frac{\ln u \ln v}{(1-uv)(1-u)(1-v)} \ du \ dv \\ &= \int_{0}^{1} \frac{\ln v}{(1-v)^{2}} \int_{0}^{1} \left(\frac{\ln u}{1-u} – \frac{v \ln u}{1-vu} \right) \ du \ dv \\ &= \int_{0}^{1} \frac{\ln v}{(1-v)^{2}} \left( \int_{0}^{1} \frac{\ln u}{1-u} \ du – v \int_{0}^{1} \frac{\ln u}{1-vu} \ du\right) \ dv \end{align}$$

where $$ \int_{0}^{1} \frac{\ln u}{1-u} \ du = \int_{0}^{1} \frac{\ln (1-w)}{w} = -\text{Li}_{2}(1) = -\zeta(2) $$

and

$$ \int_{0}^{1} \frac{\ln u}{1-vu} \ du = – \frac{1}{v} \ln(1-vu) \ln u \Bigg|^{1}_{0} + \frac{1}{v} \int_{0}^{1} \frac{\ln (1-vu)}{u} \ du = – \frac{\text{Li}_{2}(v)}{v} .$$

Therefore,

$$ \sum_{k=1}^{\infty} \big(\psi^{(1)} (k) \big)^{2} = \int_{0}^{1} \frac{\ln v}{(1-v)^2} \Big(\text{Li}_{2}(v) – \zeta(2) \Big) \ dv .$$

Then integrating by parts,

$$ \begin{align} &= \big(\text{Li}_{2}(v) – \zeta(2) \big) \left(\ln (1-v) + \frac{v \ln v}{1-v} \right)\Bigg|^{1}_{0} + \int_{0}^{1} \frac{\ln^{2}(1-v)}{v} \ dv + \int_{0}^{1} \frac{\ln(1-v) \ln v}{1-v} \ dv \\ &= \int_{0}^{1} \frac{\ln^{2}(1-v)}{v} \ dv + \int_{0}^{1} \frac{ \ln(1-v) \ln v}{1-v} \ dv \end{align}$$

where $$ \begin{align} \int_{0}^{1} \frac{\ln^{2}(1-v)}{v} \ dt &= \ln^{2}(1-v)\ln v \Bigg|^{1}_{0} + 2 \int_{0}^{1} \frac{\ln(1-v) \ln v}{1-v} \ dv \\ &= 2 \int_{0}^{1} \frac{ \ln(1-v) \ln v}{1-v} \ dv . \end{align} $$

Therefore, $$ \begin{align} \sum_{k=1}^{\infty} \big(\psi^{(1)} (k) \big)^{2} &= 3 \int_{0}^{1} \frac{\ln(1-v) \ln v}{1-v} \ dv \\ &= -3 \int_{0}^{1} \ln v \sum_{n=1}^{\infty} H_{n}v^{n} \ dv \\ & = -3 \sum_{n=1}^{\infty} H_{n} \int_{0}^{1} v^{n} \ln v \ dv \\ &= 3 \sum_{n=1}^{\infty} \frac{H_{n}}{(n+1)^{2}} \\ &= 3 \left(\sum_{n=1}^{\infty} \frac{H_{n+1}}{(n+1)^{2}} – \sum_{n=1}^{\infty} \frac{1}{(n+1)^{3}} \right) \\ &= 3 \left( \sum_{n=1}^{\infty} \frac{H_{n}}{n^{2}} -1 – \zeta(3) + 1 \right) \\ &= 3 \big(2 \zeta(3) – \zeta(3) \big) \tag{1} \\ &= 3 \zeta(3) .\end{align} $$

$ $

$(1)$ Generalized Euler sum $\sum_{n=1}^\infty \frac{H_n}{n^q}$