Consider the series
$$\mathcal{S} = \sum_{n=1}^{\infty} \frac{\left ( \psi^{(1)}(n) \right )^2}{n^2}$$
One way to evaluate it is by invoking the identity $\displaystyle \zeta(2) - \mathcal{H}_{n-1}^{(2)} = \psi^{(1)}(n)$ where $\mathcal{H}_{n}^{(2)}$ is the harmonic number of order $2$. Invoke some Euler sum identities and voila .
However, in this link we have seen the evaluation of the series $\sum \limits_{n=1}^{\infty} \left ( \psi^{(1)}(n) \right )^2$. I wonder if we can use the same technique used in that link to tackle this series. This is the main question.
Background: This particular series arose while trying to evaluate the polygamma series
$$\mathcal{S} = \sum_{n=1}^{\infty} \frac{\psi^{(1)}(n) \psi^{(1)} (n+1)}{n^2} $$
I managed to solve the latter series by not using this one but nevertheless it's quite interesting of its own. If anyone is interested the series evaluates to $\frac{9 \zeta(6)}{8}$.
