The Fermi - Pasta - Ulam - Tsingou chain

Post your questions on this hwk in the blog at the end of the page.
If you hand your Hwk by e-mail, you can send it to A. Rosso, M. Civelli or G. Roux .

Context: The ergodic hypothesis

Consider a system described by a Hamiltonian dynamics

where qn(t) represents a positions and pn(t) a momenta and qn(0), pn(0) are the initial conditions. Here the total energy

is conserved and the Statistical Physics predicts that the microcanonical ensemble is based on the equipartition principle:
the probability density to visit a point (qn, pn) in the phase space is uniform on the energy shell E.

The assumption behind statistical physics is called ergodic hypothesis: the long-time average of an observable coincides with the ensemble average of statistical physics.In practice this means that we can replace the precise simulation of the system over a long time with the average over many independent realizations. In this exercise we will see that this assumption is not always correct even for large and interacting systems.

Part I: The harmonic chain

One prerequisite for the ergodic hypothesis is that the system has to explore the whole energy surface. A simple counter-example to this expectation is given by a chain of particles on a lattice connected to their neighbors with elastic springs (see figure below). In this systems the number of constants of motion is high enough to forbid the dynamics to explore the whole energy surface.

Particle are distributed on a lattice close to integer sites 0,1,...L. Each particle interacts with its closest neighbors with an elastic potential V. We consider the case with fixed boundary conditions.

A discrete elastic chain is shown in figure: L+1 particles are interacting on a one dimensional lattice. Each particle interacts via elastic springs with the left and right neighbor. The distance between the rest position of the particle n and its actual position is denoted by un (n=0...L).
  • The particles are initially at rest:

  • Fixed boundary conditions are set for the first and the last particle.

  • The Hamiltonian is

where the velocity is the moment conjugated to the position (taking mass equal to unity).

Q1: Using equation (1), show that the equations of motion are:

Fourier representation

It is useful to introduce a descriptio in terms of Fourier modes k = 0,1...L − 1:

which are compatible with the fixed boundary conditions.

The external image arrow-10x10.png Fourier transform is given by

The Hamiltonian can be re-written as:

Q2: Count the modes of strictly positive energy for the initial conditions described above.

Q3: At t=0, we prepare the system at energy E. If the equipartition of energy applies, what is the average energy associated to each mode k?
[the answer should be the result of an explicit calculation reported in your copy. Hand-waving justifications will not be accepted].

Numerical scheme for the time evolution

At t=0, the chain is set with this initial conditions:

The total energy is concentrated in the mode k=1, whose energy is

We may wonder whether the energy stays in this mode or decays on the others modes, as predicted by statistical physics. We now study this question numerically, by discretizing the equation of evolution (2) with a time step dt (Euler integration scheme).

import pylab,math,numpy
Ntime = 1000
L = 32
dt = 0.001
def HARMONIC(u,n):
    return u[n+1]+u[n-1]-2.*u[n]
def FT(e,k):
    return sum([e[n]*math.sin(n*k*math.pi/L) for n in range(1,L)])/math.sqrt(L/2.)
#initial condition for the position u and velocity v
u = [0]+[L*math.sin(n*math.pi/L) for n in range(1,L)]+[0]
v = [0. for n in range(L+1)]
mode1 = []
#dynamics with fixed boundary conditions: u[0] = u[L-1] = 0
for itime in range(Ntime):
    oldu = u[:]
    oldv = v[:]
    v = [0]+[oldv[n]+dt*HARMONIC(oldu,n) for n in range(1,L)]+[0]
    u = [0]+[oldu[n]+dt*oldv[n] for n in range(1,L)]+[0]
    k = 1; en1 = FT(v,k)**2/2.+2.*(math.sin(k*math.pi/(2*L)))**2*FT(u,k)**2
    #k = 2; en2 = FT(v,k)**2/2.+2.*(math.sin(k*math.pi/(2*L)))**2*FT(u,k)**2
pylab.ylabel('E_1, energy of mode 1')
time = [dt*itime for itime in range(Ntime)]
exact = [L*(L*math.sin(math.pi/(2*L)))**2 for itime in range(Ntime)]
Q4: Compare the expected and the observed outcome with a plot. Comment about the precision of the Euler integration scheme.
[maximum 2 sentences]

Q5 : Propose a more efficient scheme of discretization that you will adopt in the following. Compare the expected and the observed outcome with a plot. Comment with 1 ot 2 sentences. [Very Important Question]

Q6: Set a different initial condition u=[0]+[L/2.-abs(L/2.-n) for n in range(1,L)]+[0].
Plot some snapshots of the displacement u(t) and the time evolution of the energy Ek for modes k=1 to 4.
Comment about the validity of the ergodic hypothesis.

Part II: Questions about Python

Q7: Explain the difference between a list (e.g. [1,2,3]), a tuple (e.g. (1,2,4))?

Q8: Define a list a (e.g. a=[1,2,2]). Define b=a. Modify a. What happens to b?

Q9: If a and b are two lists, what are the differences between the commands a.append(b), a.extend(b) and a+b?

Part III: A counter-example by Fermi, Pasta, Ulam and Tsingou

In Part I we have seen that the ergodic hypothesis does not hold for integrable systems. It was taught that non-integrable systems would on the contrary be ergodic, until computers were invented...In 1955, Enrico Fermi, John Pasta, Stanislaw Ulam and Mary Tsingou proposed to add a slightly non-harmonic interaction, that would allow the mixing of the energy between the modes so as to recover the microcanonical prediction of statistical physics. The equations of motions that they considered are

The numerical study of this equation was performed with one of the very first computers, the Maniac I. To repeat the same experiment, you can use the following force
def FPU(u,n):
    return u[n+1]+u[n-1]-2.*u[n]+alpha*((u[n+1]-u[n])**2-(u[n]-u[n-1])**2)
instead of HARMONIC(u,n). We start with the following parameters:
Ntime = 10000
L = 32
alpha = 0.3
dt = 0.2
and this initial state
u = [0]+[math.sin(n*math.pi/L) for n in range(1,L)]+[0]
Q10: What do you observe? Plot the energy of the first modes as a function of time and explain why Fermi, Pasta, Ulam and Tsingou were happy.
[1 or 2 sentences.]

Q11: Take a larger number of time steps Ntime=44000. Plot the energy of the first modes as a function What do you really conclude about the ergodic hypothesis?

Difficult Part : Study of solitons

This part is difficult and not compulsory.
Here we ask you to solve the Bonus by producing very nice plots (or movies) that illustrate soliton propagation.

The complicated evolution observed for the anharmonic chain corresponds to the propagation of many not-deforming profiles, called "solitons". In some limit it is possible to construct a mapping between the present problem and the Korteweg-de Vries equation. Here we sketch the main steps of the demonstration and recall the form of the traveling wave solution (the soliton) of the non-linear equation. We ask you to test the stability of these solutions with the program that you wrote in part II.

  • You first take a macroscopic limit by defining continuous variables.

  • Then you set a field Φ(x,τ) from the definition

  • At second order the field Φ(x,τ) verifies the equation

  • Assuming that variations in time are slow, we neglected the term

  • Setting X = x − τ one can show that

  • The equation above is called the Korteweg-deVries equation. Find a traveling wave solution, i.e. of the form:

Bonus 0: Find C0 and C1 and come back to the initial variable Φ(x,τ)
[Hint: when integrating Θ to obtain Φ(x,τ), ensure that the profile of Φ goes to zero at infinity].

Bonus 1: Test numerically the stability of the obtained soliton by plotting its time evolution.
[Don't forget to compute the initial velocity correctly.]

Bonus 2: How do two colliding solitons of opposite velocities interact? Illustrate this fact numerically by plotting the time evolution.


print this page