Gaussian peak behind Stirling's formula ($t=n+\sqrt{n}s,t^n{e}^{-t}/(n^n{e}^{-n})\to e^{-s^2/2}$)

Repeated integration by parts gives the Gamma integralThe logarithm of its integrand is $n\log t-t$, whose unique maximum is at $t=n$. Write $t=n+\sqrt{n}s$. For fixed $s$ and large $n$, a Taylor expansion givesThe Big O notation here is for a fixed $s$, or uniformly for $s$ in a fixed bounded interval. It does not claim a uniform approximation over the entire growing integration domain. Consequently the peak, divided by its height, approaches $e^{-s^2/2}$.

The exact substitution givesReplacing the integral by the Gaussian integral ${\int }_{\mathbb{R}}{e}^{-s^2/2}ds=\sqrt{2\pi }$ produces $S_n$. Thus the estimate is peak height $n^n{e}^{-n}$ times width $\sqrt{n}$ times the area of a standard Gaussian peak. This is Laplace's method in a concrete form; controlling the integral's tails justifies the replacement. The next terms, including their products after exponentiation, produce the Stirling correction by Gaussian moments.