In the Sommerfeld model why can’t we obtain the heat capacity via the usual method (i.e. calculating E = 2_s \left( \frac{L}{2 \pi} \right)^3 \int \varepsilon(\mathbf{k}) n_{F}(\beta(\varepsilon-\mu)) \mathrm{d} \mathbf{k} and then using C=dE/dT)? I can imagine this is because of the Fermi-Dirac distribution taking on different shapes for T=0 and T>0. However, in the Debye model we circumvented a similar problem, where the density of states dropped to zero at some point, by just integrating up to \omega_D. Cant we do something similar here and just integrate up to \epsilon_F.
We certainly can do that. The math of approximating the resulting is a bit more complicated than with the Debye model—that’s all. The whole triangle thing is only about approximating the resulting integral.