Chapter 47Mathematics for Statistical Mechanics — The n-Dimensional Ball, Stirling's Formula, and the Zeta Function

Derivation: Volume formula for the n-dimensional ball (Part 1 — Computing $I_n$ by factoring it into a product)

Method A: Regard it as a product of $n$ independent one-variable integrals. If we use the property of the exponential function $e^{a+b}=e^ae^b$ in the reverse direction, the sum $q_1^2+\cdots+q_n^2$ in the exponent can be broken up into a product of separate factors:

$$ e^{-(q_1^2+q_2^2+\cdots+q_n^2)}=e^{-q_1^2}\cdot e^{-q_2^2}\cdots e^{-q_n^2} $$

Since each of the integration variables $q_1,\ldots,q_n$ runs independently from $-\infty$ to $\infty$, this $n$-fold integral factors into a product of $n$ one-variable integrals of the same form:

$$ I_n=\left(\int_{-\infty}^{\infty}e^{-q_1^2}\dd q_1\right)\left(\int_{-\infty}^{\infty}e^{-q_2^2}\dd q_2\right)\cdots\left(\int_{-\infty}^{\infty}e^{-q_n^2}\dd q_n\right) $$

Every factor is the Gaussian integral $\displaystyle\int_{-\infty}^\infty e^{-x^2}\dd x=\sqrt\pi$, which has the same value regardless of the name of the integration variable (this was proved in University Mathematics, Section 7.5 in two ways: by converting the double integral to polar coordinates and by using the squeeze theorem). Hence all $n$ factors equal $\sqrt\pi$, so that

\begin{equation} I_n=(\sqrt\pi)^n=\pi^{n/2}\qquad\text{……①} \label{eq:47-In-way-a} \end{equation}

is obtained.

∎ (end of the derivation of ①)

Derivation: Volume formula for the n-dimensional ball (Part 2 — Computing $I_n$ by dividing space into “shells”)

Method B: Divide $n$-dimensional space into “thin shells” of constant distance $r$ from the origin, and integrate. The integrand $e^{-(q_1^2+\cdots+q_n^2)}$ is determined only by the distance $r=\sqrt{q_1^2+\cdots+q_n^2}$ from the origin and does not depend on direction (since $q_1^2+\cdots+q_n^2=r^2$, the integrand can be rewritten as $e^{-r^2}$). When we integrate such a function, “determined by the distance $r$ alone,” over all of space, the calculation is easier if we divide $n$-dimensional space into thin shells centered at the origin, running from radius $r$ to $r+\dd r$—thin rings in two dimensions, thin spherical shells in three dimensions—and integrate over them.

Let us consider the volume $\dd V_n(r)$ of the shell of radius $r$ (the region between $r$ and $r+\dd r$). This is nothing other than the volume of the $n$-dimensional ball of radius $r+\dd r$ minus the volume of the $n$-dimensional ball of radius $r$, so by the very definition of differentiation (the limit of a difference), we have

$$ \dd V_n(r)=V_n(r+\dd r)-V_n(r)=V_n'(r)\,\dd r $$

Substituting the result of Section 47.2.1, $V_n(r)=a_nr^n$, and differentiating, we get

\begin{equation} \dd V_n(r)=na_nr^{n-1}\,\dd r \label{eq:47-shell} \end{equation}

This is the generalization to $n$ dimensions of the relations “differentiating the area $\pi r^2$ of a circle of radius $r$ with respect to $r$ gives the circumference $2\pi r$ (because circumference × $\dd r$ is the area of a thin ring)” and “differentiating the volume $\frac43\pi r^3$ of a sphere of radius $r$ with respect to $r$ gives the surface area $4\pi r^2$ (because surface area × $\dd r$ is the volume of a thin spherical shell)” (it is worth checking that in these cases $na_nr^{n-1}=2\pi r$ for $n=2$ and $na_nr^{n-1}=4\pi r^2$ for $n=3$).

Using this “volume element of a shell,” the integral $I_n$ over all of space (Eq. \eqref{eq:47-In-def}) can be rewritten as a one-variable integral in which $r$ runs from $0$ to $\infty$:

$$ I_n=\int_0^\infty e^{-r^2}\,\dd V_n(r)=\int_0^\infty e^{-r^2}\cdot na_nr^{n-1}\,\dd r=na_n\int_0^\infty r^{n-1}e^{-r^2}\,\dd r $$

To evaluate this last integral, we perform the substitution $t=r^2$. Differentiating $t=r^2$ with respect to $r$ gives $\dd t=2r\,\dd r$, that is, $\dd r=\dfrac{\dd t}{2r}=\dfrac{\dd t}{2\sqrt t}$ (we used $r=\sqrt t$). The limits of integration correspond directly, $r=0\to t=0$ and $r\to\infty\to t\to\infty$, and $r^{n-1}=(\sqrt t)^{n-1}=t^{(n-1)/2}$, so

$$ \int_0^\infty r^{n-1}e^{-r^2}\,\dd r=\int_0^\infty t^{\frac{n-1}2}e^{-t}\cdot\frac{\dd t}{2\sqrt t}=\frac12\int_0^\infty t^{\frac{n-1}2-\frac12}e^{-t}\,\dd t=\frac12\int_0^\infty t^{\frac n2-1}e^{-t}\,\dd t $$

The last integral is exactly the definition of the gamma function, $\Gamma(\alpha)=\displaystyle\int_0^\infty x^{\alpha-1}e^{-x}\dd x$ (University Mathematics, Definition 32.3), with $\alpha=\dfrac n2$, so it equals $\Gamma\!\left(\dfrac n2\right)$. Therefore,

\begin{equation} I_n=na_n\cdot\frac12\Gamma\!\left(\frac n2\right)=\frac n2\,a_n\,\Gamma\!\left(\frac n2\right)\qquad\text{……②} \label{eq:47-In-way-b} \end{equation}

∎ (end of the derivation of ②)

Q.E.D.: Comparing ① and ② to determine $a_n$

Since we have computed the same integral $I_n$ in two different ways, ① (Eq. \eqref{eq:47-In-way-a}) and ② (Eq. \eqref{eq:47-In-way-b}) are equal:

$$ \pi^{n/2}=\frac n2\,a_n\,\Gamma\!\left(\frac n2\right) $$

Solving this for $a_n$ gives

$$ a_n=\frac{\pi^{n/2}}{\dfrac n2\,\Gamma\!\left(\dfrac n2\right)} $$

Now, substituting $\alpha=\dfrac n2$ into the recurrence relation of the gamma function, $\Gamma(\alpha+1)=\alpha\,\Gamma(\alpha)$ (University Mathematics, Formula 32.4), gives $\Gamma\!\left(\dfrac n2+1\right)=\dfrac n2\,\Gamma\!\left(\dfrac n2\right)$, so rewriting the denominator, we obtain

$$ a_n=\frac{\pi^{n/2}}{\Gamma\!\left(\dfrac n2+1\right)} $$

Since $V_n(R)=a_nR^n$, as we saw in Section 47.2.1, it follows that

$$ V_n(R)=\frac{\pi^{n/2}}{\Gamma\!\left(\dfrac n2+1\right)}R^n $$

which is what was to be shown.

∎

Derivation: Formula 47.2

Eq. \eqref{eq:47-gauss-a}. Make the substitution $u=\sqrt a\,x$. Since $\dd u=\sqrt a\,\dd x$ and the interval of integration remains $-\infty\to\infty$, we get $\displaystyle\int_{-\infty}^\infty e^{-ax^2}\dd x=\frac1{\sqrt a}\int_{-\infty}^\infty e^{-u^2}\dd u=\frac{\sqrt\pi}{\sqrt a}=\sqrt{\pi/a}$ ($\int e^{-u^2}\dd u=\sqrt\pi$ is the Gaussian integral used in Section 47.2.2).

Eq. \eqref{eq:47-gauss-b}. Make the substitution $t=ax^2$. From $x=\sqrt{t/a}$ we get $\dd x=\dfrac{\dd t}{2\sqrt{a\,t}}$, and $x^{2n}=(t/a)^n$; since the interval of integration remains $0\to\infty$,

$$ \int_0^\infty x^{2n}e^{-ax^2}\dd x=\int_0^\infty\left(\frac ta\right)^ne^{-t}\frac{\dd t}{2\sqrt a\,\sqrt t}=\frac1{2a^{\,n+\frac12}}\int_0^\infty t^{\,n-\frac12}e^{-t}\dd t=\frac{\Gamma\!\left(n+\frac12\right)}{2a^{\,n+\frac12}} $$

(the last integral is the definition of the gamma function with $\alpha=n+\frac12$). Further, applying the recurrence relation $\Gamma(\alpha+1)=\alpha\Gamma(\alpha)$ $n$ times until we reach $\Gamma(\frac12)=\sqrt\pi$, we get

$$ \Gamma\!\left(n+\tfrac12\right)=\left(n-\tfrac12\right)\left(n-\tfrac32\right)\cdots\tfrac12\,\Gamma\!\left(\tfrac12\right)=\frac{(2n-1)(2n-3)\cdots1}{2^n}\sqrt\pi=\frac{(2n-1)!!}{2^n}\sqrt\pi $$

Therefore, $\dfrac{\Gamma(n+\frac12)}{2a^{n+\frac12}}=\dfrac{(2n-1)!!}{2^{n+1}}\dfrac{\sqrt\pi}{a^{\,n+\frac12}}=\dfrac{(2n-1)!!}{2^{n+1}}\sqrt{\dfrac\pi{a^{2n+1}}}$ is obtained.

∎

Derivation: Deriving Formula 47.3 by the method of Riemann sums

Taking the natural logarithm of both sides of $n!=1\cdot2\cdot3\cdots n$, the product becomes a sum by the property of logarithms $\ln(ab)=\ln a+\ln b$:

$$ \ln n!=\ln(1\cdot2\cdot3\cdots n)=\ln1+\ln2+\cdots+\ln n=\sum_{k=1}^n\ln k $$

As an approximation valid when $n$ is large, we replace this sum $\displaystyle\sum_{k=1}^n\ln k$ by the integral $\displaystyle\int_1^n\ln x\,\dd x$. The reason is as follows. On the graph of the function $y=\ln x$, tile each interval $[k-1,k]$ with a rectangle of height $\ln k$ (the value of the function at the right endpoint) and width $1$; the total area of these rectangles is nothing other than $\displaystyle\sum_{k=1}^n\ln k$. Since $\ln x$ is monotonically increasing, this total area approximates the area under the curve, $\displaystyle\int_1^n\ln x\,\dd x$, from above, and the difference can be regarded as small compared with the integral itself (Figure 47.3. This is an application of the method of Riemann sums that you learned in high school mathematics—the relation between $\Sigma$ and $\displaystyle\int$. We will examine what the difference actually is later).

To evaluate the integral $\displaystyle\int_1^n\ln x\,\dd x$, we use the indefinite integral formula $\displaystyle\int\ln x\,\dd x=x\ln x-x+C$ (obtained by integration by parts with $u=\ln x,\ \dd v=\dd x$; since $\dd u=\dd x/x,\ v=x$, we have $\int\ln x\,\dd x=x\ln x-\int x\cdot\frac1x\,\dd x=x\ln x-\int\dd x=x\ln x-x+C$). Then

$$ \int_1^n\ln x\,\dd x=\Big[x\ln x-x\Big]_1^n=(n\ln n-n)-(1\cdot\ln1-1) $$

Since $\ln1=0$, the expression in the second pair of parentheses on the right-hand side is $0-1=-1$. Therefore,

$$ \int_1^n\ln x\,\dd x=(n\ln n-n)-(-1)=n\ln n-n+1 $$

Approximating the sum by the value of this integral, we obtain

$$ \sum_{k=1}^n\ln k\approx n\ln n-n+1 $$

When $n\gg1$ ($n$ is sufficiently large), the constant $+1$ at the end is negligibly small compared with $n\ln n$ and $n$, so we may drop it:

$$ \sum_{k=1}^n\ln k\approx n\ln n-n $$

We combine the right-hand side into a single logarithm so that it will later be easy to bring down from the exponent of the exponential function. Using $n=n\ln e$ (because $\ln e=1$) and the properties of logarithms (the constant-multiple formula $k\ln a=\ln a^k$ and the subtraction formula $\ln a-\ln b=\ln\frac ab$), we have

$$ n\ln n-n=n\ln n-n\ln e=\ln n^n-\ln e^n=\ln\frac{n^n}{e^n}=\ln(n^ne^{-n}) $$

Thus we have obtained $\ln n!\approx n\ln n-n=\ln(n^ne^{-n})$. Rewriting this in the form $n!\approx n^ne^{-n}$ gives Formula 47.3. However, the meaning of “$\approx$” requires care: what we have shown here is an approximation after taking logarithms (even if $\ln A\approx\ln B$, the ratio $A/B$ is not necessarily close to $1$). In fact, the ratio of $n!$ to $n^ne^{-n}$ does not approach $1$; as we will see later, the former is about $\sqrt{2\pi n}$ times as large as the latter. From now on, writing $n!\approx n^ne^{-n}$ is shorthand that presupposes “taking the logarithm before using it.”

Supplement: The actual difference is $\ln n!-(n\ln n-n)\approx\dfrac12\ln(2\pi n)$ (obtained from the refined version in University Mathematics, Section 32.3). It is about $2.07$ for $n=10$ and about $3.22$ for $n=100$; even as $n$ increases, it grows only slowly, like a logarithm, and remains small compared with the leading term $n\ln n$. (The difference from “$n\ln n-n+1$,” which includes the constant $1$ dropped in the integral approximation, is $\frac12\ln(2\pi n)-1$.)

∎

Derivation: Proof that the range of convergence of the zeta function is $s\gt1$

Let $f(x)=1/x^s$ ($s\gt0$). Then $f(x)$ is monotonically decreasing for $x\gt0$ (the denominator $x^s$ grows with $x$, so $f(x)$ shrinks). A monotonically decreasing function takes its maximum value at the left endpoint $x=k$ and its minimum value at the right endpoint $x=k+1$ of the interval $[k,k+1]$, so

$$ f(k+1)\le f(x)\le f(k)\qquad(k\le x\le k+1) $$

holds. Integrating this inequality over the interval $[k,k+1]$ (of width $1$; an inequality keeps its direction when integrated), we obtain

$$ f(k+1)\le\int_k^{k+1}f(x)\,\dd x\le f(k) $$

(Figure 47.5. The equality sign on the left holds only when $f$ is a constant function, so, since $f$ is strictly monotonically decreasing here, the inequalities are in fact strict.)

We add this inequality up for $k=1,2,\ldots,n-1$. Adding up the left-hand sides gives $\displaystyle\sum_{k=1}^{n-1}f(k+1)=\sum_{j=2}^n f(j)=\sum_{k=2}^n\frac1{k^s}$; adding up the middle terms gives $\displaystyle\sum_{k=1}^{n-1}\int_k^{k+1}f(x)\,\dd x=\int_1^n f(x)\,\dd x$ (we have merely joined the intervals of integration); and adding up the right-hand sides gives $\displaystyle\sum_{k=1}^{n-1}f(k)=\sum_{k=1}^{n-1}\frac1{k^s}$. Hence

\begin{equation} \sum_{k=2}^n\frac1{k^s}\lt\int_1^n f(x)\,\dd x\lt\sum_{k=1}^{n-1}\frac1{k^s} \label{eq:47-zeta-bounds1} \end{equation}

holds. Adding $f(1)=1$ to both sides of the left inequality gives $\displaystyle\sum_{k=1}^n\frac1{k^s}\lt\int_1^nf(x)\,\dd x+1$, and adding $f(n)=1/n^s$ to both sides of the right inequality gives $\displaystyle\int_1^nf(x)\,\dd x+\frac1{n^s}\lt\sum_{k=1}^n\frac1{k^s}$ (note that, by choosing what to add to which side, both can be brought to the form $\displaystyle\sum_{k=1}^n\frac1{k^s}$). Combining these, we obtain

\begin{equation} \int_1^nf(x)\,\dd x+\frac1{n^s}\lt\sum_{k=1}^n\frac1{k^s}\lt\int_1^nf(x)\,\dd x+1 \label{eq:47-zeta-bounds2} \end{equation}

Next, we actually compute this integral. For $s\ne1$, by the integration formula for power functions,

$$ \int_1^nf(x)\,\dd x=\int_1^nx^{-s}\,\dd x=\left[\frac{x^{1-s}}{1-s}\right]_1^n=\frac{n^{1-s}-1}{1-s}=\frac1{s-1}\left(1-\frac1{n^{s-1}}\right) $$

(The last equality merely tidies up the sign by multiplying the numerator and denominator by $-1$. We also used $n^{1-s}=1/n^{s-1}$.) Substituting this into Eq. \eqref{eq:47-zeta-bounds2}, we obtain

\begin{equation} \frac1{s-1}\left(1-\frac1{n^{s-1}}\right)+\frac1{n^s}\lt\sum_{k=1}^n\frac1{k^s}\lt\frac1{s-1}\left(1-\frac1{n^{s-1}}\right)+1 \label{eq:47-zeta-bounds3} \end{equation}

What remains is to use this inequality to examine what happens to the partial sum $\displaystyle\sum_{k=1}^n1/k^s$ as $n\to\infty$, for each range of $s$.

[Convergence] The case $s\gt1$. In the right-hand inequality of Eq. \eqref{eq:47-zeta-bounds3}, if $s\gt1$ then $1/n^{s-1}\gt0$, so, regardless of $n$,

$$ \sum_{k=1}^n\frac1{k^s}\lt\frac1{s-1}\left(1-\frac1{n^{s-1}}\right)+1\lt\frac1{s-1}+1 $$

holds, and the partial sums have an upper bound $\dfrac1{s-1}+1$ (a fixed wall independent of $n$). On the other hand, each term $1/k^s$ is positive, so the partial sum increases monotonically each time we increase $n$. If an ever-increasing sequence has a wall that it never exceeds, the sequence has no choice but to approach some value (a sequence that grew without limit would eventually exceed any wall). This is the basic theorem, treated rigorously in university analysis, that “a monotonically increasing sequence that is bounded above converges”; in this book we stay with this intuitive explanation. Hence, for $s\gt1$, $\zeta(s)$ converges.

[Divergence] The case $s\le1$. We use the bound from below (the left-hand inequality of Eq. \eqref{eq:47-zeta-bounds2}) $\displaystyle\sum_{k=1}^n\frac1{k^s}\gt\int_1^nx^{-s}\dd x$.

• The case $0\lt s\lt1$: $\displaystyle\int_1^nx^{-s}\dd x=\frac{n^{1-s}-1}{1-s}$, and since $1-s\gt0$, we have $n^{1-s}\to\infty$. Hence the partial sums diverge to $\infty$.
• The case $s=1$: since $\displaystyle\int_1^n\frac{\dd x}x=\ln n\to\infty$, the partial sums (the harmonic series) diverge.
• The case $s\le0$: since $1/k^s=k^{-s}\ge1$, the partial sum, being a sum of $n$ terms, is at least $n$, and diverges.

Therefore, $\zeta(s)$ converges only when $s\gt1$.

(Supplement) Listing the conditions under which $1/n^{s-1}$ and $1/n^s$, which appear in Eq. \eqref{eq:47-zeta-bounds3}, approach $0$ as $n\to\infty$, we have $s\ne1$ (the integration formula can be used), $s\gt0$ ($1/n^s\to0$), and $s\gt1$ ($1/n^{s-1}\to0$); the range satisfying all three simultaneously is $s\gt1$. This agrees with the conclusion of the case analysis above, and confirms that the condition $s\gt1$ emerges in the same way from both the “wall for convergence” and the “behavior of each term in the bounding expressions.”

∎

Derivation: $\zeta(4)=\pi^4/90$ (Fourier series and Parseval's identity)

By Example 28.3, the Fourier series of $f(x)=x^2\ (-\pi\lt x\le\pi)$ is

$$ x^2=a_0+\sum_{k=1}^\infty a_k\cos kx,\qquad a_0=\frac{\pi^2}3,\quad a_k=\frac{4(-1)^k}{k^2} $$

This series satisfies $\lvert a_k\cos kx\rvert\le4/k^2$, and $\sum4/k^2$ converges (Theorem 47.1, $s=2$), so the series converges uniformly and may be integrated term by term. We therefore square both sides and integrate from $-\pi$ to $\pi$.

The left-hand side is $\displaystyle\int_{-\pi}^{\pi}x^4\dd x=\left[\frac{x^5}5\right]_{-\pi}^{\pi}=\frac{2\pi^5}5$. Squaring the right-hand side and expanding gives terms of the forms $a_0^2$, $2a_0a_k\cos kx$, and $a_ja_k\cos jx\cos kx$. Since $\int_{-\pi}^{\pi}\cos kx\,\dd x=0$ (Formula 27.2 in University Mathematics, Section 27.1), the terms $2a_0a_k\cos kx$ integrate to $0$, and the terms $\cos jx\cos kx$ with $j\ne k$ also integrate to $0$ by orthogonality (Formula 27.4). All that remain are $\displaystyle\int_{-\pi}^{\pi}a_0^2\dd x=2\pi a_0^2$ and $\displaystyle\int_{-\pi}^{\pi}a_k^2\cos^2kx\,\dd x=\pi a_k^2$ (Formula 27.4 with $m=n$), so

$$ \frac{2\pi^5}5=2\pi a_0^2+\pi\sum_{k=1}^\infty a_k^2\qquad\text{(Parseval's identity)} $$

holds. Substituting $a_0=\pi^2/3$ and $a_k^2=16/k^4$, we obtain

$$ \frac{2\pi^5}5=2\pi\cdot\frac{\pi^4}9+16\pi\sum_{k=1}^\infty\frac1{k^4}=\frac{2\pi^5}9+16\pi\,\zeta(4) $$

Solving this for $\zeta(4)$: from $16\pi\,\zeta(4)=2\pi^5\left(\dfrac15-\dfrac19\right)=2\pi^5\cdot\dfrac4{45}=\dfrac{8\pi^5}{45}$ we have

$$ \zeta(4)=\frac{8\pi^5}{45\cdot16\pi}=\frac{\pi^4}{90} $$

Applying the same procedure to $f(x)=x$ (Example 28.1: $b_k=2(-1)^{k+1}/k$), we find that $\zeta(2)=\pi^2/6$ also follows from $\displaystyle\frac{2\pi^3}3=\pi\sum_{k=1}^\infty\frac4{k^2}=4\pi\,\zeta(2)$.

∎

Derivation: Theorem 47.2 (geometric series expansion + term-by-term integration)

We rewrite the denominator of the integrand using the formula for the geometric series. Since $e^{-x}\lt1$ (when $x\gt0$), applying the sum formula for an infinite geometric series with common ratio $e^{-x}$, $\displaystyle\frac1{1-r}=\sum_{k=0}^\infty r^k$ ($|r|\lt1$), with $r=e^{-x}$, we obtain

$$ \frac1{e^x-1}=\frac{e^{-x}}{1-e^{-x}}=e^{-x}\sum_{k=0}^\infty e^{-kx}=\sum_{k=0}^\infty e^{-(k+1)x}=\sum_{k=1}^\infty e^{-kx} $$

(The last equality merely shifts the summation index, $k\to k-1$.) We substitute this into the integral and interchange the order of the series sum and the integral (since each term $x^{s-1}e^{-kx}$ is nonnegative for $x\gt0$, the order of summation and integration may be interchanged. This is justified by a theorem of university-level analysis called the monotone convergence theorem. Here we accept it and move on):

$$ \int_0^\infty\frac{x^{s-1}}{e^x-1}\,\dd x=\int_0^\infty x^{s-1}\sum_{k=1}^\infty e^{-kx}\,\dd x=\sum_{k=1}^\infty\int_0^\infty x^{s-1}e^{-kx}\,\dd x $$

The integral of each term can be reduced to the gamma function by the substitution $u=kx$ ($\dd u=k\,\dd x$) (the same calculation as we learned in University Mathematics, Example 32.4):

$$ \int_0^\infty x^{s-1}e^{-kx}\,\dd x=\int_0^\infty\left(\frac uk\right)^{s-1}e^{-u}\,\frac{\dd u}k=\frac1{k^s}\int_0^\infty u^{s-1}e^{-u}\,\dd u=\frac{\Gamma(s)}{k^s} $$

Substituting this, we obtain

$$ \int_0^\infty\frac{x^{s-1}}{e^x-1}\,\dd x=\sum_{k=1}^\infty\frac{\Gamma(s)}{k^s}=\Gamma(s)\sum_{k=1}^\infty\frac1{k^s}=\Gamma(s)\,\zeta(s) $$

For the series $\sum1/k^s$ to converge, $s\gt1$ is required (Theorem 47.1) (this also corresponds to the fact that, as $x\to0$, the integrand behaves like $x^{s-2}$, and the integral diverges unless $s\gt1$).

∎

Derivation: Formula 47.5 (by integration by parts)

Using the relation in the box above, we apply integration by parts with $u=x^p,\ \dd v=-\dfrac{\dd}{\dd x}\!\left(\dfrac1{e^x-1}\right)\dd x$ (that is, $v=-\dfrac1{e^x-1}$, $\dd u=px^{p-1}\dd x$), which gives

$$ \int_0^\infty\frac{x^pe^x}{(e^x-1)^2}\,\dd x=\left[-\frac{x^p}{e^x-1}\right]_0^\infty+p\int_0^\infty\frac{x^{p-1}}{e^x-1}\,\dd x $$

We examine the boundary term. As $x\to\infty$, $x^p/(e^x-1)\to0$ (the exponential function approaches $0$ faster than any power function, for the same reason as in the Caution in University Mathematics, Section 32.2). As $x\to0$, $e^x-1\approx x$ (we use the Taylor expansion of $e^x$ up to first order), so $x^p/(e^x-1)\approx x^{p-1}\to0$ (when $p\gt1$). Hence, if $p\ge2$ is a natural number, the boundary term becomes $0$ completely, and we obtain

$$ \int_0^\infty\frac{x^pe^x}{(e^x-1)^2}\,\dd x=p\int_0^\infty\frac{x^{p-1}}{e^x-1}\,\dd x $$

The integral on the right-hand side is exactly the form obtained from Formula 47.4 (Eq. \eqref{eq:47-bose-i}) by replacing $p$ with $p-1$, so it equals $\Gamma(p)\zeta(p)$. Therefore,

$$ \int_0^\infty\frac{x^pe^x}{(e^x-1)^2}\,\dd x=p\,\Gamma(p)\,\zeta(p)=\Gamma(p+1)\,\zeta(p)=p!\,\zeta(p) $$

(The second-to-last equality uses the recurrence relation $\Gamma(p+1)=p\Gamma(p)$, and the last equality uses $\Gamma(p+1)=p!$ for natural numbers.)

∎

Derivation: Theorem 47.3

As in the proof of Theorem 47.2, we expand the denominator of the integrand into a series. For $\dfrac1{e^x+1}=\dfrac{e^{-x}}{1+e^{-x}}$, using the infinite geometric series with common ratio $-e^{-x}$ ($|-e^{-x}|\lt1$), $\dfrac1{1+r}=\dfrac1{1-(-r)}=\displaystyle\sum_{k=0}^\infty(-r)^k$, we obtain

$$ \frac1{e^x+1}=e^{-x}\sum_{k=0}^\infty(-1)^ke^{-kx}=\sum_{k=0}^\infty(-1)^ke^{-(k+1)x}=\sum_{k=1}^\infty(-1)^{k-1}e^{-kx} $$

This is a series in which the signs alternate (an alternating series). Integrating this and applying the same procedure as in the proof of Theorem 47.2 (term-by-term integration + the integral formula for the gamma function), we obtain

$$ \int_0^\infty\frac{x^{s-1}}{e^x+1}\,\dd x=\sum_{k=1}^\infty(-1)^{k-1}\int_0^\infty x^{s-1}e^{-kx}\,\dd x=\sum_{k=1}^\infty\frac{(-1)^{k-1}\Gamma(s)}{k^s}=\Gamma(s)\sum_{k=1}^\infty\frac{(-1)^{k-1}}{k^s} $$

(The term-by-term integration can be justified by the fact that, for $s\gt1$, $\sum_k\Gamma(s)/k^s$ converges, that is, the sum of the absolute values of the terms is finite (Theorem 47.1). For $0\lt s\le1$ the alternating series is only conditionally convergent, so this argument cannot be used.) Let us express the alternating series $\displaystyle\sum_{k=1}^\infty\frac{(-1)^{k-1}}{k^s}$ that appeared at the end in terms of the zeta function $\zeta(s)$. For $s\gt1$, the sum of the odd-numbered terms and the sum of the even-numbered terms each converge, so we may compute them separately. Since the odd-numbered terms carry a plus sign and the even-numbered terms a minus sign,

$$ \sum_{k=1}^\infty\frac{(-1)^{k-1}}{k^s}=\left(\sum_{k:\,\text{odd}}\frac1{k^s}\right)-\left(\sum_{k:\,\text{even}}\frac1{k^s}\right)=\underbrace{\left(\sum_{\text{odd}}+\sum_{\text{even}}\right)}_{=\,\zeta(s)}-2\sum_{\text{even}}\frac1{k^s} $$

(The sum of the odd terms is obtained by subtracting the sum of the even terms from the sum of all the terms, so we rewrote it as “sum over odd terms” $=\zeta(s)-$ “sum over even terms” and substituted.) The sum over the even terms alone is, letting $k=2m$ ($m=1,2,3,\ldots$),

$$ \sum_{\text{even}}\frac1{k^s}=\sum_{m=1}^\infty\frac1{(2m)^s}=\frac1{2^s}\sum_{m=1}^\infty\frac1{m^s}=\frac{\zeta(s)}{2^s} $$

and therefore we obtain

$$ \sum_{k=1}^\infty\frac{(-1)^{k-1}}{k^s}=\zeta(s)-2\cdot\frac{\zeta(s)}{2^s}=\left(1-\frac2{2^s}\right)\zeta(s)=\left(1-2^{1-s}\right)\zeta(s) $$

Substituting this, we obtain

$$ \int_0^\infty\frac{x^{s-1}}{e^x+1}\,\dd x=\Gamma(s)\left(1-2^{1-s}\right)\zeta(s) $$

This proves Theorem 47.3.

∎