Question: Is there a way to find the complete asymptotic expansion (by deriving the correct $O(1)$ term) via the Euler-Maclaurin formula?
Let $$s_n=\sum_{j=1}^n j\log j$$
Fix $m\ge 1$. By the Euler-Maclaurin formula
$$
\begin{align*}
s_n&=\int_1^n t\log t\ dt+\frac{n\log n}{2}+\frac{\log n}{12} -\sum_{s=2}^{m-1} \frac{B_{2s}}{(2s)!}\left(\frac{(2s-3)!}{n^{2s-2}}-\frac{(2s-3)!}{1^{2s-2}}\right)+ R_m \\
&= \frac{n^2 \log n}{2}-\frac{n^2}{4}+\frac{n\log n}{2}+\frac{\log n}{12}
+\frac{1}{4}-\sum _{s=2}^{m-1} \frac{B_{2s}}{(2s)(2s-1)(2s-2)}\left(\frac{1}{n^{2s-2}}-1\right) +R_m(n)
\end{align*}$$
where
$$R_m(n)= \int _1^n\frac{B_{2m}-\tilde B_{2m}(t)}{(2m)!}\frac{(2m-2)!}{t^{2m-1}}\ dt=\int_1^n \frac{B_{2m}-\tilde B_{2m}(t)}{(2m)(2m-1)t^{2m-1}}\ dt$$
Now our goal is to bound $R_m(n)$ for large $n$. Since $|B_{2m}-\tilde B_{2m}(t)|\le 2|B_{2m}|$ we have by the triangle inequality for integrals
$$\begin{align*}
|R_m (n)| &\le \int _1^{n} \frac{2|B_{2m}|}{(2m)(2m-1)t^{2m-1}}\ dt \\&=\frac{2|B_{2m}|}{(2m)(2m-1)(2m-2)}\left(1-\frac{1}{n^{2m-2}}\right)
\end{align*}$$
so bounded for $m\ge 1$ and in particular, $R_m(n)=O(1)$ for large $n$. So it should be that
$$s_n\overset{?}{\sim} \frac{n^2\log n}{2}-\frac{n^2}{4}+\frac{n\log n}{2}+\frac{1}{4} -\sum_{s=2}^{m-1} \frac{B_{2s}}{(2s)(2s-1)(2s-2) n^{2s-2}}$$
for $m\ge 1$ as $n\to +\infty$. However, I don't think this is a valid asymptotic expansion since if $m=1$ the $R_m(n)$ term contributes a factor of $O(1)$; so the $1/4$ term is off by a constant factor. After doing some searching, indeed, the correct constant is $C:=\log A =0.248754\dots$ where $A$ is the Glaisher–Kinkelin constant. This can be found via the Abel–Plana formula, which is a separate topic.
When $m=1$
$$R_1(n)=\frac{1}{2}\int_1^n \frac{B_{2}-\tilde B_{2}(t)}{t}\ dt=\frac{\log n}{12} -\int _1^n \frac{\tilde B_2(t)}{t}\ dt$$
after substituting $B_2=1/6$. I do not know how to evaluate second integral.
Best Answer
This is too long for a comment.
I could have asked the same question years ago but, using Mathematica, I found a way to do it for any order
Series[Unevaluated[Sum[k Log[k], {k, 1, n}]], {n, Infinity, 7}]generates before any simplifcation
$$\log \left(e^{-\frac{n^2}{4}} n^{\frac{n^2}{2}+\frac{n}{2}+\frac{1 }{12}} \left(A+\frac{A}{720 n^2}-\frac{1433 A}{7257600 n^4}+\frac{1550887 A}{15676416000 n^6}+O\left(\frac{1}{n^8}\right)\right)\right)$$ which is $$\frac{6 n^2+6 n+1}{12}\, \log (n)-\frac{n^2}{4}+\log (A)+\frac{1}{720 n^2}-\frac{1}{5040n^4}+\frac{1}{10080n^6}+O\left(\frac{1}{n^8}\right)$$