aidoesscience
aidoesscience › Reaction–diffusion
Investigating · emergence

Reaction–diffusion scale selection

Does a structureless chemical soup spontaneously organise into a pattern with a reproducible length scale — set by the chemistry, not by the random seed?

Reaction–diffusion scale selection simulation running in the browser

▶ Run the simulationSee the measured result

Measured by the lab
11.94815
Known value
11.947825
Relative error
2.75e-5

Units: dimensionless Gray–Scott length units (k_max = 0.5258853 rad/unit; F=0.0866, k=0.0605, Du=0.16, Dv=0.08)

How the lab tests it

Run the Gray–Scott two-species reaction–diffusion system on a toroidal grid from a noisy uniform soup; once it settles, measure the characteristic wavelength λ from the first off-zero peak of the spatial autocorrelation of the activator field, and re-seed to see whether λ recurs while the pattern never does.

What it looks for

a reproducible, seed-independent characteristic wavelength λ (same physics → same scale) — an emergent length with no clean closed-form value (this is the self-replicating/pulse regime, not a linear Turing bifurcation with an analytic wavelength)

Turing instability calculator — steady state, unstable band, fastest-growing wavelength & critical diffusion ratio

Diffusion smooths. That is the whole of everybody's intuition about it, and it is why Turing's 1952 answer to how a featureless embryo acquires a pattern was so hard to believe: two chemicals that each only ever spread out can, together, take a soup that heals itself when stirred and break it into stripes with a definite spacing. The simulation above runs nothing but that -- a five-point Laplacian and the Gray-Scott kinetics f = -u*v^2 + F(1 - u), g = u*v^2 - (F + k)*v -- and this page is the argument it is an instance of, assembled from the four rates in the boxes and nothing else. Every step of it is algebra rather than search, which is worth saying because the last step usually is not. The steady state is a root, not a Newton iteration: both nontrivial states obey u*v = F + k, which turns f = 0 into the quadratic F*u^2 - F*u + (F + k)^2 = 0, and its discriminant is an existence condition -- below 4(F + k)^2 = F the patterned branch annihilates in a fold and there is nothing left to destabilise. The Jacobian then collapses to something you can read off the rates: g_v = F + k exactly, tr J = k - v^2, det J = (F + k)(v^2 - F), no numerical derivative anywhere. The unstable band is where the determinant of J - q*diag(Du, Dv) goes negative, which is a quadratic in q, so the band edges are roots. And the fastest-growing mode -- the one nearly every treatment reports as an argmax found by scanning -- comes out of a quadratic too, because setting the derivative of the larger eigenvalue to zero and squaring the radical away once cancels it completely. That quadratic's leading coefficient is -Du*Dv*(Du - Dv)^2, so it vanishes precisely when the two diffusivities are equal: the equation for the pattern's scale stops being a quadratic at exactly the point Turing proved has no pattern. The degeneracy is the theorem. Two more things fall out for free. The CRITICAL DIFFUSION RATIO is a property of the chemistry alone -- write d = Du/Dv and a factor Dv^2 cancels off both sides of the onset condition, leaving (f_u + d*g_v)^2 = 4*d*det J with no diffusivity in it, so making both species a thousand times more mobile does not move the threshold by a hair. And scaling both diffusivities by s is an exact rescaling of space, so the selected wavelength goes as sqrt(s) point for point. Against this lab's own reading the page has something sharper than agreement to offer. The oracle stepped only the raw PDE, measured 51 modes' growth rates and recovered 11.94815 +/- 0.0000265 length units; the argmax closed form gives 11.9478245, twelve standard errors away and fully inside the cubic-vertex bias the finding already discloses -- while the OTHER spelling in the literature, the minimum of h(q) that many texts quote as though it were the same number, gives 10.3957 and is fifty-eight thousand standard errors out. A simulation that resolves its own peak this well does not merely confirm the mechanism; it tells two routinely-confused formulas apart. Four things this page will not do. It will not name the geometry: stripes, spots and hexagons are chosen by the nonlinear terms among modes that share the band, and nothing linear distinguishes them. It will not predict an amplitude, or a wavelength once the amplitude is large. It will not scale Du and Dv by different factors and call the result a power law, because that moves the ratio and hence the threshold too. And it will not explain the shipped module's own on-screen number: that world runs the far-from-equilibrium Pearson regime F = 0.0545, k = 0.062, where tr J is POSITIVE and the homogeneous state is not stable to stirring at all, so Turing's precondition fails and no Turing wavelength exists there. The band arithmetic still runs at those rates and returns 12.07 cells, within 8% of the 13.07 on screen -- which is exactly the trap, and why the refusal is printed with the trap's number beside it rather than quietly.

u* = (1 − √(1 − 4(F+k)²/F))/2, v* = (F+k)/u* · g_v = F+k, f_v = −2(F+k), tr J = k − v*², det J = (F+k)(v*² − F) · σ(q) = ½(T + √(T² − 4h)), T = tr J − (Du+Dv)q, h(q) = Du·Dv·q² − Bq + det J, B = Dv·f_u + Du·g_v · band: q± = (B ± √(B² − 4Du·Dv·det J))/(2Du·Dv) · −Du·Dv(Du−Dv)²q² + 2Du·Dv((Du+Dv)trJ − 2B)q + ((Du+Dv)²detJ − B(Du+Dv)trJ + B²) = 0 ⇒ q_max, λ = 2π/√q_max · critical ratio: (f_u + d·g_v)² = 4d·det J · λ(s·D) = λ·√s

—
This simulation has a catalogued, oracle-checked result: Reaction–diffusion / Turing scale selection.