Theorem: Stirling's Formula for Gamma Function

Let $\Gamma$ denote the Gamma function.
Let $z \in \C$ with a strictly positive real part and $\size {\arg z} < \dfrac \pi 2$.
Then:
:$\map \Gamma {z + 1} = \sqrt {2 \pi z} \, z^z e^{-z} \paren {1 + \dfrac 1 {12 z} + \dfrac 1 {288 z^2} - \dfrac {139} {51 \, 840 z^3} - \dfrac {571} {2\, 488 \, 320 z^4} + \dfrac {163 \, 879} {209\, 018 \, 880 z^5} + \dfrac {5\, 246\, 819} {75\, 246\, 796\, 800 z^6} + \cdots}$

Proof:
From Logarithmic Approximation of Error Term of Stirling's Formula for Gamma Function, we have:
:$\ds \map \Ln {\map \Gamma z} = \paren {z - \dfrac 1 2} \map \Ln z - z + \dfrac {\ln 2 \pi} 2 + \sum_{n \mathop = 1}^d \frac {B_{2 n} } {2 n \paren {2 n - 1} z^{2 n - 1} } + \OO \paren {\dfrac 1 {z^{2 d + 1} } }$
Taking the first $3$ terms in the sum, we get:
{{begin-eqn}}
{{eqn | l = \Ln \map \Gamma z
      | r = \paren {z - \dfrac 1 2} \map \Ln z - z + \dfrac {\ln 2 \pi} 2 + \sum_{n \mathop = 1}^d \frac {B_{2 n} } {2 n \paren {2 n - 1} z^{2 n - 1} } + \OO \paren {\dfrac 1 {z^{2 d + 1} } }
      | c = Logarithmic Approximation of Error Term of Stirling's Formula for Gamma Function
}}
{{eqn | r = \paren {z - \dfrac 1 2} \map \Ln z - z + \dfrac {\ln 2 \pi} 2 + \frac {B_2} {2 \times 1 \times z} + \frac {B_4} {4 \times 3 \times z^3}+ \frac {B_6} {6 \times 5 \times z^5}
      | c = stopping after $3$ terms
}}
{{eqn | r = \paren {z - \dfrac 1 2} \map \Ln z - z + \dfrac {\ln 2 \pi} 2 + \dfrac {\dfrac 1 6} {2 \times 1 \times z} + \dfrac {-\dfrac 1 {30} } {4 \times 3 \times z^3}+ \dfrac {\dfrac 1 {42} } {6 \times 5 \times z^5}
      | c = {{Defof|Bernoulli Numbers}}
}}
{{eqn | r = \paren {z - \dfrac 1 2} \map \Ln z - z + \dfrac {\ln 2 \pi} 2 + \dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5}
      | c = 
}}
{{end-eqn}}
Exponentiating both sides, we obtain:
{{begin-eqn}}
{{eqn | l = \map \exp {\Ln \map \Gamma z }
      | r = \map \exp {\paren {z - \dfrac 1 2} \map \Ln z - z + \dfrac {\ln 2 \pi} 2 + \dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }
      | c = 
}}
{{eqn | ll = \leadsto
      | l = \map \Gamma z
      | r = \map \exp {\paren {z - \dfrac 1 2} \map \Ln z - z + \dfrac {\ln 2 \pi} 2 + \dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }
      | c = Exponential Function is Inverse of Natural Logarithm
}}
{{eqn | r = \map \exp {\paren {\map \Ln z }^{\paren {z - \frac 1 2} } - z + \paren {\ln 2 \pi}^{\frac 1 2} + \dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }
      | c = Logarithm of Power - Natural Logarithm
}}
{{eqn | r = z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \map \exp {\dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }
      | c = Exponential Function is Inverse of Natural Logarithm, Product of Powers
}}
{{eqn | r = z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \paren {1 + \paren {\dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} } + \dfrac {\paren {\dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }^2} {2!} }
      | c = Power Series of Exponential Function
}}
{{eqn | o = 
      | ro= +
      | r = z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \paren {\dfrac {\paren {\dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }^3} {3!} + \dfrac {\paren {\dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }^4} {4!} }
      | c = 
}}
{{eqn | o = 
      | ro= +
      | r = z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \paren {\dfrac {\paren {\dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }^5} {5!} + \dfrac {\paren {\dfrac 1 {12 z} - \dfrac 1 {360 z^3} + \dfrac 1 {1260 z^5} }^6} {6!} + \cdots}
      | c = 
}}
{{eqn | r = z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \leftparen {1 + \dfrac 1 {12 z}  + \dfrac {\paren {\dfrac 1 {12 z} }^2} {2!} + \paren {-\dfrac 1 {360 z^3} + \dfrac {\paren {\dfrac 1 {12 z} }^3} {3!} } + }
      | c = grouping terms
}}
{{eqn | o = 
      | ro= +
      | r = \paren {\dfrac {2 \paren {\dfrac 1 {12 z} } \paren {- \dfrac 1 {360 z^3} } } {2!} +  \dfrac {\paren {\dfrac 1 {12 z} }^4} {4!} } + \paren {\dfrac 1 {1260 z^5} + \dfrac {3 \paren {\dfrac 1 {12 z} }^2 \paren {- \dfrac 1 {360 z^3} } } {3!} +  \dfrac {\paren {\dfrac 1 {12 z} }^5} {5!} } +
      | c =
}}
{{eqn | o = 
      | ro= +
      | r = \rightparen {\paren {\dfrac {2 \paren {\dfrac 1 {12 z} } \paren {\dfrac 1 {1260 z^5} } + \paren {-\dfrac 1 {360 z^3} }^2 } {2!} + \dfrac {4 \paren {\dfrac 1 {12 z} }^3 \paren {-\dfrac 1 {360 z^3} } } {4!} +  \dfrac {\paren {\dfrac 1 {12 z} }^6} {6!} } + \cdots }
      | c = 
}}
{{eqn | r = z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \leftparen {1 + \dfrac 1 {12 z}  + \dfrac 1 {2! \times \paren {12 z}^2} + \dfrac {-12^2 + 5} {3! \times 5 \times \paren {12 z}^3} + \dfrac {-4 \times 12^2 + 5} {4! \times 5 \times \paren {12 z}^4} + }
      | c = 
}}
{{eqn | o = 
      | ro = +
      | r = \rightparen {\dfrac {8 \times 12^4 -14 \times 12^2 + 7} {5! \times 7 \times \paren {12 z}^5} + \dfrac {254 \times 12^4 - 140 \times 12^2 + 35} {6! \times 5 \times 7 \times \paren {12 z}^6} + \cdots }
      | c = 
}}
{{eqn | r = z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \paren {1 + \dfrac 1 {12 z}  + \dfrac 1 {288 z^2} - \dfrac {139} {51 \, 840 z^3} - \dfrac {571} {2\, 488 \, 320 z^4} + \dfrac {163 \, 879} {209\, 018 \, 880 z^5} + \dfrac {5\, 246\, 819} {75\, 246\, 796\, 800 z^6} + \cdots }
      | c = 
}}
{{eqn | ll = \leadsto
      | l = z \map \Gamma z
      | r = z \times z^{\paren {z - \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \paren {1 + \dfrac 1 {12 z}  + \dfrac 1 {288 z^2} - \dfrac {139} {51 \, 840 z^3} - \dfrac {571} {2\, 488 \, 320 z^4} + \dfrac {163 \, 879} {209\, 018 \, 880 z^5} + \dfrac {5\, 246\, 819} {75\, 246\, 796\, 800 z^6} + \cdots }
      | c = multiplying both sides by z
}}
{{eqn | ll = \leadsto
      | l = \map \Gamma {z + 1}
      | r = z^{\paren {z + \frac 1 2} } \times e^{-z} \times \sqrt {2 \pi} \times \paren {1 + \dfrac 1 {12 z}  + \dfrac 1 {288 z^2} - \dfrac {139} {51 \, 840 z^3} - \dfrac {571} {2\, 488 \, 320 z^4} + \dfrac {163 \, 879} {209\, 018 \, 880 z^5} + \dfrac {5\, 246\, 819} {75\, 246\, 796\, 800 z^6} + \cdots }
      | c = Gamma Difference Equation, Product of Powers
}}
{{eqn | r = z^z \times e^{-z} \times \sqrt {2 \pi z} \times \paren {1 + \dfrac 1 {12 z}  + \dfrac 1 {288 z^2} - \dfrac {139} {51 \, 840 z^3} - \dfrac {571} {2\, 488 \, 320 z^4} + \dfrac {163 \, 879} {209\, 018 \, 880 z^5} + \dfrac {5\, 246\, 819} {75\, 246\, 796\, 800 z^6} + \cdots }
      | c = 
}}
{{end-eqn}}
{{qed}}