The phonon density of states tells us how many phonon modes there are at every frequency. This can be calculated by choosing a uniform grid of $\vec{k}$ states in the first Brillouin zone and calculating the frequencies of the phonon modes for each $\vec{k}$ state. A histogram is then constructed out of these frequencies to plot the density of states. The number of phonon modes depends on the size of the crystal; the more atoms there are, the more normal modes there are. Usually the density of states is specified per m³.

For an fcc lattice, the Brillouin zone is a truncated octahedron. There are six square faces and eight hexagonal faces.

$D(\omega)$
	$\large \sqrt{\frac{M}{C}}\omega$

At low frequencies, the density of states rises like $\omega^2$. This is the form that is expected for a linear dispersion relation, $\omega = c_s|\vec{k}|$, where $c_s$ is the speed of sound. There is a kink in the density of states at $\sqrt{\frac{M}{C}}\omega = \sqrt{2}$. At this frequency, transverse acoustic modes touch the Brillouin zone boundary at $L$. Further kinks are observed at $\sqrt{\frac{M}{C}}\omega = 2,\,\sqrt{6}, $ and $2\sqrt{2}$ where the phonon modes are at a symmetry point of the Brillouin zone boundary.

There is still some numerical noise visible in the density of states because the $\vec{k}$ state grid is not fine enough. To get a smooth density of states, it would be necessary to calculate more points than is possible in the time it takes to load a web page. Notice that the number of numerical errors is listed. The code that calculates the eigenvalues of the matrix sometimes fails. This is common for simple routines that calculate eigenvalues. As long as there are not too many errors compared to the number of $\vec{k}$ states calculated, it should not have a large influence on the result.

Matlab code for performing this calculation more precisely can be found here.