Chapter 11.3: Tracy-Widom Distribution & KPZ Universality

Traffic Fluctuations Are Not Gaussian

The central limit theorem predicts Gaussian fluctuations for sums of independent random variables. But the ASEP is a strongly correlated system: the exclusion constraint creates long-range correlations. The remarkable result, proved by Johansson (2000) and Tracy-Widom (2009), is that the current fluctuations in the ASEP belong to the KPZ universality class and follow the Tracy-Widom distribution—the same distribution that governs the largest eigenvalue of random matrices.

This has practical consequences for traffic: the distribution of vehicle throughput at a measurement point is not Gaussian but exhibits characteristic skewness and heavy tails predicted by the Tracy-Widom distribution.

11.3.1 Current Fluctuations in the ASEP

Let $Q(t)$ be the total number of particles that have crossed a fixed bond up to time $t$. The mean current is:

$$\langle Q(t) \rangle = J \cdot t, \qquad J = (p-q)\rho(1-\rho)$$

The fluctuations of $Q(t)$ grow anomalously. For the TASEP with step initial condition (all sites occupied for$x \leq 0$, empty for$x > 0$):

$$\frac{Q(t) - \langle Q \rangle t}{(\Gamma t)^{1/3}} \xrightarrow{d} F_{\text{GUE}}$$

where $F_{\text{GUE}}$ is the Tracy-Widom GUE distribution (the distribution of the largest eigenvalue of a random matrix from the Gaussian Unitary Ensemble), and$\Gamma$ is a model-dependent constant.

Key Properties of Tracy-Widom

Skewness:$\gamma_1 \approx 0.2241$ (positive skew—large fluctuations more likely than small ones)
Excess kurtosis:$\gamma_2 \approx 0.0934$ (heavier tails than Gaussian)
Left tail:$F(s) \sim \exp(-|s|^3/12)$ (faster decay than Gaussian)
Right tail:$1-F(s) \sim \exp(-2s^{3/2}/3)$ (slower decay than Gaussian)

The $t^{1/3}$ scaling of fluctuations (compared to the $t^{1/2}$ Gaussian scaling) is the signature of KPZ universality. The exponent $1/3$ is exact and universal—it does not depend on the microscopic details of the model.

11.3.2 KPZ Universality Class

The Kardar-Parisi-Zhang (KPZ) equation describes the growth of an interface$h(x,t)$:

$$\frac{\partial h}{\partial t} = \nu \nabla^2 h + \frac{\lambda}{2}(\nabla h)^2 + \eta(x,t)$$

The connection to the ASEP is through the mapping$h(x,t) = \sum_{y \leq x} [\tau_y(t) - \rho]$—the integrated density fluctuation is a height function. The ASEP particle current corresponds to interface growth. The three KPZ exponents are:

$$\text{Growth:} \quad \beta = \frac{1}{3}, \qquad \text{Roughness:} \quad \alpha = \frac{1}{2}, \qquad \text{Dynamic:} \quad z = \frac{3}{2}$$

The dynamic exponent $z = 3/2$ gives the relaxation time scaling $\tau \sim L^{3/2}$that we encountered in the TASEP spectral gap.

Quantum Group Connection

The ASEP Markov generator commutes with the action of the quantum group$U_q(\mathfrak{sl}_2)$ with deformation parameter:

$$q_{\text{deform}} = \sqrt{\frac{q}{p}}$$

When $q = 0$ (TASEP),$q_{\text{deform}} = 0$ and we are at the crystal limit of the quantum group. When$q = p$ (SSEP),$q_{\text{deform}} = 1$ and the quantum group reduces to the classical $\mathfrak{sl}_2$. The full crossover from KPZ to diffusive scaling is encoded in this deformation parameter.

11.3.3 Full Crossover: KPZ to Diffusive

The spectral gap of the ASEP generator with periodic boundary conditions and$N$ particles on$L$ sites, at density$\rho = N/L$, is:

$$\Delta = p + q - 2\sqrt{pq}\cos\!\left(\frac{\pi}{L}\right) \approx \frac{\pi^2 \sqrt{pq}}{L^2}$$

For open boundary TASEP, the gap is$\Delta \sim L^{-3/2}$. The crossover between these two regimes is controlled by:

$$\tau_{\text{relax}} \sim \begin{cases} L^{3/2} & \text{if } q/p \ll 1/L \quad \text{(KPZ regime)} \\ L^2 / \sqrt{pq} & \text{if } q/p \gg 1/L \quad \text{(diffusive regime)} \end{cases}$$

11.3.4 Model Hierarchy

The models we have studied form a natural hierarchy:

ASEP (p, q, open boundaries)

↓ q → 0

TASEP — 3 phases, KPZ universality

↓ p → q

SSEP — diffusive, Gaussian fluctuations

ASEP + Langmuir (add ω_A, ω_D)

↓ Ω > 0

5 phases, domain wall localisation

Each level adds physical complexity. The TASEP captures the essence of unidirectional traffic with boundary effects. The ASEP adds backward motion (lane changing, reversals). Langmuir coupling adds bulk sources and sinks (on/off ramps along the road). At each level, exact results from integrable systems theory provide non-perturbative predictions.

11.3.5 Simulation: Current Fluctuations and Tracy-Widom

TASEP Current Fluctuations and Tracy-Widom Comparison

Python

script.py173 lines

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
from scipy.stats import norm

np.random.seed(42)

def simulate_tasep_current(L, rho_init, T, n_trials):
    """Simulate TASEP with periodic BC and measure current fluctuations."""
    currents = []

for trial in range(n_trials):
        # Initialize with Bernoulli occupation
        lattice = (np.random.random(L) < rho_init).astype(int)
        Q = 0  # cumulative current across bond (0, 1)

for t in range(T):
            # Random sequential update
            i = np.random.randint(0, L)
            j = (i + 1) % L
            if lattice[i] == 1 and lattice[j] == 0:
                lattice[i] = 0
                lattice[j] = 1
                if i == 0:
                    Q += 1  # count crossings of bond 0->1

currents.append(Q)

return np.array(currents)

# ── Parameters ──
L = 50
rho = 0.5  # half-filled (maximal current regime)
T = 100_000
n_trials = 500

print("Running TASEP simulations for current fluctuation statistics...")
currents = simulate_tasep_current(L, rho, T, n_trials)

# Expected mean current: J = rho(1-rho) * T/L (per Monte Carlo step)
J_theory = rho * (1 - rho)
Q_mean = np.mean(currents)
Q_std = np.std(currents)

# KPZ scaling: fluctuations ~ t^{1/3}
# Rescale: chi = (Q - <Q>) / (Gamma * t)^{1/3}
# We estimate Gamma from the data
t_eff = T / L  # effective time per site
Gamma_est = (Q_std**3 / t_eff)  # rough estimate
chi = (currents - Q_mean) / Q_std  # standardised

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.patch.set_facecolor('#0a0a0a')

# Panel 1: Histogram of current fluctuations
ax1 = axes[0, 0]
ax1.set_facecolor('#0a0a0a')
ax1.tick_params(colors='white')
for spine in ax1.spines.values():
    spine.set_color('#334155')

counts, bins, _ = ax1.hist(chi, bins=60, density=True, color='#34d399', alpha=0.7,
                            edgecolor='#065f46', label='TASEP simulation')
x_gauss = np.linspace(-4, 4, 300)
ax1.plot(x_gauss, norm.pdf(x_gauss), color='#f87171', linewidth=2, linestyle='--',
        label='Gaussian')

# Approximate Tracy-Widom GUE (using shifted/scaled approximation)
# TW_GUE has mean ≈ -1.771, std ≈ 0.813, skew ≈ 0.224
# We plot a skewed approximation
from scipy.stats import skewnorm
a_skew = 2.5  # approximate skewness parameter
tw_approx = skewnorm.pdf(x_gauss, a_skew, loc=-0.5, scale=0.9)
tw_approx = tw_approx / tw_approx.max() * counts.max() * 1.1
ax1.plot(x_gauss, skewnorm.pdf(x_gauss, a_skew, loc=-0.3, scale=0.85),
        color='#fbbf24', linewidth=2, label='Tracy-Widom (approx)')

ax1.set_xlabel('Standardised current χ', color='white')
ax1.set_ylabel('Probability density', color='white')
ax1.set_title('Current Fluctuations', color='#6ee7b7', fontsize=13)
ax1.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)

# Panel 2: Q-Q plot against Gaussian
ax2 = axes[0, 1]
ax2.set_facecolor('#0a0a0a')
ax2.tick_params(colors='white')
for spine in ax2.spines.values():
    spine.set_color('#334155')

sorted_chi = np.sort(chi)
n = len(sorted_chi)
theoretical_quantiles = norm.ppf((np.arange(1, n+1) - 0.5) / n)
ax2.scatter(theoretical_quantiles[::10], sorted_chi[::10], color='#34d399', s=8, alpha=0.6)
ax2.plot([-3, 3], [-3, 3], color='#f87171', linewidth=2, linestyle='--', label='Gaussian reference')
ax2.set_xlabel('Gaussian quantiles', color='white')
ax2.set_ylabel('Sample quantiles', color='white')
ax2.set_title('Q-Q Plot (deviation from Gaussian)', color='#6ee7b7', fontsize=13)
ax2.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

# Panel 3: Cumulants
ax3 = axes[1, 0]
ax3.set_facecolor('#0a0a0a')
ax3.tick_params(colors='white')
for spine in ax3.spines.values():
    spine.set_color('#334155')

from scipy.stats import kurtosis, skew
measured_skew = skew(chi)
measured_kurt = kurtosis(chi)  # excess kurtosis
tw_skew = 0.2241
tw_kurt = 0.0934

labels = ['Skewness', 'Excess Kurtosis']
measured = [measured_skew, measured_kurt]
tw_values = [tw_skew, tw_kurt]
gaussian = [0, 0]

x_pos = np.arange(len(labels))
width = 0.25
bars1 = ax3.bar(x_pos - width, measured, width, color='#34d399', label='Measured', edgecolor='#065f46')
bars2 = ax3.bar(x_pos, tw_values, width, color='#fbbf24', label='Tracy-Widom GUE', edgecolor='#92400e')
bars3 = ax3.bar(x_pos + width, gaussian, width, color='#f87171', label='Gaussian', edgecolor='#7f1d1d')

ax3.set_xticks(x_pos)
ax3.set_xticklabels(labels, color='white')
ax3.set_ylabel('Value', color='white')
ax3.set_title('Non-Gaussian Signatures', color='#6ee7b7', fontsize=13)
ax3.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)
ax3.axhline(y=0, color='white', linestyle='-', alpha=0.2)

# Panel 4: Fluctuation scaling with time
ax4 = axes[1, 1]
ax4.set_facecolor('#0a0a0a')
ax4.tick_params(colors='white')
for spine in ax4.spines.values():
    spine.set_color('#334155')

# Measure variance at different times
T_values = [10000, 30000, 60000, 100000]
var_values = []
for T_val in T_values:
    curr = simulate_tasep_current(L, rho, T_val, 100)
    var_values.append(np.var(curr))

T_arr = np.array(T_values) / L
var_arr = np.array(var_values)

ax4.loglog(T_arr, var_arr, 'o-', color='#34d399', linewidth=2, markersize=8, label='Measured Var(Q)')

# Reference slopes
t_ref = np.array([T_values[0]/L, T_values[-1]/L])
# KPZ: Var ~ t^{2/3}
ax4.loglog(t_ref, var_arr[0] * (t_ref/t_ref[0])**(2/3), '--', color='#fbbf24',
          linewidth=2, label=r'$t^{2/3}$ (KPZ)')
# Gaussian: Var ~ t
ax4.loglog(t_ref, var_arr[0] * (t_ref/t_ref[0])**(1.0), ':', color='#f87171',
          linewidth=2, label=r'$t^1$ (Gaussian)')

ax4.set_xlabel('Time t (MC sweeps)', color='white')
ax4.set_ylabel('Var(Q)', color='white')
ax4.set_title('Fluctuation Scaling', color='#6ee7b7', fontsize=13)
ax4.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()

print(f"Current statistics: mean={Q_mean:.1f}, std={Q_std:.1f}")
print(f"Skewness: measured={measured_skew:.4f}, TW={tw_skew:.4f}")
print(f"Excess kurtosis: measured={measured_kurt:.4f}, TW={tw_kurt:.4f}")
print("Tracy-Widom analysis complete.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

ASEP current fluctuations scale as $t^{1/3}$(not $t^{1/2}$), following the Tracy-Widom GUE distribution.
Non-Gaussian signatures: skewness $\approx 0.224$, excess kurtosis $\approx 0.093$.
The KPZ universality class connects traffic flow to interface growth, random matrices, and directed polymers.
The quantum group $U_q(\mathfrak{sl}_2)$ with$q_{\text{deform}} = \sqrt{q/p}$ controls the full KPZ-to-diffusive crossover.
The hierarchy ASEP → TASEP/SSEP → ASEP+Langmuir provides increasingly realistic traffic models, each with exact results from integrable systems.

← Five-Phase Diagram HJB Equation →