The KPZ Equation

The Kardar-Parisi-Zhang equation describes the stochastic growth of interfaces—from crystal surfaces to bacterial colonies. Here we apply it to the roughening boundary of an expanding city, where the interplay of smoothing, nonlinear lateral growth, and randomness produces universal scaling laws.

1. The KPZ Equation

Let $h(x,t)$ denote the height of an interface (or the radial distance of a city boundary from its centre) at lateral position $x$ and time $t$. The KPZ equation is:

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

Each term has a physical interpretation:

$\nu\,\nabla^2 h$ — Surface tension / smoothing. In the urban context, this represents planning regulations and infrastructure that smooth the city boundary.
$\frac{\lambda}{2}(\nabla h)^2$ — Nonlinear lateral growth. Slopes grow faster than flat regions because development “fills in” concavities. This is the term that makes KPZ distinct from the linear Edwards-Wilkinson equation.
$\eta(x,t)$ — White noise with $\langle\eta(x,t)\,\eta(x',t')\rangle = 2D\,\delta(x-x')\delta(t-t')$. Stochastic fluctuations from individual land-use decisions.

2. Roughness and Family-Vicsek Scaling

The central observable is the roughness (interface width) measured over a system of size $L$:

$$W(L,t) = \left\langle \overline{\bigl(h - \bar{h}\bigr)^2}\right\rangle^{1/2}$$

where the overbar is a spatial average and angle brackets denote ensemble averaging. The Family-Vicsek scaling ansatz states:

$$W(L,t) \sim L^{\alpha}\,f\!\left(\frac{t}{L^z}\right)$$

where the scaling function behaves as:

$$f(u) \sim \begin{cases} u^{\beta} & u \ll 1 \\ \text{const} & u \gg 1 \end{cases}$$

This defines three critical exponents:

$\alpha$ — Roughness exponent. Controls the saturated roughness: $W_{\text{sat}} \sim L^{\alpha}$.
$\beta$ — Growth exponent. Controls early-time roughening: $W \sim t^{\beta}$.
$z = \alpha/\beta$ — Dynamic exponent. Sets the crossover time $t_{\times} \sim L^z$.

For the KPZ equation in 1+1 dimensions, the exact exponents are:

$$\alpha = \frac{1}{2}, \qquad \beta = \frac{1}{3}, \qquad z = \frac{3}{2}$$

These are universal—they depend only on dimensionality and symmetry, not on the microscopic details. Any growth process in the KPZ universality class shares these exponents.

3. Why These Exponents?

The exact exponents follow from a remarkable connection. Via the Cole-Hopf transformation $Z = e^{\lambda h / 2\nu}$, the KPZ equation maps to the stochastic heat equation:

$$\frac{\partial Z}{\partial t} = \nu\,\nabla^2 Z + \frac{\lambda}{2\nu}\,\eta\,Z$$

which is linear in $Z$. The partition function $Z$ connects to directed polymers in random media, and the free energy $h \propto \ln Z$ inherits the Tracy-Widom distribution for its fluctuations.

A simple scaling argument also fixes the exponent relation. If we rescale $x \to bx$, $t \to b^z t$, $h \to b^{\alpha} h$, then each term in KPZ transforms as:

$$b^{\alpha - z} = b^{\alpha - 2} = b^{2\alpha - 2} = b^{-(1+z)/2}$$

The first equality (diffusion) gives $z = 2$ only for the linear Edwards-Wilkinson case. The nonlinear term enforces the additional constraint $\alpha + z = 2$, and the noise term gives $z = 2\alpha + 1$ in $d=1$. Solving simultaneously: $\alpha = 1/2,\; z = 3/2,\; \beta = 1/3$.

4. KPZ Simulation and Roughness Measurement

We discretise the KPZ equation on a 1D lattice and measure the roughness exponent $\beta$ from the early-time growth regime. The discretisation uses a simple Euler scheme with careful treatment of the nonlinear term.

KPZ Simulation with Roughness Exponent Measurement

Python

script.py111 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- KPZ parameters ---
L = 512          # system size
dx = 1.0
dt = 0.01
nu = 1.0         # diffusion (surface tension)
lam = 3.0        # nonlinearity strength
D_noise = 0.5    # noise amplitude
n_steps = 8000
measure_every = 20

h = np.zeros(L)  # flat initial condition

# Storage for roughness vs time
times = []
roughness = []

for step in range(n_steps):
    # Periodic boundary: gradients via central differences
    hp = np.roll(h, -1)  # h_{i+1}
    hm = np.roll(h, 1)   # h_{i-1}

laplacian = (hp - 2*h + hm) / dx**2
    grad = (hp - hm) / (2*dx)

# Euler step
    noise = np.sqrt(2*D_noise*dt/dx) * np.random.randn(L)
    h += dt * (nu * laplacian + 0.5*lam * grad**2) + noise

# Measure roughness
    if step % measure_every == 0:
        W = np.std(h)
        times.append(step * dt)
        roughness.append(W)

times = np.array(times)
roughness = np.array(roughness)

# --- Fit beta from early-time power law ---
# Use region where roughness is growing (before saturation)
mask = (times > 0.5) & (times < 20)
log_t = np.log(times[mask])
log_W = np.log(roughness[mask])
beta_fit, intercept = np.polyfit(log_t, log_W, 1)

# --- Plots ---
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))

# (a) Interface snapshots
ax = axes[0]
for snap_step in [0, 2000, 5000, 7999]:
    # Re-run to get snapshot (simplified: use final state shifted)
    pass
ax.set_title('(see roughness plot)')

# Rerun to collect snapshots
h_snap = np.zeros(L)
snapshots = {0: h_snap.copy()}
snap_times_wanted = [2000, 5000, 7999]
for step in range(n_steps):
    hp = np.roll(h_snap, -1)
    hm = np.roll(h_snap, 1)
    laplacian = (hp - 2*h_snap + hm) / dx**2
    grad = (hp - hm) / (2*dx)
    noise = np.sqrt(2*D_noise*dt/dx) * np.random.randn(L)
    h_snap += dt * (nu * laplacian + 0.5*lam * grad**2) + noise
    if step in snap_times_wanted:
        snapshots[step] = h_snap.copy()

ax = axes[0]
colors = ['#065f46', '#059669', '#34d399', '#a7f3d0']
for idx, (s, hh) in enumerate(sorted(snapshots.items())):
    ax.plot(hh[:200], color=colors[idx], alpha=0.8,
            linewidth=0.7, label=f't={s*dt:.0f}')
ax.set_xlabel('x')
ax.set_ylabel('h(x,t)')
ax.set_title('Interface Snapshots')
ax.legend(fontsize=8)

# (b) Roughness vs time (log-log)
ax = axes[1]
ax.loglog(times, roughness, 'o', color='teal', markersize=2, alpha=0.6)
t_fit = np.linspace(0.5, 20, 100)
ax.loglog(t_fit, np.exp(intercept)*t_fit**beta_fit, '--', color='crimson',
          linewidth=2, label=f'fit: beta={beta_fit:.3f}')
ax.axhline(roughness[-1], color='gray', linestyle=':', alpha=0.5, label='saturation')
ax.set_xlabel('t')
ax.set_ylabel('W(t)')
ax.set_title('Roughness Growth')
ax.legend(fontsize=9)

# (c) Local slope distribution (height differences)
ax = axes[2]
dh = np.diff(h)
ax.hist(dh, bins=60, density=True, color='teal', alpha=0.7, edgecolor='black', linewidth=0.3)
ax.set_xlabel('dh/dx')
ax.set_ylabel('Density')
ax.set_title('Slope Distribution')

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight')
plt.show()

print(f"Measured growth exponent: beta = {beta_fit:.4f}")
print(f"KPZ exact (1+1d):        beta = 0.3333")
print(f"Saturated roughness:      W_sat = {roughness[-1]:.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. City Boundary Roughening

Map the KPZ variables to urban sprawl:

$h(x,t)$ is the radial extent of the built-up area at angular position $x$ around the city centre.
The $\nu\,\nabla^2 h$ term represents planning regulations that smooth the urban boundary (green belts, zoning).
The $(\nabla h)^2$ nonlinearity captures the fact that concavities in the boundary fill in faster than convexities grow out—infill development is cheaper than greenfield expansion.
The noise $\eta$ reflects the stochastic nature of individual building decisions.

Empirical studies of cities like London and Delhi show boundary roughness exponents consistent with KPZ scaling. The universality means the precise zoning rules do not matter—only the symmetry class determines the large-scale behaviour of the boundary.

Key Insight

KPZ universality tells us that city boundaries roughen according to power laws that are insensitive to microscopic planning details. The exponents $\alpha = 1/2,\; \beta = 1/3$ are as fundamental to urban growth interfaces as they are to burning paper or bacterial colonies.

← Bayesian MCMC Fortran Parameter Search →