Knowing the Stirling's approximation for the Gamma function (factorial) for integers: $$\Gamma(n+1)=n!\approx \sqrt{2\pi n}n^ne^{-n}\bigg(1+\frac{a_1}{n}+\frac{a_2}{n^2}+\cdots\bigg)$$ Using the above approximation one can write: $$(n+1)!=\sqrt{2\pi(n+1)}(n+1)^{n+1}e^{-(n+1)}\bigg(1+\frac{a_1}{n+1}+\frac{a_2}{(n+1)^2}+\cdots\bigg)$$ We know that following recursion holds: $$(n+1)!=(n+1)n!$$ One can rewrite this: $$(n+1)!=(n+1)\sqrt{2\pi n}n^ne^{-n}\bigg(1+\frac{a_1}{n}+\frac{a_2}{n^2}+\cdots\bigg)$$ All this comes from: https://www.csie.ntu.edu.tw/~b89089/link/gammaFunction.pdf (Page 8-9) Then the author gives this expansion to calculate the $a_k$ coefficients when $n$ becomes large. Comparing these two expressions for $(n+1)!$ gives $$1+\frac{a_1}{n}+\frac{a_2}{n^2}+\cdots=\bigg(1+\frac{1}{n}\bigg)^{n+1/2}e^{-1}\bigg(1+\frac{a_1}{n+1}+\frac{a_2}{(n+1)^2}+\cdots\bigg)$$ Then he says, that after "classical series expansion" this equals: $$1+\frac{a_1}{n}+\frac{a_2-a_1+\frac{1}{12}}{n^2}+\frac{\frac{13}{12}a_1-2a_2+a_3+\frac{1}{12}}{n^3}+\cdots$$ I don't understand how he got there. Only thing that came to my mind was, that as $n\to\infty$ $\big(1+\frac{1}{n}\big)^n$ goes to $e$ but then we are left with $\big(1+\frac{1}{n}\big)^{1/2}$. When i expand this into binomial series, i get: $$\bigg(1+\frac{1}{2n}-\frac{1}{8n^2}+\frac{1}{16n^3}\cdots\bigg)\cdot\bigg(1+\frac{a_1}{n+1}+\frac{a_2}{(n+1)^2}+\cdots\bigg)$$ And I'm stuck here.
Or is there any other elementary way how to compute the coefficients $a_k$ of the stirling's series expansion for factorial/Gamma function?