7.3 The Statistical Toolkit: Densities of States, Polylogarithms, and the Bose and Fermi Integrals#

Elementary Computational Physics
Volume VII — Quantum Statistical Mechanics Notebook 7.3
Before the physics, the vocabulary. We count the quantum states of a box — literally, as lattice points — and watch the count become the density of states that converts every sum into an integral. We build the polylogarithm, the master special function of quantum gases, and evaluate the Bose and Fermi integrals three independent ways. And we meet a bounded function with a vertical tangent, whose refusal to grow past 2.612 will, three movements from now, condense a gas.
Level · advanced   •   Est. · 180–220 min
Raymond Amador v1.4.0  ·  2026-07-31  ·  CC BY 4.0 (text) / MIT (code)

Notebook overview#

This is the third and final notebook of Movement 0, and it does what the last notebook of a toolkit movement should: it assembles the specific working objects that every later notebook will consume, each one built honestly and checked against something independent. The two notebooks before it were about principles — rigidity, causality, saddle points. This one is about equipment.

Four pieces of equipment, in order of use. First, the density of states: the thermodynamics of a macroscopic system is a sum over its quantum states, and the master move of this volume converts that sum into an integral, \(\sum_{\text{states}} \to \int g(\varepsilon)\,d\varepsilon\). We refuse to take \(g(\varepsilon)\) as a formula: the particle-in-a-box modes of §6.10 are integer lattice points in an octant of \(k\)-space, and we count them — literally, with numpy.meshgrid and a boolean sum — and watch the staircase \(N(E)\) converge to the continuum sphere-volume formula, with the deficit shrinking as \(E^{-1/2}\) and carrying a name (Weyl’s surface correction) rather than an apology. The count is then repeated across dimensions and dispersions, because the exponent of \(\varepsilon\) in \(g\) is where geometry enters thermodynamics — and two of the volume’s later dramas hinge on it. Second, the Gamma–zeta–polylogarithm family, the special-function vocabulary of quantum gases, with the polylogarithm \(\mathrm{Li}_s(z)\) implemented from scratch and validated three independent ways to eight digits. Third, the Bose and Fermi integrals, evaluated three ways as well — direct quadrature, series reduction to polylogarithms, and closed forms falling to the residue-derived zeta values of §7.2 — including a number met three movements before its fame: \(\int_0^\infty x^3/(e^x-1)\,dx = \pi^4/15\), the arithmetic heart of the Stefan–Boltzmann law. Fourth, the Sommerfeld expansion, the low-temperature instrument of every metal calculation, derived and then held to account: its truncation error is measured and scales as \(T^4\), as advertised.

Two bombs are planted for later movements, deliberately and in plain sight. The two-dimensional density of states is constant — verified by counting — and that innocuous fact is why an ideal gas refuses to Bose-condense in two dimensions (§7.17). And the Bose function \(g_{3/2}(z)\) is monotonically increasing on the physical fugacity range yet bounded by \(\zeta(3/2) = 2.612\dots\), a ceiling it approaches with a vertical tangent: when the physics of §7.17 demands more than \(2.612\), the system will have to do something drastic, and the drastic thing is Bose–Einstein condensation. The toolkit knows things the physics hasn’t said yet.

Conventions (this notebook). Box energies are measured in units of the ground scale \(E_1 = \hbar^2\pi^2/2mL^2\), so a mode \(\mathbf{n} \in \mathbb{Z}_+^d\) has \(\tilde E = |\mathbf{n}|^2\); hard-wall (Dirichlet) counting and periodic counting differ only in surface terms and give the same bulk density of states, so we count the box and speak for both. Spin enters as a multiplicative degeneracy factor \(g_\sigma\) and is set to \(1\) here. The bosonic fugacity is restricted to \(0 < z \le 1\) — previewed here as the domain on which the series make sense, justified physically in §7.7. Polylogarithm partial sums always state their \(k_{\max}\) and a tail bound; near \(z = 1\) the series crawls, and the integral route or mpmath.polylog is the reference instead — a stated rule, not an improvisation.

How to read the checks. Each exercise closes with a validate call against an independent fact: the counted staircases converging to the continuum formulas in \(d = 1, 2, 3\) with the \(E^{-1/2}\) Weyl deficit; the polylogarithm agreeing three ways (series, mpmath.polylog, the defining integral via scipy.integrate.quad) to eight digits; the Bose and Fermi integrals closing onto \(\pi^2/6\), \(\pi^4/15\), \(\Gamma(3/2)\zeta(3/2)\), and \(\pi^2/12\) — the residue zetas of §7.2 doing the arithmetic; the Sommerfeld remainder shrinking sixteenfold when \(T\) halves; and the bounded, vertically-tangent approach of \(g_{3/2}\) to its ceiling. A ✓ is strong evidence; a ✗ is a prompt to locate the discrepancy.

Scope. The toolkit, not yet the physics: no ensemble is invoked, no temperature is attached to a density matrix (that is §7.4), and the occupations these integrals will average arrive in §7.7. See Pathria & Beale, Statistical Mechanics (Appendices D–E); Ashcroft & Mermin, Solid State Physics (Ch. 2, the Sommerfeld expansion); Arfken, Weber & Harris (special functions). Cross-reference §6.10/§6.16 (the box modes being counted), §7.1 (the keyhole anatomy of these integrands), §7.2 (the residue zetas; \(\zeta(3/2)\) flagged), §5.3 (the classical large-\(N\) toolkit this quantum one parallels), and forward to §7.4, §7.7, §7.10, §7.14, §7.17.

Theory in brief#

From sums to integrals: the density of states#

Everything thermodynamic in this volume will begin as a sum over quantum states, and the first tool converts that sum into an integral,

(683)#\[\sum_{\text{states}} (\cdots) \;\longrightarrow\; \int_0^\infty g(\varepsilon)\,(\cdots)\, d\varepsilon, \qquad g(\varepsilon) = \frac{dN}{d\varepsilon},\]

where \(N(E)\) counts the states with energy below \(E\). For a particle in a hard-wall box of side \(L\) (§6.10), the allowed modes are \(\mathbf{k} = \pi\mathbf{n}/L\) with \(\mathbf{n}\) a positive-integer vector: lattice points in an octant of \(k\)-space. In ground-scale units (\(\tilde E = E/E_1\), \(E_1 = \hbar^2\pi^2/2mL^2\)) a mode has \(\tilde E = |\mathbf{n}|^2\), so \(N(E)\) is the number of lattice points inside the sphere octant of radius \(\sqrt{\tilde E}\) — asymptotically its volume, \(N = \tfrac{1}{8}\cdot\tfrac{4\pi}{3}\tilde E^{3/2} = \tfrac{\pi}{6}\tilde E^{3/2}\), whence \(g(\varepsilon) \propto \sqrt{\varepsilon}\) in three dimensions. The count converges to the volume from below, with a deficit that shrinks as \(E^{-1/2}\): that deficit is the surface (Weyl) correction — the boundary’s share of the count, a celebrated piece of spectral geometry (“can one hear the shape of a drum?”), not a numerical error. Spin multiplies \(g\) by a degeneracy factor \(g_\sigma\); we set it to one and say so.

Dimensions and dispersions#

The same count in \(d\) dimensions, for the two dispersion relations physics keeps handing us,

(684)#\[\varepsilon \propto k^2:\quad g(\varepsilon) \propto \varepsilon^{d/2-1}, \qquad\qquad \varepsilon \propto k:\quad g(\varepsilon) \propto \varepsilon^{d-1},\]

from the volume of the \(k\)-space ball of radius \(k(\varepsilon)\). For massive particles: \(\sqrt{\varepsilon}\) in 3D, constant in 2D, \(\varepsilon^{-1/2}\) in 1D (an integrable divergence at the band bottom). For massless, linear dispersion (photons, phonons): \(g \propto \varepsilon^{d-1}\), in 3D the famous \(\varepsilon^2\). Two forward flags, with weight. The constant 2D density of states is the fact behind the absence of ideal-gas Bose–Einstein condensation in two dimensions (§7.17). And the \(\varepsilon^2\) photon density of states is the engine of Planck’s law (§7.14). The density of states is where geometry enters thermodynamics.

The Gamma–zeta–polylogarithm family#

The volume’s special-function vocabulary is one family. \(\Gamma(s)\) is the continued factorial (\(\Gamma(1/2) = \sqrt\pi\)); \(\zeta(s)\) collects the residue-derived even values of §7.2 (\(\zeta(2) = \pi^2/6\), \(\zeta(4) = \pi^4/90\)) and the computed odd ones (\(\zeta(3/2) = 2.612\dots\), \(\zeta(3) = 1.202\dots\)); and the polylogarithm

(685)#\[\mathrm{Li}_s(z) = \sum_{k=1}^{\infty} \frac{z^k}{k^s}, \qquad z\,\frac{d\,\mathrm{Li}_s}{dz} = \mathrm{Li}_{s-1}(z), \qquad \mathrm{Li}_s(1) = \zeta(s),\]

is the master function of quantum gases. The derivative identity follows from the series by inspection (differentiating brings down a \(k\), lowering the index) and is the engine of the vertical tangent at this notebook’s climax. The series converges for \(|z| \le 1\) (for \(s > 1\)), but slowly near \(z = 1\) — at \(z = 1\), \(s = 3/2\) the remainder after \(2\times10^6\) terms is still \(\sim 1.4\times10^{-3}\) — so the working rule is: away from \(z = 1\), partial sums with a bounded tail; near \(z = 1\), the defining integral or mpmath.polylog as reference.

The Bose and Fermi integrals#

The two integral families the whole volume runs on:

(686)#\[\int_0^\infty \frac{x^{s-1}}{z^{-1}e^{x} - 1}\,dx = \Gamma(s)\,\mathrm{Li}_s(z) \equiv \Gamma(s)\,g_s(z), \qquad \int_0^\infty \frac{x^{s-1}}{z^{-1}e^{x} + 1}\,dx = -\Gamma(s)\,\mathrm{Li}_s(-z) \equiv \Gamma(s)\,f_s(z).\]

The derivation is one geometric expansion: \(1/(z^{-1}e^x - 1) = \sum_{k\ge1} z^k e^{-kx}\), integrated term by term against \(x^{s-1}\) (each term a scaled \(\Gamma(s)\)), and the alternating version gives the Fermi family. The integrand’s anatomy — \(x^{s-1}\) against one Boltzmann-like denominator — is the keyhole integrand of §7.1, one Boltzmann factor richer, exactly as promised there. At \(z = 1\) the closed forms fall to the residue zetas of §7.2: \(\int x/(e^x-1) = \zeta(2) = \pi^2/6\), \(\int x^3/(e^x-1) = \Gamma(4)\zeta(4) = \pi^4/15\) (Stefan–Boltzmann’s number, met three movements early), \(\int \sqrt x/(e^x-1) = \Gamma(3/2)\zeta(3/2)\) (the integral BEC’s critical temperature is made of), while the Fermi side closes through the eta function \(\eta(s) = (1 - 2^{1-s})\,\zeta(s)\): \(\int x/(e^x+1) = \eta(2) = \pi^2/12\). The factor \((1-2^{1-s})\) is honest arithmetic — the alternating sum subtracts the even terms twice.

The Sommerfeld expansion#

At low temperature the Fermi function \(n_F(\varepsilon) = 1/(e^{(\varepsilon-\mu)/k_BT}+1)\) is a smeared step, and \(-\partial n_F/\partial\varepsilon\) is a peak of width \(\sim k_BT\) at \(\mu\). Integrating by parts and expanding a smooth \(H(\varepsilon)\) about \(\mu\),

(687)#\[\int_0^\infty H(\varepsilon)\,n_F(\varepsilon)\,d\varepsilon \;=\; \int_0^\mu H(\varepsilon)\,d\varepsilon \;+\; \frac{\pi^2}{6}\,(k_BT)^2\,H'(\mu) \;+\; \mathcal{O}(T^4),\]

where the coefficient is set by the second moment of the thermal peak — which is, up to bookkeeping, exactly the \(\pi^2/12\) Fermi integral above. The \(\pi^2/6\) is \(\zeta(2)\), Basel’s number, and it is bound for the electronic heat capacity of metals (§7.10). We derive the expansion and then measure its remainder: on \(H = \sqrt\varepsilon\) the truncation error shrinks sixteenfold when \(T\) halves — the \(T^4\) scaling is a certificate, not a promise.

A bounded function with a vertical tangent#

The quiet climax. On the physical fugacity range \(0 < z \le 1\) the Bose function \(g_{3/2}(z) = \mathrm{Li}_{3/2}(z)\) is monotonically increasing — and bounded:

(688)#\[g_{3/2}(z) \;\le\; g_{3/2}(1) = \zeta(3/2) = 2.612\dots, \qquad \frac{d g_{3/2}}{dz} = \frac{\mathrm{Li}_{1/2}(z)}{z} \;\xrightarrow{\;z\to1^-\;}\; \infty .\]

A finite ceiling, approached with infinite slope: the derivative identity lowers the index to \(\mathrm{Li}_{1/2}\), which diverges at \(z = 1\) (its terms \(1/\sqrt k\) are not summable). We verify both facts numerically — the ceiling respected at \(z = 0.999\), the slope already \(\approx 55\) and climbing. The puzzle to carry forward: a physical quantity proportional to \(g_{3/2}(z)\) — it will be the excited-state occupancy of a Bose gas — cannot be pushed past \(2.612\) by any fugacity. When the physics demands more, the system must do something drastic, and §7.17 names it: Bose–Einstein condensation.

Setup#

Conventions, colours, and one given power law: the density-of-states shape \(g \propto \varepsilon^{d/2-1}\) (massive) or \(\varepsilon^{d-1}\) (massless), which is Eq. 684 transcribed for the comparison plot — the exponents are derived in Exercise 2, and the helper only draws them. Everything this notebook is about you build in the exercise where it is earned: the lattice count count_states (Exercise 1), the polylogarithm polylog and the Bose quadrature bose_integral (Exercise 3), its fermionic twin fermi_integral (Exercise 5), and the two-term sommerfeld estimate (Exercise 6).

The Setup below holds this notebook’s data and instruments — nothing you are asked to build. It is collapsed so the building stays yours; expand it whenever you want the details.

Hide code cell source

import matplotlib.pyplot as plt
import mpmath
import numpy as np
from scipy.integrate import quad
from scipy.special import gamma as gamma_fn
from scipy.special import zeta

from ecp import draw, validate

ACCENT, INK, SOFT = draw.ACCENT, draw.INK, draw.SOFT
RED = "#c1121f"

# Conventions: box energies in ground-scale units (a mode n has E = |n|^2);
# hard-wall counting (positive-integer n) — periodic counting differs only in surface
# terms and shares the bulk density of states; spin degeneracy factor set to 1; bosonic
# fugacity restricted to 0 < z ≤ 1 (previewed here, justified in §7.7); polylogarithm
# partial sums always carry a stated k_max and tail bound, and NEAR z = 1 the reference
# is the defining integral or mpmath.polylog, never brute summation.


# data: eq-dos-dimensions transcribed — the two given power laws, drawn for shape
# comparison. The exponents are Exercise 2's derivation; nothing here is machinery.
def dos_shape(eps, d, dispersion="massive"):
    """The density-of-states shape g(ε) up to normalization (eq-dos-dimensions).

    Massive dispersion ε ∝ k^2 gives g ∝ ε^(d/2 − 1); massless/linear ε ∝ k gives
    g ∝ ε^(d − 1) — the exponents carry the physics (constant in 2D massive, ε^2 for 3D
    photons) while the prefactors carry only units, so this helper returns the pure
    power law for shape comparisons and plots.

    Parameters
    ----------
    eps : numpy.ndarray
        Energies (positive).
    d : int
        Spatial dimension.
    dispersion : str, optional
        'massive' (ε ∝ k^2, default) or 'massless' (ε ∝ k).

    Returns
    -------
    numpy.ndarray
        The unnormalized g(ε).
    """
    if dispersion == "massive":
        return eps ** (d / 2.0 - 1.0)
    if dispersion == "massless":
        return eps ** (d - 1.0)
    raise ValueError("dispersion must be 'massive' or 'massless'")

Exercise 1 — Counting the states of a box#

The density of states is a count before it is a formula — so we count, and we watch the formula emerge from the count. Cite Eq. 683.

  1. Derive \(N(E)\) for the 3D box from the lattice-point/sphere-octant picture, obtaining \(N = (\pi/6)\tilde E^{3/2}\) with \(\tilde E\) the energy in ground-scale units.

  2. Write count_states(E, d=3): enumerate the positive-integer lattice with numpy.meshgrid up to \(n_{\max} = \lfloor\sqrt E\rfloor\) per axis and sum the boolean condition \(|\mathbf{n}|^2 \le E\). Write this one yourself — the implementation is the lesson.

  3. Tabulate the staircase \(N(E)\) for \(E = 10^2, 10^3, 10^4\).

  4. Confirm the ratio to the continuum formula converges toward \(1\) (roughly \(0.78 \to 0.93 \to 0.98\)) and show the deficit scales as \(E^{-1/2}\) — identify it as the surface (Weyl) correction, the boundary’s share of the count.

  5. Differentiate to the density of states \(g(\varepsilon) \propto \sqrt\varepsilon\) and state the master move of the volume: \(\sum_{\text{states}} \to \int g(\varepsilon)\,d\varepsilon\). (Prose part.)

the 3D staircase vs the continuum octant volume:
  E =     100:  N =      410   (π/6)E^(3/2) =       523.6   ratio = 0.7830   deficit·√E = 2.170
  E =    1000:  N =    15399   (π/6)E^(3/2) =     16557.6   ratio = 0.9300   deficit·√E = 2.213
  E =   10000:  N =   511776   (π/6)E^(3/2) =    523598.8   ratio = 0.9774   deficit·√E = 2.258
../../_images/281a7500bfa806d1c368c6cf23d0c495249e7a63095405ca10449f715bbe5b8a.png

Fig. 627 The density of states is a count. The exact staircase \(N(E)\) of a 3D hard-wall box — every riser one more quantum state, counted as a lattice point — against the continuum octant-volume formula \((\pi/6)E^{3/2}\) (amber). The staircase hugs the curve from below: the deficit, here still visible at low \(E\), shrinks as \(E^{-1/2}\) and is Weyl’s surface correction — the boundary’s share of the count — rather than a numerical error. By \(E = 10^4\) (Exercise 1’s table) the count reaches \(97.7\%\) of the volume, and the master move of the volume, \(\sum_{\text{states}} \to \int g(\varepsilon)\,d\varepsilon\) with \(g \propto \sqrt{\varepsilon}\), is quantitatively licensed.#

Validation 1#

✓  the counted N(E) converges to the continuum (π/6)E^(3/2) from below   [ratios 0.783 → 0.930 → 0.977]
✓  the deficit scales as E^(−1/2): the surface (Weyl) correction, not a numerical error   [deficit·√E = 2.17, 2.21, 2.26]
True

Exercise 2 — Dimensions and dispersions#

Geometry enters thermodynamics through the density of states — and two later notebooks hinge on the exponent. Cite Eq. 684.

  1. Repeat the count with the count_states you wrote in Exercise 1, now in 2D and 1D, and confirm \(N \propto E\) (2D: \(g\) constant) and \(N \propto \sqrt E\) (1D: \(g \propto \varepsilon^{-1/2}\)).

  2. Derive the general massive-particle result \(g(\varepsilon) \propto \varepsilon^{d/2-1}\) from the \(k\)-space ball volume.

  3. Derive the massless/linear-dispersion result \(g(\varepsilon) \propto \varepsilon^{d-1}\) and specialize to the 3D photon/phonon density of states \(\propto \varepsilon^2\).

  4. Flag the two consequences in prose: the constant 2D density of states is the fact behind the absence of ideal-gas BEC in two dimensions (§7.17), and the \(\varepsilon^2\) law is the engine of Planck’s law (§7.14).

2D:  N/( (π/4)E )  at E = 1e5:  0.9961   (g constant — N linear in E)
1D:  N/√E          at E = 1e6:  1.000000   (g ∝ ε^(−1/2))

massive:  g ∝ ε^(d/2 − 1)  →  ε^(−1/2), const, √ε  in d = 1, 2, 3
massless: g ∝ ε^(d − 1)    →  const, ε, ε^2        in d = 1, 2, 3
../../_images/f6292705496995c903a34346f95523a9f74e5ef18611575fbf84defdee299b7a.png

Fig. 628 Where geometry enters thermodynamics. The density-of-states shapes \(g(\varepsilon)\) for massive particles (\(\varepsilon \propto k^2\), solid): \(\varepsilon^{-1/2}\) in 1D (an integrable divergence at the band bottom, clipped by the axis), constant in 2D, and \(\sqrt{\varepsilon}\) in 3D — and for massless dispersion (\(\varepsilon \propto k\), dashed red) the 3D photon/phonon law \(\varepsilon^2\). Two of the volume’s dramas hang on these exponents: the constant 2D curve is why an ideal Bose gas never condenses in two dimensions (7.17), and the \(\varepsilon^2\) curve, fed to the Bose occupation, becomes Planck’s law and the \(T^4\) of Stefan–Boltzmann (7.14). Curves are shape-normalized at \(\varepsilon = 1\); prefactors carry units, exponents carry the physics.#

Validation 2#

✓  the 2D count confirms N ∝ E: the density of states is CONSTANT in two dimensions   [got 0.996093 vs expected 1 (rtol=0.01, atol=1e-09)]
✓  the 1D count confirms N ∝ √E (g ∝ ε^(−1/2), integrably divergent at the band bottom)   [got 1 vs expected 1 (rtol=0.01, atol=1e-09)]
True

Exercise 3 — The polylogarithm, three ways#

The master special function of quantum gases, built from scratch and cross-examined by three independent routes. Cite Eq. 685.

Two of the three routes are ours to write. The first is the defining series itself; the third runs Eq. 686 backwards — \(\int_0^\infty x^{s-1}/(z^{-1}e^x - 1)\,dx = \Gamma(s)\,\mathrm{Li}_s(z)\), derived in Exercise 4 and quoted here as a known identity — so the quadrature helper written below is the same one Exercises 4, 5 and 7 will run forwards. Its one numerical decision: at \(z = 1\) the integrand behaves like \(x^{s-2}\) at the origin, an integrable endpoint singularity that scipy.integrate.quad handles cleanly only when the endpoint is handed to it explicitly, so the range is split at \(x = 1\) (the §7.1 keyhole-benchmark treatment). The middle route, mpmath.polylog, is an independent arbitrary-precision implementation and is simply called.

  1. Write polylog(s, z, kmax) as a partial sum of \(\sum z^k/k^s\). Write this one yourself — the implementation is the lesson.

  2. State the \(k_{\max}\) used and evaluate the geometric tail bound \(|R| \le |z|^{k_{\max}+1}/((k_{\max}+1)^s\,(1-|z|))\) at the three test points below.

  3. Write bose_integral(s, z) for \(\int_0^\infty x^{s-1}/(z^{-1}e^x - 1)\,dx\) by scipy.integrate.quad, splitting the range at \(x = 1\). Write this one yourself — the implementation is the lesson.

  4. Validate at \((s, z) = (3/2, 0.5)\), \((2, 0.8)\), \((3, 0.3)\) against mpmath.polylog and against bose_integral\((s, z)/\Gamma(s)\) — three routes agreeing to at least eight digits.

  5. Prove the derivative identity \(z\,d\mathrm{Li}_s/dz = \mathrm{Li}_{s-1}(z)\) from the series, and verify it numerically at \((s, z) = (2, 0.6)\) by a central difference.

  6. Demonstrate the slow convergence at \(z = 1\) (partial sums of \(\zeta(3/2)\) after \(2\times10^6\) terms still \(\sim 10^{-3}\) short, matching the \(2/\sqrt{k_{\max}}\) tail estimate) and state the practical rule: near \(z = 1\), use the integral route or the mpmath reference. (Prose + one computation.)

  (s, z) = (1.5, 0.5):  tail bound after k_max = 200000: 0.0e+00
  (s, z) = (2.0, 0.8):  tail bound after k_max = 200000: 0.0e+00
  (s, z) = (3.0, 0.3):  tail bound after k_max = 200000: 0.0e+00

Li_s(z) three ways (series | mpmath | integral/Γ):
  Li_1.5(0.5):  0.624837020820 | 0.624837020820 | 0.624837020820
  Li_2.0(0.8):  1.074794600008 | 1.074794600008 | 1.074794600008
  Li_3.0(0.3):  0.312400177893 | 0.312400177893 | 0.312400177893

z dLi_2/dz at z = 0.6:  central diff = 0.9162907320
Li_1(0.6) = −ln(0.4)   = 0.9162907319   (exact 0.9162907319)

ζ(3/2) partial sum, 2e6 terms: 2.610961   (ζ(3/2) = 2.612375)
shortfall = 1.41e-03   vs the 2/√K tail estimate 1.41e-03
/tmp/ipykernel_5407/1719476134.py:66: RuntimeWarning: overflow encountered in exp
  f = lambda x: x ** (s - 1.0) / (np.exp(x) / z - 1.0)

Validation 3#

✓  the polylogarithm partial sums agree with mpmath.polylog to 8+ digits   [max|Δ| = 5.55112e-17 (rtol=1e-10, atol=1e-09)]
✓  the defining integral (quad/Γ) agrees with the other two routes   [max|Δ| = 7.21645e-15 (rtol=1e-08, atol=1e-09)]
✓  the derivative identity z dLi_s/dz = Li_(s−1)(z), with Li_1 = −ln(1−z) as the anchor   [max|Δ| = 8.18271e-11 (rtol=1e-07, atol=1e-09)]
✓  at z = 1 the series crawls: the 2e6-term shortfall matches the 2/√K zeta tail   [got 0.00141421 vs expected 0.00141421 (rtol=0.2, atol=1e-09)]
True

Exercise 4 — The Bose integrals, and two famous numbers#

The integrals the boson movements run on — with their closed forms falling to the residue zetas of §7.2. Cite Eq. 686.

  1. Derive \(\int_0^\infty x^{s-1}/(z^{-1}e^x - 1)\,dx = \Gamma(s)\,\mathrm{Li}_s(z)\) by the geometric expansion and term-by-term \(\Gamma\)-integration; note the keyhole anatomy (§7.1) of the integrand.

  2. Verify the general-\(z\) formula at \((s, z) = (3/2, 0.7)\): the bose_integral you wrote in Exercise 3 against \(\Gamma(s)\cdot\) your Exercise 3 polylog.

  3. Evaluate the \(z = 1\) closed forms and confirm \(\int x/(e^x-1) = \pi^2/6\) and \(\int x^3/(e^x-1) = \Gamma(4)\zeta(4) = \pi^4/15\) — recognizing \(\zeta(2)\) and \(\zeta(4)\) as the residue-summation values of §7.2, and naming \(\pi^4/15\) as Stefan–Boltzmann’s number, met three movements early.

  4. Evaluate \(\int \sqrt x/(e^x-1) = \Gamma(3/2)\,\zeta(3/2) \approx 2.3152\) and flag it as the integral BEC’s critical temperature (§7.17) is made of.

∫ x^(1/2)/(e^x/z − 1) dx at z = 0.7:
  quad route:        0.888994445015
  Γ(3/2)·Li_(3/2)(z): 0.888994445016

∫ x/(e^x−1) dx   = 1.644934066848   (π^2/6  = 1.644934066848)
∫ x^3/(e^x−1) dx = 6.493939402267   (π^4/15 = 6.493939402267)

∫ √x/(e^x−1) dx  = 2.315157373342   (Γ(3/2)ζ(3/2) = 2.315157373394)
/tmp/ipykernel_5407/1719476134.py:66: RuntimeWarning: overflow encountered in exp
  f = lambda x: x ** (s - 1.0) / (np.exp(x) / z - 1.0)

Validation 4#

✓  the Bose integral equals Γ(s)·Li_s(z): quadrature meets the series route at general z   [got 0.888994 vs expected 0.888994 (rtol=1e-09, atol=1e-09)]
✓  the z = 1 Bose integrals close via the residue zetas of §7.2 (π^2/6, π^4/15, Γ(3/2)ζ(3/2))   [max|Δ| = 5.25731e-11 (rtol=1e-07, atol=1e-09)]
True

Exercise 5 — The Fermi integrals and the eta function#

The fermionic family, one alternating sign away — and the arithmetic that relates the two statistics at \(z = 1\). Cite Eq. 686.

  1. Derive \(\int_0^\infty x^{s-1}/(z^{-1}e^x + 1)\,dx = -\Gamma(s)\,\mathrm{Li}_s(-z) = \Gamma(s)\,f_s(z)\) by the alternating geometric expansion.

  2. Write fermi_integral(s, z), the Exercise 3 quadrature with the denominator’s sign flipped. The \(+1\) keeps the integrand finite at the origin for every \(z > 0\), so — unlike the series route, which needs \(|z| < 1\) — this one is valid at any fugacity, which is how Exercise 8 reaches \(z = 5\) and \(50\).

  3. Verify at \((s, z) = (2, 0.7)\): fermi_integral against \(-\Gamma(2)\cdot\) your Exercise 3 polylog\((2, -0.7)\).

  4. Evaluate the \(z = 1\) case via the eta function \(\eta(s) = (1 - 2^{1-s})\,\zeta(s)\) and confirm \(\int x/(e^x+1) = \pi^2/12\).

  5. State the Bose–Fermi relation \(f_s(1) = (1 - 2^{1-s})\,g_s(1)\) and its meaning (the alternating sum removes the even terms twice) — the two statistics differ, at \(z = 1\), by an arithmetic factor. (Prose + check, the latter against your Exercise 3 bose_integral.)

∫ x/(e^x/z + 1) dx at z = 0.7:
  quad route:          0.605158402338
  −Γ(2)·Li_2(−0.7):    0.605158402338

∫ x/(e^x+1) dx = 0.822467033424   (η(2) = π^2/12 = 0.822467033424)

f_(3/2)(1)/g_(3/2)(1) by quadrature: 0.2928932188
(1 − 2^(1−s)) at s = 3/2:            0.2928932188
/tmp/ipykernel_5407/1138173491.py:33: RuntimeWarning: overflow encountered in exp
  f = lambda x: x ** (s - 1.0) / (np.exp(x) / z + 1.0)
/tmp/ipykernel_5407/1719476134.py:66: RuntimeWarning: overflow encountered in exp
  f = lambda x: x ** (s - 1.0) / (np.exp(x) / z - 1.0)

Validation 5#

✓  the Fermi integral equals −Γ(s)·Li_s(−z): quadrature meets the alternating series   [got 0.605158 vs expected 0.605158 (rtol=1e-09, atol=1e-09)]
✓  ∫ x/(e^x+1) dx = η(2) = π^2/12: the Fermi family closes via the eta function   [got 0.822467 vs expected 0.822467 (rtol=1e-07, atol=1e-09)]
✓  the Bose–Fermi relation f_s(1) = (1 − 2^(1−s)) g_s(1) at s = 3/2   [got 0.292893 vs expected 0.292893 (rtol=1e-07, atol=1e-09)]
True

Exercise 6 — The Sommerfeld expansion, with its remainder measured#

The low-temperature tool of every metal calculation — derived, then held to account: the truncation error must shrink sixteenfold when the temperature halves. Cite Eq. 687.

  1. Derive the expansion by integrating by parts and expanding \(H\) about \(\mu\), showing the second-order coefficient is set by the \(\pi^2/12\) Fermi integral (times two): the \(\zeta(2)\) of Basel, bound for the metals of §7.10.

  2. Write sommerfeld(H, Hprime, mu, T) for the leading two terms of Eq. 687, the \(T = 0\) step integral plus \((\pi^2/6)\,T^2 H'(\mu)\), with \(k_B = 1\).

  3. For \(H(\varepsilon) = \sqrt\varepsilon\) and \(\mu = 1\), compare against the exact scipy.integrate.quad evaluation of \(\int H\,n_F\,d\varepsilon\) at \(T = 0.10\) and \(T = 0.05\) (raised limit, split at \(\mu\)), recording errors near \(8\times10^{-5}\) and \(5\times10^{-6}\).

  4. Confirm the error ratio \(\approx 2^4\) (about \(17\)\(18\) measured), establishing the \(T^4\) remainder — the expansion’s honesty certificate.

Sommerfeld two-term estimate vs exact quad, H = √ε, μ = 1:
  T = 0.1:  exact = 0.6749714537   Sommerfeld = 0.6748913370   |error| = 8.01e-05
  T = 0.05:  exact = 0.6687273819   Sommerfeld = 0.6687228343   |error| = 4.55e-06

error ratio between T = 0.10 and T = 0.05: 17.6   (2^4 = 16)
/tmp/ipykernel_5407/512719189.py:58: RuntimeWarning: overflow encountered in exp
  n_F = lambda e, T=T_S: 1.0 / (np.exp((e - mu_S) / T) + 1.0)
../../_images/9c00dc0a82c915df427d9616c4ac63ef398e59c71d2b296cc80667c4a67e9301.png

Fig. 629 The anatomy behind the Sommerfeld expansion. The Fermi function \(n_F(\varepsilon)\) at \(T = 0.10\) (amber) and \(T = 0.05\) (dark) is a step smeared over a width \(\sim k_BT\) about \(\mu\); its derivative \(-\partial n_F/\partial\varepsilon\) (red, scaled) is the narrow symmetric peak that does all the work — integrating any smooth \(H(\varepsilon)\) against \(n_F\) is, after one integration by parts, integrating \(\int H\) against this peak. The peak’s zeroth moment gives the \(T=0\) step integral, its first moment vanishes by symmetry (no odd orders in the expansion), and its second moment, \((\pi^2/3)T^2\) — the \(\pi^2/12\) Fermi integral in disguise — supplies the \((\pi^2/6)T^2 H'(\mu)\) correction. The measured \(T^4\) remainder (error ratio \(\approx 17\) when \(T\) halves) is the expansion’s honesty certificate.#

Validation 6#

✓  the two-term Sommerfeld estimate lands within its promised accuracy at both temperatures   [errors 8.0e-05 (T=0.10), 4.5e-06 (T=0.05)]
✓  the Sommerfeld truncation error scales as T^4 (ratio ≈ 16 when T halves)   [measured ratio 17.6]
True

Exercise 7 — Stefan–Boltzmann’s constant, assembled early#

All the pieces of the blackbody’s most famous number already exist in this movement — put them together, three movements before the physics arrives. Cite Eq. 686, Eq. 684.

  1. Verify \(\zeta(4) = \pi^4/90\) numerically (partial sums with their fast \(1/(3N^3)\) tail) and recall its residue-summation origin (the \(\pi\cot\) method of §7.2).

  2. Assemble \(\int_0^\infty x^3/(e^x-1)\,dx = \Gamma(4)\,\zeta(4) = \pi^4/15\) and confirm with the bose_integral you wrote in Exercise 3.

  3. Combine with the 3D photon density of states \(\propto \varepsilon^2\) (Exercise 2) to argue, dimensionally, that the radiated energy density must scale as \(T^4\) — the Stefan–Boltzmann law’s shape, before its physics.

  4. State in prose what §7.14 will add (the mode counting’s prefactors, \(\mu = 0\) for photons, and the law’s derivation proper) — this exercise is the arithmetic scouting party.

ζ(4) by series + tail: 1.082323233711
π^4/90               = 1.082323233711
scipy.special.zeta(4) = 1.082323233711

∫ x^3/(e^x−1) dx = 6.493939402267   (π^4/15 = 6.493939402267)

u ∝ T^4 · ∫x^3/(e^x−1)dx: three powers of T from ε^2·ε, one from dε — the shape of Stefan–Boltzmann
/tmp/ipykernel_5407/1719476134.py:66: RuntimeWarning: overflow encountered in exp
  f = lambda x: x ** (s - 1.0) / (np.exp(x) / z - 1.0)

Validation 7#

✓  ζ(4) = π^4/90: the tail-corrected series confirms the residue value of §7.2   [got 1.08232 vs expected 1.08232 (rtol=1e-10, atol=1e-09)]
✓  ∫ x^3/(e^x−1) dx = π^4/15: Stefan–Boltzmann's number, assembled from the toolkit   [got 6.49394 vs expected 6.49394 (rtol=1e-08, atol=1e-09)]
True

Exercise 8 — A bounded function with a vertical tangent#

The notebook’s quiet climax: the mathematics that will force a phase transition, stated three movements before it detonates. Cite Eq. 688, Eq. 685.

  1. Evaluate \(g_{3/2}(z) = \mathrm{Li}_{3/2}(z)\) at \(z = 0.5, 0.9, 0.99, 0.999\) (the polylog you wrote in Exercise 3, whose geometric tail bound still holds at \(z = 0.999\); mpmath.polylog as the cross-reference) and confirm it increases monotonically toward, but never exceeds, \(\zeta(3/2) = 2.612\dots\).

  2. Using \(z\,d\mathrm{Li}_s/dz = \mathrm{Li}_{s-1}(z)\), compute the slope at \(z = 0.999\) (\(\approx 55\) and growing) and show it diverges as \(z \to 1\), since \(\mathrm{Li}_{1/2}(z) \to \infty\) there — a finite ceiling reached with a vertical tangent.

  3. Contrast with the Fermi function \(f_{3/2}(z)\), which grows without bound as \(z \to \infty\) (evaluate at \(z = 5, 50\) with the fermi_integral you wrote in Exercise 5, the series being useless past \(z = 1\)) — fermions never saturate.

  4. Pose the puzzle in prose, to be resolved in §7.17: if a physical quantity (it will be the excited-state occupancy of a Bose gas) is proportional to \(g_{3/2}(z)\), it cannot be pushed past \(2.612\) by any fugacity — so what does the system do when the physics demands more?

g_(3/2)(z) climbing toward its ceiling ζ(3/2) = 2.612…:
  z = 0.5:  series = 0.62483702   mpmath = 0.62483702   (ceiling − value = 1.9875)
  z = 0.9:  series = 1.61443853   mpmath = 1.61443853   (ceiling − value = 0.9979)
  z = 0.99:  series = 2.27166008   mpmath = 2.27166008   (ceiling − value = 0.3407)
  z = 0.999:  series = 2.50170847   mpmath = 2.50170847   (ceiling − value = 0.1107)
  slope dg_(3/2)/dz at z = 0.9:     4.47   (asymptote √(π/−ln z) = 5.46)
  slope dg_(3/2)/dz at z = 0.99:    16.39   (asymptote √(π/−ln z) = 17.68)
  slope dg_(3/2)/dz at z = 0.999:    54.63   (asymptote √(π/−ln z) = 56.04)

f_(3/2)(5)  = 2.284211
f_(3/2)(50) = 6.320456   (past the Bose ceiling 2.6124 and still climbing)
/tmp/ipykernel_5407/1138173491.py:33: RuntimeWarning: overflow encountered in exp
  f = lambda x: x ** (s - 1.0) / (np.exp(x) / z + 1.0)
/tmp/ipykernel_5407/1138173491.py:33: RuntimeWarning: overflow encountered in exp
  f = lambda x: x ** (s - 1.0) / (np.exp(x) / z + 1.0)
../../_images/f8c06cbfe2c0f869eea59c6142a87d045eb4316f6496d540bf531f31ce356285.png

Fig. 630 The mathematics that will force a phase transition. The Bose function \(g_{3/2}(z) = \mathrm{Li}_{3/2}(z)\) rises monotonically across the physical fugacity range \(0 < z \le 1\) but cannot exceed its ceiling \(\zeta(3/2) = 2.612\ldots\) (dashed), which it meets at \(z = 1\) with a vertical tangent — the slope \(\mathrm{Li}_{1/2}(z)/z\) (red points at \(z = 0.9, 0.99, 0.999\)) diverges like \(\sqrt{\pi/(-\ln z)}\) even as the function itself stays bounded. In §7.7 the excited-state occupancy of an ideal Bose gas will be proportional to this function, and no fugacity can push it past \(2.612\): when cooling demands more, the ground state must absorb a macroscopic remainder — Bose–Einstein condensation (§7.17). The Fermi companion \(f_{3/2}(z)\) (grey, continuing past \(z=1\) by the quadrature route) has no ceiling: fermions never saturate.#

Validation 8#

✓  g_(3/2)(z) by partial sums matches mpmath.polylog along the climb to z = 0.999   [max|Δ| = 4.44089e-16 (rtol=1e-09, atol=1e-09)]
✓  g_(3/2) is monotone increasing yet bounded by ζ(3/2) = 2.612…   [g(0.999) = 2.5017 < 2.6124]
✓  the slope Li_(1/2)(z)/z diverges as z → 1: the ceiling is met with a vertical tangent   [slopes 4.5 → 16.4 → 54.6]
✓  the Fermi function f_(3/2) grows without bound — past the Bose ceiling and climbing   [f(5) = 2.284, f(50) = 6.320]
True

Exercise 9 — The vocabulary, complete#

Movement 0 is finished, and it ends where toolkits should — with everything sharpened and nothing yet cut. We can count states and turn sums into integrals, in any dimension and for any dispersion, with the count’s own error term identified as boundary physics rather than noise. We own the polylogarithm and its Gamma–zeta family, each value checkable three independent ways, and we know exactly where the series route dies and what replaces it. The Bose and Fermi integrals close on demand, their famous values — \(\pi^2/6\), \(\pi^4/15\), \(\Gamma(3/2)\zeta(3/2)\), \(\pi^2/12\) — already in hand with their §7.2 pedigrees attached. And the Sommerfeld expansion stands ready for the metals, its error scaling measured rather than promised.

Two loaded facts travel forward. A density of states that is constant in two dimensions, and a Bose function that cannot exceed \(2.612\) — one will forbid a condensation and the other will force one, and both are already fully proven, waiting only for physics to walk into them. There is a particular pleasure in meeting a famous constant before its fame: \(\pi^4/15\) is, for now, just an integral we evaluated three ways; in nine notebooks it will be the reason hot things glow with a fourth-power fury. The toolkit knows things the physics hasn’t said yet.

The next notebook (§7.4) begins the physics proper: the density matrix of §6.26 acquires a temperature, \(\rho = e^{-\beta H}/Z\), and the quantum canonical ensemble opens the volume it has all been for.

Notebook summary#

Movement 0 closes with the volume’s working equipment, each piece built from scratch and checked against something independent.

  • The density of states is a count Eq. 683: box modes are lattice points in a \(k\)-space octant; the counted staircase converges to \((\pi/6)E^{3/2}\) (ratio \(0.98\) by \(E = 10^4\)) with an \(E^{-1/2}\) deficit that is Weyl’s surface correction — boundary physics, not error. The master move is licensed: \(\sum_{\text{states}} \to \int g(\varepsilon)\,d\varepsilon\).

  • Geometry decides Eq. 684: \(g \propto \varepsilon^{d/2-1}\) (massive) and \(\varepsilon^{d-1}\) (massless), each verified by counting — the constant 2D DOS (no 2D BEC, §7.17) and the \(\varepsilon^2\) photon DOS (Planck’s engine, §7.14) planted as the volume’s two bombs.

  • The polylogarithm, three ways Eq. 685: partial sums with stated \(k_{\max}\) and geometric tail bound, mpmath.polylog, and the defining integral agree to eight digits; the index-lowering derivative identity is proved from the series; and the \(z = 1\) crawl (\(1.4\times10^{-3}\) short after \(2\times10^6\) terms, matching the \(2/\sqrt K\) tail) fixes the working rule near the endpoint.

  • The Bose and Fermi integrals Eq. 686: \(\Gamma(s)\mathrm{Li}_s(\pm z)\) by one geometric expansion (the keyhole anatomy of §7.1, one Boltzmann factor richer); closed forms via the residue zetas of §7.2\(\pi^2/6\), \(\pi^4/15\) (Stefan–Boltzmann’s number, three movements early), \(\Gamma(3/2)\zeta(3/2)\) (BEC’s integral) — and the Fermi family through \(\eta(s) = (1-2^{1-s})\zeta(s)\), the two statistics differing at \(z = 1\) by pure arithmetic.

  • The Sommerfeld expansion Eq. 687: derived from the moments of the \(-\partial n_F/\partial\varepsilon\) peak (its second moment is the \(\pi^2/12\) integral), and held to account — the truncation error shrinks \(\approx 17\times\) when \(T\) halves: the \(T^4\) remainder, measured.

  • A bounded function with a vertical tangent Eq. 688: \(g_{3/2}(z)\) climbs monotonically to \(\zeta(3/2) = 2.612\dots\) and no further, with slope \(\mathrm{Li}_{1/2}(z)/z \to \infty\) (already \(55\) at \(z = 0.999\)), while the Fermi companion grows without bound. The refusal is the seed of Bose–Einstein condensation.

The vocabulary is complete; the physics begins next door.

Outlook#

  • The thermal density matrix (§7.4). The \(\rho\) of §6.26 acquires a temperature, \(\rho = e^{-\beta H}/Z\), and the quantum canonical ensemble opens the physics proper.

  • The occupations and the classical limit (§7.7, §7.8). The \(n_B\) and \(n_F\) these integrals will average — derived by ensembles, meeting the contour route of §7.2 — and the Maxwell–Boltzmann limit both statistics share.

  • The payoffs. Sommerfeld in the metals (§7.10); \(\varepsilon^2\) and \(\pi^4/15\) in the blackbody (§7.14); the bounded \(g_{3/2}\) detonating as Bose–Einstein condensation (§7.17).

  • Weyl’s law and spectral geometry — “can one hear the shape of a drum?” — the surface term of Exercise 1 grown into a field: a horizon, named.

  • Cross-reference §6.10/§6.16 (the modes counted here), §7.1 (the keyhole anatomy), §7.2 (the residue zetas; \(\zeta(3/2)\) flagged there, used here), §5.3 (the classical large-\(N\) toolkit this quantum one parallels).

Take this notebook with you
Use the download button (↓) in the toolbar above to save this notebook and run it yourself. The published notebooks ship without worked solutions; if you would like the reference solutions — to teach from or to check your own work — get in touch: hello@ramador.me.