Consider the following integral:
\[\begin{eqnarray} \int_0^\infty dx x^{n}e^{-x}= \left[(-1)^n\frac{d^n}{d\alpha^n}\int_0^\infty dx e^{-\alpha x}\right]_{\alpha=1} =\left[(-1)^n\frac{d^n}{d\alpha^n} \frac{1}{\alpha}\right]_{\alpha=1}=n!. \tag{1} \end{eqnarray}\] Taking this definition, we can do the following: \[\begin{eqnarray} n!=\int_0^\infty dx x^{n}e^{-x}=\int_0^\infty dx e^{n ln(x)-x}. \tag{2} \end{eqnarray}\] Let’s take a close look at the function in the exponent: \[\begin{eqnarray} u(x) &=&n ln(x)-x, \tag{3} \end{eqnarray}\] as shown in Fig. 1. This function has its peak value at \(x=n\). Note that this function appears in the exponent, under the integral. The dominant contribution to the integral will come from the domain around \(x=n\). we can expand \(u(x)\) around \(x=n\):
\[\begin{eqnarray} u(x) &=&n ln(x)-x =n ln(x-n +n)-x=n ln(n[1 +\frac{x-n}{n}])-x\nonumber\\ & \simeq& n \left( ln(n)+\frac{x-n}{n} -\frac{1}{2}\left[\frac{x-n}{n}\right]^2 \right)-x =n ln(n) -n -\frac{1}{2}\frac{(x-n)^2}{n}\equiv \tilde{u}(x). \tag{4} \end{eqnarray}\]
The original function and the approximated functions are plotted in Fig. 1.
From Fig. 1, we also notice that if we extended the \(x\) range to include negative values, the integral would not change much since \(e^{\frac{1}{2}\frac{(x-n)^2}{n}}\) is rapidly decaying. Therefore we can change the lower limit of the integral from \(0\) to \(-\infty\) to get:
\[\begin{eqnarray} n!&=&\int_0^\infty dx x^{n}e^{-x}=\int_0^\infty dx e^{u(x)}\simeq \int_0^\infty dx e^{\tilde{u}(x)}=n^n e^{-n} \int_0^\infty dx e^{-\frac{1}{2}\frac{(x-n)^2}{n}}\nonumber\\ &\simeq&n^n e^{-n} \int_{-\infty }^\infty dx e^{-\frac{1}{2}\frac{(x-n)^2}{n}}=n^n e^{-n} \sqrt{2\pi n}=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n. \tag{5} \end{eqnarray}\]
Figure 2 shows the comparison of \(n!\) with the Stirling’s approximation given in Eq. (5).