Statistical Mechanics of Ideal Two-Level Systems

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

Introduction

An ideal two-level system consists of N noninteracting elements which each can occupy only two states with energies (Subscript[ε, 0]) or (Subscript[ε, 0]) + Δ.  For example, a system of N noninteracting spin-((1 )/2) atoms or nuclei subject to an external magnetic field can be treated as an ideal two-level system if mutual interactions and other degrees of freedom, such as kinetic energies, can be neglected.  Alternatively, a semiconductor is a system in which at T = 0 all single-particle states in the valence band would be filled while all states in the conduction band are empty.  Although a rigorous treatment requires the Fermi-Dirac distribution, the qualitative behavior of such a system is governed by the energy gap Δ between the top of the valence band and the bottom of the conduction  band.  Since the Pauli exclusion principle prohibits two electrons from occupying the same single-particle state, a semiconductor is an electrical insulator at T = 0 because it is not possible for electrons in the conduction band to absorb an infinitesimal momentum without also absorbing the finite energy required to promote an electron across the energy gap.  At higher temperature, thermal  fluctuations lead to a finite population of the conduction band and hence to a finite conductivity which depends strongly upon temperature according to the Boltzmann factor (e^(-βΔ)).  This factor also affects the thermodynamic properties of the electronic subsystem.  Thus, the characteristic peak in the temperature dependence of the heat capacity for a two-level system is often called the Schottky anomaly after one of the inventors of the transistor, a revolutionary  semiconductor device.  This peak is not present in classical statistical mechanics and is a manifestation of the quantization of the excitation-energy  spectrum.  At low temperatures, the thermodynamics of many systems with a gap in the excitation spectrum can be understood using a two-level model.

In this notebook we derive the thermodynamic properties of ideal two-level systems using the canonical ensemble.  The free energy is

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

where

\(TextForm\`Z\_1 = \(∑\+i 
              Exp[\(-β\) ε\_i]\ \ \ \ \ \ \ \ \ \ \ β = 1\/\(\(k\_B\) T\)\)\)

is the single-particle partition function summed over the energy levels (Subscript[ε, i]) of a single particle.  The energy levels for a binary system can be labeled (Subscript[ε, 0]) and (Subscript[ε, 0]) + Δ, so that the partition function becomes

\(TextForm\`Z\_1 = \(\(Exp[\(-β\) ε\_0]\) \((1 + 
            Exp[\(-βΔ\)])\)\(\ \)\)\)

The standard thermodynamic functions are then

\(TextForm\`		S = \(-∂\_T F\)\)
\(TextForm\`		U = F\  + \ T S\)
\(TextForm\`		C\_V = ∂\_T U\)

Although the algebra is fairly simple, we employ Mathematica to derive expressions for these functions, produce plots, and evaluate limiting cases.  To express the thermodynamic functions in terms of different choices of independent variables with minimal symbolic changes, we employ different arguments.  For example, S[T] and S[y] express the entropy in terms of either the temperature, T, or the ratio between the gap and the thermal energies, \(TraditionalForm\`\(\(y\  = \ Δ/\ \((k\_B\ \ T)\)\)\(\ \)\)\), respectively.  Similarly, we use S[τ] to express entropy in terms of reduced temperature, where τ=1/y.

Glossary:
	𝒩	= number of particles
	ε	= energy of single-particle ground-state 
	ε + Δ	= energy of single-particle excited state
	Δ	= energy gap
	(Subscript[k, B])	= Boltzmann constant
	T	= temperature
	β	= 1/(Subscript[k, B]) T 
	y	=  Δ/ ((Subscript[k, B])  T) = gap/thermal energies
	τ	= 1/y = reduced temperature
	(Subscript[Z, 1])	= single-particle partition function
	(Subscript[Z, N])	= N-particle partition function
	S	= entropy
	U	= internal energy 
	F	= Helmholtz free energy

Initialization

Defaults, packages, and symbols

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

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

Needs["Utilities`Notation`"]

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

Rules for changing variables

yToT = {y -> (Δ/((Subscript[k, B]) T)),  β -> 
            (1/((Subscript[k, B]) T))};
TtoTau  =  {T -> ((Δ τ)/((Subscript[k, B]) ))};
TauToY  =  {τ -> (1/y)};

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\)""];\)\)

Canonical Ensemble

Partition function

The single-particle partition function is obtained by adding the Boltzmann factors for the two single-particle states..

(Subscript[Z, 1]) = 
      Exp[-β ε] + 
          Exp[-β (ε + Δ)] /. 
        yToT

(Subscript[Z, N])  =  ((Subscript[Z, 1])^𝒩)

Principal thermodynamic functions

F[
        T] = ((- (Subscript[k, B]) T Log[(Subscript[Z, N])]  // PowerExpand) /. 
              Exp[a_ (b_ + c_)] -> Exp[a b] Exp[a c]  // 
            Simplify) // PowerExpand

S[T] = -(∂_(T ))F[
              T] //FullSimplify

U[T]  =  ((F[T]  +  T S[T] ) // Simplify )  // 
          Collect[#,  {ε}] & // Cancel

(Subscript[C, V])[T]  =  (∂_T) U[T]  // Simplify

Examine thermodynamic functions

It is useful to express the temperature in dimensionless form using \(TraditionalForm\`y = \(τ\^\(-1\) = Δ\/\(k\_B\ T\)\)\).  In this section we plot the principal thermodynamic functions in terms of reduced temperature, τ, or inverse temperature, y.

Express thermodynamic functions in terms of τ or y

U[τ]  =  U[T]  /. TtoTau  // Simplify;
U[y]  =  U[τ]  /. TauToY  // Simplify

S[τ]  = S[T]  /. TtoTau  // Simplify;
S[y]  =  S[τ]  /. TauToY  // Simplify

F[τ]  =  F[T]  /. TtoTau  // Simplify;
F[y]  =  F[τ]  /. TauToY  // Simplify

(Subscript[C, V])[τ]  =  (Subscript[C, V])[T]  /. TtoTau  // Simplify;
(Subscript[C, V])[y]  =  (Subscript[C, V])[τ]  //. TauToY

The ground-state energy ε simply represents an offset of the energy scale and appears nowhere else in the model — it has no thermodynamical significance because it is not coupled to the temperature and hence it may be taken to zero in the plots below.

Plots versus temperature

Here the thermodynamic functions are expressed in terms of reduced  temperature using T -> (τ Δ/(Subscript[k, B])).  We set the single-particle ground-state energy to 0, for convenience.

\(\(plotUvsT = 
        Plot[Evaluate[{U[τ]\/\(𝒩\ Δ\) /. ε -> 0}], {τ, 0, 5}, 
          PlotLabel -> "Internal Energy", 
          PlotRange -> {{0, 5}, {0, 0.5}}, 
          PlotStyle -> {{}, Dashing[{0.02, 0.02}], 
              Dashing[{0.04, 0.01, 0.02, 0.01}]}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"τ", label[U]}, 
          DisplayFunction -> Identity];\)\)

\(\(plotFvsT = 
        Plot[Evaluate[{F[τ]\/\(𝒩\ Δ\) /. ε -> 0}], {τ, 0, 5}, 
          PlotLabel -> "Free Energy", 
          PlotRange -> {{0, 5}, {\(-4\), 1}}, 
          PlotStyle -> {{}, Dashing[{0.02, 0.02}], 
              Dashing[{0.04, 0.01, 0.02, 0.01}]}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"τ", \ label[F]}, 
          DisplayFunction -> Identity];\)\)

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

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

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

Plots versus inverse temperature

\(\(plotUvsY = 
        Plot[Evaluate[
            U[y]\/\(𝒩\ Δ\) /. ε -> 0], {y, 0, 5}, PlotLabel -> "Internal Energy", 
          PlotRange -> {{0, 5}, {0, 0.5}}, GridLines -> None, 
          Frame -> True, FrameLabel -> {"y", \ label[U]}, 
          DisplayFunction -> Identity];\)\)

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

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

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

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

At high temperature (small y), the system approaches its maximum possible entropy, N (Subscript[k, B]) Log[2], as each of the (2^N) microstates is equally likely when (Subscript[k, B]) T  ≫  Δ.  The internal energy is then simply the average of the minimum and maximum possible energies, namely \(TraditionalForm\`N\ \((ε\  + \ Δ\/2)\)\).  The heat capacity, (Subscript[C, V])  =  T(∂(S/∂)T), approaches zero because as the entropy asymptotically approaches its maximum value its derivative approaches zero; the system can no longer increase its entropy by absorbing heat.  This behavior is a special consequence of the assumption that the system has a finite number of microstates and, hence, a maximum possible entropy.  However, real physical systems do not satisfy this assumption because there is always an energy scale and corresponding temperature at which it becomes possible to excite internal degrees of freedom arising from the internal structure of its constituents.  Nevertheless, it is often useful to consider a limited set of excitation modes, such as the magnetic alignment energy of the atoms within a lattice, as a distinct subsystem which may exhibit this behavior.  For systems such as paramagnetic salts, the energy of the lattice excitations is not bounded; although the spin contribution to the heat capacity vanishes asymptotically at high temperature the total heat capacity does not.

At low temperature (large y) the internal energy and free energy approach N ε, which is the lowest possible energy obtained when all particles occupy the single-particle ground-state.  Similarly, at low temperature the entropy approaches zero because a system in which all particles are found in the same state is obviously highly ordered.  The heat capacity also approaches zero, which is a nonclassical effect discussed in the next  section.

Schottky anomaly

The heat capacity decreases at high temperature (small y) because as the system approaches its maximum possible entropy it can no  longer increase its entropy by absorbing heat.  The behavior is characteristic of any system for which the energy and entropy are bounded, as discussed for spin-(1/2) paramagnetism in the notebook spin-half.nb, for example.  At low temperature (large y) the heat capacity again decreases because it is impossible for the system to absorb a quantity of heat that is smaller than the single-particle  energy gap.  When (Subscript[k, B]) T  ≲  Δ is small compared to the energy needed to promote a particle to the upper  energy level, most particles are found in the single-particle ground state and the system has very low entropy because it is highly ordered.  Because the system cannot absorb arbitrarily small quantities of heat without a finite temperature increase, its heat capacity approaches zero as τ → 0.  This behavior is decidedly nonclassical — it is an essential feature of the quantization of the excitation energy spectrum.  Classical systems can receive arbitrarily small excitation energies, but if a quantized system is not able to accept a full quantum of energy then it cannot absorb any energy at all.  Therefore, the heat capacity approaches zero for both extremes of temperature and displays a characteristic maximum for some intermediate temperature.  This phenomenon is named the Schottky effect.  The location of the maximum can be determined by differentiation, but the resulting transcendental equation must be solved numerically.

dCvDy = (∂_y) (Subscript[C, V])[y]/(𝒩 (Subscript[k, B]))  // 
        Simplify

eq2 =  -dCvDy (((1 + (E^y)))^3) / ((E^y) y) == 0  // 
        Simplify

ycvMax = FindRoot[eq2, {y, 2}]

(Subscript[C, V])[y] /. ycvMax

Limiting Behaviors

Considerable insight into the behavior of various systems is often obtained from symbolic expressions for their limiting behavior.  In this section we examine the low and high temperature limits.

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[#, {(Subscript[k, B]),  𝒩,  Δ ,  ε,  y}] &

Series[S[y], {y, 0, 2}]  // Simplify // 
      Collect[#, {(Subscript[k, B]),  𝒩,  Δ ,  ε,  τ}] &

Series[F[y], {y, 0, 2}]  // 
      Collect[#, {(Subscript[k, B]),  𝒩,  Δ ,  ε,  τ}] &

Series[(Subscript[C, V])[y], {y, 0, 2}]  // 
      Collect[#, {(Subscript[k, B]),  𝒩,  Δ ,  ε,  τ}] &

Thus, at high temperature the entropy is near its maximum possible value,  (Subscript[S, max]) = N (Subscript[k, B]) Log[2]. The heat capacity is small, (Subscript[C, V]) ∝ (T^(-2)), because the state of maximum disorder cannot increase its entropy  further by absorbing heat.

Low Temperature

The essential singularity in Exp[1/τ] at τ = 0 is handled by representing the exponentials in terms of a new variable, d = Exp[-1/τ], which is small at low temperatures and hence can be used for Taylor  expansions for positive temperature.

lowTrules = {Exp[-(1/τ)] -> d, 
          Exp[(1/τ)] -> (1/d), Exp[n_(Integer/T)] -> (d^(-n)), 
          Log[d] -> -(1/τ)};

dToTau  =  {d -> Exp[-(1/τ)]};

Series[(U[τ] /. lowTrules), {d, 0, 2}]  /. dToTau // 
      Collect[#, {(Subscript[k, B]),  𝒩,  Δ ,  ε,  τ}] &

Series[(F[τ] /. lowTrules), {d, 0, 2}]  /. dToTau // 
      Collect[#, {(Subscript[k, B]),  𝒩,  Δ,  ε,  τ}] &

\(\(\(Series[\((S[τ] /. lowTrules)\), {d, 0, 2}]\  /. dToTau\  // 
        Simplify\)\(//\)\(Collect[#, {k\_B, 𝒩, \ Δ, \ ε, \ τ}] &\)\(
\)\)\)

Series[((Subscript[C, V])[τ] /. lowTrules), {d, 0, 2}]  /. dToTau  // 
        Simplify // 
      Collect[#, {(Subscript[k, B]),  𝒩,  Δ,  ε,  τ}] &

When the thermal energy is small compared with the energy gap, most particles are  found in their lowest energy state and it becomes very difficult to populate the upper level.  Thus, at low temperature we find that the entropy of the system becomes small because the system exhibits a high  degree of order.  Similarly, the heat capacity becomes quite small because it is difficult for the system to absorb a quantity of heat on the order of (Subscript[k, B]) T when (Subscript[k, B]) T  ≲  Δ is small compared to the energy needed to promote a particle to the upper  energy level.

Special cases

Spin-(1/2) paramagnetism

The paramagnetic behavior of spin-(1/2) systems was studied using the microcanonical ensemble in the notebook spin-half.nb.  For the sake of completeness, and/or compulsiveness, we demonstrate that the present results reduce to those derived there when the variables are changed appropriately and the exponential functions are transformed to hyperbolic trigonometric functions.

ChangeVariablesToParamagnetism  = {ε -> -ε, Δ -> 2  ε, 
          y -> 2  y};

U[y] /. ChangeVariablesToParamagnetism  // ExpToTrig // 
      Simplify

(Subscript[C, V])[y] /. ChangeVariablesToParamagnetism  // ExpToTrig // 
      Simplify

S[y] /. ChangeVariablesToParamagnetism  /. 
          Log[1 + (E^(-2 y))] -> -y  +  Log[2 Cosh[y] ] // 
        ExpToTrig // Simplify

F[y] /. ChangeVariablesToParamagnetism  /. 
          Log[1 + (E^(-2 y))] -> -y  +  Log[2 Cosh[y] ] // 
        ExpToTrig // Simplify

As expected, these formulas reproduce those found in spin-half.nb.