Jamming: Where a Packing First Becomes Rigid

Pour sand into a jar and shake it: at some density the grains lock and the pile holds a shovel. Nothing melted or froze — the grains just ran out of room to move, and a fluid of nothing became a solid that supports load. That is jamming, and its cleanest laboratory is the simplest one imaginable: soft disks that pay an energy penalty only when they overlap. Compress a random arrangement of them and there is a sharp packing fraction ${\phi }_J$ where the system first cannot avoid overlaps — where the nearest energy minimum stops being zero. This page runs the O'Hern–Silbert–Liu–Nagel protocol (2003) to pin that point, confirms the rigidity condition that makes it a genuine transition, and measures how stiffness grows above it — where the model teaches a lesson about finite-size systematic error that the raw exponent alone would hide.

@njit
def energy_forces(pos, radii, L):
    """Harmonic soft disks: V = (1/2)(1 - r/sigma)^2 for r < sigma, else 0."""
    F = np.zeros((N, 2)); E = 0.0; contacts = np.zeros(N)
    for i in range(N-1):
        for j in range(i+1, N):
            d = min_image(pos[i] - pos[j], L)      # periodic box
            r = norm(d); sij = radii[i] + radii[j]
            if r < sij:                            # only overlaps cost energy
                overlap = 1 - r/sij
                E += 0.5 * overlap**2
                fpair = overlap / (sij * r) * d
                F[i] += fpair; F[j] -= fpair
                contacts[i] += 1; contacts[j] += 1
    return E, F, contacts

@njit
def fire_minimize(pos, radii, L):
    """FIRE: run MD, but whenever the velocity points against the force,
    kill it — descend to the nearest energy minimum without solving a Hessian."""
    v = zeros_like(pos); dt = 0.02; alpha = 0.1
    E, F, _ = energy_forces(pos, radii, L)
    while max_force(F) > 1e-11:
        if dot(F, v) > 0:                          # going downhill: accelerate
            v = (1-alpha)*v + alpha*norm(v)/norm(F)*F
            dt = min(dt*1.1, dt_max); alpha *= 0.99
        else:                                      # fighting the force: freeze
            v[:] = 0; dt *= 0.5; alpha = 0.1
        v, pos = velocity_verlet(pos, v, F, dt)    # unit-mass Verlet step
        E, F, _ = energy_forces(pos, radii, L)
    return E / N

The onset of rigidity is an energy minimum turning on

The order parameter for jamming is not a density or a magnetization — it is whether the potential energy of the relaxed packing is zero. Below ${\phi }_J$ a random configuration can always be squeezed into a nearby arrangement with no overlaps at all: energy zero, no restoring forces, a floppy assembly of disks that happen to be near each other. Above ${\phi }_J$ there is not enough room; every minimum has some residual overlap, some stored elastic energy, and a network of contacts that pushes back. So the protocol is: take a random (infinite-temperature) configuration, drop it straight down to its nearest energy minimum, and ask whether that minimum is zero. FIRE does the dropping — it runs molecular dynamics but kills the velocity the instant it starts fighting the force, so it slides downhill to a local minimum without ever building or inverting a Hessian. Bisect the packing fraction between a floppy sample and a rigid one and you converge on that sample's threshold.

N = 512 bidisperse disks (ratio 1.4), 40 samples
phi_J = 0.8441 +/- 0.0002   (O'Hern 2D bidisperse ~ 0.84; 40 samples)
dphi = 3e-05:  Z = 4.0057 +/- 0.0006   Z - 4 = 0.0057   e/N = 1.98e-10
dphi = 1e-04:  Z = 4.0177 +/- 0.0012   Z - 4 = 0.0177   e/N = 2.14e-09
dphi = 3e-04:  Z = 4.0411 +/- 0.0015   Z - 4 = 0.0411   e/N = 1.91e-08
dphi = 1e-03:  Z = 4.0870 +/- 0.0016   Z - 4 = 0.0870   e/N = 2.13e-07
dphi = 3e-03:  Z = 4.1676 +/- 0.0024   Z - 4 = 0.1676   e/N = 1.96e-06
dphi = 1e-02:  Z = 4.3429 +/- 0.0028   Z - 4 = 0.3429   e/N = 2.27e-05
dphi = 3e-02:  Z = 4.6330 +/- 0.0033   Z - 4 = 0.6330   e/N = 2.16e-04

fit (all offsets):         exponent 0.667
fit (dphi >= 3e-3 window):  exponent 0.577   (O'Hern exponent: 0.5)
local log-slopes (near -> far from jamming):
   3e-05 -> 1e-04:  0.946
   1e-04 -> 3e-04:  0.767
   3e-04 -> 1e-03:  0.622
   1e-03 -> 3e-03:  0.597
   3e-03 -> 1e-02:  0.595
   1e-02 -> 3e-02:  0.558

Over 40 independent samples the threshold lands at ${\phi }_J=0.8441\pm 0.0002$ — squarely on the accepted value near 0.84 for 50:50 bidisperse disks at size ratio 1.4 (the ratio is there to frustrate crystallization; equal disks would rather form a triangular lattice than jam). The spread across samples is real physics, not noise: a finite random packing has no single jamming point, it has a distribution of them that narrows as $N\to \infty$. The energy per particle in the last column vanishing as the offset shrinks — $2\times 10^{-10}$ at the closest approach — is the check that we are sitting just barely inside the jammed phase at each measurement, not deep in it.

Isostaticity: exactly enough contacts, and not one more

Why is jamming sharp? Because of a counting argument that goes back to Maxwell. A packing of $N$ disks in $d$ dimensions has $dN$ positional degrees of freedom to pin down. Each contact between two disks is one constraint (their separation is fixed at touching). Rigidity needs at least as many constraints as freedoms; stability against collapse forbids many more than that. At the jamming point the two bounds meet — the packing is isostatic, carrying exactly enough contacts to be rigid and no redundant ones:

The factor of two is because each contact is shared by two disks. The measurement confirms it cleanly: extrapolated to the threshold the contact number $Z$ runs straight to 4 (it reads 4.006 at the closest offset and falls toward 4 as the offset shrinks). One subtlety makes or breaks the number — rattlers. A packing always traps a few disks with fewer than $d+1=3$ contacts: caged loosely but not actually load-bearing, free to rattle in their holes. They carry no force and must be removed before counting, iteratively, because removing one can strand another. Count them in and $Z$ comes out below 4 and the whole isostatic story dissolves. The rigidity is in the backbone, and finding the backbone is half the measurement.

Stiffness above jamming — and where the exponent hides

Push past ${\phi }_J$ and the packing stiffens: it grows extra contacts beyond the isostatic minimum, and the excess is the signature power law of the jamming transition,

A square root — the excess coordination grows infinitely steeply at the transition, the fingerprint that made jamming look like a critical point. Now the part worth the page. Fit all seven offsets and the exponent comes out 0.667, not 0.5 — off by a third. The temptation is to fix it by measuring closer to ${\phi }_J$, where the asymptotic law "should" be cleanest. That is exactly backwards, and the local slopes say so: reading them near-to-far from jamming they run 0.95, 0.77, 0.62, 0.60, 0.60, 0.56 — steepest closest in, drifting toward 0.5 as you back away. The closest offset of all, at $3\times 10^{-5}$, has a local slope of 0.95 — nearly double the true exponent. The smallest offsets are the most contaminated, not the least.

The reason is finite size. In a box of $N=512$, the jamming transition is rounded: below a crossover $\Delta {\phi }^{\ast }\sim 1/N$ the excess coordination hits a floor set by the handful of surface contacts, and the clean power law is buried under that floor. The asymptotic $1/2$ law lives at intermediate offsets — far enough out that finite-$N$ rounding is negligible, not so far that the packing has left the critical regime. Fit only the window $\Delta \phi \ge 3\times 10^{-3}$ and the exponent drops to 0.577, with the local slope at the largest offsets down to 0.56 and still falling toward the true value the thermodynamic limit would deliver. The dotted $1/2$ reference on the plot runs right through the far data and pulls away from the near data — a picture of the rounding. The lesson generalizes past jamming: when a finite system rounds a transition, the asymptotic exponent is not at the point, it is in the window you can see before finite size takes over — and a fit that swallows the rounded region reports a number that is an artifact of the box, precise and wrong.

Extensions

• Finite-size scaling done right — run $N$ from 256 to 8192, collapse $Z-Z_{iso}$ against $N(\phi -{\phi }_J)$, and read the exponent from the data collapse rather than a single-$N$ fit — the same machinery as the Ising finite-size scaling flagship.
• The soft modes — an isostatic packing sits at the edge of stability, so its vibrational spectrum has an anomalous excess of low-frequency modes (the "boson peak"): diagonalize the dynamical matrix and watch the density of states pile up at $\omega \to 0$ as $\phi \to {\phi }_J$.
• Pressure and bulk modulus — both vanish as power laws above ${\phi }_J$ too; the contact-value machinery from the hard-sphere fluid connects the overlaps here to the pressure there.
• From jamming to glass — jamming is the $T=0$ endpoint; turn temperature back on and the same disks fall out of equilibrium into a glass. The hard-sphere fluid compressed fast enough never reaches its crystal and jams instead — the granular and the glassy are two faces of the same rigidity.

Reproduce it

scripts/gen_jamming.py runs the full protocol — numba-compiled FIRE, per-sample bisection, iterative rattler removal — in a couple of minutes. Its history is the page's own lesson: an earlier version reported the exponent as high as 0.75 by fitting all offsets and chasing the number with tighter bisection tolerances, treating a finite-size artifact as a convergence problem. The fix was not a better solver but a better window — the realization that the box, not the algorithm, sets how close to ${\phi }_J$ the asymptotic law is still visible.