Chapter 10.2: TASEP Phase Diagram

The Totally Asymmetric Simple Exclusion Process

The TASEP is perhaps the most important exactly solvable model in non-equilibrium statistical mechanics. Particles on a one-dimensional lattice hop to the right with rate 1, subject to exclusion: a hop is blocked if the target site is occupied. With open boundaries—particles injected at the left at rate $\alpha$ and extracted at the right at rate $\beta$—the model exhibits a genuine non-equilibrium phase transition with three distinct phases.

In urban terms, the TASEP is a single-lane road segment with an on-ramp (left boundary) and off-ramp (right boundary). The phase diagram maps directly onto the regimes observed in real traffic: free flow, congested flow, and maximum capacity.

10.2.1 Model Definition and Dynamics

The TASEP consists of $L$ sites on a one-dimensional lattice. Each site $i \in \{1, 2, \ldots, L\}$ is either occupied ($\tau_i = 1$) or empty ($\tau_i = 0$). The dynamics in continuous time are:

Bulk hopping: For$i = 1, \ldots, L-1$: if site$i$ is occupied and site$i+1$ is empty, the particle hops right at rate 1.
Left boundary (injection): If site 1 is empty, a particle is inserted at rate $\alpha$.
Right boundary (extraction): If site$L$ is occupied, the particle is removed at rate $\beta$.

Master Equation

Let $P(\{\tau\}, t)$ be the probability of configuration $\{\tau_1, \ldots, \tau_L\}$at time $t$. The master equation is:

$$\frac{dP}{dt} = \sum_{\{\tau'\}} \bigl[W(\{\tau'\} \to \{\tau\}) P(\{\tau'\}) - W(\{\tau\} \to \{\tau'\}) P(\{\tau\})\bigr]$$

The transition rates $W$ are: for bulk hops,$W = \tau_i(1-\tau_{i+1})$; for injection,$W = \alpha(1-\tau_1)$; for extraction,$W = \beta\tau_L$.

10.2.2 Mean-Field Theory: Three Phases

Taking the mean-field (product measure) approximation$\langle \tau_i \tau_{i+1} \rangle \approx \langle \tau_i \rangle \langle \tau_{i+1} \rangle$, we write $\rho_i = \langle \tau_i \rangle$. In steady state, the current $J$ is uniform across the lattice:

$$J = \rho_i(1 - \rho_{i+1}) = \text{const} \quad \forall\, i$$

For a spatially uniform bulk ($\rho_i = \rho$), we get $J = \rho(1-\rho)$. The boundary conditions constrain the allowed bulk density:

Low-Density (LD) Phase

When $\alpha < \beta$ and$\alpha < 1/2$, the left boundary is the bottleneck. The bulk density and current are controlled by the injection rate:

$$\rho_{\text{LD}} = \alpha, \qquad J_{\text{LD}} = \alpha(1-\alpha)$$

High-Density (HD) Phase

When $\beta < \alpha$ and$\beta < 1/2$, the right boundary is the bottleneck:

$$\rho_{\text{HD}} = 1-\beta, \qquad J_{\text{HD}} = \beta(1-\beta)$$

Maximal-Current (MC) Phase

When $\alpha > 1/2$ and$\beta > 1/2$, neither boundary is a bottleneck. The system self-organises to maximise throughput:

$$\rho_{\text{MC}} = \frac{1}{2}, \qquad J_{\text{MC}} = \frac{1}{4}$$

Phase Boundaries

The three phase boundaries are:

$$\text{LD/HD:} \quad \alpha = \beta \quad (\alpha, \beta < \tfrac{1}{2})$$

$$\text{LD/MC:} \quad \alpha = \tfrac{1}{2} \quad (\beta > \tfrac{1}{2})$$

$$\text{HD/MC:} \quad \beta = \tfrac{1}{2} \quad (\alpha > \tfrac{1}{2})$$

The LD/HD transition is first order (discontinuous density). The LD/MC and HD/MC transitions are second order (continuous density, diverging correlation length).

10.2.3 Exact Solution: Matrix Product Ansatz

The remarkable feature of the TASEP is that its steady-state distribution can be computed exactly using the Derrida-Evans-Hakim-Pasquier matrix product ansatz (1993). The steady-state probability of a configuration$\{\tau_1, \ldots, \tau_L\}$ is:

$$P(\tau_1, \ldots, \tau_L) = \frac{1}{Z_L} \langle W | \prod_{i=1}^{L} \bigl(\tau_i D + (1-\tau_i) E\bigr) | V \rangle$$

Here $D$ and$E$ are (infinite-dimensional) matrices, and$\langle W|$,$|V\rangle$ are boundary vectors. The key algebraic relations are:

$$DE - ED = D + E$$

$$\langle W | E = \frac{1}{\alpha} \langle W |, \qquad D | V \rangle = \frac{1}{\beta} | V \rangle$$

The normalisation (partition function) is:

$$Z_L = \langle W | (D + E)^L | V \rangle$$

The algebra $DE - ED = D + E$ ensures that the steady-state condition (probability current conservation) is automatically satisfied for every bond. The exact current is:

$$J = \frac{Z_{L-1}}{Z_L}$$

In the thermodynamic limit $L \to \infty$, the partition function can be evaluated by contour integral methods, confirming the mean-field phase diagram exactly—a rare case where mean-field gives the correct phase boundaries.

10.2.4 Domain Wall Random Walk

On the coexistence line $\alpha = \beta < 1/2$, a sharp interface (domain wall) separates a low-density region (left) from a high-density region (right). This domain wall performs a random walk on the lattice:

$$D_+ = J_{\text{LD}} = \alpha(1-\alpha), \qquad D_- = J_{\text{HD}} = \beta(1-\beta)$$

When $\alpha = \beta$, the drift vanishes and the wall performs an unbiased random walk. The wall position is uniformly distributed over$[0, L]$, leading to a linear density profile interpolating between $\rho_- = \alpha$ and$\rho_+ = 1-\beta$.

When $\alpha \neq \beta$, the wall is biased toward one boundary, creating the pure LD or HD phase. The domain wall picture provides an intuitive understanding of the phase transition: LD has the wall pinned at the right, HD has it pinned at the left, and coexistence has it delocalised.

10.2.5 Bethe Ansatz and KPZ Connection

The TASEP is related to the XXZ spin chain via the mapping$\tau_i \to S_i^z + 1/2$. The Bethe ansatz yields the exact spectrum of the generator (Markov matrix). The spectral gap—which determines the relaxation time—scales as:

$$\Delta \sim L^{-3/2} \quad \Longrightarrow \quad \tau_{\text{relax}} \sim L^{3/2}$$

The exponent $3/2$ is the hallmark of the KPZ universality class (Kardar-Parisi-Zhang). This connects the TASEP to a vast class of non-equilibrium phenomena: interface growth, random matrices, directed polymers, and last-passage percolation. The current fluctuations in the TASEP follow the Tracy-Widom distribution—we will explore this in Module 11.

For traffic, the $L^{3/2}$ scaling means that a road segment of length $L$ requires time $\sim L^{3/2}$ to reach steady state after a perturbation—significantly longer than the$L^2$ diffusive scaling one might naively expect.

10.2.6 Simulation: TASEP Phases and Phase Diagram

We simulate the TASEP using a random-sequential update algorithm and visualise the density profiles in each phase, along with the full phase diagram.

TASEP Simulation: Density Profiles in All Phases

Python

script.py88 lines

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

np.random.seed(42)

def simulate_tasep(L, alpha, beta, T_steps, T_measure):
    """Simulate TASEP with random-sequential update."""
    lattice = np.zeros(L, dtype=int)
    density_profile = np.zeros(L, dtype=float)
    measurements = 0

for t in range(T_steps):
        # Random-sequential update: pick a random site/event
        r = np.random.random()
        if r < alpha / (L + 2):
            # Left boundary injection
            if lattice[0] == 0:
                lattice[0] = 1
        elif r < (alpha + beta) / (L + 2):
            # Right boundary extraction
            if lattice[L-1] == 1:
                lattice[L-1] = 0
        else:
            # Bulk hop
            i = np.random.randint(0, L-1)
            if lattice[i] == 1 and lattice[i+1] == 0:
                lattice[i] = 0
                lattice[i+1] = 1

# Measure after equilibration
        if t > T_measure:
            density_profile += lattice
            measurements += 1

return density_profile / max(measurements, 1)

L = 200
T_steps = 2_000_000
T_measure = 500_000

# ── Three phases ──
params = [
    (0.2, 0.8, 'LD Phase (α=0.2, β=0.8)'),
    (0.8, 0.2, 'HD Phase (α=0.8, β=0.2)'),
    (0.8, 0.8, 'MC Phase (α=0.8, β=0.8)'),
    (0.3, 0.3, 'Coexistence (α=β=0.3)'),
]

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.patch.set_facecolor('#0a0a0a')
colors_phase = ['#34d399', '#f87171', '#60a5fa', '#fbbf24']

for idx, (alpha, beta, title) in enumerate(params):
    ax = axes[idx // 2, idx % 2]
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white')
    for spine in ax.spines.values():
        spine.set_color('#334155')

profile = simulate_tasep(L, alpha, beta, T_steps, T_measure)
    x = np.arange(L) / L

ax.plot(x, profile, color=colors_phase[idx], linewidth=1.5, alpha=0.8)

# Mean-field predictions
    if alpha < beta and alpha < 0.5:
        ax.axhline(y=alpha, color='white', linestyle='--', alpha=0.5, label=f'MF: ρ={alpha}')
    elif beta < alpha and beta < 0.5:
        ax.axhline(y=1-beta, color='white', linestyle='--', alpha=0.5, label=f'MF: ρ={1-beta}')
    elif alpha >= 0.5 and beta >= 0.5:
        ax.axhline(y=0.5, color='white', linestyle='--', alpha=0.5, label='MF: ρ=0.5')
    else:
        ax.axhline(y=alpha, color='#fbbf24', linestyle='--', alpha=0.4, label=f'α=β={alpha}')
        ax.axhline(y=1-beta, color='#f472b6', linestyle=':', alpha=0.4, label=f'1-β={1-beta}')

ax.set_xlabel('Position x/L', color='white')
    ax.set_ylabel('Density ρ', color='white')
    ax.set_title(title, color='#6ee7b7', fontsize=12)
    ax.set_ylim(-0.05, 1.05)
    ax.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("TASEP density profiles plotted for all three phases + coexistence.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

TASEP Phase Diagram

Python

script.py66 lines

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

# ── TASEP Phase Diagram ──
fig, ax = plt.subplots(figsize=(8, 8))
fig.patch.set_facecolor('#0a0a0a')
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
for spine in ax.spines.values():
    spine.set_color('#334155')

# Phase regions
alpha_range = np.linspace(0, 1, 500)
beta_range = np.linspace(0, 1, 500)
A, B = np.meshgrid(alpha_range, beta_range)

# Classify phases
phase = np.zeros_like(A)
# LD: alpha < beta, alpha < 0.5
ld_mask = (A < B) & (A < 0.5)
# HD: beta < alpha, beta < 0.5
hd_mask = (B < A) & (B < 0.5)
# MC: alpha > 0.5, beta > 0.5
mc_mask = (A >= 0.5) & (B >= 0.5)
# Coexistence line patches
coex_mask = (np.abs(A - B) < 0.02) & (A < 0.5)

phase[ld_mask] = 1
phase[hd_mask] = 2
phase[mc_mask] = 3

from matplotlib.colors import ListedColormap
cmap = ListedColormap(['#1a1a2e', '#065f46', '#7f1d1d', '#1e3a5f'])
ax.pcolormesh(A, B, phase, cmap=cmap, shading='auto', alpha=0.6)

# Phase boundaries
ax.plot([0, 0.5], [0, 0.5], color='#fbbf24', linewidth=3, label='LD/HD (1st order)')
ax.plot([0.5, 0.5], [0.5, 1.0], color='#34d399', linewidth=3, label='LD/MC (2nd order)')
ax.plot([0.5, 1.0], [0.5, 0.5], color='#f472b6', linewidth=3, label='HD/MC (2nd order)')
ax.plot([0.5], [0.5], 'o', color='white', markersize=10, zorder=5, label='Multicritical point')

# Labels
ax.text(0.15, 0.4, 'LD', fontsize=24, color='#34d399', fontweight='bold', ha='center')
ax.text(0.4, 0.15, 'HD', fontsize=24, color='#f87171', fontweight='bold', ha='center')
ax.text(0.75, 0.75, 'MC', fontsize=24, color='#60a5fa', fontweight='bold', ha='center')

# Current in each region
ax.text(0.15, 0.35, r'$J = \alpha(1-\alpha)$', fontsize=11, color='#6ee7b7', ha='center')
ax.text(0.4, 0.10, r'$J = \beta(1-\beta)$', fontsize=11, color='#fca5a5', ha='center')
ax.text(0.75, 0.70, r'$J = 1/4$', fontsize=11, color='#93c5fd', ha='center')

ax.set_xlabel('α (injection rate)', color='white', fontsize=14)
ax.set_ylabel('β (extraction rate)', color='white', fontsize=14)
ax.set_title('TASEP Phase Diagram', color='#6ee7b7', fontsize=16)
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_aspect('equal')
ax.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=10, loc='upper right')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("TASEP phase diagram plotted.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

The TASEP with open boundaries has three phases: LD ($\rho = \alpha$), HD ($\rho = 1-\beta$), and MC ($\rho = 1/2$).
The matrix product ansatz provides the exact steady-state distribution through the algebra $DE - ED = D + E$.
On the LD/HD coexistence line, a domain wall random walks along the lattice, giving a linear density profile.
The relaxation time scales as $\tau \sim L^{3/2}$, placing the TASEP in the KPZ universality class.
Mean-field gives exact phase boundaries but misses correlations visible in the density profile near boundaries.

← Langmuir Fundamental Diagram LWR & Spinodal Decomposition →