Asymptotic expansion of $\sum _{k=1}^n \left(\frac{k}{n}\right)^k$

How to calculate the asymptotic expansion of these two sums w.r.t $n$ to arbitrary precision?

$$\sum _{k=1}^n \left(\frac{k}{n}\right)^k,\sum _{k=1}^n \Gamma\left(\frac{k}{n}\right)^{-k}$$

According to O.Furdui's problem book Limits, Series and fractional part integrals it might be a open problem. I have no idea other than knowing the first sum tends to $\frac{e}{e-1}$ as $n\rightarrow \infty$ (Tannery theorem, for instance) and would like you to help. Note that similar sums $\small\sum\limits _{k=1}^n \left(\frac{k}{n}\right)^n$, $\small\sum\limits _{k=1}^n \Gamma\left(\frac{k}{n}\right)^{-n}$ can be approximated using Taylor expansion but not so helpful for this case. Thanks in advance!

My try: Split the sum into three parts $$\sum_{k=1}^n \left(\frac{k}{n}\right)^k = \sum_{1\leq k\leq K} \left(\frac{k}{n}\right)^k + \sum_{K+1 \leq k < n-n^\epsilon} e^{k\log\left(\frac{k}{n}\right)} + \sum_{n-n^\epsilon \leq k \leq n} e^{k\log\left(\frac{k}{n}\right)} $$ for some integer $K$ which defines the order and some small $\epsilon >0$ (say $\epsilon=1/2$). It is easy to see that $k\log\left(\frac{k}{n}\right)$ has a unique minimum at $k=\frac{n}{e}$ somewhere in the range of the middle term for large $n$. Therefore we evaluate the boundary terms of the middle term for some estimate $$k=K+1: \quad \left(\frac{K+1}{n}\right)^{K+1} \\ k=n-n^\epsilon: \quad e^{n(1-n^{\epsilon-1})\log(1-n^{\epsilon-1})} \leq e^{-n^\epsilon + n^{2\epsilon -1}} \, .$$ For fixed $K$ and sufficiently large $n$ the right boundary obviously vanishes exponentially (the optimal $\epsilon$ is $1-\frac{\log 2}{\log n}$ so that $n^\epsilon=n/2$) and the largest value in that range is that for $k=K+1$. Hence the middle term is of order ${\cal O}(n^{-K})$.

For the last term substitute $k\rightarrow k-n$ so that it becomes $$\sum_{0\leq k \leq n^\epsilon} e^{-k +\left[(n-k)\log\left(1-\frac{k}{n}\right) + k\right]} \, .$$ It remains to estimate the square bracket $$(n-k)\log\left(1-\frac{k}{n}\right) + k = -(n-k) \sum_{m=1}^\infty \frac{1}{m}\left(\frac{k}{n}\right)^m + k \\ = \frac{k^2}{n} - (n-k) \sum_{m=2}^\infty \frac{1}{m}\left(\frac{k}{n}\right)^m \\ = \frac{k^2}{2n} + \sum_{m=2}^\infty \frac{k}{m(m+1)} \left(\frac{k}{n}\right)^m = \sum_{m=1}^\infty \frac{k}{m(m+1)} \left(\frac{k}{n}\right)^m$$ which vanishes vor large $n$. For an order $K$ approximation we can thus write $$\sum_{0\leq k \leq n^\epsilon} e^{-k +\left[(n-k)\log\left(1-\frac{k}{n}\right) + k\right]} \\ = \sum_{0\leq k \leq n^\epsilon} e^{-k} \left\{ 1 + \sum_{l=1}^\infty \frac{1}{l!} \sum_{m_1=1}^\infty \cdots \sum_{m_l=1}^\infty \frac{k^{l+m_1+\dots+m_l}}{m_1(m_1+1)\cdots m_l(m_l+1)} \frac{1}{n^{m_1+\dots+m_l}} \right\} \\ = \sum_{0\leq k \leq n^\epsilon} e^{-k} \left\{ 1 + \sum_{p=1}^\infty \frac{k^p}{n^p} \sum_{l=1}^p \frac{k^{l}}{l!} \substack{ \sum_{m_1=1}^\infty \cdots \sum_{m_l=1}^\infty \\ m_1+\dots+m_l \, \stackrel{!}{=} \, p }\frac{1}{m_1(m_1+1)\cdots m_l(m_l+1)} \right\} \, .$$

When evaluating the moments $$ \sum_{0\leq k \leq n^\epsilon} e^{-k} \, k^{p+l} $$ for $p=0,1,2,...,K-1$, the range of summation can be extended up to infinity, because that only introduces an exponentially suppressed error term ${\cal O}\left(n^{(p+l)\epsilon} \, e^{-n^\epsilon}\right)$.

Collecting terms, it is seen that $$\sum_{k=1}^n \left(\frac{k}{n}\right)^k = a_0 + \sum_{k=1}^{K-1} \frac{k^k + a_k}{n^k} + {\cal O}\left(n^{-K}\right) $$ where $$a_0 = \frac{e}{e-1} \\ a_k = \sum_{l=1}^k \frac{\sum_{q=0}^\infty q^{k+l} \, e^{-q}}{l!} \substack{ \sum_{m_1=1}^\infty \cdots \sum_{m_l=1}^\infty \\ m_1+\dots+m_l \, \stackrel{!}{=} \, k }\frac{1}{m_1(m_1+1)\cdots m_l(m_l+1)} \, .$$

For $k\geq 2$ the $a_k$ are extremely close to $k^k$, that is less than $0.04\%$ relative error, so that the total coefficient for $k\geq 2$ is in good approximation $2k^k$.

One term beyond leading order we have for $K=2$ $$\sum_{k=1}^n \left(\frac{k}{n}\right)^k = \frac{e}{e-1} + \frac{1+\frac{e(e+1)}{2(e-1)^3}}{n} + {\cal O}(n^{-2}) \, .$$

Increasing the order $K$ also shifts the range of validity to higher $n$, i.e. it is an asymptotic series. The zero, first and fifth order approximations are shown below. The fifth order is visually not distinguishable from the approximation where $a_k=k^k$ has been used for $k\geq 1$.

Since @Diger's answer captures the main idea, this answer merely amends it to cover the case of $\Gamma$, and provides some computations. First let's restate the result: for $n\to\infty$ $$\sum_{k=1}^{n}(k/n)^k\asymp A_0+\sum_{j=1}^{(\infty)}(j^j+A_j)n^{-j},\qquad A_j=\sum_{k=0}^{\infty}a_j(k),$$ where $a_j(k)$ are the expansion coefficients of $(1-k/n)^{n-k}$ in powers of $1/n$ (for fixed $k$): $$\sum_{j=0}^{\infty}a_j(k)x^j:=(1-kx)^{(1-kx)/x}=\exp\left[-k\left(1-\sum_{j=1}^\infty\frac{(kx)^j}{j(j+1)}\right)\right].$$

Similarly, the contribution to the asymptotics of $$\sum_{k=1}^{n}\big(\Gamma(k/n)\big)^{-k}\asymp\sum_{j=0}^{(\infty)}B_j n^{-j},$$ up to $n^{-j}$, is the one from the first $j$ terms of the defining sum, plus the one from a handful of the last terms, with the "handful" tending to infinity. Explicitly, $B_j=\sum_{k=1}^{j}b_j(k)+\sum_{k=0}^{\infty}c_j(k)$, where $$\big(\Gamma(kx)\big)^{-k}=:\sum_{j=k}^{\infty}b_j(k)x^j,\qquad\big(\Gamma(1-kx)\big)^{-(1-kx)/x}=:\sum_{j=0}^{\infty}c_j(k)x^j.$$ For computations, we use the known expansion $$\log\Gamma(1-x)=\gamma x+\sum_{j=2}^{\infty}\frac{\zeta(j)}{j}x^j$$ from which one deduces $1/\Gamma(x)=\sum_{j=1}^{\infty}g_j x^j$ with $$g_1=1,\quad j g_{j+1}=\gamma g_j-\sum_{k=2}^j(-1)^k\zeta(k)g_{j-k+1}.$$

The first few values of $a_j(k)$ are \begin{align*} a_0(k)&=e^{-k} \\a_1(k)&=\frac{e^{-k}}{2} k^2 \\a_2(k)&=\frac{e^{-k}}{24} (3 k^4 + 4 k^3) \\a_3(k)&=\frac{e^{-k}}{48} (k^6 + 4 k^5 + 4 k^4) \end{align*} The corresponding values of $A_j$ are then \begin{align*} A_0&=\frac{e}{e-1} \\A_1&=\frac{e(e + 1)}{2(e-1)^3} \\A_2&=\frac{e(7 e^3 + 45 e^2 + 21 e - 1)}{24(e-1)^5} \\A_3&=\frac{e(9 e^5 + 193 e^4 + 422 e^3 + 102 e^2 - 7 e + 1)}{48(e-1)^7} \end{align*} Denoting $c:=e^\gamma$, the first three values of $B_j$ are \begin{align*} B_0&=\frac{c}{c-1}, \\B_1&=1-\left(\frac{\pi^2}{12}-\gamma\right)\frac{c(c+1)}{(c-1)^3}, \\B_2&=4+\gamma+\frac{c}{(c-1)^5} \\&\times\left[\left(\frac{\pi^2}{12}-\frac{\zeta(3)}{3}\right)(c^3+3c^2-3c-1)\right. \\&+\left.\frac{1}{2}\left(\frac{\pi^2}{12}-\gamma\right)^2(c^3+11c^2+11c+1)\right]. \end{align*} [The expression for $B_3$ looks too cumbersome to put here.]