Stirling's formula: proof?

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}$$

Inspired by the two references below, here's a short proof stripped of motivation and details.

For $t>0$, define

$$ g_t(y) = \begin{cases} \displaystyle \left(1+\frac{y}{\sqrt{t}}\right)^{\!t} \,e^{-y\sqrt{t}} & \text{if } y>-\sqrt{t}, \\ 0 & \text{otherwise}. \end{cases} $$

It is not hard to show that $t\mapsto g_t(y)$ is decreasing for $y\geq 0$ and increasing for $y\leq 0$. The limit is $g_\infty(y)=\exp(-y^2/2)$, so dominated convergence gives

$$ \lim_{t\to\infty}\int_{-\infty}^\infty g_t(y)\,dy = \int_{-\infty}^\infty \exp(-y^2/2)\,dy = \sqrt{2\pi}. $$

Use the change of variables $x=y\sqrt{t}+t$ to get

$$ t! = \int_0^\infty x^t e^{-x}\,dx = \left(\frac{t}{e}\right)^t \sqrt{t} \int_{-\infty}^\infty g_t(y)\,dy. $$

References:

[1] J.M. Patin, A very short proof of Stirling's formula, American Mathematical Monthly 96 (1989), 41-42.

[2] Reinhard Michel, The $(n+1)$th proof of Stirling's formula, American Mathematical Monthly 115 (2008), 844-845.