e-mail: [email protected]
Keywords: Reversible simulation, Hard disk fluid, Collisions of polygons, Absolute age, Computer algebra methods
PACS: 05.45.Pq Ln Er Mn
APPENDIX: Animation of Hard Disk Gas
The evolution of many particle systems cannot be solved analytically and that is why the only option left was numerical simulation, called molecular dynamics. Already in the beginnings of the latter a question arose about the possibility to preserve the reversibility of its solutions. For example, N. Margolus [1] gave a description of universal reversible cellular automata thus outlining the essentials of reversible computation. However, it is only after the appearance of the Levesque-Verlet bit reversible algorithm [2], which represents apparently the most efficient method used to obtain reversible trajectories, that the reversible computer simulation of dynamical systems has become an important field of computational physics [3].
Several authors have compared the behavior of exactly reversible simulations with standard molecular dynamics, demonstrating that even a negligibly small deliberate noise can cause the loss of reversibility [4]. Nevertheless, not too much attention is being paid to possible interrelations existing between numerical ir/reversibility and corresponding system's dynamics.
I have recently studied the model of Generalized Baker Map (GBM), treated as a system displaying microscopic reversibility combined with macroscopic irreversibility. The goal of the study was to analyze the interrelation between both these aspects of evolution and to arrive — with the help of the analysis — at an understanding of the origin of irreversibility in reversible systems. It seems plausible that many of the results concerning GBM described in the paper [5], and some properties observed in the animation of GBM1, should apply also to models mimicking more realistic systems. To corroborate this view, I have turned to the simulation of 2D hard sphere fluid which would be reversible to essentially unlimited precision. The most serious obstacle on the way to achieving this goal is the collisional dynamics of disks, the simulation of which necessarily involves goniometric functions. The latter return irrational values and in such a way prevent attaining the "numerical reversibility". To avoid this drawback, I have chosen a method different from that of Levesque-Verlet — actually a standard molecular dynamics simulation [6] which conserves energy — and I approximate "particles" by regular polygons adjusting appropriately the dynamics of collisions. All the calculations are then performed with rational numbers, i.e. with fractions of integers.
Very interesting aspect of the GBM model turned to be the "absolute age" t of its states, which can be assigned in an unambiguous way and with reasonable precision to any point [5]. It is derived from the length of integers expressing positions of points as rational numbers. It increases in equal steps from the past (being negative there) to the future (where it is positive) — except the small neighborhood of the "present" state where its behavior is less predictable. This age is absolute in the sense that it does depend only on the distance from the present, i. e., on the number of iteration steps. Since this notion is intimately connected to the origin of irreversibility, it was interesting to find out whether it will be possible to implement the latter into the hard polygon model.
In all simulations described here the particles move within the unit square and are being reflected by its boundaries.2 They are supposed to be regular polygons with 4n sides, n ³ 1, four of which are parallel to the sides of the unit square.
The particles interact via elastic collisions only when in direct contact with each other and they are not allowed to rotate (their motion is restricted to translations). It is supposed that in collisions they exchange components of velocities perpendicular to the sides in contact, irrespective of positions of their centers. Hence the dynamics can be treated as a sequence of two-body elastic collisions.
The system starts with all particles having the same speed with pseudorandom uniform distribution of directions over the interval [0; 2p), and placed on a regular square lattice covering the whole unit square or just a quarter of it. Such conditions were chosen to readily recognize the return of the system to initial state after velocity reversal.
The inner diameter of the polygons is chosen to be either 1/25 or 1/50 of the unit square side, so that the number of particles, corresponding to the maximum density (to avoid ternary collisions), is bounded from above by 400 and 900, respectively. The lower bound is four particles which corresponds to the minimum density.
Due to the hard disk potential, the calculations evidently reduce to consideration of all pairs of particles i and j and finding the minimum of times tij for their next collision. One then moves all particles forward in time until the collision occurs and calculates the postcollision velocities of the colliding pair.
Evidently, in the case of nonrotating squares one can compute in this way the dynamics with absolute precision if all the initial values are rational. The collisional time tij is, namely, a linear combination of rational values in such a case, and velocity components are just exchanged. With the growing n the situation is less simple but can be handled in a way to be explained later.
In order to compare the model with real fluid, we suppose that the chosen diameter of particles corresponds to about 1 nm in reality, so that to simulate the gas at normal pressure and temperature we should have about four particles in the box — the lowest number used in the simulations. Choosing the initial speed of particles so as to correspond to about 500 m/s we then obtain mean free paths comparable to those of real gases.
The simulation was programmed in the Java programming language that enables the program to be run as an application as well as an applet that can be embedded into an HTML page (good for demonstration purposes on the web). However, when used as an applet, the program was very slow for n > 1, and therefore it was necessary to make it faster without losing reversibility. Simulations, using rational (double precision) numbers for both positions and velocities, confirmed that the contributions of positions to irreversibility are negligible for small n so that one can obtain reasonable approximation to reversible dynamics using rational values for velocities and letting positions as floats. Such are the simulations given at my web site3 (denoted there as "animations"), representing a compromise between the number of particles and the speed of the system's evolution. These animations can serve consequently only demonstrative purposes, but the full program is free of such limitations.
I have started with n = 1, i. e. with hard squares. Such system is free from one of the essential problems of numerical simulation: the appearance of irrational numbers which would be otherwise introduced by goniometric functions in calculations of projections. The squares namely collide only in the x and y directions so that they exchange only horizontal or vertical components of their velocities, but do not change their values. Under such conditions, components of velocities are just exchanged between particles, so that speeds rapidly attain a distribution which will not change thereafter and which is determined by the initial one. This means that there is no universal distribution function, particularly, that the Maxwell-Boltzmann distribution, expressed for a 2D system of N particles with mean velocity v0 as
| (1) |
cannot be reached if the system has not started with it — see an example of the system starting with a "singular" initial distribution, Fig. 1. Moreover, the damped oscillatory motion of the center of mass is the first (and the only) sign of macroscopic "hydrodynamic" behavior caused by rigid boundaries.
FIG. 1: The histogram of initially identical velocities of 400 squares after 5000 collisions. The velocities are divided into 20 intervals and p denotes the percentage of particles with corresponding values. The figure demonstrates that the Maxwell-Boltzmann distribution (represented by continuous line) is not being approached.
The exclusive use of rational numbers in all calculations should guarantee exact reversibility. This was actually confirmed for 500,000 collisions with 400 particles and nothing suggested that more particles and/or longer runs should violate the reversibility.
One cannot assign absolute age t [5] to the states of hard squares system, since this variable is meaningful only if positions and/or velocities are fractions of integers which grow with growing time (number of iteration steps). One cannot expect such behavior with squares, since if we start with rational positions and velocities, the latter will belong to a finite set of combinations of starting values and positions will be linear functions of time. It is true that the positions are being represented by fractions of growing integers but this growth is not monotonous (it is actually chaotic) so that there is no way how to extract information about time elapsed from such data.
The system of hard squares is evidently ergodic and mixing with respect to positions, but these properties are not shared by velocities which are essentially preserved due to the above mentioned reasons.
We now go over to the next approximation, n = 2, to the 2D hard disk gas — the gas composed of regular octagons. The octagons are again not allowed to rotate so that they can collide in four different directions only, two of which are identical with those found in the case of squares and the two new ones are inclined with respect to vertical and horizontal directions by the angle of j = p/4, see Fig. 2. To calculate the projections of velocities onto the two new directions requires to multiply corresponding values by the irrational factor of cos(p/4) = Ö2/2 = 1/Ö2. To prevent the latter from causing irreversibility, we can think of cos(p/4) as Ö2/2 in the forward evolution and as 1/Ö2 in the backward one. This has the effect that factors Ö2, accumulated during the forward evolution, are being cancelled out in reversed evolution so that the velocities are perfectly recovered.
FIG. 2: Collisional dynamics of octagons. The octagons collide in the direction perpendicular to the line of contact so that components of velocities in this direction are exchanged and components parallel to it are being preserved.
In numerical calculations we must, of course, approximate Ö2 by a rational number. Such approximation evidently introduces an error into positions so that the longer the time required for the reversible evolution, the more places are to be kept in decimal approximation of the number. However, any finite time of reversible evolution can be achieved by appropriate choice of the approximation. Practical tests show that 1.4142 represents sufficient rational approximation of Ö2 to obtain reversibility for reasonable number (at least thousands) of collisions.
We therefore see that the described simulation of octagons cannot be absolutely reversible for an unlimited time. But since the deviation from reversibility can be made arbitrarily small, and the time of reversible evolution consequently arbitrarily long, the simulation retains the most salient property of reversible systems with macroscopic irreversible behavior — irreversibility is practically not observable during any finite run. We can call such a property practical reversibility.
It should be moreover noted that the behavior of positions is further influenced by the specific character of polygon scattering. Since the effective cross section is in this case discrete, the dependence of positions on initial conditions is "sensitive" only at a discrete set of angles (corresponding to the vertices of the polygon). One could call such sensitivity "singular" as it is defined on the set with zero measure. With growing n the set becomes larger so that the sensitivity is the more pronounced the greater the n. The simulation with octagons demonstrated that the irreversibility of positions is in this specific case still practically negligible.
The absolute age t [5] cannot be again brought into relation with positions due to their linear dependence on time (number of iterations). However, the velocities of particles are multiplied, after each collision occurring in the direction not parallel to x or y axes, by (rational approximation to) cos(p/4) so that a possibility arises that their rational expressions will consist of fractions of growing natural numbers. Actually, the simulation confirms that such situation is typical. Therefore, it should be possible to construct an explicit expression for the absolute age based on rational components of velocities of all N particles which will be growing both in the future and in the past — one could make t even approximately linear function of the iteration index k, t = const |k|. Denoting numerators and denominators in rational expressions of numbers by superscripts A and B, respectively, one of possible choices is represented by the formula
| (2) |
where vlx and vly denote x and y components of the velocity of the l-th particle and c is the scaling factor which adapts the age to be approximately equal to the iteration index k (we keep meanwhile s = 1). This may seem to be a bit artificial and unphysical but, nevertheless, it turns out to be very reliable.
Such age will grow both in the future and in the past. To make it positive in the future and negative in the past one should calculate s as follows. If the sum
| (3) |
of differences of absolute values of numerators and denominators of particles i and j is greater after the collision than before it, then s = +1, if it is smaller, then s = -1; otherwise it does not change. This rather speculative way to keep the sign of the age monotonous has, however, the drawback that the age is then not the function of the state of the system — which is our primary goal — but depends on two consecutive states.
FIG. 3: The histogram of velocities after 1200 collisions of 400 octagons demonstrates that the Maxwell-Boltzmann distribution (represented by the continuous line) is being approached quickly.
There is a distant resemblance of the age to the entropy in that both depend on the distance of the current state of the system from the initial (non equilibrium) state. However, the age is monotonously growing with time whereas entropy approaches a limit value. Entropy is the measure of equilibrium, the age is roughly proportional to collision number. This difference is best seen in the fact that age grows even if the system starts in equilibrium.
Since the change of velocities is now limited by less constraints, the equilibrium Maxwell-Boltzmann distribution could be attained in principle. Observations, actually, confirm that the latter is readily approached by this system, see Fig. 3. The above mentioned influence of n on sensitivity of positions explains why this approach depends on n. We see that already octagons represent, in most respects, qualitatively sufficient approximation to hard disk fluid. It is maybe worth mentioning, that the irreversible dynamics would have still higher rate of convergence.
It is evident that all the results described in this section will be seen with regular 4n-gons also for n > 2. The only difference expected will be the growing speed of approach to the equilibrium distribution of positions and velocities with growing n. This, in my view, illustrates the dependence of interparticle interactions on the ir/reversibility and approach to the universal Maxwell-Boltzmann distribution of velocities, the idea mentioned already in [7].
Numerical molecular dynamics simulation of fluids is of great theoretical and practical importance. It enables to replace real experiments with computational ones and allows to calculate relevant parameters of flows — such as, e. g., transport coefficients — departing from theoretical assumptions.
But the reversible microscopic basis of fluid behavior is not taken into account in the usual case of simulations which are based on double precision calculations. The inherent rounding errors, namely, generate the irreversibility of the models, which is just a numerical artifact with no physical counterpart. The errors, making such simulations microscopically irreversible, are hiding one of the important aspects of the problem: the observed macroscopic irreversibility should ensue from underlying reversible dynamics and under no conditions should it result from the lack of numerical precision or from the presence of numerical noise.
The fully reversible simulation, on the contrary, enables to observe irreversibility in a model with granted microscopic reversibility which otherwise would be impossible to attain in real experiments. That is why only such simulation can serve the basis for the exact explanation of origin of ubiquitous macroscopic irreversibility of microscopically reversible dynamics.
The reversible simulation of many-particle systems was undertaken by several authors who arrived at many interesting results, concerning mostly the approach to equilibrium [1,2] and influence of noise on reversibility [4,8]. My stress was more on the influence of system dynamics on the properties of reversible evolution. That is why I have chosen to simulate particles by polygons with 4n sides and tried to see how the growing n, entailing the approach of polygons to disks, will influence the behavior of the systems. The case of low n shows rapid appearance of thermodynamic behavior, the case of high n will be the subject of a separate work.
In the extreme case of n = 1, the system is less sensitive to changes of positions than those of velocities, since small changes of impact parameter (i. e. changes smaller than the particle diameter) do not influence the direction of particles after collision. This restricts the variability of speeds to the extent that the attainment of the Maxwell-Boltzmann distribution is being excluded. In the next step, n = 2, the discrete set of collisional angles is less restricted, and the octagons model described in the paper demonstrated that the reversible dynamics of octagons has already most macroscopic properties expected of gases with hard particle potential. As a result, no assumptions about irreversibility on microscopic level are needed in order to obtain macroscopic irreversibility.
The simulation using rational numbers for coordinates and velocities enabled to introduce the notion of absolute age t which proved to be rather fruitful for understanding the origin of irreversibility. The possibility to define t is due to the fact that velocities of particles contain "traces" of their evolution, which enables us, at least in principle, to recover the iteration index k from arbitrary state. Then the question, why we are not able to prepare states which would evolve contrary to the second law of thermodynamics, can be answered trivially by stating that we can prepare states with "simple" particle velocities (corresponding to "present" states) more easily then those with "complex" ones (corresponding to the future states which upon velocity reversal become the past ones). From this point of view, the irreversibility is related to substantially different requirements for preparation of the present and past states and this seems to be the full explanation.
It should be noted that the most salient property of the described simulation is its ability to regain the positions and velocities of the hard disks in the reversed trajectory. This makes the method comparable to that of the Levesque-Verlet bit reversible algorithm [2] where only positions are exact and energy is not being preserved.
The author is thankful to William G. Hoover for his many suggestions. This work was partly supported by the Scientific Grant Agency (VEGA) of the Slovak Academy of Science and Ministry of Education of the Slovak Republic under the Grant No. 1/2198/05.
1http://www.geocities.com/kumicak/BakAnima.html
2Periodic boundaries could be handled easier but they seem to be less realistic.
3http://www.geocities.com/kumicak/HSqGas.html