Is there a similar formula like the Stirling one on the sum over $\ln(n)$ (take logarithms on its factorial representation), $$ \sum_{n=1}^N \ln(n) = N\cdot \ln(N)−N+\frac{\ln(N)}{2}+\frac{\ln(2π)}{2}+O\left(\frac{\ln(N)}{N}\right), $$ but for the following sum?
$$\sum_{n=1}^N H_n \cdot \ln(n)$$
I already advanced on getting good approximation on asymptotics using Euler-McLaurin approximation on summation terms till $O(\ln(N))$ order. But further advance is becoming hard for me in $O(1)$ term. This is often referred as the Ramanujan summation of this series.
The constant term is difficult to construct. I'll present the terms to $O(1/n)$ which includes the constant. This problem was asked in MSE 2891159, in which case I answered one of three; this problem is that problem's $\sigma_c(n).$ The problem can be setup exactly as that previous problem, by differentiating with respect to $s$ the following, and taking limits when needed: $$ \sum_{k=1}^n k^{-s}=\zeta{(s)}-\frac{1}{(s-1)n^{s-1}}+\frac{1}{2n^s}- \frac{1}{12}\frac{s}{n^{s+1}} +.... $$ Higher terms are neglected because what is present is sufficient to get to $O(1/n)$ terms. The formula follows from the Euler-McLaurin summation. As in the other problem, break the sum you want into primitives consisting of sums over $\log{k}/k^m$ or $\log^2{k}/k^m.$ That is, $$\sum_{k=1}^n H_k \log{k} = \sum_{k=1}^n \Big(\gamma + \log{k} + \frac{1}{2k} -\frac{1}{12k^2} \Big) \log{k}$$ where the asymptotic formula for $H_k = \gamma + \psi(k+1)$ has been used and a sufficient number of terms have been taken, $except$ for those need to derive the constant term (we'll come back to that later). Letting $L=\log{n}$ we have the following: $$ v_0=\sum_{k=1}^n \log^2{k}=n\big(L^2+2(1-L)) + \frac{L^2}{2}+\frac{L}{6n}+\frac{\gamma^2}{2}-\frac{\pi^2}{24}-\frac{1}{2}\log^2{(2\pi)}+\gamma_1$$ $$ \quad v_1=\sum_{k=1}^n \log{k}=n(L-1)+\frac{L}{2}-\zeta'(0)+\frac{1}{12n}$$ $$v_2= \sum_{k=1}^n \frac{\log{k}}{k} =\frac{L^2}{2}+\gamma_1+\frac{L}{2n}$$ $$ v_3=\sum_{k=1}^n \frac{\log{k}}{k^2} =-\zeta'(2)-\frac{L+1}{n}$$ Let $\tilde{v_k}=(v_k$ with constant term set to 0). Then $$\sigma_c(n)=\sum_{k=1}^n H_k \log{k} =\tilde{v_0}+\gamma\, \tilde{v_1} + \frac{1}{2} \tilde{v_2}-\frac{1}{12} \tilde{v_3} + C $$ $$=n\big(L^2+(\gamma-2)(L-1)\big) + \frac{3L^2}{4}+\gamma\frac{L}{2}+\frac{1}{2n}(L+\frac{\gamma+1}{6}) + C $$ where $C$ is the unknown constant term. A correct way to determine C is simply by $$ C=\lim_{n \to \infty} \Big(\sum_{k=1}^n H_k \log{k} - \Big(n\big(L^2+(\gamma-2)(L-1)\big) + \frac{3L^2}{4}+\gamma\frac{L}{2} \Big)\, \Big)$$ Proceed further only if you are comfortable with formal mathematics. Now I'm going to find $C$ in a different manner. It is well-known that $$ \gamma = \lim_{n \to \infty}\Big( \sum_{k=1}^n\frac{1}{k} - \log{n} \Big) = \int_0^\infty \Big(\frac{1}{e^x-1} - \frac{e^{-x}}{x} \Big)dx.$$ The first expression is similar in spirit to how we defined $C,$ but the second equation is useful in that many digits of $\gamma$ can be extracted quite easily with your favorite numerical integration routines.
In derving $v_0$ through $v_3$ I stopped because this is sufficient to get terms to $O(1/n).$ If we kept going, there would be a $v_5 = \sum_{k=1}^n \log{k}/k^4= -\zeta'(4) + o(1/n)$, etc. The natural thing to do is add all the $-\zeta'(2n)$ terms up with the appropriate weights to get a new constant that, when added to the ones in the $v_k$, constitute $C.$ This was done in MSE 2891159 and very fortunately the series converged. No such luck here. The weights are from the asymptotic expansion of the harmonic numbers, $$H_k = \gamma + \log{k} + \frac{1}{2k}-\sum_{m=1}^\infty \frac{B_{2m}}{2m} k^{-2m} $$ so the constant we want to assign meaning to is $$ \kappa \,\,\dot{=} \sum_{m=1}^\infty \frac{B_{2m}}{2m} \zeta'(2m) $$ where the dotted equals means 'representation' instead of 'equals.' The idea is akin to saying $$\sum_{k=1}^\infty k \,\, \dot{=} -1/12$$ which can be given rigor in terms of zeta regularization. First differentiate the well-known Euler-like integral for the zeta function to find $$\zeta'(2s) =\frac{1}{\Gamma(2s)} \int_0^\infty \frac{t^{2s-1}}{e^t-1} \Big( \log{t} - \psi(2s) \Big) .$$ Insert into the definition for $\kappa$ and interchange $\int$ and $\sum$
$$ (K) \quad \kappa \,\,\dot{=} \int_0^\infty \frac{dt/t}{e^t-1}\sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \big( \log{t} - \psi(2m) \big). $$ The inner sum converges (its been Borel transformed), we just need to find an expression that's not a power series. We need some formulas. It is well-known that, and it will be used repeatedly, $$\sum_{m=1}^\infty \frac{B_m}{m!} t^m = \frac{t}{e^t-1} -1, \quad \sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m}= \frac{t}{e^t-1} -1+t/2 $$ An integration of the previous equation leads to $$\sum_{m=1}^\infty\frac{B_{2m}}{(2m)(2m)!} t^{2m}= -\log\big(\frac{t}{e^t-1}\big) -t/2 $$ Let's begin on the 'psi' term: $$ \Psi(t):=\sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \,\psi(2m)= \sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \,\big(\psi(2m)+\frac{1}{2m} - \frac{1}{2m} \big)=$$ $$= \sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \,\psi(2m+1)+\Big(t/2+\log{\big(\frac{t}{e^t-1}\big)} \Big) =$$ $$= \sum_{m=1}^\infty \frac{B_{m}}{m!} t^m \,\psi(m+1)+ \Big( \frac{1}{2}(1-\gamma)t \Big) + \Big(t/2+\log{\big(\frac{t}{e^t-1}\big)} \Big) .$$ Again use $H_m = \gamma + \psi(m+1)$ to find $$\Psi(t)=\sum_{m=1}^\infty \frac{B_m}{m!} t^m \,H_m - \gamma\Big(\frac{t}{e^t-1} - 1\Big)+ \big(1-\frac{\gamma}{2}\big)t + \log{\big(\frac{t}{e^t-1}\big) }$$ The harmonic number is used because there is the integral relation $$H_m = -m \int_0^1 dx x^{m-1} \log{(1-x)} $$ Insert this, switch $\int$ and $\sum$, sum up the series in closed form, and $finally$ $$\Psi(t)=-\int_0^t \frac{\log{(1-u/t)}}{e^u-1}\Big(1-\frac{u\,e^u}{e^u-1}\Big)\,du- \gamma\Big(\frac{t}{e^t-1} - 1\Big)+ \big(1-\frac{\gamma}{2}\big)t + \log{\big(\frac{t}{e^t-1}\big) }$$ In eq. $(K)$ the sum before the $\log{t}$ is one of the known formulas. Therefore a have a (double) integral relationship that is perfectly well-behave for $\kappa$ $$\kappa=\int_0^\infty \frac{dt/t}{e^t-1}\Big(\, \log{t}\big(\frac{t}{e^t-1}-1+t/2\big) - \Psi(t)\Big) = -0.077596...$$ When using $\kappa$ the answer can be stated as $$ C=\kappa + \frac{\gamma^2}{2}-\frac{\pi^2}{24}+\frac{1}{2}\gamma \, \log(2\pi) - \frac{1}{2} \log^2{(2\pi)} + \frac{3}{2} \gamma_1 $$ Six digits of agreement are obtained by using the asymptotic formula and comparing to the brute force summation for $n=$ 20.