The Fermi-Pasta-Ulam problem


The Fermi-Pasta-Ulam problem is named after the numerical experiments first performed by Enrico Fermi, John Pasta, and Stanislaw Ulam in 1954–1955 on the Los Alamos MANIAC computer, one of the first electronic computers. They wanted to understand how a crystal evolves toward thermal equilibrium by simulating a chain of particles coupled by springs.

    Suppose that we start with a simple model of a crystal, a system of particles of mass m connected by harmonic springs – a model inspired by the Debye theory of a crystal. The energy of such a system can be expressed as:

where ui is the displacement of particle i from its equilibrium position, pi is its momentum, and k is the usual spring constant. The two ends of the chain are assumed to be fixed and hence u0 = uN = 0.

    The usual way of understanding a system of masses connected by springs is to express the energy in terms of the normal modes. The latter are related to the displacement by

The energy of the system can be expressed in terms of the normal mode amplitudes as


    The equipartition theorem in statistical mechanics says that the average contribution of each quadratic contribution to the energy is kT/2, where T is the absolute temperature and k is Boltzmann's constant. However, the energy in each normal mode remains constant because the normal modes are independent of each other. Hence, if we solve Newton's equations of motion for a system of masses connected by linear springs, the energy in each mode will be the same as it was initially. That is, if we initially put all of the energy into a single mode, the energy will stay in this mode indefinitely, and the equipartition theorem will not be satisfied.

    A system that is ergodic (or more accurately quasi-ergodic) is one for which all accessible microstates are equally probable after a sufficiently long time. If the system were ergodic, the energy would be equally distributed among the various Fourier modes. Hence, we conclude that a system of particles connected by linear springs is non-ergodic.

    What happens if we make the springs nonlinear? For example, suppose we take the potential energy to be

The last term leads to a coupling between the modes. Fermi, Pasta, and Ulam thought that this coupling would lead to the energy that was initially in a single mode to eventually drift into the other modes until the equipartition theorem is satisfied. However, when they did the simulation, they found to their great surprise that this energy was shared by only a few other modes; the remaining modes were hardly excited. Moreover, after a long time the initial state was almost completely recovered. This result, which is known as the FPU problem or the FPU paradox, shows that nonlinearity is not sufficient to guarantee the equipartition of energy. It turns out that ergodic behavior is only observed when the magnitude of the nonlinear term is more than a certain critical value.

The Program

The program uses the Verlet algorithm to solve Newton's equations of motion numerically. The user can select the initial mode (k ≥ 1) of the system. The energy in this mode (red) and the two neighboring modes k + 1 (blue) and k - 1 (green) are plotted. The mean energy in these three modes is given when the simulation is stopped.

    If the system is ergodic, then all particles will see the same average environment, and the time average of any physical quantity associated with individual particles will be the same for each particle if the time t is sufficiently long. The time average of this quantity f for particle i is defined as

The time average of f averaged over all particles is


    We will choose the physical quantity of interest to be the potential energy of a particle, ei, defined as

ei = Σi ≠ j u(rij),

where u(rij) is the potential energy of interaction between particle i and j.

    A measure of the difference between the mean value of ei compared to its average over all particles is the energy metric Ω(t), defined as


If the system is ergodic over the time interval t, then it can be shown that Ω(t) decreases as 1/t. The program plots 1/Ω(t) versus t.


  1. Consider the behavior of a system of harmonic springs by setting b = 0. Describe some of the normal modes for N = 4. Does the energy that is initially in a given mode change with time? What is the behavior of the metric?
  2. Simulate the system using the default parameters. What is the behavior of the metric for t < 2 × 105? What is the behavior of the energy in the three modes that are plotted? What is the behavior of the energies and metric for longer times? (The unit of time is such that k = m = 1.) The simulation can be speeded up by choosing the number of steps per display to be at least 1000 and closing the window that displays the displacements of the particles.
  3. Increase the value of b and determine the approximate value of b for which the energy in the various modes becomes approximately equal and the inverse metric increases approximately linearly.


Updated 2 March 2009.