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

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. $$

The upper bound can be improved using asymptofic series :