Part IV: Statistical Thermodynamics | Ensemble Theory

Ensemble Theory

The systematic framework connecting microscopic dynamics to macroscopic thermodynamics through statistical ensembles — from Boltzmann's entropy to fluctuation theory and the grand canonical formalism

1. Introduction: Statistical Ensembles

At the heart of statistical mechanics lies a deceptively simple idea: rather than tracking the trajectory of every particle in a macroscopic system, we consider a vast collection of hypothetical copies of the system — an ensemble — each consistent with the known macroscopic constraints. Observable thermodynamic quantities are then identified with ensemble averages.

This ensemble approach, introduced by J. Willard Gibbs in his foundational 1902 treatise Elementary Principles in Statistical Mechanics, transforms the intractable problem of solving equations of motion for $10^{23}$ particles into a tractable exercise in probability theory. The central question becomes: what is the probability $P_i$ that the system occupies microstate $i$?

The Three Fundamental Ensembles

Different macroscopic constraints give rise to different ensembles, each suited to particular physical situations:

Microcanonical Ensemble (NVE): Isolated system with fixed energy $E$, volume $V$, and particle number $N$.
Canonical Ensemble (NVT): System in thermal contact with a heat bath at temperature $T$, with fixed $V$ and $N$.
Grand Canonical Ensemble ($\mu$VT): System exchanging both energy and particles with a reservoir at temperature $T$ and chemical potential $\mu$.

In the thermodynamic limit ($N \to \infty$, $V \to \infty$ with $N/V$ fixed), all three ensembles yield identical predictions for intensive thermodynamic quantities. This equivalence is one of the deepest results in statistical mechanics, rooted in the sharpness of probability distributions for macroscopic systems.

2. Derivation 1: The Microcanonical Ensemble

Consider an isolated system with fixed energy $E$, volume $V$, and particle number $N$. The fundamental question is: how do we assign probabilities to the accessible microstates?

2.1 The Equal A Priori Probability Postulate

The foundational postulate of statistical mechanics states that for an isolated system in equilibrium, all accessible microstates are equally probable:

Equal A Priori Probability Postulate

$$P_i = \begin{cases} \frac{1}{\Omega(E,V,N)} & \text{if } E_i = E \\ 0 & \text{otherwise} \end{cases}$$

where $\Omega(E,V,N)$ is the number of microstates with energy $E$.

This postulate cannot be derived from mechanics alone — it is the bridge between dynamics and thermodynamics. Its justification rests on ergodic theory: for a sufficiently chaotic system, the time average of any observable equals the ensemble average over the microcanonical distribution.

2.2 Boltzmann Entropy: S = k ln $\Omega$

We now derive the connection between the microcanonical ensemble and thermodynamics. The entropy of an isolated system must satisfy three requirements: (i) it is extensive, (ii) it is maximized at equilibrium, and (iii) it is additive for independent subsystems. Consider two independent subsystems $A$ and $B$:

$$\Omega_{\text{total}} = \Omega_A \cdot \Omega_B$$

The total number of microstates is multiplicative for independent systems.

For entropy to be additive ($S_{\text{total}} = S_A + S_B$) while $\Omega$ is multiplicative, $S$ must be a logarithmic function of $\Omega$:

Boltzmann Entropy

$$S = k_B \ln \Omega(E, V, N)$$

The proportionality constant $k_B = 1.381 \times 10^{-23}$ J/K is Boltzmann's constant, chosen to match the thermodynamic temperature scale.

2.3 Temperature from the Microcanonical Ensemble

We derive the thermodynamic temperature by considering two subsystems $A$ and $B$ in thermal contact, sharing total energy $E = E_A + E_B$. The combined density of states is:

$$\Omega_{\text{total}}(E) = \sum_{E_A} \Omega_A(E_A)\,\Omega_B(E - E_A)$$

Equilibrium corresponds to the energy partition $E_A^*$ that maximizes $\Omega_{\text{total}}$, equivalently maximizing $\ln \Omega_A + \ln \Omega_B$. Setting the derivative to zero:

$$\frac{\partial \ln \Omega_A}{\partial E_A} = \frac{\partial \ln \Omega_B}{\partial E_B}$$

This equality of derivatives at equilibrium is precisely the condition of equal temperature. We identify:

Microcanonical Temperature

$$\frac{1}{T} = \left(\frac{\partial S}{\partial E}\right)_{V,N} = k_B \frac{\partial \ln \Omega}{\partial E}$$

Similarly, pressure and chemical potential follow from $P/T = (\partial S/\partial V)_{E,N}$ and $\mu/T = -(\partial S/\partial N)_{E,V}$.

2.4 Application: Ideal Gas in the Microcanonical Ensemble

For $N$ non-interacting classical particles in a box of volume $V$ with total energy $E$, we must count the volume of the $3N$-dimensional momentum-space hypersphere $\sum_{i=1}^{3N} p_i^2 = 2mE$. The result, using the surface area of a$d$-dimensional sphere of radius $R$ is $S_d(R) = 2\pi^{d/2} R^{d-1}/\Gamma(d/2)$:

$$\Omega(E,V,N) = \frac{V^N}{N!\,h^{3N}} \cdot \frac{2\pi^{3N/2}}{(3N/2 - 1)!\,} (2mE)^{(3N-1)/2}\,\delta E$$

The $N!$ accounts for indistinguishability; $h^{3N}$ discretizes phase space; $\delta E$ is the energy shell width.

Taking the logarithm using Stirling's approximation ($\ln N! \approx N \ln N - N$) for large $N$:

$$S = Nk_B \left[\ln\!\left(\frac{V}{N}\right) + \frac{3}{2}\ln\!\left(\frac{4\pi m E}{3Nh^2}\right) + \frac{5}{2}\right]$$

This is the Sackur-Tetrode equation for the entropy of an ideal gas.

From this entropy, all ideal gas thermodynamics follows. Applying the temperature relation:

$$\frac{1}{T} = \frac{\partial S}{\partial E} = \frac{3Nk_B}{2E} \quad \Longrightarrow \quad E = \frac{3}{2}Nk_BT$$

And applying the pressure relation $P/T = (\partial S/\partial V)_{E,N} = Nk_B/V$:

Ideal Gas Law from Statistical Mechanics

$$PV = Nk_BT$$

The equation of state emerges purely from counting microstates — no kinetic theory assumptions needed.

3. Derivation 2: The Canonical Ensemble

The canonical ensemble describes a system in thermal equilibrium with a heat bath at temperature $T$. Unlike the microcanonical ensemble, the system's energy is not fixed but fluctuates around a mean value $\langle E \rangle$.

3.1 Derivation via Maximum Entropy

We derive the canonical distribution by maximizing the Gibbs entropy subject to the constraint of fixed average energy. The Gibbs entropy for a discrete probability distribution $\{P_i\}$ is:

$$S = -k_B \sum_i P_i \ln P_i$$

We maximize $S$ subject to two constraints: (1) normalization $\sum_i P_i = 1$and (2) fixed mean energy $\sum_i P_i E_i = \langle E \rangle$. Using Lagrange multipliers$\alpha$ and $\beta$:

$$\frac{\partial}{\partial P_i}\left[-k_B \sum_j P_j \ln P_j - \alpha\!\left(\sum_j P_j - 1\right) - \beta\!\left(\sum_j P_j E_j - \langle E\rangle\right)\right] = 0$$

Carrying out the differentiation for each $P_i$:

$$-k_B(\ln P_i + 1) - \alpha - \beta E_i = 0 \quad \Longrightarrow \quad P_i = \exp\!\left(-\frac{\alpha + k_B}{k_B}\right)\exp\!\left(-\frac{\beta E_i}{k_B}\right)$$

Imposing normalization and identifying $\beta = 1/k_BT$ by comparison with thermodynamics:

Canonical (Boltzmann) Distribution

$$P_i = \frac{e^{-E_i / k_BT}}{Z} \qquad \text{where} \qquad Z = \sum_i e^{-E_i / k_BT}$$

$Z$ is the canonical partition function — the single most important quantity in statistical mechanics.

3.2 Thermodynamic Quantities from the Partition Function

The partition function $Z$ encodes all equilibrium thermodynamic information. Introducing$\beta = 1/k_BT$, we derive each quantity systematically:

Helmholtz Free Energy:

$$A = -k_BT \ln Z$$

Proof: From the definition $A = E - TS$ and $S = -k_B\sum_i P_i \ln P_i$, substitute $P_i = e^{-\beta E_i}/Z$:

$$S = -k_B \sum_i P_i(-\beta E_i - \ln Z) = \frac{\langle E\rangle}{T} + k_B \ln Z$$

$$\Longrightarrow \quad A = \langle E\rangle - TS = -k_BT \ln Z \qquad \square$$

Average Energy:

$$\langle E \rangle = -\frac{\partial \ln Z}{\partial \beta}$$

Proof: $\langle E\rangle = \sum_i E_i P_i = \frac{1}{Z}\sum_i E_i e^{-\beta E_i} = -\frac{1}{Z}\frac{\partial Z}{\partial \beta} = -\frac{\partial \ln Z}{\partial \beta}$.

Entropy:

$$S = k_B\!\left(\ln Z + \beta \langle E \rangle\right) = k_B \ln Z + \frac{\langle E \rangle}{T}$$

3.3 Equivalence of Ensembles

A crucial question is whether the canonical and microcanonical ensembles give the same thermodynamic predictions. The partition function can be written as:

$$Z = \sum_i e^{-\beta E_i} = \int dE\,\Omega(E)\,e^{-\beta E}$$

The integrand $f(E) = \Omega(E)e^{-\beta E} = e^{S(E)/k_B - \beta E}$ has an extremely sharp maximum at the energy $E^*$ satisfying $\partial S/\partial E = 1/T = k_B\beta$. For a macroscopic system, the width of this peak scales as $\sim\sqrt{N}$ while $E^*$ scales as $N$, so the relative width vanishes as $1/\sqrt{N}$:

Ensemble Equivalence in the Thermodynamic Limit

$$\ln Z \approx \ln \Omega(E^*) - \beta E^* = \frac{S(E^*)}{k_B} - \beta E^* = -\beta A$$

In the limit $N \to \infty$, the canonical free energy equals the microcanonical free energy exactly.

4. Derivation 3: The Grand Canonical Ensemble

When a system can exchange both energy and particles with a reservoir, we need the grand canonical ensemble. The reservoir is characterized by temperature $T$ and chemical potential $\mu$, which controls the average particle number.

4.1 The Grand Partition Function

Following the same maximum-entropy logic as Section 3, but now with an additional constraint on the average particle number $\langle N \rangle$, we introduce a third Lagrange multiplier. The resulting probability for microstate $i$ with energy $E_i$ and particle number $N_i$ is:

$$P_i = \frac{1}{\Xi}\,e^{-\beta(E_i - \mu N_i)}$$

The grand partition function $\Xi$ (also denoted $\mathcal{Z}$) sums over all possible particle numbers and all microstates for each $N$:

Grand Partition Function

$$\Xi = \sum_{N=0}^{\infty} e^{\beta\mu N}\,Z_N = \sum_{N=0}^{\infty} e^{\beta\mu N} \sum_{\{i\}_N} e^{-\beta E_i^{(N)}}$$

where $Z_N$ is the canonical partition function for exactly $N$ particles. The fugacity $z = e^{\beta\mu}$ controls particle number fluctuations.

4.2 Thermodynamic Quantities

The grand potential $\Phi_G = A - \mu N = -PV$ connects directly to $\Xi$:

Grand Potential and Equation of State:

$$PV = k_BT \ln \Xi$$

This remarkable relation gives the equation of state directly from the grand partition function.

The average particle number follows from differentiating with respect to $\mu$:

$$\langle N \rangle = k_BT \frac{\partial \ln \Xi}{\partial \mu} = \frac{1}{\beta}\frac{\partial \ln \Xi}{\partial \mu}$$

4.3 Application: Ideal Gas Chemical Potential

For a single-species ideal gas, the canonical partition function is $Z_N = z_1^N / N!$ where $z_1 = V/\Lambda^3$ is the single-particle partition function and $\Lambda = h/\sqrt{2\pi mk_BT}$ is the thermal de Broglie wavelength. The grand partition function becomes:

$$\Xi = \sum_{N=0}^{\infty} \frac{(e^{\beta\mu}\,z_1)^N}{N!} = \exp\!\left(e^{\beta\mu}\,\frac{V}{\Lambda^3}\right)$$

From $\langle N \rangle = e^{\beta\mu}\,V/\Lambda^3$, we immediately solve for the chemical potential:

Ideal Gas Chemical Potential

$$\mu = k_BT \ln\!\left(\frac{N\Lambda^3}{V}\right)$$

Since $N\Lambda^3/V \ll 1$ for a classical gas, $\mu < 0$ — adding a particle to an ideal gas always lowers the free energy.

We can also verify the equation of state: $PV = k_BT \ln \Xi = k_BT \cdot e^{\beta\mu}V/\Lambda^3 = Nk_BT$, recovering the ideal gas law as expected.

5. Derivation 4: Fluctuations in the Canonical Ensemble

One of the most powerful aspects of the canonical ensemble is its ability to predict not just average quantities, but fluctuations around those averages. These fluctuation relations provide a deep connection between microscopic dynamics and macroscopic response functions.

5.1 Energy Fluctuations

The variance of the energy in the canonical ensemble is computed from the partition function. Starting from $\langle E \rangle = -\partial \ln Z / \partial \beta$:

$$\langle E^2 \rangle = \frac{1}{Z}\sum_i E_i^2\,e^{-\beta E_i} = \frac{1}{Z}\frac{\partial^2 Z}{\partial \beta^2}$$

$$\langle E^2 \rangle - \langle E \rangle^2 = \frac{\partial^2 \ln Z}{\partial \beta^2} = -\frac{\partial \langle E \rangle}{\partial \beta}$$

Since $\beta = 1/k_BT$, we have $\partial/\partial\beta = -k_BT^2\,\partial/\partial T$, and recalling $C_V = \partial\langle E\rangle/\partial T$:

Energy Fluctuation-Response Relation

$$\langle (\Delta E)^2 \rangle = k_BT^2\,C_V$$

Energy fluctuations are directly proportional to the heat capacity — a fluctuation-dissipation relation.

5.2 Relative Fluctuations and the Thermodynamic Limit

For an ideal gas, $C_V = \frac{3}{2}Nk_B$ and $\langle E\rangle = \frac{3}{2}Nk_BT$, so:

$$\frac{\sqrt{\langle(\Delta E)^2\rangle}}{\langle E\rangle} = \frac{\sqrt{k_BT^2 \cdot \frac{3}{2}Nk_B}}{\frac{3}{2}Nk_BT} = \sqrt{\frac{2}{3N}}$$

Relative Energy Fluctuation

$$\frac{\Delta E_{\text{rms}}}{\langle E \rangle} \sim \frac{1}{\sqrt{N}} \;\xrightarrow{N\to\infty}\; 0$$

For $N \sim 10^{23}$, the relative fluctuation is $\sim 10^{-12}$ — energy is effectively sharp, and the canonical ensemble becomes indistinguishable from the microcanonical.

This result is general: for any extensive quantity, relative fluctuations scale as $1/\sqrt{N}$. This is precisely why macroscopic thermodynamics works — fluctuations are unobservably small for bulk matter, and different ensembles yield identical predictions. The exception occurs at phase transitions, where $C_V$ diverges and fluctuations become macroscopic.

6. Applications of Ensemble Theory

Modern Computational Applications

Monte Carlo Simulations: The Metropolis algorithm generates configurations distributed according to the canonical (Boltzmann) distribution. At each step, a trial move is accepted with probability $\min(1, e^{-\beta\Delta E})$, where $\Delta E$ is the energy change. This enables efficient sampling of high-dimensional configuration spaces without computing $Z$ explicitly.
Molecular Dynamics: By coupling the equations of motion to a thermostat (Nose-Hoover, Langevin), molecular dynamics trajectories sample the canonical ensemble. The ergodic hypothesis ensures that time averages equal ensemble averages, allowing computation of thermodynamic properties from dynamical simulations.
Free Energy Calculations: Techniques like thermodynamic integration, free energy perturbation, and the Wang-Landau algorithm exploit ensemble theory to compute free energy differences — quantities inaccessible to simple averaging. These are central to drug design and materials science.
Phase Transitions: The grand canonical ensemble naturally describes first-order phase transitions where particle number fluctuations diverge. Near critical points, the fluctuation-response relations of Section 5 predict divergent susceptibilities, connecting ensemble theory to the renormalization group and universality.

7. Historical Context

The Architects of Statistical Mechanics

The ensemble concept has a rich intellectual history spanning several decades of the late 19th and early 20th centuries.

Ludwig Boltzmann (1870s–1900s): Pioneered the statistical interpretation of entropy through his $H$-theorem (1872) and the relation $S = k \ln W$ (inscribed on his tombstone). Boltzmann's combinatorial approach to counting microstates established the microcanonical foundation, though he faced fierce opposition from the anti-atomist school led by Ernst Mach and Wilhelm Ostwald.
J. Willard Gibbs (1902): In his masterwork Elementary Principles in Statistical Mechanics, Gibbs introduced the ensemble concept with full mathematical rigor. He defined the microcanonical, canonical, and grand canonical ensembles, proved their equivalence in the thermodynamic limit, and derived the fluctuation formulas. Gibbs's abstract, axiomatic approach transcended the debates about atomism and remains the foundation of modern statistical mechanics.
Albert Einstein (1902–1904): Independently developed the foundations of statistical mechanics in a series of three papers, unaware of Gibbs's work. His 1904 paper on energy fluctuations derived the relation $\langle(\Delta E)^2\rangle = k_BT^2 C_V$ and applied it to blackbody radiation, discovering the particle-like (photon) contribution to energy fluctuations — a key step toward quantum mechanics.

8. Python Simulation: 1D Ising Model Monte Carlo

The Ising model is the simplest system exhibiting cooperative phenomena. In one dimension, it has an exact solution (no finite-temperature phase transition), but Monte Carlo simulation reveals the crossover from disordered to ordered behavior as temperature decreases. We implement the Metropolis algorithm sampling from the canonical ensemble.

1D Ising Model: Monte Carlo Simulation with Metropolis Algorithm

Python

script.py135 lines

#!/usr/bin/env python3
"""
ising_1d_monte_carlo.py
Monte Carlo simulation of the 1D Ising model using the Metropolis algorithm.
Computes energy and magnetization vs temperature, showing the crossover
from paramagnetic to ferromagnetic-like behavior.
"""
import numpy as np

# ============================================
# 1D Ising Model Parameters
# ============================================
N_spins = 100          # Number of spins in the chain
J = 1.0                # Coupling constant (ferromagnetic)
N_equil = 5000         # Equilibration sweeps
N_measure = 10000      # Measurement sweeps
N_temps = 20           # Number of temperature points

# Temperature range (in units of J/k_B)
T_min = 0.2
T_max = 5.0
temperatures = np.linspace(T_min, T_max, N_temps)

# Initialize random number generator
rng = np.random.default_rng(42)

def initialize_spins(N):
    """Initialize random spin configuration (+1 or -1)."""
    return rng.choice([-1, 1], size=N)

def compute_energy(spins, J):
    """Compute total energy with periodic boundary conditions."""
    return -J * np.sum(spins * np.roll(spins, 1))

def compute_magnetization(spins):
    """Compute total magnetization."""
    return np.sum(spins)

def metropolis_sweep(spins, beta, J):
    """Perform one full sweep of Metropolis updates."""
    N = len(spins)
    for i in range(N):
        # Compute energy change for flipping spin i
        left = spins[(i - 1) % N]
        right = spins[(i + 1) % N]
        dE = 2.0 * J * spins[i] * (left + right)

# Metropolis acceptance criterion
        if dE <= 0.0 or rng.random() < np.exp(-beta * dE):
            spins[i] *= -1
    return spins

# ============================================
# Main Simulation Loop
# ============================================
print(f"1D Ising Model Monte Carlo Simulation")
print(f"N_spins = {N_spins}, J = {J}")
print(f"Equilibration sweeps: {N_equil}, Measurement sweeps: {N_measure}")
print(f"{'='*65}")
print(f"{'T/J':>8s} {'<E>/N':>12s} {'<|M|>/N':>12s} {'C_V/N':>12s} {'chi/N':>12s}")
print(f"{'-'*65}")

energies_avg = []
mag_avg = []
cv_list = []
chi_list = []

for T in temperatures:
    beta = 1.0 / T

# Initialize spins
    spins = initialize_spins(N_spins)

# Equilibration
    for _ in range(N_equil):
        spins = metropolis_sweep(spins, beta, J)

# Measurement
    E_samples = []
    M_samples = []
    for _ in range(N_measure):
        spins = metropolis_sweep(spins, beta, J)
        E = compute_energy(spins, J)
        M = compute_magnetization(spins)
        E_samples.append(E)
        M_samples.append(abs(M))

E_arr = np.array(E_samples)
    M_arr = np.array(M_samples)

# Compute averages per spin
    E_mean = np.mean(E_arr) / N_spins
    M_mean = np.mean(M_arr) / N_spins

# Fluctuation relations (per spin)
    E2_mean = np.mean(E_arr**2)
    M2_mean = np.mean(M_arr**2)

# Heat capacity: C_V = (<E^2> - <E>^2) / (k_B T^2)
    C_V = (E2_mean - np.mean(E_arr)**2) / (T**2 * N_spins)

# Susceptibility: chi = (<M^2> - <|M|>^2) / (k_B T)
    chi = (M2_mean - np.mean(M_arr)**2) / (T * N_spins)

energies_avg.append(E_mean)
    mag_avg.append(M_mean)
    cv_list.append(C_V)
    chi_list.append(chi)

print(f"{T:8.3f} {E_mean:12.6f} {M_mean:12.6f} {C_V:12.6f} {chi:12.6f}")

# ============================================
# Exact Results for Comparison
# ============================================
print(f"\n{'='*65}")
print(f"Exact 1D Ising Results (Periodic BC, N -> infinity)")
print(f"{'='*65}")
print(f"{'T/J':>8s} {'E_exact/N':>12s} {'C_exact/N':>12s}")
print(f"{'-'*40}")

for T in temperatures:
    beta = 1.0 / T
    # Exact energy per spin: e = -J * tanh(beta*J)
    e_exact = -J * np.tanh(beta * J)
    # Exact heat capacity per spin
    c_exact = (beta * J / np.cosh(beta * J))**2
    print(f"{T:8.3f} {e_exact:12.6f} {c_exact:12.6f}")

print(f"\nKey observations:")
print(f"  - No true phase transition in 1D (Ising, Peierls argument)")
print(f"  - Crossover at T ~ J: spins become correlated below T ~ J")
print(f"  - C_V peak near T ~ 0.8J marks maximum energy fluctuations")
print(f"  - <|M|>/N -> 1 as T -> 0 (ordered); -> 0 as T -> inf (disordered)")
print(f"  - Fluctuation-response: C_V/N = (<E^2> - <E>^2) / (NkT^2)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Fortran Simulation: Quantum Harmonic Oscillator Thermodynamics

We compute the canonical partition function and all thermodynamic properties for a system of $N$ independent quantum harmonic oscillators. The single-oscillator partition function is $z_1 = [2\sinh(\beta\hbar\omega/2)]^{-1}$, and the $N$-oscillator partition function is $Z = z_1^N$. We compare the quantum results to the classical limit $E_{\text{classical}} = Nk_BT$ (equipartition).

Canonical Partition Function for Quantum Harmonic Oscillators

Fortran

program.f9097 lines

program qho_ensemble
  implicit none
  ! ============================================================
  ! Compute canonical partition function and thermodynamic
  ! properties for N quantum harmonic oscillators.
  !
  ! Z_1 = 1 / (2 * sinh(beta*hbar*omega/2))
  ! Z_N = Z_1^N
  ! A = -kT * ln(Z_N)
  ! E = -d(ln Z)/d(beta) = N*(hbar*omega/2)*coth(beta*hbar*omega/2)
  ! S = (E - A) / T
  ! C_V = dE/dT
  ! ============================================================
  integer, parameter :: n_temps = 25
  integer :: i
  double precision :: T, beta, x, hw
  double precision :: z1, ln_z, A_free, E_avg, S_ent, C_v
  double precision :: E_classical, C_classical
  double precision :: kB, hbar, omega
  integer :: N_osc
  double precision :: T_min, T_max, dT

! Physical constants (SI)
  kB    = 1.380649d-23       ! Boltzmann constant [J/K]
  hbar  = 1.054571817d-34    ! Reduced Planck constant [J s]
  omega = 1.0d13             ! Angular frequency [rad/s] (typical molecular)
  N_osc = 100                ! Number of oscillators

hw = hbar * omega          ! Energy quantum

! Temperature range
  T_min = 10.0d0             ! [K]
  T_max = 1000.0d0           ! [K]
  dT = (T_max - T_min) / dble(n_temps - 1)

write(*,'(A)') 'Quantum Harmonic Oscillator Ensemble Thermodynamics'
  write(*,'(A,I6)')    '  Number of oscillators N = ', N_osc
  write(*,'(A,ES12.4)') '  hbar*omega             = ', hw
  write(*,'(A,ES12.4)') '  Characteristic T (hw/kB)= ', hw/kB
  write(*,'(A)') ''
  write(*,'(A)') repeat('=', 90)
  write(*,'(A8, A14, A14, A14, A14, A14, A14)') &
       'T [K]', 'E/N [hw]', 'E_cl/N [hw]', 'A/N [hw]', &
       'S/NkB', 'Cv/NkB', 'Cv_cl/NkB'
  write(*,'(A)') repeat('-', 90)

do i = 1, n_temps
    T = T_min + dble(i - 1) * dT
    beta = 1.0d0 / (kB * T)
    x = beta * hw / 2.0d0        ! dimensionless parameter

! Single-oscillator partition function
    z1 = 1.0d0 / (2.0d0 * sinh(x))

! N-oscillator ln(Z)
    ln_z = dble(N_osc) * log(z1)

! Helmholtz free energy per oscillator (in units of hw)
    A_free = -kB * T * ln_z / (dble(N_osc) * hw)

! Average energy per oscillator (in units of hw)
    ! E/N = (hw/2) * coth(x)
    E_avg = 0.5d0 * cosh(x) / sinh(x)

! Classical energy per oscillator: kT (equipartition)
    E_classical = kB * T / hw

! Entropy per oscillator (in units of kB)
    ! S/NkB = x*coth(x) - ln(2*sinh(x))
    S_ent = x * cosh(x) / sinh(x) - log(2.0d0 * sinh(x))

! Heat capacity per oscillator (in units of kB)
    ! Cv/NkB = x^2 / sinh^2(x)
    C_v = x**2 / (sinh(x)**2)

! Classical heat capacity
    C_classical = 1.0d0

write(*,'(F8.1, F14.6, F14.6, F14.6, F14.6, F14.6, F14.6)') &
         T, E_avg, E_classical, A_free, S_ent, C_v, C_classical
  end do

write(*,'(A)') repeat('=', 90)
  write(*,'(A)') ''
  write(*,'(A)') 'Analysis:'
  write(*,'(A,F8.2,A)') '  Einstein temperature T_E = hw/kB = ', hw/kB, ' K'
  write(*,'(A)') '  For T >> T_E: quantum -> classical (Cv -> kB per oscillator)'
  write(*,'(A)') '  For T << T_E: quantum freezing (Cv -> 0 exponentially)'
  write(*,'(A)') '  Zero-point energy: E(T=0) = N*hw/2 (quantum, nonzero!)'
  write(*,'(A)') '  Classical limit misses zero-point energy and freezeout'
  write(*,'(A)') ''
  write(*,'(A)') 'Ensemble theory connection:'
  write(*,'(A)') '  Z = [2*sinh(beta*hw/2)]^(-N)  (canonical partition function)'
  write(*,'(A)') '  All thermodynamics derived from A = -kT*ln(Z)'
  write(*,'(A)') '  Fluctuations: <(dE)^2> = kT^2 * Cv  (Section 5 result)'

end program qho_ensemble

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Part IV Overview Physical Chemistry Home →