Why doesn't the phonon theory reproduce a Debye-like curve for $C_V$ vs. $T$? 
*

*Debye's theory is a model of specific heat in which a high-frequency cut-off $\omega_D$ is put in by hand. The resulting curve of $C_V$ versus temperature $T$ grows as $T^3$ at low $T$ and saturates to a constant value at high $T$. The ad hoc parameter $\omega_D$ is then fitted experimentally. 

*The phonon theory also enables us to calculate the contribution of specific heat coming from phonons. It uses the dispersion relation $\omega=c_sk$. But that only gives rise to the $T^3$ behaviour which is only the low-temperature behaviour. See Ashcroft & Mermin, page 457 for a derivation.
Why does the phonon theory of solids fail to reproduce a Debye-like curve? In case I am ill-informed and it can reproduce a Debye-like curve, please suggest me a reference (preferably a book/lecture note not some research article).
 A: At the top of page 455 in Ashcroft and Mermin, they note that one can expand the Bose-Einstein distribution in powers of $1/T$ to get the high-temperature expansion.  The first term yields the Law of Dulong and Petit, and the rest of the terms decay as a function of $T$. Due to the natural high-frequency cut-off caused by the finite distance between atoms, the integrals in each of these terms is also finite. These two facts render the heat capacity finite at high temperatures and that the heat capacity matches the Law of Dulong and Petit as $T\to\infty$. For a specific case, consider the following.
1D harmonic model
To the extent that we can treat the vibrations of the lattice as harmonic (i.e. no non-linear oscillator terms), the phonon theory for a 1D solid yields a Debye-like curve. 
Dispersion relation and density of states
The dispersion relation for a 1D chain of atoms of mass $m$ separated by a distance $a$ and connected by springs of spring-constant $c$ is
$$
\omega = 2\sqrt{\frac{c}{m}}\left|\sin\left(\frac{k_xa}{2}\right)\right|,~~~~~~~\mbox{where}~~~~~~-\frac{\pi}{a}\leq k_x\leq \frac{\pi}{a}.
$$
Accordingly, the (1D) density of states is given by
$$
D(\omega) = \frac{L}{\pi}\frac{1}{d\omega/dk},
$$
where $L=Na$ is the total length of the solid, and $N$ is the number of atoms (also the number of primitive unit cells, since there is one atom per unit cell). Using the dispersion relation above, this becomes$$
D(\omega)=
\frac{N}{\omega_{\textrm{max}}}\frac{2}{\pi}
\begin{cases}
\frac{1}{\sqrt{1-(\omega/\omega_{\textrm{max}})^2}} & \omega \leq \omega_{\textrm{max}} \\
0 & \textrm{otherwise}
\end{cases},
$$
where $\omega_{\textrm{max}}=2\sqrt{c/m}$. Importantly, the cutoff frequency $\omega_{\textrm{max}}$ is built in automatically: the maximum frequency (minimum wavelength) emerges due to the finite spacing between atoms.
Integral expression for the internal energy
Upon converting the integral over $k$ to one over $\omega$, the internal energy becomes
$$
U = \int_0^{\infty}d\omega~D(\omega)\frac{\hbar\omega}{e^{\hbar\omega/k_{\textrm{B}}T}-1}
= \int_0^{\omega_{\textrm{max}}}d\omega~\frac{N}{\omega_{\textrm{max}}}\frac{2}{\pi}\frac{1}{\sqrt{1-(\omega/\omega_{\textrm{max}})^2}}\frac{\hbar\omega}{e^{\hbar\omega/k_{\textrm{B}}T}-1}.
$$
Upon changing variables to $x=\omega/\omega_{\textrm{max}}$, this becomes
$$
U= \frac{2N}{\pi}\hbar\omega_{\textrm{max}}\int_0^{1}dx~\frac{1}{\sqrt{1-x^2}}\frac{x}{e^{x(\hbar\omega_{\textrm{max}}/k_{\textrm{B}}T)}-1}.
$$
Heat capacity
The heat capacity is given by
$$
\frac{\partial U}{\partial T} = k_{\textrm{B}}\frac{2N}{\pi}\left(\frac{\hbar\omega_{\textrm{max}}}{k_{\textrm{B}}T}\right)^2\int_0^{1}dx~\frac{1}{\sqrt{1-x^2}}\frac{x^2e^{x(\hbar\omega_{\textrm{max}}/k_{\textrm{B}}T)}}{\left(e^{x(\hbar\omega_{\textrm{max}}/k_{\textrm{B}}T)}-1\right)^2}.
$$
This integral is finite for all positive values of $T$. Below, I've plotted this heat capacity (solid) along with the corresponding Debye interpolation (dashed) as a function of $T$. We can see that both follow the Law of Dulong and Petit at high temperatures, that they are linear at low temperatures (and agree), and that the phonon model yields a higher heat capacity at intermediate temperatures due to the rounding-off of the dispersion relation (i.e. more modes at lower frequencies means the higher frequency modes "turn on" at lower temperatures than in the Debye model with the linear dispersion relation. (Note that in order to compare the two, I matched the long-wavelength speed of sound of the two models, which results in $\omega_{\textrm{D}} = \frac{\pi}{2}\omega_{\textrm{max}}$.)

