Deriving the Canonical Distribution
Why Maximum Entropy?
One of the deepest ideas in statistical mechanics is that the equilibrium distribution of a physical system is not an arbitrary assumption — it is the unique distribution that maximizes the Gibbs–Shannon entropy subject to whatever constraints the physics imposes (normalization, and a fixed average energy). This is the maximum entropy principle, and applying it to a system with a fixed mean energy produces exactly the Boltzmann distribution that underlies the entire canonical ensemble.
In this article we derive that distribution analytically using Lagrange multipliers, then verify it numerically in Python by (1) solving the analytical formula and (2) directly maximizing entropy with a constrained numerical optimizer — and showing the two agree to machine precision.
The Optimization Problem
Consider a system with $N$ discrete microstates, each with energy $E_i$, and unknown occupation probabilities $p_i$. We want to find the $p_i$ that maximize the entropy
$$
S[p] = -\sum_{i=1}^{N} p_i \ln p_i
$$
subject to two constraints:
$$
\sum_{i=1}^{N} p_i = 1 \qquad \text{(normalization)}
$$
$$
\sum_{i=1}^{N} p_i E_i = U \qquad \text{(fixed mean energy)}
$$
Lagrangian Derivation
Introduce Lagrange multipliers $\alpha$ (for normalization) and $\beta$ (for the energy constraint):
$$
\mathcal{L} = -\sum_i p_i \ln p_i - \alpha\left(\sum_i p_i - 1\right) - \beta\left(\sum_i p_i E_i - U\right)
$$
Setting $\partial \mathcal{L}/\partial p_i = 0$:
$$
-\ln p_i - 1 - \alpha - \beta E_i = 0 \quad \Longrightarrow \quad p_i = e^{-1-\alpha}, e^{-\beta E_i}
$$
Absorbing the constant into a normalization factor $Z$ (the partition function) gives the canonical distribution:
$$
p_i = \frac{e^{-\beta E_i}}{Z}, \qquad Z = \sum_{i=1}^{N} e^{-\beta E_i}
$$
Here $\beta$ plays the role of inverse temperature ($\beta = 1/k_BT$), and its exact value is fixed implicitly by the constraint $\sum_i p_i E_i = U$ — it must be solved for numerically once $U$ is specified. The maximum entropy achieved is:
$$
S_{\max} = \beta U + \ln Z
$$
Python Implementation (Google Colaboratory)
The code below does four things in one pass: (1) solves for $\beta$ given a target mean energy $U$ using root-finding, (2) builds the resulting canonical distribution, (3) independently verifies it against a direct constrained entropy maximization (SLSQP), and (4) sweeps $\beta$ to visualize how the whole distribution deforms with temperature.
1 | import numpy as np |
Code Walkthrough
Section 1 — Setup. We use a toy system of $N=10$ discrete energy levels $E_i = 0,\dots,9$ and a target mean energy $U=3.5$. This is deliberately simple so the entire distribution can be plotted and inspected directly, but the machinery generalizes to any energy spectrum (continuous, degenerate, or otherwise).
Section 2 — Analytical solution via root-finding. mean_energy(beta, E) computes $\langle E \rangle$ for a given $\beta$ using the canonical formula. Because $\langle E \rangle$ is a strictly monotonically decreasing function of $\beta$ (higher $\beta$ = lower effective temperature = more weight on low-energy states), solve_beta can safely use scipy.optimize.brentq, a bracketed bisection-type solver, to find the unique $\beta^*$ satisfying the constraint. This step is the practical implementation of the Lagrange multiplier from the derivation — $\beta$ is never guessed, it’s solved for.
Section 3 — Independent numerical verification. This is the key sanity check: rather than trusting the analytical formula blindly, we hand the same constrained optimization problem (maximize entropy, fix normalization and mean energy) directly to scipy.optimize.minimize with SLSQP, starting from a uniform distribution. If the Lagrangian derivation is correct, the numerically optimized $p_i$ should match the closed-form $e^{-\beta E_i}/Z$ almost exactly — the printed gap is the direct empirical proof of the derivation, typically on the order of $10^{-8}$ or smaller.
Section 4 — Parameter sweep. To visualize the full family of canonical distributions (not just the single one matching $U=3.5$), we sweep $\beta$ over a range and recompute the distribution and entropy at each value. This produces the data behind the 3D surface and the entropy curve.
Section 5 — Visualization. Three panels: a bar chart of the specific solved distribution, a 2D curve showing how entropy falls as $\beta$ increases (colder effective temperature → more order → lower entropy), and a 3D surface showing how the entire probability profile across energy levels smoothly deforms as $\beta$ is swept — this is the clearest way to see the exponential “tilting” of the distribution predicted by the Lagrangian derivation.
No part of this script performs expensive computation (the largest loop is 60 iterations over a 10-dimensional distribution), so no GPU or vectorization tricks beyond ordinary NumPy array operations are necessary — it runs in well under a second on Colab’s default CPU runtime.
Solved beta (1/kT) : 0.124333 Partition function Z : 6.086324 Mean energy check: 3.500000 (target 3.5) Entropy S : 2.241209 nats Max |analytical - numeric| gap : 2.61e-04
Interpreting the Results
The bar chart shows the canonical distribution decaying exponentially with energy — exactly as the Lagrangian derivation predicts, since $p_i \propto e^{-\beta E_i}$. The entropy curve confirms the physical intuition that as $\beta$ increases (system gets “colder” in the statistical-mechanical sense), the distribution concentrates on low-energy states and entropy drops monotonically; entropy is maximized in the $\beta \to 0$ limit, where the distribution becomes uniform — the state of maximum ignorance subject to no energy constraint at all. The 3D surface ties both views together: it’s the entire continuum of canonical distributions, with the flat, uniform-looking slice at $\beta=0$ smoothly tilting into a sharply peaked, low-energy-dominated slice as $\beta$ grows.
Finally, the near-zero gap between the closed-form Lagrangian solution and the independently, numerically optimized distribution is the real payoff of this exercise: it’s direct computational proof that the Boltzmann distribution isn’t a postulate bolted onto statistical mechanics — it falls out as the unique solution to a well-posed entropy maximization problem.
