# The Fermi - Pasta - Ulam - Tsingou chain

DM2 The Fermi - Pasta - Ulam - Tsingou chain

Don't hesitate to ask questions and make remarks on this wiki page.

# Context: The ergodic hypothesis

Consider a system described by a Hamiltonian dynamics

$\dot p_n = - \frac{\partial\mathcal H}{\partial q_n}\quad$
$\dot q_n = \frac{\partial\mathcal H}{\partial p_n} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(1)$
where qn(t) represents a positions and pn(t) a momenta and qn(0), pn(0) are the initial conditions. Here the total energy
$E=\mathcal H(q_n(0),p_n(0))$
is conserved and Statistical physics predicts the microcanonical ensemble 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 approaches 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 [Group 1]

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:
$\quad v_n(0)=0$

• Fixed boundary conditions are set for the first and the last particle.
$u_0=u_L=0, \quad v_0= v_L=0$

• The Hamiltonian is
$\mathcal H(\{u_n,v_n\}) = \frac 12 \sum_{n=0}^{L-1} \Big[ v_n^2 + (u_{n+1}-u_n)^2\Big]$
where the velocity is the moment conjugated to the position (taking unit mass).

Q1: Using equation (1), show that the equations of motion are:
$\ddot u_n = u_{n+1}+u_{n-1}-2u_n \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(2)$

### Fourier representation

It is useful to introduce a descriptio in terms of Fourier modes k = 0,1...L − 1:
$\hat u_k = \sqrt{\tfrac 2L}\sum_{n=0}^{L-1} u_n \sin\frac{nk\pi}{L}\quad ; \quad\quad\quad \hat v_k = \sqrt{\tfrac 2L}\sum_{n=0}^{L-1} v_n \sin\frac{nk\pi}{L}$
which are compatible with the boundary conditions.

The Fourier transform is given by
$u_n = \sqrt{\tfrac 2L}\sum_{k=0}^{L-1} \hat u_k \sin\frac{nk\pi}{L}\quad ; \quad\quad\quad v_n = \sqrt{\tfrac 2L}\sum_{k=0}^{L-1} \hat v_k \sin\frac{nk\pi}{L}$

The Hamiltonian can be re-written as:
$\mathcal H(\{\hat u_k,\hat v_k\}) = \sum_{k=0}^{L-1} E_k = \sum_{k=0}^{L-1} \frac 12 \Big[ \hat v_k^2 + \omega_k^2 \hat u_k^2\Big] \quad\text{where}\quad \omega_k=2\sin\frac{k\pi}{2L}$

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 prepared as:
$u_n(t=0) ~=~ L\sin\frac{n\pi}{L} \quad\quad v_n(t=0) ~=~ 0$

The energy is concentrated in the mode k=1, of energy
$E(t=0)=E_1=L^3 \sin^2 \frac{\pi}{2L}$

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.
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
mode1.append(en1)
pylab.xlabel('time')
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)]
pylab.plot(time,exact,'b-')
pylab.plot(time,mode1,'r-')
pylab.axis([0,Ntime*dt,75.,82.])
pylab.show()

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 [Group 1]

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 [Group 1]

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 is

$\ddot u_n = u_{n+1}+u_{n-1}-2u_n + \alpha \Big[(u_{n+1}-u_n)^2-(u_{n}-u_{n-1})^2\Big] \quad \text{with} \; \alpha \ll 1$
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)

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 required.
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 undeforming 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 travelling wave solution (the soliton) of the non-linear equation. We ask you to test the stability of these solutions with the program you wrote in part II.

### From Fermi Pasta Ulam Tsingou to the Korteweg-de Vries equation

• You first take a macroscopic limit by defining continuous variables.
$x=\varepsilon n \quad \text{and} \quad\tau=\varepsilon , \quad \text{with } \varepsilon \ll 1$
• Then you set a field φ(x,τ) from the definition
$u_n(t)=\varepsilon \phi(\varepsilon n, \epsilon t)$
• At second order the field φ(x,τ) verifies the equation
$\partial_\tau^2\phi = \partial_x^2\phi+\bigg[2\alpha\partial_x\phi\partial_x^2\phi +\frac 1{12} \partial_x^4\phi\bigg]\varepsilon^2$
• assuming that variations in time are slow we neglected the term
$\partial_\tau^2 \phi$
• One sets X = x − τ and shows that
$\partial_\tau\Phi +\bigg[\alpha\Phi\partial_X\Phi +\frac 1{24} \partial_X^3\Phi\bigg]\varepsilon^2=0 \quad \text{where} \Phi = \partial_x \phi$
where
• This previous equation is called the Korteweg-de Vries equation. Find a traveling wave solution, i.e. of the following form:
$\Phi(X,t)= C_0 \Phi_0\big(C_1(X-v\tau)\big)$
where the invariant shape is
$\Phi_0(\bar X)=\frac 1{\cosh^2\bar X}$

Bonus 0: Find C_0 and C_1 and come back to the initial variable φ(x,τ) [Hint: when integrating from Φ to φ, ensure that the profile φ 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.