Series expansion of incomplete gamma function ratio

If we start with a little rearranging,

$$ \begin{align} \frac{\Gamma(k+1,a)}{k!} &= \frac{1}{k!} \int_a^\infty t^k e^{-t}\,dt \\ &= \frac{1}{k!} \left( \int_0^\infty t^k e^{-t}\,dt - \int_0^a t^k e^{-t}\,dt \right) \\ &= 1 - \frac{1}{k}\int_0^a t^k e^{-t}\,dt, \end{align} $$

then make the change of variables $t = ae^{-s}$, we end up with

$$ 1 - \frac{\Gamma(k+1,a)}{k!} = \frac{a^{k+1}}{k!} \int_0^\infty e^{-ks} \exp\{-s-ae^{-s}\}\,ds. $$

If we expand

$$ \begin{align} \exp\{-s-ae^{-s}\} &= \sum_{n=0}^{\infty} c_n s^n \\ &= e^{-a} + e^{-a}(a-1) s + \tfrac{1}{2} e^{-a} (a^2-3a+1) s^2 + \cdots \end{align} $$

then by Watson's lemma we can exchange the order of integration and summation to obtain an asymptotic series for the integral;

$$ \begin{align} 1 - \frac{\Gamma(k+1,a)}{k!} &\approx \frac{a^{k+1}}{k!} \sum_{n=0}^{\infty} c_n \int_0^\infty e^{-ks} s^n \,ds \\ &= \frac{a^{k+1}}{k!} \sum_{n=0}^{\infty} \frac{n! c_n}{k^{n+1}} \\ &= \frac{a^{k+1} e^{-a}}{k \cdot k!} \left(1 + \frac{a-1}{k} + \frac{a^2 - 3a + 1}{k^2} + \cdots \right). \end{align} $$

In particular, to first order we have

$$ 1 - \frac{\Gamma(k+1,a)}{k!} \sim \frac{a^{k+1} e^{-a}}{k \cdot k!}. $$

Using a CAS, what it found as expansion for large values of $k$ is $$S(k)=\frac{\Gamma(k+1,a)}{k!}\simeq 1-\frac{ a^{k+1} e^{k \left(\log \left(\frac{1}{k}\right)+1\right)-a}}{\sqrt{2 \pi }k^{3/2}}$$ I hope this could be of some help.