Toys

Reaction Diffusion

Click the canvas to start/pause the simulation.

Every simulation here uses periodic boundary conditions, so the domain is topologically a torus. The View button chooses the geometry: Flat view uses the flat Laplacian, while Torus view renders the surface in 3D and swaps in the torus Laplace-Beltrami operator (in conservative flux form), so diffusion feels the donut's surface geometry. Drag the torus to rotate it. Toggle views mid-simulation and watch the pattern re-equilibrate. Vector and tensor models (Oseen-Frank, Q-Tensor) get the scalar operator applied componentwise; proper tangent-plane transport is a future upgrade.

These simulations are based on reaction-diffusion models of the form

$$\frac{\partial \mathbf{u}}{\partial t} = \gamma \Delta \mathbf{u} + g(\mathbf{u})$$

where $\mathbf{u}$ represents one or more chemical concentrations or phase fields evolving in a spatial domain. Different choices of the reaction term $g$ and the number of fields give rise to qualitatively different patterns.

The Gray-Scott model tracks two chemicals $u$ and $v$ with the stoichiometry

$$\begin{array}{rcl}U + 2V &\rightarrow& 3V,\\ V &\rightarrow& P.\end{array}$$

Species $U$ is fed into the system at rate $f$ and $V$ is removed at kill rate $k$, giving the PDE

$$\left\{\begin{array}{rcl}\dfrac{\partial u}{\partial t} &=& \gamma_u \Delta u - uv^2 + f(1-u),\\[6pt]\dfrac{\partial v}{\partial t} &=& \gamma_v \Delta v + uv^2 - (f+k)v.\end{array}\right.$$

Different $(\gamma_u,\gamma_v,f,k)$ produce spots, stripes, spirals, and coral-like branching -- all called Turing patterns. The display shows $|u - v|$.

Numerical scheme: Forward Euler, 9-point Laplacian, periodic BCs.

References: MIT GrayScott, Karl Sims, Pearson (1993).

Cahn-Hilliard (Cahn and Hilliard, 1958) describes phase separation in a binary mixture through a conserved order parameter $\phi \in [-1,1]$ that labels the two phases:

$$\left\{\begin{array}{rcl} \dfrac{\partial \phi}{\partial t} &=& \Delta\mu,\\[8pt] \mu &=& -\varepsilon^2\Delta\phi + W'(\phi), \end{array}\right. \quad W(\phi) = \frac{(1-\phi^2)^2}{4}.$$

The chemical potential $\mu$ is the variational derivative of the Ginzburg-Landau energy, so Cahn-Hilliard is the $H^{-1}$ gradient flow of $$E[\phi] = \int_\Omega \frac{\varepsilon^2}{2}|\nabla\phi|^2 + W(\phi)\,d\Omega.$$ Unlike Allen-Cahn, the total mass $\int \phi$ is conserved -- the two phases can only rearrange, not appear or disappear. Starting from near-uniform noise (Spinodal preset), the mixture spontaneously separates into interconnected bicontinuous domains, then slowly coarsens: smaller features dissolve and feed the growth of larger ones. In the Off-Critical preset the minority phase instead forms isolated droplets. The Ostwald Ripening preset quenches an off-critical mixture straight into a dense emulsion of small $+1$ droplets, then lets it coarsen: larger droplets grow while smaller ones dissolve and vanish, all while the total droplet area stays fixed -- the conserved-dynamics counterpart of the Allen-Cahn coarsening above, where domains simply shrink and disappear. Dark = phase $-1$, bright = phase $+1$.

Numerical scheme: Two-pass forward Euler split -- pass 1 computes $\mu$ from $\phi$ via the chemical potential equation and stores $(\phi, \mu)$ in an intermediate texture; pass 2 samples that texture to advance $\phi$ via $\phi_t = \Delta\mu$. Periodic BCs. Only conditionally stable -- large $\Delta t$ or small $\varepsilon^2$ can cause blow-up.

Research connection: Cahn-Hilliard is a central model in Justin's dissertation on energy-stable schemes for gradient flow problems. The scheme here is forward Euler -- only conditionally stable and not guaranteed to decrease $E$. The dissertation develops semi-implicit schemes that are unconditionally energy-stable. See /p/research/.

References: Cahn and Hilliard (1958), J. W. Cahn (1961), L.-Q. Chen (2002).

The Oseen-Frank model describes a nematic liquid crystal using a director field $\mathbf{n}$ -- a unit vector giving the local average molecular orientation (with $\mathbf{n} \equiv -\mathbf{n}$ since molecules have no head or tail). In the one-constant approximation all elastic modes cost the same energy, giving

$$E[\mathbf{n}] = \frac{K}{2}\int_\Omega |\nabla\mathbf{n}|^2\,d\Omega.$$

The overdamped (no-flow) gradient flow of this energy subject to $|\mathbf{n}|=1$ is

$$\frac{\partial \mathbf{n}}{\partial t} = K\bigl(\nabla^2\mathbf{n} - (\mathbf{n}\cdot\nabla^2\mathbf{n})\,\mathbf{n}\bigr).$$

The projection term $-(\mathbf{n}\cdot\nabla^2\mathbf{n})\,\mathbf{n}$ enforces the unit-norm constraint. The display maps director angle $\theta$ to hue with period $\pi$ (reflecting $\mathbf{n}\equiv-\mathbf{n}$); defect cores appear as points where the hue cycles rapidly.

Topological defects are points where $\mathbf{n}$ cannot be continuously defined. In 2D nematics the fundamental defects have half-integer winding numbers $\pm\tfrac{1}{2}$. Integer ($\pm 1$) defects also appear; the $+1$ radial ("hedgehog" or "aster") and the $-1$ hyperbolic ("saddle") are shown in the Hedgehog preset as an attracting pair -- watch them collide and annihilate. Annihilation shows two $+\tfrac{1}{2}$ and two $-\tfrac{1}{2}$ defects in a symmetric cross pattern; each $+\tfrac{1}{2}$ is nearest to a $-\tfrac{1}{2}$, so two simultaneous annihilation events occur. Coarsening starts from random directors -- many defect pairs nucleate and annihilate until only a few survive.

Numerical scheme: Forward Euler, 9-point Laplacian on each director component, with projection and renormalization each step. Stable for $K \cdot \Delta t \le 0.5$ (stencil max eigenvalue $-2$); current parameters give margin $1 - K\Delta t \cdot 2 = 0.2$.

Research connection: The Oseen-Frank model breaks down inside defect cores where $|\mathbf{n}|\to 0$. The full Landau-de Gennes Q-tensor model regularizes the core by allowing the order parameter to melt there -- that is the model in Justin's dissertation work, and you can run it here: select Q-Tensor (LC) in the model dropdown above. See /p/research/.

References: Oseen (1933), Frank (1958), de Gennes and Prost "The Physics of Liquid Crystals" (1993).

The Landau-de Gennes Q-tensor model describes a nematic liquid crystal with a symmetric, traceless tensor order parameter instead of a unit director. In 2D the tensor has two independent components,

$$Q = \left(\begin{array}{cc} q_1 & q_2 \\ q_2 & -q_1 \end{array}\right), \qquad q_1 + i\,q_2 = \tfrac{S}{2}\,e^{2i\theta},$$

where $\theta$ is the director angle and $S$ is the scalar order parameter -- how aligned the molecules are. Crucially, $S$ is free to vary: $S \to 0$ means the material locally melts into the isotropic phase. Because the director angle enters doubled, $\mathbf{n} \equiv -\mathbf{n}$ is built into the representation and half-integer defects are native -- no branch cuts. In strict 2D, $\mathrm{tr}(Q^3) = 0$ identically, so the one-constant Landau-de Gennes energy reduces to

$$F[Q] = \int_\Omega \frac{L}{2}|\nabla Q|^2 + \frac{a}{2}\,\mathrm{tr}(Q^2) + \frac{c}{4}\bigl(\mathrm{tr}\,Q^2\bigr)^2\,d\Omega,$$

and its $L^2$ gradient flow, written in the components $q = (q_1, q_2)$, is just a two-species reaction-diffusion system:

$$\frac{\partial q_i}{\partial t} = L\Delta q_i - a\,q_i - 2c\,|q|^2 q_i.$$

With $a < 0$ the isotropic state is unstable and the system orders with $|q|_{eq} = \sqrt{-a/2c}$ -- this is exactly the vector Ginzburg-Landau equation, the tensorial cousin of Allen-Cahn above.

Quench drops the system from the isotropic phase into the nematic: order develops locally, incompatible regions meet, and a gas of $\pm\tfrac{1}{2}$ defects nucleates and coarsens by pair annihilation. +1 Splits in Two is the showpiece: a $+1$ radial defect costs elastic energy proportional to its winding number squared, so splitting into two $+\tfrac{1}{2}$ defects (each costing a quarter) halves the energy. Watch the core fission and the halves repel. The vector Oseen-Frank model above cannot do this -- a unit vector field has no half-integer windings, so its $+1$ is artificially stable. Pair Annihilation shows a $+\tfrac{1}{2}$ and $-\tfrac{1}{2}$ attracting and annihilating.

Displays: Director maps director angle to hue (period $\pi$) and dims by $S$, so defect cores appear as dark points with hue pinwheels around them. Order shows the $S$ field directly -- watch the cores melt. Schlieren simulates the view between crossed polarizers, $I \propto S^2\sin^2(2\theta)$: dark brushes meet at defects, two brushes for half-integer charges and four for integer ones, just like a real microscope image.

Numerical scheme: Forward Euler, 9-point Laplacian on each component, periodic BCs. Conditionally stable ($\Delta t \, L \lesssim 1.25$ for this stencil); defect core radius $\sim\sqrt{L/|a|}$ is about 2 pixels at the default parameters.

Research connection: This is the model from Justin's dissertation -- there in 3D, with energy-stable semi-implicit schemes that unconditionally dissipate $F$ instead of the conditionally stable forward Euler used here. See /p/research/.

References: de Gennes and Prost "The Physics of Liquid Crystals" (1993), Mottram and Newton "Introduction to Q-tensor theory" (2014), Schopohl and Sluckin (1987).

The complex Ginzburg-Landau equation governs a complex field $A = u + iv$ and is the universal amplitude equation: near any Hopf bifurcation, an oscillatory medium reduces to it, which is why so many of the spiral patterns on this page have it lurking underneath. In one-constant form,

$$\frac{\partial A}{\partial t} = A + D(1 + i\alpha)\,\nabla^2 A - (1 + i\beta)\,|A|^2 A.$$

The linear term grows any small perturbation; the cubic term saturates it at $|A| = 1$. The two real parameters $\alpha$ (linear dispersion) and $\beta$ (nonlinear dispersion) control everything. The field is stored as its real and imaginary parts in the same two channels the Q-tensor uses -- and indeed this is the oscillatory, non-variational cousin of that gradient flow: same complex field, completely different dynamics.

The whole phase diagram turns on the sign of $1 + \alpha\beta$, the Benjamin-Feir-Newell criterion. When $1 + \alpha\beta > 0$, plane waves are stable and the field freezes into a Spiral Lattice: a tiling of rotating spirals, each a phase defect where $A = 0$. When $1 + \alpha\beta < 0$, plane waves are unstable and the medium turns chaotic. Just past the line you get Phase Turbulence -- the amplitude stays near $1$ while the phase wanders, so a few large spirals persist in a restless background. Push further and amplitude defects proliferate into full Defect Turbulence: a spatiotemporal chaos of spirals endlessly nucleating and annihilating. All three grow from the same tiny random seed.

Displays: Phase maps $\arg A$ to hue and dims by the amplitude $|A|$, so defect cores show as dark points. Amplitude shows $|A|$ directly -- the holes are the defects. Pure Phase drops the dimming for an undimmed phase portrait.

Numerical scheme: Forward Euler, 9-point Laplacian, periodic BCs. The dispersive $(1+i\alpha)$ Laplacian tightens stability to $\Delta t\, D \, |\lambda| < 2/(1+\alpha^2)$, so the turbulent presets ($\alpha = 2$) step with a smaller $\Delta t$ than the lattice ($\alpha = 0$).

References: Aranson and Kramer, "The world of the complex Ginzburg-Landau equation," Rev. Mod. Phys. (2002); Cross and Hohenberg (1993).

Several models here connect directly to Justin's dissertation on energy-stable schemes for gradient flow problems. Allen-Cahn and (upcoming) Cahn-Hilliard are canonical gradient flows; the versions in this toy use forward Euler for simplicity. The energy-stable semi-implicit schemes that unconditionally preserve the free energy appear in the dissertation -- see /p/research/.

Source code: github

Toys

Reaction Diffusion

Navigation