I'm trying to derive the Stirling Series from the gamma function
\begin{equation} \begin{split} \Gamma (N+1) &= \int_0^\infty x^N e^{-x} dx \\ &= \int_0^\infty e^{N (\ln x - \frac{x}{N})} dx \\ \end{split} \end{equation}
Thus the $\phi(x)$ for the saddle point approximation here is $\ln{x} - \frac{x}{N}$. Expanding this around it's maximum at $x_{max} = N$
\begin{equation} \begin{split} \phi(x) &\approx \ln{N}-1 + \phi'(x_{max}) + \frac{1}{2!}\left(\frac{-1}{N^2}\right)\left(x-N\right)^2 + \frac{1}{3!} \left( \frac{2}{N^3} \right) \left(x-N\right)^3 - \frac{6}{4!N^4} (x-N)^4 + \frac{24}{5! N^3} (x-N)^5 + ...\\ &\approx \ln{N}-1 - \frac{1}{2N^2}\left(x-N\right)^2 + \frac{1}{3N^3} \left(x-N\right)^3 - \frac{1}{4N^4} (x-N)^4 + \frac{1}{5N^3} (x-N)^5 + ...\\ \end{split} \end{equation}
Back to the gamma function
\begin{equation} \begin{split} \Gamma (N+1) &= \int_0^\infty e^{N (\ln N - \frac{x}{N})} dx \\ &= \int_0^\infty exp \left(N\left(\ln{N}-1 - \frac{1}{2N^2}\left(x-N\right)^2 + \frac{1}{3N^3} \left(x-N\right)^3 - \frac{1}{4N^4} (x-N)^4 + \frac{1}{5N^3} (x-N)^5 + ... \right)\right) dx \\ \\ &= \left( \frac{N}{e} \right)^N \int_0^\infty exp \left(- \frac{1}{2N}\left(x-N\right)^2 + \frac{1}{3N^2} \left(x-N\right)^3 - \frac{1}{4N^3} (x-N)^4 + \frac{1}{5N^4} (x-N)^5 + ... \right) dx \\ &= \left( \frac{N}{e} \right)^N \int_0^\infty \prod_{n=0}^\infty e^{\:\frac{(-1)^n}{nN^{n-1}} (x-N)^n} dx \\ \end{split} \end{equation}
I am lost as to how to proceed from here. Maybe expanding all of the exponentials? I have a few other ideas, but as to not clutter I'll leave it there.