In 1955 Enrico Fermi, John Pasta, Stanislaw Ulam (and programmer Mary Tsingou) ran one of the first computer simulations in physics. They expected a chain of weakly nonlinear oscillators to gradually reach thermal equilibrium — energy spreading equally across all normal modes. Instead, after visiting a few higher modes, the energy returned almost entirely to the first mode: the FPU recurrence. This violated the prevailing assumption that weak nonlinearity alone is sufficient for thermalization.
Set β (nonlinearity) to zero for a perfectly linear chain — mode 1 keeps all its energy forever. Increase β to watch energy spill into modes 2, 3, 4 … then recur. Raise the initial energy E₀ for stronger nonlinearity and faster (but less exact) recurrences. Track more modes to see the full spectrum. The yellow dashed line marks the equipartition level — the height each bar would reach in true thermalization. Click Mode 2 to start with mode 2 excited instead.
The FPU paradox was resolved in the 1960s by the discovery of solitons: localized nonlinear waves that pass through each other without breaking up. In the FPU chain, energy concentrates into near-soliton structures that periodically return it to the initial mode. At high enough energy the system eventually does thermalize — but over timescales far longer than originally expected. The problem spurred the entire field of nonlinear dynamics and integrable systems.
This simulation integrates the FPU β-model: a chain of N=32 particles connected by springs that have both linear and cubic (β) restoring forces. The equations of motion are ẍⅰ = (xᵢ₋₁ − xᵢ)(1 + β(xᵢ₋₁ − xᵢ)²) + (xᵢ₋₁ − xᵢ)(1 + β(xᵢ₋₁ − xᵢ)²), integrated with a 4th-order Runge-Kutta scheme for accuracy. Energy is initially loaded into a single normal mode by displacing particles in its sinusoidal eigenfunction pattern. The modal energies are computed via a sine-transform of the particle displacements.
Three visualizations update simultaneously: the particle chain shows actual displacements colored by speed; the bar spectrum shows how energy is distributed across the first N modes (dashed yellow line = equipartition level); and the time-history traces show each mode’s energy over the last 300 time-samples. When β = 0 only mode 1 has energy and it stays there. With β > 0 energy leaks into modes 2, 3, 4 … and then returns — the recurrence. Recurrence is detected when mode-1 carries >70% of total energy.
What is the Fermi-Pasta-Ulam recurrence?
It is the surprising phenomenon where energy injected into the lowest normal mode of a nonlinear oscillator chain spreads to higher modes but then returns nearly completely to the first mode, rather than distributing evenly (thermalizing) as statistical mechanics predicts.
What does the β parameter control?
β sets the strength of the cubic nonlinearity in the spring force. At β=0 the chain is perfectly linear and energy never leaves mode 1. As β increases, coupling between modes grows, energy spreads faster, recurrences become less exact, and at very high β the system may eventually thermalize over long timescales.
What is modal energy and how is it computed?
A normal mode is a collective oscillation pattern of the whole chain. Its energy is computed by projecting the particle displacements and velocities onto the mode’s sinusoidal eigenfunction via a discrete sine transform, then evaluating the kinetic plus potential energy in that mode coordinate.
Classical statistical mechanics and the equipartition theorem predict that weakly coupled systems should reach thermal equilibrium, with energy equally shared among all degrees of freedom. Fermi expected this to happen quickly. Instead the energy kept returning to the first mode, suggesting that weak nonlinearity does not guarantee thermalization on accessible timescales.
It marks the equipartition energy level — the height each bar would reach if the total energy were equally distributed among all tracked modes. Real thermalization would bring all bars to this line. The FPU recurrence keeps bar 1 periodically much taller than the equipartition level.
In 1965 Zabusky and Kruskal showed that the continuous limit of the FPU chain is the Korteweg-de Vries (KdV) equation, which supports soliton solutions. These localized waveforms travel without dispersion and pass through each other elastically. In the discrete chain, near-soliton structures carry energy without sharing it, explaining the recurrence.
With strong nonlinearity the system eventually escapes the near-integrable soliton regime. Recurrences become incomplete and irregular, and over very long simulation times the system does approach equipartition — true thermalization. The crossover energy is called the FPU threshold or stochasticity threshold.
RK4 (fourth-order Runge-Kutta) is used with a small timestep, conserving the total Hamiltonian energy to high precision. You can verify this by watching the Total energy readout — it should remain nearly constant throughout the simulation.
The original 1955 FPU calculation used N=32 or N=64. With 32 interior particles the recurrence period is short enough to be visible in real time, and the modal structure is easy to display. Fewer particles make modes fewer and the spectrum cruder; more particles slow down the computation.
The α-model adds a quadratic nonlinear term (αΔx²) to the spring force; the β-model used here adds a cubic term (βΔx³). The β-model is simpler to analyze because it preserves the reflection symmetry x → −x of the chain. Both models show the recurrence phenomenon, but with different recurrence periods and mode-coupling patterns.
The FPU problem launched the study of integrable systems, Hamiltonian chaos, and computational physics as a discipline. It influenced the KAM theorem (Kolmogorov-Arnold-Moser) on persistence of quasi-periodic orbits, and connected to studies of discrete breathers, rogue waves, and nonlinear lattice dynamics still active today.