A closed form of $\sum_{k=1}^\infty \psi^{(1)} (k+a)\psi^{(1)} (k+b)$? The following result
$$
\sum_{k=1}^\infty\left(\psi^{(1)} (k)\right)^2 = 3\zeta(3)
$$ where $\psi^{(1)}$ is the polygamma function makes me think there is a nice sum for the series

$$
\sum_{k=1}^\infty \left(\psi^{(1)} (k+a)\right)^2 
$$

or

$$
\sum_{k=1}^\infty \psi^{(1)} (k+a)\psi^{(1)} (k+b)
$$

where $a$ and $b$ are any real numbers such that $a >-1, b>-1.$
Could you help me to find it?

We have
$$
\int^{1}_{0}\int^{1}_{0}\frac{\log u\log v}{(1-uv)(1-u)(1-v)}dudv=3\zeta(3)
$$
Also
$$
\int^{1}_{0}\int^{1}_{0}\frac{(uv)^k\log u\log v}{(1-u)(1-v)}dudv=Z(2,k+1)\psi'(k+1)\textrm{, where }Re(k)>-1
$$
The function $Z(s,k)=\sum^{\infty}_{l=0}\frac{1}{(l+k)^s}$ is Hurwitz zeta function.
Hence
$$
\sum^{a-1}_{k=0}\int^{1}_{0}\int^{1}_{0}\frac{(uv)^k\log u\log v}{(1-u)(1-v)}dudv
=\int^{1}_{0}\int^{1}_{0}\frac{(1-(uv)^a)\log u\log v}{(1-uv)(1-u)(1-v)}dudv.
$$
Hence
$$
\sum^{\infty}_{k=1}\left(\psi'(k+a)\right)^2=\int^{1}_{0}\int^{1}_{0}\frac{(uv)^a\log u \log v}{(1-u v)(1-u)(1-v)}dudv=3\zeta(3)-\sum^{a-1}_{k=0}Z(2,k+1)\psi'(k+1)
$$
By trying to generalize the problem we have
$$
C_{\nu}=\underbrace{\int^{1}_{0}\ldots\int^{1}_{0}}_{\nu}\frac{\log u_1\ldots \log u_{\nu}du_1\ldots du_{\nu}}{(1-u_1\ldots u_{\nu})(1-u_1)\ldots(1-u_{\nu})}=\sum^{\infty}_{l=0}\left(\int^{1}_{0}\frac{u^l\log u}{1-u}du\right)^{\nu}=
$$
$$
=\sum^{\infty}_{l=0}\left(-Z(2,l+1)\right)^{\nu}
$$
Also
$$
\underbrace{\int^{1}_{0}\ldots\int^{1}_{0}}_{\nu}\frac{u_1^k\ldots u_{\nu}^k \log u_1\ldots\log u_{\nu}du_1\ldots du_{\nu}}{(1-u_1)\ldots (1-u_{\nu})}=(-1)^{\nu}Z(2,k+1)(\psi'(k+1))^{\nu-1}
$$
Hence
$$
\sum^{a-1}_{k=0}\underbrace{\int^{1}_{0}\ldots\int^{1}_{0}}_{\nu}\frac{u_1^k\ldots u_{\nu}^k \log u_1\ldots\log u_{\nu}du_{1}\ldots du_{\nu}}{(1-u_1)\ldots (1-u_{\nu})}=(-1)^{\nu}\sum^{a-1}_{k=0}Z(2,k+1)(\psi'(k+1))^{\nu-1}=
$$
$$
=\underbrace{\int^{1}_{0}\ldots\int^{1}_{0}}_{\nu}\frac{(1-u_1^a\ldots u_{\nu}^a) \log u_1\ldots\log u_{\nu}du_1\ldots du_{\nu}}{(1-u_1\ldots u_{\nu})(1-u_1)\ldots (1-u_{\nu})}.
$$
Hence
$$
\underbrace{\int^{1}_{0}\ldots\int^{1}_{0}}_{\nu}\frac{(u_1^a\ldots u_{\nu}^a) \log u_1\ldots\log u_{\nu}du_1\ldots du_{\nu}}{(1-u_1\ldots u_{\nu})(1-u_1)\ldots (1-u_{\nu})}=
$$
$$
=C_{\nu}+(-1)^{\nu-1}\sum^{a-1}_{k=0}Z(2,k+1)(\psi'(k+1))^{\nu-1}
$$
And finaly
$$
\sum^{\infty}_{k=1}\left(\psi'(k+a)\right)^{\nu}=C_{\nu}+(-1)^{\nu-1}\sum^{a-1}_{k=0}Z(2,k+1)(\psi'(k+1))^{\nu-1}
$$
