Stirling’s formula: proof?

Suppose we want to show that $$n! \sim \sqrt{2 \pi} n^{n+(1/2)}e^{-n}$$

Instead we could show that $$\lim_{n \to \infty} \frac{n!}{n^{n+(1/2)}e^{-n}} = C$$ where $C$ is a constant. Maybe $C = \sqrt{2 \pi}$.

What is a good way of doing this? Could we use L’Hopital’s Rule? Or maybe take the log of both sides (e.g., compute the limit of the log of the quantity)? So for example do the following $$\lim_{n \to \infty} \log \left[\frac{n!}{n^{n+(1/2)}e^{-n}} \right] = \log C$$

A proof I found a while ago entirely relies on creative telescoping. Since $$\frac{1}{n^2}-\frac{1}{n(n+1)}=\frac{1}{n^2(n+1)}$$,

$$\begin{eqnarray*} \sum_{n\geq m}\frac{1}{n^2}&=&\sum_{n\geq m}\left(\frac{1}{n}-\frac{1}{(n+1)}\right)+\frac{1}{2}\sum_{n\geq m}\left(\frac{1}{n^2}-\frac{1}{(n+1)^2}\right)\\&+&\frac{1}{6}\sum_{n\geq m}\left(\frac{1}{n^3}-\frac{1}{(n+1)^3}\right)-\frac{1}{6}\sum_{n\geq m}\frac{1}{n^3(n+1)^3}\tag{1}\end{eqnarray*}$$
hence, by the series representation for $$\psi(z)=\frac{d}{dz}\log\Gamma(z)$$ (where $$\Gamma(z)$$ is the analytic continuation of $$\int_{0}^{+\infty}t^{z-1}e^{-t}\,dt$$, defined for $$\text{Re}(z)>0$$):
$$\psi'(m)=\sum_{n\geq m}\frac{1}{n^2}\leq \frac{1}{m}+\frac{1}{2m^2}+\frac{1}{6m^3}\tag{2}$$
and in a similar fashion:
$$\psi'(m) \geq \frac{1}{m}+\frac{1}{2m^2}+\frac{1}{6m^3}-\frac{1}{30m^5}.\tag{3}$$
Integrating twice, we have that $$\log\Gamma(m)$$ behaves like:
$$\log\Gamma(m)\approx\left(m-\frac{1}{2}\right)\log(m)-\color{red}{\alpha} m+\color{blue}{\beta}+\frac{1}{12m}\tag{4}$$
where $$\color{red}{\alpha=1}$$ follows from $$\log\Gamma(m+1)-\log\Gamma(m)=\log m$$.

That gives Stirling’s inequality up to a multiplicative constant.

$$\color{blue}{\beta=\log\sqrt{2\pi}}$$ then follows from Legendre’s duplication formula and the well-known identity:

$$\Gamma\left(\frac{1}{2}\right)=2 \int_{0}^{+\infty}e^{-x^2}\,dx = \sqrt{\pi}.\tag{5}$$

Addendum: if we apply creative telescoping like in the second part of this answer, i.e. by noticing that $$k(x)=\frac{60x^2-60x+31}{60x^3-90x^2+66x-18}$$ gives $$k(x)-k(x+1)=\frac{1}{x^2}+O\left(\frac{1}{x^8}\right)$$, we arrive at

$$\begin{eqnarray*} m!&\approx& 2^{\frac{37-32m}{42}}e^{\frac{1}{84} \left(42-\sqrt{35} \pi -84 m+2 \sqrt{35} \arctan\left[\sqrt{\frac{5}{7}} (2m-1)\right]\right)} \\ &\cdot&\sqrt{\pi}\, m\,(2m-1)^{\frac{8}{21}(2m-1)}\left(m^2-m+\frac{3}{5}\right)^{\frac{5}{84}(2m-1)}\tag{6}\end{eqnarray*}$$
that is much more accurate than the “usual” Stirling’s inequality, but also way less “practical”.
However, it might be fun to plug in different values of $$m$$ in $$(6)$$ to derive bizarre approximate identities involving $$e,\pi,\sqrt{\pi}$$ and values of the arctangent function, like

$$\sqrt{\pi}\,\exp\left[-\frac{1}{42}\left(147+\sqrt{35} \arctan\frac{1}{\sqrt{35}}\right)\right]\approx 2^{\frac{19}{6}}3^{\frac{1}{6}}5^{\frac{5}{12}} 7^{-\frac{37}{12}}.\tag{7}$$