Approximation of gamma function via Riemann sums at integer points

Question

I found something curious. We know that the gamma function is defined as $$ \Gamma(n+1) := \int_{t=0}^\infty t^n \exp(-t) dt,$$ and it has the property that $\Gamma(n+1) = n!$ for non-negative integer $n$. I recently came across the following Riemann sum approximation of the gamma function: $$ S(n) := \sum_{k=1}^\infty k^n \exp(-k). $$ A priori there is no reason to suspect that $S(n)$ should approximate $\Gamma(n+1)$ at all, since the points are integer, so they are not getting closer.

However, I plugged in a couple of values of this sum into Wolfram Alpha (e.g., https://www.wolframalpha.com/input?i=sum+exp%28-n%29+n%5E6) and saw that, to my surprise, $S(n)$ has a closed form (some messy formula), and numerically $S(n)$ is pretty close to $n!$, seemingly improving as $n$ gets large.

Is this known? How do we derive the closed form for the sums in general? I tried searching but couldn't find anything. Any pointers would be super helpful! Ultimately I want to prove that $S(n)$ approximates $n!$ as $n$ gets large so any ideas for showing that would also be helpful.

Answer 1

Here is a variation of Gonçalo's argument. We will consider the exponential generating function of the $S(n)$, namely

$$\begin{eqnarray*} f(z) = \sum_{n \ge 0} \frac{S(n)}{n!} z^n &=& \sum_{k \ge 0} \left( \sum_{n \ge 0} \frac{k^n z^n}{n!} \right) \exp(-k) \\ &=& \sum_{k \ge 0} \exp(kz) \exp(-k) \\ &=& \frac{1}{1 - \exp(z-1)}. \end{eqnarray*}$$

The asymptotic behavior of the coefficients of this generating function can be described using meromorphic singularity analysis, as described e.g. in Flajolet and Sedgewick. The dominant pole occurs at $z = 1$ with residue $1$ which is responsible for the dominant term of the asymptotic $\frac{S(n)}{n!} \approx 1$. The error is dominated by the poles next closest to the origin, which occur at $z = 1 \pm 2 \pi i$; this gives (with a bit more work)

$$\boxed{ \frac{S(n)}{n!} = 1 + O \left( \frac{1}{|1 + 2 \pi i|^n} \right) }.$$

Here $|1 + 2 \pi i| = \sqrt{1 + 4 \pi^2} = 6.362 \dots $ so, for example, for $n = 6$ the multiplicative error is on the order of $1.5 \times 10^{-5}$. However, this exponential does grow more slowly than the factorial so the additive error does eventually increase.

Closed forms for any particular $n$ can also be deduced by expanding the generating function. They are probably related to some well-known sequence of polynomials but I'm not sure which. I also don't know what bounds you get from applying Euler-Maclaurin here.

Answer 2

The difference may be expressed explicitly using the known expansion of the polylogarithm: \begin{align*} \sum\limits_{k = 1}^\infty {k^n \mathrm{e}^{ - k} } & = \operatorname{Li}_{ - n} (\mathrm{e}^{ - 1} ) = \Gamma (n + 1) + \sum\limits_{m = 0}^\infty {( - 1)^m \frac{{\zeta ( - n - m)}}{{m!}}} \\& = \Gamma (n + 1) - 2\sum\limits_{m = 1}^\infty {( - 1)^m \frac{{\Gamma (n + m)}}{{\Gamma (m)}}\frac{{\zeta (n + m)}}{{(2\pi )^{n + m} }}\cos \left( {\frac{{\pi (n + m)}}{2}} \right)} \end{align*} for $n\ne 0,-1,-2,\ldots$. Here $\zeta$ denotes the Riemann zeta function.

Addendum. I derive an estimate for the absolute error when $n>0$. From the above, we obtain \begin{align*} \frac{1}{{\Gamma (n + 1)}}\left( {\sum\limits_{k = 1}^\infty {k^n {\rm e}^{ - k} } } \right) - 1 & = 2\sum\limits_{m = 1}^\infty {( - 1)^{m - 1} \binom{n+m-1}{m-1}\frac{{\zeta (n + m)}}{{(2\pi )^{n + m} }}\cos \left( {\frac{{\pi (n + m)}}{2}} \right)} \\ & = \frac{{( - 1)^{n} }}{\pi }{\mathop{\rm Im}\nolimits} \sum\limits_{m = 0}^\infty {\binom{n+m}{m}\frac{{\zeta (n + m + 1)}}{{( - 2\pi {\rm i})^{n + m} }}} \\ & = \frac{{( - 1)^{n} }}{\pi }{\mathop{\rm Im}\nolimits} \sum\limits_{k = 1}^\infty {\frac{1}{k}\sum\limits_{m = 0}^\infty {\binom{n+m}{m}\frac{1}{{( - 2\pi {\rm i}k)^{n + m} }}} } \\ & = 2\,{\mathop{\rm Re}\nolimits} \sum\limits_{k = 1}^\infty {\frac{1}{{(2\pi {\rm i}k + 1)^{n + 1} }}} . \end{align*} Thus $$ \left| {\frac{1}{{\Gamma (n + 1)}}\left( {\sum\limits_{k = 1}^\infty {k^n {\rm e}^{ - k} } } \right) - 1} \right| \le 2\sum\limits_{k = 1}^\infty {\frac{1}{{\left| {2\pi {\rm i}k + 1} \right|^{n + 1} }}} = 2\sum\limits_{k = 1}^\infty {\frac{1}{{(4\pi ^2 k^2 + 1)^{(n + 1)/2} }}} . $$