<s>[INST] A closed form of the series $\sum_{n=1}^{\infty} \frac{H_n^2-(\gamma + \ln n)^2}{n}$ I have found a closed form of the following new series involving non-linear harmonic numbers.     

Proposition. $$\sum_{n=1}^{\infty} \dfrac{H_n^2-(\gamma + \ln n)^2}{n} = \dfrac{5}{3}\zeta(3)-\dfrac{2}{3}\gamma^3-2\gamma \gamma_{1}-\gamma_{2} $$ 

where      \begin{align} & H_{n}: =\sum_{k=1}^{n}\frac{1}{k} \\ &\gamma: =\lim_{n\to\infty} \left(\sum_{k=1}^n \frac{1}{k}-\ln n\right) \\ & \gamma_{1}:=\lim_{n\to\infty} \left(\sum_{k=1}^n \frac{\ln k}{k}-\frac{1}{2}\ln^2 n\right)\\& \gamma_{2}: =\lim_{n\to\infty} \left(\sum_{k=1}^n \frac{\ln^2 k}{k}-\frac{1}{3}\ln^3 n\right), \end{align} 
$\gamma_1, \gamma_2$ being Stieltjes constants.
What tool would use to prove it?
 [/INST] 
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.
</s>