Part I: Ensembles | Chapter 3

The Grand Canonical Ensemble

Systems exchanging both energy and particles with a reservoir: the grand partition function, chemical potential, grand potential, and particle number fluctuations

Historical Context

Gibbs introduced the grand canonical ensemble to handle systems with variable particle number. The chemical potential \(\mu\), first defined by Gibbs in his 1876 thermodynamic work, controls the flow of particles just as temperature controls the flow of energy. The grand canonical formalism became indispensable for quantum statistical mechanics, where particle creation and annihilation are natural.

The fugacity \(z = e^{\beta\mu}\) introduced by R.H. Fowler provides an elegant parameterization that simplifies many calculations, especially for quantum gases.

1. Derivation of the Grand Canonical Distribution

Derivation 1: From the Microcanonical Ensemble

Consider a system S in contact with a reservoir R that can exchange both energy and particles. The total energy \(E_{\text{tot}}\) and particle number \(N_{\text{tot}}\)are fixed. When S is in a microstate with energy \(E_s\) and particle number \(N_s\):

\[p(E_s, N_s) \propto \Omega_R(E_{\text{tot}} - E_s, N_{\text{tot}} - N_s)\]

Taylor expanding the reservoir entropy to first order:

\[S_R \approx S_R(E_{\text{tot}}, N_{\text{tot}}) - \frac{E_s}{T} + \frac{\mu N_s}{T}\]

where we used \(\partial S/\partial E = 1/T\) and \(\partial S/\partial N = -\mu/T\). Exponentiating:

\[\boxed{p_i(E_i, N_i) = \frac{1}{\Xi} e^{-\beta(E_i - \mu N_i)}}\]

2. The Grand Partition Function

Derivation 2: Structure of \(\Xi\)

The grand partition function is the normalization constant:

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

where \(z = e^{\beta\mu}\) is the fugacity and\(Z_N\) is the canonical partition function for \(N\) particles. This beautifully separates the particle-number sum from the energy sum.

Thermodynamic Relations from \(\Xi\)

Grand potential: \(\Phi_G = -k_BT \ln \Xi = -PV\)

Average particle number: \(\langle N \rangle = z\frac{\partial \ln \Xi}{\partial z}\bigg|_{\beta,V} = k_BT\frac{\partial \ln \Xi}{\partial \mu}\bigg|_{\beta,V}\)

Average energy: \(\langle E \rangle = -\frac{\partial \ln \Xi}{\partial \beta}\bigg|_{z,V} + \mu\langle N\rangle\)

Pressure: \(PV = k_BT \ln \Xi\)

Entropy: \(S = k_B\left(\ln \Xi + \beta\langle E\rangle - \beta\mu\langle N\rangle\right)\)

Derivation 3: Grand Potential and the Euler Relation

The grand potential is the Legendre transform of the Helmholtz free energy:

\[\Phi_G = F - \mu N = E - TS - \mu N\]

From the Euler relation for extensive quantities (\(E = TS - PV + \mu N\)):

\[\Phi_G = -PV\]

This remarkable result means that \(k_BT\ln\Xi = PV\), connecting the grand partition function directly to the equation of state.

3. Particle Number Fluctuations

Derivation 4: Compressibility Sum Rule

The variance in particle number is:

\[\langle (\Delta N)^2 \rangle = \frac{1}{\beta}\frac{\partial \langle N \rangle}{\partial \mu}\bigg|_{T,V} = k_BT\left(\frac{\partial \langle N \rangle}{\partial \mu}\right)_{T,V}\]

This can be related to the isothermal compressibility \(\kappa_T = -\frac{1}{V}\frac{\partial V}{\partial P}\big|_T\):

\[\boxed{\frac{\langle(\Delta N)^2\rangle}{\langle N \rangle} = \frac{k_BT}{V}\langle N \rangle \kappa_T = \frac{\rho k_BT \kappa_T}{1}}\]

For an ideal gas, \(\kappa_T = 1/P = V/(Nk_BT)\), giving Poisson-like fluctuations:\(\langle(\Delta N)^2\rangle = \langle N \rangle\). Near a critical point where\(\kappa_T\) diverges, particle fluctuations become anomalously large, leading to critical opalescence.

4. Applications

Ideal Gas in the Grand Canonical Ensemble

For the classical ideal gas with \(Z_N = V^N/(N!\lambda_{dB}^{3N})\):

\[\Xi = \sum_{N=0}^{\infty} \frac{z^N}{N!}\left(\frac{V}{\lambda_{dB}^3}\right)^N = \exp\left(\frac{zV}{\lambda_{dB}^3}\right)\]

\[\langle N \rangle = z\frac{\partial}{\partial z}\ln\Xi = \frac{zV}{\lambda_{dB}^3}, \qquad PV = k_BT\ln\Xi = \langle N \rangle k_BT\]

Adsorption: The Langmuir Isotherm

Consider \(M\) independent adsorption sites, each either empty or occupied by one particle with binding energy \(-\epsilon\). The grand partition function for a single site is:

\[\xi = 1 + ze^{\beta\epsilon}\]

For \(M\) independent sites: \(\Xi = \xi^M\). The mean occupation fraction (coverage) is:

\[\theta = \frac{\langle N \rangle}{M} = \frac{ze^{\beta\epsilon}}{1 + ze^{\beta\epsilon}} = \frac{P/P_0}{1 + P/P_0}\]

where \(P_0 = k_BT\lambda_{dB}^{-3}e^{-\beta\epsilon}\). This is the celebrated Langmuir isotherm.

5. Preview: Quantum Occupation Numbers

Derivation 5: Single-Mode Grand Partition Function

For a single quantum mode with energy \(\epsilon\), the grand partition function depends on the quantum statistics:

Bosons (\(n = 0,1,2,\ldots\))

\[\xi_{\text{BE}} = \sum_{n=0}^{\infty} (ze^{-\beta\epsilon})^n = \frac{1}{1 - ze^{-\beta\epsilon}}\]

\[\langle n \rangle_{\text{BE}} = \frac{1}{z^{-1}e^{\beta\epsilon} - 1}\]

Fermions (\(n = 0,1\))

\[\xi_{\text{FD}} = 1 + ze^{-\beta\epsilon}\]

\[\langle n \rangle_{\text{FD}} = \frac{1}{z^{-1}e^{\beta\epsilon} + 1}\]

The total grand partition function is a product over all single-particle states:\(\Xi = \prod_k \xi_k\), giving \(\ln\Xi = \sum_k \ln\xi_k\).

6. Computational Exploration

This simulation explores the grand canonical ensemble, demonstrating particle number fluctuations, the Langmuir isotherm, and quantum occupation numbers.

Grand Canonical Ensemble: Particle Fluctuations, Langmuir Isotherm, and Quantum Statistics

Python

script.py120 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# Grand Canonical Ensemble Explorations
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Grand Canonical Ensemble', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#1a0a0a')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('#fca5a5')
    for spine in ax.spines.values():
        spine.set_color('#333')

# --- Panel 1: Particle number distribution for ideal gas ---
ax1 = axes[0, 0]
for avg_N in [10, 30, 100]:
    N_range = np.arange(max(0, avg_N - 4 * int(np.sqrt(avg_N))),
                         avg_N + 4 * int(np.sqrt(avg_N)) + 1)
    # Poisson distribution
    from math import lgamma
    log_P = N_range * np.log(avg_N) - avg_N - np.array([lgamma(n + 1) for n in N_range])
    P = np.exp(log_P)
    ax1.plot(N_range, P, linewidth=1.5, label='<N> = ' + str(avg_N))
ax1.set_xlabel('N')
ax1.set_ylabel('P(N)')
ax1.set_title('Particle Number Distribution (Ideal Gas)')
ax1.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Langmuir isotherm ---
ax2 = axes[0, 1]
P_over_P0 = np.linspace(0, 10, 300)
theta_langmuir = P_over_P0 / (1 + P_over_P0)
ax2.plot(P_over_P0, theta_langmuir, color='#f87171', linewidth=2, label='Langmuir')
# BET (multilayer) for comparison
c_BET = 5.0
theta_BET = c_BET * P_over_P0 / ((1 - P_over_P0 / 10) * (1 - P_over_P0 / 10 + c_BET * P_over_P0 / 10))
valid = P_over_P0 < 9
ax2.plot(P_over_P0[valid], np.clip(theta_BET[valid], 0, 3), '--', color='#a78bfa',
         linewidth=1.5, label='BET (multilayer)', alpha=0.7)
ax2.axhline(y=1, color='white', linestyle=':', alpha=0.3)
ax2.set_xlabel('P / P_0')
ax2.set_ylabel('Coverage theta')
ax2.set_title('Langmuir Adsorption Isotherm')
ax2.set_ylim(0, 1.5)
ax2.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Quantum occupation numbers ---
ax3 = axes[1, 0]
eps_minus_mu = np.linspace(-3, 5, 500)
for kT in [0.2, 0.5, 1.0, 2.0]:
    beta_val = 1.0 / kT
    n_BE = np.where(eps_minus_mu > 0.01,
                    1.0 / (np.exp(beta_val * eps_minus_mu) - 1), np.nan)
    n_FD = 1.0 / (np.exp(beta_val * eps_minus_mu) + 1)
    if kT == 1.0:
        ax3.plot(eps_minus_mu, n_FD, color='#60a5fa', linewidth=2, label='FD (kT=1)')
        ax3.plot(eps_minus_mu[eps_minus_mu > 0.01],
                 n_BE[eps_minus_mu > 0.01], color='#f87171', linewidth=2, label='BE (kT=1)')
    elif kT == 0.2:
        ax3.plot(eps_minus_mu, n_FD, color='#60a5fa', linewidth=1, alpha=0.5, linestyle='--', label='FD (kT=0.2)')
ax3.set_xlabel('(epsilon - mu) / eps_0')
ax3.set_ylabel('<n>')
ax3.set_title('Quantum Occupation Numbers')
ax3.set_ylim(0, 5)
ax3.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Compressibility and fluctuations near critical point ---
ax4 = axes[1, 1]
# Van der Waals gas: kappa_T diverges near T_c
T_over_Tc = np.linspace(0.85, 2.0, 300)
# Approximate kappa_T/kappa_0 near critical point ~ |T - Tc|^{-gamma}
gamma_exp = 1.0  # mean-field
kappa_ratio = np.where(T_over_Tc > 1.0,
                       1.0 / np.abs(T_over_Tc - 1.0)**gamma_exp,
                       1.0 / np.abs(T_over_Tc - 1.0)**gamma_exp)
kappa_ratio = np.clip(kappa_ratio, 0, 100)
ax4.semilogy(T_over_Tc, kappa_ratio, color='#fb923c', linewidth=2)
ax4.axvline(x=1.0, color='#fbbf24', linestyle='--', alpha=0.5, label='T = T_c')
ax4.set_xlabel('T / T_c')
ax4.set_ylabel('kappa_T / kappa_0 (log scale)')
ax4.set_title('Compressibility Divergence Near T_c')
ax4.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")

# Numerical results
print("\n" + "=" * 60)
print("GRAND CANONICAL ENSEMBLE: NUMERICAL RESULTS")
print("=" * 60)

print("\n--- Ideal Gas Particle Fluctuations ---")
for avg_N in [10, 100, 1000, 10000]:
    rel_fluct = 1.0 / np.sqrt(avg_N)
    print("  <N> = " + str(avg_N) + ": sigma_N/<N> = " + str(round(rel_fluct, 6)) + ", sigma_N = " + str(round(np.sqrt(avg_N), 2)))

print("\n--- Langmuir Isotherm ---")
for p_ratio in [0.1, 0.5, 1.0, 5.0, 10.0]:
    theta = p_ratio / (1 + p_ratio)
    print("  P/P0 = " + str(p_ratio) + ": theta = " + str(round(theta, 4)))

print("\n--- Occupation Numbers at kT = 1 ---")
for eps_mu in [0.1, 0.5, 1.0, 2.0, 5.0]:
    n_be = 1.0 / (np.exp(eps_mu) - 1)
    n_fd = 1.0 / (np.exp(eps_mu) + 1)
    n_mb = np.exp(-eps_mu)
    print("  (eps-mu) = " + str(eps_mu) + ": n_BE = " + str(round(n_be, 4)) + ", n_FD = " + str(round(n_fd, 4)) + ", n_MB = " + str(round(n_mb, 4)))
print("\n" + "=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Summary and Key Results

Core Formulas

Distribution: \(p \propto e^{-\beta(E-\mu N)}\)
Grand potential: \(\Phi_G = -k_BT\ln\Xi = -PV\)
Fluctuations: \(\langle(\Delta N)^2\rangle = k_BT\partial\langle N\rangle/\partial\mu\)
Fugacity expansion: \(\Xi = \sum_N z^N Z_N\)

Physical Insights

Chemical potential controls particle flow
Ideal gas: Poisson particle statistics
Compressibility links to density fluctuations
Natural framework for quantum gases

← Canonical Ensemble Quantum Statistics →