Reversible Simulation of Hard Sphere Gas

Juraj Kumičák

Remark: If your browser does not support Java virtual machine, necessary to view the animations, you can download it from this site.

Since the system of hard spheres (discs) is chaotic, many of the results given in the paper, and some properties observed in the animation of generalized baker map should apply to it. My primary goal was to develop the program simulating hard sphere gas which would be reversible to unlimited precision. The most serious obstacle on this way is the collisional dynamics of spheres, the simulation of which necessarily involves goniometric functions. The latter return irrational values and in such a way prevent attaining the "numerical reversibility". Therefore, I was forced to approximate "particles" by regular polygons and to adjust appropriately the mechanics of collisions. Moreover, to obtain reversibility, I had to perform also all the calculations with the use of rational numbers, i.e. with fractions of integers. Programs performing such precise calculations are, unfortunately, extremly slow and therefore I give here reversible animations serving only demonstrative purposes, representing a compromise between the precision and the speed. The stress was on reliable reversibility of the animation.

In all the presented simulations the particles move in the unit square and are being reflected by its boundaries. They are supposed to be regular polygons with 4n sides, four of which are parallel to the sides of the unit square. The particles are not allowed to rotate and they interact via elastic collisions only when in direct contact with each other. Their initial positions are equidistant and they can occupy initially the full unit square or its upperleft quarter. The initial velocities of all particles have the same absolute values and their directions are pseudorandom. Such initial distribution was chosen to recognize immediately if the initial state is reached after reversal. Four of the particles are colored to facilitate observation of their motion: the first three are painted with easy to remember RGB sequence, the last one (positioned approximately in the center) is orange. The black cross shows the position of the system's center of mass and the user is advised to watch its damped oscillatory motion which is the first sign of macroscopic "hydrodynamic" behavior.

To use the applets, enter first the number of particles in the text field (the number will be automatically rounded to the nearest smaller square) and click Set (do not use Enter instead). You can next optionally choose between small and big squares (by checking the Make particles bigger box) and between Large grid and Small grid. Choose Printout if you can open the Java console of your browser and wish to see the output there. You begin the animation by pressing the Run button. The particles start moving and the speed of animation can be changed any time using the slider. After stopping the motion you can reverse the velocities by clicking Reverse, and pressing afterwards the Run button, you will observe the motion you can never see in reality, since it is violating the second law of thermodynamics. The reversed state will at last return to the original positions (with reversed velocities) where it stops. During forward and backward run you can also pause the motion any time.

The graphical user interface shows the number of collisions and time elapsed. The latter is, however, not proportional to the physical time, rather to the computational one, since the speed of animation is roughly proportional to the number of collisions to be calculated. Typically, the dense initial states take essentially more time than the "diluted" ones. The small graph on the right shows the distribution of speeds (histogram) as the collisions take place. The blue line is the theoretical result for an ideal gas in thermal equilibrium – a two dimensional Maxwell-Boltzmann distribution.

Opening the Java console – if available – you get futher information about the system as collisions proceed. The printout of the data slows down the animation, so that you can switch it off by unchecking the Printout box. This can be done repeatedly during forward and backward run so that you can adapt the printout to your requirements.

It is recommeded that you start experimenting with the lowest number of particles, N = 4, at the moderate animation speed. You will see how the dynamics of collisions looks like and you will notice the time lag between stopping and pausing the run. This is due to the requirement of reversibility which led me to stop the run at well defined states. The mutual collisions of particles and their collisions with walls turned out to be such events.

Animation of squares with double precision numbers

The system of hard squares lacks one of the essential problems of numerical simulation: irrational numbers which are being introduced by goniometric functions in calculations of projections. The squares namely collide only in the x and y directions so that they only exchange horizontal or vertical components of their velocities, but do not change their values. The exclusive use of rational numbers in all calculations, however, slows down the calculations even under these favorable conditions. Fortunately, this does not prevent to have reversibility after thousands of collisions. Therefore, I present at first the animation of hard squares, using double precision numbers, just to have a reversible simulation with acceptable speed. The user is advised to play with this model at first.

The above properties of the hard squares dynamics imply also that the Maxwell-Boltzmann distribution cannot be reached, as can be seen by observing the histogram. The reversibility of this model is then guaranteed for moderate number of particles and for fairly large number of collisions. To have an idea about what one can expect, I note that with 400 small particles starting in small grid I have observed perfect reversibility after 500 000 collisions, corresponding to 303 time units as shown in Time elapsed. If there are, however, too many particles, multiple collisions occur which may cause irreversibility, but sometimes even the multiple collisions are restored during the backward run. When using the small grid, the reasonable number of big and small particles should not exceed 11 x 11 = 121, and 23 x 23 = 529, respectively. For the large grid the numbers may be higher, but they are limited in any case by 30 x 30 = 900 particles.

The Java console shows futher information about the system as collisions proceed: number of (remaining) collisions, coordinates and velocities of the "central" (orange) particle and coordinates of the center of mass.

Animation of octagons with rational velocities

The reversible animation demonstrable by an applet is necessarily only approximate but is expected to retain – at least in certain respects – the most salient properties of reversible systems with macroscopic irreversible behavior. Some of such properties can be observed in the next approximation to the hard sphere 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 of the squares. To calculate the projections of velocities onto the two new directions requires to multiply corresponding values by the factor of cos(π/4) = √2/2 = 1/√2. Simulations, using rational (or double precision) numbers for both positions and velocities, confirmed that the contributions of positions to irreversibility are negligible so that one can obtain reasonable approximation of reversible dynamics using rational values for velocities only. To this effect one can approximate √2 by a rational number and use for cos(π/4) in the forward evolution √2/2 and in the backward one 1/√2. This has the effect that √2 is crossed away in reversed evolution so that the velocities are perfectly recovered. Practical tests show that 1.4142 represents sufficient rational approximation of √2 to obtain reversibility after reasonable number (at least thousands) of collisions.

Large octagons are replaced in the animation by circles, small ones by squares, within which the octagons are inscribed. Java console output gives positions and rational velocity components of the colliding particles. The absolute age τ (defined, however, differently) is also given, derived from numerators in rational expression of velocities of all particles. One observes that the Maxwell-Boltzmann distribution is readily approached by this system.

The rational numbers are handled in calculations with the help of the class BigInteger from the Java.math package. Note the growth of integers representing the rational velocities with growing number of iterations and their fall after reversal. The growing numbers are slowing down the calculations limiting the number of usable iterations.

It is evident that the same behavior will be seen with all regular 4n-gons. The only difference expected will be the speed of approach to the equilibrium distribution of positions and velocities.


The inspiration to develop the applets comes from the nice piece of work of Julio R. Gea-Banacloche. I have made use of the author's graphical user interface and of his overall approach to molecular dynamics simulation. The idea of perfectly reversible simulation using polygons is, however, mine.

(August 22, 2005)

Back to contents

Author's official and private home pages.

Hosted by www.Geocities.ws

1