Chapter 11.1: General ASEP with Langmuir Kinetics

Beyond TASEP: Two-Way Hopping and Bulk Reactions

The TASEP is a special case of the Asymmetric Simple Exclusion Process (ASEP), which allows hopping in both directions. Moreover, on real roads, vehicles enter and exit not just at boundaries but also along the road (on-ramps, off-ramps, parking). This bulk attachment/detachment is precisely the Langmuir kinetics from Chapter 10.1. Coupling ASEP with Langmuir yields a rich model that interpolates between transport-dominated and reaction-dominated regimes.

11.1.1 The Full ASEP Master Equation

In the general ASEP, particles hop to the right with rate$p$ and to the left with rate$q$, subject to exclusion. The bulk update rules are:

$$10 \xrightarrow{p} 01, \qquad 01 \xrightarrow{q} 10$$

The mean-field current in the bulk is:

$$J = (p - q)\,\rho(1-\rho)$$

The effective drift velocity is $v_{\text{eff}} = p - q$. When $q = 0$, this reduces to the TASEP. When $p = q$, we get the Symmetric Simple Exclusion Process (SSEP) with zero net current.

General Boundary Conditions

The most general open boundary conditions involve four parameters:

α: Injection rate at the left boundary ($0 \to 1$ at site 1)
γ: Extraction rate at the left boundary ($1 \to 0$ at site 1)
δ: Injection rate at the right boundary ($0 \to 1$ at site L)
β: Extraction rate at the right boundary ($1 \to 0$ at site L)

In traffic terms: $\alpha$ is the main on-ramp entrance rate, $\beta$ is the main off-ramp exit rate, $\delta$ represents vehicles entering from a merging road at the end, and $\gamma$ represents vehicles making U-turns near the beginning.

The boundary densities imposed by the left and right reservoirs are:

$$\rho_L = \frac{\alpha}{\alpha + \gamma}, \qquad \rho_R = \frac{\delta}{\delta + \beta}$$

11.1.2 Langmuir Coupling: Bulk Attachment and Detachment

In addition to boundary-driven transport, we now allow particles to attach to and detach from the lattice at every site in the bulk:

$$0 \xrightarrow{\omega_A} 1 \quad \text{(attachment)}, \qquad 1 \xrightarrow{\omega_D} 0 \quad \text{(detachment)}$$

These are exactly the Langmuir adsorption/desorption processes from Chapter 10.1. In traffic, they represent vehicles entering from side streets (attachment) or exiting to side streets (detachment).

The Competition Parameter

The crucial insight is the scaling of the Langmuir rates with system size. For a lattice of $L$ sites, the total Langmuir contribution scales as $\omega L$(summed over all bulk sites), while the boundary contribution is$O(1)$. For both mechanisms to compete, we set:

$$\omega_A = \frac{\Omega_A}{L}, \qquad \omega_D = \frac{\Omega_D}{L}$$

The dimensionless competition parameter $\Omega = (\Omega_A + \Omega_D)$ controls the relative importance of bulk kinetics versus boundary-driven transport:

$\Omega \ll 1$: Transport-dominated — boundary conditions control the density profile (TASEP-like behavior)
$\Omega \gg 1$: Reaction-dominated — Langmuir kinetics control the bulk density toward the equilibrium$\rho_{\text{eq}} = \Omega_A / (\Omega_A + \Omega_D)$
$\Omega \sim 1$: Both mechanisms compete — new phases emerge

Mean-Field Equation

The continuum mean-field equation for the ASEP+Langmuir system is:

$$\frac{\partial \rho}{\partial t} + (p-q)\frac{\partial}{\partial x}\bigl[\rho(1-\rho)\bigr] = \Omega_A(1-\rho) - \Omega_D \rho$$

The left side is the ASEP conservation law; the right side is the Langmuir source/sink. In steady state, the density profile satisfies the ODE:

$$(p-q)(1-2\rho)\frac{d\rho}{dx} = \Omega_A(1-\rho) - \Omega_D \rho$$

11.1.3 The Langmuir Equilibrium Density

Setting the Langmuir source term to zero gives the bulk equilibrium density:

$$\rho_{\text{eq}} = \frac{\Omega_A}{\Omega_A + \Omega_D} = \frac{K}{1 + K}$$

where $K = \Omega_A / \Omega_D$ is the Langmuir equilibrium constant. This is a fixed point of the bulk dynamics, attracting the local density. When transport is weak, the entire bulk converges to$\rho_{\text{eq}}$ regardless of boundary conditions.

The competition between boundary-imposed densities and the Langmuir equilibrium density creates the rich phase structure we will explore in the next chapter. The key question is: what happens when $\rho_L \neq \rho_{\text{eq}} \neq \rho_R$? The answer is a five-phase diagram with domain wall localisation.

11.1.4 Simulation: ASEP + Langmuir

We simulate the general ASEP with Langmuir kinetics using a continuous-time Monte Carlo algorithm and compare transport-dominated versus reaction-dominated regimes.

ASEP + Langmuir: Transport vs Reaction Regimes

Python

script.py152 lines

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

np.random.seed(42)

def simulate_asep_langmuir(L, p_hop, q_hop, alpha, beta, omega_A, omega_D, T_steps, T_eq):
    """Simulate general ASEP + Langmuir with random-sequential update."""
    lattice = np.zeros(L, dtype=int)
    density = np.zeros(L, dtype=float)
    measurements = 0
    total_current = 0.0

# Total rate for event selection
    for t in range(T_steps):
        # Select random event
        site = np.random.randint(0, L)

if site == 0:
            # Left boundary
            r = np.random.random()
            if lattice[0] == 0 and r < alpha:
                lattice[0] = 1
            elif lattice[0] == 1 and site < L-1 and lattice[1] == 0:
                if np.random.random() < p_hop:
                    lattice[0] = 0
                    lattice[1] = 1
        elif site == L - 1:
            # Right boundary
            r = np.random.random()
            if lattice[L-1] == 1 and r < beta:
                lattice[L-1] = 0
            elif lattice[L-1] == 1 and lattice[L-2] == 0:
                if np.random.random() < q_hop:
                    lattice[L-1] = 0
                    lattice[L-2] = 1
        else:
            # Bulk
            r = np.random.random()
            if lattice[site] == 1 and site < L-1 and lattice[site+1] == 0:
                if r < p_hop:
                    lattice[site] = 0
                    lattice[site+1] = 1
            elif lattice[site] == 1 and site > 0 and lattice[site-1] == 0:
                if r < q_hop:
                    lattice[site] = 0
                    lattice[site-1] = 1

# Langmuir attachment/detachment
            r2 = np.random.random()
            if lattice[site] == 0 and r2 < omega_A:
                lattice[site] = 1
            elif lattice[site] == 1 and r2 < omega_D:
                lattice[site] = 0

if t > T_eq:
            density += lattice
            measurements += 1

return density / max(measurements, 1)

L = 100
T_steps = 500_000
T_eq = 100_000

# Parameters
p_hop = 1.0
q_hop = 0.0  # TASEP limit
alpha = 0.2
beta = 0.8

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

# ── Panel 1: Pure TASEP (no Langmuir) ──
ax = axes[0, 0]
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
for spine in ax.spines.values():
    spine.set_color('#334155')
profile = simulate_asep_langmuir(L, p_hop, q_hop, alpha, beta, 0, 0, T_steps, T_eq)
x = np.arange(L) / L
ax.plot(x, profile, color='#34d399', linewidth=1.5)
ax.axhline(y=alpha, color='#fbbf24', linestyle='--', alpha=0.6, label=f'ρ_LD = {alpha}')
ax.set_xlabel('Position x/L', color='white')
ax.set_ylabel('Density ρ', color='white')
ax.set_title('Pure TASEP (Ω=0)', color='#6ee7b7', fontsize=12)
ax.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)
ax.set_ylim(-0.05, 1.05)

# ── Panel 2: Weak Langmuir ──
ax = axes[0, 1]
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
for spine in ax.spines.values():
    spine.set_color('#334155')
omega_A_weak = 0.5 / L
omega_D_weak = 0.5 / L
rho_eq = omega_A_weak / (omega_A_weak + omega_D_weak)
profile = simulate_asep_langmuir(L, p_hop, q_hop, alpha, beta, omega_A_weak, omega_D_weak, T_steps, T_eq)
ax.plot(x, profile, color='#2dd4bf', linewidth=1.5)
ax.axhline(y=alpha, color='#fbbf24', linestyle='--', alpha=0.6, label=f'ρ_LD = {alpha}')
ax.axhline(y=rho_eq, color='#f472b6', linestyle=':', alpha=0.6, label=f'ρ_eq = {rho_eq:.2f}')
ax.set_xlabel('Position x/L', color='white')
ax.set_ylabel('Density ρ', color='white')
ax.set_title('Weak Langmuir (Ω=1)', color='#6ee7b7', fontsize=12)
ax.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)
ax.set_ylim(-0.05, 1.05)

# ── Panel 3: Strong Langmuir ──
ax = axes[1, 0]
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
for spine in ax.spines.values():
    spine.set_color('#334155')
omega_A_str = 5.0 / L
omega_D_str = 3.0 / L
rho_eq2 = omega_A_str / (omega_A_str + omega_D_str)
profile = simulate_asep_langmuir(L, p_hop, q_hop, alpha, beta, omega_A_str, omega_D_str, T_steps, T_eq)
ax.plot(x, profile, color='#a78bfa', linewidth=1.5)
ax.axhline(y=rho_eq2, color='#f472b6', linestyle=':', alpha=0.6, label=f'ρ_eq = {rho_eq2:.2f}')
ax.set_xlabel('Position x/L', color='white')
ax.set_ylabel('Density ρ', color='white')
ax.set_title(f'Strong Langmuir (Ω={5+3})', color='#6ee7b7', fontsize=12)
ax.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)
ax.set_ylim(-0.05, 1.05)

# ── Panel 4: Asymmetric Langmuir (ρ_eq ≠ 1/2) ──
ax = axes[1, 1]
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
for spine in ax.spines.values():
    spine.set_color('#334155')
omega_A_asym = 6.0 / L
omega_D_asym = 2.0 / L
rho_eq3 = omega_A_asym / (omega_A_asym + omega_D_asym)
profile = simulate_asep_langmuir(L, p_hop, q_hop, 0.3, 0.3, omega_A_asym, omega_D_asym, T_steps, T_eq)
ax.plot(x, profile, color='#f87171', linewidth=1.5)
ax.axhline(y=rho_eq3, color='#f472b6', linestyle=':', alpha=0.6, label=f'ρ_eq = {rho_eq3:.2f}')
ax.axhline(y=0.3, color='#fbbf24', linestyle='--', alpha=0.6, label='α = β = 0.3')
ax.set_xlabel('Position x/L', color='white')
ax.set_ylabel('Density ρ', color='white')
ax.set_title('Asymmetric Langmuir (ρ_eq=0.75)', color='#6ee7b7', fontsize=12)
ax.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)
ax.set_ylim(-0.05, 1.05)

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("General ASEP + Langmuir simulation complete.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

11.1.5 Transport-Dominated vs Reaction-Dominated

The steady-state ODE can be solved by separation of variables. Defining$\xi = x/L$ (rescaled position on$[0,1]$):

$$(1-2\rho)\frac{d\rho}{d\xi} = \Omega_A(1-\rho) - \Omega_D \rho$$

This has a fixed point at $\rho = \rho_{\text{eq}}$. The flow direction depends on whether$\rho > \rho_{\text{eq}}$ or$\rho < \rho_{\text{eq}}$:

When $\rho < 1/2$: The factor $(1-2\rho) > 0$, so$d\rho/d\xi > 0$ when$\rho < \rho_{\text{eq}}$ (density increases rightward toward equilibrium).

When $\rho > 1/2$: The factor $(1-2\rho) < 0$, so the flow reverses. This creates the possibility of matching LD and HD solutions through a shock.

Mean-Field ODE: Transport vs Reaction Regimes

Python

script.py60 lines

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
from scipy.integrate import solve_ivp

# ── Solve the mean-field ODE for different Omega values ──
def rhs(xi, rho, Omega_A, Omega_D):
    denom = 1 - 2 * rho
    if abs(denom) < 1e-10:
        return 0.0
    return (Omega_A * (1 - rho) - Omega_D * rho) / denom

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.patch.set_facecolor('#0a0a0a')

# Different regimes
regimes = [
    (0.1, 0.1, 'Transport-dominated (Ω=0.2)'),
    (1.0, 1.0, 'Crossover (Ω=2.0)'),
    (5.0, 3.0, 'Reaction-dominated (Ω=8.0)'),
]

alpha_vals = [0.2, 0.3, 0.15]
colors_line = ['#34d399', '#2dd4bf', '#fbbf24', '#f87171']

for idx, (Omega_A, Omega_D, title) in enumerate(regimes):
    ax = axes[idx]
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white')
    for spine in ax.spines.values():
        spine.set_color('#334155')

rho_eq = Omega_A / (Omega_A + Omega_D)

for j, alpha in enumerate(alpha_vals):
        xi_span = [0, 1]
        xi_eval = np.linspace(0, 1, 500)
        sol = solve_ivp(rhs, xi_span, [alpha], t_eval=xi_eval,
                       args=(Omega_A, Omega_D), method='RK45',
                       max_step=0.001)
        if sol.success:
            rho_sol = np.clip(sol.y[0], 0, 1)
            ax.plot(sol.t, rho_sol, color=colors_line[j], linewidth=2,
                   label=f'ρ(0)={alpha}')

ax.axhline(y=rho_eq, color='#f472b6', linestyle=':', linewidth=2, alpha=0.7,
              label=f'ρ_eq={rho_eq:.2f}')
    ax.axhline(y=0.5, color='white', linestyle='--', alpha=0.2)
    ax.set_xlabel('Position ξ = x/L', color='white')
    ax.set_ylabel('Density ρ', color='white')
    ax.set_title(title, color='#6ee7b7', fontsize=11)
    ax.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=8)
    ax.set_ylim(-0.05, 1.05)

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Transport vs reaction regime comparison plotted.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

The general ASEP allows bidirectional hopping with rates$p$ (forward) and$q$ (backward), giving current$J = (p-q)\rho(1-\rho)$.
Langmuir kinetics (attachment $\omega_A$, detachment $\omega_D$) model bulk on/off-ramps.
The scaling $\omega \sim 1/L$ ensures competition between boundary transport and bulk reactions.
The parameter $\Omega = (\Omega_A + \Omega_D)L$ controls the crossover from TASEP-like to Langmuir-like behavior.
The Langmuir equilibrium density $\rho_{\text{eq}} = \Omega_A/(\Omega_A + \Omega_D)$ acts as a bulk attractor, competing with boundary-imposed densities.

← LWR & Spinodal Decomposition Five-Phase Diagram →