The purpose of this post is to present a derivation of the Planck distribution which governs the phenomenon of so-called blackbody radiation. I feel a better (though longer) name would have been universal quantum thermal radiation (UQTR), because ultimately any object, independent of its material composition, glows when it bathes in the canonical ensemble with an ambient quantum gas of photons \(\gamma\) at some temperature \(T\).
Degeneracy of Accessible Microstates
To start, forget about whatever objects happen to be immersed in the photon gas, and just focus on the photon gas itself. It is in some ways similar to a classical ideal gas, because photons do not interact with each other \(\gamma+\gamma\to\gamma+\gamma\). If we imagine more precisely that the photon gas is enclosed within a cubic cavity of side lengths \(L\) and volume \(V=L^3\) (essentially a three-dimensional infinite potential well), then each photon \(\gamma\) lives in some momentum eigenstate \(|\textbf k\rangle\) with normalized position space plane wavefunction \(\langle\textbf x|\textbf k\rangle=\frac{1}{\sqrt{V}}e^{i\textbf k\cdot\textbf x}\). We then enforce \(L\)-periodic boundary conditions \(\langle\textbf x+L\hat{\textbf e}_i|\textbf k\rangle=\langle\textbf x|\textbf k\rangle\) in all \(3\) directions \(i=1,2,3\), thereby gluing together opposite faces of the cubic cavity as if topologically connected by a portal (the motivation for this is either that we don’t want to deal with edge effects or that we can think of this as a control subvolume within a larger sea of photons). This quantizes the allowed wavevectors of the photons in the gas as \(\textbf k_{\textbf n}=\frac{2\pi}{L}\textbf n\) where here \(\textbf n\in\textbf Z^3\) is some lattice point. The allowed energies of individual photons are therefore also quantized according to the relativistic dispersion relation for \(m=0\) massless particles \(E_{\textbf n}=\hbar c|\textbf k_{\textbf n}|=\frac{2\pi\hbar c}{L}|\textbf n|\) which by the way exhibits considerable degeneracy (e.g. \(E_{-\textbf n}=E_{\textbf n}, E_{(0,1,4)}=E_{(2,2,3)}\), etc.). For fun, here’s a ChatGPT rendering of the photon gas in the cubic cavity:
Now we want to learn something about the macroscopic observables of the photon gas. This means chugging the system through the machine of statistical mechanics, which in turn means computing the canonical partition function \(Z=\sum_{|\mu_{\gamma\text{ gas}}\rangle}e^{-\beta E_{|\mu_{\gamma\text{ gas}}\rangle}}\). Now, what are the accessible microstates \(|\mu_{\gamma\text{ gas}}\rangle\) of the photon gas? One might initially think that for a gas of \(N\) photons we simply need to assign some lattice point \(\textbf n_1,\textbf n_2,…,\textbf n_N\) to each of the \(N\) photons, and that would be a complete specification of the entire photon gas microstate \(|\mu_{\gamma\text{ gas}}\rangle\). However, there are \(2\) things wrong with this. The first is that in addition to each photon’s translational degrees of freedom (as described by \(\textbf n\)), each photon also has \(2\) orthogonal polarization states, analogous to the spin angular momentum qubit degree of freedom in spin-\(1/2\) quantum particles such as electrons \(e^-\). The second point is deeper and more subtle; the statement “for a gas of \(N\) photons” implicitly assumes that the number of photons \(N\) is conserved \(\dot N=0\) as one would expect in the canonical/\(NVT\) ensemble. But a fundamental property of photons is that they are not conserved \(\dot N\neq 0\) (through various absorption/emission interactions with matter like the walls of the cubic cavity, although this is a minor detail). Put another way, the chemical potential associated with creation or annihilation of photons vanishes \(\mu=0\) in the grand canonical ensemble (although note this is not about exchanging photons with a surrounding particle bath the way one might exchange atoms for instance; the total number of photons \(N\) is simply not conserved, whether in the photon gas itself, or in the composite system of the photon gas with the surrounding particle bath).
The upshot of this is that the space of accessible microstates of the photon gas is much larger than initially thought (though perhaps not strictly so in the cardinality sense!), consisting of all possible numbers of photons \(N\in\textbf N\) together with, for each possible photon population \(N\), some specification of \(N\) lattice points \(\textbf n_1,\textbf n_2,…,\textbf n_N\in\textbf Z^3\) and \(N\) polarizations \(|\psi_1\rangle,|\psi_2\rangle,…,|\psi_N\rangle\in\{|0\rangle,|1\rangle\}\) for each of the \(N\) photons (note that photon and boson rhyme so photons are bosons and there is no restriction on the accessible microstates coming from the Pauli exclusion principle which applies only to fermions). Thus, the canonical partition function for the photon gas is actually (remembering also the \(1/N!\) Gibbs factor):
\[Z=\sum_{N=0}^{\infty}\frac{1}{N!}\sum_{(\textbf n_1,\textbf n_2,…,\textbf n_N)\in\textbf Z^{3N}}\sum_{(|\psi_1\rangle,|\psi_2\rangle,…,|\psi_N\rangle)\in\{|0\rangle,|1\rangle\}^N}e^{-\beta E_{(\textbf n_1,\textbf n_2,…,\textbf n_N)}}\]
where the total energy of the photon gas in that particular microstate is degenerate with respect to their polarizations, instead depending only on the momenta \(E_{(\textbf n_1,\textbf n_2,…,\textbf n_N)}=\frac{2\pi\hbar c}{L}\sum_{i=1}^N|\textbf n_i|\). As a result, we can move the exponential outside the sum over polarization states \(\sum_{(|\psi_1\rangle,|\psi_2\rangle,…,|\psi_N\rangle)\in\{|0\rangle,|1\rangle\}^N}1=2^N\). We can also convert the exponential of the series of energies into a product \(e^{-\beta E_{(\textbf n_1,\textbf n_2,…,\textbf n_N)}}=\prod_{i=1}^Ne^{-2\pi\beta\hbar c|\textbf n_i|/L}\).
The series of products \(\sum_{(\textbf n_1,\textbf n_2,…,\textbf n_N)\in\textbf Z^{3N}}\prod_{i=1}^Ne^{-2\pi\beta\hbar c|\textbf n_i|/L}\) is tricky to calculate but can be accurately estimated with the aid of a useful approximation.
Digression on the Euler-Maclaurin Summation Formula
The Euler-Maclaurin summation formula asserts that, for any smooth function \(f\in C^{\infty}(\textbf R\to\textbf C)\), if we sum the values of \(f\) evaluated on a bunch of consecutive integers \(\sum_{n=a}^bf(n)\) with \(a,b\in\textbf Z\cup\{-\infty,\infty\}\), then this is often well-approximated by the integral \(\int_a^bf(n)dn\) of the function \(f\) on the interval \([a,b]\), treating \(n\) as though it were a continuous variable \(n\in\textbf R\) rather than an integer \(n\in\textbf Z\) (of course this is a two-way street, meaning that many integrals can also be approximated by suitable series). The canonical example of the Euler-Maclaurin summation formula concerns the harmonic series \(f(n)=1/n\), with the following diagram given to rationalize why the Euler-Maclaurin approximation should be a good approximation:
However, the full Euler-Maclaurin formula also acknowledges that the approximation is imperfect by explicitly giving a formula for the error/discrepancy between the series and its integral approximation:
\[\sum_{n=a}^bf(n)=\int_a^bf(n)dn+\frac{f(a)+f(b)}{2}+\sum_{n=1}^{\infty}\frac{B_{2n}}{(2n)!}\left(f^{[2n-1]}(b)-f^{[2n-1]}(a)\right)\]
where the asymptotic series on the right features the non-zero even Bernoulli numbers \(B_2=1/6,B_4=-1/30\), etc. This interesting sequence of rational numbers were historically first discovered in the context of Faulhaber’s formula \(1^p+2^p+…+n^p=\frac{1}{p+1}\left(B_0{{p+1}\choose{0}}n^{p+1}+B_1{{p+1}\choose{1}}n^{p}+…+B_p{{p+1}\choose{p}}n\right)\) for sums of powers, but arguably the real utility of the Bernoulli numbers lies in their appearance in the Euler-Maclaurin summation formula (indeed Faulhaber’s formula is just the special case of the Euler-Maclaurin summation formula with \(f(n)=n^p\)). The quickest way to compute these Bernoulli numbers is to use a Pascal-like triangle construction (details can be found in this paper of Kawasaki and Ohno). The idea is to first write down a couple terms of the harmonic sequence \(1, 1/2,1/3,…\) in the first row. Then, whereas in Pascal’s triangle one would compute the number “below” two numbers by simply adding them, here instead one first has to subtract the two (left number minus right number, or if you forget then bigger number minus smaller number to always obtain a positive difference), and then multiply the result by the denominator of whichever term in the original harmonic sequence \(1, 1/2,1/3,…\) lies on the same downward-sloping “diagonal column” so to speak. A picture is worth \(10^3\) words:
The Bernoulli numbers \(B_0=1,B_1=1/2,B_2=1/6\), etc. are then simply read off as the elements on the first (left-most) downward-sloping diagonal column (SHOULD INCLUDE SOME COMMENTS ABOUT WHEN EULER-MACLAURIN SUMMATION IS VALID AND RELATE IT TO THE STATISTICAL MECHANICAL APPLICATION FUNCTION WHICH IS EXPONENTIALLY DECAYING!).
Back to the Physics
All this is to say that we expect we can approximate the series of products earlier \(\sum_{(\textbf n_1,\textbf n_2,…,\textbf n_N)\in\textbf Z^{3N}}\prod_{i=1}^Ne^{-2\pi\beta\hbar c|\textbf n_i|/L}\) over the discrete lattice \(\textbf Z^{3N}\) by an integral over the smooth manifold \(\textbf R^{3N}\) instead, obtaining a \(3N\)-dimensional analog of the Euler-Maclaurin summation formula:
\[\sum_{(\textbf n_1,\textbf n_2,…,\textbf n_N)\in\textbf Z^{3N}}\prod_{i=1}^Ne^{-2\pi\beta\hbar c|\textbf n_i|/L}\approx \int_{(\textbf n_1,\textbf n_2,…,\textbf n_N)\in\textbf R^{3N}}\prod_{i=1}^Ne^{-2\pi\beta\hbar c|\textbf n_i|/L}d^3\textbf n_i\]
Notice that the \(N\) integrals for each photon completely decouple from each other (as expected because the photons are non-interacting):
\[\int_{(\textbf n_1,\textbf n_2,…,\textbf n_N)\in\textbf R^{3N}}\prod_{i=1}^Ne^{-2\pi\beta\hbar c|\textbf n_i|/L}d^3\textbf n_i=\prod_{i=1}^N\iiint_{\textbf n_i\in\textbf R^3}e^{-2\pi\beta\hbar c|\textbf n_i|/L}d^3\textbf n_i\]
Moreover, each of these integrals are identical (because each photon is an identical boson) so the product becomes an \(N\)-fold exponentiation:
\[\prod_{i=1}^N\iiint_{\textbf n_i\in\textbf R^3}e^{-2\pi\beta\hbar c|\textbf n_i|/L}d^3\textbf n_i=\left(\iiint_{\textbf n\in\textbf R^3}e^{-2\pi\beta\hbar c|\textbf n|/L}d^3\textbf n\right)^N\]
Expressing the volume element \(d^3\textbf n=4\pi|\textbf n|^2d|\textbf n|\) in spherical coordinates, the integral comes out as:
\[\iiint_{\textbf n\in\textbf R^3}e^{-2\pi\beta\hbar c|\textbf n|/L}d^3\textbf n=4\pi\int_0^{\infty}|\textbf n|^2e^{-2\pi\beta\hbar c|\textbf n|/L}d|\textbf n|=4\pi\left(\frac{L}{2\pi\beta\hbar c}\right)^3\Gamma(3)=\frac{V}{\pi^2\beta^3\hbar^3c^3}\]
So assembling everything together, the canonical partition function \(Z\) for the photon gas becomes:
\[Z=\sum_{N=0}^{\infty}\frac{(2V/\pi^2\beta^3\hbar^3c^3)^N}{N!}=e^{2V/\pi^2\beta^3\hbar^3c^3}\]
Unfortunately, I believe the correct answer is actually \(Z=e^{\pi^2 V/45\beta^3\hbar^3c^3}\). My answer is correct up to a numerical factor of \(\zeta(4)=\pi^4/90\) which I cannot for the life of me figure out where I forgot to incorporate.