Debye model

The Debye model is a solid-state equivalent of Planck's law of black body radiation, where one treats electromagnetic radiation as a photon gas confined in a vacuum space. The Debye model treats atomic vibrations as phonons confined in the solid’s volume. Most of the calculation steps are identical as both are examples of a massless Bose gas with linear dispersion relation.

Consider a cube of side ${\displaystyle L}$. From the particle in a box article, the resonating modes of the sonic disturbances inside the box (considering for now only those aligned with one axis) have wavelengths given by

${\displaystyle \lambda _{n}={2L \over n}\,,}$

where ${\displaystyle n}$ is an integer. The energy of a phonon is

${\displaystyle E_{n}\ =h\nu _{n}\,,}$

where ${\displaystyle h}$ is the Planck constant and ${\displaystyle \nu _{n}}$ is the frequency of the phonon. Making the approximation that the frequency is inversely proportional to the wavelength, we have

${\displaystyle E_{n}=h\nu _{n}={hc_{\rm {s}} \over \lambda _{n}}={hc_{s}n \over 2L}\,,}$

in which ${\displaystyle c_{s}}$ is the speed of sound inside the solid. In three dimensions we will use

${\displaystyle E_{n}^{2}={p_{n}^{2}c_{\rm {s}}^{2}}=\left({hc_{\rm {s}} \over 2L}\right)^{2}\left(n_{x}^{2}+n_{y}^{2}+n_{z}^{2}\right)\,,}$

in which ${\displaystyle p_{n}}$ is the magnitude of the three-dimensional momentum of the phonon.

The approximation that the frequency is inversely proportional to the wavelength (giving a constant speed of sound) is good for low-energy phonons but not for high-energy phonons (see the article on phonons). This disagreement is one of the limitations of the Debye model. It produces incorrect results at intermediate temperatures, whereas the results are exact at the low and high temperatures limits.

Let's now compute the total energy in the box,

${\displaystyle E=\sum _{n}E_{n}\,{\bar {N}}(E_{n})\,,}$

where ${\displaystyle {\bar {N}}(E_{n})}$ is the number of phonons in the box with energy ${\displaystyle E_{n}}$. In other words, the total energy is equal to the sum of energy multiplied by the number of phonons with that energy (in one dimension). In 3 dimensions we have:

${\displaystyle U=\sum _{n_{x}}\sum _{n_{y}}\sum _{n_{z}}E_{n}\,{\bar {N}}(E_{n})\,.}$

Here, the Debye model and Planck's law of black body radiation differ. Unlike electromagnetic photon radiation in a box, there is a finite number of phonon energy states because a phonon cannot have arbitrarily high frequencies. Its frequency is bounded by the medium of its propagation—the atomic lattice of the solid. Consider an illustration of a transverse phonon below.

It is reasonable to assume that the minimum wavelength of a phonon is twice the atom separation, as shown in the lower figure. There are ${\displaystyle N}$ atoms in a solid. Our solid is a cube, which means there are ${\displaystyle {\sqrt[{3}]{N}}}$ atoms per edge. Atom separation is then given by ${\displaystyle L/{\sqrt[{3}]{N}}}$, and the minimum wavelength is

${\displaystyle \lambda _{\rm {min}}={2L \over {\sqrt[{3}]{N}}}\,,}$

making the maximum mode number ${\displaystyle n}$ (infinite for photons)

${\displaystyle n_{\rm {max}}={\sqrt[{3}]{N}}\,.}$

This number bounds the upper limit of the triple energy sum

${\displaystyle U=\sum _{n_{x}}^{\sqrt[{3}]{N}}\sum _{n_{y}}^{\sqrt[{3}]{N}}\sum _{n_{z}}^{\sqrt[{3}]{N}}E_{n}\,{\bar {N}}(E_{n})\,.}$

For slowly varying, well-behaved functions, a sum can be replaced with an integral (also known as Thomas–Fermi approximation)

${\displaystyle U\approx \int _{0}^{\sqrt[{3}]{N}}\int _{0}^{\sqrt[{3}]{N}}\int _{0}^{\sqrt[{3}]{N}}E(n)\,{\bar {N}}\left(E(n)\right)\,dn_{x}\,dn_{y}\,dn_{z}\,.}$

So far, there has been no mention of ${\displaystyle {\bar {N}}(E)}$, the number of phonons with energy ${\displaystyle E\,.}$ Phonons obey Bose–Einstein statistics. Their distribution is given by the famous Bose–Einstein statistics formula

${\displaystyle \langle N\rangle _{BE}={1 \over e^{E/kT}-1}\,.}$

Because a phonon has three possible polarization states (one longitudinal, and two transverse which approximately do not affect its energy) the formula above must be multiplied by 3,

${\displaystyle {\bar {N}}(E)={3 \over e^{E/kT}-1}\,.}$

Actually one uses an effective sonic velocity ${\displaystyle c_{s}:=c_{\rm {eff}}}$, i.e. the Debye temperature ${\displaystyle T_{\rm {D}}}$ (see below) is proportional to ${\displaystyle c_{\rm {eff}}}$, more precisely ${\displaystyle T_{\rm {D}}^{-3}\propto c_{\rm {eff}}^{-3}:=(1/3)c_{\rm {long}}^{-3}+(2/3)c_{\rm {trans}}^{-3}}$, where one distinguishes longitudinal and transversal sound-wave velocities (contributions 1/3 and 2/3, respectively). The Debye temperature or the effective sonic velocity is a measure of the hardness of the crystal.

Substituting into the energy integral yields

${\displaystyle U=\int _{0}^{\sqrt[{3}]{N}}\int _{0}^{\sqrt[{3}]{N}}\int _{0}^{\sqrt[{3}]{N}}E(n)\,{3 \over e^{E(n)/kT}-1}\,dn_{x}\,dn_{y}\,dn_{z}\,.}$

The ease with which these integrals are evaluated for photons is due to the fact that light's frequency, at least semi-classically, is unbound. As the figure above illustrates, this is not true for phonons. In order to approximate this triple integral, Debye used spherical coordinates.

${\displaystyle \ (n_{x},n_{y},n_{z})=(n\sin \theta \cos \phi ,n\sin \theta \sin \phi ,n\cos \theta )}$

and approximated the cube by an eighth of a sphere

${\displaystyle U\approx \int _{0}^{\pi /2}\int _{0}^{\pi /2}\int _{0}^{R}E(n)\,{3 \over e^{E(n)/kT}-1}n^{2}\sin \theta \,dn\,d\theta \,d\phi \,,}$

where ${\displaystyle R}$ is the radius of this sphere, which is found by conserving the number of particles in the cube and in the eighth of a sphere. The volume of the cube is ${\displaystyle N}$ unit-cell volumes,

${\displaystyle N={1 \over 8}{4 \over 3}\pi R^{3}\,,}$

so we get

${\displaystyle R={\sqrt[{3}]{6N \over \pi }}\,.}$

The substitution of integration over a sphere for the correct integral introduces another source of inaccuracy into the model.

The energy integral becomes

${\displaystyle U={3\pi \over 2}\int _{0}^{R}\,{hc_{s}n \over 2L}{n^{2} \over e^{hc_{\rm {s}}n/2LkT}-1}\,dn}$.

Changing the integration variable to ${\displaystyle x={hc_{\rm {s}}n \over 2LkT}}$,

${\displaystyle U={3\pi \over 2}kT\left({2LkT \over hc_{\rm {s}}}\right)^{3}\int _{0}^{hc_{\rm {s}}R/2LkT}{x^{3} \over e^{x}-1}\,dx.}$

To simplify the appearance of this expression, define the Debye temperature ${\displaystyle T_{\rm {D}}}$

${\displaystyle T_{\rm {D}}\ {\stackrel {\mathrm {def} }{=}}\ {hc_{\rm {s}}R \over 2Lk}={hc_{\rm {s}} \over 2Lk}{\sqrt[{3}]{6N \over \pi }}={hc_{\rm {s}} \over 2k}{\sqrt[{3}]{{6 \over \pi }{N \over V}}}}$

where ${\displaystyle V}$ is the volume of the cubic box of side ${\displaystyle L}$.

Many references[2][3] describe the Debye temperature as merely shorthand for some constants and material-dependent variables. However, as shown below, ${\displaystyle kT_{\rm {D}}}$ is roughly equal to the phonon energy of the minimum wavelength mode, and so we can interpret the Debye temperature as the temperature at which the highest-frequency mode (and hence every mode) is excited.

Continuing, we then have the specific internal energy:

${\displaystyle {\frac {U}{Nk}}=9T\left({T \over T_{\rm {D}}}\right)^{3}\int _{0}^{T_{\rm {D}}/T}{x^{3} \over e^{x}-1}\,dx=3TD_{3}\left({T_{\rm {D}} \over T}\right)\,,}$

where ${\displaystyle D_{3}(x)}$ is the (third) Debye function.

Differentiating with respect to ${\displaystyle T}$ we get the dimensionless heat capacity:

${\displaystyle {\frac {C_{V}}{Nk}}=9\left({T \over T_{\rm {D}}}\right)^{3}\int _{0}^{T_{\rm {D}}/T}{x^{4}e^{x} \over \left(e^{x}-1\right)^{2}}\,dx\,.}$

These formulae treat the Debye model at all temperatures. The more elementary formulae given further down give the asymptotic behavior in the limit of low and high temperatures. As already mentioned, this behavior is exact, in contrast to the intermediate behavior. The essential reason for the exactness at low and high energies, respectively, is that the Debye model gives (i) the exact dispersion relation ${\displaystyle E(\nu )}$ at low frequencies, and (ii) corresponds to the exact density of states ${\textstyle (\int g(\nu )\,d\nu \equiv 3N)}$, concerning the number of vibrations per frequency interval.

Debye derived his equation somewhat differently and more simply. Using continuum mechanics, he found that the number of vibrational states with a frequency less than a particular value was asymptotic to

${\displaystyle n\sim {1 \over 3}\nu ^{3}VF\,,}$

in which ${\displaystyle V}$ is the volume and ${\displaystyle F}$ is a factor which he calculated from elasticity coefficients and density. Combining this formula with the expected energy of a harmonic oscillator at temperature T (already used by Einstein in his model) would give an energy of

${\displaystyle U=\int _{0}^{\infty }\,{h\nu ^{3}VF \over e^{h\nu /kT}-1}\,d\nu \,,}$

if the vibrational frequencies continued to infinity. This form gives the ${\displaystyle T^{3}}$ behaviour which is correct at low temperatures. But Debye realized that there could not be more than ${\displaystyle 3N}$ vibrational states for N atoms. He made the assumption that in an atomic solid, the spectrum of frequencies of the vibrational states would continue to follow the above rule, up to a maximum frequency ${\displaystyle \nu _{m}}$chosen so that the total number of states is

${\displaystyle 3N={1 \over 3}\nu _{m}^{3}VF\,.}$

Debye knew that this assumption was not really correct (the higher frequencies are more closely spaced than assumed), but it guarantees the proper behaviour at high temperature (the Dulong–Petit law). The energy is then given by

${\displaystyle {\begin{aligned}U&=\int _{0}^{\nu _{m}}\,{h\nu ^{3}VF \over e^{h\nu /kT}-1}\,d\nu \,,\\&=VFkT(kT/h)^{3}\int _{0}^{T_{\rm {D}}/T}\,{x^{3} \over e^{x}-1}\,dx\,.\end{aligned}}}$

Substituting ${\displaystyle T_{\rm {D}}}$ for ${\displaystyle h\nu _{m}/k}$,

${\displaystyle {\begin{aligned}U&=9NkT(T/T_{\rm {D}})^{3}\int _{0}^{T_{\rm {D}}/T}\,{x^{3} \over e^{x}-1}\,dx\,,\\&=3NkTD_{3}(T_{\rm {D}}/T)\,,\end{aligned}}}$

where ${\displaystyle D_{3}}$ is the function later given the name of third-order Debye function.

First we derive the vibrational frequency distribution; the following derivation is based on Appendix VI from.[4] Consider a three-dimensional isotropic elastic solid with N atoms in the shape of a rectangular parallelepiped with side-lengths ${\displaystyle L_{x},L_{y},L_{z}}$. The elastic wave will obey the wave equation and will be plane waves; consider the wave vector ${\displaystyle \mathbf {k} =(k_{x},k_{y},k_{z})}$ and define ${\displaystyle l_{x}={\frac {k_{x}}{|\mathbf {k} |}},l_{y}={\frac {k_{y}}{|\mathbf {k} |}},l_{z}={\frac {k_{z}}{|\mathbf {k} |}}}$. Note that we have

${\displaystyle l_{x}^{2}+l_{y}^{2}+l_{z}^{2}=1.}$

(1)

Solutions to the wave equation are

${\displaystyle u(x,y,z,t)=\sin(2\pi \nu t)\sin \left({\frac {2\pi l_{x}x}{\lambda }}\right)\sin \left({\frac {2\pi l_{y}y}{\lambda }}\right)\sin \left({\frac {2\pi l_{z}z}{\lambda }}\right)}$

and with the boundary conditions ${\displaystyle u=0}$ at ${\displaystyle x,y,z=0,x=L_{x},y=L_{y},z=L_{z}}$, we have

${\displaystyle {\frac {2l_{x}L_{x}}{\lambda }}=n_{x};{\frac {2l_{y}L_{y}}{\lambda }}=n_{y};{\frac {2l_{z}L_{z}}{\lambda }}=n_{z}}$

(2)

where ${\displaystyle n_{x},n_{y},n_{z}}$ are positive integers. Substituting (2) into (1) and also using the dispersion relation ${\displaystyle c_{s}=\lambda \nu }$, we have

${\displaystyle {\frac {n_{x}^{2}}{(2\nu L_{x}/c_{s})^{2}}}+{\frac {n_{y}^{2}}{(2\nu L_{y}/c_{s})^{2}}}+{\frac {n_{z}^{2}}{(2\nu L_{z}/c_{s})^{2}}}=1.}$

The above equation, for fixed frequency ${\displaystyle \nu }$, describes an eighth of an ellipse in "mode space" (an eighth because ${\displaystyle n_{x},n_{y},n_{z}}$ are positive). The number of modes with frequency less than ${\displaystyle \nu }$ is thus the number of integral points inside the ellipse, which, in the limit of ${\displaystyle L_{x},L_{y},L_{z}\to \infty }$ (i.e. for a very large parallelepiped) can be approximated to the volume of the ellipse. Hence, the number of modes ${\displaystyle N(\nu )}$ with frequency in the range ${\displaystyle [0,\nu ]}$ is

${\displaystyle N(\nu )={\frac {1}{8}}{\frac {4\pi }{3}}\left({\frac {2\nu }{c_{\mathrm {s} }}}\right)^{3}L_{x}L_{y}L_{z}={\frac {4\pi \nu ^{3}V}{3c_{\mathrm {s} }^{3}}},}$

(3)

where ${\displaystyle V=L_{x}L_{y}L_{z}}$ is the volume of the parallelepiped. Note that the wave speed in the longitudinal direction is different from the transverse direction and that the waves can be polarised one way in the longitudinal direction and two ways in the transverse direction; thus we define ${\displaystyle {\frac {3}{c_{s}^{3}}}={\frac {1}{c_{\text{long}}^{3}}}+{\frac {2}{c_{\text{trans}}^{3}}}}$.

Following the derivation from,[5] we define an upper limit to the frequency of vibration ${\displaystyle \nu _{D}}$; since there are N atoms in the solid, there are 3N quantum harmonic oscillators (3 for each x-, y-, z- direction) oscillating over the range of frequencies ${\displaystyle [0,\nu _{D}]}$. Hence we can determine ${\displaystyle \nu _{D}}$ using

${\displaystyle 3N=N(\nu _{\rm {D}})={\frac {4\pi \nu _{\rm {D}}^{3}V}{3c_{\rm {s}}^{3}}}}$.

(4)

By defining ${\displaystyle \nu _{\rm {D}}={\frac {kT_{\rm {D}}}{h}}}$, where k is the Boltzmann constant and h is the Planck constant, and substituting (4) into (3), we get

${\displaystyle N(\nu )={\frac {3Nh^{3}\nu ^{3}}{k^{3}T_{\rm {D}}^{3}}},}$

(5)

this definition is more standard. We can find the energy contribution for all oscillators oscillating at frequency ${\displaystyle \nu }$. Quantum harmonic oscillators can have energies ${\displaystyle E_{i}=(i+1/2)h\nu }$ where ${\displaystyle i=0,1,2,\dotsc }$ and using Maxwell-Boltzmann statistics, the number of particles with energy ${\displaystyle E_{i}}$ is

${\displaystyle n_{i}={\frac {1}{A}}e^{-E_{i}/(kT)}={\frac {1}{A}}e^{-(i+1/2)h\nu /(kT)}.}$

The energy contribution for oscillators with frequency ${\displaystyle \nu }$ is then

${\displaystyle dU(\nu )=\sum _{i=0}^{\infty }E_{i}{\frac {1}{A}}e^{-E_{i}/(kT)}}$.

(6)

By noting that ${\displaystyle \sum _{i=0}^{\infty }n_{i}=dN(\nu )}$ (because there are ${\displaystyle dN(\nu )}$ modes oscillating with frequency ${\displaystyle \nu }$), we have

${\displaystyle {\frac {1}{A}}e^{-1/2h\nu /(kT)}\sum _{i=0}^{\infty }e^{-ih\nu /(kT)}={\frac {1}{A}}e^{-1/2h\nu /(kT)}{\frac {1}{1-e^{-h\nu /(kT)}}}=dN(\nu ).}$

From above, we can get an expression for 1/A; substituting it into (6), we have

${\displaystyle {\begin{aligned}dU&=dN(\nu )e^{1/2h\nu /(kT)}(1-e^{-h\nu /(kT)})\sum _{i=0}^{\infty }h\nu (i+1/2)e^{-h\nu (i+1/2)/(kT)}\\\\&=dN(\nu )(1-e^{-h\nu /(kT)})\sum _{i=0}^{\infty }h\nu (i+1/2)e^{-h\nu i/(kT)}\\&=dN(\nu )h\nu \left({\frac {1}{2}}+(1-e^{-h\nu /(kT)})\sum _{i=0}^{\infty }ie^{-h\nu i/(kT)}\right)\\&=dN(\nu )h\nu \left({\frac {1}{2}}+{\frac {1}{e^{h\nu /(kT)}-1}}\right).\end{aligned}}}$

Integrating with respect to ν yields

${\displaystyle U={\frac {9Nh^{4}}{k^{3}T_{\rm {D}}^{3}}}\int _{0}^{\nu _{D}}\left({\frac {1}{2}}+{\frac {1}{e^{h\nu /(kT)}-1}}\right)\nu ^{3}d\nu .}$

The temperature of a Debye solid is said to be low if ${\displaystyle T\ll T_{\rm {D}}}$, leading to

${\displaystyle {\frac {C_{V}}{Nk}}\sim 9\left({T \over T_{\rm {D}}}\right)^{3}\int _{0}^{\infty }{x^{4}e^{x} \over \left(e^{x}-1\right)^{2}}\,dx.}$

This definite integral can be evaluated exactly:

${\displaystyle {\frac {C_{V}}{Nk}}\sim {12\pi ^{4} \over 5}\left({T \over T_{\rm {D}}}\right)^{3}.}$

In the low-temperature limit, the limitations of the Debye model mentioned above do not apply, and it gives a correct relationship between (phononic) heat capacity, temperature, the elastic coefficients, and the volume per atom (the latter quantities being contained in the Debye temperature).

The temperature of a Debye solid is said to be high if ${\displaystyle T\gg T_{\rm {D}}}$. Using ${\displaystyle e^{x}-1\approx x}$ if ${\displaystyle |x|\ll 1}$ leads to

${\displaystyle {\frac {C_{V}}{Nk}}\sim 9\left({T \over T_{\rm {D}}}\right)^{3}\int _{0}^{T_{\rm {D}}/T}{x^{4} \over x^{2}}\,dx}$

which upon integration gives

${\displaystyle {\frac {C_{V}}{Nk}}\sim 3\,.}$

This is the Dulong–Petit law, and is fairly accurate although it does not take into account anharmonicity, which causes the heat capacity to rise further. The total heat capacity of the solid, if it is a conductor or semiconductor, may also contain a non-negligible contribution from the electrons.

Debye vs. Einstein. Predicted heat capacity as a function of temperature.

The Debye and Einstein models correspond closely to experimental data, but the Debye model is correct at low temperatures whereas the Einstein model is not. To visualize the difference between the models, one would naturally plot the two on the same set of axes, but this is not immediately possible as both the Einstein model and the Debye model provide a functional form for the heat capacity. As models, they require scales to relate them to their real-world counterparts. One can see that the scale of the Einstein model is given by ${\displaystyle \epsilon /k}$:

${\displaystyle C_{V}=3Nk\left({\epsilon \over kT}\right)^{2}{e^{\epsilon /kT} \over \left(e^{\epsilon /kT}-1\right)^{2}}.}$

The scale of the Debye model is ${\displaystyle T_{\rm {D}}}$, the Debye temperature. Both are usually found by fitting the models to the experimental data. (The Debye temperature can theoretically be calculated from the speed of sound and crystal dimensions.) Because the two methods approach the problem from different directions and different geometries, Einstein and Debye scales are not the same, that is to say

${\displaystyle {\epsilon \over k}\neq T_{\rm {D}}\,,}$

which means that plotting them on the same set of axes makes no sense. They are two models of the same thing, but of different scales. If one defines the Einstein condensation temperature as

${\displaystyle T_{\rm {E}}\ {\stackrel {\mathrm {def} }{=}}\ {\epsilon \over k}\,,}$

then one can say

${\displaystyle T_{\rm {E}}\neq T_{\rm {D}}\,,}$

and, to relate the two, we must seek the ratio

${\displaystyle {\frac {T_{\rm {E}}}{T_{\rm {D}}}}\,.}$

The Einstein solid is composed of single-frequency quantum harmonic oscillators, ${\displaystyle \epsilon =\hbar \omega =h\nu }$. That frequency, if it indeed existed, would be related to the speed of sound in the solid. If one imagines the propagation of sound as a sequence of atoms hitting one another, then it becomes obvious that the frequency of oscillation must correspond to the minimum wavelength sustainable by the atomic lattice, ${\displaystyle \lambda _{min}}$, where

${\displaystyle \nu ={c_{\rm {s}} \over \lambda }={c_{\rm {s}}{\sqrt[{3}]{N}} \over 2L}={c_{\rm {s}} \over 2}{\sqrt[{3}]{N \over V}}}$,

which makes the Einstein temperature ${\displaystyle T_{\rm {E}}={\epsilon \over k}={h\nu \over k}={hc_{\rm {s}} \over 2k}{\sqrt[{3}]{N \over V}}\,,}$ and the sought ratio is therefore

${\displaystyle {T_{\rm {E}} \over T_{\rm {D}}}={\sqrt[{3}]{\pi \over 6}}\ =0.805995977...}$

Using the ratio, both models can be plotted on the same graph. This ratio is the cube root of the ratio of the volume of one octant of a 3-dimensional sphere to the volume of the cube that contains it, which is just the correction factor used by Debye when approximating the energy integral above. Alternatively, the ratio of the two temperatures can be seen to be the ratio of Einstein's single frequency at which all oscillators oscillate and Debye's maximum frequency. Einstein's single frequency can then be seen to be a mean of the frequencies available to the Debye model.

Even though the Debye model is not completely correct, it gives a good approximation for the low temperature heat capacity of insulating, crystalline solids where other contributions (such as highly mobile conduction electrons) are negligible. For metals, the electron contribution to the heat is proportional to ${\displaystyle T}$, which at low temperatures dominates the Debye ${\displaystyle T^{3}}$ result for lattice vibrations. In this case, the Debye model can only be said to approximate the lattice contribution to the specific heat. The following table lists Debye temperatures for several pure elements[2] and sapphire:

The Debye model's fit to experimental data is often phenomenologically improved by allowing the Debye temperature to become temperature dependent;[6] for example, the value for water ice increases from about 222 K[7] to 300 K[8] as the temperature goes from absolute zero to about 100 K.

For other bosonic quasi-particles, e.g., for magnons (quantized spin waves) in ferromagnets instead of the phonons (quantized sound waves), one can derive analogous results. In this case at low frequencies one has different dispersion relations of momentum and energy, e.g., ${\displaystyle E(\nu )\propto k^{2}}$ in the case of magnons, instead of ${\displaystyle E(\nu )\propto k}$ for phonons (with ${\displaystyle k=2\pi /\lambda }$). One also has different density of states (e.g., ${\displaystyle \int g(\nu ){\rm {d}}\nu \equiv N\,}$). As a consequence, in ferromagnets one gets a magnon contribution to the heat capacity, ${\displaystyle \Delta C_{\,{\rm {V|\,magnon}}}\,\propto T^{3/2}}$, which dominates at sufficiently low temperatures the phonon contribution, ${\displaystyle \,\Delta C_{\,{\rm {V|\,phonon}}}\propto T^{3}}$. In metals, in contrast, the main low-temperature contribution to the heat capacity, ${\displaystyle \propto T}$, comes from the electrons. It is fermionic, and is calculated by different methods going back to Sommerfeld's free electron model.

It was long thought that phonon theory is not able to explain the heat capacity of liquids, since liquids only sustain longitudinal, but not transverse phonons, which in solids are responsible for 2/3 of the heat capacity. However, Brillouin scattering experiments with neutrons and with X-rays, confirming an intuition of Yakov Frenkel,[9] have shown that transverse phonons do exist in liquids, albeit restricted to frequencies above a threshold called the Frenkel frequency. Since most energy is contained in these high-frequency modes, a simple modification of the Debye model is sufficient to yield a good approximation to experimental heat capacities of simple liquids.[10] More recently, it has been shown that instantaneous normal modes associated with relaxations from saddle points in the liquid energy landscape, which dominate the frequency spectrum of liquids at low frequencies, may determine the specific heat of liquids as a function of temperature over a broad range.[11]

The Debye frequency (Symbol: ${\displaystyle \omega _{\rm {Debye}}}$ or ${\displaystyle \omega _{\rm {D}}}$) is a parameter in the Debye model. It refers to a cut-off angular frequency for waves of a harmonic chain of masses, used to describe the movement of ions in a crystal lattice and more specifically, to correctly predict the heat capacity in such crystals to be constant for high temperatures (Dulong–Petit law). The term was first introduced by Peter Debye in 1912.[12]

Throughout this whole article periodic boundary conditions are assumed.

Relation to Debye's temperature

The Debye temperature ${\displaystyle \theta _{\rm {D}}}$, another parameter in Debye model, is related to the Debye frequency by the relation

$${\displaystyle \theta _{\rm {D}}={\frac {\hbar }{k_{\rm {B}}}}\omega _{\rm {D}},}$$

where ${\displaystyle \hbar }$ is the reduced Planck constant and ${\displaystyle k_{\rm {B}}}$is the Boltzmann constant.

Derivation with the actual dispersion relation

Because only the discretized points matter, two different waves could render the same physical manifestation (see Phonon).

This problem could be made more insightful by making it more complex. Instead of using the dispersion relation ${\displaystyle \omega =v_{\rm {s}}k}$, the correct dispersion relation is now going to be assumed. From classical mechanics it is known that for an equidistant chain of masses which interact harmonically with each other the dispersion relation reads as follows[15]

$${\displaystyle \omega (k)=2{\sqrt {\frac {\kappa }{m}}}\left|\sin \left({\frac {ka}{2}}\right)\right|.}$$

After plotting this relation, it is clear that Debye's estimation of the cut-off wavelength was right after all. Because for every wavenumber bigger than ${\displaystyle \pi /a}$ (that is: ${\displaystyle \lambda }$ is smaller than ${\displaystyle 2a}$) a wavenumber that is smaller than ${\displaystyle \pi /a}$ could be found with the same angular frequency. This means the resulting physical manifestation for the mode with the larger wavenumber is indistinguishable from the one with the smaller wavenumber. Thereby, the study of the dispersion relation can be limited to the first brillouin zone[16] i.e. for ${\textstyle k\in \left[-{\frac {\pi }{a}},{\frac {\pi }{a}}\right]}$.This is possible because the system consists of discretized points, as is demonstrated in the animated picture. Dividing the dispersion relation by ${\displaystyle k}$ and inserting ${\displaystyle \pi /a}$ for ${\displaystyle k}$, we find the speed of a wave with ${\displaystyle k=\pi /a}$ to be

$${\displaystyle v_{\rm {s}}(k=\pi /a)={\frac {2a}{\pi }}{\sqrt {\frac {\kappa }{m}}}.}$$

By simply inserting ${\displaystyle k=\pi /a}$ in the original dispersion relation we find

$${\displaystyle \omega (k=\pi /a)=2{\sqrt {\frac {\kappa }{m}}}=\omega _{\rm {D}}.}$$

Combining these results the same result is once again found

$${\displaystyle \omega _{\rm {D}}={\frac {\pi v_{\rm {s}}}{a}}.}$$

However, for diatomic chains (and more complex chains) the associated cut-off frequency (and wavelength) is not very accurate, since the cut-off wavelength is twice as big and the dispersion relation consists of two branches (for a diatomic chain). It is also not certain from this whether for more dimensional systems the cut-off frequency was accurately predicted by Debye.

Alternative derivation

The physical result of two waves can be identical when at least one of them has a wavelength that is bigger than twice the initial distance between the masses (taken from Nyquist–Shannon sampling theorem).

For a one dimensional chain this result could also be reproduced using theory on aliasing. The Nyquist–Shannon sampling theorem is used in the following derivation; the main difference being that in the following derivation the discretization is not in time, but in space. If we use the correct dispersion relation from last paragraph, it will be clear in another insightful way why the cut-off frequency has the value previously (twice) derived. So again,

$${\displaystyle \omega (k)=2{\sqrt {\frac {\kappa }{m}}}\left|\sin \left({\frac {ka}{2}}\right)\right|}$$

is assumed.

This derivation is completely equivalent to the previous one, that is: the same assumptions are made to retrieve the result. It is not more or less accurate, it is just a different approach.

To determine where the cut-off frequency should be, it is useful to first determine where the cut-off of the wavelength should be. From the dispersion relation we know that for ${\displaystyle k>\pi /a}$ every mode is repeated, so the cut-off wavelength would be at ${\displaystyle \lambda _{\rm {D}}=2a}$. From this and the periodic boundary conditions you can immediately see that the total number of modes per polarization would be ${\displaystyle N}$. As seen in the gif of the previous paragraph this is because every wave with a wavelength shorter than ${\displaystyle 2a}$ could be replaced by a wave with a wavelength longer than ${\displaystyle 2a}$ to regain the same physical result.

However, the dispersion relation from previous paragraph (the correct one) is not even necessary in reasoning as to why the cut-off should be at ${\displaystyle \lambda =2a}$. Because, as is depicted, only waves with a longer wavelength than ${\displaystyle 2a}$ could render the same physical result as another one. So this is another way to correctly predict the cut-off wavelength of phonons without using the correct dispersion relation (or even knowledge from classical mechanics as Debye did). However, using the wrong dispersion relation which Debye assumed, waves with a smaller wavelength would have a higher frequency, but the relative movement of the masses would be the same, so this does not render new modes.

This results again in ${\displaystyle k_{\rm {D}}=\pi /a}$, rendering

$${\displaystyle \omega _{\rm {D}}={\frac {\pi v_{\rm {s}}}{a}}.}$$

Also here it does not matter which dispersion relation is used (the correct one or the one Debye used), the same cut-off frequency would be found.

Unfortunately, the same method could not be used (as easily) for a two- or three-dimensional crystal, because diagonal waves would have a larger cut-off wavelength, which are also difficult to predict.

• Bose gas
• Gas in a box
• Grüneisen parameter
• Bloch–Grüneisen temperature
• Electrical resistivity and conductivity#Temperature dependence

• CRC Handbook of Chemistry and Physics, 56th Edition (1975–1976)
• Schroeder, Daniel V. An Introduction to Thermal Physics. Addison-Wesley, San Francisco (2000). Section 7.5.

• Experimental determination of specific heat, thermal and heat conductivity of quartz using a cryostat.
• Simon, Steven H. (2014) The Oxford Solid State Basics (most relevant ones: 1, 2 and 6)