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?

The leading behavior may be understood by rendering
$\Sigma_{k=0}^n k^n = n^n\Sigma_{k=0}^n(k/n)^n = n^n\Sigma_{m=0}^n(1-m/n)^n.$
Here $m=n-k$. Then the terms in the last sum successively approach $e^{-m}$ as $n$ increases without bound, rendering the asymptotic result
$\Sigma_{k=0}^n k^n\text{ ~ }  n^n\Sigma_{m=0}^\infty e^{-m} = (n^n)[e/(e-1)].$
