A closed form of the series $\sum_{n=1}^{\infty} \frac{H_n^2-(\gamma + \ln n)^2}{n}$

Here is a nice solution from R. Tauraso, for those interested in this problem:

Solution proposed by Roberto Tauraso, Roma, Italy.

The Multiple Harmonic Sum is defined by $$H_{n}(s_{1}, \dots, s_{l}) == \sum_{0 < k_{1} < k_{2} < \dots < k_{l} \leq n}\prod_{i = 1}^{l}\frac{1}{k_{i}^{s_{i}}}$$ Then, by the known properties of the stuffle product, \begin{align} H_{n}^{2}(1) &= 2H_{n}(1, 1) + H_{n}(2)\notag\\ H_{n}^{3}(1) &= 6H_{n}(1, 1, 1) + 3H_{n}(1, 2) + 3H_{n}(2, 1) + H_{n}(3)\notag \end{align} Therefore, for $N > 0$, \begin{align} \sum_{n = 1}^{N}\frac{H_{n}^{2}}{n} &= 2\sum_{n = 1}^{N}\frac{H_{n}(1, 1)}{n} + \sum_{n = 1}^{N}\frac{H_{n}(2)}{n}\notag\\ &= 2H_{N}(1, 1, 1) + 2H_{N}(1, 2) + H_{N}(3) + H_{N}(2, 1)\notag\\ &= \frac{1}{3}H_{N}^{3}(1) + \frac{2}{3}H_{N}(3) + H_{N}(1, 2)\notag\\ &= \frac{1}{3}H_{N}^{3}(1) + \frac{2}{3}\zeta(3) + \zeta(2, 1) + o(1)\notag\\ &= \frac{1}{3}H_{N}^{3}(1) + \frac{5}{3}\zeta(3) + o(1)\tag{1} \end{align} because $\zeta(2, 1) = \zeta(3)$ (see for example Thirty-two Goldbach Variations by J. M. Borwein and D. M. Bradley).

Moreover, $H_{N}(1) = \ln N + \gamma + O(1/N)$ implies that \begin{align} \frac{1}{3}H_{N}^{3}(1) - \sum_{n = 1}^{N}\frac{(\gamma + \ln n)^{2}}{n} &= \frac{1}{3}H_{N}^{3}(1) - \gamma^{2}H_{N}(1) - 2\gamma\sum_{n = 1}^{N}\frac{\ln n}{n} - \sum_{n = 1}^{N}\frac{\ln^{2}n}{n}\notag\\ &= \frac{\ln^{3}N}{3} + \frac{\gamma^{3}}{3} + \gamma\ln^{2}N + \gamma^{2}\ln N - \gamma^{2}(\ln N + \gamma)\notag\\ &\,\,\,\,\,\,\,\, -2\gamma\left(\gamma_{1} + \frac{\ln^{2}N}{2}\right) - \left(\gamma_{2} + \frac{\ln^{3}N}{3}\right) + o(1)\notag\\ &= -\frac{2}{3}\gamma^{3} - 2\gamma\gamma_{1} - \gamma_{2} + o(1)\tag{2} \end{align} By adding together $(1)$ and $(2)$, and by taking the limit as $N \to \infty$, we get the result.