Debye Model and the T-cubed Heat Capacity

Derivation

Debye Model and the T-cubed Heat Capacity

D-257 Home PU-303 Threads energy · waves · chance Depends on Quantization of Lattice Vibrations, bose-einstein-distribution

Statement

Treating a monatomic crystal of \(N\) atoms as an isotropic elastic continuum whose \(3N\) normal modes are phonons with a linear dispersion \(\omega = v_s k\) truncated at a Debye cutoff \(\omega_D\), the lattice internal energy is \(U = 9Nk_BT\,(T/\theta_D)^3\!\int_0^{\theta_D/T} x^3/(e^x-1)\,dx\), and its temperature derivative yields the molar heat capacity \(C_V = \tfrac{12\pi^4}{5}Nk_B\,(T/\theta_D)^3\) for \(T\ll\theta_D\) (the Debye \(T^3\) law) and \(C_V \to 3Nk_B\) for \(T\gg\theta_D\) (Dulong–Petit), where \(\theta_D = \hbar\omega_D/k_B\).

Why it matters

Measured heat capacities of insulating crystals fall as \(T^3\) at low temperature, not the constant \(3Nk_B\) that classical equipartition predicts. The Debye model captures this with a single material parameter, the Debye temperature \(\theta_D\), and simultaneously reproduces the classical high-temperature plateau. It is the canonical example of how quantizing collective modes and applying Bose–Einstein statistics rescues thermodynamics from the ultraviolet catastrophe of continuum elasticity.

The same phonon density of states and the same \(T^3\) counting reappear across condensed-matter and astrophysics: in the lattice contribution that must be subtracted to expose the electronic \(\gamma T\) term in metals, in the phonon energy density of a solid, and — with photons replacing phonons — in the Stefan–Boltzmann law, which shares the \(\int x^3/(e^x-1)\,dx\) integral.

Assumptions

Linear (acoustic) dispersion \(\omega = v_s k\) with a single effective sound speed \(v_s\). If dropped, \(g(\omega)\) is no longer \(\propto\omega^2\); optical branches and dispersion curvature enter, and the clean \(\omega^2\) density of states — and hence the exact \(T^3\) coefficient — is lost. Three phonon polarizations per wavevector, giving exactly \(3N\) modes. If dropped, the normalization \(\int_0^{\omega_D} g\,d\omega = 3N\) that fixes \(\omega_D\) fails, and the Dulong–Petit limit no longer comes out to \(3Nk_B\). A sharp Debye cutoff \(\omega_D\) replacing the true Brillouin-zone boundary. If dropped, the upper limit of every integral becomes anisotropic and branch-dependent; the model's virtue — that low-\(T\) results are independent of the cutoff — survives, but the crossover region shifts. Modes are independent harmonic oscillators (non-interacting phonons) obeying Bose–Einstein statistics with zero chemical potential. If dropped, anharmonic phonon–phonon interactions add temperature-dependent corrections and thermal expansion, so measured \(C_p\) drifts above the harmonic \(C_V\). The crystal is isotropic, so a scalar \(v_s\) suffices. If dropped, one longitudinal and two transverse speeds must be averaged as \(3/v_s^3 = 1/v_L^3 + 2/v_T^3\); using a single number misestimates \(\theta_D\).

Derivation

\[ dN_{\text{modes}} = 3\cdot\frac{V}{(2\pi)^3}\,4\pi k^2\,dk = \frac{3V}{2\pi^2}\,k^2\,dk \]

Count allowed wavevectors in a spherical shell of radius \(k\) under periodic boundary conditions: each mode occupies \((2\pi)^3/V\) of \(k\)-space, and the factor \(3\) is the three polarizations. A

\[ g(\omega)\,d\omega = \frac{3V}{2\pi^2}\,\frac{\omega^2}{v_s^2}\,\frac{d\omega}{v_s} = \frac{3V}{2\pi^2 v_s^3}\,\omega^2\,d\omega \]

Change variables with the linear dispersion \(k = \omega/v_s\), \(dk = d\omega/v_s\). This defines the Debye density of states \(g(\omega)\propto\omega^2\). A

\[ \int_0^{\omega_D} g(\omega)\,d\omega = \frac{3V}{2\pi^2 v_s^3}\cdot\frac{\omega_D^3}{3} = 3N \;\;\Longrightarrow\;\; \omega_D^3 = \frac{6\pi^2 N v_s^3}{V} \]

Fix the cutoff by demanding the total mode count equal \(3N\); this is the physical content that a discrete lattice has finitely many modes. A

\[ U = \int_0^{\omega_D} g(\omega)\,\hbar\omega\,\langle n(\omega)\rangle\,d\omega,\qquad \langle n(\omega)\rangle = \frac{1}{e^{\hbar\omega/k_BT}-1} \]

Each mode is a quantized oscillator whose mean occupation is the Bose–Einstein distribution (prior result); measuring energy from the zero-point reference, the thermal energy per mode is \(\hbar\omega\langle n\rangle\). A

\[ U = \frac{3V}{2\pi^2 v_s^3}\int_0^{\omega_D}\frac{\hbar\,\omega^3}{e^{\hbar\omega/k_BT}-1}\,d\omega \]

Insert \(g(\omega)\) from step 2. The integrand is purely a function of \(\omega\) and \(T\); all geometry sits in the prefactor. A

\[ x \equiv \frac{\hbar\omega}{k_BT},\qquad \omega = \frac{k_BT}{\hbar}x,\qquad \omega^3\,d\omega = \left(\frac{k_BT}{\hbar}\right)^{4} x^3\,dx \]

Nondimensionalize so temperature is scaled out of the integrand and into the limits and prefactor. The upper limit becomes \(x_D = \hbar\omega_D/k_BT = \theta_D/T\). B

\[ U = \frac{3V\,k_B^4 T^4}{2\pi^2 v_s^3 \hbar^3}\int_0^{\theta_D/T}\frac{x^3}{e^x-1}\,dx \]

Substitute step 6 into step 5, collecting the constants \(\hbar(k_BT/\hbar)^4 = k_B^4T^4/\hbar^3\). Legal because the map \(\omega\mapsto x\) is smooth and monotonic. B

\[ \frac{V}{v_s^3} = \frac{6\pi^2 N}{\omega_D^3},\qquad \theta_D^3 = \frac{\hbar^3\omega_D^3}{k_B^3} \;\Longrightarrow\; \frac{3V k_B^4 T^4}{2\pi^2 v_s^3\hbar^3} = 9Nk_BT\left(\frac{T}{\theta_D}\right)^3 \]

Eliminate \(V/v_s^3\) using step 3 and rewrite the prefactor through \(\theta_D\). The explicit \(v_s\) and \(V\) disappear in favour of the single measurable \(\theta_D\). B

\[ \boxed{\,U = 9Nk_BT\left(\frac{T}{\theta_D}\right)^{3}\int_0^{\theta_D/T}\frac{x^3}{e^x-1}\,dx\,} \]

Combine steps 7–8. This is the exact Debye internal energy for all \(T\); the remaining physics is in the two limits of the dimensionless integral. A

\[ C_V = \left(\frac{\partial U}{\partial T}\right)_V = 9Nk_B\left(\frac{T}{\theta_D}\right)^3\int_0^{\theta_D/T}\frac{x^4 e^{x}}{(e^x-1)^2}\,dx \]

Differentiate step 9. The \(T\)-dependence in both the prefactor and the upper limit is differentiated; the boundary term from the limit cancels against part of the prefactor derivative, leaving this compact form. C

\[ T\ll\theta_D:\quad \int_0^{\theta_D/T}\!\to\int_0^{\infty}\frac{x^3}{e^x-1}\,dx = \frac{\pi^4}{15} \]

At low temperature the upper limit \(\theta_D/T\to\infty\); the exponential cutoff makes the integrand's tail negligible, so the finite limit is replaced by the standard \(\Gamma(4)\zeta(4)\) integral, independent of \(\theta_D\). C

\[ U \to \frac{3\pi^4}{5}Nk_BT\left(\frac{T}{\theta_D}\right)^3,\qquad C_V = \frac{dU}{dT} = \frac{12\pi^4}{5}Nk_B\left(\frac{T}{\theta_D}\right)^{3} \]

Insert \(\pi^4/15\) into step 9 and differentiate; the low-\(T\) integral is now temperature-independent so \(C_V\propto T^3\) follows directly. B

\[ T\gg\theta_D:\quad e^x-1\approx x \Rightarrow \int_0^{\theta_D/T}\frac{x^3}{e^x-1}\,dx \approx \int_0^{\theta_D/T}\! x^2\,dx = \frac{1}{3}\left(\frac{\theta_D}{T}\right)^3 \]

At high temperature \(x\le\theta_D/T\ll1\), so linearize the Bose factor. Each mode reaches its classical energy \(k_BT\). B

\[ U \to 9Nk_BT\left(\frac{T}{\theta_D}\right)^3\cdot\frac{1}{3}\left(\frac{\theta_D}{T}\right)^3 = 3Nk_BT,\qquad C_V \to 3Nk_B \]

The \((T/\theta_D)^3\) prefactor cancels the \((\theta_D/T)^3\) from the integral, recovering equipartition over \(3N\) modes. A

Result

\[ C_V \;=\; \frac{12\pi^4}{5}\,Nk_B\left(\frac{T}{\theta_D}\right)^{3}\;\;(T\ll\theta_D),\qquad C_V \;\to\; 3Nk_B\;\;(T\gg\theta_D) \]

Reading. Deep in the cold regime only the long-wavelength modes with \(\hbar\omega\lesssim k_BT\) are thermally excited; their number grows as \((k_BT/\hbar v_s)^3 \propto T^3\), and each carries \(\sim k_B\) of heat capacity, so \(C_V\propto T^3\). As \(T\) rises past \(\theta_D\) every mode is fully excited and the capacity saturates at the classical Dulong–Petit value \(3Nk_B\) (per mole, \(3R\approx 24.9\ \mathrm{J\,mol^{-1}K^{-1}}\)). \(\theta_D\) is the single knob setting the crossover.

Units check. \(Nk_B\) has units \(\mathrm{J\,K^{-1}}\); \(12\pi^4/5\) and \((T/\theta_D)^3\) are dimensionless, so \(C_V\) is \(\mathrm{J\,K^{-1}}\) as required. In \(U=9Nk_BT(\cdots)\), \(Nk_BT\) is \(\mathrm{J}\) and the integral is dimensionless, giving energy. \(\theta_D=\hbar\omega_D/k_B\) has \(\mathrm{(J\,s)(s^{-1})/(J\,K^{-1})} = \mathrm{K}\).

Limiting cases

\(T\to 0\): \(C_V\to 0\) as \(T^3\), satisfying the third law (\(S\to0\)); no classical model achieves this.
\(T\ll\theta_D\): \(C_V = \frac{12\pi^4}{5}Nk_B(T/\theta_D)^3\); a plot of \(C_V/T\) versus \(T^2\) is a straight line through the origin (for an insulator).
\(T\gg\theta_D\): \(C_V\to 3Nk_B\), the Dulong–Petit law, with leading correction \(C_V = 3Nk_B[1-\tfrac{1}{20}(\theta_D/T)^2+\cdots]\).
Continuum limit \(\omega_D\to\infty\): the cutoff drops out entirely and only the \(T^3\) law survives at all \(T\) — the phonon analogue of blackbody radiation.
Single mode \(N=1\): reduces to one Einstein oscillator only if the \(\omega^2\) spectrum is collapsed to a delta function; Debye and Einstein agree at high \(T\) but differ at low \(T\) (Einstein gives exponential, not \(T^3\), suppression).

Breaks when

Metals at very low temperature. Below a few kelvin the conduction electrons contribute a linear \(\gamma T\) term that dominates the \(\propto T^3\) lattice part; the pure Debye law fails and one must fit \(C_V = \gamma T + \beta T^3\).
Anharmonic / high-temperature regime near melting. Phonon–phonon interactions, vacancy formation, and thermal expansion push the measured \(C_p\) above \(3Nk_B\); the harmonic independent-mode assumption breaks and \(C_p-C_V = \alpha^2 VT/\kappa_T\) is no longer negligible.
Layered or low-dimensional solids. In quasi-2D materials (graphite, layered crystals) the density of states is not \(\propto\omega^2\) over the relevant range; flexural modes give \(C_V\propto T\) or \(T^2\) behaviour instead of \(T^3\) at intermediate temperatures.
Glasses and amorphous solids. Two-level tunnelling states add a near-linear term and a "boson peak" excess over the Debye prediction below \(\sim 1\,\mathrm{K}\).

Failure modes

Dropping the polarization factor of 3, giving \(Nk_B\) instead of \(3Nk_B\) at high \(T\) and a \(T^3\) coefficient a factor of 3 too small.
Using \(\int_0^\infty x^4e^x/(e^x-1)^2\,dx\) inconsistently: the \(C_V\) integral equals \(4\pi^4/15\), not \(\pi^4/15\); mixing it with the \(U\) integral gives the wrong low-\(T\) coefficient.
Forgetting that the upper limit depends on \(T\) when differentiating \(U\), producing a spurious boundary term and the wrong \(C_V\).
Treating \(\theta_D\) as a fundamental constant rather than a fitted average of \(v_s\); quoting a single \(\theta_D\) as if it were exact across all temperatures.
Applying the \(T^3\) law to a metal at 1 K without subtracting the electronic \(\gamma T\) term, then concluding "Debye is wrong."
Confusing \(C_V\) and \(C_p\): experiments measure \(C_p\); at high \(T\) the two differ by the thermal-expansion term.

Discussion

The physical engine behind the \(T^3\) law is a phase-space count. At temperature \(T\) only phonons with \(\hbar\omega \lesssim k_BT\) are thermally populated, i.e. those inside a sphere of radius \(k_T \sim k_BT/\hbar v_s\) in wavevector space. The number of such modes scales as \(k_T^3\propto T^3\), and because each active mode contributes of order \(k_B\) to the heat capacity while the frozen-out modes contribute nothing, \(C_V\propto T^3\) follows without evaluating a single integral. The \(\pi^4/15\) is merely the precise numerical prefactor from doing the count carefully with Bose statistics.

Debye's construction is best seen against Einstein's earlier model, which placed all \(3N\) oscillators at one frequency. Einstein correctly predicted the drop below \(3Nk_B\) but got exponential freeze-out, \(C_V\sim e^{-\theta_E/T}\), because a single gapped mode has no low-energy excitations. Debye's insight was that acoustic modes are gapless — \(\omega\to0\) as \(k\to0\) — so there are always arbitrarily soft modes to excite, and it is precisely this gaplessness, encoded in the linear dispersion, that produces a power law rather than an exponential. The two models are the same physics with different densities of states; they agree at high \(T\) where only mode counting matters.

The connection to radiation is exact in form: replace phonons by photons, set \(v_s\to c\), remove the cutoff (\(\omega_D\to\infty\)) and drop the longitudinal polarization (photons have two), and \(U\propto T^4\) with the identical \(\int_0^\infty x^3/(e^x-1)\,dx=\pi^4/15\) becomes the Stefan–Boltzmann law. A phonon gas is a photon gas in a box with a finite number of modes and a finite sound speed; the \(T^3\) heat capacity and the \(T^4\) energy density are two readings of the same Bose integral.

The independence of the low-\(T\) coefficient from the cutoff \(\omega_D\) is a renormalization-style statement: infrared physics (long-wavelength acoustic phonons) governs the leading low-temperature thermodynamics, and the ultraviolet details of the true dispersion near the Brillouin-zone edge enter only through subleading corrections and through the overall value of \(\theta_D\) that sets the crossover scale. This is why the crude replacement of the real zone boundary by a sphere of equal mode count is so successful for \(C_V(T\to0)\) yet unreliable in the crossover region, where the actual density of states — with its van Hove singularities — matters.

Common misconceptions. The \(T^3\) law is not a property of "three dimensions plus quantum mechanics" alone — it requires gapless linear dispersion; optical phonons and any energy gap give Einstein-like exponential behaviour. \(\theta_D\) is not a phase-transition temperature and nothing special happens at \(T=\theta_D\); it is only the scale marking the smooth crossover from \(T^3\) rise to Dulong–Petit saturation.

Worked examples

\[ \text{Molar } C_V \text{ of copper at } T=10\ \mathrm{K},\quad \theta_D = 343\ \mathrm{K}\ \text{(lattice part only)} \]

Since \(T/\theta_D = 10/343 \approx 0.029 \ll 1\), the low-temperature \(T^3\) law applies. Set \(N=N_A\) so \(Nk_B=R\). A

\[ C_V = \frac{12\pi^4}{5}\,R\left(\frac{T}{\theta_D}\right)^3,\qquad \frac{12\pi^4}{5} = \frac{12(97.409)}{5} = 233.78 \]

Substitute symbols first; evaluate the numerical prefactor. A

\[ \left(\frac{10}{343}\right)^3 = (0.029155)^3 = 2.478\times10^{-5} \]

Insert the numbers into the dimensionless ratio. A

\[ C_V = 233.78\times 8.314\ \mathrm{J\,mol^{-1}K^{-1}}\times 2.478\times10^{-5} \]

Multiply prefactor, \(R\), and ratio. A

\[ C_V \approx 4.8\times10^{-2}\ \mathrm{J\,mol^{-1}K^{-1}} \]

Reading. At 10 K the lattice holds only about \(0.2\%\) of its Dulong–Petit capacity (\(3R=24.9\)). Measured copper \(C_V\) at 10 K is a few times larger because the electronic term \(\gamma T\) (with \(\gamma\approx 0.69\ \mathrm{mJ\,mol^{-1}K^{-2}}\), giving \(\sim 6.9\times10^{-3}\)) is comparable — illustrating exactly why metals need both terms.

Units check. \(R\) in \(\mathrm{J\,mol^{-1}K^{-1}}\) times dimensionless factors gives \(\mathrm{J\,mol^{-1}K^{-1}}\).

\[ \text{Estimate } \theta_D \text{ of copper from an average sound speed } v_s = 3810\ \mathrm{m\,s^{-1}} \]

Use \(\theta_D = (\hbar/k_B)\,v_s\,(6\pi^2 n)^{1/3}\) with number density \(n=N/V\), from \(\omega_D^3 = 6\pi^2 n v_s^3\). B

\[ n = \frac{\rho N_A}{M} = \frac{8960\ \mathrm{kg\,m^{-3}}\times 6.022\times10^{23}\ \mathrm{mol^{-1}}}{0.06355\ \mathrm{kg\,mol^{-1}}} = 8.49\times10^{28}\ \mathrm{m^{-3}} \]

Compute atomic number density from mass density and molar mass. A

\[ \omega_D = v_s\,(6\pi^2 n)^{1/3} = 3810\,\big(6\pi^2\times 8.49\times10^{28}\big)^{1/3} \]

Symbols before numbers; \(6\pi^2 = 59.22\), so \(6\pi^2 n = 5.03\times10^{30}\ \mathrm{m^{-3}}\) and its cube root is \(1.714\times10^{10}\ \mathrm{m^{-1}}\). A

\[ \omega_D = 3810\times 1.714\times10^{10} = 6.53\times10^{13}\ \mathrm{rad\,s^{-1}} \]

Multiply to get the Debye angular frequency. A

\[ \theta_D = \frac{\hbar\omega_D}{k_B} = \frac{1.055\times10^{-34}\times 6.53\times10^{13}}{1.381\times10^{-23}} \]

Convert cutoff frequency to a temperature. A

\[ \theta_D \approx 4.99\times10^{2}\ \mathrm{K} \approx 500\ \mathrm{K} \]

Reading. The single-\(v_s\) estimate (\(\sim500\ \mathrm{K}\)) overshoots the calorimetric value (\(343\ \mathrm{K}\)) because copper's transverse and longitudinal speeds differ substantially; using the proper average \(3/v_s^3 = 1/v_L^3 + 2/v_T^3\) (transverse modes, being slower, dominate) lowers \(\theta_D\) toward the measured value. The exercise shows both the model's predictive reach and the limits of a scalar sound speed.

Units check. \((\mathrm{m\,s^{-1}})(\mathrm{m^{-1}}) = \mathrm{s^{-1}}\) for \(\omega_D\); \((\mathrm{J\,s})(\mathrm{s^{-1}})/(\mathrm{J\,K^{-1}}) = \mathrm{K}\) for \(\theta_D\).

Problems

Show from the low-temperature Debye energy \(U = \tfrac{3\pi^4}{5}Nk_BT(T/\theta_D)^3\) that the phonon entropy is \(S = \tfrac{4\pi^4}{5}Nk_B(T/\theta_D)^3 = \tfrac{1}{3}C_V\), and confirm \(S\to0\) as \(T\to0\).
Solution
Since \(C_V = T\,(\partial S/\partial T)_V\) and \(C_V = \tfrac{12\pi^4}{5}Nk_B T^3/\theta_D^3\), we have \(dS = C_V\,dT/T = \tfrac{12\pi^4}{5}Nk_B\,T^2/\theta_D^3\,dT\). Integrating from 0: \(S = \tfrac{12\pi^4}{5}Nk_B\,\tfrac{T^3}{3\theta_D^3} = \tfrac{4\pi^4}{5}Nk_B(T/\theta_D)^3\). Comparing, \(S = C_V/3\). As \(T\to0\), \(S\propto T^3\to0\), satisfying the third law.
For silicon, \(\theta_D = 645\ \mathrm{K}\). Compute the molar \(C_V\) at \(T=20\ \mathrm{K}\) and state what fraction of Dulong–Petit this is.
Solution
\(T/\theta_D = 20/645 = 0.03101\); cube \(=2.981\times10^{-5}\). \(C_V = 233.78\times 8.314\times 2.981\times10^{-5} = 0.0580\ \mathrm{J\,mol^{-1}K^{-1}}\). Fraction of \(3R=24.94\): \(0.0580/24.94 = 2.3\times10^{-3}\), i.e. \(0.23\%\). (Silicon is a good insulator so no electronic term is needed.)
A crystal is measured to have \(C_V/T = a + bT^2\) with \(a = 0.70\ \mathrm{mJ\,mol^{-1}K^{-2}}\) and \(b = 4.8\times10^{-2}\ \mathrm{mJ\,mol^{-1}K^{-4}}\). Identify the two contributions and extract \(\theta_D\).
Solution
The linear-in-\(T\) piece \(aT\) is the electronic heat capacity \(\gamma T\); the \(bT^3\) piece is the Debye lattice term with \(b = \tfrac{12\pi^4}{5}R/\theta_D^3\). Solve \(\theta_D^3 = \tfrac{12\pi^4}{5}R/b = 233.78\times 8.314/(4.8\times10^{-5}\ \mathrm{J\,mol^{-1}K^{-4}}) = 1943.6/4.8\times10^{-5} = 4.049\times10^{7}\ \mathrm{K^3}\). Thus \(\theta_D = (4.049\times10^7)^{1/3} = 344\ \mathrm{K}\) (note \(b\) converted to \(\mathrm{J}\): \(4.8\times10^{-5}\ \mathrm{J\,mol^{-1}K^{-4}}\)). This is copper.
Evaluate the leading high-temperature correction: show \(C_V = 3Nk_B\big[1 - \tfrac{1}{20}(\theta_D/T)^2 + \cdots\big]\), and find the temperature at which \(C_V\) is \(1\%\) below \(3Nk_B\) for \(\theta_D = 300\ \mathrm{K}\).
Solution
Expand the exact \(C_V\) integrand for small \(x\): \(\tfrac{x^4e^x}{(e^x-1)^2} = x^2\big(1 - \tfrac{x^2}{12} + \cdots\big)\). Then \(C_V = 9Nk_B(T/\theta_D)^3\int_0^{\theta_D/T}(x^2 - x^4/12)\,dx = 9Nk_B(T/\theta_D)^3\big[\tfrac{1}{3}(\theta_D/T)^3 - \tfrac{1}{60}(\theta_D/T)^5\big] = 3Nk_B\big[1 - \tfrac{1}{20}(\theta_D/T)^2\big]\). Set \(\tfrac{1}{20}(\theta_D/T)^2 = 0.01\Rightarrow (\theta_D/T)^2 = 0.2\Rightarrow T = \theta_D/\sqrt{0.2} = 300/0.4472 = 671\ \mathrm{K}\).
Photon analogue: by replacing \(v_s\to c\), removing the cutoff, and keeping only 2 polarizations, adapt step 7 to derive the blackbody energy density \(u = U/V = \tfrac{\pi^2 k_B^4}{15\hbar^3 c^3}T^4\).
Solution
With 2 polarizations the density of states is \(g(\omega) = \tfrac{V\omega^2}{\pi^2 c^3}\) (two-thirds of the phonon prefactor \(\tfrac{3V}{2\pi^2}\) becomes \(\tfrac{V}{\pi^2}\), with \(v_s\to c\)). Then \(U = \tfrac{V\hbar}{\pi^2 c^3}\int_0^\infty \tfrac{\omega^3}{e^{\hbar\omega/k_BT}-1}\,d\omega\). Substituting \(x=\hbar\omega/k_BT\): \(U = \tfrac{V k_B^4 T^4}{\pi^2 c^3\hbar^3}\int_0^\infty\tfrac{x^3}{e^x-1}\,dx = \tfrac{V k_B^4 T^4}{\pi^2 c^3\hbar^3}\cdot\tfrac{\pi^4}{15}\). Dividing by \(V\): \(u = \tfrac{\pi^2 k_B^4}{15\hbar^3 c^3}T^4\), the Stefan–Boltzmann energy density, sharing the same \(\pi^4/15\) Bose integral as the phonon result.