Instantaneous Pairing of Lyapunov Exponents in Chaotic Hamiltonian Dynamics and the 2017 Ian Snook Prize
William Graham Hoover Carol Griswold Hoover














The time-averaged Lyapunov exponents, { λ_i  }, support a mechanistic description of
the chaos generated in and by nonlinear dynamical systems. The exponents are ordered from largest
to smallest with the largest one describing the exponential growth rate of the ( small ) distance
between two neighboring phase-space trajectories.  Two exponents, λ_1 + λ_2, describe
the rate for areas defined by three nearby trajectories.  λ_1 + λ_2 + λ_3 is
the rate for volumes defined by four nearby trajectories, and so on.  Lyapunov exponents for
Hamiltonian systems are symmetric. The time-reversibility of the motion equations links the growth
and decay rates together in pairs.  This pairing provides  a more detailed explanation than
Liouville's for the conservation of phase volume in Hamiltonian mechanics. Although correct for
long-time averages, the dependence of trajectories on their past is responsible for the observed 
lack of detailed pairing for the instantaneous “local” exponents, { λ_i(t)  } . The
2017 Ian Snook Prizes will be awarded to the author(s) of an accessible and pedagogical discussion
of local Lyapunov instability in small systems. We desire that this discussion build on the two
nonlinear models described here, a double pendulum with Hooke's-Law links and a periodic chain of
Hooke's-Law particles tethered to their lattice sites.  The latter system is the ϕ^4 model
popularized by Aoki and Kusnezov.  A four-particle version is small enough for comprehensive
numerical work and large enough to illustrate ideas of general validity.




Instantaneous Pairing of Lyapunov Exponents in Chaotic Hamiltonian Dynamics and the 2017 Ian Snook Prizes ;

Short Running-Head Title for CMST :
 2017 Snook Prizes : How Lyapunov Exponents Pair


    
William Graham Hoover and Carol Griswold Hoover       

Ruby Valley Research Institute  

HC 60 Box 601    

Ruby Valley, NV 89833                             


    December 30, 2023
====================================================================================================================================================================================================





§ INTRODUCTION


The elucidation of Hamiltonian chaos and Lyapunov instability by Poincaré and Lorenz is familiar
textbook material. Models which capture aspects of complexity, the Logistic and Baker Maps, the
Lorenz attractor and the Mandelbrot Set, combine visual appeal with mechanistic understanding in
the bare minimum of spatial dimensions, two for maps and three for flows. Mechanical models with
only three- or four-dimensional phase spaces are simple enough that the entire phase space can be
explored exhaustively.“Small Systems” can augment our understanding of nature in terms of numerical
models by introducing more complexity. Just a few more degrees of freedom make an ergodic exhaustive
sampling impossible.  For the small systems we treat here we take on the more difficult task of
defining and analyzing the time-dependent convergence of “typical” trajectories. 

Chaos involves the exponential growth of perturbations.  Joseph Ford emphasized the consequence that
the number of digits required in the initial conditions is proportional to the time for which an
accurate solution is desired. Accordingly a “typical” nonexhaustive trajectory or history is the
best that we can do.  To go beyond the simplest models to those which elucidate macroscopic phenomena,
like phase transitions and the irreversibility described by the Second Law of Thermodynamics, we like
Terrell Hill's idea of small-system studies ( in the 1960s he wrote a prescient book,
Thermodynamics of Small Systems. ) In what follows we describe two small-system models which
are the foci of the Ian Snook Prize Problem for 2017.  These models are Hamiltonian, both with four
degrees of freedom so that their motions are described in eight-dimensional phase spaces.






 §.§ The Springy Pendulum and the Springy Double Pendulum

The double pendulum with rigid links is an excellent model for the table-top demonstration of chaos.  Bill
saw one in action at an all-day Stanford lecture given by James Yorke.  An even simpler mathematical model
for chaos can be obtained with a single pendulum. For chaos the single pendulum needs a spring rather than
a rigid link.  The single springy pendulum moves in a four-dimensional phase space, just as does the
double pendulum with rigid links.  Along with Harald Posch<cit.> we investigated mathematical models
for chaos based on chains of pendula, both rigid and springy.  We studied many-body instabilities by
characterizing the form of the detailed description of many-dimensional chaos, the Lyapunov spectrum.  We
considered two kinds of model Hamiltonians describing chains in a gravitational field :
[ 1 ] chains composed of particles
with equal masses, as in a physical length of chain; [ 2 ] chains in which only the bottom mass was affected
by gravity, as in a light chain supporting a heavy weight. Figure 1 shows five snapshots, equally spaced in
time, from a chaotic double-pendulum trajectory.  Initially the motionless chain was placed in the
horizontal configuration appearing at the top right of Figure 1.  If gravity affects only the lower of the
two masses ( as in the type-2 models supporting a heavy weight ) the  corresponding Hamiltonian is

    H = [  p_1^2 + p_2^2  ]/2 + (κ/2)[  (r_1-1)^2 + (r_12-1)^2  ] + y_2  .

where r_1 and r_12 are the lengths of the upper and lower springs. To enhance the coupling between
the springs and gravity we choose the force constant κ = 4 here.





 §.§ The Spectrum of Time-Averaged Lyapunov Exponents, { λ } 


The Lyapunov exponents making up the spectrum are conventionally numbered in the descending order of their
long-time-averaged values. We begin with the largest, λ_1 . λ_1 describes
the long-time-averaged rate at which the distance between the trajectories of two nearby phase-space
points increases. That rate,
λ_1 ≡⟨ λ_1(t)  ⟩≡⟨ (dlnδ/dt)  ⟩ ,
is necessarily positive in a chaotic system.  A more detailed description of rates of change of lengths and areas, and
volumes, and hypervolumes of dimensionality up to that of the phase space itself, leads to definitions
of additional Lyapunov exponents.  The next exponent, λ_2 , is needed to describe the rate at which
a typical phase-space  area, defined by three nearby points, increases ( or decreases ) with increasing time,
λ_1 + λ_2 ≡⟨ (dln A/dt)  ⟩ =
⟨ λ_1(t) + λ_2(t)  ⟩.
Again an average over a sufficiently long time for convergence is required.  Likewise the time-averaged
rate of change of a three-dimensional phase volume defined by four neighboring trajectories is λ_1
+ λ_2 + λ_3 . This sequence of rates and exponents continues for the rest of the spectrum.
There are D exponents for a D-dimensional phase-space description.



 §.§ Local and Global Lyapunov-Exponent “Pairing” for Hamiltonian Systems


The time-reversibility of Hamiltonian mechanics implies that all the rates of change change sign if the
direction of time is reversed.  This suggests, for instance, that all the exponents, { λ } and
{ λ(t)  } , are “paired”, with the rates forward in time opposite to those backward in
time.  This turns out to be “true” for the long-time-averaged exponents but could be “false” for the local
exponents.
Local exponents depend upon the recent past history of neighboring trajectories.  The global exponents,
which describe the growth and decay of the principal axes of comoving hyperellipsoids in phase space are paired,
though the time required to show this through numerical simulation can be long.  This exponent pairing is
the focus of the 2017 Snook Prize, as we detail in what follows.  There is a vast literature describing and
documenting the numerical evaluation and properties of Lyapunov spectra.  The theoretical treatments are
sometimes abstruse and lacking in numerical verification.  This year's Prize Problem seeks to help remedy this
situation.  The numerical foundation for the study of Lyapunov
exponents is an algorithm developed by Shimada and
Nagashima in Sapporo<cit.> and Benettin in Italy, along with his colleagues Galgani, Giorgilli,
and Strelcyn<cit.>, beginning in the late 1970s.  Google indicates hundreds of thousands of internet
hits for “Lyapunov Spectrum”.  We mention only a few other references<cit.> here. The
internet makes these and most of the rest readily available.




 §.§ The ϕ^4 Model for Chaos and Heat Conduction in Solids





Aoki and Kusnezov popularized the ϕ^4 model as a prototypical atomistic lattice-based model leading to
Fourier heat conduction<cit.>.  In addition to a nearest-neighbor Hooke's-Law potential the
model incorporates quartic tethers binding each particle to its own lattice site.  

Here we denote the
displacements of the particles from their sites as { q_i  } . In our one-dimensional
case the spacing between the lattice sites does not appear in the Hamiltonian or in the equations of
motion.  In numerical work it is convenient to choose the spacing equal to zero while setting the
particle masses, force constants for the pairs, and those for the tethers all equal to unity.  For a
four-particle problem in an eight-dimensional phase space the three-part Hamiltonian is :

    H = ∑_i=1^4 [  (p_i^2/2) + (q_i^4/4)  ] + ∑_4^springs (q_i,j^2/2)  .

The periodic boundary condition includes the spring linking particles 1 and 4 :

    q̈_1 = -q_1^3 + q_2 + q_4 - 2q_1  ;  q̈_4 = -q_4^3 + q_1 + q_3 - 2q_4  .

See Figure 2 for two ways of visualizing the periodic boundary conditions of the ϕ^4 chain.

The energy range over which chaos is observed in the ϕ^4 model includes about nine orders of
magnitude<cit.>.  The chaotic range for a four-body chain includes the two cases we discuss in
the present work, { E = 8,  288  ;  (E/N) = 2, 72  } . With both the springy pendulum and
the ϕ^4 models in mind we turn next to a description of their chaotic properties.



§ THE CHAOTIC DYNAMICS OF THE SPRINGY DOUBLE PENDULUM

 
Like most smoothly-differentiable Hamiltonian systems the double springy pendulum has infinitely
many periodic  or quasiperiodic phase-space solutions surrounded by a chaotic sea. Dynamics in the
sea is exponentially sensitive to perturbations.  The dynamics occurs in an eight-dimensional phase
space. Perturbations oriented along the trajectory or perpendicular to the energy surface, where there
is no longtime growth at all, give two zeroes, so that the maximum number of nonzero Lyapunov
exponents is six.

Each positive exponent is necessarily paired with its negative twin, with the two changing roles if
the direction of time is reversed. It is often stated that this time-reversible pairing links not
only the time-averaged rates of the dynamics, but also the “local” or “instantaneous” rates<cit.>. 
Because chaotic pendulum problems give different local exponents if Cartesian and polar coordinates
are used one might think that pairing could be hindered by using a mixture of these coordinates.
To check on this idea we considered a mixed-coordinate Hamiltonian for the model of Figure 1 with
polar coordinates for the “inside” Particle 1 :

    H = (1/2)[  p_r^2 + (p_θ/r)^2 + p_x^2 + p_y^2  ] + y_2
    + (κ/2)[  (r-1)^2 + (r_12-1)^2  ]  ;
                                                                                                                 

    r_12 = √( x_2^2 + y_2^2 + r_1^2 - 2r_1x_2sin(θ_1) + 2r_1y_2cos(θ_1)) ;  κ = 4  .

Formulating and solving the motion equations in mixed Cartesian and polar coordinates is an
intricate error-prone task.  It is useful first to solve the problem in Cartesian coordinates.
That solution then provides a check for the more complicated mixed-coordinate case. Energy
conservation is a nearly-infallible check of the programming.  We computed spectra of
Lyapunov exponents averaged over one billion fourth-order and one billion  fifth-order Runge-Kutta
timesteps, dt = 0.001 . This ensures that the numerical truncation errors of order (dt^5/120) or
(dt^6/720) are of the same order as the double-precision roundoff error. We chose the initial
condition of Figure 1 with both masses motionless at the support level, { x_1,y_1,x_2,y_2 } =
{ 1,0,2,0  } , so that the initial potential, kinetic, and total energies all vanished.
Only the outer Cartesian mass interacts with the gravitational field.

The simplest numerical method for obtaining Lyapunov spectra<cit.> is first to generate a
D-dimensional “reference trajectory” in the D-dimensional phase space.  Then a set of D
similar “offset” trajectories, an infinitesimal distance away, δ , are generated in
the same space with numerical offset vectors of length δ = 0.00001 or 0.000001.  While advancing the
resulting D(D+1) D-dimensional differential equations the local Lyapunov exponents are
obtained by “Gram-Schmidt” orthonormalization. This process rescales the vectors to their original
length and rotates all but the first of them in order to maintain their orthonormal arrangement.
The rescaling operation portion of the Gram-Schmidt process gives local values for the D Lyapunov
exponents :

    λ_i(t) ≡ (-1/dt)ln(δ_i^after/δ_i^before)  ;  λ_i ≡⟨ λ_i(t)  ⟩ .

For the type-2 double pendulum of Figure 1 the time-averaged Lyapunov spectrum is :

    { λ } = { +0.143, +0.076. +0.034, 0.000, 0.000, -0.034, -0.076, -0.143  } .

The rms fluctuations in these rates are typically orders of magnitude larger than the rates themselves. 
The uncertainty in the exponents as well as the differences between exponents using fourth-order
or fifth-order Runge-Kutta integrators with dt = 0.001 are both of order ± 0.001 . Our numerical work
shows that the pairing of the exponents is maintained if one of the pendula is described by polar
coordinates with the other pendulum Cartesian.  The local exponents are different but still paired.








§ CONVERGENCE AND ORDERING OF LOCAL LYAPUNOV EXPONENTS


The algorithm for generating the Lyapunov exponents<cit.> requires the ordering of D offset
vectors in the vicinity of a reference trajectory.  The first vector follows exactly the same motion
equations with the proviso that its length is constant.  The second vector, also of constant length,
is additionally required to remain orthogonal to the first so that the combination of the two gives
the rate of expansion or contraction of two-dimensional areas in the vicinity of the reference
trajectory. In general the nth offset vector satisfies n constraints in all, keeping its own
length constant while also maintaining its orthogonality to the preceding n-1 vectors.

Although the local rates { λ(t)  } associated with the vectors are necessarily ordered
when time-averaged  over a sufficiently long time to give the { λ } , this
ordering is regularly violated, locally, as Figures 3 and 4 show. Offhand one would
expect that increasing the Lyapunov exponents or decreasing the accuracy of the simulation would
lead to more rapid convergence of the ordering of the vectors.  For this reason we consider a
model which is as simple as possible, with a relatively large chaotic range, and is easy to simulate.
This ϕ^4 model, named for its quartic tethering potential, has proved particularly useful in
the simulation of heat flow.  We consider the equilibrium version of the model here, an isolated system.



§ THE DYNAMICS OF ONE-DIMENSIONAL PERIODIC Φ^4 MODELS






The simplest Lyapunov algorithm for the ϕ^4 model is exactly that used with the springy pendula.  We
follow D+1 trajectories in the D-dimensional phase space, rescaling them at every timestep
to obtain the complete spectrum of D = 8 instantaneous Lyapunov exponents. This phase-space integration
of nine trajectories, followed by Gram-Schmidt orthonormalization, can be modified
by using Lagrange multipliers to impose the eight constant-length constraints and the (1/2)(8·7) =
28 orthogonality constraints. A third approach, particularly simple to implement for the ϕ^4
model with its power-law equations of motion, is to linearize the motion equations so that the
offset vectors, rather than being small, can be taken as unit vectors in “tangent space”.  By
using separate integrators for the “reference trajectory” and for the eight unit vectors the
programming is at about the same level of difficulty as is that of the straightforward phase-space 
approach.  We implemented both approaches for the ϕ^4 problems and found good agreement for
the Lyapunov spectra at a visual level, even for calculations using a billion timesteps.  This is
because the reference trajectories for the phase-space and tangent-space algorithms are identical.



§ USEFUL INTEGRATION TECHNIQUES






Fourth-order and fifth-order Runge-Kutta integrators are particularly useful algorithms for small
systems.  First, these integrators  are easy to program. These integrators are also explicit, a real
simplification whenever a variable timestep is desirable. Their errors are typically opposite in sign.
For the simple harmonic oscillator the fourth-order energy decays while the fifth-order energy diverges.
By choosing a sufficiently small timestep, for which the two algorithms agree, one can be confident in
the accuracy of the trajectories.  Another useful technique is adapative integration : comparing solutions
with a single timestep dt to those from two successive half steps with (dt/2).  The timestep is then
adjusted up or down
by a factor of two whenever it is necessary to keep the root-mean-squared error in a prescribed band,
10^-12 >  error > 10^-14 for instance.<cit.>

At the expense of about a factor of fifty in computer time, FORTRAN makes it possible to carry out
quadruple-precision simulations with double-precision programming by changing the gnu compiler command :

    gfortran  - O  - o  xcode  code.f⟶ gfortran  - O  - o  xcode - freal - 8 - real - 16  code.f

Here the FORTRAN program is code.f and the executable is xcode .



§ THE 2017 IAN SNOOK PRIZE PROBLEM


The springy pendula and ϕ^4 problems detailed here show that “pairing” is typically present after
sufficient time, with that time sensitive to the largest Lyapunov exponent as well as to the initial
conditions. There are several features of these introductory problems that merit investigation :

[ 1 ] To what extent is there an unique chaotic sea ? Can the symmetry of the initial conditions limit the portion
of phase space visited when the dynamics is chaotic ?

[ 2 ] Within the ϕ^4 model's chaotic sea do
the time-averaged kinetic temperatures  { T_i = ⟨ p_i^2  ⟩ } , agree for all the particles ?
( If not, a thermal cycle applying heat and extracting work from the chain could be developed so as to
violate the Second Law.<cit.> )

[ 3 ] Is the pairing time simply related to the Lyapunov exponents and the chain length ?

[ 4 ] Is the accuracy of the pairing simply related to the accuracy of the integrator ?

The next and last question, which motivated this year's Prize Problem seems just a bit more difficult : [ 5 ]
Can relatively-simple autonomous Hamiltonian systems be devised for which long-time local pairing is
absent ? Our exploratory work has suggested that dynamical disturbances induced by collisions, with those
collisions separated by free flight, could lead to repeated violations of pairing<cit.>. On the other
hand Dettmann and Morriss have published a proof of pairing for isokinetic systems<cit.>.  A simple gas of
several diatomic or triatomic molecules is likely to be enough to settle that question.

The 2017 Ian Snook Prize will be awarded to the most interesting paper discussing and elucidating these
questions.  Entries should be submitted to Computational Methods in Science and Technology, cmst.eu, prior
to 1 January 2018. The Prize Award of 500 United States dollars sponsored by ourselves, and the Additional
Ian Snook Prize Award, also 500, will be awarded to the author(s) of the paper best addressing
this Prize Problem. 




§ ACKNOWLEDGMENTS


We are grateful to the Poznan Supercomputing and Networking Center for their support of these prizes
honoring our late Australian colleague Ian Snook (1945-2013).  We also appreciate useful comments,
suggestions, and very helpful numerical checks of our work furnished by Ken Aoki, Carl Dettmann, Clint
Sprott, Karl Travis, and Krzysztof Wojciechowski.  We particularly recommend Aoki's reference 10 for a
comprehensive study of the dynamics of one-dimensional equilibrium ϕ^4 systems.

99

b1  W. G. Hoover, C. G. Hoover, and H. A. Posch, “Lyapunov Instability of Pendula, Chains and
              Strings”, Physical Review A 41, 2999-3004 (1990).

b2  H. A. Posch, “Symmetry Properties of Orthogonal and Covariant Lyapunov Vectors and Their
              Exponents”, Journal of Physics A 46, 254006 (2013).

b3  I. Shimada and T. Nagashima, “A Numerical Approach to Ergodic Problems of Dissipative                 
              Dynamical Systems”, Progress of Theoretical Physics 61, 1605-1616 (1979).

b4  G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, “Lyapunov Characteristic                   
              Exponents for Smooth Dynamics Systems and for Hamiltonian Systems; a Method for Computing              
              All of Them, Parts I and II: Theory and Numerical Application”, Meccanica 15, 9-20
              and 21-30 (1980).

b5  J.-P Eckmann and D. Ruelle, “Ergodic Theory of Chaos and Strange Attactors”, Reviews of
              Modern Physics 57, 617-56 (1985).

b6  B. A. Bailey, “Local Lyapunov Exponents; Predictability Depends on Where You Are”, in
              Nonlinear Dynamics and Economics, W. A. Barnett, A. P. Kirman, and M. Salmon, editors
              (Cambridge University Press, 1996) pages 345-359.

b7  Hong-Liu Yang and Günter Radons, “Comparison of Covariant and Orthogonal Lyapunov Vectors”,
              Physical Review E 82, 046204 (2010) = arχiv:1008.1941.

b8  H. A. Posch and R. Hirschl, “Simulation of Billiards and of Hard-Body Fluids” in Hard           
              Ball Systems and the Lorentz Gas, Encyclopedia of the Mathematical Sciences 101, edited
              by D. Szász (Springer Verlag, Berlin, 2000), pages 269-310.

b9  K. Aoki and D. Kusnezov, “Nonequilibrium Statistical Mechanics of Classical Lattice ϕ^4
              Field Theory”, Annals of Physics 295, 50-80 (2002).

b10  K. Aoki, “Stable and Unstable Periodic Orbits in the One-Dimensional Lattice ϕ^4 Theory”,
              Physical Review E 94, 042209 (2016) .

b11  W. G. Hoover and K. Aoki, “Order and Chaos in the One-Dimensional ϕ^4 Model :
               N-Dependence and the Second Law of Thermodynamics”, Communications in Nonlinear Science
               and Numerical Simulation (in press, 2017) = arχiv 1605.07721.

b12  W. G. Hoover, J. C. Sprott, and C. G. Hoover, “Adaptive Runge-Kutta Integration for Stiff
               Systems: Comparing Nosé  and Nosé-Hoover Dynamics for the Harmonic Oscillator”,
               American Journal of Physics 84, 786-794 (2016).

b13  Wm. G. Hoover and C. G. Hoover, “Time-Symmetry Breaking in Hamiltonian Mechanics”,
               Computational Methods in Science and Technology 19, 77-87 (2013) = arχiv 1302.2533.

b14  Wm. G. Hoover and C. G. Hoover, “What is Liquid? Lyapunov Instability Reveals Symmetry-Breaking
               Irreversibilities Hidden Within Hamilton's Many-Body Equations of Motion”, Condensed Matter
               Physics 18, 1-13 (2015) = arχiv 1405.2485.

b15  C. P. Dettmann and G. P Morriss, “Proof of Lyapunov Exponent Pairing for Systems at Constant  
               Kinetic Energy”, Physical Review E 53, R5545-R5548 (1996).



