Pedagogy, polylogarithms and particles
April 30, 2015. The polylogarithm series often appears as the butt of convergence tests; here, we go further, evaluating it and discussing the application to quantum statistics.
Introduction
Prerequisites: basic real analysis, recurrence relations, power series, quantum statistics (last section).
A few years back, James Sellers wrote a nice article for teachers of first analysis courses, which calculates the exact value of
\[\sum_{k\geq1} \frac{k^n}{m^k}\]for $n \in \mathbb{N}$ and $m > 1$. This series usually appears only to have the ratio test applied to it, and is promptly forgotten. Sellers points out that engaged students, armed only with power series or recurrence relations, can actually do the summation, so these “toy” series are something of a missed pedagogical opportunity. We discuss Sellers’s approach, then the connection to the polylogarithm function and quantum statistical mechanics.
Recurrence relation
Sellers’s first method for finding exact values uses a recurrence relation. Let’s define
\[a_m(n) := \sum_{k\geq1} \frac{k^n}{m^k}.\]Applying the convergence test shows that this converges if $m > 1$. To seed our recurrence, we start off with the geometric series ($n = 0$):
\[a_m(0) = \sum_{k\geq1} \frac{1}{m^k} = \frac{1}{m1}.\]Using the binomial theorem, we get an equation for $a_m(n)$:
\[\begin{align*} a_m(n) & = \sum_{k\geq1} \frac{k^n}{m^k} \\ & = \frac{1}{m} + \sum_{k\geq2} \frac{k^n}{m^k}\\ & = \frac{1}{m} + \frac{1}{m}\sum_{k\geq1} \frac{(k+1)^n}{m^k}\\ & = \frac{1}{m}\left[1 + \sum_{k\geq1} \frac{1}{m^k}\sum_{j=0}^n \binom{n}{j}k^j\right] \\ & = \frac{1}{m}\left[1 + \sum_{j=0}^n \binom{n}{j}\sum_{k\geq1} \frac{k^j}{m^k}\right] = \frac{1}{m}\left[1 + \sum_{j=0}^n \binom{n}{j}a_m(j)\right]. \end{align*}\]Solving for $a_m(n)$ gives us the recurrence relation we want:
\[a_m(n) = \frac{1}{m1}\left[1 + \sum_{j=0}^{n1} \binom{n}{j}a_m(j)\right].\]In theory, we can use this to calculate $a_m(n)$ for whatever values we like. Let’s set $m=2$ and calculate the first few values of $n$ just to see how it works:
\[\begin{align*} a_2(0) & = 1,\\ a_2(1) & = 1 + a_2(0) = 2,\\ a_2(2) & = 1 + a_2(0) + 2a_2(1) = 6,\\ a_2(3) & = 1 + a_2(0) + 3a_2(1) + 3a_3(1) = 26. \end{align*}\]We can use a computer to calculate the values for large $n$. Here’s a quick and dirty implementation in Haskell:
mySumRecur :: Double > Integer > Double
mySumRecur m n = foldr step 1 [0..(n1)] * 1/(m  1)
where step j accum = accum + (choose n j)*(mySumRecur m j)
(Supply your own choose function here.) A closed form would be nice. Solving the recurrence is tricky, so we turn to a different method based on power series.
Exercise 1. Using the recurrence, show that $a_m(n)$ is rational if $m$ is rational.
Power series
Sellers’s second approach uses a power series
\[f_n(x) := \sum_{k\geq1} k^nx^k.\]Note that $a_m(n) = f_n(m^{1})$. We’re going to express $f_n(x)$ in terms of the geometric power series
\[\frac{1}{1x} = 1 + x + x^2 + \cdots = \sum_{k\geq0}x^k.\]This is analogous to seeding the recurrence for $a_m(n)$ using $a_m(0)$. We can do this by applying the operator $x\frac{d}{dx}$, which simply pulls down the power of each term. Applying $n$ times gives
\[\left(x\frac{d}{dx}\right)^n\sum_{k\geq 0}x^k = \sum_{k\geq 1}k^nx^k = f_n(x).\]To get something useful, we can apply the same operator $n$ times to $(1x)^{1}$ and evaluate at $x = m^{1}$. This is like solving the recurrence, so apply the operator in steps first:
\[\begin{align*} \left(x\frac{d}{dx}\right)\frac{1}{1x} & = \frac{x}{(1x)^2}, \\ \left(x\frac{d}{dx}\right)^2\frac{1}{1x} & = \frac{x^2+x)}{(1x)^3}, \\ \left(x\frac{d}{dx}\right)^3\frac{1}{1x} & = \frac{x^3+4x^2 +x}{(1x)^4}, \\ \left(x\frac{d}{dx}\right)^4\frac{1}{1x} & = \frac{x^4+11x^3+11x^2 +x}{(1x)^5}, \\ \left(x\frac{d}{dx}\right)^5\frac{1}{1x} & = \frac{x^5+26x^4+66x^3+26x^2 +x}{(1x)^6}. \end{align*}\]We can see the pattern is $f_n(x) = g_n(x)(1x)^{(n+1)}$ for some polynomial $g_n(x)$ of degree $n$. But what is $g_n(x)$? The sequence of polynomial coefficients in the numerator (including $0$ for the constant terms) is
\[1, 0, 1, 1, 0, 1, 4, 1, 0, 1, 11, 11, 1, 0, 1, 26, 66, 26, 1, \ldots.\]Searching the Online Encyclopedia of Integer Sequences reveals this is A123125, the sequence of Eulerian numbers. These are denoted $e(n, j)$ for $n \geq 1$ and $0 \leq j \leq n$, and defined by $e(1,1) = 1$, $e(1, k) = 0$ for $k \leq 0$, and
\[e(n, j) = je(n1,j)+(nj+1)e(n1,j1).\]So, it looks like $g_n(x)$ has coefficients given by the Eulerian numbers, but let’s check using induction. We’ve already done the base case, so assume it works for $n$ and try to prove it works for $n+1$:
\[\begin{align*} f_{n+1}(x) & = \left(x\frac{d}{dx}\right)f_n(x) \\ & = \left(x\frac{d}{dx}\right) \frac{g_n(x)}{(1x)^{n+1}} \\ & = \left(x\frac{d}{dx}\right) \frac{1}{(1x)^{n+1}} \sum_{j=0}^n e(n, j)x^j\\ & = \frac{1}{(1x)^{n+1}} \sum_{j=0}^n je(n, j)x^j + \frac{(n + 1)}{(1x)^{n+2}} \sum_{j=0}^n e(n, j)x^{j+1}\\ & = \frac{1}{(1x)^{n+2}} \sum_{j=0}^n \big[je(n, j)x^j + (n  j + 1)e(n, j)x^{j+1}\big]\\ & = \frac{1}{(1x)^{n+2}} \sum_{j=0}^{n+1} \big[je(n, j) + (n  j + 2)e(n, j1)\big]x^j\\ & = \frac{1}{(1x)^{n+2}} \sum_{j=0}^{n+1} e(n+1, j)x^j. \end{align*}\]So it’s true for all $n$.
We now have a method for computing $f_n(x)$, but another recurrence to solve, this time for $e(n, j)$. We can make our lives easier a bit easier here by using a closedform expression for $e(n, j)$,
\[4e(n, j) = \sum_{\ell=0}^j(1)^\ell (j\ell)^n\binom{n+1}{\ell}.\]Here is some more quick Haskell code for computing $f_n(x)$:
eulerian :: Integer > Integer > Double
eulerian n j = sum $ map f [0..j]
where f l = (1)^l * (choose (n + 1) l) * fromIntegral (j  l)^n
mySumPower :: Integer > Double > Double
mySumPower n x = (sum $ map f [1..n]) / (1  x)^(n + 1)
where f j = x^j * eulerian n j
Exercise 2. Prove the closedform expression for $e(n, j)$.
Polylogarithms and quantum statistics
The polylogarithm is a special function defined for complex $z<1$ by the series
\[\mbox{Li}_s(z) = \sum_{k\geq 1} \frac{z^k}{k^s}.\]It can be analytically continued to $z \geq 1$. In terms of the power series introduced earlier, $f_n(x) = \mbox{Li}_{n}(x)$. For $z < 1$, it’s possible to rewrite $\mbox{Li}_s(z)$ in a slightly surprising way:
\[\begin{align*} \frac{1}{\Gamma(s)}\int_0^\infty dx\, \frac{x^{s1}}{z^{1}e^x  1} & = \frac{1}{\Gamma(s)}\int_0^\infty dx\, ze^{x}\frac{x^{s1}}{1  ze^{x}} \\ & = \frac{1}{\Gamma(s)}\int_0^\infty dx\, x^{s1} \sum_{k\geq 1}z^{k}e^{kx} \\ & = \frac{1}{\Gamma(s)}\sum_{k\geq 1}z^{k}\int_0^\infty \frac{dy}{k}\left(\frac{y}{k}\right)^{s1} e^{y} \\ & = \frac{1}{\Gamma(s)}\sum_{k\geq 1}\frac{z^{k}}{k^s}\int_0^\infty dy\,y^{s1} e^{y} \\ & = \frac{1}{\Gamma(s)}\sum_{k\geq 1}\frac{z^{k}}{k^s}\Gamma(s) = \sum_{k\geq 1}\frac{z^{k}}{k^s} = \mbox{Li}_s(z). \end{align*}\]In other words, the polylogarithm can be represented by an integral,
\[\mbox{Li}_s(z) = \frac{1}{\Gamma(s)}\int_0^\infty dx\, \frac{x^{s1}}{z^{1}e^x  1}.\]This integral turns out to be very important in quantum statistics: it is used to find expected values for observables in the grand canonical ensemble for a free gas of bosons. In less highfalutin terms, it is used to find the average energy, pressure, density, colour, etc, of a gas of quantum particles which can be in the same state.
It would take us too far afield to discuss the subject in any detail. We’ll content ourselves with a handwaving explanation of each term in the integral:
 The integration variable $x$ is a dimensionless variable related to distribution of particle energies by $x = \beta E = E/kT$, where $k \approx 1.381 \times 10^{23} \mbox{ J K}^{1}$ is the Boltzmann constant and $T$ is the (fixed) temperature of the gas.

The factor which does not depend on $s$,
\[\frac{dx}{z^{1}e^x + 1}\]is the expected number of particles per unit volume in a “slice” of energy space at $x = \beta E$ of width $dx$.
 The power in the numerator $x^{s1}$ tells us how the observable depends on the energy.
 Finally, $z$ is a mysterious quantity called the fugacity, related to the chemical potential $\mu$ by $z = e^{\mu \beta}$. The potential $\mu$ tells us the energy cost of adding a particle to the gas at fixed entropy and volume.
So, it’s nice to see how the $n$ and $m$ can take on physical meaning. In fact, for vanilla cases (a 3D, nonrelativistic gas) the observables have halfinteger $s$ values we can’t handle with Sellers’s methods. However, for ultrarelativistic particles they will be integers and we’re in business! We also get extremely interesting behaviour in the limit $z \to 1$: BoseEinstein condensation and an associated phase transitions. For more detail on quantum statistics and applications of the polylogarithm, I refer interested readers to David Tong’s excellent notes on statistical physics, which in turn owe a great deal to Mehran Kardar’s approach.
I should also mention that the polylogarithm does dual service as an expectation integral for fermi gases. These are gases consisting of particles which are forbidden to occupy the same state. The integral we use is
\[\mbox{Li}_s(z) = \frac{1}{\Gamma(s)}\int_0^\infty dx\, \frac{x^{s1}}{z^{1}e^x + 1},\]with the terms interpreted as in the bosonic case. Once again, observables for nonrelativistic gases have halfinteger values of $s$, while ultrarelativistic gases have integer values. Although nonrelativistic fermi gases are interesting, stellar remnants like white dwarfs and neutron stars are made of ultrarelativistic fermi gases. In conclusion, our original “toy” series is not only a good opportunity to learn some summation techniques, but a good opportunity to see what this summation can tell us about the interior of a neutron star!