Nonergodic dynamics of force-free granular gases

We study analytically and by event-driven molecular dynamics simulations the nonergodic and aging properties of force-free cooling granular gases with both constant and velocity-dependent (viscoelastic) restitution coefficient $\varepsilon$ for particle pair collisions. We compare the granular gas dynamics with an effective single particle stochastic model based on an underdamped Langevin equation with time dependent diffusivity. We find that both models share the same behavior of the ensemble mean squared displacement (MSD) and the velocity correlations in the small dissipation limit. However, we reveal that the time averaged MSD of granular gas particles significantly differs from this effective model due to ballistic correlations for systems with constant $\varepsilon$. For velocity-dependent $\varepsilon$ these corrections become weaker at longer times. Qualitatively the reported non-ergodic behavior is generic for granular gases with any realistic dependence of $\varepsilon$ on the impact velocity.

Ergodicity is a fundamental concept of statistical mechanics. Starting with Boltzmann, the ergodic hypothesis states that long time averages of physical observables are identical to their ensemble averages [8]. In the wake of modern microscopic techniques such as single particle tracking [9], in which individual trajectories of single molecules or submicron tracers are routinely measured, knowledge of the ergodic properties of the system is again pressing: while time averages are measured in single particle assays or massive simulations, generally ensemble averages are more accessible theoretically. How measured time averages can be interpreted in terms of ensemble approaches and diffusion models is thus imminent [10,11].
Here we demonstrate in quantitative detail how ergodicity is violated even in simple mechanical systems such as force-free granular gases. We analytically derive the time and ensemble averaged mean squared displacements (MSDs) and show that for both constant and viscoelastic restitution coefficients the time average of the MSD is fundamentally different from the corresponding ensemble MSD. Moreover, the amplitude of the time averaged MSD is shown to be a decaying function of the length of FIG. 1: Collisions in a free granular gas with a below-unity restitution coefficient lead to a cooling of the gas. Along with the reduced kinetic energy of the gas particles, the diffusion coefficient in a free cooling granular gas decreases with time.
the measured trajectory (aging). Comparison to the effective single particle underdamped scaled Brownian motion (SBM) demonstrates that this behavior is due to the non-stationarity invoked by the time dependence of the granular temperature. Ballistic correlations of the granular gas are shown to be relevant beyond this effective SBM description. Our results for generic granular gases are relevant both from a fundamental statistical mechanical point of view and for the practical analysis of time series of granular gas particles from observations (in particular, of granular gases in Space) and simulations.
Granular gas particles collide inelastically and a fraction of their kinetic energy is transformed into heat stored in internal degrees of freedom. The dissipative nature of granular gases gives rise to many interesting physical properties [2]. In absence of external forces the gas gradually cools down. During the first stage of its evolution, the granular gas is in the homogeneous cooling state characterized by uniform density and absence of macroscopic fluxes [2], realized, e.g., in microgravity environments [12]. Eventually instabilities occur and vertexes develop [2]. Here we focus on spatially uniform systems.
Collisions are quantified by the restitution coefficient ε = |v ′ 12 /v 12 |, the ratio of the relative speed of two granu-lar particles after and before a collision event [13]; ε = 1 denotes perfectly elastic collisions, while ε = 0 reflects the perfectly inelastic case [14]. For 0 < ε < 1 the granular temperature T (t) = m v 2 /2 given by the mean kinetic energy of particles with mass m continuously decreases according to Haff's law for granular gases [15], Here is the inverse characteristic time of the granular temperature decay, involving the initial value of the inverse mean collision time τ −1 c (t) ∝ T (t)/m. Weak dissipation (ε ≃ 1) thus implies τ 0 ≫ τ c . Due to the temperature decrease the self-diffusion coefficient D(t) of the gas is time dependent [2,[16][17][18], (Fig. 1). For ε = 1 we recover normal diffusion with constant diffusivity. Most studies of granular gases assume that ε is constant. Different approaches consider the relative collision speed dependence of ε as ε(v 12 ) ≃ 1 − C 1 v 1/5 12 + C 2 v 2/5 12 [19]. The coefficients C 1 and C 2 depend on material properties and the sizes of gas particles [20]. The granular temperature of the viscoelastic gas scales as T (t) ∼ t −5/3 [21,22] implying D(t) ∼ t −5/6 [17]. Especially in mixtures of viscoelastic granular particles a crossover from superdiffusion to a giant diffusion occurs [23,24].
Simulations. We perform event-driven Molecular Dynamics (MD) simulations [25] of a gas of hard-sphere granular particles of unit mass and radius, colliding with constant (ε = 0.8, Fig. 2) and viscoelastic (ε (v 12 ), Fig. 3) restitution coefficients. The particles move freely between pairwise collisions, while during the instantaneous collisions the particle velocities are updated according to certain rules: the post-collision velocities are given by the velocities before collision and the restitution coefficient ε. We simulate N = 1000 particles in a three-dimensional cubic box with edge length L = 40 and periodic boundary conditions. The volume density is φ ≈ 0.065.
We evaluate the gas dynamics in terms of the standard ensemble MSD R 2 (t) obtained from averaging over all gas particles at time t, as well as the time averaged MSD for a time series R(t) of length t as function of the lag time ∆. Here the angular brackets denote the average δ 2 (∆) = (1/N ) N i δ 2 i (∆) over all N particle traces. For an ergodic system, such as an ideal gas with unit restitution coefficient corresponding to normal particle diffusion, ensemble and time averaged MSDs are equivalent, R 2 (∆) = δ 2 (∆) [10,11]. In contrast, several systems characterized by anomalous diffusion with powerlaw MSD R 2 (t) ≃ t α (α = 1) or a corresponding logarithmic growth of the MSD, are nonergodic and display the disparity R 2 (∆) = δ 2 (∆) [10,11,[26][27][28][29][30]. While the ensemble MSD crosses over from ballistic motion R 2 (t) ∼ t 2 for t ≪ τ0 to the logarithmic law R 2 (t) ∼ log(t) for t ≫ τ0, the time averaged MSD starts ballistically and crosses over to the scaling δ 2 (∆) ∼ ∆/t given by Eq. (7). Fig. 2 shows the results of our simulations of a granular gas with constant ε = 0.8. The ensemble MSD shows initial ballistic particle motion, R 2 (t) ∼ t 2 . Eventually, the particles start to collide and gradually lose kinetic energy. The ensemble MSD of the gas in this regime follows the logarithmic law R 2 (t) ∼ log(t) (red line in Fig. 2, top) [2] (see also Supplementary Material (SM) [31]). The time averaged MSD at short lag times ∆ preserves the ballistic law δ 2 (∆) ∼ ∆ 2 . At longer lag times, we observe the linear growth δ 2 (∆) ∼ ∆ (black symbols in Fig. 2, top). In addition to this nonergodic behavior, the time averaged MSD decreases with increasing length t of the recorded trajectory, δ 2 (∆) ∼ 1/t. This highly non-stationary behavior is also referred to as aging, the dependence of the system on its time of evolution [30]. It implies that the system is becoming progressively slower. We observe the convergence lim ∆→t δ 2 (∆) → R 2 (t) . Fig. 3 depicts the MD simulations results for a granular gas with viscoelastic restitution coefficient ε(v 12 ). In this case the ensemble MSD scales as R 2 (t) ∼ t 1/6 for t ≫ τ 0 . The time averaged MSD does not seem to follow a universal scaling law but appears to transiently change from the power-law δ 2 ∼ ∆ 7/6 at intermediate lag times to δ 2 ∼ ∆ at longer ∆, see the bounds derived in SM [31]. As function of the length t of particle traces, we observe the crossover from δ 2 ∼ t −5/6 to δ 2 ∼ 1/t. Granular gas with constant ε. Let us explore this behavior in more detail. The dynamics of a granular gas can be mapped to that of a molecular gas by a rescaling of time from t to τ as dτ = T (t)/T (0)dt [32]. Using Haff's law (1), it follows that τ = τ 0 log (1 + t/τ 0 ). As function of τ the granular temperature remains con-stant, and the velocity-velocity correlation function for each spatial component decays exponentially [2,32], The MSD then has the exponential approach for the velocity correlator (β = τ 0 /τ v (0)). The MSD is At short times the particles move ballistically, R 2 (t) ∼ 3D 0 t 2 /τ v (0), crossing over to the logarithmic growth R 2 (t) ∼ 6D 0 τ 0 log (t/τ 0 ) as seen in Fig. 2 valid in the range τ 0 ≪ ∆ ≪ t, where τ 0 is the characteristic decay time of the temperature law (1). This result indeed explains the observed behavior of Fig. 2: the time averaged MSD scales linearly with the lag time and inverse-proportionally with the length t of the time traces. Comparison of Eqs. (6) and (7) demonstrates the nonergodicity and aging properties of the gas particles. Viscoelastic granular gas. For a velocity-dependent ε(v 12 ) the temperature decays like T (t) ≃ T 0 (t/τ 0 ) −5/3 , and the time transformation reads τ = 6τ 5/6 0 t 1/6 . The MSD in this case exhibits the long time scaling R 2 (t) ∼ 36D 0 τ 5/6 0 t 1/6 seen in Fig. 3, top. For the time averaged MSD we analytically obtain the bounds δ 2 (∆) ∼ ∆ 7/6 /t and δ 2 (∆) ∼ ∆/t 5/6 , compare the details in the SM [31]. These bounds are given by the dashed lines in Fig. 3, top. Concurrent to this change of slopes as function of the lag time, Fig. 3, bottom, shows the change of slope of δ 2 (∆) as function of the trajectory length t from the slope −5/6 to −1 at a fixed lag time.
Scaled Brownian motion. For unit restitution coefficient individual gas particles at long times perform Brownian motion at a fixed temperature defined by the initial velocity distribution. For the dissipative granular gases considered herein, the temperature scales like T (t) ≃ 1/t 2 and ≃ 1/t 5/3 , respectively. Single-particle stochastic processes with power-law time-varying temperature or, equivalently, time dependent diffusivity D(t), are well known. Such SBM is described in terms of the overdamped Langevin equation with diffusivity D(t) ∼ t α−1 for 0 < α < 2 [33,34]. SBM is highly non-stationary and known to be nonergodic and aging [11,34,35].
To study whether SBM provides an effective singleparticle description of the diffusion in dissipative granular gases we extend SBM to the underdamped case, driven by white Gaussian noise ξ(t) with correlation function ξ i (t 1 )ξ j (t 2 ) = δ i,j δ(t 1 − t 2 ) for the components. A granular gas with constant ε < 1 corresponds to the limiting case α = 0, giving rise to the logarithmic form R 2 (t) ≃ D 0 log(t) and the velocity correlation [36] v This result for the ultraslow SBM formally coincides with the velocity correlation function (5) for granular gases in the limit β ≫ 1, in which the velocity correlation time τ v (0) is much smaller than the characteristic decay time τ 0 of the granular temperature. This is achieved for sufficiently small dissipation in the system (ε 1). A more careful analysis shows that there exists a difference of the SBM model to the full result for the granular gas due to ballistic correlations [31]: The restitution coefficient affects only the normal component of the velocity of the colliding particles, while the tangential components remain unchanged. Therefore, in a granular gas with constant restitution coefficient the trajectories of particles become more and more aligned, effecting longtermed correlations. In the case of viscoelastic particles the restitution coefficient tends to unity when the relative velocities of colliding particles decrease with time. Therefore, at long times, the trajectories become more chaotic, and the ballistic correlations disappear. In this viscoelastic case ballistic correlations play a significant role only for intermediate lag times, while they become suppressed at longer lag times, effecting a gradual crossover between the bounding behaviors δ 2 (∆) ≃ ∆ 7/6 /t and ≃ ∆/t 5/6 ( Fig. 3, top). In contrast, for a constant ε the ballistic correlations are relevant during the entire evolution of the gas, canceling the leading term of the SBM result for the time averaged MSD [31].
Conclusions. The occurrence of nonergodicity in the form of the disparity between (long) time and ensemble averages of physical observables and aging, is not surprising in strongly disordered systems described by the prominent class of continuous time random walk models involving divergent time scales of the dynamics [10,11,26,27]. Examples include diffusive motion in amorphous semiconductors, structured disordered environments, or living biological cells [11].
In contrast to such complex systems, we here demonstrated how nonergodicity arises in force-free granular gases. This may seem surprising for such simple mechanistic systems. However, the physical reason for the nonergodicity is due to the strong non-stationarity brought about by the continuous decay of the gas temperature.
For a constant restitution coefficient the MSD R 2 (t) grows logarithmically, while the time averaged MSD δ 2 (∆) is linear in the lag time and decays inverse proportionally with the length of the analyzed time traces (aging). We derived the observed nonergodicity and the aging behavior of granular gases from the velocity autocorrelation functions. We note that aging in the homogeneous cooling state of granular gases was reported previously [37], however, it was not put in context with the diffusive dynamics of gas particles.
The decaying temperature of the dissipative force-free granular gas corresponds to an increase of the time span between successive collisions of gas particles, a feature directly built into the SBM model [34]. As we showed here, SBM and its ultraslow extension with the logarithmic growth of the MSD indeed captures certain features of the observed motion and may serve as an effective single-particle model for the granular gas, which is particularly useful when more complex situations are considered, such as the presence of external force fields.
Granular gases represent a fundamental physical system in statistical mechanics, extending the ideal gas model to include dissipation on particle collisions. Granular gases are a reference model in granular matter physics [38] with applications ranging from interstellar clouds or planetary rings to technologies in food and construction industries [2,38]. Our results shed new light on the physics of granular gases with respect to their violation of ergodicity in the Boltzmann sense. Moreover, we obtained the detailed influence of ballistic correlations in the granular gas dynamics. Both results are important for a better understanding of dissipation in free gases as well as the analysis of experimental observations and MD studies of granular gases.
We thank F. Spahn  Under the assumptions (S22) the contribution of the terms (S23) can be neglected. Finally, assuming that the second term in Eq. (S21) is small enough compared to Eqs. (S24) for the range of parameters k 1,2 chosen, we get to leading order δ 2 (∆) ≃ 36k Then we get in the limit ∆ ≪ t that In addition to these analytical estimates, we computed numerically the full expression for the time averaged MSD δ 2 (∆) = δ 2 0 (∆) + Ξ(∆). It agrees well with our MD simulation data, compare the curves in Fig. S.4 where we explicitly plot δ 2 /∆. It shows that in the range of parameters τ 0 , τ v and D 0 consistent with the MD simulations presented in the main text, the transient scaling behavior δ 2 (∆) ∼ ∆ 7/6 is realized in a limited range of ∆. Note also that in this range the linear SBM scaling for δ 2 (∆) as prescribed by Eq. (S27) is no longer valid. The reason is that for large values of ∆, when ∆ → t for any length of the trajectory, in Eq. (S19) the evolution of the time averaged MSD with the lag time ∆ becomes inherently nonlinear.