Polylogarithms of negative integer order The polylogarithms of order $s$ are defined by
$$\mathrm{Li}_s (z) = \sum_{k \geqslant 1} \frac{z^k}{k^s}, 
\quad |z| < 1.$$
From the above definition, derivatives for the polylogarithms immediately follow. Here we have
$$\frac{d}{dz} \mathrm{Li}_s (z) = \frac{\mathrm{Li_{s - 1}}(z)}{z}.$$
When $s = 0$ the following series which is geometric results
$$\mathrm{Li}_0(z) = \sum_{k \geqslant 1} z^k = \frac{z}{1 - z}.$$
Starting with the expression for $\mathrm{Li}_0(z)$, repeated application of the derivative rule allows explicit expressions for the polylogarithms to all negative integer orders to be found. The first few are:
\begin{align*}
\mathrm{Li}_{-1} (z) &= \frac{z}{(1 - z)^2}\\
\mathrm{Li}_{-2}(z) &= \frac{z(1 + z)}{(1 - z)^3}\\
\mathrm{Li}_{-3}(z) &= \frac{z(1 + 4z +z^2)}{(1 - z)^4}\\
& \vdots
\end{align*}
All are rational functions in $z$. A closed form expression for $\mathrm{Li}_{-n}(z)$ where $n \in \mathbb{N}$ in terms of either Stirling numbers or Eulerian numbers can be found. In terms of Eulerian numbers $\left\langle n \atop k\right\rangle$ the result is
$$\mathrm{Li}_{-n} (z) = \frac{1}{(1 - z)^{n + 1}} \sum^n_{k = 0} \left\langle n \atop k\right\rangle z^{n - k}, \quad n \geqslant 1.$$
What I wish to know is how can such a result be proved. Of course once one knows what one is looking for, the result can be readily guessed, and I suppose induction used to prove the claim is true for all $n$, but is there a way to prove the result directly?    
 A: Suppose we seek to show that
$$\mathrm{Li}_{-n}(w) = \frac{1}{(1-w)^{n+1}}
\sum_{k=0}^{n} \left\langle {n\atop  k}\right\rangle w^{n-k}.$$
We are done if we can prove that
($\left\langle {n\atop  n}\right\rangle$ is zero)
$$[w^m] \frac{1}{(1-w)^{n+1}}
\sum_{k=0}^{n-1} \left\langle {n\atop  k}\right\rangle w^{n-k}
= m^n.$$
Extracting coefficients we get
$$\sum_{k=0}^{n-1} \left\langle {n\atop k}\right\rangle
{m+k-n + n\choose n}
= \sum_{k=0}^{n-1} \left\langle {n\atop k}\right\rangle
{m+k\choose n}.$$
Start from the bivariate generating function of the Eulerian numbers,
which seems like a reasonable starting point and which is
$$\frac{u-1}{u-\exp((u-1)z)}.$$
Now introduce
$${m+k\choose n} = [v^n] (1+v)^{m+k}$$
which gives for the sum
$$[v^n] (1+v)^m
\sum_{k=0}^{n-1} \left\langle {n\atop k}\right\rangle
(1+v)^k.$$
Substitute the  generating function of the Eulerian  numbers into this
to get
$$[v^n] (1+v)^m
n! [z^n] \frac{v}{1+v-\exp(vz)}
\\ = [v^n] (1+v)^m
n! [z^n]
\frac{1}{1-(\exp(vz)-1)/v}.$$
Observe that $(\exp(vz)-1)/v$ expanded  about zero with respect to $z$
starts at $z.$ Since we  are extracting the $[z^n]$ coefficient we may
write
$$[v^n] (1+v)^m n! [z^n]
\sum_{p=0}^n \frac{(\exp(vz)-1)^p}{v^p}.$$
Now the coefficient extraction from the term in $z$ is
$$n! [z^n] (\exp(vz)-1)^p =
v^n  \times n! \times \frac{1}{v^n} [z^n] (\exp(vz)-1)^p
\\ = v^n \times n! [z^n] (\exp(z)-1)^p
= v^n \times p! \times {n\brace p}.$$
Returning to the main sum we thus obtain
$$[v^n] (1+v)^m
\sum_{p=0}^n v^{n-p} \times p! \times {n\brace p}
\\ = \sum_{p=0}^n [v^p] (1+v)^m \times p! \times {n\brace p}
= \sum_{p=0}^n {m\choose p} \times p! \times {n\brace p}.$$
To conclude observe that
$$m^n = n! [z^n] \exp(mz)
= n! [z^n] \sum_{p=0}^m {m\choose p} (\exp(z)-1)^p
\\ = \sum_{p=0}^m {m\choose p} \times p! \times {n\brace p}.$$
Since  the  effective  upper  limit  is always  $\min(n,m)$  which  is
symmetric we may use  either $m$ or $n$ as the upper  limit and we are
done.
Alternatively when $n \lt m$ we may change the upper limit from $m$ to
$n$ because the Stirling number is zero for $n\lt p \le m.$ If $n \gt
m$ we may also set the upper limit to $n$ because the binomial
coefficient is zero for $n\ge p \gt m.$

Remark.  This  last sum  may  also  be evaluated  combinatorially.
Count all functions from $[1,n]$ to  $[1,m].$ To do this first fix the
number  $p$ of different  values that  will occur,  next choose  a set
partition of $n$ into $p$ sets  the constituents of each of which will
be mapped  to the same value,  and finally choose the  $p$ values from
$m$ to use and some ordering thereof, giving the result
$$m^n =
\sum_{p=0}^n {n\brace p} \times {m\choose p} \times p!.$$
A: The Wiki page presenting polylogarithms states identities with Stirling numbers of the second kind, one of them being
\begin{align*}
\mathrm{Li_{-n}}(z)=\sum_{k=0}^{n}k!
\begin{Bmatrix}n+1\\k+1\end{Bmatrix}
\left(\frac{z}{1-z}\right)^{k+1}\qquad\qquad n\geq 0\tag{1}
\end{align*}
and an identity with Eulerian numbers
\begin{align*}
\mathrm{Li_{-n}}(z)=\frac{1}{(1-z)^{n+1}}\sum_{k=0}^{n-1}k!
\left\langle n \atop k \right\rangle
z^{n-k}\qquad\qquad n\geq 1\tag{2}
\end{align*}
According to OPs question I'll try to provide a plausible reasoning how these formulae could be derived.

First part: Stirling numbers of the Second kind


In order to derive the identity (1) a crucial observation is the recurrence relation of the polylogarithm regarding the differential operator $\mathrm{D}=\frac{d}{dz}$. We can write for $n\geq 1$
\begin{align*}
\mathrm{Li_{-n}}(z)=(z\mathrm{D})\mathrm{Li_{-n+1}}(z)
\end{align*}


From this relation we obtain by successive application of the operator $z\mathrm{D}$ as already shown in OPs question


\begin{align*}
\mathrm{Li_{-n}}(z)=(z\mathrm{D})\mathrm{Li_{-n+1}}(z)=\ldots
=(z\mathrm{D})^n\mathrm{Li_{0}}(z)=(z\mathrm{D})^n\frac{z}{1-z}\tag{3}
\end{align*}

The connection with Stirling numbers of the second kind is given via the operator $z\mathrm{D}$, since the following well known operator relation is valid
\begin{align*}
(z\mathrm{D})^n=\sum_{k=1}^n\begin{Bmatrix}n\\k\end{Bmatrix}z^k\mathrm{D}^k\qquad\qquad n\geq 1\tag{4}
\end{align*}
Information about this relation can be found e.g. in Stirling numbers of the second kind by R. Milson.

Applying the relation (4) to $\frac{z}{1-z}$ we obtain for $n\geq 1$
\begin{align*}
(z\mathrm{D})^n\frac{z}{1-z}&=\sum_{k=1}^n\begin{Bmatrix}n\\k\end{Bmatrix}z^k\mathrm{D}^k\frac{z}{1-z}\\
&=\sum_{k=1}^n\begin{Bmatrix}n\\k\end{Bmatrix}z^k\mathrm{D}^k\sum_{j=1}^{\infty}z^{j}\\
&=\sum_{k=1}^n\begin{Bmatrix}n\\k\end{Bmatrix}z^k\sum_{j=k}^{\infty}k!\binom{j}{k}z^{j-k}\tag{5}\\
&=\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}z^k\sum_{j=0}^{\infty}\binom{j+k}{k}z^{j}\\
&=\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}z^k\sum_{j=0}^{\infty}\binom{-(k+1)}{j}(-z)^j\tag{6}\\
&=\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\frac{z^k}{(1-z)^{k+1}}\tag{7}
\end{align*}

Comment:

*

*In (5) we differentiate $k$ times and use $j(j-1)\cdots(j-k+1)=k!\binom{j}{k}$


*In (6) we use $\binom{-n}{k}=\binom{n+k-1}{k}(-1)^k$


*In (7) we apply the binomial series

Combining the polylogarithm representation (3) with (7) we obtain


\begin{align*}
\mathrm{Li_{-n}}(z)=\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\frac{z^k}{(1-z)^{k+1}}
\qquad\qquad n\geq 1\tag{8}
\end{align*}

We have found a relationship beween the polylogarithm and the Stirling numbers of the second kind. Note, this is not the representation (1) stated in the Wiki page, but we will show equality. In order to do so we recall the recurrence relation of Stirling numbers of the second kind
\begin{align*} 
\begin{Bmatrix}n+1\\k\end{Bmatrix}=k\begin{Bmatrix}n\\k\end{Bmatrix}+\begin{Bmatrix}n\\k-1\end{Bmatrix}
\qquad\qquad n\geq 0, k\geq 1
\end{align*}
and claim

The following holds true
\begin{align*}
\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\frac{z^k}{(1-z)^{k+1}}=
\sum_{k=0}^{n}k!\begin{Bmatrix}n+1\\k+1\end{Bmatrix}\left(\frac{z}{1-z}\right)^{k+1}
\qquad\qquad n\geq 1
\end{align*}


We start with the RHS and get
\begin{align*}
\sum_{k=0}^nk!&\begin{Bmatrix}n+1\\k+1\end{Bmatrix}\left(\frac{z}{1-z}\right)^{k+1}\\
&=\sum_{k=0}^nk!\left((k+1)\begin{Bmatrix}n\\k+1\end{Bmatrix}+
\begin{Bmatrix}n\\k\end{Bmatrix}\right)\left(\frac{z}{1-z}\right)^{k+1}\tag{9}\\
&=\sum_{k=0}^n(k+1)!\begin{Bmatrix}n\\k+1\end{Bmatrix}\left(\frac{z}{1-z}\right)^{k+1}
+\sum_{k=0}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\left(\frac{z}{1-z}\right)^{k+1}\\
&=\sum_{k=1}^{n+1}k!\begin{Bmatrix}n\\k\end{Bmatrix}\left(\frac{z}{1-z}\right)^{k}
+\sum_{k=0}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\left(\frac{z}{1-z}\right)^{k+1}\tag{10}\\
&=\sum_{k=1}^{n}k!\begin{Bmatrix}n\\k\end{Bmatrix}\left[\left(\frac{z}{1-z}\right)^{k}
+\left(\frac{z}{1-z}\right)^{k+1}\right]\\
&=\sum_{k=1}^{n}k!\begin{Bmatrix}n\\k\end{Bmatrix}\frac{z^k}{(1-z)^{k+1}}\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\Box
\end{align*}

Comment:

*

*In (9) we apply the recurrence relation to $\begin{Bmatrix}n+1\\k+1\end{Bmatrix}$


*In (10) we shift the index in the left sum by one and note that in the left sum the upper limit $k=n+1$ gives a factor $\begin{Bmatrix}n\\n+1\end{Bmatrix}=0$ and in the right sum the lower limit $k=0$ gives a factor $\begin{Bmatrix}n\\0\end{Bmatrix}=0$. This is respected in the following line.

Conclusion: The differential operator $z\mathrm{D}$  provides a convenient connection between the polylogarithm $\mathrm{Li_{-n}}(z)$ of negative integer order and the Stirling numbers of the second kind.

The Eulerian numbers are strongly connected to Stirling numbers of the second kind. So, it's plausible that there is also a nice representation of the polylogarithm in terms of Eulerian numbers.

Second step: Eulerian numbers


In order to show the relation (2) stated in the Wiki page, we consider the Eulerian polynomials $A_n(z)$ defined as
\begin{align*}
A_n(z)=\sum_{k=0}^{n-1}\left\langle n \atop k \right\rangle z^k\qquad\qquad\qquad n\geq 1
\end{align*}
These polynomials are linked with Stirling number of the second kind by a Theorem by Frobenius. See e.g. the Identity section in Eulerian polynomials. The theorem states
\begin{align*}
A_n(z)=\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}(z-1)^{n-k}\qquad\qquad n\geq 1
\end{align*}

With this theorem it's easy to derive the relation (2). We claim

The following holds true
\begin{align*}
\frac{1}{(1-z)^{n+1}}\sum_{k=0}^{n-1}\left\langle n \atop k \right\rangle z^{n-k}=
\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\frac{z^k}{(1-z)^{k+1}}
\qquad\qquad n\geq 1
\end{align*}


This is valid since we obtain for $n\geq 1$
\begin{align*}
\frac{1}{(1-z)^{n+1}}\sum_{k=0}^{n-1}\left\langle n \atop k \right\rangle z^{n-k}
&=\frac{z^n}{(1-z)^{n+1}}\sum_{k=0}^{n-1}\left\langle n \atop k \right\rangle \left(\frac{1}{z}\right)^k\\
&=\frac{z^n}{(1-z)^{n+1}}\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\left(\frac{1}{z}-1\right)^{n-k}\tag{11}\\
&=\frac{z^n}{(1-z)^{n+1}}\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\left(\frac{1-z}{z}\right)^{n-k}\\
&=\sum_{k=1}^nk!\begin{Bmatrix}n\\k\end{Bmatrix}\frac{z^k}{(1-z)^{k+1}}\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\Box
\end{align*}

Comment:

*

*In (11) we apply Frobenius Theorem.


*Note also that we used here as a matter of convenience the derived identity (8) and not the identity stated in the Wiki page.
