Quantum Statistics of Ideal Gases: Occupation Numbers

An ideal gas of identical quantum particles is completely described by one list: how many particles sit in each single-particle energy level. Because identical particles cannot be labelled, a microstate is not "which particle is where" but only the set of occupation numbers $\{n_\varepsilon\}$ — and this single shift in bookkeeping is what makes quantum statistics tractable. This article opens the ideal-quantum-gas chapter by unifying the two occupation laws into one formula and setting up the master equation that every later computation — the Fermi sea, Bose–Einstein condensation, Planck's law, the Debye solid — specializes. We take the microscopic input from the quantum-mechanics track. Which particles obey which law is owned by aq-qm-bosons-fermions: "integer spin means boson, half-integer spin means fermion," with fermions exclusive because "put two in the same state and the antisymmetric combination cancels itself to zero." The hard cap on fermion occupancy is owned by aq-qm-pauli-exclusion-principle, where the antisymmetric state "vanishes when the labels coincide" so the occupancy is "forced to $n\in\{0,1\}$ ." This article does not re-derive antisymmetry or exclusion; it takes those rules as the allowed range of the integers $n_\varepsilon$ and builds the thermodynamic consequence. The payoff is a single unified result. Working in the grand canonical ensemble, the grand partition function factorizes into an independent product over single-particle levels, and each level's mean occupation comes out as $\langle n_\varepsilon\rangle=1/(e^{\beta(\varepsilon-\mu)}\mp1)$ — the lower sign ( $-$ ) for bosons, the upper ( $+$ ) for fermions. The two celebrated distributions differ by a single sign. The full derivations of each live in their own homes (aq-sm-bose-einstein-distribution, aq-sm-fermi-dirac-distribution); here we assemble both from one calculation, then write the master sum $N=\sum_\varepsilon\langle n_\varepsilon\rangle$ that fixes the chemical potential and launches the chapter.

Key takeaways

A microstate of an ideal quantum gas is a set of occupation numbers $\{n_\varepsilon\}$ ; quantum statistics enters only through which integers are allowed — $n\in\{0,1\}$ (fermions) or $n\in\{0,1,2,\dots\}$ (bosons).
In the grand canonical ensemble the grand partition function factorizes over single-particle levels, so each level is an independent open subsystem.
One single-level sum, evaluated with the two allowed ranges, gives both laws at once: $\langle n_\varepsilon\rangle=1/(e^{\beta(\varepsilon-\mu)}\mp1)$ , lower sign bosons, upper sign fermions.
The individual derivations are owned by aq-sm-bose-einstein-distribution and aq-sm-fermi-dirac-distribution; this article unifies them and sets up the chapter.
The master equation $N=\sum_\varepsilon\langle n_\varepsilon\rangle\to\int_0^\infty g(\varepsilon)\langle n(\varepsilon)\rangle\,\mathrm{d}\varepsilon$ fixes $\mu(T,n)$ and feeds every SM5 result.
Key idea: because the grand sum factorizes over levels, one single-level partition function produces both statistics through a single sign flip in the denominator — $-1$ for bosons (occupation can diverge as $\varepsilon\to\mu^+$ ), $+1$ for fermions (occupation capped at $1$ ).

The core insight

Treat each single-particle level as its own open subsystem exchanging particles with a reservoir at $(\mu,T)$ . The grand partition function then factorizes, $\Xi=\prod_\varepsilon\Xi_\varepsilon$ , and the mean occupation of a level is $\langle n_\varepsilon\rangle=z\,\partial_z\ln\Xi_\varepsilon$ . Evaluating the single-level sum over its allowed occupancies — two terms for fermions, a geometric series for bosons — yields both quantum distributions in one stroke: $\langle n_\varepsilon\rangle=1/(e^{\beta(\varepsilon-\mu)}\mp1)$ . The whole difference between Bose and Fermi statistics is the $\mp1$ .

Where this lives today

The occupation-number picture is the working language of every quantum many-body calculation. In condensed matter it underlies the theory of metals, semiconductors, and superconductors — one writes $N=\int g(\varepsilon)\langle n(\varepsilon)\rangle\,\mathrm{d}\varepsilon$ and reads off carrier densities, Fermi surfaces, and heat capacities. In quantum optics the boson occupation of field modes is what a photon-counting detector measures, and its super-Poissonian variance is the signature of thermal light. In cold-atom laboratories both signs are realized on demand: trapped $^{40}$ K and $^6$ Li form ideal Fermi gases whose momentum-space step is imaged directly, while $^{87}$ Rb and $^{23}$ Na condense into a single macroscopically occupied mode — the same $\langle n_\varepsilon\rangle$ , dialled through its sign. Second quantization, the operator form of this representation, is the starting point of quantum field theory itself.

The physical picture

Picture the single-particle spectrum as a long rack of levels, each a tollbooth wired to a giant reservoir that hands out particles at a price set by the chemical potential $\mu$ . A level at energy $\varepsilon$ accepts a particle with grand-canonical weight $e^{-\beta(\varepsilon-\mu)}$ each. The only difference between the two kinds of matter is the booth's capacity: a fermion booth holds at most one occupant (full or empty), while a boson booth has no ceiling, so its average filling can grow without bound as the price falls toward the level energy. That single structural difference — capacity one versus capacity infinite — is the entire content of quantum statistics, and it shows up as one sign. The first figure is the hero: the same reservoir, the same energy axis, three occupation curves. Fermi–Dirac never climbs above one (the Pauli ceiling); Maxwell–Boltzmann is the classical exponential; Bose–Einstein diverges as $\varepsilon\to\mu^+$ , the seed of condensation. Read left to right they are one formula with the $\mp1$ dialled through $+1,0,-1$ .

Figure

Three occupation laws from one formula: the sign flip splits Bose Maxwell and Fermi

The punchline: the $\mp1$ is the only difference. Dial it to $+1$ and you get a ceiling; to $-1$ and you get a divergence; drop it and you recover the classical exponential.

The second idea is why the grand canonical ensemble is the right tool. If we insisted on fixing the total particle number $N$ (canonical ensemble), the levels would be coupled — every added particle steals from the shared budget — and the sum would not factorize. Trading the constraint $N=\text{fixed}$ for a fixed chemical potential $\mu$ decouples the levels: each may independently hold whatever it likes, and we recover the true mean $N$ afterward by summing. The figure below shows that factorization pictorially.

Figure

Why the grand ensemble factorizes: independent tollbooths trading with one reservoir

Fixing $\mu$ instead of $N$ is what buys the factorization; the true $N$ is recovered as a mean at the end.

The core equations

The starting point is the grand partition function summed over every many-body microstate. Because a microstate is just the list $\{n_\varepsilon\}$ with energy $E=\sum_\varepsilon\varepsilon\,n_\varepsilon$ and number $N=\sum_\varepsilon n_\varepsilon$ , the grand-canonical weight $e^{-\beta(E-\mu N)}$ splits into a factor per level.

\Xi=\sum_{\{n_\varepsilon\}}e^{-\beta\sum_\varepsilon(\varepsilon-\mu)n_\varepsilon}=\prod_\varepsilon\Xi_\varepsilon,\qquad \Xi_\varepsilon=\sum_{n_\varepsilon}\big(z\,e^{-\beta\varepsilon}\big)^{n_\varepsilon} \qquad (1)

where: $\Xi$ — grand partition function; $\Xi_\varepsilon$ — single-level factor; $z=e^{\beta\mu}$ — fugacity; $\beta=1/k_BT$ . Physical reading: the exponential of a sum is a product, and relaxing the total- $N$ constraint in favour of $\mu$ makes the level occupancies independent — the one algebraic fact that this entire chapter rests on. The two statistics are two choices for the range of $n_\varepsilon$ . Fermions stop at $n=1$ ; bosons run to infinity as a geometric series that converges only when $z\,e^{-\beta\varepsilon}<1$ .

\Xi_\varepsilon^{\text{FD}}=1+z\,e^{-\beta\varepsilon},\qquad \Xi_\varepsilon^{\text{BE}}=\frac{1}{1-z\,e^{-\beta\varepsilon}}\quad(z\,e^{-\beta\varepsilon}<1) \qquad (2)

where: the upper labels FD/BE denote fermions/bosons; the boson series requires $\mu<\varepsilon$ for every level. Physical reading: the single stopped term at $n=1$ is the whole arithmetic difference between a two-term sum and an infinite one — everything downstream is a consequence of that. The mean occupation of a level is the logarithmic fugacity-derivative of its factor, the standard grand-canonical relation for a single independent mode.

\langle n_\varepsilon\rangle=z\left(\frac{\partial\ln\Xi_\varepsilon}{\partial z}\right)=\frac{1}{e^{\beta(\varepsilon-\mu)}\mp1} \qquad (3)

where: the lower sign ( $-$ ) is Bose–Einstein, the upper ( $+$ ) is Fermi–Dirac. Physical reading: this is the unified occupation law. The full step-by-step derivation of each case is owned by aq-sm-bose-einstein-distribution ( $-1$ ) and aq-sm-fermi-dirac-distribution ( $+1$ ); the value of this article is that one calculation yields both. Summing the grand potential over all levels gives the pressure, and summing the occupation gives the particle number — the equation that fixes $\mu$ .

\ln\Xi=\mp\sum_\varepsilon\ln\big(1\mp z\,e^{-\beta\varepsilon}\big),\qquad pV=k_BT\ln\Xi \qquad (4)

where: $pV=-\Omega=k_BT\ln\Xi$ is the grand potential relation. Physical reading: the sign in front and inside is the same $\mp$ ; bosons ( $-$ ) and fermions ( $+$ ) share one expression, differing only in that sign. Finally, the sums become integrals over the single-particle spectrum weighted by the density of states — the master equation of the whole chapter.

N=\sum_\varepsilon\langle n_\varepsilon\rangle=\int_0^\infty\!g(\varepsilon)\,\frac{\mathrm{d}\varepsilon}{e^{\beta(\varepsilon-\mu)}\mp1} \qquad (5)

where: $g(\varepsilon)$ — single-particle density of states, owned by aq-sm-density-of-states-free-particles (including the spin degeneracy $g_s=2s+1$ ). Physical reading: this implicit equation determines $\mu(T,n)$ ; feeding the same $g(\varepsilon)$ and the two signs into it produces the Fermi sea, Bose–Einstein condensation, and the photon and phonon gases in turn.

Derivations

Derivation ABoth quantum distributions from one factorized grand partition function

The goal is to obtain the unified occupation law $\langle n_\varepsilon\rangle=1/(e^{\beta(\varepsilon-\mu)}\mp1)$ for bosons and fermions together, from the single fact that the grand partition function factorizes over single-particle levels. We take as given only the allowed occupancy ranges owned by quantum mechanics — $n\in\{0,1\}$ for fermions (aq-qm-pauli-exclusion-principle: occupancy "forced to $n\in\{0,1\}$ ") and $n\in\{0,1,2,\dots\}$ for bosons (aq-qm-bosons-fermions: bosons in the "same state" carry no penalty). The route is: factorize, evaluate the single-level sum for each range, then differentiate. This card owns the unification; the individual algebra of each case is owned by the BE and FD articles.

Step 1. Write the grand sum in the occupation-number basis. A microstate is $\{n_\varepsilon\}$ with $E-\mu N=\sum_\varepsilon(\varepsilon-\mu)n_\varepsilon$ , so $\Xi=\sum_{\{n_\varepsilon\}}\prod_\varepsilon(z\,e^{-\beta\varepsilon})^{n_\varepsilon}=\prod_\varepsilon\big(\sum_{n_\varepsilon}(z\,e^{-\beta\varepsilon})^{n_\varepsilon}\big)$ . Why: the exponential of a sum factorizes, and once $\mu$ replaces the fixed- $N$ constraint the levels are independent, so the unconstrained multi-sum splits into a product of single-level sums.

Step 2. Evaluate the single-level sum for each allowed range. For fermions the sum has two terms, $\Xi_\varepsilon^{\text{FD}}=1+z\,e^{-\beta\varepsilon}$ ; for bosons it is a geometric series, $\Xi_\varepsilon^{\text{BE}}=1/(1-z\,e^{-\beta\varepsilon})$ , convergent when $z\,e^{-\beta\varepsilon}<1$ . Why: the exclusion cap stops the fermion sum after the $n=1$ term — exactly where the boson series would keep going — so the two ranges are the only input that distinguishes the statistics.

Step 3. Take the logarithm of both cases at once. Writing $x=z\,e^{-\beta\varepsilon}$ , the fermion factor gives $\ln\Xi_\varepsilon=+\ln(1+x)$ and the boson factor gives $\ln\Xi_\varepsilon=-\ln(1-x)$ , i.e. $\ln\Xi_\varepsilon=\pm\ln(1\pm x)$ with the upper sign for fermions. Why: the two logarithms have the same shape up to the sign of $x$ and an overall sign, which is what lets one formula carry both.

Step 4. Differentiate to get the mean. Using $\langle n_\varepsilon\rangle=x\,\partial_x\ln\Xi_\varepsilon$ , both cases give $\langle n_\varepsilon\rangle=x/(1\pm x)=1/(x^{-1}\pm1)=1/(e^{\beta(\varepsilon-\mu)}\pm1)$ , upper sign fermions. Relabelling the sign as $\mp$ (lower for bosons) gives Eq. (3). Why: the algebra collapses to a single rational function whose $\pm1$ is the entire fingerprint of the statistics — $+1$ enforces the ceiling $\langle n\rangle<1$ , $-1$ permits the divergence at $\varepsilon\to\mu^+$ . Verified: the two-term and geometric single-level sums, the derivative $z\,\partial_z\ln\Xi_\varepsilon$ reproducing $1/(e^{\beta(\varepsilon-\mu)}\mp1)$ for each sign, and the denominator difference of exactly $2$ are checked symbolically in verify_quantum-statistics-occupation-numbers.py.

Go deeper · Derivation B · 1 minThe master particle-number equation and the low-temperature Fermi closure

Derivation BThe master particle-number equation and the low-temperature Fermi closure

The goal is to turn the level-by-level occupation into the single equation that fixes the chemical potential and thereby launches every ideal-quantum-gas result. We convert the sum over discrete levels into an energy integral weighted by the density of states $g(\varepsilon)$ , write $N=\int g(\varepsilon)\langle n(\varepsilon)\rangle\,\mathrm{d}\varepsilon$ , and then check that in the fermion $T\to0$ limit this closes onto the standard Fermi energy — a concrete confirmation that the master equation is correct before it is used elsewhere. We do not re-derive $g(\varepsilon)$ here; that object is owned by aq-sm-density-of-states-free-particles.

Step 1. Convert sum to integral. For a macroscopic box the levels are dense, so $\sum_\varepsilon\to\int_0^\infty g(\varepsilon)\,\mathrm{d}\varepsilon$ with $g(\varepsilon)=(g_sV/4\pi^2)(2m/\hbar^2)^{3/2}\sqrt{\varepsilon}$ for free particles in 3D. Then $N=\int_0^\infty g(\varepsilon)\,\langle n(\varepsilon)\rangle\,\mathrm{d}\varepsilon$ . Why: the density of states packages the mode counting so the same integral serves both statistics — only the $\mp1$ in $\langle n\rangle$ changes.

Step 2. Specialize to fermions at $T=0$ . The Fermi function becomes a step $\theta(\mu-\varepsilon)$ , so $N=\int_0^{\varepsilon_F}g(\varepsilon)\,\mathrm{d}\varepsilon$ with $\mu(0)\equiv\varepsilon_F$ . Why: at zero temperature the exclusion ceiling forces every level below $\varepsilon_F$ to be exactly full and every level above exactly empty — the master integral truncates cleanly.

Step 3. Do the integral. With $g\propto\sqrt\varepsilon$ , $\int_0^{\varepsilon_F}A\sqrt\varepsilon\,\mathrm{d}\varepsilon=\tfrac23A\,\varepsilon_F^{3/2}$ , and setting this equal to $N$ and solving gives $\varepsilon_F=(\hbar^2/2m)(3\pi^2 n)^{2/3}$ with $n=N/V$ and $g_s=2$ . Why: the $\tfrac23$ from the $\sqrt\varepsilon$ integral and the $(3\pi^2n)^{2/3}$ from inverting are exactly the numbers that reappear in the Fermi-energy article — the master equation reproduces them, confirming its normalization.

Step 4. Read the physics. The same equation with the boson sign and $\mu\to\varepsilon_0^-$ produces a divergent ground-state term (Bose–Einstein condensation), and with $\mu=0$ produces the photon gas. Why: one master integral, three signs of physics — the unification is not cosmetic; it is the computational spine of the whole chapter. Verified: the $N=\tfrac23A\varepsilon_F^{3/2}$ integral, the closure to $\varepsilon_F=(\hbar^2/2m)(3\pi^2n)^{2/3}$ , and the numeric copper/sodium Fermi energies are checked in verify_quantum-statistics-occupation-numbers.py.

Limiting cases

The single formula $\langle n_\varepsilon\rangle=1/(e^{\beta(\varepsilon-\mu)}\mp1)$ has four diagnostic limits.

Case	Result	Meaning
Classical, $e^{\beta(\varepsilon-\mu)}\gg1$	$\langle n_\varepsilon\rangle\approx e^{-\beta(\varepsilon-\mu)}$	The $\mp1$ is negligible; both statistics reduce to Maxwell–Boltzmann when $n\lambda_T^3\ll1$
Degenerate fermions, $T\to0$	$\langle n_\varepsilon\rangle\to\theta(\varepsilon_F-\varepsilon)$	Sharp step: all levels below $\varepsilon_F$ filled, above empty — the Fermi sea
Bosons near threshold	$\langle n_{\varepsilon_0}\rangle\to\infty$ as $\mu\to\varepsilon_0^-$	Macroscopic ground-state occupation — the seed of Bose–Einstein condensation
Photons, $\mu=0$	$\langle n_\varepsilon\rangle=1/(e^{\beta\hbar\omega}-1)$	Number not conserved, so $\mu=0$ ; gives the Planck distribution

Worked examples

Example 1 — Count the microstates three ways and see indistinguishability at work.

Step 1. Setup. Put three particles into four single-particle states and count the allowed microstates for distinguishable (classical), boson, and fermion particles. This is the finite version of "which integers are allowed."

Step 2. Distinguishable. Each of the three labelled particles independently picks one of four states: $4^3=\mathbf{64}$ microstates. Labels matter, so all assignments are distinct.

Step 3. Bosons. A boson microstate is a set of occupation numbers summing to three over four states, i.e. the number of ways to place three identical balls in four boxes: $\binom{4+3-1}{3}=\binom{6}{3}=\mathbf{20}$ . Multiple occupancy is allowed, labels are gone.

Step 4. Fermions. Now no state may hold more than one particle, so we choose which three of the four states are occupied: $\binom{4}{3}=\mathbf{4}$ microstates. The drop from $64\to20\to4$ is indistinguishability plus exclusion, made concrete — and it is exactly why the occupation sums (two-term vs geometric) differ.

Example 2 — Fermi energy and temperature of copper from the master equation.

Step 1. Setup. Copper contributes one conduction electron per atom, $n=8.49\times10^{28}\,\mathrm{m^{-3}}$ . The $T=0$ master integral $N=\int_0^{\varepsilon_F}g\,\mathrm{d}\varepsilon$ closes to $\varepsilon_F=(\hbar^2/2m)(3\pi^2 n)^{2/3}$ for a free spin-½ gas.

Step 2. Fermi wavevector. $3\pi^2 n=3\pi^2(8.49\times10^{28})=2.514\times10^{30}\,\mathrm{m^{-3}}$ , so $k_F=(2.514\times10^{30})^{1/3}=\mathbf{1.36\times10^{10}\,m^{-1}}$ .

Step 3. Fermi energy. $\varepsilon_F=\dfrac{(1.0546\times10^{-34})^2(1.36\times10^{10})^2}{2(9.109\times10^{-31})}=1.13\times10^{-18}\,\mathrm{J}=\mathbf{7.0\ eV}$ , matching the measured value for copper.

Step 4. Fermi temperature. $T_F=\varepsilon_F/k_B=(1.13\times10^{-18})/(1.381\times10^{-23})=\mathbf{8.2\times10^{4}\,K}$ . Since room temperature is $\sim300\,\mathrm{K}\ll T_F$ , only the sliver of levels within $\sim k_BT$ of $\varepsilon_F$ ever changes occupation — the observation the rest of the chapter builds on. (Sodium, $n=2.65\times10^{28}\,\mathrm{m^{-3}}$ , gives $\varepsilon_F\approx\mathbf{3.2\ eV}$ by the same route.)

Origins · 1 minoccupation number and how this was discovered

occupation number

representation matured 1926–1928 · The occupation-number description of a many-body state grew out of the 1924–1926 work of Bose, Einstein, Fermi and Dirac on quantum gases, and was placed on a formal footing by Dirac and Jordan in 1927–1928 with the creation and annihilation operators of second quantization. Its power is exactly what this article uses: because a microstate is a list of occupancies rather than a labelled arrangement, the grand partition function factorizes over levels and the two quantum distributions fall out of one calculation.

History

1924

— Bose derives Planck's law by counting photon states, introducing the new statistics

1924

— Einstein extends Bose's counting to a gas of massive atoms and predicts condensation

1926

— Fermi formulates the statistics of an ideal gas obeying Pauli exclusion (Rendiconti Lincei)

1926

— Dirac independently derives the same statistics from antisymmetric wavefunctions (Proc. Roy. Soc.)

1927

— Dirac introduces creation and annihilation operators for the radiation field

1928

— Jordan and Wigner cast fermion occupation numbers into an operator algebra, formalizing the representation

Going further

Texts: Pathria & Beale, Statistical Mechanics (3rd ed.), Ch. 6 (The Theory of Quantum Statistics); Kardar, Statistical Physics of Particles, Ch. 7; Reif, Fundamentals of Statistical and Thermal Physics, Ch. 9; Sethna, Statistical Mechanics, Ch. 7; Landau & Lifshitz, Statistical Physics Part 1, §§53–56.
Microscopic owner: the allowed occupancy ranges are owned by the quantum-mechanics track — $n\in\{0,1\}$ for fermions by aq-qm-pauli-exclusion-principle (occupancy "forced to $n\in\{0,1\}$ "), and the boson/fermion census by aq-qm-bosons-fermions ("integer spin means boson, half-integer spin means fermion"). This article assembles the thermodynamic occupation from those rules.
Builds on: Density of States for Free Particles · Grand Canonical Ensemble · Indistinguishability and Symmetry of Wavefunctions.
Leads to: Degenerate Fermi Gas at Zero Temperature · Bose-Einstein Condensation in the Ideal Gas · Sommerfeld Expansion · Photon Gas and Planck Radiation Law · Chemical Potential of a Fermi Gas vs Temperature. The individual distributions are owned by aq-sm-bose-einstein-distribution and aq-sm-fermi-dirac-distribution.
Provenance: AllQFT `sm` export (Pathria & Beale, Ch. 6), aq-sm-quantum-statistics-occupation-numbers, ref STAT.IQGAS.QSOCC.1217. All factorization, sign, limit, variance, and Fermi-closure claims python-verified (24 checks, all PASS). The copper $\varepsilon_F\approx7.0$ eV, $T_F\approx8.2\times10^4$ K and sodium $\varepsilon_F\approx3.2$ eV reproduce the measured metal values.

Why does one formula give both distributions — is that just a notational trick?

It is more than notation. The grand partition function factorizes over single-particle levels for any ideal gas, so each level's occupation is governed by the same single-level sum $\Xi_\varepsilon=\sum_n(z\,e^{-\beta\varepsilon})^n$ . The only physical input that changes is the range of $n$ : two values for fermions, all non-negative integers for bosons. That single change turns the two-term sum into a geometric series, which is exactly what flips the $+1$ to a $-1$ in the denominator. So the unification reflects a real structural fact — one calculation, two allowed ranges — not a coincidence of symbols.

Where does this article stop and the Bose–Einstein / Fermi–Dirac articles begin?

This article owns the unification and the setup: the factorization of $\Xi$ , the single formula $\langle n_\varepsilon\rangle=1/(e^{\beta(\varepsilon-\mu)}\mp1)$ obtained from both allowed ranges at once, and the master equation $N=\int g\langle n\rangle\,\mathrm{d}\varepsilon$ that launches the chapter. The detailed algebra, structural identities, and physics of each case — the bunching and $\mu<\varepsilon_0$ constraint for bosons, the particle-hole symmetry and Fermi step for fermions — are owned by aq-sm-bose-einstein-distribution and aq-sm-fermi-dirac-distribution respectively. Think of this as the chapter's map; those two are the destinations.

Why must we use the grand canonical ensemble rather than the canonical one?

Because only in the grand ensemble do the single-particle levels become independent. If you fix the total number $N$ (canonical ensemble), then filling one level uses up part of a shared budget, so the occupancies are correlated and the grand sum cannot factorize into a product over levels. Trading the constraint $N=\text{const}$ for a fixed chemical potential $\mu$ removes that coupling: each level may hold whatever it likes, the product form holds, and the mean occupation follows from a single factor. You then recover the physical $N$ afterward as the ensemble average $\sum_\varepsilon\langle n_\varepsilon\rangle$ .

What fixes the chemical potential in practice?

The master equation $N=\int_0^\infty g(\varepsilon)/(e^{\beta(\varepsilon-\mu)}\mp1)\,\mathrm{d}\varepsilon$ is an implicit equation for $\mu$ at given density $n=N/V$ and temperature $T$ . You know $N$ , $T$ , and $g(\varepsilon)$ ; you solve for the $\mu$ that makes the integral come out to $N$ . For fermions at low temperature this gives $\mu\approx\varepsilon_F$ ; for bosons on cooling it drives $\mu$ up toward the ground-state energy until condensation intervenes; for photons, whose number is not conserved, there is no such constraint and $\mu=0$ identically.

Does the occupation being a "mean" mean the levels are not really quantized?

No — each level's actual occupancy $n_\varepsilon$ is always an integer (0 or 1 for fermions; any non-negative integer for bosons). The quantity $\langle n_\varepsilon\rangle$ is the ensemble average of that integer, which is generally fractional. For a fermion level, $\langle n_\varepsilon\rangle$ literally equals the probability the level is occupied, since $n$ takes only the values 0 and 1. For a boson level the average can exceed one because the level can genuinely hold many particles. The fractional value is a statement about probabilities over many copies of the system, not about a fractional particle.

How is this the same object as the density of states from quantum mechanics?

The density of states $g(\varepsilon)$ that appears in the master integral is exactly the one owned by aq-sm-density-of-states-free-particles, which in turn is the same $\rho(E)$ that the quantum-mechanics track uses to count final states — "both count quantum states per unit energy." This article does not re-derive it; it plugs it into the occupation sum. The occupation formula supplies the "how full each level is," and the density of states supplies the "how many levels there are per unit energy." Multiply and integrate, and you have the thermodynamics.