Wave 109 — 20 New Physics & Science Simulations (940 Total, −60 to 1000)

Every simulation on this site began as a question: what does this equation actually look like? The Debye length is taught in every plasma physics course as λ_D = √(ε₀k_BT / nq²), a formula that predicts how far an electrostatic perturbation penetrates a plasma before being screened. But until you watch a heatmap of the electric potential collapse exponentially from a test charge, the shielding remains abstract. Wave 109 asks that question twenty times across as many domains of physics and science.

Physics

Debye Length — λ_D Shielding Heatmap

In a plasma, free electrons rearrange themselves around any excess charge to neutralise it over the characteristic Debye length. The simulation places a test charge at the centre of a 2D plasma grid and renders the resulting electrostatic potential as a false-colour heatmap, showing the exponential decay φ(r) = (q / 4πε0r) exp(−r / λD). Sliders control electron temperature T and number density n, and the heatmap redraws in real time as λ_D scales with √(T/n). The transition from a diffuse, far-reaching potential at low density to a tightly confined hotspot at high density is immediately legible in colour. A companion plot shows the radial profile on a log scale, confirming the exponential envelope.

Topological Defects — 2D XY Model BKT Vortices

The Berezinskii-Kosterlitz-Thouless transition is one of the most unusual phase transitions in statistical mechanics: it occurs in two-dimensional systems where true long-range order is forbidden by the Mermin-Wagner theorem, yet something measurably sharp still happens at a critical temperature TBKT. The mechanism is the unbinding of topological defects — vortex-antivortex pairs in the 2D XY model, where each spin is a unit vector on the plane. Below TBKT, vortices and antivortices remain bound in pairs and the system has quasi-long-range order with algebraically decaying correlations. Above it, free vortices proliferate and correlations decay exponentially.

The simulation runs a Monte Carlo sweep of the 2D XY model on a square lattice, rendering spin orientations as arrows coloured by angle. Vortices are detected by computing the winding number around each plaquette and highlighted with red (positive) and blue (negative) markers. Dragging the temperature slider through TBKT produces a visible proliferation of unbound defect pairs.

Stochastic Resonance — Langevin Bistable SNR Curve

Stochastic resonance is the counterintuitive phenomenon in which adding noise to a system can improve its ability to detect a weak signal. A particle in a double-well potential driven by a weak sub-threshold periodic force cannot cross the barrier on its own — but at an optimal noise level, thermal fluctuations kick it across the barrier in synchrony with the forcing, producing a sharp peak in the signal-to-noise ratio as a function of noise amplitude. The Langevin equation governing the dynamics is:

dx/dt = −dV/dx + A cos(ωt) + η(t)
V(x) = −x²/2 + x⁴/4
⟨η(t)η(t′)⟩ = 2Dδ(t−t′)

The simulation integrates this equation numerically for a large ensemble of particles, sweeping the noise intensity D and computing the output SNR at the driving frequency using a windowed DFT. The resulting bell-shaped SNR curve is plotted in real time, with the optimal noise level marked. The stochastic resonance effect is relevant to sensory neuroscience, where neurons appear to exploit background noise to detect weak stimuli at threshold.

Granular Jamming — Bidisperse Packing φ_J ≈ 0.84

When a dense granular material such as sand, glass beads or wet coffee grounds is compressed or sheared, it undergoes a jamming transition: it suddenly acquires rigidity and resists deformation like a solid. The transition occurs at a critical packing fraction φ_J that depends on particle shape and size distribution. For a bidisperse mixture of discs (two particle sizes in a 1:1 number ratio chosen to suppress crystallisation), the jamming point lies near φ ≈ 0.84 in two dimensions.

The simulation uses a molecular dynamics engine with soft-sphere repulsive contacts to compress a bidisperse disc system at constant rate, tracking the pressure P and coordination number z as the packing fraction rises. The jamming transition is visible as a sharp upturn in pressure and a jump in z toward the Maxwell isostatic value z = 2d = 4. A live packing visualisation colours particles by their contact count, revealing the emergence of a rigid, force-chain-spanning network at φ_J.

Bragg Diffraction — nλ = 2d sinθ, FCC/BCC Powder Pattern

X-ray diffraction is the primary tool for determining crystal structure. When a monochromatic beam of wavelength λ strikes a crystalline lattice with plane spacing d at glancing angle θ, constructive interference occurs whenever Bragg's law is satisfied: nλ = 2d sinθ. Different families of planes (characterised by Miller indices hkl) contribute peaks at different angles, and the intensity of each peak depends on the structure factor, which encodes the atomic positions within the unit cell.

The simulation computes powder diffraction patterns for FCC, BCC and simple cubic structures, summing over all allowed reflections up to a user-defined maximum 2θ. The pattern is displayed as a stick diagram overlaid on a simulated Debye-Scherrer ring pattern. Switching between FCC and BCC changes the systematic absences immediately: FCC forbids mixed-parity hkl reflections, while BCC suppresses peaks where h+k+l is odd.

Spring Network — Maxwell Criterion z_c ≈ 4

James Clerk Maxwell showed in 1864 that a network of springs is mechanically rigid if and only if the mean coordination number z — the average number of springs per node — satisfies z ≥ 2d, where d is the spatial dimension. In two dimensions the critical threshold is zc ≈ 4. Below this value the network has floppy modes: large-scale deformations that cost no energy. Above it, the network is rigid and every deformation requires elastic energy.

The simulation generates random spring networks on a 2D grid with tunable connectivity and computes the normal modes using sparse linear algebra. Floppy modes are highlighted in red and rigid modes in blue, and the count of each is tracked against the Maxwell prediction. Clicking a node and dragging shows whether the network deforms rigidly or accommodates the displacement through zero-energy floppy rearrangements. The transition from floppy to rigid as bonds are added one by one is a vivid realisation of the Maxwell criterion.

Quantum

Quantum Dot Emission — Brus Equation, CdSe Spectra

Semiconductor quantum dots are nanoscale crystals small enough that quantum confinement shifts their band gap, making their emission colour a direct function of size. The Brus equation gives the size-dependent band gap of a spherical quantum dot of radius R:

E(R) = Eg,bulk + ħ²π² / (2μR²) − 1.8e² / (4πεε0R)

For CdSe quantum dots, the bulk band gap is 1.74 eV (near-infrared) and confinement blue-shifts the emission into the visible as the radius decreases from 5 nm to 1.5 nm, sweeping from red through orange, yellow and green. The simulation renders the emission spectrum as a Gaussian peak whose centre wavelength is determined by the Brus equation and whose width models phonon broadening. A palette of coloured circles represents dots of different sizes, glowing in the colour of their emission. Dragging a radius slider produces an immediate colour shift that exactly reproduces the rainbow palette seen in photographs of CdSe quantum dot solutions.

Spin Lattice — Heisenberg Ferromagnet Metropolis MC, T_c ≈ 2.27J/k

The Heisenberg model is the quantum-mechanical generalisation of the Ising model: each lattice site carries a spin vector S rather than a scalar ±1, and the Hamiltonian is H = −J Σ⟨ij⟩ Si · Sj. The 3D Heisenberg model has a ferromagnetic transition at Tc ≈ 2.27J/kB for the simple cubic lattice.

The simulation runs Metropolis Monte Carlo on a 3D Heisenberg lattice, rendering a 2D cross-section as a colour field (spin orientation mapped to hue). The magnetisation, specific heat CV and magnetic susceptibility χ are plotted as functions of temperature. Sweeping through Tc shows the order parameter vanish, the specific heat peak, and the susceptibility diverge — the classic signatures of a continuous phase transition.

Optics

Airy Disk — Bessel J₁, Rayleigh Criterion, 2D ImageData

When light from a point source passes through a circular aperture it does not form a perfect point image. Diffraction spreads the light into an Airy pattern whose intensity profile is:

I(θ) = I0 [2J1(ka sinθ) / (ka sinθ)]²

Rayleigh's resolution criterion states that two point sources are just resolved when the central maximum of one falls on the first zero of the other, at angular separation θ = 1.22λ/D. The simulation renders the 2D Airy pattern by evaluating the Bessel function on a pixel grid using the ImageData API, achieving sub-millisecond frame times for a 512×512 canvas. Two adjustable point sources can be placed side by side; their combined intensity pattern is displayed alongside the Rayleigh criterion line.

Chemistry

Reaction Front — Fisher-KPP PDE, v = 2√(Dr)

The Fisher-KPP equation couples logistic local growth with spatial diffusion:

∂u/∂t = D ∂²u/∂x² + ru(1−u)

For sharp initial conditions the asymptotic wavefront speed is v = 2√(Dr), independent of the initial profile. Both the monostable Fisher-KPP front and a bistable travelling wave variant are implemented with user-adjustable parameters. The 1D spatial profile is displayed as a moving intensity strip, the front position versus time is plotted to verify the asymptotic speed formula, and a phase-plane portrait in the co-moving frame shows the heteroclinic orbit connecting the two equilibria.

Thermodynamics

Entropy of Mixing — ΔS_mix = −nRΣx_i ln x_i, Particle Simulation

When ideal gases mix, the entropy increases by ΔSmix = −nRΣxi ln xi, always positive, reaching its maximum for equal proportions. The simulation places particles of two or three species in separate compartments and removes the dividing partition. Each particle performs a random walk on a 2D grid, and the instantaneous entropy of mixing is computed from the empirical composition in each cell. The approach to maximum entropy is plotted in real time, with the theoretical asymptote overlaid. Users can compare mixing with different diffusion coefficients to see how transport rates affect the timescale but not the thermodynamic endpoint.

Mathematics

Soap Bubble — Gauss-Seidel Laplace Relaxation, Minimal Surface

Soap films minimise their area subject to boundary constraints, satisfying the Laplace equation ∇²h = 0 for the height function. The Gauss-Seidel relaxation scheme iterates:

hi,j ← (hi+1,j + hi−1,j + hi,j+1 + hi,j−1) / 4

until convergence. The simulation lets users draw arbitrary boundary curves on a grid and watch the minimal surface relax into place, rendered as a height-coloured 2D map and optionally as a 3D mesh. Convergence speed with successive over-relaxation (SOR) is compared against standard Gauss-Seidel. The connection to electrostatics — where the same equation governs electric potential between conductors — is demonstrated by switching the boundary interpretation from height to voltage.

3D Cellular Automaton — Three.js InstancedMesh, 20³ Grid, 5 Rules

The simulation runs a 20×20×20 3D cellular automaton using five distinct rule sets drawn from the extended Birth/Survival/Neighbourhood notation. The live grid is rendered using Three.js InstancedMesh for 8000 cubes at interactive frame rates, with each cell's state mapped to opacity and colour. Rules include a 3D Game of Life analogue, a crystal-growth rule that produces dendritic forms, a chaos rule that fills the volume with aperiodic turbulence, and two rules that generate self-similar fractal structures. Internal structure is revealed by adjusting an opacity threshold that makes low-activity cells transparent.

Algorithms

Boolean Network — Kauffman NK, K = 2 Critical Phase Boundary

Stuart Kauffman's NK Boolean network model places N binary nodes, each receiving input from exactly K randomly chosen other nodes and updating via a random Boolean function. At the critical connectivity K = 2 the network sits at the boundary between order and chaos: attractors scale as √N, information propagates globally, and single-bit perturbations neither die out immediately nor cascade forever — properties Kauffman proposed as signatures of biological gene regulatory networks.

The simulation implements NK networks with adjustable N (up to 64) and K, visualising state space as a network of nodes coloured by binary state with attractor cycles highlighted. The Derrida annealed approximation for the phase boundary Kc = 1/(2p(1−p)) is overlaid on a K–p phase diagram, and single-node perturbation experiments directly probe the order-chaos boundary.

Life Science

Population Wave — Fisher-KPP 2D, Float32Array, Allee Effect

The 2D Fisher-KPP equation describes the spatial invasion of a species across a landscape. This simulation solves it on a Float32Array grid using a finite-difference scheme, with population density u(x,y,t) rendered as a false-colour map. The Allee effect — reduced per-capita growth at low density — is incorporated via a modified growth term:

f(u) = ru(u−a)(1−u),   0 < a < 1

Below the Allee threshold a, local populations decline rather than grow, creating a minimum viable patch size and a critical wave speed below which the invasion fails. Comparison between simple logistic growth and the bistable Allee model shows how the critical threshold creates sharp invasion boundaries and extinction debt in fragmented landscapes.

Networks

Network Percolation — Erdős-Rényi, Union-Find, Giant Component

For the Erdős-Rényi random graph G(N,p), a giant connected component of size O(N) emerges sharply at mean degree ⟨k⟩ = pN = 1. The simulation adds edges one at a time and tracks component sizes using Union-Find with path compression and union by rank, achieving near-O(1) amortised merge time. The giant component fraction S = |Cmax|/N is plotted against ⟨k⟩ in real time, with the mean-field self-consistency prediction S = 1 − exp(−⟨k⟩S) overlaid. Nodes are sized proportionally to their component, making the moment the giant component bootstraps itself visually dramatic.

Generative Art

Random Growth — Eden + DLA Fractal, D ≈ 1.71 Box-Counting

The Eden model grows a compact cluster by randomly selecting a perimeter site; the result is roughly circular with a rough but non-fractal surface. Diffusion-limited aggregation (DLA) launches random walkers from far away and adds each where it first touches, producing a tenuous branching structure with fractal dimension D ≈ 1.71 in two dimensions. The simulation runs both models side by side and computes the fractal dimension by the box-counting method in real time, plotting log N(ε) against log(1/ε) with a regression fit that converges to the theoretical fractal dimension as the cluster grows. The visual contrast between the compact Eden cluster and the spidery DLA aggregate makes the connection between local growth rules and large-scale fractal geometry immediately apparent.

Space Science

Hall Thruster — E×B Drift, Isp, Thrust

Hall-effect thrusters are the dominant propulsion technology for satellite station-keeping and orbit-raising. Xenon ions are accelerated by an electric field E to exhaust velocities of 15–30 km/s (Isp 1500–3000 s), while electrons are trapped in closed Hall drift orbits v = E×B/B² by a transverse magnetic field that prevents them from neutralising the accelerating potential.

The simulation models particle trajectories of both species using a Boris pusher in the crossed fields of the discharge channel, computing time-averaged thrust and specific impulse. Animated electron drift rings visualise the Hall current. Sliders control discharge voltage, propellant mass flow rate, and channel geometry; thrust F = ṁve and Isp = ve/g0 update continuously, making concrete how small voltage changes produce large Isp changes at modest thrust cost.

Cosmic Ray Shower — Heitler Cascade, e⁻/γ/μ/π, NKG Profile

When a high-energy cosmic ray proton strikes the upper atmosphere it initiates an extensive air shower. The Heitler binary cascade model describes the electromagnetic component: each electron radiates a bremsstrahlung photon and each photon pair-produces an electron-positron pair, doubling the particle count every radiation length X0 until energies fall below the critical energy Ec and the shower reaches its maximum depth Xmax.

The simulation generates a branching cascade tree from a primary energy of 1–100 PeV, tracking electrons, photons, muons and pions in distinct colours. The shower profile (particle number versus atmospheric depth) follows the Nishimura-Kamata-Greisen (NKG) function, overlaid on the Monte Carlo result. The depth of shower maximum Xmax shifts with primary mass, demonstrating why it serves as a composition probe: proton showers are deeper and more fluctuating than iron showers of equal total energy.

Electronics

Kirchhoff's Laws — MNA Matrix, Gaussian Elimination, 4 Circuit Presets

Modified Nodal Analysis (MNA) encodes Kirchhoff's current and voltage laws simultaneously as a matrix system Ax = b:

[G  B] [v]   [i]
[C  D] [j] = [e]

where G encodes conductances, B and C encode voltage source constraints, v contains node voltages, j branch currents through voltage sources, and e source values. The simulation builds the MNA matrix for four preset circuits — voltage divider, Wheatstone bridge, ladder network, and a multi-source circuit with dependent sources — and solves by Gaussian elimination with partial pivoting. Node voltages are annotated on the schematic; branch currents are shown as directed arrows whose width is proportional to magnitude. Adjusting any component value updates the solution immediately, making the dependence of circuit behaviour on component ratios tangible rather than algebraic.

Try It Yourself

Six simulations to start with if you are new to this wave:

The full Wave 109 catalogue is browsable via the category index. All twenty simulations are also available in Ukrainian via the /uk/ language prefix.

What Comes Next

The 1000-simulation target is 60 away. Wave 110 will reach it. Current planning includes simulations of Rydberg atom blockade, atmospheric boundary layer turbulence, topological insulator edge states, optimal transport (Earth-mover distance), protein folding energy landscapes, and a revisit of the Lorenz attractor with full Lyapunov spectrum computation. The emphasis throughout is on simulations that make a specific quantitative prediction visible at interactive speed: not merely animations, but tools for building physical intuition about the governing equations.

If you have a simulation you would like to see built, or a model from your research area that is underrepresented in interactive browser science, the contact page is open. Suggestions from researchers and educators have shaped several waves already, and Wave 110 planning is not yet closed.