Sum with squared trigamma function: $\sum_{n=1}^{\infty}\left(n\psi_{1}^{2}(n)-\frac{1}{n}\right)$ Here is a question I found on web
$$S=\sum_{n=1}^{\infty}\left(n\psi_{1}^{2}(n)-\frac{1}{n}\right)$$
The closed form is: $\frac{3}{2}\zeta(3)+\frac{3}{2}-\frac{\pi^2}{12}$. I know how to handle $\psi_{1}^{2}(n)$ by using double integrals and integral representation of the trigamma function. But with the term $\frac{1}{n}$ I don't know how to get rid of it.
Clearly, the sum of $n\psi_1^2(n)$ is divergent, so I doubt there is a way to eleminate the term $\frac{1}{n}$ but I can't find it yet.
Any help is welcome, and I am grateful for that, thank you so much.
 A: Perhaps this has something to do with the theory of multiple zeta values. Here is a self-contained solution using basic properties of $\psi_1(n)=\sum_{m\geqslant n}1/m^2$, namely $0<\psi_1(n)-1/n<1/n^2$: $$\small\psi_1(n)-\frac1n=\sum_{m\geqslant n}\left(\frac1{m^2}-\frac1{m(m+1)}\right)=\sum_{m\geqslant n}\frac1{m^2(m+1)}<\sum_{m\geqslant n}\frac1{nm(m+1)}=\frac1{n^2}$$ (which implies $\lim\limits_{n\to\infty}n\psi_1(n)=1$ and the convergence of $S$) and $\sum_{n\geqslant 1}\psi_1(n)/n=2\zeta(3)$, equivalent to the simplest Euler sum a.k.a. $\zeta(2,1)$; an easy proof is as follows:
$$\small\Xi=\sum_{n\geqslant 1}\frac1n\sum_{m\geqslant n}\frac1{m^2}=\sum_{m\geqslant 1}\frac1{m^2}\sum_{1\leqslant n\leqslant m}\frac1n=\sum_{m\geqslant 1}\frac1{m^2}\sum_{n\geqslant 1}\left(\frac1n-\frac1{n+m}\right)=\sum_{n,m\geqslant 1}\frac1{nm(n+m)}\\\small=\sum_{n,m\geqslant 1}\left(\frac1n+\frac1m\right)\frac1{(n+m)^2}=2\sum_{n,m\geqslant 1}\frac1{n(n+m)^2}=2\sum_{n\geqslant 1}\frac1n\sum_{m>n}\frac1{m^2}=2\big(\Xi-\zeta(3)\big).$$

We compute (in two ways) the limit $\lim\limits_{N\to\infty}(S_N-H_N)$, where $$\color{gray}{H_N=\sum_{n=1}^N\frac1n,}\qquad S_N=\sum_{n=1}^N n\left(\sum_{m=n}^N\frac1{m^2}\right)^2.$$
First, writing $S_{N-1}=\sum_{n=1}^{N-1}n\big(\psi_1(n)-\psi_1(N)\big)^2$ and expanding the square, we get $$S_{N-1}-H_{N-1}=\sum_{n=1}^{N-1}\left(n\psi_1^2(n)-\frac1n\right)-2\psi_1(N)\sum_{n=1}^{N-1}n\psi_1(n)+\psi_1^2(N)\frac{N(N-1)}{2}.$$
Now $\lim\limits_{n\to\infty}n\psi_1(n)=1$ implies $\lim\limits_{N\to\infty}(1/N)\sum_{n=1}^{N-1}n\psi_1(n)=1$ by Stolz–Cesàro, hence $$\lim_{N\to\infty}(S_N-H_N)=\color{blue}{S-\frac32}.$$
Second, fix $N$ and put $A_n=n(n+1)/2$ and $B_n=\sum_{m=n}^N 1/m^2$ (with $B_{N+1}=0$), then $$S_N=\sum_{n=1}^N(A_n-A_{n-1})B_n^2=\sum_{n=1}^N A_n(B_n^2-B_{n+1}^2)\\\small=\sum_{n=1}^N A_n(B_n-B_{n+1})(B_n+B_{n+1})=\sum_{n=1}^N\frac{n+1}{2n}\left(2B_n-\frac1{n^2}\right)$$ and, since $\sum_{n=1}^N B_n=H_N$, and $B_n\to\psi_1(n)$ monotonically (as $N\to\infty$), we get $$\lim_{N\to\infty}(S_N-H_N)=\sum_{n=1}^\infty\left(\frac{\psi_1(n)}{n}-\frac{n+1}{2n^3}\right)=\color{blue}{\frac32\zeta(3)-\frac12\zeta(2)}.$$
