# Fermi-Pasta-Ulam nonlinear lattice oscillations

Post-publication activity

Curator: Stefano Ruffo

The Fermi-Pasta-Ulam problem bears the name of the three scientists who were looking for a theoretical physics problem suitable for an investigation with one of the very first computers, the Maniac. They decided to study the thermalization process of a solid. The Fermi-Pasta-Ulam (FPU) problem, first written up in a Los Alamos report in May 1955 (Fermi, Pasta, Ulam 1955), marked the beginning of both a new field, nonlinear physics (this problem is of central importance in the theories of solitons and chaos), and the age of computer simulations of scientific problems. Mary Tsingou (now Mary T. Menzel), who took part in the numerical study, is not an author of the report, but her contribution is recognized by two lines of acknowledgment. It is time for a proper recognition of her contribution. Let us refer from now on to the Fermi-Pasta-Ulam-Tsingou problem (Dauxois 2008). Figure 1: Schematic picture of the FPU model: masses that can move only in one dimension are coupled by nonlinear springs. $$u_n$$ is the relative displacement with respect to the equilibrium position of the $$n$$-th mass. The two ends of the chain were assumed to be fixed, i.e., $$u_0 = u_{N} = 0\ .$$ Figure 2: FPU recurrence for a FPU-$$\alpha$$ model with $$N=32$$ masses and fixed ends. The plot shows the time evolution of the energy (kinetic + potential) $$E_k=(\dot A_k^2+\omega_k^2A_k^2)/2$$ of each of the three lowest normal modes, related to the displacements through $$A_k = \sqrt{2/(N+1)}\sum_{n=1}^N u_n \sin\left({nk\pi}/(N+1)\right)$$ with the frequencies $$\omega_k^2 = 4 \sin^2\left(k\pi/(2N+2)\right)\ .$$ Initially, only mode $$k=1$$ (blue) is excited. After flowing to other modes, $$k=2$$ (green), $$k=3$$ (red), etc., the energy almost fully returns to mode $$k=1\ :$$ this was a surprise! This picture might be easily reproduced using the MATLAB code provided below.

## Contents

The original idea, proposed by Enrico Fermi, was to simulate the one-dimensional analogue of atoms in a crystal: a long chain of particles linked by springs that obey Hooke’s law (a linear interaction), but with a weak nonlinear correction (quadratic for the FPU-$$\alpha$$ model or cubic for the FPU-$$\beta$$ model), see Figure 1.

A purely linear law for the springs guarantees that energy given to a single 'normal' mode always remains in that mode (see caption of Figure 2 for the definition of normal modes in terms of atom displacements from their equilibrium positions).

Fermi, Pasta and Ulam thought that, due to the nonlinear correction, the energy introduced into the lowest frequency mode $$k=$$1 should have slowly drifted to the other modes, until the equipartition of energy, a consequence of ergodicity, would have been reached. The beginning of the calculation indeed suggested that this was the case. Modes $$k=$$2, $$k=$$3,…, were successively excited, reaching a state close to equipartition, as shown in Figure 2. However, by accident, one day, they let the program run longer. When they realized their oversight and came back to the computer room, they noticed that the system, after remaining in the near equipartition state for a while, had then departed from it. To their great surprise, after 157 periods of the mode $$k=$$1, almost all the energy (all but 3%) was back to this mode. Further calculations, performed later with faster computers, showed that the same phenomenon repeats many times, and that a super-recurrence exists, at which the initial state is recovered with an even higher accuracy (see Ford, 1992 and Weissert, 1997 for a historical account).

Therefore, the system behaved in a surprising way: contrary to the predictions of statistical mechanics when the number of particles$$N$$ is going to infinity, the energy equipartition state was not reached and energy was periodically returning to the initially excited $$k=$$1 mode. This highly remarkable result, known as the FPU paradox, shows that nonlinearity is not enough to guarantee the equipartition of energy.

The solution of the paradox is two-fold. It requires the consideration of two important physical phenomena: the existence and stability of solitons in one dimension and the presence of deterministic chaos.

### In real space: generation of solitons Figure 3: Space-time evolution of the initial condition $$w(x,0)=\sin(\pi x)$$ in a system with periodic boundary conditions, obeying the Korteweg-de Vries equation (2). Time and space are respectively plotted vertically and horizontally. The amplitude of the local velocity $$w(x,t)$$ is shown in color scale. Insets on the right present successive time snapshots of $$w(x,t)\ .$$ The top one corresponds to the initial condition and the second shows the development of a shock wave. One distinguishes on the third the formation of a train of solitons. Courtesy by N. Zabusky.

In 1965, using the so-called continuum limit, Zabusky and Kruskal (1965) were able to relate the periodic behavior observed in Figure 2 to the dynamics of localized excitations, nowadays known as solitons. Soliton generation was found for periodic boundary conditions, a situation which shares some features with the original fixed boundary conditions problem (Chirikov et al. 1973). Starting from the FPU-$$\alpha$$ equations of motion, $$\tag{1} \ddot u_n=(u_{n+1}+u_{n-1}-2u_{n})+\alpha \left[(u_{n+1}-u_{n})^2-(u_{n}-u_{n-1})^2\right] \ ,$$

and restricting the investigation to long wavelength modes and right-propagating waves, they derived the Korteweg-de Vries (KdV) equation (Korteweg and De Vries 1895) $$\tag{2} w_\tau+\frac{1}{24}w_{\xi\xi\xi}+\alpha \,ww_\xi=0\ ,$$

where the field $$w$$ is a conveniently rescaled spatial derivative of the displacement field $$u_n\ ,$$ while $$\tau$$ is a rescaled time and $$w_x=\partial w/\partial x\ .$$ If the nonlinear force is cubic, i.e. for the FPU-$$\beta$$ model, one ends up with a modified KdV equation (mKdV).

If the amplitude of the field $$w$$ is large enough, the sinusoidal initial condition develops a sharp front (see Figure 3), which then breaks into a series of pulses, which are solitons. They preserve their shapes and velocities and, during their motion in the finite system with periodic boundary conditions, from time to time, they come back to the positions they had initially, restoring the initial condition. This is why FPU recurrence is observed. Quantitative estimates of the recurrence time were derived using soliton theory for the Toda lattice, which approximates well the FPU-$$\alpha$$ model for small nonlinearities (Toda 1978).

It is important to emphasize that the FPU paradox would probably not have been a mystery for more than 10 years if, before Zabusky and Kruskal, somebody had got the idea to look carefully at the dynamics of the nonlinear lattice as a function of the space coordinate.

### In normal mode space: chaotic properties and Chirikov's resonance overlap criterion

Early attempts to explain the inefficient energy transfer among normal modes using linear resonance theory failed. A partial success was obtained using an adapted Krylov-Bogoliubov averaging method, which provided an estimate of the recurrence time (Jackson, 1963). Kolmogorov-Arnold-Moser (KAM) theorem had a strong impact on the debate around the FPU problem. KAM theorem states that most orbits of a slightly perturbed integrable Hamiltonian remain quasi-periodic. It looked natural to interpret the recurrent motion of normal modes in the FPU model as a consequence of the presence of such regular orbits. However, up to now, no quantitative estimate of the FPU recurrence time has been obtained using KAM theory. Moreover, the applicability of KAM theorem to the FPU Hamiltonian has been assessed only recently for the fixed end FPU-$$\beta$$ model (Rink 2007), by proving that Birkhoff normal forms constitute an integrable approximation of the FPU Hamiltonian, the integrals being the linear energies of the normal modes. It has been also questioned whether the original FPU initial condition really lies or not on a KAM torus (Casetti et al 1997).

Compared with the rather inconclusive application of KAM theorem, a remarkable result was obtained by Izrailev and Chirikov (1966) applying the concept of overlap of nonlinear resonances. When the nonlinearity is small, one can separately consider each nonlinear resonance and apply perturbation theory. Besides KAM regular orbits, chaotic orbits with a positive Lyapunov exponent are present, but they are confined in thin layers and the Lyapunov exponent is small. On the contrary, for stronger nonlinearity, resonances are no more separated in frequency space and their overlap gives rise to strongly chaotic orbits, which have a larger Lyapunov exponent and diffuse in phase space. Using the resonance overlap criterion, Izrailev and Chirikov predicted that, as the nonlinearity of the FPU model is increased, by e.g. increasing energy, a transition should occur to a regime where relaxation to equipartition should be fast, as a consequence of the strongly chaotic nature of the orbits. This theoretical prediction was successfully tested numerically (Izrailev et al. 1970). Independendly, Bocchieri and collaborators (Bocchieri et al. 1970) showed the presence of the transition for a Lennard-Jones potential in one dimension. The strong stochasticity threshold (a term first used in Pettini and Landolfi 1990) exhibited in the numerical experiments of Livi and coworkers (Livi et al., 1985) is nothing but the transition from weak to strong chaos predicted by Chirikov. Nowadays we know that even below the strong stochasticity threshold, relaxation to equipartition occurs on longer and longer time scales as the energy is reduced, due to the overlap of higher order nonlinear resonances (De Luca et al. 1995, Shepelyansky, 1997).

The infinite-dimensional version of KAM-theory (KAM-tori are infinite-dimensional and stable) is found in the equation with FPU-nonlinearity by L. D. Pustyl'nikov in 1994 (See Reference below).

## Relaxation to equipartition for short wavelengths: chaotic breathers Figure 4: Time evolution of the $$\pi$$-mode, $$u_n(0)=a(-1)^n\ ,$$ for the FPU-$$\beta$$ model. In panel (a), the horizontal axis indicates lattice sites and the vertical axis is time. The grey scale goes from $$E_n=0$$ (white) to the maximum $$E_n$$-value (black) of the local energy defined in Eq.(3). The lower rectangle corresponds to $$0<t<3000$$ and the upper one to $$5.994\ 10^5<t<6.10^5\ .$$ Panels (b), (c) and (d) show the instantaneous local energy $$E_n$$ along the chain at three different times. Remark the difference in vertical amplitude in panel (c), when the chaotic breather is present, showing that it concentrates a large part of the energy initially given to the lattice.

If instead of preparing a long wavelength (low-frequency) initial state, one now puts the energy in the short wavelength (high-frequency) part of the normal mode spectrum (Zabusky and Deem 1967), the pathway to equipartition may lead to the creation of highly localized excitations that have an oscillating amplitude and bear most of the energy: they have been called chaotic breathers (Cretegny et al. 1998). Contrary to exact breathers (Flach and Willis 1998), which are periodic localised solutions whose existence has been proved also for FPU lattices (James 2001), these excitations display a chaotic evolution in space and time. This behavior is generic, because it is intimately related to modulational instability, a self-induced modulation of the steady state resulting from a balance between nonlinear and dispersive effects. This phenomenon appears in many branches of physics: fluid dynamics, nonlinear optics and plasma physics.

Figure 4 shows the time evolution of the $$\pi$$-mode, $$u_n(0)=a(-1)^n\ ,$$ above the modulation instability critical amplitude $$a>a_c \approx \pi/(\sqrt{6}N)\ ,$$ where $$N$$ is the number of masses in a chain with periodic boundary conditions. Interestingly, this threshold amplitude coincides with the threshold estimate obtained using the Chirikov overlap criterion for short wavelength modes (Dauxois et al 2005a). The energy residing on site $$n$$ is $$\tag{3} E_n = {1\over 2} \dot{u}_n^2 + {1\over 2}V(u_{n+1}-u_n) + {1\over 2} V(u_n-u_{n-1}),$$

where the FPU-$$\beta$$ potential is $$V(x) = {1\over 2}x^2 + {1\over 4} x^4\ .$$

Figure 4 displays the full evolution, from the initial state to the generation of the moving chaotic breather. Figure 4(b), Figure 4(c) and Figure 4(d) are three successive snapshots of the local energy $$E_n$$ along the chain. At short time, a slight modulation of the energy in the system appears (see Figure 4(b)) as a consequence of the destabilization of the $$\pi$$-mode. Later on, as shown in Figure 4(a), only a few breathers emerge. Inelastic collisions of breathers have a systematic tendency to favour the growth of the big breathers at the expense of small ones. Hence, in the course of time, the breather number decreases and only one, of very large amplitude, survives (see Figure 4(c)): this is the localized excitation called chaotic breather. The chaotic breather moves along the lattice with an almost ballistic motion: sometimes it stops or reflects back. During its motion the chaotic breather collects energy and its amplitude increases. It is important to note that the chaotic breather is never at rest and that it propagates with a given subsonic speed. Finally, the chaotic breather decays and the system reaches energy equipartition, as illustrated in Figure 4(d).

## Recent approaches

It has been observed that relaxation to equipartition for long wavelength initial conditions proceeds in two stages. On a short time scale a packet of low frequency normal modes is formed, with the higher modes cut-off exponentially (Fucito et al. 1982): this phase is associated with the formation of the soliton train of Figure 3. On longer times, energy is slowly flowing to high frequency modes. This scenario has been recently advocated by Galgani and coworkers (Berchialla et al 2004) on the basis of numerical simulations and of a resonant normal form approach to the KdV equation (Bambusi and Ponno 2006). Grist to this mill is also brought by the discovery of periodic orbits (so called q-breathers) that are exponentially localised in mode space (Flach 2005, Penati and Flach 2007).

An important open question is what remains of the FPU phenomenology in two-dimensional and three-dimensional lattices. Here results are scarce due to the numerical difficulty to simulate large lattices. Worth to mention is the work of Benettin (Benettin 2004), who argues against the two-stage relaxation picture.

## MATLAB Code

The following simple MATLAB code solves the equations for the dynamics of the FPU model and computes the normal modes (Dauxois et al. 2005b). It allows for the reproduction of different aspects of the FPU problem.

N=32;         %Number of particles must be a power of 2
alpha=0.25;   %Nonlinear parameter
TMAX=10000;
DT=20;
tspan=[0:DT:TMAX];
options=odeset('Reltol',1e-4,'OutputFcn','odeplot','OutputSel',[1,2,N]); % Test different tolerances, changing Reltol
for I=1:N,
a=1; b(I)=a*sin(pi*I/(N+1)); b(I+N)=0;  % FPU initial condition
%a=1; b(I)=a*sin(pi*N*I/(N+1)); b(I+N)=0; % Zabusky-Deem init. cond.
%k=0.8; sk=(sinh(k))\verb!^!2; ek=exp(-k); i1=I-N/4; i2=i1-N/2; %Solitons
%b(I)=-0.5/alpha*log((1+exp(2*k*(i1-1)))/(1+exp(2*k*i1)));  % Kink
%b(I)=b(I)+0.5/alpha*log((1+exp(2*k*(i2-1)))/(1+exp(2*k*i2)));  % Anti-kink
%b(I+N)= sk*ek/alpha/cosh(k*i1)/(exp(-k*i1)+exp(k*i1)*exp(-2*k));
%b(I+N)=b(I+N)-sk*ek/alpha/cosh(k*i2)/(exp(-k*i2)+exp(k*i2)*exp(-2*k));
omegak2(I)=4*(sin(pi*I/2/N))\verb!^!2;   % Mode Frequencies
end
[T,Y]=ode45('fpu1',tspan,b',options,N);  % Time integration
for IT=1:(TMAX/DT),
TIME(IT)=IT*DT*sqrt(omegak2(1))/2/pi;  % Time iteration loop
YX(IT,1:N+1)=[0 Y(IT,1:N)];
YV(IT,1:N+1)=[0 Y(IT,N+1:2*N )];
sXF(IT,:)=imag(fft([YX(IT,1:N+1) 0  -YX(IT,N+1:-1:2)]))/sqrt(2*(N+1));
sVF(IT,:)=imag(fft([YV(IT,1:N+1) 0  -YV(IT,N+1:-1:2)]))/sqrt(2*(N+1));
Energ(IT,1:N)=(omegak2(1:N).*(sXF(IT,2:N+1).\verb!^!2)+sVF(IT,2:N+1).\verb!^!2)/2;
for J=2:N-1,   % Space loop
DifY(IT,J)=Y(IT,J+1)-Y(IT,J);
end
end
plot(TIME,Energ(:,1),TIME,Energ(:,2),TIME,Energ(:,3),TIME,Energ(:,4));
surf(DifY);  % Space derivative field to show the soliton dynamics

fpu1 function
function dy=fpu1(t,y);
N=32;alpha=0.25;
D(N+1)=y(2)  -2*y(1)+alpha*((y(2)-y(1))\verb!^!2-y(1)\verb!^!2);D(1)=y(N+1);
D(2*N)=y(N-1)-2*y(N)+alpha*(y(N)\verb!^!2-(y(N)-y(N-1))\verb!^!2);D(N)=y(2*N);
for I=2:N-1,
D(N+I)=y(I+1)+y(I-1)-2*y(I)+alpha*((y(I+1)-y(I))\verb!^!2-(y(I)-y(I-1))\verb!^!2);
D(I)=y(N+I);
end
dy=D';