Asymptotic for $\sum_{k=1}^n k^n$ Consider the OEIS sequence A031971, which is defined as:
$$a_n=\sum\limits_{k=1}^n k^n\quad\color{gray}{(1,\,5,\,36,\,354,\,4425,\,67171,\,1200304,\,.\!.\!.\!)}\tag{1}$$ I'm interested in the asymptotic behavior of $a_n$ for $n\to\infty$.
Empirically, it appears that
$$a_n\stackrel{\color{gray}?}\sim\frac{e}{e-1}\,n^n\cdot\left(1-\frac{e+1}{2\,(e-1)^2}\,n^{-1}+c\,n^{-2}+\mathcal O\!\left(n^{-3}\right)\right),\tag{2}$$
where $c\approx0.6310116...$ (I haven't found a plausible closed form it). The leading term $\frac{e}{e-1}\,n^n$ is given in the OEIS.
How can we prove this formula and find higher terms in it?

Let us write
$$
\sum\limits_{k = 1}^n {k^n }  = n^n \sum\limits_{k = 0}^{n - 1} {\left( {1 - \frac{k}{n}} \right)^n } .
$$
Now by power series expansions
\begin{align*}
\left( {1 - \frac{k}{n}} \right)^n & = \exp \left( {n\log \left( {1 - \frac{k}{n}} \right)} \right) = \exp \left( { - n\sum\limits_{j = 1}^\infty  {\frac{{k^j }}{j}\frac{1}{{n^j }}} } \right)
\\ & = e^{ - k}  - e^{ - k} \frac{{k^2 }}{2}\frac{1}{n} + e^{ - k} \left( {\frac{{k^4 }}{8} - \frac{{k^3 }}{3}} \right)\frac{1}{{n^2 }} -  \ldots \,.
\end{align*}
Substituting back to the original sum, extending that sum to $n=\infty$, and using closed forms of $\sum_{k=0}^\infty k^px^k$ (with $x=e^{-1}$), we deduce
\begin{align*}
\sum\limits_{k = 1}^n {k^n } & \sim n^n \left( {\frac{e}{{e - 1}} - \frac{{e^2  + e}}{{2(e - 1)^3 }}\frac{1}{n} - \frac{{e(5e^3  - 9e^2  - 57e - 11)}}{{24(e - 1)^5 }}\frac{1}{{n^2 }} +  \ldots } \right)
\\ & =\frac{e}{{e - 1}}n^n \left( {1 - \frac{{e + 1}}{{2(e - 1)^2 }}\frac{1}{n} - \frac{{5e^3  - 9e^2  - 57e - 11}}{{24(e - 1)^4 }}\frac{1}{{n^2 }} +  \ldots } \right)
\end{align*}
as $n\to +\infty$. You can obtain more terms of the expansion if you like, the process should be clear.
