Rotations and Vibrations for Ideal Gases

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

Introduction

The thermodynamics of rotations and vibrations is investigated.  We consider heteronuclear diatomic molecules and omit the degeneracy of the nuclear wave function.  The rotational energy is expressed in terms of a temperature scale, (Subscript[T, rot])  =  (ℏ^2)/(2 (Subscript[I, rot]) (Subscript[k, B])) where (Subscript[I, rot]) is the moment of inertia.  Only a single vibrational mode is considered in the harmonic approximation.  The vibrational energy is also expressed in terms of a temperature scale, (Subscript[T, vib])  =  ℏ ω / (Subscript[k, B]).  Electronic excitations are omitted.  Some of these restrictions are generalized in the problems at the end.

Glossary:
	𝒩	= number of atoms
	(Subscript[k, B])	= Boltzmann constant
	j	= angular momentum
	(Subscript[I, rot])	= rotational inertia
	T	= temperature
	(Subscript[T, rot])	= rotational temperature	= (ℏ^2)/(2 (Subscript[I, rot]) (Subscript[k, B]))
	(Subscript[T, vib])	= vibrational temperature 	= ℏ ω / (Subscript[k, B])
	y	= inverse temperature (in rotational or vibrational units)
	(Subscript[Z, rot])	= single-particle partition function for rotational motion
	(Subscript[Z, vib])	= single-particle partition function for vibrational motion
	F	= Helmholtz free energy
	U	= internal energy 
	S	= entropy
	(Subscript[C, V])	= isochoric heat capacity

Initialization

Defaults,  packages, and symbols

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

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

Needs["Miscellaneous`PhysicalConstants`"];
Needs["Miscellaneous`Units`"];
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])_];Symbolize[(Subscript[F, x])_];Symbolize[(Subscript[S, x])_];Symbolize[(Subscript[U, x])_];Symbolize[(Subscript[CV, x])_];Symbolize[(Subscript[I, x])_]

Conversion rules

Tscales = {(Subscript[T, rot]) -> ((ℏ^2)/(2 (Subscript[I, rot]) (Subscript[k, B]))), 
          (Subscript[T, vib]) -> ((ℏ ω)/(Subscript[k, B])), ℏ -> 
            PlanckConstantReduced, 
          (Subscript[k, B]) -> BoltzmannConstant};

Error messages

Several annoying error messages can be suppressed, if desired, but should be enabled when developing or debugging code.  The first message occurs during the construction of the rotational partition function when the Boltzmann factor becomes incredibly small for very large angular momenta.  The second occurs when the plotting routines uses interpolation functions slightly outside their defined ranges.

Off[NSum::"nsnum"];
Off[InterpolatingFunction::"dmval"]

Rotational contribution

Partition function

Since the partition function cannot be summed analytically, we must perform numerical evaluations.  It is useful to divide the temperature range into low, medium, and high ranges with different step sizes.  The ranges are then joined and an interpolation function is constructed which can later be used to evaluate the thermodynamic functions.  The temperature and energy are expressed here in rotational units.

Note that care must be taken to avoid overlap in the temperature ranges because the interpolation function performs poorly when points are duplicated.

(Subscript[Z, rot])[y_] := 
      NSum[(2 j + 1) Exp[-j (j + 1) y], {j, 
          0, ∞}]

TableZ1RotLowT = 
        Table[{T, (Subscript[Z, rot])[y] /. y -> (1/T)}, {T, 0.05, 0.5, 
            0.05}];

TableZ1RotMidT = 
        Table[{T, (Subscript[Z, rot])[y] /. y -> (1/T)}, {T, 0.6, 2.9, 
            0.1}];

TableZ1RotHighT = 
        Table[{T, (Subscript[Z, rot])[y] /. y -> (1/T)}, {T, 3, 10, 
            0.5}];

TableZ1Rot = 
        Join[TableZ1RotLowT, TableZ1RotMidT, TableZ1RotHighT];

Z1RotInterpolatingFunction = 
        Interpolation[TableZ1Rot, InterpolationOrder -> 5];

\(\(DisplayTogether[
        ListPlot[TableZ1Rot, PlotStyle -> AbsolutePointSize[4]], 
        Plot[Z1RotInterpolatingFunction[T], {T, 0, 10}], Frame -> True, 
        FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""(Subscript[Z, 1])""}, 
        PlotLabel -> "Rotational Partition Function"];\)\)

Thermodynamic functions

(Subscript[F, rot])[
          T_] = -𝒩 (Subscript[k, B]) T Log[Z] /. {Z -> 
              Z1RotInterpolatingFunction[T]};
(Subscript[S, rot])[
          T_] = -(∂_T) (Subscript[F, rot])[T];
(Subscript[U, rot])[T_] = (Subscript[F, rot])[T] + T (Subscript[S, rot])[T]; 
    (Subscript[CV, rot])[T_] = (∂_T) (Subscript[U, rot])[T];

\(\(plotFrot\  = \ 
        Plot[Evaluate[F\_rot[T]\/\(𝒩\ k\_B\ T\)], {T, 0.05, 
            10}, PlotLabel -> "Rotational Free Energy", 
          PlotRange -> {{0, 10}, Automatic}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""F/(𝒩(Subscript[k, B])T)""}, 
          DisplayFunction -> Identity];\)\)

\(\(plotUrot\  = \ 
        Plot[Evaluate[U\_rot[T]\/\(𝒩\ k\_B\ T\)], {T, 0.05, 
            10}, PlotLabel -> "Rotational Energy", 
          PlotRange -> {{0, 10}, {0, 1}}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""U/(𝒩(Subscript[k, B])T)""}, 
          DisplayFunction -> Identity];\)\)

\(\(plotSrot\  = \ 
        Plot[Evaluate[S\_rot[T]\/\(𝒩\ k\_B\)], {T, 0.05, 10}, 
          PlotLabel -> "Rotational Entropy", 
          PlotRange -> {{0, 10}, {0, 4}}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""S/(𝒩(Subscript[k, B]))""}, 
          DisplayFunction -> Identity];\)\)

\(\(plotCVrot\  = \ 
        Plot[Evaluate[CV\_rot[T]\/\(𝒩\ k\_B\)], {T, 0.05, 
            10}, PlotLabel -> "Rotational Heat Capacity", 
          PlotRange -> {{0, 10}, {0, 1.2}}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""(Subscript[C, v])/(𝒩(Subscript[k, B]))""}, 
          DisplayFunction -> Identity];\)\)

Show[
        GraphicsArray[{{plotUrot, plotFrot}, {plotSrot, plotCVrot}}, 
          GraphicsSpacing -> {-0.1, 0.05}, 
          DisplayFunction -> $DisplayFunction ]];

Numerical estimates

To obtain an estimate of the characteristic rotational temperature for a diatomic molecule, we need to express the rotational inertia about the center-of-mass in terms of the nuclear masses and the distance between them.  The masses are obtained from the atomic weight and the interatomic separation is typically 1 - 2 angstroms,  but more precise values for the separation requires experimental data or a more detailed model.

(Subscript[r, cm])  = (( (Subscript[m, 1]) (Subscript[r, 1]) +  (Subscript[m, 2]) (Subscript[r, 2]))/((Subscript[m, 1]) + (Subscript[m, 2]))) /. {(Subscript[r, 2]) -> (Subscript[r, 1]) + d}  // Simplify

Irule  =  
      (Subscript[I, rot]) ->  (Subscript[m, 1]) ((((Subscript[r, 1]) - (Subscript[r, cm])))^2)  +  (Subscript[m, 2]) ((((Subscript[r, 2]) - (Subscript[r, cm])))^2)  /. {(Subscript[r, 2]) -> (Subscript[r, 1]) + d}  // Simplify

values[HCl]  =  {(Subscript[m, 1]) -> AtomicMassUnit, 
        (Subscript[m, 2]) -> 35 AtomicMassUnit, d -> 1.5 Angstrom}

(Subscript[I, rot])  /. Irule /. values[HCl]

((Subscript[T, rot])  //. Tscales  /. Irule  /. values[HCl]  // SI) /. 
      Joule -> Kilogram (Meter^2) (Second^(-2))

Notice that the rotational inertia for diatomic molecules is smallest when the atoms are light, such as (Subscript[H, 2]), and the corresponding rotational temperature is then relatively large.  For HCl, (Subscript[T, rot]) is only about 10 kelvin and will be smaller for most nonhydrogenic molecules.  Therefore, the rotational contribution to the heat capacity usually saturates at temperatures of only a few tens of kelvins, very much below room temperature.

Vibrational contribution

Thermodynamic functions

The temperature and energy are expressed here in vibrational units.  The zero-point energy does not contribute to the heat capacity and is omitted.

(Subscript[Z, vib]) = (∑_(n = 0)^∞) Exp[-n y]

(Subscript[F, vib])[
          T_] = -𝒩 (Subscript[k, B]) T Log[(Subscript[Z, vib])]  /. {y -> 
                1/T} // Simplify;
(Subscript[S, vib])[T_] = -(∂_T) (Subscript[F, vib])[T] // Simplify;
(Subscript[U, vib])[T_] = (Subscript[F, vib])[T] + T (Subscript[S, vib])[T] // Simplify;
(Subscript[CV, vib])[T_] = (∂_T) (Subscript[U, vib])[T] // Simplify;

{(Subscript[U, vib])[T], (Subscript[CV, vib])[T]}

\(\(plotFvib\  = \ 
        Plot[Evaluate[F\_vib[T]\/\(𝒩\ k\_B\ T\)], {T, 0.05, 
            10}, PlotLabel -> "Vibrational Free Energy", 
          PlotRange -> {{0, 10}, Automatic}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""F/(𝒩(Subscript[k, B])T)""}, 
          DisplayFunction -> Identity];\)\)

\(\(plotSvib\  = \ 
        Plot[Evaluate[S\_vib[T]\/\(𝒩\ k\_B\)], {T, 0.05, 10}, 
          PlotLabel -> "Vibrational Entropy", 
          PlotRange -> {{0, 10}, {0, 4}}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""S/(𝒩(Subscript[k, B]))""}, 
          DisplayFunction -> Identity];\)\)

\(\(plotUvib\  = \ 
        Plot[Evaluate[U\_vib[T]\/\(𝒩\ k\_B\ T\)], {T, 0.05, 
            10}, PlotLabel -> "Vibrational Energy", 
          PlotRange -> {{0, 10}, {0, 1}}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""U/(𝒩(Subscript[k, B])T)""}, 
          DisplayFunction -> Identity];\)\)

\(\(plotCVvib\  = 
        Plot[Evaluate[CV\_vib[T]\/\(𝒩\ k\_B\)], {T, 0.05, 
            10}, PlotLabel -> "Vibrational Heat Capacity", 
          PlotRange -> {{0, 10}, {0, 1.2}}, GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, rot])"", \*""(Subscript[C, V])/(𝒩(Subscript[k, B]))""}, 
          DisplayFunction -> Identity];\)\)

Show[
        GraphicsArray[{{plotUvib, plotFvib},  {plotSvib, plotCVvib}}, 
          GraphicsSpacing -> {-0.1, 0.05}, 
          DisplayFunction -> $DisplayFunction ]];

Numerical estimates

Reliable numerical estimates for the characteristic vibrational temperature are more difficult to make without recourse to experimental data for the interatomic potential.  Typically one finds that the vibrational quantum is on the order of 0.5 eV.

(Subscript[T, vib])  //. 
          Tscales  /. {ω ->  
            0.5 ElectronVolt/PlanckConstantReduced} // SI

Therefore, the vibrational temperature is typically a few thousand kelvins and the vibrational contribution to the heat capacity near room temperature is quite small.  However, for some diatomic molecules the vibrational energies can be as much as a factor of 10 smaller and hence can be relevant near room temperature.

Examine temperature dependence of heat capacity

First we construct a function which combines the translational, rotational, and vibrational contributions with proper attention to their scale parameters and the interpolation range for the rotational contribution.

(Subscript[C, V])[T_, (Subscript[T, rot]) : _, (Subscript[T, vib]) : _] := 
      Module[{(Subscript[τ, rot]) = (T/(Subscript[T, rot])), (Subscript[τ, vib]) = (T/(Subscript[T, vib])), 
          rotCV, τ, x}, rotCV[τ_ /; τ < 0.05] := 0; 
        rotCV[τ_ /; τ > 0 && τ < 10] := 
          (Subscript[CV, rot])[x] /. x -> τ; 
        rotCV[τ_ /; τ > 10] := 𝒩 (Subscript[k, B]); 
        N[((3 𝒩 (Subscript[k, B]))/2) + rotCV[(Subscript[τ, rot])] + 
              (Subscript[CV, vib])[x] /. x -> (Subscript[τ, vib])]]

The following parameters for the HCl molecule are take from the literature (e.g. Rosser, Sec. 15.10).

values = {(Subscript[T, rot]) -> 15, (Subscript[T, vib]) -> 4130};

A log-linear plot can be used to display several orders of magnitude variation in temperature, but because Mathematica does not choose its points wisely for logarithmic plots, we first construct a list of values uniformly distributed in Log[T].

\(\(LogLinearListPlot[
        Table[{10\^x, 
            C\_V[10\^x, T\_rot, T\_vib]\/\(𝒩\ k\_B\) /. 
              values}, {x, 0, 4, 0.05}], 
        PlotRange -> {Automatic, {0, 4}}, PlotJoined -> True, 
        PlotLabel -> "Heat Capacity", GridLines -> None, 
        Frame -> True, 
        FrameLabel -> {"T", \*""(Subscript[C, V])/(𝒩(Subscript[k, B]))""}, 
        DisplayFunction -> $DisplayFunction];\)\)

For very low temperatures, well below (Subscript[T, rot]), the heat capacity displays the exponential behavior characteristic of systems with an energy gap.  Near (Subscript[T, rot]), the rotational  contribution overshoots the equipartition limit.  Near room temperature, the molar heat capacity for this gas is well approximated by the equipartition value of (5/2) R obtained by assuming that two rotational degrees of freedom are active while the vibrational mode is frozen out.  Only for temperatures near (Subscript[T, vib]), here about 5000 Kelvin, is the full classical vibrational contribution obtained.

Problems

(Subscript[Z, r]^FD)  =  (Subscript[s, A]) (2  (Subscript[s, A]) + 1) 
          (Subscript[Z, even])  +  ((Subscript[s, A]) + 1) (2  (Subscript[s, A]) + 1) 
          (Subscript[Z, odd])

(Subscript[Z, even]) = (∑_(ℓ even))(2  ℓ + 1) 
          exp[-ℓ (ℓ + 1) 
              ((Subscript[T, rot])/T)]

(Subscript[Z, odd])  =  (∑_(ℓ odd)) (2  ℓ + 1) 
          exp[-ℓ (ℓ + 1) 
              ((Subscript[T, rot])/T)]

(Subscript[C, r]^FD)  =  ((Subscript[s, A])/(2  (Subscript[s, A]) + 1)) 
          (Subscript[C, even])  +    (((Subscript[s, A]) + 1)/(2  (Subscript[s, A]) + 1)) 
          (Subscript[C, odd])