Series expansion of $\sum_{k=1}^{n} \frac{\ln(k+1)}{k}$ to order $o(1)$ I am interested in the asymptotic development of $a_n\stackrel{\rm def}{=} \sum_{k=1}^{n} \frac{\ln(k+1)}{k}$, when $n\to \infty$. Via a standard comparison series/integral, I got that
$$
\int_1^{n+1} dx\,\frac{\ln(x+1)}{x} \leq a_n \leq \int_0^{n} dx\,\frac{\ln(x+1)}{x}
$$
which, conditioning on my not messing up the computations, leads to
$$a_n = \frac{1}{2}{\ln^2 n} + O(1)
$$
where the $O(1)$ term is, ignoring lower-order terms, between $\frac{\pi^2}{6}$ and $\frac{\pi^2}{12}$.
For my application, however, it is crucial that I get $a_n$ to order $o(1)$ (and the higher the better: $o(\frac{1}{n})$, for instance, would be great). I tried a bit this morning, but couldn't figure out out to do so: by starting the comparison at a higher indice, I can reduce the range of uncertainty about the constant term, but didn't get to completely get it.
 A: Let $f(x) = \dfrac{\ln (x+1)}{x}$. Taking a first step towards the Euler-Maclaurin formula, we have
$$\sum_{k = 1}^n f(k) = \frac{f(1) + f(n)}{2} + \int_1^n f(x)\,dx + \int_1^n p_1(x) f'(x)\,dx,$$
where $p_1(x) = \lbrace x\rbrace - \frac{1}{2}$. Generally, I let $p_n(x) = \frac{1}{n!} \mathscr{B}_n(\lbrace x\rbrace)$, where $\mathscr{B}_n$ is the $n^{\text{th}}$ Bernoulli polynomial. Evaluating the first integral to any desired precision is rather easy:
\begin{align}
\int_1^n \frac{\ln (x+1)}{x}\,dx
&= \int_1^n \frac{\ln x}{x}\,dx + \int_1^n \frac{\ln (1 + x^{-1})}{x}\,dx \\
&= \frac{1}{2}(\ln n)^2 + \int_{1/n}^1 \frac{\ln (1+u)}{u}\,du \tag{$u = x^{-1}$}\\
&= \frac{1}{2}(\ln n)^2 + \sum_{k = 1}^{\infty} \frac{(-1)^{k-1}}{k^2}(1 - n^{-k}) \\
&= \frac{1}{2}(\ln n)^2 + \frac{\pi^2}{12} + \sum_{k = 1}^{\infty} \frac{(-1)^k}{k^2n^k}.
\end{align}
Now we can decide whether
$$K = \int_1^{\infty} p_1(x) f'(x)\,dx$$
can be evaluated with the desired accuracy easily enough to stop there. If so, we have
$$\sum_{k = 1}^n f(k) = \frac{1}{2}(\ln n)^2 + \biggl(\frac{\pi^2}{12} + \frac{\ln 2}{2} + K\biggr) + \frac{\ln(n+1)}{2n} + \sum_{k = 1}^{\infty} \frac{(-1)^k}{k^2n^k} - \int_n^{\infty} p_1(x)f'(x)\,dx,$$
and due to the cancellations because of the sign-changes of $p_1$ we have
$$\int_n^{\infty} p_1(x)f'(x)\,dx = \int_n^{\infty} p_1(x) \biggl(\frac{1}{x(x+1)} - \frac{\ln(x+1)}{x^2}\biggr)\,dx \in O\biggl(\frac{\ln n}{n^2}\biggr),$$
thus
$$\sum_{k = 1}^n \frac{\ln (k+1)}{k} = \frac{1}{2}(\ln n)^2 + \biggl(\frac{\pi^2}{12} + \frac{\ln 2}{2} + K\biggr) + \frac{\ln n}{2n} - \frac{1}{n} + O\biggl(\frac{\ln n}{n^2}\biggr).\tag{$\ast$}$$
If we don't consider the integral defining $K$ to be nice enough, we can integrate by parts to find
$$\int_1^n p_1(x)f'(x)\,dx = p_2(x)f'(x)\biggr\rvert_1^n - \int_1^n p_2(x)f''(x)\,dx = \frac{1}{12}\bigl(f'(n) - f'(1)\bigr) - \int_1^n p_2(x)f''(x)\,dx.$$
Then we face the same decision as above, do we consider
$$\int_1^{\infty} p_2(x)f''(x)\,dx$$
easy enough to evaluate? If not, we can once more integrate by parts. Since $p_{2m+1}(1) = p_{2m+1}(n) = 0$ for $m \geqslant 1$, this time we don't pick up boundary terms and find
$$-\int_1^n p_2(x)f''(x)\,dx = \int_1^n p_3(x)f'''(x)\,dx.$$
One can continue integrating by parts, but I doubt whether the faster decay of the higher-order derivatives of $f$ makes the numerical evaluation of the integral
$$\int_1^{\infty} p_n(x) f^{(n)}(x)\,dx$$
easier for long, since the $p_n$ become more complicated.
