How might one make this Gaussian integral derivation of Stirling's approximation more airtight? I will be following an argument presented by the YouTuber Flammable Maths - which is in fairness not their own, they're just giving simplified content to their audience!
Their derivation is appealingly intuitive, but I definitely feel there is some rigour lacking - however I am not able to follow the only other proof I have seen, so I will be trying to improve upon the argument they make.
We want to show: $$z!\sim z^ze^{-z}\cdot\sqrt{2\pi z}$$
We begin with: $$\begin{align}z!=\Gamma(z+1)&=\int_0^\infty e^{-t}t^z\,dt=\int_0^\infty\exp(z\ln(t)-t)\,dt\\&=\int_{-1}^\infty\exp(z\ln(z(1+x))-z(1+x))\cdot z\,dx\\&=z^{z+1}e^{-z}\int_{-1}^\infty\exp(z\cdot(\ln(1+x)-x))\,dz\end{align}$$
So far, no problems at all. However, these next lines of argument (paraphrased by me) feel... wishy-washy:

As the expression $\ln(1+x)-x\lt0$ for all non-zero $x$, the integrand is $\exp(-s)$ for some large quantity $s$, as $z\to\infty$. On the interval $(-1,0)$, $\ln(1+x)-x$ diverges rapidly to $-\infty$, even more so in multiplication by large $z$, and hence the exponential decays quickly to $0$. Likewise, on the interval $(0,\infty)$ the integrand goes sharply to $0$. We may then consider the integrand as non-negligible for $x$ in a small neighbourhood of $0$, and as such the Maclaurin expansion of $\ln(1+x)$ may be used to get the integrand as $\exp(-z\cdot\frac{x^2}{2})\cdot\exp(-z\cdot O(x^3))$ as $x\to0$. As $z$ gets larger, the neighbourhood in which the integrand is non-negligible gets smaller, and thus the $O(x^3)$ term vanishes. We now obtain: $$z!\approx z^{z+1}e^{-z}\int_{-\infty}^\infty\exp\left(-z\cdot\frac{x^2}{2}\right)\,dx$$ For large $z$, extending the lower bound to $-\infty$ as that extra area is again negligible. The Gaussian integral now returns $z!\approx z^{z+1}e^{-z}\sqrt{2\pi/z}=z^ze^{-z}\sqrt{2\pi z}$

There are many problems here, in my view. We first are discarding the majority of the range of integration, invoking a Maclaurin approximation but discarding the $O(x^3)$ term without proof that this does not affect things (for example, could I not simply include the $x^3$ term as well, losing the beloved Gaussian integral, and then truncate at $x^4$ instead), and then re-instating a much larger range of integration - the entire real line - of a different function. We also have not shown that the negligible portions of the integral remain negligible even under multiplication by $z^{z+1}e^{-z}$. We are also missing the important aspect of asymptotic equivalence: we have only (heuristically) shown that $z!\in\Theta(z^ze^{-z}\sqrt{2\pi z})$. The curious thing is, all the steps I balk at are in fact, as visually corroborated by Desmos, fairly good approximations anyway - but I have no idea how to show that, as there are too many moving parts and shortcuts in this argument for me to handle.
My question: is the argument rigorously salvageable, or should I give up and look to another approach?
Many thanks.
 A: Allow me to use $t$ instead of $z$ since $z$ is usually thought of as a complex variable. The argument can be made rigorous; here is a precise claim that we will prove using the outline you gave:
$$
\lim_{t\to\infty}\frac{\int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx}{\sqrt{2\pi/t}} = 1.
$$
In particular, we will prove the limit exists. As an immediate consequence, we will prove Stirling's approximation:
$$
\lim_{t\to\infty}\frac{t!}{t^te^{-t}\sqrt{2\pi t}} = 1.
$$
One way to make this kind of argument rigorous is to introduce a parameter $\newcommand{\eps}{\varepsilon} \eps > 0$ and split the region of integration into $(-1,-\eps)\cup(-\eps,\eps)\cup(\eps,\infty)$, and estimate each piece separately. We will always assume $\eps\in(0,1/10)$ is fixed.
(i) For all $x\in(-\eps,\eps)$, by considering the Maclaurin series of $\ln(1+x)$ one can show
$$(-x^2/2+x^3/3)(1+\eps)\le\ln(1+x)-x\le -x^2/2+x^3/3.\tag{1}$$
We also have the immediate inequality
$$-\frac{\eps x^2}{3}\le \frac{x^3}{3} \le \frac{\eps x^2}{3},\tag{2}$$
so we can estimate for an upper bound
$$
\int_{-\eps}^\eps e^{t[\ln(1+x)-x]}\,dx \le \int_{-\eps}^\eps e^{-tx^2/2}e^{tx^3/3}\,dx \le \int_{-\eps}^\eps e^{-t\frac{x^2}{2}(1-\frac{2\eps}{3})}\,dx\le\sqrt{2\pi (1/t)(1-2\eps/3)^{-1}}.
$$
We proceed analogously for the lower bound, but we use the other sides of the inequalities $(1)$ and $(2)$:
\begin{align*}
\int_{-\eps}^\eps e^{t[\ln(1+x)-x]}\,dx &\ge \int_{-\eps}^\eps e^{-(1+\eps)tx^2/2}e^{(1+\eps)tx^3/3}\,dx\\
&\ge \int_{-\eps}^\eps e^{-(1+\eps)tx^2/2}e^{-\eps(1+\eps)tx^2/3}\,dx\\
&=\int_{-\eps}^\eps \exp[{-(1+\eps)(1+\frac{2\eps}{3})\frac{tx^2}{2}}]\,dx\\
&\ge \int_{-\eps}^\eps\exp[-(1+\eps)^2\frac{tx^2}{2}]\,dx\\
&=(1+\eps)^{-1}\sqrt{2\pi/t} - 2\int_\eps^\infty\exp[-(1+\eps)^2\frac{tx^2}{2}]\,dx\\
&\ge (1+\eps)^{-1}\sqrt{2\pi/t} - 2\int_\eps^\infty\exp[-\frac{tx^2}{2}]\,dx.
\end{align*}
In the second-to-last line, we used the fact that $e^{ax^2}$ is an even function for any $a$. We need an upper bound for the last integral to continue bounding the expression we care about from below, so we estimate:
\begin{align*}
\int_\eps^\infty 1\cdot e^{-tx^2/2}\,dx \le \int_{\eps}^\infty \frac{x}{\eps}\cdot e^{-tx^2/2}\,dx,
\end{align*}
and by a routine $u$-substitution $u = tx^2/2$, we get
$$ 2\int_\eps^\infty\exp[-\frac{tx^2}{2}]\,dx\le \frac{2}{t\eps}e^{-t\eps^2/2}.$$
Therefore, our lower bound is
$$
\int_{-\eps}^\eps e^{t[\ln(1+x)-x]}\,dx \ge (1+\eps)^{-1}\sqrt{2\pi/t} - \frac{2}{t\eps}e^{-t\eps^2/2}.
$$
(ii) For $x\in (\eps,\infty)$, I find it somewhat more convenient to rewrite the integrand as
$$ \int_\eps^\infty (\frac{1+x}{e^x})^t\,dx.$$
For $x\in(\eps,\infty)$, we have $1+x\le e^{(1-\frac{\eps}{10})x}$.  We deduce
$$
0\le \int_\eps^\infty (\frac{1+x}{e^x})^t\,dx \le \int_\eps^\infty e^{-t\eps x/10}\,dx = \frac{10 e^{-t\eps^2/10}}{t\eps}.
$$
(iii) To handle $x\in (-1,-\eps)$, note that the upper bound of (1) holds unconditionally for all $x$, so we can estimate
$$
0\le \int_{-1}^{-\eps}e^{t[\ln(1+x)-x]}\,dx\le \int_{-1}^{-\eps} e^{-t(x^2/2-x^3/3)}\,dx \le e^{-t(\frac{\eps^2}{2}-\frac{\eps^3}{3})}.
$$
Putting these three estimates together, we see that for any $\eps\in(0,1/10)$ and all $t>0$,
\begin{align*}
(1+\eps)^{-1}\sqrt{2\pi/t} - \frac{2}{t\eps}e^{-t\eps^2/2} &\le \int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx\\
&\le \sqrt{2\pi (1/t)(1-2\eps/3)^{-1}} + \frac{10 e^{-t\eps^2/10}}{t\eps} + e^{-t(\frac{\eps^2}{2}-\frac{\eps^3}{3})}
\end{align*}
Dividing by $\sqrt{2\pi/t}$ and sending $t\to\infty$, we get
$$
(1+\eps)^{-1}\le \liminf_{t\to\infty}\frac{\int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx}{\sqrt{2\pi/t}}\le\limsup_{t\to\infty}\frac{\int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx}{\sqrt{2\pi/t}} \le \sqrt{(1-2\eps/3)^{-1}}.
$$
As $\eps\in(0,1/10)$ was arbitrary, we deduce that $\lim_{t\to\infty}\frac{\int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx}{\sqrt{2\pi/t}}$ exists and equals $1$.
Since it was shown that
$$
t! = t^{t+1}e^{-t}\int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx,
$$
we deduce that
\begin{align*}
\frac{t!}{t^te^{-t}\sqrt{2\pi t}} &= \frac{t^{t+1}e^{-t}\int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx}{t^te^{-t}\sqrt{2\pi t}}\\
&= \frac{\int_{-1}^\infty e^{t[\ln(1+x)-x]}\,dx}{\sqrt{2\pi/t}}\xrightarrow{t\to\infty} 1,
\end{align*}
as desired.
