What is special about polylogarithms that leads to so many interesting identities and applications?
The reason why polylogarithm are so important/interesting/ubiquitous is they are the simplest non trivial examples of analytical functions underlying variations of mixed Hodge structure. This goes back to Beilinson and Deligne.
A variation of mixed Hodge structure is a very sophisticated gadget. You can think of it as
- a nice differential equation (the underlying connexion)
- its solutions (the underlying local system of horizontal sections) with a $\mathbb{Q}$-structure
- some additional data that make the structure very rigid.
Typical examples of VMHS on $X$ come from the cohomology of families of varieties parametrized by $X$. They can be used to encode the interaction between between topological and arithmetical properties of $X$.
For example, there is a rank 2 variation of mixed Hodge structure $K \in Ext^1_{VMHS(\mathbb{C}^\times)}(\mathbb{Q},\mathbb{Q}(1))$ over $\mathbb{C}^\times$ known as the Kummer sheaf. The underlying local system $K_{\mathbb{Q}}$ has fiber $H_1(\mathbb{C}^\times,\{1,z\};\mathbb{Q})$ at $z$. The underlying connexion is a trivial vector bundle of rank 2 with nilpotent connexion $$ \nabla = d - \begin{pmatrix} 0 & 0 \cr % \frac{dt}{t} & 0 \end{pmatrix} $$ The "periods" are obtained by integrating the coefficient of the matrix over the paths $\gamma \in H_1(\mathbb{C}^\times,\{1,z\};\mathbb{Q})$. So we get a the non trivial period by integrating $\frac{dt}{t}$ over paths $[1,z]$, i.e. we get determinations of $\log(z)$. Conclusion: we have an object $K$ in $VMHS(\mathbb{C}^\times)$ "categorizing" the classical logarithm function. On the arithmetic side, the transcendance of $\log(z)$ for generic $z$ mirrors the fact that $H_1(\mathbb{C}^\times) = \mathbb{Z}$.
The same can be done for polylogarithm functions. The Logarithmic sheaf is the symmetric algebra $Log := Sym(K)$ (it corresponds to the whole family of $\log^n(z)$, $n\in \mathbb{N}$). The Polylogarithm sheaf is a canonical extension $Pol$ of $\mathbb{Q}$ by the restriction of $Log(1)$ to $\mathbb{P}^1\setminus \{0,1,\infty\}$. Its periods encodes the monodromy of the polylogarithm functions in the same way the Kummer sheaf does for the logarithm function. These are the most elementary unipotent variations of mixed Hodge structures.
Now we have "categorized" the classical polylogarithm functions. In fact, this can actually be done on a more fundamental level using only algebraic cycles defined over $\mathbb{Z}$ (this is the motivic story, the variation of mixed Hodge structure being just a realization of the motivic object). This has very interisting arithmetic consequences. For example, specializing to 1, this implies that we have motivic cohomology classes in $H^1(\mathbb{Z},\mathbb{Q}(n))$ whose images under the Hodge regulator corresponds to $\zeta(n)$. Using this picture, the period conjecture then implies that the $\zeta(2n+1)$ are algebraicly independent over $\mathbb{Q}$. To give you an idea of how powerful this intuition is: we can't even prove $\zeta(5)$ is irrational!
In conclusion, the polylogarithm functions are interesting because they correspond to non trivial algebraic cycles. This leads to interactions between the analytical properties of the functions, arithmetics of special values and algebraic geometry. Lots of classical functions should have similar interpretations. For example, there is a similar picture for Euler's Beta and Gamma functions.
One of the reasons is that it satisfies $$\Theta_z f = T_s f $$ where $\Theta_z = z\frac{\partial}{\partial z}$ and $T_s g(s) = g(s-1)$ (where you take $f(s,z) = Li_s(z)$, and is in some sense the simplest non-trivial function which satisfies this. This first-order mixed (Euler) differential-shift equation was studied some in the 1950s, but then ignored for a long time. That, of course, is not the only reason, but a number of 'weird' identities have simpler proofs when looked at this way.
Polylogarithms have interesting connections in the Theory of Partitions. This is mainly because a class of generating functions for bivariate partition statistics can be approximated by polylogarithms.
Consider for example the class of generating functions with $s>1,$
$\displaystyle \sum_{n=0}^\infty \sum_{m=0}^\infty p(n,m)x^n u^m =\prod_{k=1}^\infty \frac{1}{(1-xu^k)^{k^{s-2}}}. $
When $s=2,$ $p(n,m)$ is the number of partitions of $m$ with $n$ parts. When $s=3,$ $p(n,m)$ represents the number of plane partitions of $m$ with trace $n.$
For simplicity assume $1>x>0,$ as $u\to ^-1,$
$\displaystyle \prod_{k=1}^\infty \frac{1}{(1-xu^k)^{k^{s-2}}} \thicksim\sqrt[ - \frac{1}{\zeta(2-s)}]{1-x}\exp(\Gamma(s-1)Li_{s}(x)/(\ln u^{-1})^{s-1}). $
When used in conjunction with other techniques like the Circle Method, this estimate serves as a foundation for many facts about $p(n,m).$
Off the top of my head here are three papers that use this kind of estimate: http://www.math.drexel.edu/~rboyer/papers/partitions_experimental.pdf http://plms.oxfordjournals.org/content/s2-36/1/117.full.pdf http://www.springerlink.com/content/d65348128235x7m1/