Statistical Mechanics of Lattice Vibrations

from Statistical Physics using Mathematica
© James J. Kelly, 1996-2002

Introduction

The contribution of lattice vibrations to the heat capacity of a crystalline solid is determined by its spectrum of normal modes.  For a system of N atoms there will be 3  N normal modes.  The simplest approximation, due to Einstein, is to assume that all modes have the same frequency.   A better approximation due to Debye assumes that the density of states is proportional to the square of the frequency up to a maximum frequency determined by the total number of modes.  This assumption for the excitation spectrum (density of states) is often called the Debye interpolation formula because it provides a method for interpolating the temperature dependence between low and high temperature limits which can be derived from rather general considerations.  This model of the excitation spectrum is fairly accurate for long wavelengths, but the behavior of real crystals at short wavelengths is much more complicated.  Nevertheless, the Debye interpolation formula remains useful because many thermodynamic properties, such as the specific heat, are not terribly sensitive to details of the excitation spectrum.

The quantized normal modes of vibration are represented by phonons which obey the dispersion relation ω  =  k v, where k is the wave number and v is the propagation velocity.  If we assume that the sound velocity is independent of wave length and that the system is isotropic, the density of states g[k]  dk  =  
        V 4  π (k^2) d(k/(((2  π))^3)) can be expressed in terms of frequency as  g[ω]  dω  =  
        (Subscript[A, ω]) (ω^2) dω, where  (Subscript[A, ω]) is a normalization factor.  The normalization of the density of states is determined by the condition that the total number of states below the cut-off, (Subscript[k, D]) or (Subscript[ω, D]), must be equal to the total number of normal modes, 3  N, such that

\(∫\_0\%\(k\_D\)V\ \(\(k\^2\) dk\)\/\((2  π)\)\^3\  = \ \(3  N\ \ ⟹\ \ k\_D\^3 = 
        18 \( π\^2\) N\/V\)\)
\(∫\_0\%\(ω\_D\)\(A\_ω\) \(ω\^2\) dω\  = \ \(3 
           N\ \ ⟹\ \ A\_ω = \(\(9  N\)\/ω\_D\^3 = \(9  N\)\/\((\(k\_D\) v)\)\^3\)\)\)

For a crystal with a simple, nearly isotropic, structure, we can classify the vibrational modes into two types: longitudinal waves with density variations parallel to the direction of propagation or transverse waves with density variations perpendicular to the direction of propagation.  We assume that the sound velocity is independent of frequency but permit different velocities for longitudinal and transverse modes so that the dispersion relations for these types become (Subscript[ω, ℓ])  =  k (Subscript[v, ℓ]) and (Subscript[ω, t])  =  k (Subscript[v, t]).  Since there are two independent transverse polarization directions but only one longitudinal polarization direction, transverse modes should be normalized to 2  N and longitudinal to N.  Provided that the propagation speeds are not drastically different, it is useful to interpret  \(TraditionalForm\`\(\(ω\_D\  = \ \(k\_\(\(D\)\(\ \)\)\) 
          v\&_\)\(\ \)\)\) as a weighted average over polarizations, such that

\(3\/v\&_\^3 = 
      1\/v\&_\_ℓ\^3 + \ 2\/v\&_\_t\^3\)

where \(TraditionalForm\`v\&_\) is the weighted average propagation speed.  Therefore, combining these results, we find that the Debye frequency,

\(ω\_D = \(\((18 \( π\^2\) N\/V)\)\^\(1/3\)\) 
        v\&_\)

is proportional to the cube root of the density or the average interatomic spacing.  The density dependence of (Subscript[ω, D]) can be used to estimate the expansivity of the crystal.

Glossary:
	𝒩	= number of atoms
	V	= volume of crystal
	k	= wave number
	ω	= vibrational frequency
	v	= speed of sound (averaged over polarization)
	n[y]	= phonon occupation probability
	g[k] 	= density of states in terms of wave number
	g[ω] 	= density of states in terms of frequency
	(Subscript[ω, D])	= Debye frequency (max. ω)
	(Subscript[k, B])	= Boltzmann constant
	T	= temperature
	y	= ℏ  ω / ((Subscript[k, B]) T) = ratio of phonon/thermal energies
	(Subscript[T, E])	= Einstein temperature
	(Subscript[T, D])	= Debye temperature
	U	= internal energy 
	(Subscript[C, V])	= heat capacity

Initialization

Defaults,  packages, and labels

ClearAll["Global`*"];
Off[General::spell, General::spell1];

$DefaultFont = {"Times", 12};
$TextStyle = {FontFamily -> "Times", FontSize -> 12, 
          FontSlant -> "Italic"};

Needs["Utilities`Notation`"];
Needs["Graphics`Master`"]

Symbolize[(Subscript[k, B])];Symbolize[(Subscript[C, V])];Symbolize[(Subscript[T, x])_];Symbolize[(Subscript[ω, x])_];Symbolize[(Subscript[Z, m])_];

SetAttributes[{(Subscript[k, B]), ℏ }, Constant]

\(\(label[C\_V]\  = \ 
        StyleForm[\*""\!\(C\_V\)/(3\!\(𝒩k\_B\))""];\)\)

Simplification assumptions

MyAssumptions = {m > 0, V > 0, (Subscript[k, B]) > 0, 
          T > 0, ℏ > 0, β > 0, 𝒩 > 0, 
          (Subscript[k, D]) > 0, (Subscript[ω, D]) > 0, (Subscript[T, D]) > 0, (Subscript[T, E]) > 0};

MySimplify = 
        Simplify[#, MyAssumptions] &;
MyFullSimplify = FullSimplify[#, MyAssumptions] &;

Rules for changing variables

βToT  =  {β -> 
            (1/((Subscript[k, B]) T))}; 
Ttoβ  =  {T -> (1/((Subscript[k, B]) β))};
yToω  =  {y -> ((ℏ ω)/((Subscript[k, B]) T))};
ωTox  =  {ω -> x (Subscript[ω, D])};
ωToT  =  {(Subscript[ω, D]) ->  (((Subscript[k, B]) (Subscript[T, D]))/ℏ) };
TDtoτ   =  {T  ->  τ (Subscript[T, D])}; 
TDtoy  =  {T  ->  (Subscript[T, D])/y};
TEtoτ   =  {T  ->  τ (Subscript[T, E])};  
    TEtoy  =  {T  ->  (Subscript[T, E])/y};

Density of states

The phonon occupation probability was derived and studied in the notebook hotherm.nb.  The density of states assumes that wave numbers are distributed isotropically in momentum space.  The normalization condition assumes that there is a maximum wave number that is based upon the interatomic spacing of the crystal.

n[y_] := (1/(Exp[y] - 1))

g[k]  =  V ((4  π (k^2))/(((2 π))^3)) ;

kDsol  =  Solve[(∫_0^((Subscript[k, D]))) 
                g[k] \[DifferentialD]k ==  3  𝒩,  
            (Subscript[k, D])][[3]]  // MySimplify

ωDsol  =  (Subscript[ω, D])  →  (Subscript[k, D]) v  /. 
        kDsol

g[ω]  =  
      A (ω^2)  /. Solve[(∫_0^((Subscript[ω, D])))A (ω^2) \[DifferentialD]ω  ==  3  𝒩,  
            A][[1]]

Einstein model

To better understand the Debye interpolation formula, it will be helpful to compare with the Einstein model.  The relevant results are quoted here.  For further information, refer to the notebook hotherm.nb.

uEinstein = 3 𝒩 (Subscript[k, B]) (Subscript[T, E]) n[((Subscript[T, E])/T)]

cvEinstein = (∂_T) uEinstein

Heat capacity for isotropic Debye crystal

The internal energy is obtained by integrating the energy for each mode times its occupation probability against the density of states.  Although this integral can be performed in closed form, the PolyLog function is probably not familiar.  Therefore, we compare the high and low temperature behaviors of the Debye model with those of the simpler Einstein model below.  The following section then compares the Debye and Einstein models graphically.

uDebye = ((∫_0^((Subscript[ω, D])))n[((ℏ ω)/((Subscript[k, B]) T))] ℏ ω g[ω] \[DifferentialD]ω)  /. ωToT // Simplify // 
        Collect[#, {𝒩, (Subscript[k, B]), T}] &

cvDebye = (∂_T) uDebye  // Simplify // 
        Collect[#, {𝒩, (Subscript[k, B]), T}] &

The prediction that the specific heat capacity depends only upon a simple reduced temperature variable, τ, is an example of the law of corresponding states in which the thermodynamic functions assume universal forms.  This principle often applies, to a reasonable approximation, when the dynamics of a system are governed by a single identifiable scale, in this case the Debye temperature.  However, because the excitation spectrum for most crystals is more complicated than the Debye interpolation formula, real systems exhibit departures from this simple model.  For example, the normal modes for layered systems may have quite different propagation speeds for vibrations within a layer compared with vibrations between layers.  In such cases there will be at least two temperature scales and the heat capacity will be a more complicated function of temperature.  Nevertheless, the Debye interpolation formula remains useful because the heat capacity is not particularly sensitive to the details of the excitation spectrum.

cvDebye  /. TDtoτ  // Simplify

Limiting behaviors

High temperature limits

Series[cvEinstein /. TEtoy, {y, 0, 2}]  // 
      Collect[#1, {𝒩, (Subscript[k, B])}] &

Series[cvDebye /. TDtoy, {y, 0, 2}]  // 
        Simplify//Collect[#1, {𝒩, 
            (Subscript[k, B])}] &

Low temperature limits

To evaluate the low temperature behavior of the Einstein heat capacity, we must handle the essential singularity in Exp[1/T] explicitly.

lowTrules = {Exp[((x_.  (Subscript[T, E]))/T)] -> (d^(-x)), 
          Exp[((x_.  (Subscript[T, D]))/T)] -> (d^(-x))};

cvEinsteinLowT = cvEinstein /. 
            lowTrules /. {d -> Exp[-((Subscript[T, E])/T)]} // 
        Series[#1, {d, 0, 1}] &

To evaluate the low temperature behavior of the Debye heat capacity, we return to the energy integral and recognize that the temperature appears only in the occupation number and that for low T the range of integration can be extended to infinity.   Therefore, we can differentiate under the integral and the resulting integral can then be evaluated in closed form.  A change of variables is used to isolate the temperature dependence.

integrand = (∂_T) 
              n[((ℏ ω)/((Subscript[k, B]) T))] ℏ ω g[ω] (ω/x) /. {ω -> x (Subscript[k, B]) T/ℏ}

cvDebyeLowT = Integrate[integrand, {x, 0, ∞}]

Compare Einstein and Debye heat capacities

Plot[
        Evaluate[{(cvEinstein/(3 𝒩 (Subscript[k, B]))), 
              (cvDebye/(3 𝒩 (Subscript[k, B])))} /. {(Subscript[T, E]) -> 1, 
              (Subscript[T, D]) -> 1}], {T, 0.015, 3}, Frame -> True, 
        GridLines -> Automatic, 
        PlotLabel -> "Heat Capacity for Lattice Vibrations", 
        FrameLabel -> {"reduced temperature",  label[(Subscript[C, V])]}, 
        PlotStyle -> {Dashing[{0.02, 0.02}], {}}, 
        PlotLegend -> {"Einstein", "Debye"}, 
        LegendPosition -> {0.2, -0.2}, 
        LegendSize -> {0.5, 0.2}];

Although the Einstein model does predict that the heat capacity decreases at low temperature, and hence is an important improvement over the equipartition prediction of Dulong and Petit, its exponential behavior is too precipitous.  This exponential behavior occurs because the excitation spectrum of the Einstein crystal has an energy gap which leads to a heat capacity with the characteristic two-state form when T ≪ (Subscript[T, E]).  The low temperature behavior of the Debye model is more realistic because its excitation spectrum is continuous — although the range of excitations which contribute strongly decreases as temperature decreases, modes are always available at energies less than (Subscript[k, B]) T such that the heat capacity exhibits a power law instead of an exponential suppression at low temperatures.  In particular, its prediction that at low temperature (Subscript[C, V])  ∝  (T^3) is in good agreement with the data for most solid electrical insulators.  Furthermore, the value of (Subscript[T, D]) obtained by fitting the temperature dependence of the specific heat is usually in good agreement with that obtained from elastic properties (sound speed and density).

Similar considerations can also be applied to liquids, with several important differences.  First, fluids do not support sheer stress and hence do not support transverse modes of vibration; therefore, sound waves in fluids are longitudinal and only that contribution to the normalization of the Debye frequency should be considered.  At low temperature the heat capacity is then \(TraditionalForm\`C\_V\  ≈ \(4  π\^4\)\/\(\(5\)\(\ \)\)\ N\ k\_B\ \((T\/T\_D)\)\^3\).  Second, there are other modes of excitations, such as vortices, that can mask the phonon contribution at modest temperatures.  Only liquid helium-4 remains liquid at sufficiently low temperature for the phonon contribution to clearly dominate such that the  (T^3) behavior becomes evident.  Nevertheless, one finds that the specific heat for liquid-4 predicted from its density and sound velocity is in good agreement with calorimetry for T  ≲  0.6 kelvin.

Grüneisen model

In addition to the lattice vibrations described by the Debye model, the energy of a crystal also depends upon the mean interatomic spacing.  This effect can be represented by a potential energy Φ[V] which depends upon the volume of the crystal.  Furthermore, the frequency for each mode also depends upon volume.  The density dependence of phonon energies is parametrized by

\(γ = \(\(-\(\(∂ln\ ω\)\/\(∂ln\ V\)\)\)\  = \ \(-\(V\/ω\)\) ∂ω\/∂V\  ≈ \ 1\/3\)\)

which is assumed to be constant over the relevant ranges of temperature and pressure and is known as the Grüneisen constant.  Thus, we can obtain the mechanical equation of state from the free energy

F[T, 
        V]  =  Φ  -  (Subscript[k, B]) 
          T ∫dω g[ω] 
              Log[(Subscript[Z, 1])[T, ω]]

where (Subscript[Z, 1])[T, ω] is the partition function for a single mode.  Furthermore, we assume that Φ[V] can be expanded around the equilibrium density, such that

\(Φ[
        V]\  ≈ \ Φ\_0\  + \ \((V - V\_0)\)\^2\/\(2\ κ\ V\_0\)\)

where (N/(Subscript[V, 0]))is the equilibrium density and κ is an elasticity parameter.  Therefore, this extension to the Debye model allows the temperature dependence of the expansivity of the solid to be related to its equilibrium elastic properties.

Initialization

Symbolize[(Subscript[F, 1])];Symbolize[(Subscript[U, 1])];Symbolize[(Subscript[p, 1])];Symbolize[(Subscript[V, 0])];Symbolize[(Subscript[Φ, 0])];

SetAttributes[{γ,  (Subscript[Φ, 0]),  (Subscript[V, 0]),   κ }, 
      Constant]

DtRules  =  {Dt[ω, 
              V] -> -γ ω / V,  
          Dt[β, T] -> -β/T, 
          Dt[β, V] -> Dt[β, T] Dt[V, T], 
          Dt[ω, T] -> Dt[ω, V] Dt[V, T]};

GrüneisenRules ={(Subscript[ω, D]) -> (Subscript[ω, 0]) (((((Subscript[V, 0])/V)))^γ), Φ[
                V] -> (Subscript[Φ, 0]) + ((((V - (Subscript[V, 0])))^2)/(2 κ (Subscript[V, 0]))),Φ'[V]->(((V - (Subscript[V, 0])))/( κ (Subscript[V, 0]))), Φ''[V]->(1/( κ (Subscript[V, 0])))};

Single mode quantities

(Subscript[Z, 1])[ω_]  = ( (∑_(n = 0)^∞) 
                Exp[-β (n  +  
                      (1/2)) ℏ ω])  // ExpToTrig // 
        Simplify

(Subscript[F, 1])[ω_]  =   -(Subscript[k, B]) T Log[ 
          (Subscript[Z, 1])[ω] ]

(Subscript[U, 1])[ω_]  =  -D[
          Log[ (Subscript[Z, 1])[ω] ],  β]

The pressure is obtained from the partial derivative of F wrt V for constant T.

\(p\_1[ω_]\  = \ \(\((\(-Dt[F\_1[ω], V]\)\  //. 
                DtRules)\) /. {Dt[T, V] -> 0, 	Dt[V, T] -> 0}\)\  /. 
          Ttoβ // Simplify\)

Thus, we find that the pressure and the energy density contributions for each mode are simply proportional to each other, such that (Subscript[p, 1])  =   γ  ((Subscript[U, 1])/V).  Given the assumption that the Grüneisen constant, γ, is independent of frequency, this same relationship applies to the entire phonon contribution.  In addition, there will be the contribution from Φ such that the equation of state becomes

\(p = \ \(-\(∂Φ\/∂V\)\)\  + \ γ\ U\/V\)

where U is the total phonon energy.  We need to differentiate this equation with respect to V in order to obtain the expansivity and for that we need compare derivatives of the phonon energy wrt V and T.

dU1dT  = (Dt[(Subscript[U, 1])[ω], T]  //. DtRules)  // 
        Simplify

dU1dV  =  (Dt[(Subscript[U, 1])[ω], V]  //. DtRules)  // 
        Simplify

CV1  =  dU1dT  /. {Dt[V, T] -> 0}

dU1dVconstantT  =  dU1dV  /. {Dt[V, T] -> 0}

dU1dVconstantT  -  (γ/V) (CV1 T  -  
                (Subscript[U, 1])[ω])  /. βToT  // Simplify

Therefore, we can establish the following rule.

dUdVrule  =  
        Dt[U, V, 
            Constants -> {N, 
                T}] -> (γ/V) (T (Subscript[C, V])  -  
                U);

Grüneisen-Debye equation of state and thermal expansivity

EquationOfState  =  
      p  ==  -(∂_V) Φ[
                V]  +  ((γ U)/V)

alphaEquation  = Dt[EquationOfState, T, 
              Constants -> {p, N}] //. {Dt[U, T, 
                  Constants -> {N, p}] -> 
                (Subscript[C, V])  +  α V Dt[U, V, Constants -> {N, T}], 
              Dt[V, T, Constants -> {N, p}]  -> α V} /. 
          dUdVrule // Simplify

alphaSolution  =  Solve[
            alphaEquation, α][[1]]  // Simplify

pRule  =  Solve[EquationOfState /. GrüneisenRules, 
          p][[1]]

URule  = Solve[EquationOfState /. GrüneisenRules, 
          U][[1]]

expansivity  = (α  /. alphaSolution  /. URule /. 
            GrüneisenRules)  // Simplify

compressibility  =  -1/ (V Dt[p /. pRule, V])  /.  
            Dt[U, V] -> (γ/V) (T (Subscript[C, V])  -  U)  /. 
          URule // Simplify

expansivity / compressibility  // Simplify

Therefore, we obtain a very simple relationship between the isobaric expansivity and the isothermal compressibility.  Finally, we employ a general thermodynamic relationship between the principal heat capacities to obtain the isobaric heat capacity for this model.

(Subscript[C, p])  =  
      (Subscript[C, V])  +  T V ((expansivity^2)/compressibility)  // 
        Simplify

In the limit that the pressure is small (compared to the internal forces which determine (Subscript[V, 0])), the isobaric and isochoric heat capacities are nearly equal.

(Subscript[C, p]) / (Subscript[C, V ]) /. {p -> 0, V -> (Subscript[V, 0])}  // 
        Series[#, {(Subscript[C, V]), 0, 1}] & // Simplify

(Subscript[C, p])  -  (Subscript[C, V ]) /. {p -> 0, V -> (Subscript[V, 0])}  // 
        Series[#, {(Subscript[C, V]), 0, 2}] & // Simplify

Therefore, recognizing that (Subscript[C, V])  ∝  (T^3) for small T, we conclude that (Subscript[C, p])  -  (Subscript[C, V])  ∝   (T^7) and is quite small.

Problems

ℏ (Subscript[ω, roton])  = Δ + ((((ℏ^2)(k-(Subscript[k, 0])))^2)/(2μ))