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 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.


Introduction
Doug Henderson and John Barker helped to set the stage for our own Nonequilibrium developments through their equilibrium work on Thermodynamic Perturbation Theory [1]. This novel approach solved the problem of calculating accurate liquid-state thermodynamics by approximating the structure of a liquid 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 Nonequilibrium Systems has just been published by World Scientific Publishers in Singapore and released in March 2015.

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 either 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 mathematician with roots in Saint Petersburg, Alexander Lyapunov . Around 1979-1980, Shimada and 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 satellite 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 orthonormalization 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 timeaveraged rate at which two nearby trajectories separate -can be determined from the growth rate of the first Lyapunov vector δ 1 (t ) e λ 1 t . 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 multiplier 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 twodimensional 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: 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): 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 procedure, 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) Lyapunov exponents [6].

13003-2
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: 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: 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 exponents 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 ideas mentioned 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 shocked material. 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 backward 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 imperfect. Though these calculations were bit-reversible, they spanned only tens of thousands of timesteps. We consider much longer simulations in the present work.  The leading vectors going forward in time emphasize the leading edges of the crystals, where the collision 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 making 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 describes 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 velocities, the colliding hexagons can bounce or coalesce. At higher velocities, several smaller crystallites or 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: In order to propagate the particles reversibly, we use Levesque and Verlet's bit-reversible algorithm [11] which we detail in the following section.

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: 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.

Combining Newtonian and Hamiltonian mechanics
In our molecular dynamics work we combine the Newtonian and the Hamiltonian forms of mechanics. The coordinate-based Newtonian Leapfrog algorithm advances pairs of coordinate configurations {q t , q t ±dt } while Hamilton's first-order motion equations advance from one {q t , p t } phase point to the next {q t +dt , p t +dt }: Combining the two forms of mechanics requires a definition of momentum based on the coordinate information generated by the Leapfrog algorithm. The unimaginative first-order choice, with errors can and should be improved upon by using instead the third-order definition [10], The formal error in this last definition is − (1/30) ..... q dt 4 . Figure 5. Pairs of 7-particle crystallites and their collision products for two choices of the initial velocities, ±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.

Numerical liquid models for the collision process
To minimize numerical errors, it is useful to choose very smooth force laws. As demonstration problems 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 w i 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: The normalization constant (5/h 2 π) is chosen so that the integral of the weight function is unity: 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:

Crystallite collisions and their time-reversed twins
Continuing to pursue simplicity we choose initial conditions with symmetric coordinates and momenta for the two colliding crystallites. Sample geometries are shown in figures 4 and 5.  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, 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 symmetry in a macroscopic way after approximately 50 000 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 symmetrized 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, magnifying them by a factor of 10 16 , 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: likewise taking care to compute symmetrized densities for each particle: The pairing of local (instantaneous) Lyapunov exponents has been much discussed [13][14][15][16][17][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 vectors 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 uncertainty, 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. Originally, at time zero, the center-to-center separation of the two crystallites was 50 with all particle velocities (0, +0.1) in the leftmost hexagon and (0, −0.1) in the rightmost hexagon. The least-energy nearestneighbor separation is 0.8611 2127 0463, 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.   edges is h = 3.5. The original center-to-center distance at time zero was 50. This system was followed for 10 000 000 timesteps alternating segments with 500 000 forward steps followed by 500 000 backward steps in order to ensure the convergence of the offset vectors both forward and backward in time. 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 235 000, 1235 000, 2 235 000, . . . 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.

13003-8
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 successors, 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,  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 supercomputer 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 anywhere in the world on a laptop computer. It is interesting that despite all the change in computers and computation, the basic principles of mechanics, and the conclusions emerging from them, are unchanged.

Conclusions
Many facets of classical mechanics still remain to be illuminated. Some of them directly concern liquids and dynamic instability. Today, computational advances make it possible to evaluate Lyapunov spectra 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 Thermodynamics. 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-  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][20][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]. The motivation for using periodic rather than rigid boundaries in molecular 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 possibility 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 identifying 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 Hamiltonian mechanics, it is clear that the present results do violate pairing in the forward time direction. In the reversed 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.