On the Arrow of Time

We believe the following three ingredients are enough to explain the mystery of the arrow of time: (1). equations of dynamics of gas molecules, (2). chaotic instabilities of the equations of dynamics, (3). unavoidable perturbations to the gas molecules. The level of physical rigor or mathematical rigor that can be reached for such a theory is unclear.


The Theory
There have been a lot of writings on the arrow of time; nevertheless, the problem still remains a mystery. In the equilibrium thermodynamics context, the arrow of time refers to the second law of thermodynamics. Consider the example of gas in a box, each gas molecule moves according to the equation of dynamics (free flying and collision) which is reversible in the coordinate time; on the other hand, the whole body of gas can only move in an entropy increasing direction among thermodynamic equilibria, e.g. moving from half box to whole box. The entropy increasing direction can be viewed as an intrinsic time direction -the arrow of time. In fact, there are many similar irreversible features in nature; for instances, living things can only grow older, sugar and water do not spontaneously un-mix, broken glass does not spontaneously fix itself. Explaining all these arrows of time is a challenging task. According to L. Boltzmann [1]:"The second law (of thermodynamics) can never be proved mathematically by means of the equations of dynamics alone." Here we shall approach the problem with the following three ingredients: (1) equations of dynamics of gas molecules, (2) chaotic instabilities of the equations of dynamics, (3) unavoidable perturbations to the gas molecules. Trying to use chaos to explain irreversibility is not a new attempt, see for example I. Prigogine's writings [7] [9] [10] [8]. Of course, there are many other attempts from different perspectives. The most recent one is that of Maccone [4] from quantum perspective, see also its critiques [2] [3] [6] [5].
As Boltzmann said, the equations of dynamics alone will not be able to explain irreversibility. Unavoidable perturbations to the gas molecules have to be taken into account. These perturbations are inside the isolated system, and they can be for examples physical noises and chemical fluctuations. If the dynamics of the gas molecules is non-chaotic, then these small perturbations will not generate any substantial effect; and the dynamics of the whole gas body will be reversible at least partially. The reality is that the system of equations of gas molecules is a non-integrable system and its dynamics is chaotic. Under such chaotic instabilities, the small unavoidable perturbations totally alter and liberate the dynamics of the whole gas body during the relaxation time from one thermodynamic equilibrium to another. In fact, the final thermodynamic equilibrium is independent of the initial conditions of the gas molecules, under such chaotic instabilities and unavoidable perturbations. For the example of gas in half box moving to full box, even the coordinate time and velocity of every gas molecule are reversed at any time during the relaxation time from half box to full box, the final thermodynamic equilibrium will still be the full box. Thus we arrive at the diagram for the arrow of time, Figure 1. Starting from the coordinate time 'Now', in both forward and backward coordinate time, the thermodynamic equilibrium entropy is increasing. The entropy time is folded into only positive direction -the arrow of time. In a more generalized sense, for a living thing at 'Now', even its every molecule's coordinate time is reversed, the living thing will not get younger, rather continue to grow older! Next we will explain a little more on chaotic instabilities and the effect of unavoidable perturbations. The main character of chaotic instabilities is 'sensitive dependence on initial data'. This is commonly characterized by amplifications of initial perturbations. For the traditional lower dimensional chaos, such amplifications are exponential in time. For the dynamics of gas moleculesm such amplifications are actually linear in time as shown later in the numerical simulation of one dimensional gas dynamics. The dynamics of the gas molecules is chaotic in both forward and backward coordinate time. Due to the unavoidable perturbations, the real physical orbits of the gas molecules are not the mathematical orbits prescribed by the equations of dynamics and the initial conditions. Only between two consecutive perturbations among all gas molecules, the gas molecules follow a segment of the mathematical orbits. The duration between two consecutive perturbations is extremely small. So the real physical orbits quickly diverge from the mathematical ones. For an N particles system on the interval [0, 1] for example, the phase space has 2N dimensions, coordinated by the particles' locations and velocities. The region where all particles are in the half interval [0, 1/2], has measure 2 −N . The main theme of our study is to investigate the possibility of return to this small region when the time and velocities are reversed at some moment on the orbit originating from a point in this small region.

A Popular Challenger to Irreversibility -The Poincaré Recurrence Theorem
The unwarranted argument here is that due to the Poincaré recurrence theorem, full box gas will return near half box gas after enough time. Physicists making such an argument apply the Poincaré recurrence theorem without its condition. The crucial condition is that one can find a bounded invariant region in the phase space of all the gas molecules coordinated by their locations and momenta. Then inside the bounded invariant region, Poincaré recurrence theorem holds. The system of equations of gas molecules is a non-integrable system for which such a bounded invariant region does not exist. Even for near-integrable Hamiltonian systems, the well-known mechanism of Arnold diffusion can lead to unbounded drifting of the momenta.

Ergodicity
There are two types of chaotic dynamics of gas molecules when the whole system moves from one thermodynamic equilibrium to another. Type I chaotic dynamics is the 'transitional chaotic dynamics' during the relaxation time, while type II chaotic dynamics is the 'saturated chaotic dynamics' when the system is in a thermodynamic equilibrium. The transitional chaotic dynamics has a non-zero mean motion, and it is the non-zero mean motion that eliminates the possibility of ergodicity. The saturated chaotic dynamics has no mean motion and may have ergodicity. Since the system is in a thermodynamic equilibrium, its entropy maintains constant and is compatible with ergodicity of the gas molecules' dynamics.

Numerical Simulations
The main setup is as follows: Particles move in the interval [0, 1], particles collide with each other and the walls (the two boundaries 0 and 1 of the interval [0, 1]) elastically. In particular, after they collide, two particles exchange their speeds to maintain the conservation of momentum assuming all particles have the same mass; on the other hand, after it collides with a wall, a particle just reverses its velocity. 4.1.1. Reversibility. In Figure 2, initially all the particles are randomly positioned in the half interval [0, 0.5], and the initial velocities of the particles are randomly chosen from [−0.5, 0.5]. After each collision, a random perturbation of magnitude less than 10 −6 is introduced to the velocities, while no perturbation is introduced to the positions of the particles. Special care needs to be taken when a speed is small after a collision. As time increases, the particles will evolve into the full interval [0, 1]. After 10, 000 collisions, the particles velocities are reversed and the time is reversed too. Then after each new collision, check if the particles are all inside the half interval [0, 0.5]. The return probability P r is given by the quotient of the total number of collisions N r after which the particles are all inside the half interval [0, 0.5] and 10, 000: For any fixed number of particles N p , 10 such numerical experiments are conducted, and the average return probability P r is ploted in Figure 2.
In Figure 3, initially all the particles are randomly positioned in the full interval [0, 1], and the initial velocities of the particles are randomly chosen from [−0.5, 0.5]. No random perturbation is introduced. After each collision, check if the particles are all inside the half interval [0, 0.5]. The numerical experiment ends after 10, 000 collisions. The return probability P r is defined the same as above. The average return probability P r of 10 such numerical experiments is ploted in Figure 3.
In  collisions. The return probability P r is defined the same as above (with 10, 000 replaced by 50, 000). The average return probability P r of 10 such numerical experiments is ploted in Figure 4.
From our preliminary numerical experiments, we conclude that • Without perturbation, no matter how many particles are initially positioned in the half interval, and no matter how many total collisions, reversibility back to the half interval is always achieved as long as the computer accuracy is sufficient. This is due to the reversibility of the mathematical orbits. • With perturbation, reversibility back to the half interval is possible for small number of particles (N p < 16). When the number of particle N p ≥ 16, extremely tiny perturbations (of magnitude less than 10 −6 ) will lead to irreversibility. This is due to the sensitive dependence on initial data of the orbits in the phase space. Take  in the 32-dimensional phase space is chosen, from which a mathematical orbit exists. This rigorous mathematical orbit is reversible. But such a mathematical orbit has sensitive dependence on initial data. Even a tiny perturbation can lead to another diverging orbit. This causes the irreversibility of the numerical orbit when tiny random perturbations are introduced. We believe that the specifics of the perturbations and the specific ways in which they are introduced are important to irreversibility. What is crucial for irreversibility is the existence of perturbation.
In Figure 5, initially all the particles are randomly positioned in the half interval [0, 0.5], and the initial velocities of the particles are randomly chosen from [−0.5, 0.5]. After each collision, a random perturbation of magnitude less than 10 −6 is introduced to the velocities, while no perturbation is introduced to the positions of the particles. Special care needs to be taken when a speed is small after a collision. As time increases, the particles will evolve into the full interval [0, 1]. After different relaxation time (10000, 20000, 30000), the particles velocities are reversed and the time is reversed too. Then the particles are tracked the same amount of Due to the size of the computation, we usually run the number of experiments N around 50 to a few hundreds. 4.1.2. Chaotic Dynamics. First we look at the time series plots for two particles case (N p = 2). In Figure 6, for four different initial conditions, the location time series of the right particle are plotted. These time series look chaotic. But such chaotic dynamics may not be in the chaos in the classical sense. Bear in mind that at a collision, the velocity has a discontinuity. Next we explore further the chaotic dynamics along the line of Liapunov exponents.
If we introduce small perturbations only on the initial locations of the particles, this will cause the perturbed collisions and the unperturbed collisions happen at a little different time, then the velocity deviation can be order one during a small time interval inbetween the perturbed and the unperturbed collisions. On the other Figure 6. The location time series of the right particle in the two particles system. hand, such a velocity deviation does not contribute substantially to the location deviation. So the location deviation maintains the same order as the initial location perturbation. Our numerical experiments verified this. Therefore, the main contribution for irreversibility comes from velocity perturbations.
If we introduce small perturbations only on the initial velocities (or both on the initial velocities and on the initial locations) of the particles, then the initial velocity deviation will generates a location deviation which is linear in time, with speed given by the initial velocity deviation. Near collisions, the location deviation will diverge from this main line, resulting in a rough oscillating continous nondifferentiable location deviation function in time. The velocity deviation mainly follows the horizontal line in time given by the initial velocity deviation. Near collisions, it will just to other horizontal lines of order one away from the horizontal axis. Overall, the velocity deviation is a piecewise continuous function in time, with discontinuities at the the collisions, see  close to the horizontal line, given by the initial velocity deviation of order 10 −6 : where (v 1 , v 2 ) and (v 1 , v 2 ) are the unperturbed and perturbed velocities of the two particles. The other four lines are near:  In the case of Figure 11(a), right after the collision, both particle 1 and both particle 2's location deviations will decrease. In the case of Figure 11(b), right after the collision, particle 1 (the circle)'s location deviations will increase, while particle 2 (the square)'s location deviations will decrease. More specifically, right after the collision, particle 1's velocity deviation is v 2 − (v 1 +v 1 ), while particle 2's velocity deviation is v 1 − (v 2 +v 2 ).
Therefore, particle 1's velocity has an extra deviation (due to the collision) of  Thus the extra gain (due to the collision) of particle 1's location deviation is the same with the extra loss of particle 2's location deviation. In term of 1 -norm (sum of every particle's location deviation), in the case of Figure 11(b), the total location deviation of particles 1 and 2 is unchanged in lieu of the collision. On the other hand, in term of 2 -norm (Euclidean norm), the total location deviation may increase or decrease. So far, we never observed the 'increase' case. If a particle collides with a boundary, it is easy to see that the particle's location deviation will decrease right after the earlier collision. Scattered points below the main line in the location deviation always appear when the experiment time is long enough, see Figure 12 for the two particles case. If we increase the size of the initial perturbation, scattered points appear earlier. In  Figures 13 and 14, we use the initial perturbation size of 10 −4 and the experiment time of 5000 and 10000, we can see a global picture in time of the location deviation for the two particles case. Bear in mind that the blue region represents a continuous non-differentiable function in time.
In conclusion, what we have observed is that the particle dynamics has sensitive dependence upon initial data, but instead of "exponential" sensitive dependence like for traditional chaos, we observed here "linear" sensitive dependence; that is, the particle location deviation induced by an initial perturbation grows linearly   Limited numerical experiments indicates that the return probability P r (defined above) from the half to the whole rectangle follows the law of P r = 2 −N where N is the number of disks; while in the 1D case, the return probability does not exactly follow this law due to particles blocking each other. 4.2.2. Chaotic Dynamics. In the 2D case, the dynamics of the disks is completely chaotic.
First (Figure 15) we choose the number of disks N = 6, and the radius of each disk r = 0.05. The initial velocities of the disks are randomly chosen with speed bounded by 10. The initial locations of the disks are also randomly chosen within the rectangle. By adding a perturbation of size 10 −15 to the initial locations, we observe the location deviations of the disks in time. Figure 15 shows the time series of the location deviations in Euclidean norm.
Next we compute the Liapunov exponent. For a perturbation of size 10 −15 to the initial locations, we stop the computation when the Euclidean norm of the where ∆p represents the Euclidean norm of the location deviations, and ∆t represents the time step. Then we average these 100 quantities to get an estimated Liapunov exponent for this numerical experiment. Finally, we average these estimated Liapunov exponents for 200 such numerical experiments to obtain the final estimated Liapunov exponent. The dependence of the final estimated Liapunov exponent upon the size and the number of the disks is shown in Figure 16. Finally the slide shows on the dynamics of six disks are shown in Figures 17  and 18.

Epilogue
We are convinced that the three ingredients mentioned in the section on Theory are enough to explain the mystery of the arrow of time. The level of physical rigor or mathematical rigor that can be reached for this theory is unclear. The numerical simulations above (especially in the 2D case) offer positive support to our Theory.      Figure 18. The same setup as in Figure 17 except that the perturbation size is 10 −12