Solution (source code)

= Solution

The harmonic <phonon> contribution to the <internal energy> is
$$
U=\int_0^{\omega_D}g(\omega)\hbar\omega
\left[\frac12+\frac1{e^{\hbar\omega/(k_BT)}-1}\right]d\omega.
$$
For $T\gg\Theta_D$, every retained mode has $x=\hbar\omega/(k_BT)\ll1$. Expanding $1/(e^x-1)=1/x-1/2+O(x)$ cancels the explicit half quantum and gives $k_BT$ per mode to leading order. Therefore
$$
U=3Nk_BT+O\!\left(\frac{N(k_B\Theta_D)^2}{k_BT}\right),
\qquad\boxed{C_V=\left(\frac{\partial U}{\partial T}\right)_V\longrightarrow3Nk_B.}
$$
Equivalently, differentiating before integrating gives
$$
C_V=9Nk_B\left(\frac{T}{\Theta_D}\right)^3\int_0^{\Theta_D/T}\frac{x^4e^x}{(e^x-1)^2}\,dx.
$$
The integrand tends to $x^2$, and its integral is $(\Theta_D/T)^3/3$, producing the same <Dulong-Petit law>. In classical language, the <equipartition theorem> assigns $k_BT/2$ each to the kinetic and elastic parts of every harmonic <normal mode>. Fixed <zero-point energy> does not contribute to the <heat capacity>.