Time-Symmetry Breaking in Hamiltonian Mechanics. III.
A Memoir for Douglas James Henderson [1934-2020]

Doug (Figure 1) had the good fortune to do his “Theory of Liquids” PhD work with Henry Eyring at the University of Utah in the company of Francis Ree. Francis completed his own PhD work on “Random Walks” with Eyring a year earlier. Meanwhile Bill was completing his “Virial Series” PhD with Andy De Rocco at the University of Michigan. All three budding chemical physicists, Doug, Francis, and Bill, got together repeatedly over the next forty years. The details of Doug’s research career are spelled out in a 59-page Curriculum Vitae at Brigham Young University’s “Emeritus Faculty” portion of the Chemistry and Biochemistry website.

After nearly a decade of teaching jobs at the University of Utah, Arizona State University, and the University of Waterloo, Doug settled into a 20-year research career at IBM’s Almaden Research Center in San José California. It took only a little time for him to contact like-minded physicists at the nearby Lawrence Radiation Laboratory. At LRL about a dozen such, including Francis and Bill, had been brought together by the efforts of Berni Alder. All of us had the benefit of working with Doug’s colleagues at IBM, John Barker and Fareed Abraham, both of them in the forefront of theoretical and computational explorations of fluid properties.

This IBM-LRL cross-fertilization led to joint workb3 , and to a series of West Coast Statistical Mechanics Meetings, held at various University of California locations and the Almaden Center. These activities continued through the personnel cutbacks of the early 1970s up to the 17 October 1989 Loma Prieta earthquake, which killed dozens and injured thousands of Bay Area residents. As a consequence of the earthquake Doug relocated to the Latter-Day Saints’ Brigham Young University at Provo Utah, 600 miles east, over 800 by road. A majority of his 500 publications were generated with hundreds of colleagues all over the world following his 1990 move. His collaborations and travels continued worldwide until his final illness this past Summer and Autumn of 2020.

Doug’s commitment to his Latter Day Saints’ faith was strongb4 . In our experience most scientists’ faith, if any at all, makes up only a small portion of their day-to-day activities. By way of explaining this Feynman gave a 1956 talk noting that science progresses through scepticism and doubtb5 . Scientists try to reduce the uncertainty inherent in models of the real world by comparing theories, experiments, and computation. Feynman thought it difficult to accept the simultaneous realities of religious faith and scientific discovery because any “metaphysical” [miraculous] aspects of faith, if definite enough to be checked, might later be found to contradict scientific inquiry. Feynman’s doubt that religious metaphysics can withstand scientific investigation brings the example of evolution to mind. Doug’s belief was clearb4 : “Life is [not] possible without God”. Many, probably most, nonspecialist scientists would instead accept Richard Dawkins’ Darwin-based explanations as reasonableb6 .

Doug’s seminal work with John Barker [1925-1995] on the understanding of liquid properties from a knowledge of hard-sphere structural and thermodynamic properties was the product of decades of research. We can do no better here than to quote from a 1998 review b7 of their joint contribution to the Reviews of Modern Physicsb8 , “What is Liquid?” :

it was a lifelong fascination: finding the interatomic and intermolecular forces, models, and simulation techniques which could reproduce first qualitatively, but ultimately quantitatively, liquid properties. Condensed rare gases were a recurring subject in [their] studies. [Their] work spanned the transition from crude mechanical liquid-structure models (Bernal’s aggregates of steel ball bearings and Hildebrand’s coloured gelatin spheres, suspended in a fish tank, come to mind) to the computer simulations of liquid structure which are commonplace today [1998]. In their early work on liquid perturbation theory, Barker and Henderson discovered that the van der Waals hard-sphere picture provided a workable structural model for computing static liquid properties. Liquid free energies, suitable for thermodynamic calculations, could be estimated by finding the density and temperature dependence of the optimum hard sphere diameter.

Onward from the 1970s hard-sphere perturbation theory has provided a useful route to the thermodynamic properties of simple liquids. By now this approach, for which Google returns over six million results, has been superceded by computer simulation. Today it is feasible to compute hundreds, or even thousands, of thermodynamic state points with straightforward Monte Carlo or molecular dynamics simulations. Pressure and temperature, rather than energy and volume, can be made to serve as the independent variables for atomistic simulations by using Nosé-Hoover dynamics, which was developed in the 1980sb9 ; b10 ; b11 .

Nonequilibrium versions of thermostated dynamics have provided accurate estimates of transport coefficients, the shear and bulk viscosities along with the thermal conductivity. These linear transport coefficients can also be obtained from equilibrium Green-Kubo correlation functions so that both the equilibrium and hydrodynamic constitutive relations follow from relatively simple algorithms. Wikipedia articles and Amazon’s book offerings provide up-to-date expositions of the details. Through the work of Doug and others the understanding of equilibrium thermodynamic properties is no longer a mystery. Today it is nonequilibium systems which provide interesting and challenging research areas. They include the scientific subject of this memoir, the tension between the irreversibility of real life and the idealized time-reversible models with which we describe irreversible processes.

“Time marches on”. The steady, unstoppable advancement of time is all too apparent in writing a Memoir. Rudolf Clausius recognized it in his famous 1865 version of the thermodynamic laws:

“Die Energie der Welt ist constant.” [ First Law ]

“Die Entropie der Welt strebt [strives] einem Maximum zu.” [Second Law],

The Second Law describes the unstoppable aspect of real life, stating that an isolated system’s entropy increases to its equilibrium maximum. In support of that law Ludwig Boltzmann’s H Theorem showed that a probabilistic description of dilute-gas collisions in an isolated system provides an (approximate) entropy model which cannot decrease b12 . Thus it appears that this unidirectional law rules out the time-reversible nature of atomistic mechanical models, both classical and quantum. The classical model leads to Liouville’s [exact] Theorem, which deals with an entropy defined by Gibbs’ phase-volume formulation. For an isolated system the Theorem shows that this entropy cannot change with time.

In principle, but definitely not in practice, given initial conditions, the entire past and present of a Hamiltonian system can be obtained by solving the motion equations. This oversimplified picture–the very existence of longtime solutions–is subject to a serious caveat: in principle it is necessary to specify coordinates and momenta with infinitely fine precision in order to begin the work of an “exact” solution. No matter how precise these conditions exponential Lyapunov instability guarantees an ultimate frustration. In view of, or even despite, this nit picking, it is certainly reasonable to wonder why time-reversible dynamics is such a good model for the observed unidirectional evolution of the world around us. In our view this interesting subject remains ripe for a better understanding based on manybody computer simulations of natural processes. The inelastic collision of two crystallites or drops to make a single larger drop is an exampleb1 ; b13 ; b14 . The organized motion required to reverse such a process looks entirely implausible. No doubt it is. On a more macroscopic scale the spreading of Rayleigh waves on the surface of a pond, launched by a thrown rock, is another such example. Reversing a strongly irreversible process requires an inordinate amount of information.

Conventional Navier-Stokes fluid mechanics supports the law. The instability of a violation, a reversed flow, is built into the model. When the pressure-tensor component $P_{xx}$ includes the viscous stress $\eta(\partial u/\partial x)$ (with $\eta$ the shear viscosity) the Navier-Stokes motion equation is intrinsically irreversible. In a reversed Navier-Stokes flow viscous stresses, proportional to the velocity gradient $(\partial u/\partial x)$ change sign, leading to exponential instability. This instability is the result of poor modeling rather than Lyapunov’s amplified roundoff error.

Simulations can be decisive in understanding instabilities. Figure 2 is an illustration, showing both the stable progress of a steady shockwave and the irreversibility of its reversed motion in a manageable small simulation of 8192 two-dimensional particles, a simulation requiring only a few minutes of computer time.

In our experience irreversible processes, when simulated numerically, invariably do evolve in the direction of increasing entropy, even though entropy resists a useful definition for nonequilibrium states. Simulations provide the best means for understanding the details enabling the Second Law. In Figure 2, with the advantage of considerable simulations with varying velocities, timesteps, initial conditions, and numbers of particles we adopted a short-ranged purely-repulsive pair potential for the dynamicsb15 :

Initially, in the uppermost frame, the entire 8192-particle cold crystallite of Figure 2 travels to the right at velocity $(0.97,0.00)$. This choice of initial velocity yields the uniform equilibrium hot fluid shown in the third configuration counting down from the top. The hot fluid is now twice as dense as the initial cold crystallite, with the conversion complete at a time $t=114.6$. Evidently the velocity of the left-moving shockwave is $(-221.7025/2)/114.6\simeq-0.97$, just as is required for a homogeneous twofold compression. The molecular dynamics simulation in the Figure was continued onward for an additional time of 57.3, to the configurational snapshot at the bottom of Figure 2. Rather than retracing its recent history (as might well be expected from the symmetry of the motion equations) the hot fluid generates something entirely new, the rarefaction fan shown at the bottom.

There are three different plausible futures from the maximum-compression state at time 114.6 shown in Figure 2. All of them, at time 171.9, are shown in Figure 3: [1] at the top the straightforward continuation of the molecular dynamics; [2] in the middle the evolution to that same time, 171.9, with a Maxwell-Boltzmann distribution of velocities at the “hot” kinetic temperature of 0.115; [3] at the bottom the result of changing the signs of all the velocities at the maximum-density configuration at $t=114.6$. The failure of this last operation to recapture the initial condition at the top of Figure 2 illustrates the Lyapunov instability of the reversed shockwave. This last future is not accessible to continuum mechanics due to that mechanics lack of time reversibility.

By now both atomistic and continuum simulations of shockwaves are familiar. Our own work goes back to 1967. For a few references see Bill’s 1979 workb16 comparing Klimenko and Dremin’s atomistic shockwave profiles to the predictions of Navier-Stokes hydrodynamics. Back in 1967 Berni Alder had termed shockwaves the “most-irreversible” of processesb17 , suggesting that they are well worth studying. We continue that irreversibility work here. In the atomistic snapshots displayed in Figure 2 it is noteworthy that the underlying equations of motion are Hamilton’s, equivalent to Newton’s $\{\ F=ma\ \}$. Both mechanics models are precisely time-reversible. Suppose we play a movie of the simulated shock process backwards. What would we see? An unstable impossible flow. The visual reversal is only made possible by previously storing the forward flow. In the projectile simulation of Figure 2 we found that the largest Lyapunov exponent is of order unity in either time direction. The exponent describes the exponential growth rate of small perturbations. Thus the exponential growth of roundoff errors, of order $e^{115}\simeq 10^{50}$ for the time interval illustrated in Figure 2, is a Herculean barrier to longtime reversibility.

The shockwave stagnation of this figure is taken from our recent Memoir to Berni Alderb1 . In this computational model, entropy has clearly changed in the melting and heating of the once-cold projectile. Despite the time-reversibility of the underlying dynamical simulation the second law of thermodynamics, entropy increase, is satisfied. As shocks are irreversible we were curious what would happen when an atomistic reversible simulation was reversed, either by changing the signs of all the velocities or by changing the sign of the timestep $dt$ in the integration algorithm. What we found was surprising. A reversed steady shockwave, with a width of about two particle diameters, was soon replaced by a nonsteady wave, growing in width linearly in the time, and soon spanning all densities between that of the hot shocked fluid and a zero-density gas. We soon learned that such a structure, well known in fluid mechanics, is a “rarefaction wave” and that its expansive behavior is easily understood for simple model fluidsb18 , as we will soon see.

At the start of the present work we guessed that a detailed picture of the reversed-shock instability could be obtained by defining and characterizing local one-particle Lyapunov exponents. Repetition of our earlier simulationsb1 , along with analyses of single-particle Lyapunov exponents was disappointing. Though conventional Lyapunov analysis places the instability squarely at the shock, single-particle analyses provided no useful information. We abandoned their study and took up the analyses of shocks and rarefactions which makes up the balance of this paper.

Figure 4 illustrates three methods for generating shockwaves with projectiles and boundary forces, as in our previous work. At the bottom we illustrate a fourth approach, developed for the present work. It avoids boundary potentials by modeling the collision of two mirror-image projectiles. Using the same short-ranged soft-disk potential model as beforeb15 ; b1 , with projectile velocities of $\pm 0.97$ we obtained twofold compression of the original stress-free triangular lattice. Figure 5 illustrates the collision of two 4096-particle projectiles, with a fourth-order Runge-Kutta timestep of 0.01, and an initial kinetic temperature (to break the perfect-lattice symmetry) of 0.0001. This simulation shows that the stress-free cold projectiles at density $\sqrt{(4/3)}$ and temperature 0.0001, are steadily converted to hot fluid, the same fluid state, with density $2\sqrt{(4/3)}$ and kinetic temperature 0.115, as in our earlier workb1 , but without the need of any boundary forces.

Because two shocks rather than one are involved in the conversion the projectile method requires only half the time of our previous stagnation-based work. Both projectiles are twofold compressed at a time of 58. Figure 6 shows the steady decline of the total kinetic energy $K$ during compression, followed by its recovery during the rarefaction [expansion] which begins at maximum compression. It is noteworthy that the time variation is symmetric about the minimum within an accuracy of about one percent.

Figure 5 shows five snapshots of the 8192-particle two-projectile system, equally spaced in time with an interval of 29. The middle image is accurately compressed twofold, but its subsequent expansion is not at all like that of its past history. The pair of rarefaction waves visible in the fourth and fifth images were a complete surprise. Let us focus next on what happens to the hot twofold-compressed fluid in the center of the figure when it expands. In order to analyze the results we introduce smooth-particle averaging, a useful technique for generating continuously differentiable field variables from discrete atomistic information.

Just as in the single-projectile simulations of Figures 2 and 3, we consider three different continuations of molecular dynamics simulations, in Figure 7 from the middle view of Figure 5: [1] Conventional molecular dynamics; [2] Molecular dynamics with the initial velocities chosen from an equilibrium Maxwell-Boltzmann distribution; [3] Molecular dynamics with all velocities reversed. To make quantitative comparisons of the densities in these three different scenarios beginning with the maximum-density configuration we will evaluate a local smoothly-continuous density function, $\rho(x)$ for each of them. We use a “smooth-particle” algorithm based on Lucy’s weight functionb1 with a smoothing length of $h=3$. The weight function $w(|dx|)$ is zero for distances beyond 3 and has two vanishing space derivatives there, guaranteeing a smoothly differentiable density profile throughout its range:

With 8192 particles and original projectile dimensions $[L_{x}=256\sqrt{(3/4)}=221.70]\times[L_{y}=32]$, densities in the range $|x_{g}|<(L_{x}/2)-h$ represent averaged contributions from an area of $32\times 6=192$. The density at any such gridpoint location $x_{g}$ is simply given by the sum of all such nearby particles’ smoothed contributions,

Figure 8 compares the three future density profiles: molecular dynamics, Maxwell-Boltzmann, and reversal. The detailed evolution appears in Figure 9 for the two-projectile method. The cold material density, $\sqrt{(4/3)}=1.1547$, is doubled to 2.3094 in the shocked fluid.

Continuing the simulation over the interval $58<t<116$ with either the old velocities or with new ones from a $T=0.115$ Maxwell-Boltzmann distribution makes no significant difference. We conclude from this comparison that the maximum-compression configuration is effectively at equilibrium and that reversing the velocities indicates a nonequilibrium correlation between coordinates and velocities that recaptures a part, but far from all, of the original time-zero configuration. Evidently the molecular dynamics, due to Lyapunov instability of order $e^{t}$, is unable to recapture the past. No doubt that a sufficiently precise simulation (well beyond research budgets) could reproduce the initial configuration within visual accuracy. Trials using timesteps $\{\ 0.04,0.01,0.0025,0.000625\ \}$ as well as the (orders of magnitude slower) quadruple precision version of the Runge-Kutta integrator confirmed that storing the old configurations or using a bit-reversible algorithm are the only feasible approaches to capturing the past. Recently we considered the simulation of a similar single-projectile rarefaction problemb1 . We found that reversing the velocities at maximum compression gave a largest Lyapunov exponent with its largest contributors localized near the shockwave front. That finding indicates that the shockwave itself is the main barrier to time reversibility.

The steady nature of shockwaves is well knownb18 and leads to a simple numerical method for finding their structure from a known model for the constitutive equations. Here we choose simple model mechanical and thermal equations of state suggested by linear dependences of the energy on density and temperature and satisfying the Rankine-Hugoniot relation equating the work done in shock compression to the energy change:

This relation follows from the “Symmetric” compression picture from Figure 4. Choices of pairs of mechanical and thermal variables satisfying the Rankine-Hugoniot relation for the model include

To complete the constitutive equations needed for continuum simulations we choose to use a constant unit shear viscosity, setting the scale of the model’s shock width at unity. The pressure tensor is $P_{xx}=P-(du/dx)\ ;\ P_{yy}=P+(du/dx)$, where $P$ is the equilibrium pressure. For simplicity we ignore heat conductivity. The steady nature of the wave implies that the mass, momentum, and energy fluxes are all constant. In the case we have chosen here they are

Knowing the values of the three constants the result can be written as an ordinary differential equation for the velocity gradient, $(du/dx)$. A Runge-Kutta solution for a similar problemb19 appears in Figure 10.

Another (more challenging) numerical approachb20 solves the partial differential evolution equations on a fixed “Eulerian” grid. The resulting partial differential equations are as follows:

Notice that the viscous term in $P_{xx}$ contributes to a diffusion equation, $(\partial u/\partial t)\propto(\partial^{2}u/\partial x^{2})$, leading to decay, $\propto e^{-x^{2}/4t}$, rather than divergent  growth, $\propto e^{+x^{2}/4t}$, which would follow from changing the direction of time . We were able to find solutions by following the description given on page 360 of our Kharagpur Lectures bookb20 . It is interesting that the time-dependent equations require a much smaller step than do the steady state equations. See Figure 12 for details. We verified that the two methods, steady-state and Eulerian, are consistent and turned to considering the detailed differences between shock and rarefaction waves.

The usual textbook explanation of such a wave’s behavior can be applied to the existence of shock and rarefaction waves. In a typical fluid the sound velocity increases with density. Imagine the initial density distribution shown in the caricature of Figure 11. A wave traveling rightward must spread while a wave traveling leftward must steepen, reaching a finite width due to the diffusive nature of viscosity. The rightmoving wave’s width grows without bound in the absence of physical boundaries. This simple explanation, illustrated in Figure 11, is responsible for the extreme instability of the shockwave to time reversal. A molecular dynamics simulation with $dt=0.000625$ keeps the energy constant to over a dozen digits throughout a fourth-order Runge-Kutta simulation, but still, in a relatively short time the reversed shockwave becomes a rarefaction.

Figure 12 illustrates the transformation, using 3000 spatial intervals and a timestep of 0.00005. Even so we found it necessary to remove an even/odd instability periodically (every 20,000 steps) by averaging the density, velocity, and energy according to the following schemeb20 :

The reversed shockwave broadens roughly linearly in the time, as can be seen in the Figure. With this successful simulation of the conversion of a shockwave to a rarefaction wave following velocity reversal we see that our molecular dynamics simulations are quite consistent with the predictions of simple computational fluid dynamics.

In this work we set out to avoid the flow perturbations introduced by boundary potentials and found that simulations using mirror-image pairs of projectiles provide high-quality simulations of shock and rarefaction waves. Smooth-particle averaging makes it possible to compare atomistic and continuum profiles. One can also formulate Lucy function recipes for the local thermodynamic properties. Local averages of kinetic temperature and the density can provide local entropies, from which entropy production could be analyzed. Although the future choices of Lyapunov-unstable dynamical simulations must depend upon the local density of accessible phase-space states, accessible through roundoff error and its amplification, we have not attempted such a project. No doubt a simpler few-body nonequilibrium system, such as a short conducting $\phi^{4}$ chainb21 or a single heat-conducting oscillator in a temperature gradientb22 would provide a better starting point.

The hydrodynamic shock/rarefaction instability we have discussed and modeled here is just one among the many instabilities found with continuum constitutive relations. The Rayleigh-Taylor instability, with a denser fluid pressing upon a less-dense one, is similar. The reason why the shockwave/rarefaction instability discussed here has no classic name is that laboratory experiments (as opposed to computer experiments) do not yet have the capability of reversing all the velocities in a fluid system!

We are grateful to the referee for his careful reading of the manuscript.