<s>[INST] How to expand $\sum_{k=0}^\infty \frac{k^k}{k!}e^{-k\lambda}$ I know that  $\sum_{k=0}^\infty \dfrac{k^k}{k!}\dfrac{e^{-k-k\tfrac{s}{n}}}{\sqrt{n}}$ can be expanded in powers of $\frac{1}{\sqrt{n}}$ to yield:
$\frac1{\sqrt{2s}}+\frac{1}{3\sqrt{n}}+\frac{\sqrt{2s}}{12n}+\dotso$
My question is twofold: (i) is there an actual proof of this (beyond numerical),
                        (ii) how does this expansion continue (general formula, or 
                             at least the next few terms)?
 [/INST] We will assume $s > 0$ or the series diverges.
Recall for $|z| \le \frac1e$, the $0$-branch of the Lambert W function has following expansion:
$$W_0(z) = \sum_{k=1}^\infty \frac{(-k)^{k-1}}{k!} z^k$$
Pick any $\lambda_0 > 1$ and substitute $z$ by $-e^{-\lambda}$ for $\lambda > \lambda_0$. Above expansion converges uniformly for such $\lambda$. So does the series of its derivatives for each individual term. In this situation, it is legal to exchange the order of summation and taking derivatives. This implies
$$f(\lambda) \stackrel{def}{=}\sum_{k=0}^\infty \frac{k^k}{k!} e^{-\lambda k}
= 1 + \sum_{k=1}^\infty \frac{k^k}{k!} e^{-\lambda k}
= 1 + \frac{d}{d\lambda} W_0(-e^{-\lambda})
$$ 
Let $u = W_0(-e^{-\lambda})$, we have $u e^u = -e^{-\lambda}$. Differentiate against $\lambda$, we get
$$(u+1)e^u \frac{du}{d\lambda} = e^{-\lambda} = -u e^u 
\implies \frac{du}{d\lambda} = - \frac{u}{1+u}$$
This leads to
$$f(\lambda) = 1 + \frac{du}{d\lambda} = \frac{1}{1+u} = \frac{1}{1 + W_0(-e^{-\lambda})}$$
It is easy to see the series at hand evaluates to $\displaystyle\;\frac{1}{\sqrt{n}}f\left(1 + \frac{s}{n}\right)$
For simplicity, let's assume $n = 1$. Let $\lambda = 1 + s$ and $w = 1 + u$. We have 
$$f(1 + s) = \frac{1}{w}\quad\text{ and }\quad ue^u = -e^{-\lambda} 
\iff (1-w) e^w = e^{-s}$$ This leads to
$$s = -(\log(1-w) + w) = \frac{w^2}{2} + \frac{w^3}{3} + \frac{w^4}{4} + \cdots
$$
For small $s$, the leading approximation of $w$ is clearly $\sqrt{2s}$ and hence
$f(1+s) \sim \frac{1}{\sqrt{2s}}$.
To process further, write $s$ as $\frac{t^2}{2}$ and expanse $w$ as a power series in $t = \sqrt{2s}$. By matching coefficients of different powers of $t$ using a CAS, we get
$$\begin{align}w &= 
t-\frac{t^2}{3}+\frac{t^3}{36}+\frac{t^4}{270}+\frac{t^5}{4320}-\frac{t^6}{17010} + \cdots\\
\implies
\frac1w &= 
\frac{1}{t}+\frac{1}{3}+\frac{t}{12}+\frac{2t^2}{135}+\frac{t^3}{864}-\frac{t^4}{2835}-\frac{139 t^5}{777600} + \cdots
\end{align}
$$
Putting $n$ back, we get
$$\frac{1}{\sqrt{n}}f\left(1 + \frac{s}{n}\right)
= \frac{1}{\sqrt{2s}} + \frac{1}{3n^{\frac12}}
+ \frac{\sqrt{2s}}{12n}
+\frac{4\,s}{135\,{n}^{\frac{3}{2}}}
+\frac{\sqrt{2s}^3}{864 n^2}
-\frac{4\,{s}^{2}}{2835\,{n}^{\frac{5}{2}}} + \cdots
$$
It looks like there is a typo in the third term of the expansion you get.
</s>