Asymptotic expansion of $\small\sum_{k=1}^n H_k^{-1}$, where $\small H_k$ are harmonic numbers I'm interested in an asymptotic expansion of $\,\sum_{k=1}^n H_k^{-1}$ for $\,n\to\infty$, where $H_k=\sum_{m=1}^km^{-1}$ are harmonic numbers. I would like all terms in the expansion to use $H_n$ rather than $\log n$ to avoid irrational numbers where possible.
So far, I have figured the two initial terms by replacing the harmonic numbers with their asymptotic expansion $H_n\sim\log n+\gamma+O\!\left(n^{-1}\right)$ and approximating the sum by the corresponding integral. If I haven't made any mistakes, then
$$\sum_{k=1}^n H_k^{-1}\sim n\,H_n^{-1}+n\,H_n^{-2}+O\!\left(n\,H_n^{-3}\right).$$
But computations become more complicated after that, so I couldn't make any further progress. Could you propose an approach that allows to compute more terms, or find a general formula for them?
 A: Denote your sum by $S(n)$. In terms of the digamma function, $H_n=\psi(n+1)+\gamma$. Thus, by the Euler–Maclaurin formula
$$
S(n) = \int_1^n {\frac{{\mathrm{d}t}}{{\psi (t + 1) + \gamma }}}  + \mathcal{O}(1).
$$
Taking $s = \psi (t + 1)$ gives
$$
S(n) = \int_{\psi (2)}^{\psi (n + 1)} {\frac{1}{{s + \gamma }}\frac{{\mathrm{d}t}}{{\mathrm{d}s}}\mathrm{d}s}  + \mathcal{O}(1).
$$
It is known that
$$
t =  - \frac{1}{2} + \mathrm{e}^s  + \mathcal{O}(\mathrm{e}^{ - s} ),
$$
whence
\begin{align*}
S(n) & = \int_{\psi (2)}^{\psi (n + 1)} {\frac{{\mathrm{e}^s }}{{s + \gamma }}\mathrm{d}s}  + \mathcal{O}(1)\int_{\psi (2)}^{\psi (n + 1)} {\frac{{\mathrm{e}^{ - s} }}{{s + \gamma }}\mathrm{d}s}  + \mathcal{O}(1) \\ & = \int_{\psi (2)}^{\psi (n + 1)} {\frac{{\mathrm{e}^s }}{{s + \gamma }}\mathrm{d}s}  + \mathcal{O}(1) = \mathrm{e}^{ - \gamma } \int_{\psi (2) + \gamma }^{H_n } {\frac{{\mathrm{e}^v }}{v}\mathrm{d}v}  + \mathcal{O}(1) = \mathrm{e}^{ - \gamma } {\mathop{\rm Ei}\nolimits} (H_n ) + \mathcal{O}(1).
\end{align*}
Here ${\mathop{\rm Ei}\nolimits}$ is the exponential integral function. Using the known asymptotic expansion of this function, we find
\begin{align*}
S(n) & \sim \mathrm{e}^{H_n  - \gamma } \sum\limits_{k = 0}^\infty  {\frac{{k!}}{{H_n^{k + 1} }}}  + \mathcal{O}(1)= \mathrm{e}^{\log n + \mathcal{O}(1/n)} \sum\limits_{k = 0}^\infty  {\frac{{k!}}{{H_n^{k + 1} }}} + \mathcal{O}(1)\\ & = n\sum\limits_{k = 0}^\infty  {\frac{{k!}}{{H_n^{k + 1} }}}  + \mathcal{O}(1)\sum\limits_{k = 0}^\infty  {\frac{{k!}}{{H_n^{k + 1} }}} + \mathcal{O}(1)\sim n\sum\limits_{k = 0}^\infty  {\frac{{k!}}{{H_n^{k + 1} }}} ,
\end{align*}
i.e.,
$$\boxed{
S(n) \sim n\sum\limits_{k = 0}^\infty  {\frac{{k!}}{{H_n^{k + 1} }}} }
$$
as $n\to +\infty$.
