Part IV: Advanced Topics | Chapter 3

Non-Equilibrium Statistical Mechanics

Master equation, Fokker-Planck equation, Langevin dynamics, detailed balance, and linear response theory

Historical Context

While equilibrium statistical mechanics was largely complete by the early 20th century, the non-equilibrium theory remains an active frontier. Albert Einstein's 1905 theory of Brownian motion and Paul Langevin's stochastic equation (1908) launched the field. Adriaan Fokker (1914) and Max Planck (1917) developed the partial differential equation for probability evolution, and Lars Onsager's reciprocal relations (1931, Nobel Prize 1968) established the first general results for systems near equilibrium.

Ryogo Kubo's fluctuation-dissipation theorem (1966) and linear response theory provided the bridge between equilibrium fluctuations and non-equilibrium transport, showing that transport coefficients can be computed from equilibrium correlation functions. Modern developments include the Jarzynski equality and Crooks theorem, extending thermodynamics to far-from-equilibrium processes.

1. The Master Equation

Derivation 1: Probability Evolution for Markov Processes

Consider a system that can be in discrete states \(\{n\}\) with transition rates \(W_{mn}\) (probability per unit time to go from state \(n\)to state \(m\)). The probability \(P_n(t)\) evolves according to:

\[\boxed{\frac{dP_n}{dt} = \sum_{m \ne n}\left[W_{nm}P_m(t) - W_{mn}P_n(t)\right]}\]

The first term (gain) represents transitions into state \(n\); the second term (loss) represents transitions out. This can be written in matrix form:

\[\frac{d\mathbf{P}}{dt} = \mathbf{W}\mathbf{P}, \qquad W_{nm} = \begin{cases} W_{nm} & n \ne m \\ -\sum_{k \ne n}W_{kn} & n = m \end{cases}\]

The steady state satisfies \(\mathbf{W}\mathbf{P}^{ss} = 0\). If the transition rates satisfy detailed balance, the equilibrium distribution\(P_n^{eq} = e^{-\beta E_n}/Z\) is the unique steady state.

2. Detailed Balance

Derivation 2: Equilibrium Condition

Detailed balance requires that at equilibrium, each transition is individually balanced by its reverse:

\[\boxed{W_{nm}P_m^{eq} = W_{mn}P_n^{eq} \quad \text{for all } n, m}\]

For a thermal system: \(W_{nm}/W_{mn} = e^{-\beta(E_n - E_m)}\). This is stronger than the steady-state condition (which only requires net zero flow for each state, not each pair). Detailed balance guarantees time-reversal symmetry of the dynamics and ensures relaxation to the Boltzmann distribution.

The Metropolis algorithm satisfies detailed balance with the choice:

\[W_{nm} = \min\left(1, e^{-\beta(E_n - E_m)}\right)\]

Systems that violate detailed balance include driven systems (e.g., molecular motors, active matter) and systems with non-conservative forces. These exhibit persistent probability currents even in steady state.

3. The Fokker-Planck Equation

Derivation 3: Continuum Limit of the Master Equation

For continuous variables \(x\), the master equation becomes an integral equation. Expanding for small jumps (Kramers-Moyal expansion) and keeping the first two terms gives the Fokker-Planck equation:

\[\boxed{\frac{\partial P(x,t)}{\partial t} = -\frac{\partial}{\partial x}\left[A(x)P\right] + \frac{1}{2}\frac{\partial^2}{\partial x^2}\left[B(x)P\right]}\]

where \(A(x) = \langle\Delta x\rangle/\Delta t\) is the drift coefficient and\(B(x) = \langle(\Delta x)^2\rangle/\Delta t\) is the diffusion coefficient. For a particle in a potential \(V(x)\) with friction \(\gamma\)at temperature \(T\):

\[A(x) = -\frac{1}{\gamma}\frac{dV}{dx}, \qquad B = \frac{2k_BT}{\gamma}\]

The steady-state solution is the Boltzmann distribution:\(P^{ss}(x) \propto e^{-V(x)/(k_BT)}\), provided the Einstein relation\(D = k_BT/\gamma\) holds. This connects the fluctuation strength (diffusion) to the dissipation (friction) -- a manifestation of the fluctuation-dissipation theorem.

4. The Langevin Equation

Derivation 4: Stochastic Dynamics

The Langevin equation describes the trajectory of a single Brownian particle:

\[\boxed{m\ddot{x} = -\gamma\dot{x} - \frac{dV}{dx} + \sqrt{2\gamma k_BT}\,\xi(t)}\]

where \(\xi(t)\) is Gaussian white noise: \(\langle\xi(t)\rangle = 0\),\(\langle\xi(t)\xi(t')\rangle = \delta(t - t')\). In the overdamped limit (\(m\ddot{x} \approx 0\)):

\[\gamma\dot{x} = -\frac{dV}{dx} + \sqrt{2\gamma k_BT}\,\xi(t)\]

The Langevin and Fokker-Planck descriptions are equivalent: the Langevin equation gives individual trajectories, while the Fokker-Planck equation gives the probability distribution. The connection is that the ensemble of Langevin trajectories generates the Fokker-Planck probability density.

Einstein Relation and Brownian Motion

Einstein showed in 1905 that the diffusion coefficient of a Brownian particle is\(D = k_BT/\gamma\). Combined with Stokes' law \(\gamma = 6\pi\eta a\)for a sphere of radius \(a\) in a fluid of viscosity \(\eta\), this gives the Stokes-Einstein relation: \(D = k_BT/(6\pi\eta a)\). Perrin's experimental verification of \(\langle x^2\rangle = 2Dt\) (1908) provided definitive evidence for the existence of atoms.

5. Linear Response Theory

Derivation 5: Kubo Formula

Consider a system in equilibrium perturbed by a time-dependent external field\(f(t)\) coupled to observable \(A\):\(H(t) = H_0 - f(t)A\). The linear response of observable \(B\) is:

\[\langle\delta B(t)\rangle = \int_{-\infty}^{t}\chi_{BA}(t - t')\,f(t')\,dt'\]

The response function (susceptibility) is given by the Kubo formula:

\[\boxed{\chi_{BA}(t) = \frac{1}{k_BT}\frac{d}{dt}\langle B(t)A(0)\rangle_0\,\Theta(t)}\]

where \(\langle\cdots\rangle_0\) denotes the equilibrium average and\(\Theta(t)\) is the step function ensuring causality. In frequency space:

\[\chi(\omega) = \chi'(\omega) + i\chi''(\omega)\]

The fluctuation-dissipation theorem (FDT) relates the dissipative part to equilibrium fluctuations:

\[\boxed{\chi''(\omega) = \frac{\omega}{2k_BT}\tilde{C}(\omega)}\]

where \(\tilde{C}(\omega)\) is the Fourier transform of the equilibrium correlation function \(C(t) = \langle A(t)A(0)\rangle_0\). This remarkable result means transport coefficients (conductivity, viscosity, diffusion) can be computed from equilibrium simulations via Green-Kubo relations.

6. Applications

Key Applications

Electrical conductivity: The Kubo formula gives the conductivity as an integral of the current-current correlation function:\(\sigma = \frac{1}{k_BTV}\int_0^{\infty}\langle\mathbf{J}(t)\cdot\mathbf{J}(0)\rangle\,dt\).

Kramers escape problem: The rate of escape from a metastable potential well is \(r \sim \omega_0 e^{-\Delta V/(k_BT)}\) (Arrhenius), with prefactor corrections from the Fokker-Planck equation.

Molecular dynamics: Langevin thermostats (adding friction and noise) are used in molecular dynamics simulations to sample the canonical ensemble.

Active matter: Self-propelled particles (bacteria, active colloids) are described by modified Langevin equations that violate detailed balance, leading to phenomena like motility-induced phase separation.

7. Computational Exploration

This simulation demonstrates Langevin dynamics, Fokker-Planck probability evolution, Brownian motion statistics, and linear response.

Non-Equilibrium: Langevin Dynamics, Fokker-Planck, Brownian Motion, and Autocorrelation

Python

script.py185 lines

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

# ============================================================
# Non-Equilibrium Stat Mech: Langevin, Fokker-Planck, Response
# ============================================================
np.random.seed(42)

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Non-Equilibrium Statistical Mechanics', 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: Langevin dynamics in a double-well potential ---
ax1 = axes[0, 0]

# V(x) = x^4 - 2x^2
gamma = 1.0
kBT = 0.5
dt = 0.01
n_steps = 50000

x = np.zeros(n_steps)
x[0] = -1.0
noise_amp = np.sqrt(2 * gamma * kBT * dt)

for i in range(1, n_steps):
    dVdx = 4 * x[i-1]**3 - 4 * x[i-1]
    x[i] = x[i-1] + (-dVdx / gamma) * dt + noise_amp * np.random.randn()

t_arr = np.arange(n_steps) * dt
ax1.plot(t_arr[:10000], x[:10000], color='#60a5fa', linewidth=0.3, alpha=0.7)

# Overlay potential
ax1_twin = ax1.twinx()
x_pot = np.linspace(-2, 2, 200)
V_pot = x_pot**4 - 2 * x_pot**2
ax1_twin.plot(x_pot, V_pot, color='#f87171', linewidth=2, alpha=0.5)
ax1_twin.set_ylabel('V(x)', color='#f87171')
ax1_twin.tick_params(axis='y', colors='#f87171')
ax1_twin.set_ylim(-2, 5)

ax1.set_xlabel('Time')
ax1.set_ylabel('x(t)', color='#60a5fa')
ax1.tick_params(axis='y', colors='#60a5fa')
ax1.set_title('Langevin in Double Well (kBT=' + str(kBT) + ')')
ax1.set_xlim(0, 100)

# --- Panel 2: Fokker-Planck steady state vs Boltzmann ---
ax2 = axes[0, 1]

# Histogram of Langevin trajectory = steady-state P(x)
x_bins = np.linspace(-2.5, 2.5, 100)
hist, bin_edges = np.histogram(x[5000:], bins=x_bins, density=True)
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])

# Boltzmann prediction
V_arr = bin_centers**4 - 2 * bin_centers**2
P_boltz = np.exp(-V_arr / kBT)
P_boltz /= np.trapezoid(P_boltz, bin_centers)

ax2.bar(bin_centers, hist, width=bin_centers[1]-bin_centers[0], alpha=0.5, color='#60a5fa', label='Langevin histogram')
ax2.plot(bin_centers, P_boltz, color='#f87171', linewidth=2.5, label='Boltzmann: exp(-V/kBT)')
ax2.set_xlabel('x')
ax2.set_ylabel('P(x)')
ax2.set_title('Steady State = Boltzmann Distribution')
ax2.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Mean squared displacement (Brownian motion) ---
ax3 = axes[1, 0]

# Free Brownian motion: <x^2> = 2Dt
n_particles = 200
n_time = 5000
dt_bm = 0.01
D_vals = [0.1, 0.5, 1.0, 2.0]

for D, color, label in zip(D_vals,
                            ['#f87171', '#60a5fa', '#34d399', '#a78bfa'],
                            ['D = 0.1', 'D = 0.5', 'D = 1.0', 'D = 2.0']):
    # Simulate many particles
    positions = np.zeros((n_particles, n_time))
    noise = np.sqrt(2 * D * dt_bm) * np.random.randn(n_particles, n_time - 1)
    for i in range(1, n_time):
        positions[:, i] = positions[:, i-1] + noise[:, i-1]

msd = np.mean(positions**2, axis=0)
    t_bm = np.arange(n_time) * dt_bm
    ax3.plot(t_bm, msd, color=color, linewidth=2, label=label)
    ax3.plot(t_bm, 2 * D * t_bm, '--', color=color, linewidth=1, alpha=0.5)

ax3.set_xlabel('Time')
ax3.set_ylabel('<x^2>')
ax3.set_title('Mean Squared Displacement (solid: MC, dashed: 2Dt)')
ax3.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Autocorrelation and relaxation ---
ax4 = axes[1, 1]

# Ornstein-Uhlenbeck process: dx = -x/tau * dt + sqrt(2D) * dW
tau = 5.0
D_ou = 1.0
n_ou = 100000
dt_ou = 0.05
x_ou = np.zeros(n_ou)

for i in range(1, n_ou):
    x_ou[i] = x_ou[i-1] - x_ou[i-1]/tau * dt_ou + np.sqrt(2*D_ou*dt_ou) * np.random.randn()

# Compute autocorrelation
max_lag = 2000
lags = np.arange(max_lag)
x_mean = np.mean(x_ou[10000:])
x_fluct = x_ou[10000:] - x_mean
C0 = np.mean(x_fluct**2)

autocorr = []
for lag in lags:
    if lag == 0:
        autocorr.append(1.0)
    else:
        c = np.mean(x_fluct[:-lag] * x_fluct[lag:]) / C0
        autocorr.append(c)

autocorr = np.array(autocorr)
t_lags = lags * dt_ou

ax4.plot(t_lags, autocorr, color='#f87171', linewidth=2, label='Numerical C(t)/C(0)')
ax4.plot(t_lags, np.exp(-t_lags / tau), 'w--', linewidth=1.5, alpha=0.6, label='exp(-t/tau), tau=' + str(tau))
ax4.set_xlabel('Time lag t')
ax4.set_ylabel('C(t) / C(0)')
ax4.set_title('Autocorrelation (Ornstein-Uhlenbeck)')
ax4.set_xlim(0, 40)
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("NON-EQUILIBRIUM STAT MECH: NUMERICAL RESULTS")
print("=" * 60)

print("\n--- Brownian Motion ---")
print("  Einstein relation: D = kBT / gamma")
print("  Mean squared displacement: <x^2> = 2Dt (1D)")
for D in D_vals:
    print("  D = " + str(D) + ": <x^2>(t=10) = " + str(round(2*D*10, 2)))

print("\n--- Ornstein-Uhlenbeck Process ---")
print("  tau = " + str(tau))
print("  D = " + str(D_ou))
print("  Steady-state variance: <x^2> = D*tau = " + str(D_ou*tau))
print("  Measured variance: " + str(round(np.var(x_ou[10000:]), 3)))

print("\n--- Double-Well Kramers Escape ---")
V_barrier = 1.0  # V(0) - V(-1) = 0 - (-1) = 1
rate_theory = np.sqrt(4 * 12) / (2 * np.pi) * np.exp(-V_barrier / kBT)
print("  Barrier height: Delta V = " + str(V_barrier))
print("  kBT = " + str(kBT))
print("  Kramers rate estimate: r ~ " + "{:.4e}".format(rate_theory))

print("\n--- Fluctuation-Dissipation Theorem ---")
print("  chi''(omega) = omega/(2kBT) * C_tilde(omega)")
print("  Green-Kubo: sigma = (1/kBTV) * integral <J(t).J(0)> dt")
print("  Einstein: D = lim (t->inf) <x^2>/(2t)")

print("\n--- Key Relations ---")
print("  Fokker-Planck: dP/dt = -d/dx[A*P] + (1/2)d^2/dx^2[B*P]")
print("  Langevin: gamma*dx = -dV/dx*dt + sqrt(2*gamma*kBT)*dW")
print("  Detailed balance: W_nm * P_m^eq = W_mn * P_n^eq")
print("\n" + "=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

8. Summary and Key Results

Core Formulas

Master eq: \(dP_n/dt = \sum_m(W_{nm}P_m - W_{mn}P_n)\)
Fokker-Planck: \(\partial_t P = -\partial_x(AP) + \frac{1}{2}\partial_x^2(BP)\)
FDT: \(\chi''(\omega) = \frac{\omega}{2k_BT}\tilde{C}(\omega)\)
Einstein: \(D = k_BT/\gamma\)

Physical Insights

Detailed balance ensures relaxation to Boltzmann
Fluctuation = dissipation (noise strength tied to friction)
Transport coefficients from equilibrium correlations
Langevin and Fokker-Planck are equivalent descriptions

← Universality Fluctuation Theorems →