SWIRL

What Is Swirl?

Swirl simulates the Belousov-Zhabotinsky (BZ) chemical reaction — one of the most striking examples of spontaneous self-organization in chemistry. In a shallow dish, a carefully chosen mixture of malonic acid, bromate, and cerium catalyst oscillates rhythmically between oxidized and reduced states, and when left undisturbed forms rotating spiral waves that sweep outward at several millimetres per minute. No computer, no external input — just chemistry, far from equilibrium, organizing itself.

The simulation uses the Barkley model (1991), a minimal two-variable reaction-diffusion system that captures the essential excitable-media dynamics of the BZ reaction. A 200×200 periodic grid evolves under partial differential equations; eight numerical steps run per animation frame. The amber wavefronts are the oxidized activator; the dark teal wake is the refractory inhibitor that prevents a region from re-exciting too soon, which is why the spiral rotates rather than collapsing or exploding.

The Belousov-Zhabotinsky Reaction

Boris Belousov discovered the oscillating reaction in 1951 while trying to model the Krebs cycle biochemically. He found that a mixture of citric acid, cerium sulfate, and potassium bromate in sulfuric acid would spontaneously alternate between yellow (oxidized Ce⁴⁺) and colourless (reduced Ce³⁺) every few minutes. His manuscript was rejected twice by Soviet journals — reviewers refused to believe a chemical reaction could oscillate, since thermodynamics seemed to forbid periodic behaviour in a single-phase homogeneous system.

He was right. Anatol Zhabotinsky extended the work in the 1960s, confirmed the oscillation mechanism, and — crucially — observed the spiral waves that form when the reaction is run in a thin layer rather than a well-stirred beaker. Ilya Prigogine won the 1977 Nobel Prize in Chemistry in part for explaining how such dissipative structures — ordered patterns sustained by a continuous throughput of energy and matter — are thermodynamically consistent. They are not exceptions to the second law; they are the second law doing something spectacular.

The Barkley Model

The Barkley model (Dwight Barkley, 1991) distils excitable-media dynamics into two coupled partial differential equations. Despite having only four parameters, it reproduces the full phenomenology of the BZ reaction: plane waves, ring waves, spiral waves, and spatiotemporal chaos (spiral breakup).

∂u/∂t = (1/ε) · u · (1−u) · [u − (v + b)/a]  +  Du · ∇²u
∂v/∂t = u − v

u  — activator (HBrO₂, oxidized state) — diffuses rapidly, drives wavefronts
v  — inhibitor (catalyst Ce⁴⁺)         — evolves slowly, suppresses re-excitation

a = 0.75   — excitability threshold (free end of wavefront)
b = 0.06   — rest-state offset (sets the refractory threshold)
ε = 0.020  — time-scale separation (small ε → fast fronts, narrow arms)
Du = 0.01  — normalised diffusion constant for u
Dv ≈ 0     — v doesn't diffuse appreciably (ionic, larger molecule)

The kinetic term (1/ε)·u·(1−u)·(u−(v+b)/a) is a cubic that creates bistability between the resting state (u=0) and the excited state (u=1), gated by the inhibitor level v. When v is low, any small perturbation above the threshold (v+b)/a drives u rapidly to 1 — this is the wavefront. Behind the front, v rises (since dv/dt = u), which drives u back to 0 and makes the region refractory: unable to re-excite until v decays again.

Why Spirals?

A free end on a wavefront is unstable. In an excitable medium, the tip of a broken wavefront curves: the newly exposed resting tissue ahead of the tip begins to be excited by diffusion, while the refractory tissue behind the tip cannot follow. The net effect is that the tip rotates around the core — a small region of tissue that cycles periodically but never fully excites.

Once established, a spiral wave is a topologically stable object. Its rotation period T is set by the kinetics (inversely related to ε); its wavelength λ and the speed c of the arms are related by c = λ/T. The core position drifts slowly if the medium is heterogeneous. Two counter-rotating spirals annihilate when their arms collide; two co-rotating spirals maintain distance. This is the dynamics on display every time a BZ dish is left undisturbed in a lab.

Spiral period:    T ≈ 40–120 simulation time units  (depends on ε)
Wavefront speed:  c ≈ Du · (kinematics) — set by diffusion + excitability
Core radius:      r_core → 0 as ε → 0  (tip tightens to a point)

At small ε: thin, fast, tightly coiled spirals (high excitability)
At large ε: wide, slow, loosely wound spirals (reduced excitability)
At ε > ε_c ≈ a/b: spirals break up into spatiotemporal chaos
             (Winfree's "cardiac fibrillation" regime)

Where This Appears in Nature

Cardiac tissue        — re-entrant spiral waves in the heart are the
                         mechanism of ventricular fibrillation. The spiral
                         tip drifts, fragments, and transitions to chaos —
                         the same spiral breakup visible when ε is pushed high.

Neural cortex         — spreading cortical depression (migraines with aura)
                         propagates as a slow BZ-like wave at ~3 mm/min,
                         leaving a refractory wake. The visual aura is the front.

cAMP signalling       — Dictyostelium discoideum (slime mould) aggregates
                         using cyclic-AMP spiral waves to coordinate thousands
                         of cells into a single fruiting body.

CO oxidation          — platinum catalysis in a gas-phase reaction
                         produces spiral patterns at the nanoscale,
                         observable by photoemission electron microscopy.

Zebra fish pigment    — some skin pattern formation follows a BZ-like
                         activator-inhibitor mechanism (Turing + excitable).

Galaxy arms           — density waves in spiral galaxies share mathematical
                         structure with BZ spiral arms (though the physics differs).

Implementation

The simulation uses explicit Euler integration on a 200×200 toroidal grid with periodic boundary conditions. Two Float32Arrays hold the current u and v fields; two more hold the next time step. After computing the update, the buffers swap — no allocation per frame.

Grid:         200×200 periodic (toroidal), Float32Array × 4
Timestep:     dt = 0.01  (stable: Du · dt / dx² = 0.0001 ≪ 0.5)
Substeps:     8 per animation frame → 0.08 time units per frame
Display:      200×200 → 480×480 CSS (image-rendering: pixelated)

Inner loop per cell (8 × 40 000 = 320 000 cells per frame):
  lap(u) = u[↑] + u[↓] + u[←] + u[→] − 4·u   (5-point stencil, dx=1)
  f = (1/ε) · u · (1−u) · (u − (v+b)/a)
  u_new = clamp(u + dt · (Du · lap + f),  0, 1)
  v_new = clamp(v + dt · (u − v),         0, 1)

Rendering (per pixel, after each frame):
  ex = u                       — activator brightness
  rf = v · (1 − u)             — refractory signature (v without wavefront)
  R  = min(255, ex·235 + rf·14 + 10)   → amber wavefront
  G  = min(255, ex·175 + rf·36 +  5)   → gold
  B  = min(255, ex·60  + rf·96 +  4)   → teal refractory

Spiral seeding (initial condition):
  Three half-disk stimuli: excited hemisphere adjacent to refractory hemisphere.
  The free wavefront tip curls around the inhibitor boundary → stable rotor.
  Click to place additional ring waves (radius-5 excited disc, v=0).