(The presentation that follows is a very light introduction. For more austere results, consult Paper 1).

Let's see if we can find a nice expansion for Lambert's W function.

We know: W(z) has an infinite number of branches. Exactly one of these branches is analytic at 0. We then need a series expansion for W(z) around 0. Maple has the series for W built in. It is:
>series(W(z),z,8);

(If one wants to evaluate the coefficients of this series up to a finite degree, one can use Robert Israel's Maple code from his book: "Calculus the Maple Way" from Lab 16, page 135. The following code calculates the coefficients for n=5:

> restart;
> W:=LambertW;
> ys:=sum(a[n]*x^n,n=1..5);  #five terms
> ys*exp(ys)=x;     #implicit definition of the W
> g:=series(op(1,")-op(2,"),x=0,6);
> coeff(g,x,2);
> eqs:={seq(coeff(g,x,nn)=0,nn=1..5)};
> solve(eqs,{a[1],a[2],a[3],a[4],a[5]});
> subs(",ys);

However, both one's knowledge of the explicit a[n]'s above, either from Maple's expansion or from the explicit evaluation, are not very useful for the Lemma. We need the general term of the expansion.

So we tackle the problem by using some smart tricks in Maple:
The first obvious similarity comes from looking at the series expansion of the series z*e^-z, which I by accident stumped into while working on the infinite exponentials analysis:

Let's see the first series in Maple:
>series(z*exp(-z),z,8);

Now look at:
>series(W(z),z,8);

It is obvious that the series are suspiciously similar. We pick immediately the basic terms of the W expansion, by using some judicious observation, as:

>sum((-1)^(n-1)*z^n/(n-1)!, n=1..7);

Can we now compensate for the differences?
Yes, we can. Let us look at some more Maple code:

>s1:=series(W(z),z,15);
>s1:=convert(s1,polynom);
>s2:=sum((-1)^(n-1)*z^n/(n-1)!,n=1..14);
>s2:=convert(s2,polynom);

Now, let us extract the extra terms. For this, we compare terms of same degree in both expansions, and in fact, we can have Maple look for a pattern, by factoring the relevant exponents:

>for n from 2 to 14 do
>print('iterate',n,' factor',ifactor(op(n,s1)/op(n,s2)));
>od:

from the generated list, we can now pick up the extra factors, again by some simple observation, since they are already factored for us by Maple. Those extra factors will then be:
n^n-2.

Therefore, the general MacLaurin series expansion for the W, will be given by:

>sum((-1)^(n-1)*n^(n-2)*z^n/(n-1)!,n=1..oo);

We comb it a little bit to make it prettier by multiplying numerator and denominator by n, and it becomes:

>sum((-1)^(n-1)*n^(n-1)*z^n/n!,n=1..oo);

Here then, is the infamous Lambert's W function, pretty formatted:

Now we can easily find the radius of convergence for the expansion using the ratio test:

lim_n->oo|a_n+1/a_n| =

(and after some simple calculations...)

lim_n->oo|(1+1/n)^n-1 z|

= |e z|

<1, provided |z|<1/e.

We are done. The disk of convergence where the above expansion of the W is valid, is therefore: R = {z: |z|<1/e} [*].

(Naturally, the radius of convergence could not be greater than 1/e, since the W has a branch point at -1/e}

If |z| = 1/e then
|(-1)^(n-1)*n^(n-1)*zⁿ/n!| =
n^(n-1)/[eⁿ*n!],

and using Stirling's approximation:
Sqrt(2*Pi)*n^(n+1/2)/eⁿ < n!,

we get:

n^(n-1)/[eⁿ*n!] < n^(n-1)/[Sqrt(2*Pi)*n^(n+1/2)]
= 1/[Sqrt(2*Pi)*n^3/2].

Since the series sum(1/n^3/2, n=1..oo) converges,

the original series converges absolutely on all points of the circle of convergence: C = {z: |z|=1/e} as well. QED.

(Thanks to Paul D. Hanna for verifying this)