I would like to obtain a closed form for the following limit: $$I_2=\lim_{k\to \infty} \left ( - (\ln k)^2/2 +\sum_{n=1}^{k} \psi(n) \, \psi'(n) \right)$$
Here $\psi(n)$ is digamma function.
Using the method detailed in this answer, I was able to compute simpler, related series: $$\lim_{k\to \infty} \sum_{n=1}^{k} \left (\psi'(n) -1/n \right) =1 $$ $$\sum_{n=1}^{\infty} \psi''(n) =-\frac{\pi^2}{3} $$ But $I_2$ seems to be tougher because of the product of two digamma's. The divergence of $(\ln k)^2$ is matched by the first terms of asymptotic series for $\psi(n) \psi'(n)$, via the definition of the Stieltjes number, $$ \lim_{k\to \infty} \left ( \sum_{n=1}^{k} \frac{\ln n}{n} - (\ln k)^2/2 \right ) =\gamma_1 $$ but I am stuck with the reminder term.
Side remark: the problem originates in physics, see my older question. In particular, I was able to show that $\langle x \rangle \approx -0.251022$ defined in that question actually equals exactly $-(1+\gamma_0)/(2 \pi)$ where $\gamma_0$ is Euler's constant. This answer I seek here is the only piece missing on my path to closed form $\langle x^2 \rangle$.