This is question I could have been asking for years. I wonder if the basis of this magnificent approximation has been published anywhere.
What I have seen (the problem is that I do not remember where) are approximations built as
$$n!\sim \sqrt{\pi}\left(\frac ne\right)^n\sqrt [2k]{P_k(n)}$$ where remains the mystery of the $\sqrt [2k]{.}$. The first case I saw is Gosper approximation.
The coefficients of the polynomials were obtained taking the logarithms of both sides and identified with the Stirling series. So, what was obtained are
$$P_1(n)=2n+12$$
$$P_2(n)=4n^2+12n+18$$
$$\color{red}{P_3(n)=8 n^3+4 n^2+n+\frac{1}{30}}$$ and for sure, with the tools we have today, we could continue for ever
$$P_4(n)=16 n^4+\frac{32 n^3}{3}+\frac{32 n^2}{9}+\frac{176 n}{405}-\frac{128}{1215}$$
$$P_5(n)=32 n^5+\frac{80 n^4}{3}+\frac{100 n^3}{9}+\frac{178 n^2}{81}-\frac{95
n}{972}+\frac{2143}{40824}$$
$$P_6(n)=64 n^6+64 n^5+32 n^4+\frac{128 n^3}{15}+\frac{8 n^2}{15}+\frac{8
n}{105}+\frac{596}{1575}$$ for more and more accuracy.
For example, for $n=5$, the magic formula given by the great Ramanujan gives $120.000147066$ while the last given here leads to $120.000000406$.
There is one thing interesting to notice : up to $k=3$ the coefficients of powers of $n >0$ are all integer numbers.
Edit
Staying with the power $\frac 16$ from Ramanujan, we could extend it as
$$P_3(n)=8 n^3+4 n^2+n+\frac{1}{30}\left(1+\sum_{i=1}^\infty \frac {a_i}{n^i}\right)$$ and the sequence of the first $a_i$'s is
$$\left\{-\frac{11}{8},\frac{79}{112},\frac{3539}{6720},-\frac{9511}{13440},-\frac{30
153}{71680},\frac{233934691}{212889600},\frac{3595113569}{5960908800},\cdots\right\}$$