Small and Simple Systems That Favor the Arrow of Time

The 2nd law of thermodynamics yields an irreversible increase in entropy until thermal equilibrium is achieved. This irreversible increase is often assumed to require large and complex systems to emerge from the reversible microscopic laws of physics. We test this assumption using simulations and theory of a 1D ring of N Ising spins coupled to an explicit heat bath of N Einstein oscillators. The simplicity of this system allows the exact entropy to be calculated for the spins and the heat bath for any N, with dynamics that is readily altered from reversible to irreversible. We find thermal-equilibrium behavior in the thermodynamic limit, and in systems as small as N=2, but both results require microscopic dynamics that is intrinsically irreversible.

The second law of thermodynamics is the fundamental barrier that prevents us from going back in time.
Indeed, this "2 nd law" causes the unidirectional increase in entropy needed for the "arrow of time" [1][2][3][4][5][6] and the maximum in entropy for thermal equilibrium.Because other basic laws of physics are reversible, the 2 nd law dominates our daily lives.However, most (but not all [7][8][9]) scientists believe that irreversible behavior emerges from reversible microscopic dynamics only for systems that are large and complex.As an early example, Boltzmann often sought to obtain the irreversible 2 nd law from reversible Newtonian kinetics, with mixed success.His goal was to justify the maximum in entropy needed for thermal equilibrium, now embodied by Boltzmann's factor   ∝  −  / , where Boltzmann's constant () connects internal energies (  ) to an ideal heat bath at temperature .(An ideal heat bath is effectively infinite, with weak but essentially instantaneous thermal coupling to the system [10].)Although experiments must ultimately decide the range of 2 nd -law behavior [11][12][13][14], the models we describe here have the advantage of dynamics that is easily changed from reversible to irreversible, system sizes that are precisely fixed, and entropy of the entire system (including heat bath) that can be calculated exactly at every step.
Computer simulations provide unique insight into how macroscopic behavior emerges from microscopic dynamics [15][16][17].In fact, the first molecular dynamics (MD) simulations on a digital computer surprised Enrico Fermi and his team by showing that small anharmonic systems governed by Newton's laws fail to reach thermal equilibrium [18].Although MD simulations of larger systems usually exhibit thermalequilibrium averages [19,20], when Fermi's ideas are extended to study dynamics, many systems have energy fluctuations that diverge from standard statistical mechanics as  → 0 [21].This disconnect between Boltzmann's factor and Newton's laws arises when conservation of local energy overwhelms the weak coupling to the heat bath.Here we study simpler systems as a function of size and coupling to the bath.We find that an intrinsically irreversible step in the microscopic dynamics is needed for the 2 nd law.
All systems studied here are based on the Ising model, consisting of interacting binary degrees of freedom ("spins") on a periodic lattice [22].We start with the one-dimensional (1D) Ising model having  spins in a 1D "ring" (to avoid anomalous endpoints).The potential energy is  = − ∑      +1  =1 . Here,  is an interaction energy, and   governs whether this interaction ("bond") occurs between neighboring spins at lattice sites  and  + 1, with   between spins at  =  and  = 1.Thus, bonded spins have interaction energy − if they are aligned or + if they are anti-aligned, or 0 if they are unbonded.In the standard Ising model, all spins are bonded to all nearest neighbors,   ≡ 1.Here we usually assume that broken bonds (  = 0) occur intermittently, allowing a thermal-equilibrium distribution of bonds  = ∑    =1 < .These broken bonds are needed to maximize entropy for the 2 nd law, plus they add complexity via a third energy value between neighboring spins.The 1D Ising model can be treated analytically for systems of any size using Boltzmann's entropy [23] Eq. ( 1) Here   is the number of high-energy bonds between anti-aligned spins and  0 is the number of noninteracting bonds, leaving  −  0 −   low-energy bonds between aligned spins.Note that the term 2  0 is a consequence of  0 separate segments of spins (having non-bonded endpoints) in a continuous ring, which is replaced by 2  0 +1 for free boundaries, or when  0 = 0.However, the special case of  0 = 0 is negligible if  ≫ 10 and  is not too low.Then, the system rarely fluctuates far from the condition that maximizes the entropy of mixing,  0 / = 1/2.At the energies investigated here, equilibrium values are: 0 / = 0.51 − 0.57.In terms of the bond types given in Eq. ( 1), the potential energy becomes: Although simulations of the Ising model usually utilize algorithms based on Boltzmann's factor, we argue that such simulations yield both a conceptual and physical mismatch by assuming that an effectivelyinfinite heat bath instantly couples to a small set of spins.A more balanced approach, introduced by Creutz [24], couples the spins to a similarly-sized set of "demons" (Einstein oscillators, all having the same spacing between discrete energy levels).In Creutz's model, each demon serves as a local heat bath that exchanges energy with nearby spins; thus, the demons serve as a source of kinetic energy, .This model enhances the realism of simple simulations by requiring that total energy ( =  + ) is always explicitly conserved, and by ensuring that the spin system and its heat bath have similar sizes and impact on each other.Of special interest is when the Creutz model is made reversible [25], where simulations can be run for an unlimited sequence of steps, then if reversed, after the same number of steps the system returns to its exact initial state.Such precise reversibility is rare in simulations of complex systems, where round-off error and sensitivity to initial conditions often yield divergent trajectories [26,27].Although many types of cellular automata show reversibility [28], the Creutz model also has discrete states for its heat bath, yielding thermal statistics for small systems and thermodynamic behavior for large systems.
Unlike most simulations of Ising spins that utilize a fixed temperature from an ideal heat bath, simulations of the Creutz model conserve energy, so that  is found from thermal averages.Consider  sources of , each with energy-level spacing  to match the spins.If each  was coupled to an ideal heat bath, its average energy would be given by Bose-Einstein statistics:  ∞ = /( / − 1).Inverting this equation yields the average temperature: / = 1/ ln(1 + / ∞ ).Although finite-size effects in small systems corrode the concept of , the exact Boltzmann's entropy for  sources of  can be written for systems of any size [29]: In fact, we find accurate thermal averages for  ≥ 2 from counting the ways that  =  +  can be distributed between  and .Specifically, the entropies, Eq. ( 1) and Eq. ( 3), yield the multiplicities,   =    / and   =    / , and the net probabilities: , = (  *   )/ ∑ (  *   ) , .
Eq. ( 4) Table 1 gives key properties of the  = 2 system having total energy / ≡ 1.Even in a relatively short simulation, all possible spin states and bonds are likely to occur, with weighted averages given in the bottom row of Table 1.Similar tables can be constructed for systems having other sizes and energies, but such tables become unwieldy for large systems.
We simulated Creutz-like models as detailed in the Supplemental Material.The main part of Fig. 1 shows the resulting total entropy density as a function of time for various simulations of a system having  = 210 6 and / ≡ 1.Four lines in Fig. 1 A (black, green, blue, and red) show that the entropy during each initial simulation rises sharply from an initial state, with rates that increase with increasing randomness.Specifically, for the entropy to reach 95% of the maximum value, irreversible dynamics (red) requires a single sweep, while reversible dynamics requires seven sweeps with global  (blue) and eleven sweeps with local  (green); whereas systems with no broken bonds (black) never approach this maximum.Two lines in Figs.
1 A and C (red and magenta) show that systems with intrinsic randomness are irreversible.All other lines are reversible, precisely following every step back in time, even after more than 2.6210 11 steps.
Figure 1 B displays these same simulations over an intermediate interval of times on a greatly expanded scale.Although initial simulations (green, blue, and red) have considerable overlap in their fluctuations, their time-averaged entropies (solid horizontal lines) increase significantly with increasing randomness (standard errors are less than the line thickness).Furthermore, subsequent simulations exhibit a large jump in entropy from dynamics that is reversible (grey and cyan) to irreversible (magenta), which is even clearer when time-averaged (dashed horizontal lines).Thus, Fig. 1 establishes that the total entropy from irreversible dynamics is significantly and persistently higher than that from reversible dynamics.
Moreover, if the spins are initially in an equilibrium state due to irreversible dynamics, this entropy evolves to a lower value if the dynamics becomes reversible.
The inset in Fig. 1   .The total entropy in Fig. 2 A tracks the entropy of  in Fig. 2 B, but mirrors the entropy of  in Fig. 2 C. Thus, reversible dynamics fails to reach maximum entropy because of the explicit heat bath, which is absent in simulations utilizing an ideal heat bath where total entropy is not calculated.Our results are consistent with MD simulations showing deviations from standard statistical mechanics [21], attributable to the explicit conservation of energy and intrinsic reversibility of Newton's laws.Now focus on the behavior during the middle third of Fig. 2.During these middle times, all simulations have the average rate of bond-change attempts reduced to 1/10 th the rate for spin-change attempts.Note that the total entropy is significantly altered only for reversible dynamics.Interestingly, this total entropy is reduced during the initial simulation, but increased during subsequent simulations.Thus, reversible dynamics yields non-equilibrium steady states with entropy that depends on the relative time scales of the dynamics, and on the initial conditions [30,31].Whereas, entropy is significantly and consistently higher for irreversible dynamics.128, green symbols), but such finite-size effects are expected when a fixed total energy is insufficient to thermally occupy higher energy levels.Thus, Boltzmann's factor applies only in the thermodynamic limit, and only if the dynamics is intrinsically irreversible.

Symbols in
Figure 3 displays various entropy densities, and their differences, as a function of .The main plot has logarithmic axes with 2 ≤  ≤ 210 6 , while the inset has linear axes with 2 ≤  ≤ 16.Line and symbol colors identify the total energy during each simulation (see legends).The uppermost lines (dash-dotted) show that   / increases with increasing , approaching a constant value at large .Symbols in the inset show that   / from simulations (circles) coincide with  ℎ / from theory (squares), found from summing the entropies of all states weighted by their multiplicities, similar to those given in Table 1.The inset also shows (  −  ℎ )/ ℎ (triangles) in percent.Although Boltzmann's factor cannot describe such small systems, for intrinsically irreversible dynamics, thermal equilibrium from the entropy-weighted sum over all states remains valid down to  = 2.Of course, thermal equilibrium fails for  = 1 because there can be no randomness in choosing the single .
Symbols in the main part of Fig. 3 show the size dependence of the difference between entropies, (  −   )/ = ∆/.Symbol color and type identify the total energy and heat bath, given in the legend.
Lines, from fits to the data using ∆/ = / +  with s and c constants, show general agreement with the behavior from local (dashed) and global (solid) .Irreversible dynamics increases the entropy, ∆ > 0, except for three blue triangles missing from Fig. 3 where  = (1  4)10 4 with / ≡ 2. Having ∆ < 0 can be understood from Fig. 2 D by the number of values of effective  needed to describe the dynamics, which changes from many values for small systems (green symbols), to two values for large reversible systems (black symbols) but only one value for irreversible dynamics (red symbols).More importantly, even as  → 2 where ∆/ > 0.01, the inset shows that only   matches calculated multiplicities.Most importantly, when  > 10 4 , ∆/ tends to increase with increasing , opposite to the behavior needed for   to approach   in the thermodynamic limit.
We identify oscillations in local energy as the primary cause of non-thermal behavior during reversible dynamics.Grey and cyan lines in Fig. 1 manifest these oscillations in the entropy at the start (A) and end (C) of the simulations, and in the inset from the peak at highest frequencies.These oscillations occur as fluctuations in local energy are first traded back-and-forth between  and  before dissipating via bonded spins.Similar fast oscillations in energy during MD simulations cause analogous deviations from standard statistical mechanics [21].Thus, Boltzmann's factor fails to describe fluctuations in large and complex systems governed by Newton's laws, and in our simple systems governed by reversible dynamics.In both cases, the failure of standard statistical mechanics can be attributed to conservation of local energy overwhelming the weak coupling to an explicit heat bath.However, when our simple systems have intrinsically irreversible dynamics, red symbols and lines in Figs. 1 and 2 A show significantly higher entropy, and a well-defined  as  → ∞ (Fig. 2 D).Furthermore, the inset of Fig. 3 shows thermal-equilibrium behavior for systems as small as  = 2, but only if the dynamics is intrinsically irreversible.
We speculate that intrinsic irreversibility in real systems may come from wavefunction collapse when coupling to a bath.In our case, each step involves one spin coupling to one , so that  = 2 may be related to a double-slit experiment.For reversible dynamics, the choice of  follows a prescribed sequence, similar to always knowing which slit passes the particle.Whereas irreversible dynamics involves intrinsic randomness in the choice of , similar to the intrinsic randomness of wavefunction collapse when a particle is measured as it passes through a slit.Thus, our results suggest that the measurment processes may extend to microscopic systems, such as when a single spin senses a heat bath of at least two particles.Related ideas have previously been proposed for the 2 nd law [1,3,8,9].
Conclusions -Our simulations reveal some unanticipated results in the thermal behavior of simple systems.An intrinsically irreversible step is needed for maximum entropy in systems of all sizes.
Specifically, intrinsically irreversible dynamics is needed for equilibrium thermal statistics from entropyweighted sums in systems as small as  = 2, and for Boltzmann's factor to apply in the thermodynamic limit,  ≫ 10  1. Properties of  = 2 system having total energy / ≡ 1.The first column shows representative configurations for low energy (•), high energy (×), or no energy (⊙) bonds; with spins that may be up (⇑) or down (⇓).The next four columns give the number of no-energy ( 0 ) and high-energy (  ) bonds, / (Eq.( 2)) and   / (Eq.( 1)).The next two columns give  =  −  and   / (Eq.( 3)).The final two columns give   =   +   and  , (Eq.( 4)).The bottom row gives weighted averages.  and / ≡ 1. Line color identifies the type of dynamics and thermal bath for initial simulations (main axis), and when averaged over subsequent simulations on an expanded scale of 10 2 in A and C, and 10 5 in B, with a common offset in each case.Specifically, red (magenta) lines show initial (averaged) behavior of irreversible simulations.Initial (averaged) reversible dynamics using local  is shown by the green (grey) lines, whereas dynamics using global  is shown by the blue (cyan) lines.Horizontal lines in B show timeaveraged values of each type.Black lines in A and C are also from reversible dynamics using global , but with no broken bonds.The inset shows power-spectral densities as a function of relative frequency on a double-logarithmic plot, with line color matching that in the main figures.We simulate Creutz-like models as follows.Each simulation utilizes a 1D ring of  spins and  sources of , with one source at each spin site.The spins are initialized into a low-energy configuration, with the first half of the spins in the ring oriented up and the second half oriented down.This configuration allows simulations at relatively low energies without the risk of freezing all spins into their fully-aligned state. is initialized by adding units of energy to randomly-chosen sites until  reaches the set value for the simulation.A single step in the simulation proceeds by first choosing a site, then attempting to change the spin and bond associated with that site.Changes that may occur include inverting the spin and/or making or breaking the bond with a neighboring spin.Such changes occur if-and-only-if energy can be conserved.
For example, if making the change increases , there must be sufficient energy at the site to keep  ≥ 0 after the change, consistent with demons acting as sources of kinetic energy.Alternatively, if the attempted change does not increase  the step always proceeds, with any reduction in  given to .Such steps are repeated  times to yield a single sweep.The sweeps are repeated many times, often yielding ~10 11 steps per simulation.For post-processing evaluation, a set of several quantities is recorded in a file.
Each recorded set includes a moving average of /, /, /, and   of the lowest eight levels of .
For convenience, file size is limited to a total of 2 17 = 131,072 sets of data per simulation.For short simulations (or large systems), the data set is recorded after every sweep, whereas for long simulations of smaller systems, each quantity is averaged over 10, 100, 1000, or 10,000 sweeps before recording.
We intermix reversible and irreversible simulations to study both types of dynamics.Specifically, the first simulation is reversible, so that the system returns to its far-from-equilibrium starting state at the end of the simulation.The second simulation continues from the final state of the first simulation, redistributes , then proceeds irreversibly.The third simulation continues from the final state of the second simulation, redistributes , then proceeds using the reversible dynamics of the first simulation.Such simulations are repeated to yield 6-10 simulations.Initial evolutions of reversible and irreversible dynamics come from the 1 st and 2 nd simulations, respectively.Equilibrium thermal behavior comes from averaging subsequent simulations, e.g. 3, 5, 7, and 9 th simulations for reversible dynamics and 4, 6, 8, and 10 th simulations for irreversible dynamics.
For irreversible dynamics we utilize the Mersenne Twister pseudo-random number generator [32] that is seeded randomly for effectively random behavior.First, a spin site is chosen randomly from the lattice.
Then, the  that will govern the dynamics is chosen randomly.Often, this single source of  governs attempts to change the spin and its bond, in random order.Other simulations utilize two randomlyordered sources of , one for attempts to change the spin and the other for its bond.Because all such simulations show similar results, we identify the crucial ingredient to be the continually varying choice of which source of  is used for each attempt.
All algorithms used for simulating reversible dynamics involve creating arrays of randomly-ordered (but fixed) sequences of sites, with each site occurring once (and only once) in each array.The sequences are chosen randomly to reduce correlations, but fixed so that the dynamics can be reversed by reversing the sequence.Here we focus on two algorithms.The "local " algorithm utilizes a single array for the site to be updated each step, with the single  at that site governing the attempt to first change the spin and then its bond.The "global " algorithm utilizes a second randomly-ordered array for the spin-change  and a third randomly-ordered array for the bond-change , with alternating spin-first or bond-first attempts.
The entropy increases slightly by tripling the number of randomly-ordered arrays, but the increase is small indicating that no reversible algorithm can match the maximum entropy of irreversible dynamics.
Figures 1 A, B, and Cshow, respectively, the dynamics at the beginning, near the middle (vertical scale 105 ), and end of the simulations.The upper three sets of data in Figs.1 A and Calso have an expanded vertical scale (10 2 ), with a common offset in both cases.Two sets of simulations are reversible, with the dynamics of each step governed by either the local  (green and grey), or global  (blue and cyan) chosen from the entire bath using randomly-ordered (but fixed) arrays.The third set of simulations (red and magenta) are irreversible, with global  chosen using a new random number each step.Main colors (green, blue, and red) come from the initial simulation, with spins always starting in the same highly-aligned state.Secondary colors (grey, cyan, and magenta) come from averaging three subsequent simulations of each type, alternating reversible then irreversible dynamics, with spins starting in the final state of the previous simulation.All reversible simulations (green, grey, blue, and cyan) start and end in the same state, while irreversible simulations (red and magenta) show no tendency to return to their initial state.Horizontal lines in Fig.1B show that when averaged over all intermediate steps, net entropy increases slightly with increasing randomness for the initial simulations (solid), while subsequent simulations (dashed) show a sharp increase from dynamics that is reversible (grey and cyan), to irreversible (magenta).Black lines in Figs.1 Aand C come from the initial simulation of a similar system, with global  and reversible dynamics, but no broken bonds.The peak entropy densities of all other simulations are about  ln(2) higher, as expected from the entropy of mixing when broken bonds are allowed.Such large increases in entropy indicate that real systems must also intermittently break their interactions, if physically feasible.
shows power spectral density (PSD) as a function of relative frequency () using the same sets of simulations and line colors given in the main part of the figure.Each PSD is found by taking the absolute value squared of the Fourier transform of √.The two sets of irreversible simulations (magenta) often overlap, giving a measure of the uncertainty.These PSD decrease monotonically with increasing f for 10 log() > 30, consistent with the overdamped relaxation shown in Fig.1A.Lines from reversible simulations having local  (grey) and global  (cyan) overlap at high frequencies.Both of these PSD increase monotonically with increasing f for 10 log() > 30, reaching a peak at maximum f, consistent with the fast oscillations shown in Figs.1 A and C. With decreasing f, reversible dynamics utilizing local  (grey) shows a 1/f-like divergence, indicative of slow energy diffusion at long times.

Figure 2 A
Figure 2 A shows moving averages of the total entropy from the simulations of Fig. 1, with Figs. 2 B and C from the separate entropies of  and , respectively.Each symbol comes from averaging 10 4 sweeps, positioned at their median time.Global  is used for both reversible dynamics (black squares) and irreversible dynamics (red circles).Open symbols come from initial simulations, with closed symbols from averaging three pairs of subsequent simulations.Error bars (visible when larger than the symbol size) give the standard error from averaging the subsequent simulations.Total entropies (Fig. 2 A) of the initial reversible simulation (open squares) are nearly as high as the initial irreversible simulation (open

Fig. 2 (
D) show time dependences of the ratio of probabilities from adjacent levels in .Specifically: ln(  / +1 ) with  = 0 (squares),  = 1 (circles),  = 2 (up triangles), and  = 3 (down triangles).Thus, if Boltzmann's factor can be used for , each symbol gives ∆/, where ∆ =  +1 −   = .Coinciding red symbols indicate that a single  applies to all levels only for irreversible dynamics of large systems.Black symbols show that reversible dynamics requires one value of effective  for  = 0 and  = 2, another value for  = 1 and  = 3, with additional values for middle times when bond-change rates are reduced.The concept of a single  also fails for irreversible dynamics in smaller systems ( =

FIG 1 .
FIG 1. Main figures show / as a function of time at the beginning (A), ending (C), and intermediate (B) times of simulations having  = 2106 and / ≡ 1. Line color identifies the type of dynamics and thermal bath for initial simulations (main axis), and when averaged over subsequent simulations on an expanded scale of 10 2 in A and C, and 10 5 in B, with a common offset in each case.Specifically, red (magenta) lines show initial (averaged) behavior of irreversible simulations.Initial (averaged) reversible dynamics using local  is shown by the green (grey) lines, whereas dynamics using global  is shown by the blue (cyan) lines.Horizontal lines in B show timeaveraged values of each type.Black lines in A and C are also from reversible dynamics using global , but with no broken bonds.The inset shows power-spectral densities as a function of relative frequency on a double-logarithmic plot, with line color matching that in the main figures.

FIG 2 .FIG 3 .
FIG 2. Moving averages of / from the simulations shown in Fig. 1 for spins (C), bath (B), and combined (A).Symbol style identifies dynamics as reversible (black squares) or irreversible (red circles), from initial simulations (open) or from averaging over subsequent simulations (filled), with error bars from the standard error.During the middle third of each simulation the average rate of bond-change attempts is reduced to 1/10 th the rate of spin-change attempts.(D) Moving averages of the ratio of probabilities between adjacent energy levels in .Symbol shape identifies the levels:  = 0 (squares),  = 1 (circles),  = 2 (up triangles), and  = 3 (down triangles).Black and red symbols are from the simulations shown in A-C having reversible and irreversible dynamics, respectively.Green symbols are from irreversible dynamics in a system having  = 125.
4. The failure of standard statistical mechanics to describe reversible dynamics can be attributed to conservation of local energy causing excess fluctuations that are not fully constrained by weak coupling to a realistic heat bath.Because maximum entropy behavior are often observed in our daily lives, our results suggest that intrinsic randomness may be common in nature.I gratefully acknowledge insightful comments from Sumiyoshi Abe, Ruth E. Kastner, Vladimiro Mujica, and George H. Wolf.I also thank Research Computing at Arizona State University for use of their