Statistical Mechanics of Noninteracting Oscillators

Einstein model of lattice vibrations in a crystal

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

Introduction

There are many physical systems which can be represented as a collection of noninteracting oscillators.  For example, the Einstein model of lattice vibrations represents the quantized normal modes of vibration for a crystal  (phonons) using a single mean frequency.  In that case the number of normal  modes is 3 times the number of atoms in the crystal because each atom has three independent modes of vibration.  However, the normal modes of a real crystal exhibit a more complex frequency spectrum than postulated by the Einstein model.  Nevertheless, the Einstein model contains enough of the physics to explain some of the most important features of the thermodynamics of lattice vibrations.  The most important feature of this model is that the energy for lattice vibrations is quantized.  Consequently, one finds that the heat capacity approaches zero rapidly at very low temperature, in contradiction with the classical model for which the equipartition theorem predicts that the heat capacity is simply 3  N (Subscript[k, B]).  This result helped to gain support for the early quantum theory and hence the Einstein crystal is of considerable historical importance.  Furthermore, the model is both simple and robust enough to employ as a starting point for analyzing many other situations.  For example, more sophisticated models of lattice vibrations can be constructed by representing the spectrum of lattice vibrations as a collection of Einstein crystals.  Similarly, vibrations of linear molecules can be analyzed in terms of normal  modes.  Furthermore, the Planck model of electromagnetic radiation in an enclosure is also based upon quantized harmonic modes.  Therefore, we start with the simplest model of this type.  Consider a system of n noninteracting harmonic oscillators in thermal equilibrium.  It is important to recognize that n is the number of oscillator modes, not the number of atoms.

Glossary:
	n	= number of oscillator modes
	ϵ	= ℏ ω = energy of a single quantum
	(Subscript[k, B])	= Boltzmann constant
	T	= temperature
	m	= total number of quanta
	x	= average excitation per oscillator mode
	y	= ϵ / (Subscript[k, B]) T = ratio of excitation/thermal energies
	(Subscript[T, E])	= ℏ ω / (Subscript[k, B]) = Einstein temperature = characteristic temperature scale
	τ	= 1/y = T/(Subscript[T, E]) = reduced temperature
	N	= number of atoms
	S	= entropy
	F	= Helmholtz free energy
	U	= (m + (n/2)) ℏ ω = internal energy, including zero-point energy 
	(Subscript[C, V])	= heat capacity
	(Subscript[Z, 1])	= single-mode partition function
	Z	= canonical partition function

To express the thermodynamic functions in terms of different choices of independent variables with minimal symbolic changes, we employ different arguments.  For example, S[x] and S[y]  express the entropy in terms of either the relative excitation, x, or the inverse temperature, y, respectively.  Similarly, we use S[T] to express entropy in terms of temperature where T is expressed in dimensionless form as 1/y.

Initialization

Defaults, packages, and symbols

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

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

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

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

Establish  transformation rules

Stirling's formula

StirlingApprox = {Log[x_!] -> x Log[x] - x};

Rules for changing variables

yToT = {y -> (ϵ/((Subscript[k, B]) T))};
yToTau = {y -> (1/τ)};
TtoY = {T -> (ϵ/((Subscript[k, B]) y)), β ϵ -> y};

dydT = -(y/T);
dydB = (y/B);

Plot labels

label[U]  =   StyleForm["U/nϵ"];
label[F]  =   StyleForm["F/nϵ"];
\(\(label[S]\  = \ \ StyleForm[\*""S/\!\(nk\_B\)""];\)\)
\(\(label[C\_V]\  = \ 
        StyleForm[\*""\!\(C\_V\)/\!\(nk\_B\)""];\)\)

Microcanonical Ensemble

Multiplicity and entropy for large n

The multiplicity function for this system is determined by the number of ways that m indistinguishable objects (quanta) can be distributed among n distinguishable boxes (modes).  Suppose that the boxes are represented by vertical lines marking their boundaries and the objects are represented by circles, as sketched below for a particular configuration in which 15 quanta are distributed among 7 modes with the first mode containing 2 quanta, the second 4, the third none, the fourth 3, the fifth 5, the sixth none, and the seventh 1.

oo|oooo||ooo|ooooo||o

Note that the outermost walls need not be included.  Thus, the number of combinations of m objects and n - 1 interior partitions is simply (m + n - 1)!.  However, because the m! permutations of the balls or the (n - 1)! permutations of the partitions among themselves leave the system invariant, the degeneracy is simply the binomial coefficient (({{m + n - 1}, {m}})).

multiplicity = (((n + m - 1)!)/(m! (n - 1)!));

logMultiplicity = ((Log[multiplicity] // PowerExpand) /. 
                StirlingApprox //. {m -> x n, n - 1 -> n, 
                n + n x -> (1 + x) n}) // PowerExpand // 
        Simplify

S[x] = (Subscript[k, B]) logMultiplicity

plotSvsX  =  
        Plot[Evaluate[S[x]/(n (Subscript[k, B]))], {x, 0.01, 5}, 
          PlotLabel -> "Entropy", 
          PlotRange -> {{0, 5}, {0, 3}}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"x",  label[S]}, 
          DisplayFunction -> Identity];

Determine relationship between temperature and energy

eq1 = β == 
        Simplify[(∂_x) S[x]/(n (Subscript[k, B]) ϵ)]

xsol = Solve[eq1, x][[1]] //. 
        TtoY

U[y] = (x + (1/2)) n ϵ /. xsol  // ExpToTrig  // 
        Simplify

plotUvsY = 
        Plot[Evaluate[U[y]/(n ϵ)], {y, 0.01, 5}, 
          PlotLabel -> "Energy", PlotRange -> {{0, 5}, {0, 5}}, 
          GridLines -> None, Frame -> True, 
          FrameLabel -> {"y",  label[U]}, 
          DisplayFunction -> Identity];

Heat capacity

(Subscript[C, V])[y] = (∂_y) U[y] dydT /. TtoY

plotCVvsY = 
        Plot[Evaluate[(Subscript[C, V])[y]/(n (Subscript[k, B]))], {y, 0.01, 5}, 
          PlotLabel -> "Heat Capacity", 
          PlotRange -> {{0, 5}, {0, 1}}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"y",  label[(Subscript[C, V])]}, 
          DisplayFunction -> Identity];

Entropy and free energy as functions of inverse temperature

S[y]  = (S[x] /. 
                xsol /. {Log[
                    (1/(-1 + (E^y)))] -> -(y/2) - Log[2  Sinh[(y/2)]], 
                Log[1 + (1/(-1 + (E^y)))] -> 
                  (y/2) - Log[2  Sinh[(y/2)]]}) // ExpToTrig // 
        Simplify

plotSvsY = 
        Plot[Evaluate[S[y]/(n (Subscript[k, B]))], {y, 0.01, 5}, 
          PlotLabel -> "Entropy", 
          PlotRange -> {{0, 5}, {0, 3}}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"y",  label[S]}, 
          DisplayFunction -> Identity];

Helmholtz free energy

F[y] =  U[y] - T S[y] /. TtoY  // Simplify

plotFvsY = 
        Plot[Evaluate[F[y]/(n ϵ)], {y, 0.01, 5}, 
          PlotLabel -> "Free Energy", 
          PlotRange -> {{0, 5}, {-5, 1}}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"y",  label[F]}, 
          DisplayFunction -> Identity];

Display principal thermodynamic functions

Show[
        GraphicsArray[{{plotUvsY, plotFvsY}, {plotCVvsY, plotSvsY}}, 
          GraphicsSpacing -> {-0.1, 0.05}, 
          DisplayFunction -> $DisplayFunction]];

Examine thermodynamic functions of T

Here the thermodynamic functions are expressed in terms of reduced temperature using y -> (τ^(-1)).

U[τ] = U[y] /. yToTau;
S[τ] = S[y] /. yToTau;
F[τ] = F[y] /. yToTau;
(Subscript[C, V])[τ] = (Subscript[C, V])[y] /. yToTau;

plotUvsT = 
        Plot[Evaluate[U[τ]/(n ϵ)], {τ, 0.01, 5}, 
          PlotLabel -> "Energy", PlotRange -> {{0, 5}, {0, 4}}, 
          GridLines -> None, Frame -> True, 
          FrameLabel -> {"τ",  label[U]}, 
          DisplayFunction -> Identity];

plotCVvsT = 
        Plot[Evaluate[(Subscript[C, V])[τ]/(n (Subscript[k, B]))], {τ, 0.01, 5}, 
          PlotLabel -> "Heat Capacity", 
          PlotRange -> {{0, 5}, {0, 1}}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {"τ",  label[(Subscript[C, V])]}, 
          DisplayFunction -> Identity];

plotSvsT = 
        Plot[Evaluate[S[τ]/(n (Subscript[k, B]))], {τ, 0.01, 5}, 
          PlotLabel -> "Entropy", 
          PlotRange -> {{0, 5}, {0, 3}}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"τ",  label[S]}, 
          DisplayFunction -> Identity];

plotFvsT = 
        Plot[Evaluate[F[τ]/(n ϵ)], {τ, 0.01, 5}, 
          PlotLabel -> "Free Energy", 
          PlotRange -> {{0, 5}, {-4, 1}}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"τ",  label[F]}, 
          DisplayFunction -> Identity];

Show[
        GraphicsArray[{{plotUvsT, plotFvsT}, {plotCVvsT, plotSvsT}}, 
          GraphicsSpacing -> {-0.1, 0.05}, 
          DisplayFunction -> $DisplayFunction]];

Limiting Behaviors

Considerable insight into the behavior of physical systems is often obtained from analytic expressions for their limiting behaviors under extreme conditions.  In this section we examine the low and high temperature limits for the principle thermodynamic functions.

High Temperature

The high temperature limit can be obtained easily using Taylor series  expansion with respect to the inverse temperature, y, about zero.  We can easily transform between y and T, if desired.

Series[U[y], {y, 0, 2}] // 
      Collect[#1, {n, (Subscript[k, B]), ϵ}] &

Series[S[y], {y, 0, 2}]  // 
        Simplify//Collect[#1, {n, 
            (Subscript[k, B]), ϵ}] &

Series[F[y], {y, 0, 2}]  // Simplify // 
      Collect[#1, {n, (Subscript[k, B]), ϵ}] &

Series[(Subscript[C, V])[y], {y, 0, 2}] // 
      Collect[#1, {n, (Subscript[k, B]), ϵ}] &

Thus, we find that the heat capacity approaches the classical equipartition  limit, (Subscript[C, V])  ->  n (Subscript[k, B]), in the high temperature (small y) limit, with the correction term being of order (((((ℏ ω)/((Subscript[k, B]) T))))^2).

Low Temperature

The low temperature limit is more difficult because the Boltzmann factor, \(TraditionalForm\`Exp[\(±\(\(ℏ\ ω\)\/\(k\_B\ T\)\)\)]\), has an essential singularity for T -> 0.  Nevertheless, the limiting behavior can be obtained by using explicit  rules for isolating the singularity in the Boltzmann factor as manifested by various exponential and hyperbolic functions.  Since we expect all of the thermodynamic functions to be well behaved at low temperature, the strategy in formulating these rules is to factor the part that appears divergent, namely \(TraditionalForm\`Exp[\(+\(\(ℏ\ ω\)\/\(k\_B\ T\)\)\)]\), leaving behind an expansion in powers of the small quantity \(TraditionalForm\`Exp[\(-\(\(ℏ\ ω\)\/\(k\_B\ T\)\)\)]\).

lowTrules = {
	Cosh[a_. 1/τ]-> Exp[a/τ] (1 + Exp[-2  a/τ])/2,
	Sinh[a_. 1/τ]-> Exp[a/τ] (1 - Exp[-2 a/τ])/2,
 	Sech[a_. 1/τ]-> 2 Exp[-a/τ] (1 - Exp[-2 a/τ]),
   	Csch[a_. 1/τ]-> 2 Exp[-a/τ] (1 + Exp[-2 a/τ]),
   	Tanh[a_. 1/τ]-> 1 - 2 Exp[-2 a/τ],
  	Coth[a_. 1/τ]-> 1 + 2 Exp[-2 a/ τ],
  	Log[1-Exp[-1/τ]]-> -Exp[-1/τ]
};

U[τ]//.lowTrules

(F[τ] //.lowTrules//PowerExpand)//.lowTrules//Simplify

(S[τ] //. lowTrules // PowerExpand) //. lowTrules  // 
        Expand  // 
      Collect[#1, {n, (Subscript[k, B]), (E^(-1/τ))}] &

(Subscript[C, V])[τ] //. lowTrules  // Expand  // 
      Collect[#1, {n, (Subscript[k, B]), τ}] &

Comparison with classical statistical mechanics

The Einstein model of lattice vibrations in a crystal (the Einstein  crystal) approximates the spectrum of normal modes for the motion of N atoms by a single average frequency ω.  Since each atom can move in three dimensions, the number of normal  modes available to the crystal is n  =  3  N.  Each normal mode is quantized, such that the energy can be expressed as  U =  (m + (n/2)) ℏ ω where m is the number of quanta of excitation.  It is useful to introduce an Einstein temperature,  (Subscript[T, E])  =  ℏ ω / 
            (Subscript[k, B]) , such that y = (Subscript[T, E])/T is simply the inverse of the reduced temperature is T/(Subscript[T, E]).  Hence, this model contains one parameter, (Subscript[T, E]), characterizing the stiffness of the lattice, which may be fitted to  experimental observations.  Typically one finds that (Subscript[T, E]) is a few hundred kelvin.

Heat capacity

At high temperatures (y ≪ 1), we obtain the classical behavior predicted by the equipartition theorem  in which each of the n kinetic and n potential energies contribute (1/2) (Subscript[k, B]) T to the internal energy and (1/2) (Subscript[k, B]) to the heat capacity.  At low temperature (y ≫ 1), the heat capacity declines as \(TraditionalForm\`y\^2\ e\^\(-y\) = \((T\_E\/T)\)\^2\ Exp[\(-\(T\_E\/\(\(\ \)\(T\)\)\)\)]\) because it becomes very difficult to excite quantized oscillations when the thermal energy is small compared with the oscillator spacing.  This  suppression of the heat capacity at low temperature is characteristic of any system for which there is a gap in the excitation spectrum and hence is obtained whenever the energy levels are discrete (quantized).  The most  important difference between the classical and quantum models is that classically it is possible to excite a vibrational mode with arbitrarily  small amplitude and energy, whereas quantum mechanically that mode cannot be  excited with less than one full quantum.  Thus, when the thermal energy (Subscript[k, B]) T becomes comparable to or smaller than the basic quantum of excitation ℏ ω, that mode is frozen out and the classical equipartition theorem fails.   Therefore, even though the Einstein model is a very simplistic description of the thermodynamics of lattice vibrations in crystals, it nevertheless captures an essential aspect of the physics of quantized systems.

Show[Graphics[ plotCVvsT], PlotRange -> {{0, 2}, {0, 1}}, 
        DisplayFunction -> $DisplayFunction];

Note that when the reduced temperature reaches 1, the heat capacity is already close to its classical limit.  Comparing this heat capacity with experimental values, one finds that the characteristic temperature in the Einstein crystal, (Subscript[T, E])  =  ℏ ω / 
            (Subscript[k, B]) , is typically a few hundred degrees Kelvin, depending upon the stiffness of the crystal.  Thus, significant departures from classical behavior are observed at room temperature for many materials.  The Einstein model provides a good description of experimental data except at very low temperatures where the use of a single average frequency breaks down.  The success of this model in explaining such an elementary thermodynamic quantity for which classical reasoning failed contributed to the acceptance of the early quantum theory.

Zero-point energy

It is also interesting to compare the classical energy to quantum mechanical energy with and without the zero-point contribution.  In the classical model the internal energy is directly proportional to the temperature, whereas in the quantum mechanical model the internal energy approaches a constant at low temperature, namely the zero-point energy.   It is interesting that the zero-point energy is needed to obtain agreement with the classical limit at high temperatures.  Although such an offset has no dynamical consequences for this system, it is in principle observable via mass measurements (gravitational or inertial) and hence it is comforting to recover the classical limit when zero-point motion is included.

\(\(Plot[
        Evaluate[{U[y]\/\(n\ ϵ\), U[y]\/\(n\ ϵ\) - 0.5, 
              T} /. {y -> 1\/T}], {T, 0.01, 3}, 
        PlotStyle -> {{}, Dashing[{0.04, 0.02}], Dashing[{0.02, 0.02}]}, 
        PlotLabel -> "Energy", PlotRange -> {{0, 3}, {0, 3}}, 
        GridLines -> None, Frame -> True, 
        FrameLabel -> {"τ", label[U]}, 
        PlotLegend -> {"Quantum", "no zero-point", "Classical"}, LegendPosition -> {0.3, \(-0.3\)}, 
        LegendSize -> {1.0, 0.3}];\)\)

Canonical Ensemble

Finally, we demonstrate that the same thermodynamic functions can be obtained more efficiently using the canonical ensemble.  The single-particle partition function is denoted (Subscript[Z, 1]).

(Subscript[Z, 1])  =  (∑_(m = 0)^∞) 
              Exp[-(m + (1/2)) y]  // ExpToTrig  // 
        Simplify

fcanonical = -n (Subscript[k, B]) T Log[(Subscript[Z, 1])] //. TtoY  // 
        Simplify

scanonical =  -(∂_y) fcanonical dydT  /. TtoY  // 
        Simplify

ucanonical =  (∂_y)(y fcanonical)  // 
        Simplify

Since these results all agree with those obtained from the microcanonical ensemble, further analysis would be redundant.  However, for most systems we will find that the canonical ensemble provides simpler methods of calculation, either exact or approximate.