Solution (source code)

= Solution

The <density of states> $g(E)$ counts one-electron states per energy interval, including both <electron spin> orientations. It may be a total or a per-volume quantity; $U$ must have the matching normalization. The <Fermi-Dirac distribution>
$$
f(E,T)=\frac1{e^{(E-\mu(T))/(k_BT)}+1}
$$
is the occupation probability of such a state. In the fixed-volume <free electron model>, the integration variable $E$ and the <density of states> do not depend on $T$. The occupation depends explicitly on $T$ and implicitly through the <chemical potential> $\mu(T)$ needed to preserve electron number. Thermal expansion or a changing <electronic band structure> would require additional effects outside this model.

At $T\ll T_F$, only states within order $k_BT$ of the <Fermi energy> change their occupation. There are order $g(E_F)k_BT$ such states, each carrying an excitation energy of order $k_BT$. Thus $U(T)-U(0)\propto g(E_F)(k_BT)^2$ and the <electronic heat capacity> is proportional to $g(E_F)k_B^2T$. Using the supplied coefficient,
$$
\boxed{C_V=\gamma T,\qquad\gamma=\frac{\pi^2}{3}k_B^2g(E_F),\qquad
U(T)-U(0)=\frac12\gamma T^2.}
$$
The shell argument fixes the scaling, not the numerical coefficient; the controlled <Sommerfeld expansion> fixes that coefficient at fixed electron number. Filled states far below $E_F$ cannot accept arbitrary thermal excitations because of the <Pauli exclusion principle>.