What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion
Typical Hamiltonian liquids display exponential "Lyapunov instability", also called "sensitive dependence on initial conditions". Although Hamilton's equations are thoroughly time-reversible, the forward and backward Lyapunov instabilities can differ, qualitatively. In numer...
Збережено в:
| Дата: | 2015 |
|---|---|
| Автори: | , |
| Формат: | Стаття |
| Мова: | English |
| Опубліковано: |
Інститут фізики конденсованих систем НАН України
2015
|
| Назва видання: | Condensed Matter Physics |
| Онлайн доступ: | http://dspace.nbuv.gov.ua/handle/123456789/153582 |
| Теги: |
Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Цитувати: | What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion / Wm.G. Hoover, C.G. Hoover // Condensed Matter Physics. — 2015. — Т. 18, № 1. — С. 13003:1-13. — Бібліогр.: 22 назв. — англ. |
Репозитарії
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
irk-123456789-153582 |
|---|---|
| record_format |
dspace |
| spelling |
irk-123456789-1535822019-06-15T01:29:57Z What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion Hoover, Wm.G. Hoover, C.G. Typical Hamiltonian liquids display exponential "Lyapunov instability", also called "sensitive dependence on initial conditions". Although Hamilton's equations are thoroughly time-reversible, the forward and backward Lyapunov instabilities can differ, qualitatively. In numerical work, the expected forward/backward pairing of Lyapunov exponents is also occasionally violated. To illustrate, we consider many-body inelastic collisions in two space dimensions. Two mirror-image colliding crystallites can either bounce, or not, giving rise to a single liquid drop, or to several smaller droplets, depending upon the initial kinetic energy and the interparticle forces. The difference between the forward and backward evolutionary instabilities of these problems can be correlated with dissipation and with the Second Law of Thermodynamics. Accordingly, these asymmetric stabilities of Hamilton's equations can provide an "Arrow of Time". We illustrate these facts for two small crystallites colliding so as to make a warm liquid. We use a specially-symmetrized form of Levesque and Verlet's bit-reversible Leapfrog integrator. We analyze trajectories over millions of collisions with several equally-spaced time reversals. 2015 Article What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion / Wm.G. Hoover, C.G. Hoover // Condensed Matter Physics. — 2015. — Т. 18, № 1. — С. 13003:1-13. — Бібліогр.: 22 назв. — англ. 1607-324X DOI:10.5488/CMP.18.13003 arXiv:1405.2485 PACS: 05.10.-a, 05.45.-a http://dspace.nbuv.gov.ua/handle/123456789/153582 en Condensed Matter Physics Інститут фізики конденсованих систем НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
English |
| description |
Typical Hamiltonian liquids display exponential "Lyapunov instability", also called "sensitive dependence on initial conditions". Although Hamilton's equations are thoroughly time-reversible, the forward and backward Lyapunov instabilities can differ, qualitatively. In numerical work, the expected forward/backward pairing of Lyapunov exponents is also occasionally violated. To illustrate, we consider many-body inelastic collisions in two space dimensions. Two mirror-image colliding crystallites can either bounce, or not, giving rise to a single liquid drop, or to several smaller droplets, depending upon the initial kinetic energy and the interparticle forces. The difference between the forward and backward evolutionary instabilities of these problems can be correlated with dissipation and with the Second Law of Thermodynamics. Accordingly, these asymmetric stabilities of Hamilton's equations can provide an "Arrow of Time". We illustrate these facts for two small crystallites colliding so as to make a warm liquid. We use a specially-symmetrized form of Levesque and Verlet's bit-reversible Leapfrog integrator. We analyze trajectories over millions of collisions with several equally-spaced time reversals. |
| format |
Article |
| author |
Hoover, Wm.G. Hoover, C.G. |
| spellingShingle |
Hoover, Wm.G. Hoover, C.G. What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion Condensed Matter Physics |
| author_facet |
Hoover, Wm.G. Hoover, C.G. |
| author_sort |
Hoover, Wm.G. |
| title |
What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion |
| title_short |
What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion |
| title_full |
What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion |
| title_fullStr |
What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion |
| title_full_unstemmed |
What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion |
| title_sort |
what is liquid? lyapunov instability reveals symmetry-breaking irreversibilities hidden within hamilton's many-body equations of motion |
| publisher |
Інститут фізики конденсованих систем НАН України |
| publishDate |
2015 |
| url |
http://dspace.nbuv.gov.ua/handle/123456789/153582 |
| citation_txt |
What is liquid? Lyapunov instability reveals symmetry-breaking irreversibilities hidden within Hamilton's many-body equations of motion / Wm.G. Hoover, C.G. Hoover // Condensed Matter Physics. — 2015. — Т. 18, № 1. — С. 13003:1-13. — Бібліогр.: 22 назв. — англ. |
| series |
Condensed Matter Physics |
| work_keys_str_mv |
AT hooverwmg whatisliquidlyapunovinstabilityrevealssymmetrybreakingirreversibilitieshiddenwithinhamiltonsmanybodyequationsofmotion AT hoovercg whatisliquidlyapunovinstabilityrevealssymmetrybreakingirreversibilitieshiddenwithinhamiltonsmanybodyequationsofmotion |
| first_indexed |
2025-07-14T05:02:06Z |
| last_indexed |
2025-07-14T05:02:06Z |
| _version_ |
1837597276448489472 |
| fulltext |
Condensed Matter Physics, 2015, Vol. 18, No 1, 13003: 1–13
DOI: 10.5488/CMP.18.13003
http://www.icmp.lviv.ua/journal
What is liquid? Lyapunov instability reveals
symmetry-breaking irreversibilities hidden within
Hamilton’s many-body equations of motion
Wm.G. Hoover, C.G. Hoover
Ruby Valley Research Institute, Highway Contract 60, Box 601, Ruby Valley, Nevada 89833
Received June 17, 2014, in final form November 25, 2014
Typical Hamiltonian liquids display exponential “Lyapunov instability”, also called “sensitive dependence on
initial conditions”. Although Hamilton’s equations are thoroughly time-reversible, the forward and backward
Lyapunov instabilities can differ, qualitatively. In numerical work, the expected forward/backward pairing of
Lyapunov exponents is also occasionally violated. To illustrate, we consider many-body inelastic collisions in
two space dimensions. Two mirror-image colliding crystallites can either bounce, or not, giving rise to a sin-
gle liquid drop, or to several smaller droplets, depending upon the initial kinetic energy and the interparticle
forces. The difference between the forward and backward evolutionary instabilities of these problems can be
correlated with dissipation and with the Second Law of Thermodynamics. Accordingly, these asymmetric sta-
bilities of Hamilton’s equations can provide an “Arrow of Time”. We illustrate these facts for two small crys-
tallites colliding so as to make a warm liquid. We use a specially-symmetrized form of Levesque and Verlet’s
bit-reversible Leapfrog integrator. We analyze trajectories over millions of collisions with several equally-spaced
time reversals.
Key words: Lyapunov instability, exponent pairing, chaotic dynamics, irreversibility
PACS: 05.10.-a, 05.45.-a
1. Introduction
Doug Henderson and John Barker helped to set the stage for our own Nonequilibrium developments
through their equilibriumwork on Thermodynamic Perturbation Theory [1]. This novel approach solved
the problem of calculating accurate liquid-state thermodynamics by approximating the structure of a liq-
uid with hard-sphere or soft-sphere pair distribution functions. In our 2004 contribution, we described
Nonequilibrium Molecular Dynamics, the offshoot of classical mechanics designed to treat mechanical
and thermal gradients according to generalizations of Gibbs’ statistical mechanics. Since then, we have
published a book summarizing these ideas [2]. Its successor, Simulation and Control of Chaotic Nonequi-
librium Systems has just been published by World Scientific Publishers in Singapore and released in
March 2015.
2. Lyapunov instability and Lyapunov spectra
A key finding of the nonequilibrium work was that steady-state distribution functions are singular
and fractal rather than Gibbsian and smooth, emphasizing the rarity of nonequilibrium states [3]. In ei-
ther of these cases, equilibrium or nonequilibrium, the necessary mixing in n-dimensional phase space
is facilitated by Lyapunov instability, the exponential growth of small perturbations. This instability is the
focus of our present work. Lyapunov instability is named for a Russian, a gifted and prolific mathemati-
cian with roots in Saint Petersburg, Alexander Lyapunov (1857–1918). Around 1979–1980, Shimada and
© Wm.G. Hoover, C.G. Hoover, 2015 13003-1
http://dx.doi.org/10.5488/CMP.18.13003
http://www.icmp.lviv.ua/journal
Wm.G. Hoover, C.G. Hoover
Nagashima [4] as well as Benettin, Galgani, Giorgilli, and Strelcyn [5] developed numerical methods for
evaluating the spectrum of all n Lyapunov exponents.
The spectrum describes the n-dimensional nature of Lyapunov instability in n-dimensional space.
The resulting orthogonal description of instabilities is much like the orthogonal description of vibrations
making up the solid-phase frequency distributions. The basic idea is to follow the motion of n satel-
lite trajectories in the neighborhood of an n-dimensional reference trajectory. The orthogonality of the
n-dimensional vectors separating the satellites from the reference can be enforced by Gram-Schmidt or-
thonormalization or by an equivalent set of constraining Lagrange multipliers [6]. When the motion is
Lyapunov unstable, the largest of the n exponents describing the instability — λ1 ≡ 〈 λ1(t) 〉 — the time-
averaged rate at which two nearby trajectories separate — can be determined from the growth rate of
the first Lyapunov vector δ1(t)≃ eλ1t .
Just as with the whole spectrum, this determination of λ1(t) can be done in either of two ways: (i)
rescale the distance between a satellite trajectory and the reference trajectory at each discrete timestep
or, (ii) constrain the length of the offset vector |δ1| with a Lagrange multiplier, λ1(t). The Lagrange mul-
tiplier approach [6] entails one multiplier for each of the n(n −1)/2 angles defined by a pair of vectors,
plus n additional multipliers for the lengths of the vectors. Since the constrained problem, along with all
of its Lagrange multipliers, is time-reversible, the complete set of multipliers going forward needs only to
change sign to maintain the orthonormality constraints in the reversed direction. Numerical work shows
that this reversibility is illusory (as is quite well known to the experts). The reversed set of vectors is
unstable, as we will see presently.
A second Lyapunov exponent, λ2(t), can be added to the first to describe the rate at which a two-
dimensional area in n space changes with time. Continuing this process to the third, fourth, . . .multiplier,
the sum of all n Lyapunov exponents gives the rate at which the n-dimensional hypervolume in the (q, p)
phase space changes with time:
⊗̇(t)/⊗ (t)≡
n
∑
1
λi (t) .
Figure 1 illustrates the relationship of the Lyapunov exponents to the orthogonalized vectors separating
the n satellite trajectories from the n-dimensional reference trajectory.
For simplicity, we develop all of our manybody models in two space dimensions, using particles of
unit mass. The dimensionality of the corresponding phase space is four times the number of particles,
n = 4N . There is a separate phase-space direction for each particle coordinate and momentum (velocity,
for particles of unit mass):
{q, p} ≡ {xi , ẋi , yi , ẏi }.
Typically, in many-body systems, the Lyapunov exponents are of the same order as the collision rate, and
the “spectra” of all the exponents resemble the Debye spectra of solid-state physics.
The time-reversibility of Hamilton’s equations of motion extends also to the time-reversibility of the
differential equations governing the Lyapunov vectors {δ(t)} and their associated exponents {λ(t)}. This
Figure 1. The Lyapunov exponents {λ1,λ2,λ3} are respectively the growth rates (δ̇/δ) of small orthogonal
vectors {δ1,δ2,δ3} in one-, two-, three-dimensional subspaces of n-dimensional phase space. If the vectors
are allowed to grow during each timestep, then they are rescaled in length by the Gram-Schmidt proce-
dure, which also maintains their orthogonality. If the vectors are instead constrained to a constant length
with Lagrange multipliers{λ1 ,λ2,λ3}, those multipliers are identical to the “local” (time-dependent) Lya-
punov exponents [6].
13003-2
What is liquid? Symmetry: Lyapunov + Hamilton
has two interesting consequences: (i) Directly from Hamilton’s equations of motion one would (naïvely)
expect that for every positive exponent and vector there is a time-reversed pair, with all the momenta
reversed:
{
+λi (+δq,+δp)t ←→−λn+1−i (+δq,−δp)t
}
.
Although this is true, it turns out that only one of the vectors in each pair is “stable”. The observed vectors
going forward in time are quite distinct from those going backward. (ii) For every Lyapunov vector of
the form δ= (δq,δp) there is also another paired orthogonal vector with an oppositely-signed Lyapunov
exponent and with the coordinate and velocity components of the original vector switched:
{
+λi (+δq,+δp)t ←→−λn+1−i (−δp,+δq)t
}
.
By permuting the components of the vectors in this way, orthogonality is guaranteed. This pairing is
mostly true. Typically, there is a simple relationship between vectors corresponding to Lyapunov expo-
nents with opposite signs. But we will see that this is not always the case. The occasional exceptions
brought about the present work.
We previously investigated the first and simpler of the two pairing ideasmentioned above, comparing
the forward and backward phase-space offset vectors δf
1 and δb
1 for planar shockwaves generated by
two colliding crystals. These mirror-image crystals moved toward each other in the x direction. In the
y direction, the boundary conditions were periodic [7]. See figure 2. Forward in time, the “important
particles”, making above-average contributions to δf
1, were concentrated within the hot shockedmaterial.
Backward in time, and with the very same configurations with opposite velocities, the above-average
contributions were less spatially concentrated [7].
In order to eliminate transients in such calculations, we hit upon the idea of cycling a bit-reversible
(exactly reversible, to the last bit) simulation forward and backward in time until the forward and back-
ward vectors δf
1 and δb
1 had converged to machine accuracy. Again, the important particles going forward
and backward in time were qualitatively different [8, 9]. Similar effects were found for binary collisions
of two crystallites in the absence of any spatial periodicity [10]. All these simulations, with or without
spatial or temporal periodicities, agreed in finding qualitative differences between the Lyapunov vectors
forward and backward in time. In reference [10], where the full Lyapunov spectrum for two colliding
37-particle hexagons was computed, the only vector and exponent pairing observed was quite imper-
fect. Though these calculations were bit-reversible, they spanned only tens of thousands of timesteps. We
consider much longer simulations in the present work.
Figure 2. Important particles, shown as open circles, in the collision of two periodic 240-particle
triangular-lattice crystallites. The forward and backward times are 2, 4, and 6, with reversal at 12. See
reference [7] for more details, especially figures 20 and 21 of that reference.
13003-3
Wm.G. Hoover, C.G. Hoover
Figure 3. (Color online) Important particles, shown in black, in the collision of two 400-particle balls
reversed at time t = 100. The snapshots correspond to three different times with identical particle co-
ordinates in both the forward and the backward directions of time. The δ vectors giving the forward
and backward values of λ1 are quite different in the two time directions. See figures 12 and 13 of refer-
ence [10].
The leading vectors going forward in time emphasize the leading edges of the crystals, where the col-
lision is taking place. In the absence of periodic boundary conditions, the leading vectors going backward
in time (and we will soon describe the best way to go backward) instead emphasize the “necking” region,
where the compound liquid drop formed by the colliding crystals relives its past history as two separate
bodies. See figure 3 for the collision of two 400-particle crystalline balls. In that figure, the particles mak-
ing above-average contributions to the largest Lyapunov exponent going both forward and backward in
time are emphasized.
In the present work, we reduce the intricacies of the Lyapunov spectrum and the instabilities it de-
scribes by considering smaller systems for longer times. These are all Hamiltonian systems with two
crystallites colliding to form one or more fragments. These smaller simpler systems make it possible to
study Lyapunov instability and the pairing of vectors with greater precision.
Figure 4 shows sample 74-particle snapshots for two different initial velocities. At relatively low ve-
locities, the colliding hexagons can bounce or coalesce. At higher velocities, several smaller crystallites or
Figure 4. Pairs of 37-particle crystallites and their collision products. The initial condition for all of these
simulations corresponds to the upper left-hand illustration. At the low velocity of ±0.10, the result is
the “liquid” ball shown at the lower left. At the higher velocity of ±0.50, the results for bit-reversible
simulations with and without symmetrization (described in section 6) are shown at the right where the
74 particles have separated into six fragments. The two bit-reversible simulations differ as a result of the
Lyapunov-unstable amplification of computer roundoff errors, as discussed in section 6.
13003-4
What is liquid? Symmetry: Lyapunov + Hamilton
drops are formed. To simplify both the dynamics and the analysis for corresponding pairs of particles in
the two 37-particle hexagons, we choose inversion-symmetric initial conditions:
{
xLeft + xRight = yLeft + yRight = 0; ẋLeft + ẋRight = ẏLeft + ẏRight = 0
}
.
In order to propagate the particles reversibly, we use Levesque and Verlet’s bit-reversible algorithm [11]
which we detail in the following section.
3. Levesque-Verlet bit-reversible simulations
The study of Lyapunov instabilities requires special numerical methods. Since our aim here is to
compare the stabilities of forward and backward motion equations for millions of timesteps, we begin
with Levesque and Verlet’s observation that the “Leapfrog” algorithm for solving Newton’s equations of
motion can be made precisely time-reversible by restricting the particle coordinates to (large) integer
values:
{
qt+d t ≡ 2qt −qt−d t +
[
(dt 2/m)F (qt )
]
int
}
←→
{
qt−d t ≡ 2qt −qt+d t +
[
(dt 2/m)F (qt )
]
int
}
.
Apparently, this identity ensures reversibility. All that has to be done is to compute and round off the force
terms, as indicated by the brackets [. . .]int, in precisely the same way whether going forward in time, to
t+dt , or backward, to t−dt . By using “long” 16-digit integers, a satisfactory precision can be obtained. An
alternative is to store the billions of trajectory coordinates describing the forward trajectory. With either
method, a strictly “bit-reversible” reference trajectory can be generated forward in time and can then
be used in reversed order to describe its backward relative. Since the bit-reversible Leapfrog technique
relies on a “conservative” one-to-one mapping of successive pairs of configurations, it cannot be applied
to “dissipative” motion equations. Dissipative motions cause the phase volume to shrink. Ultimately, that
chaotic shrinking generates fractal strange attractors [2, 3, 7].
The application of the bit-reversible algorithm to the computation of Lyapunov spectra was pioneered
by M. Romero-Bastida, D. Pazó, J.M. López, and M.A. Rodríguez [12]. With a strictly reversible reference
trajectory, the motion of the n nearby satellite trajectories can then be generated with straightforward
Runge-Kutta integration. The corresponding numerical method is described in section 4. These ideas are
then applied to the idealized liquid Hamiltonian described in section 5. The results, for typical collisions,
are detailed in section 6. The final section 7 sums up the connection of these differences to irreversibility,
as described by the Second Law of Thermodynamics.
4. Combining Newtonian and Hamiltonian mechanics
In our molecular dynamics work we combine the Newtonian and the Hamiltonian forms of mechan-
ics. The coordinate-based Newtonian Leapfrog algorithm advances pairs of coordinate configurations
{qt , qt±dt } while Hamilton’s first-order motion equations advance from one {qt , pt } phase point to the
next {qt+dt , pt+dt }:
{
q̈ = F (q)
}
versus
{
q̇ = p; ṗ = F (q)
}
.
Combining the two forms of mechanics requires a definition of momentum based on the coordinate infor-
mation generated by the Leapfrog algorithm. The unimaginative first-order choice, with errors
(1/2)
...
q dt 2,
pt ≃
[
qt+dt −qt−dt
]
/(2dt),
can and should be improved upon by using instead the third-order definition [10],
pt ≡ (4/3)
[
qt+dt −qt−dt
]
/(2dt)− (1/3)
[
qt+2dt −qt−2dt
]
/(4dt).
The formal error in this last definition is −(1/30)
.....
q dt 4.
13003-5
Wm.G. Hoover, C.G. Hoover
Figure 5. Pairs of 7-particle crystallites and their collision products for two choices of the initial veloc-
ities, ±0.1 (above) and ±0.5 (below). Notice that in the higher-speed collision, the crystallites exchange
one particle before separating. The inversion symmetry has been maintained by symmetrization at each
timestep but the original fourfold symmetry of the hexagons has been lost. The symmetrization process
is described in section 6.
5. Numerical liquid models for the collision process
To minimize numerical errors, it is useful to choose very smooth force laws. As demonstration prob-
lems we consider here collisions of two hexagonal crystallites. Snapshots appear in figure 4 and 5 for
collisions of both 37-particle and 7-particle crystallites. We follow reference [10] and use a many-body
attractive binding energy, (1/2)(ρi −1)2 for each Particle i . This potential, when differentiated gives the
force on Particle i due to Particle j as (2−ρi −ρ j )∇i wi j . Each Particle’s personal density is computed
using Lucy’s weight function [2]. All particles { j } (including i ) within a distance h make a contribution to
the density there:
ρi ≡
14 or 74
∑
j=1
w(r < h), w = (5/h2π)(1+3z)(1− z)3, z ≡ r /h.
The normalization constant (5/h2π) is chosen so that the integral of the weight function is unity:
h
∫
0
w(r )2πr dr ≡ 1.
For the model systems discussed here, with N = 2×37 = 74 or N = 2×7 = 14, we have used h = 3.00 and
h = 3.50, respectively. A typical nearest-neighbor distance in all of these problems is of order unity.
In order to minimize the occurrence of very closely-spaced pairs of particles, it is expedient to include
a short-ranged repulsive potential in the Hamiltonian. For our demonstration problems we choose a very
smooth “soft-disk” short-ranged pair potential:
φpair(r < 1) = (1− r 2)4.
With these specially smooth choices for the attractive binding potential and the repulsive pair potential
(and with our third-order definition of the momentum), the energy throughout a run remains constant
with six-figure accuracy.
The satellite solutions for each timestep are launched basing the offset vectors on the current leapfrog
values of the coordinates and their associated momenta {q, p}. Then, a fourth- or fifth-order Runge-Kutta
integrator provides close to machine accuracy in incrementing the satellite motions with a timestep
of 0.001.
13003-6
What is liquid? Symmetry: Lyapunov + Hamilton
6. Crystallite collisions and their time-reversed twins
Continuing to pursue simplicity we choose initial conditions with symmetric coordinates and mo-
menta for the two colliding crystallites. Sample geometries are shown in figures 4 and 5.
{
qi =−qN+1−i ; pi =−pN+1−i
}
, i = 1. . . N .
Figure 6 shows the evolution of the summed-up largest and smallest Lyapunov exponents, λ1(t)+λ56(t),
both forward and backward in time, but for 14 particles rather than 74 or 800 and for one million
timesteps, half a million forward and half a million reversed. The initial orientations of the offset vec-
tors,
{
δ≡ (q, p)sat − (q, p)ref
}
, can either be chosen “randomly” or as rows of a unit matrix. A convenient
length for the vectors in these numerical simulations is 0.0001. For our collision problems it takes about
50 000 timesteps (with dt = 0.001) for the time-dependent vectors to converge stably to a set independent
of the initial conditions.
Figure 6. Deviation from exponent pairing going forward in time (left of center) and backward (right
of center). The Lyapunov sum, λ1 +λ56 is plotted in both time directions for a 7+ 7 particle collision
with velocities ±0.10. The data are plotted for the third repetition, for iterations between 3 million and 4
million, where dt =±0.001. The deviation from pairing going forward in time is at its maximum at a time
near 243, corresponding to iteration numbers from 3238000 to 3247000. The insets show details for time
windows of width 9. In the reversed motionwith dt negative and at time 243, the pairing is nearly perfect,
invisible on the main plot. The reversed-time deviation from pairing is three orders of magnitude smaller
than the forward-time deviation. The huge backward spike near iteration 3510000 corresponds to the
conversion of the forward set of δ vectors into the backward set following time reflection at iteration
3500000 where t = 500.
We reverse the direction of time periodically, by changing the sign of dt every half million timesteps,
with nineteen changes in all in the course of a ten-million-timestep run. The reference trajectories for
the forward and backward segments are each repeated ten times and are identical to machine accuracy
using the Levesque-Verlet bit-reversible integrator. The satellite trajectories rotate about the reference
trajectory, constrained to remain orthogonal and of constant length. After the first repeated forward
and backward segments, these trajectories hardly change in subsequent segments. The only significant
changes occur near the middle of the collision process, forward in time. For typical details see the inset
of figure 6.
One would expect that the left-right and up-down symmetries imposed on the initial conditions would
persist throughout any bit-reversible calculation. Wrong! Because the order in which the summed up
forces give the total force on a particle are not necessarily symmetric, the sums can, and eventually do,
differ in the last decimal, with the resulting difference amplifying exponentially in time.
Take a toy-model example in which Particles 1, 2, 3, and 4 are arranged in order on a line with the
pair interactions added up in the usual (i < j ) order. The force on the leftmost particle, Particle 1, is a sum
of first, second, and third neighbor forces, in that order. The force on the rightmost particle, Particle 4,
13003-7
Wm.G. Hoover, C.G. Hoover
is instead the sum of third, second, and first neighbor forces, the opposite order. On a finite-precision
machine, the totals are likely different. Since the average force is typically close to zero, it is often the
case that significant figures are lost in the process of summing the forces on particles.
It is annoying to learn that an initially perfectly symmetrized bit-reversible problem loses its sym-
metry in a macroscopic way after approximately 50000 timesteps, about the same time as is required
for Lyapunov exponent pairs to pair up. Just as in the simple four-particle toy-model example, the sym-
metrized problem loses its symmetry due to the capricious nature of the ordering of the force sums.
Lyapunov instability can quickly move these tiny errors from the last decimal place to the first, mag-
nifying them by a factor of 1016, and so producing configurations which visibly lack symmetry. For an
inadvertent example see the asymmetry shown in the bit-reversible figure 5 of reference [10]. In our
“symmetric” demonstration problems here, the inversion symmetry can be maintained by the expedient
of symmetrizing the forces on each member of the particle pairs at every timestep:
Fi = [Fi −Fi+7]/2; Fi+7 = [Fi+7 −Fi ]/2,
likewise taking care to compute symmetrized densities for each particle:
ρi = [ρi +ρi+7]/2; ρi+7 = [ρi+7 +ρi ]/2.
The pairing of local (instantaneous) Lyapunov exponents has beenmuch discussed [13–18], mostly for
very small systems with only a few degrees of freedom. Our 74-particle simulations indicated exponent
pairing most of the time. But the complexity of the Gram-Schmidt calculations, orthogonalizing 296 vec-
tors at each timestep, left unclear the reason for the occasional loss of pairing. A numerical source seemed
likely because the largest Lyapunov exponent was closely reproduced (visually there was no detectable
change) from one million-timestep sequence to the next while variations in the smallest exponent, λ296
were visible and wholly responsible for the nonzero values of the sum, λ1(t)+λ296(t). To avoid any un-
certainty, we chose to concentrate on the simpler 14-particle problems with their reduced demands on
the Gram-Schmidt orthonormalization step.
The moment of collision for the problem that we choose to study in detail is shown in figure 7. Orig-
inally, at time zero, the center-to-center separation of the two crystallites was 50 with all particle veloc-
ities (0,+0.1) in the leftmost hexagon and (0,−0.1) in the rightmost hexagon. The least-energy nearest-
neighbor separation is 0.861121270463, which minimizes the potential energy, and there is no thermal
motion. Without any interaction, the two crystallites would have overlapped perfectly at a time of 250.
As shown in the figure, the collision actually begins at time [50−3.5−2×0.8611]/0.2 ≃ 223.9. By a time
of 250, the collision has converted two cold crystallites to a single warm drop. All of the results described
here for this collision use an offset vector length |δ| = 0.0001 and a timestep dt = 0.001.
Figure 6 showed the gross features of the sum of the first and last Lyapunov exponents over the fourth
interval of one million steps. On the scale shown, all of the intervals from the second through the tenth
Figure 7. Two 7-particle cold hexagons at the moment of collision, where the separation of their leading
edges is h = 3.5. The original center-to-center distance at time zero was 50. This system was followed
for 10000000 timesteps alternating segments with 500000 forward steps followed by 500000 backward
steps in order to ensure the convergence of the offset vectors both forward and backward in time.
13003-8
What is liquid? Symmetry: Lyapunov + Hamilton
3rd
collision
230 < time < 250
O O
3rd
collision
230 < time < 250
O O
4th
collision
O O
4th
collision
O O
5th
collision
O O
5th
collision
O O
6th
collision
O O
6th
collision
O O
Figure 8. Transition region details are illustrated for four different collisions, separated by one million
timesteps. We show the time-dependence of the first and last Lyapunov exponents, +λ1(t) ≃−λ56(t) for
the inelastic collision of two 7-particle hexagons with opposite initial velocities (±0.1,0.0). Apart from the
transition region between the two open circles, the pairing of the exponents is nearly perfect.
behaved identically, with the Lyapunov exponents paired during the reversed motion but consistently
undergoing an episode of nonpairing, going forward in time, at a time of 235 (corresponding to 235000,
1235000, 2235000, . . . timesteps). A detailed investigation shows that the lack of pairing occurs during
a shearing motion of the central rows of particles relative to the upper and lower ones. It is this same
shearing motion which was associated with the dynamics of hard disks at melting in 1963 [19]. For most
of the simulation, the first and the last exponents nearly sum to zero, λ1 +λ56 ≃ 0. The second half of the
run retraces the configurations beginning at the maximum time of 500. Figure 6 shows that λ1 reverses
visually to a time of 488 where the set of δ vectors begins its change from the unstable reversed forward
vectors to the stable backward ones. For the remainder of the run, back to the initial configuration at time
zero, no further disturbances to pairing are observed. Indeed, as the figure shows, the pairing is mostly
quite good (as the sum is very close to zero) forward in time too, in the interval 0 < t < 500, except where
the collisional effects are a maximum near t = 243. It is interesting to see that both before and after
the collision, the pairing is nearly perfect. Within the collision, there are some brief but quite significant
differences.
In figure 8 we show the neighborhood of the transition region before and after the time of maximum
shear. Four different collisions, separated by one million timesteps from their predecessors and succes-
sors, are detailed there. Between the open circles, the lack of pairing is noticeable, and changes from one
collision to the next. Outside the transition region and throughout, the reversed collision pairing is closely
satisfied. Let us next analyze the contributions of the individual particles to the collision process within
the transition region.
Figure 9 shows the individual-particle details of the change from the original hexagonal crystallite
shape during the collision. In this low-speed example, two of the particles in each hexagon shear over
and under their neighbors in adjacent rows.
For each of the particles we can quantify its importance to the dynamical instability by computing the
set of all fourteen particle amplitudes,
{
δx2
i +δy2
i +δẋ2
i +δẏ2
i
}
−→
N
∑
δ2
i ≡ 1.
Figure 10 shows the time-dependent amplitudes for Particles 3 and 10 during the maximum shear
period with time increasing: 240 < t < 245. The second and third repetitions of this time period are
both plotted here but nearly coincide. Among the fourteen particles, Particle 3 (and its inverted image
Particle 10) stand out with amplitudes near the maximum possible (1/2) throughout the shearing motion.
Evidently, the lack of pairing during the collision (as shown in figure 8) is linked to the change in stability
of the motions of these two particles as they pass by their neighbors, indicating the importance of shear
to the stability of fluid deformation.
13003-9
Wm.G. Hoover, C.G. Hoover
-2
-1
0
1
2
-3 < x < 3
Before and After Shear
Times = 241 and 243
Figure 9. Snapshots at times of 241 and 243, before (filled circles) and after (open circles) the cooperative
shearing motion of the central row of particles.
Fifty years ago we helped point out that the hard-disk solid melts when it becomes possible for shear
of exactly this same kind to occur [19]. See figure 11. This observation came from watching supercom-
puter movies of hard-disk molecular dynamics at the Radiation Laboratory in Livermore. The number of
disks was 870 and the boundary conditions were periodic. Now this same mechanism can be seen any-
where in the world on a laptop computer. It is interesting that despite all the change in computers and
computation, the basic principles ofmechanics, and the conclusions emerging from them, are unchanged.
7. Conclusions
Many facets of classical mechanics still remain to be illuminated. Some of them directly concern liq-
uids and dynamic instability. Today, computational advances make it possible to evaluate Lyapunov spec-
tra for systems with a few thousand exponents. Just the first of them, λ1(t), when compared to the last,
λn(t), shows that the melting mechanism identified from computer movies fifty years ago is still active
today in the irreversibly unstable collisional processes responsible for the Second Law of Thermodynam-
ics. The irreversibility of plastic flow in solids corresponds to the hysteresis of dislocation motions. In
that flow process: stored energy −→ heat. The shearing motion of dislocations is a close relative of the
instability found in the present work. Ergodicity and the packing of hard disks has highlighted the im-
0
0.1
0.2
0.3
0.4
0.5
240 < time < 245
Amplitude for Particle 3
( as time Increases )
Figure 10. The contributions of the two “important” particles (the two filled circles closest to x = 0 in
figure 9) are identical and together account for about 90% of the amplitude of δf1. The variation is shown
in the time window of maximum shear, 240 < t < 245. Data for the second and third repetitions of the
window are nearly identical and are shown as solid and dashed lines, respectively. The initial velocities
of the two 7-particle hexagons are ẋ =±0.1.
13003-10
What is liquid? Symmetry: Lyapunov + Hamilton
Figure 11. Figures from fifty years ago showing that the mechanism for hard-disk melting is the same as
that found for the irreversible aspect of an inelastic crystal collision. The 2-particle van der Waals loop
at the right of the figure comes from the cooperative shearing motion shown at the left for a periodic
two-disk system. Reprinted from reference [19].
portance of cooperative shearing to transitional behavior [19–21]. The “free volumes” accessible to liquid
particles were of interest in Doug Henderson and John Barker’s work in 1976 [1], in our work [20, 21],
and in others’ [22]. Themotivation for using periodic rather than rigid boundaries inmolecular dynamics
simulations becomes very clear on examining the “zoo” of hard-disk structures which can be “locked in”
in the rigid case [22] where shear flow is prevented.
Nonequilibrium steady states generated with time-reversible thermostats [3] are relatively simple to
“understand”. When the long-time-averaged change in phase-space density is nonzero, the only possibil-
ity consistent with stability (associated with a bounded phase-space distribution) is a strange attractor,
with the comoving change of density positive, 〈(d f /dt〉)> 0, while the fractal density itself is unchanged,
(∂ f /∂t) ≡ 0. Familiarity with this seemingly paradoxical state of affairs has led to its acceptance. A similar
understanding of Hamiltonian irreversibility is not yet so clear.
An illuminating fringe benefit of our efforts was discovering that bit-reversible dynamics may not
retain the symmetry of its initial conditions. Whether or not symmetry is maintained is dependent upon
the manner in which the particle forces are computed. Retention cannot be guaranteed unless the forces
are summed up in an explicitly symmetric manner. Long-time reversibility studies need to make this
choice explicit.
We have emphasized that the local growth and coalescence rates in phase space are identical going
forward or backward in time [7]. It is only through the influence of the “past” that there can be a lack
of symmetry in the local Lyapunov spectra. This hidden source of irreversibility is well worth mining
as it is the most obvious property distinguishing one direction of time travel from the other. The lack of
symmetry between the past and the future leaves its signature in the vectors and the exponents iden-
tifying important particles. The localization of these particles, which jumps from spot to spot in larger
systems [20], is an explicitly time-irreversible property, a tantalizing hint toward the understanding of
Hamiltonian irreversibility in terms of dynamical instabilities.
Though the pairing of exponents seems to be a widely-accepted consequence of Hamiltonianmechan-
ics, it is clear that the present results do violate pairing in the forward time direction. In the reversed
13003-11
Wm.G. Hoover, C.G. Hoover
direction, fluctuations in pairing are at least two orders of magnitude smaller. Why isn’t pairing violated
in the backward direction? Evidently, there is still more work to be done.
Acknowledgment
We thank the anonymous referee for several useful suggestions. We address most of them here, in
order to clarify some of the underlying concepts and details of our work.
The referee wished for evidence as to the Debye-like nature of manybody Lyapunov spectra. This
became clear for manybody simulations during the 1980s. In addition to figure 1 of reference [6], two
nonequilibrium spectra appear as figure 1 in Wm.G. Hoover, “The Statistical Thermodynamics of Steady
States”, Physics Letters A, 1999, 255, 37–41.
That paper also illustrates the fact that low-temperature “stochastic” thermostats, represented by a
drag force −(p/τ), can lead to fractal distributions. This fractal nature, as judged by the Kaplan-Yorke
dimension, occurs because the phase-space offset vectors rotate very rapidly. Vectors linking a reference
trajectory to an orthogonal array of rapidly-rotating satellite trajectories mix the dimensionality loss as-
sociated with drag among all phase-space directions.
It needs to be emphasized that the offset vectors are in phase space rather than tangent space. The
insensitivity of our results to the length of these vectors and to the timestep was carefully checked.
The inversion symmetry of the balls and crystallites in figures 3–5 implies that the top of the leftmost
projectile corresponds to the bottom of its rightmost image. This inversion symmetry is particularly no-
ticeable in the “important particles” shown here in figure 3.
As the referee kindly points out, it has been stated that Hamiltonian offset vectors are paired so that
the sum of the first and last local Lyapunov exponents should vanish exactly. Here, we simply state our
finding that this pairing is sometimes violated, not only initially, but also at long times.
The reason why we have used a collective density-dependent attractive potential for these collision
problems is that this reduces the yield strength so that the flow canmore easily be observed and analyzed.
References
1. Barker J.A., Henderson D., Rev. Mod. Phys., 1976, 48, 587; doi:10.1103/RevModPhys.48.587.
2. Hoover Wm.G., Hoover C.G., Time Reversibility, Computer Simulation, Algorithms, Chaos, World Scientific, Sin-
gapore, 2012.
3. Holian B.L., Hoover Wm.G., Posch H.A., Phys. Rev. Lett., 1987, 59, 10; doi:10.1103/PhysRevLett.59.10.
4. Shimada I., Nagashima T., Prog. Theor. Phys., 1979, 61, 1605; doi:10.1143/PTP.61.1605.
5. Benettin G., Galgani L., Giorgilli A., Strelcyn J.M., Meccanica, 1980, 15, 9; doi:10.1007/BF02128236.
6. Hoover W.G., Posch H.A., Phys. Lett. A, 1985, 113, 82; doi:10.1016/0375-9601(85)90659-0.
7. Hoover Wm.G., Hoover C.G., AIP Conf. Proc., 2011, 1332, 23; doi:10.1063/1.3569486; Preprint arXiv:1008.4947,
2010.
8. Hoover Wm.G., Hoover C.G., Computational Methods in Science and Technology, 2013, 19, 69;
doi:10.12921/cmst.2013.19.02.69-75; Preprint arXiv:1112.5491, 2011.
9. Hoover Wm.G., Hoover C.G., Computational Methods in Science and Technology, 2013, 19, e1;
doi:10.12921/cmst.2013.19.03.e1.
10. Hoover Wm.G., Hoover C.G., Computational Methods in Science and Technology, 2013, 19, 77;
doi:10.12921/cmst.2013.19.02.77-87; Preprint arXiv:1302.2533, 2013.
11. Levesque D., Verlet L., J. Stat. Phys., 1993, 72, 519; doi:10.1007/BF01048022.
12. Romero-Bastida M., Pazó D., López J.M., Rodríguez M.A., Phys. Rev. E, 2010, 82, 036205;
doi:10.1103/PhysRevE.82.036205.
13. Hoover Wm.G., Hoover C.G., Posch H.A., Phys. Rev. A, 1990, 41, 2999; doi:10.1103/PhysRevA.41.2999 [see sec-
tion III].
14. Ramaswamy R., Eur. Phys. J. B, 2002, 29, 339; doi:10.1140/epjb/e2002-00313-8 [see section 3].
15. Bosetti H., Posch H.A., Dellago Ch., Hoover Wm.G., Phys. Rev. E, 2010, 82, 046218;
doi:10.1103/PhysRevE.82.046218; Preprint arXiv:1004.4473, 2010.
13003-12
http://dx.doi.org/10.1103/RevModPhys.48.587
http://dx.doi.org/10.1103/PhysRevLett.59.10
http://dx.doi.org/10.1143/PTP.61.1605
http://dx.doi.org/10.1007/BF02128236
http://dx.doi.org/10.1016/0375-9601(85)90659-0
http://dx.doi.org/10.1063/1.3569486
http://arxiv.org/abs/1008.4947
http://dx.doi.org/10.12921/cmst.2013.19.02.69-75
http://arxiv.org/abs/1112.5491
http://dx.doi.org/10.12921/cmst.2013.19.03.e1
http://dx.doi.org/10.12921/cmst.2013.19.02.77-87
http://arxiv.org/abs/1302.2533
http://dx.doi.org/10.1007/BF01048022
http://dx.doi.org/10.1103/PhysRevE.82.036205
http://dx.doi.org/10.1103/PhysRevA.41.2999
http://dx.doi.org/10.1140/epjb/e2002-00313-8
http://dx.doi.org/10.1103/PhysRevE.82.046218
http://arxiv.org/abs/1004.4473
What is liquid? Symmetry: Lyapunov + Hamilton
16. Hoover Wm.G., Hoover C.G., In: Proceedings of the 40th Conference “Advanced Problems in Mechanics” (St. Pe-
tersburg, 2012), Institute of Problems of Mechanical Engineering, St. Petersburg, 2012; Preprint arXiv:1205.1276,
2012.
17. Posch H.A., J. Phys. A: Math. Theor., 2013, 46, 254006; doi:10.1088/1751-8113/46/25/254006; Preprint
arXiv:1107.4032, 2012.
18. Waldner F., Hoover Wm.G., Hoover C.G., Chaos, Solitons & Fractals, 2014, 60, 68; doi:10.1016/j.chaos.2014.01.006.
19. Alder B.J., Hoover W.G., Wainwright T.E., Phys. Rev. Lett., 1963, 11, 241; doi:10.1103/PhysRevLett.11.241.
20. Hoover Wm.G., Boercker K., Posch H.A., Phys. Rev. E, 1998, 57, 3911; doi:10.1103/PhysRevE.57.3911.
21. Hoover W.G., Hoover N.E., Hanson K., J. Chem. Phys., 1979, 70, 1837; doi:10.1063/1.437660.
22. Carlsson G., Gorham J., Kahle M., Mason J., Phys. Rev. E, 2012, 85, 011303; doi:10.1103/PhysRevE.85.011303;
Preprint arXiv:1108.5719, 2011.
Що таке рiдина? Нестiйкiсть Ляпунова розкриває
необоротнiсть порушення симетрiї, яка прихована у
багаточастинкових рiвняннях руху Гамiльтона
В.Г. Гувер, К.Г. Гувер
Дослiдницький iнститут м. Рубi Валей, Рубi Валей, Невада 89833, США
Типовi рiдини Гамiльтона проявляють експоненцiальну “нестiйкiсть Ляпунова”, яку також називають “за-
лежнiстю, чутливою до початкових умов”. Хоча рiвняння Гамiльтона є цiлком оборотнi по часу, прямi i
зворотнi необоротностi Ляпунова можуть якiсно рiзнитися. При числових обрахунках очiкуване спарю-
вання вперед/назад експонент Ляпунова також iнколи порушується. Наприклад, розглянемо непружнi
зiткнення багатьох тiл у двовимiрному просторi. Два дзеркально вiдображенi зiштовхуванi кристалiти
можуть пiдстрибнути або нi, в результатi чого виникає єдина крапелька рiдини, або декiлька менших
краплин в залежностi вiд початкової кiнетичної енегiї та мiжчастинкових сил. Вiдмiннiсть мiж прямими
i зворотнiми еволюцiйними нестiйкостями може корелювати з дисипацiєю та з другим законом термо-
динамiки. Отже, цi асиметричнi стiйкостi рiвнянь Гамiльтона можуть забезпечити “Стрiлу Часу”. У данiй
статтi цi факти проiлюстровано для двох невеличких кристалiтiв, якi зiштовхуються з тим, щоб утворити
теплу рiдину. Тут використано спецiально симетризовану форму бiт-оборотнього iнтегратора “Лiп-фрог”
Левека i Верле. Проаналiзовано траєкторiї понад мiльйони зiткнень при декiлькох рiвновiддалених часо-
вих оберненнях.
Ключовi слова: нестiйкiсть Ляпунова, експонентне спарювання, хаотична динамiка, необоротнiсть
13003-13
http://arxiv.org/abs/1205.1276
http://dx.doi.org/10.1088/1751-8113/46/25/254006
http://arxiv.org/abs/1107.4032
http://dx.doi.org/10.1016/j.chaos.2014.01.006
http://dx.doi.org/10.1103/PhysRevLett.11.241
http://dx.doi.org/10.1103/PhysRevE.57.3911
http://dx.doi.org/10.1063/1.437660
http://dx.doi.org/10.1103/PhysRevE.85.011303
http://arxiv.org/abs/1108.5719
Introduction
Lyapunov instability and Lyapunov spectra
Levesque-Verlet bit-reversible simulations
Combining Newtonian and Hamiltonian mechanics
Numerical liquid models for the collision process
Crystallite collisions and their time-reversed twins
Conclusions
|