Clustering impact regime with shocks in freely evolving granular gas

. A freely cooling granular gas without any external force evolves from the initial homogeneous state to the inhomogeneous clustering state, at which the energy decay deviates from the Ha ﬀ ’s law. The asymptotic behavior of energy in the inelastic hard sphere model have been predicted by several theories, which are based on the mode coupling theory or extension of inelastic hard rods gas. In this study, we revisited the clustering regime of freely evolving granular gas via large-scale molecular dynamics simulation with up to 16 . 7 million inelastic hard disks. We found novel regime regarding on collisions between “clusters” spontaneously appearing after clustering regime, which can only be identiﬁed more than a few million particles system. The volumetric dilatation pattern of semicircular shape originated from density shock propagation are well characterized on the appearing of “cluster impact” during the aggregation process of clusters.


Introduction
Clustering in the granular gas is one of the fascinating phenomena in the granular physics [1][2][3]. In general, a freely cooling granular gas with an inelastic hard sphere (IHS) model is known to evolve from the initial homogeneous equilibrium state to the inhomogeneous clustering state, at which the energy decay deviates from so-called "Haff's law" [4]. Due to simplicity of this model, there have been much investigated by various methodologies with a wide coarse-graining levels from the microscopic (molecular dynamics (MD)) and the mesoscopic (Enskog-Boltzmann equation), to the macroscopic (Navier-Stokes (NS)-like continuum equations and Ginzburg-Landau theory). However, there are few studies based on theoretical approach after clustering regime.
Previously [5][6][7][8], we have investigated the macroscopic statistical properties on the freely evolving quasi-elastic IHS model by performing a large-scale event-driven molecular dynamics with mainly 512 2 [5], 4096 2 [7] in two-dimensional system (2D) and 128 3 [8] in three dimensional system (3D). These extensive numerical simulations have revealed a quite similarity with fluid Navier-Stokes turbulence, such as enstrophy cascade (Kraichnan-Leith-Bachelor theory), Kolmogorov scaling and Bose-Einstein condensation, which are well-known concept in the fluid turbulence [9]. Reynolds number with Taylor's microscale before clustering grow up and was estimated as R λ ∼ 327 (2D) and ∼ 1470 (3D) at the onset of clustering regime [7,8] respectively, which are sufficiently high and reasonable that the ideal granular gas system shows the specific laws of turbulent behaviors.
At the initial stage of the cooling process, density and velocity filed remain homogeneous. However, the kinetic e-mail: isobe@nitech.ac.jp energy (or granular temperature) decreases via inelastic collisions in time. In this homogeneous cooling state (HCS), the granular temperature T (t) decays according to the Haff's law, where τ is collision number and τ 0 is a fitting parameter, which can be estimated explicitly by the kinetic theory in dilute systems. The reduced temperature T * (t) for freely cooling IHS model is equivalently the reduced total energy in the hard sphere system since there is no potential energy. In the clustering regime, due to inhomogeneity in density, T * (t) deviate from Haff's law [4]. The asymptotic behavior of T * (t) after clustering regime has been predicted by several theories, which are based on the mode coupling theory [10] and the extension of inelastic (sticky) hard rods gas [11]. Those theoretical predictions show the asymptotic behavior as T * (τ) ∼ Aτ −d/2 (d ≥ 2) and ∼ Bt −d/2 (2 ≥ d ≥ 4), respectively, where A = A(r, ν), B and r are independent constant, d is dimension. In 2D case, both theories predict T * ∼ τ −1 and ∼ t −1 , respectively [12]. In this paper, we performed large-scale and long-time event-based molecular dynamics simulation [13,14] in IHS model with up to 16.7 million particles to clarify the validity of asymptotic behaviors predicted by these theories, and especially focus on inhomogeneous clustering regimes in the dense systems, at which the theoretical approach cannot be resolved.

Model and numerical method
We address the simple IHS model without any external force (such as gravity), i.e., so-called "freely cooling granular gas". At microscopic molecular dynamics level, the IHS model can completely be characterized by three independent parameters: the total number of disks N, the packing fraction ν, and the restitution coefficient r. In case of r = 1, the system is govern by only solid-fluid phase transition, so-called "Alder transition" (See, recent progress in Refs. [15,16]). The whole systems with a side length L in the two dimensional (2D) square systems include monodisperse hard disks up to N = 4096 2 with a radius of disk σ, in which we adopted the periodic boundary conditions. The relations between L and packing fraction ν are L/(2σ) = √ πN/ν/2. Initially, the system is prepared for the equilibrium states by a long preparation run in an elastic mono-disperse system with r = 1. For the dissipation, the normal restitution coefficient r < 1 of sequential binary instantaneous inelastic collisions is introduced (i.e., parallel to the relative position of the two colliding disks in contact), where the momentum is conserved and energy is dissipated. The collisons are evolved by an efficient eventbased algorithms (See, in Refs. [13,14]).
Because of monotonically decreasing total energy in the evolving process, the physical time scale is slowing down due to dissipation between inelastic collisions. The time rescaled system is also often considered, in which the particle velocities are rescaled so as to keep the energy of the whole system being constant. This procedure has a direct relation to the thermostat in the the methodology of molecular dynamics. In general, thermostat methods are categorized by deterministic or stochastic ones. This time rescaling procedure resembles "velocity scaling method". This procedure do not affect the particle trajectories in the hard sphere system, which just affect the time scale of relative velocities between future collides particles, i.e., next collision time. In the rescaled time system, the unit of evolving time in the original system t is modified as new where T (t) is the granular temperature as a function of t [17,18].

Asymptotic energy decay after Haff's law
It is known that the freely cooling granular gas systems become unstable and several spatiotemporal scales appear. The standard evolution scenario is that spatial instabilities occur from the initial homogeneous cooling state (HCS: 1st stage) to the velocity field (shearing regime: 2nd stage), and then, to the density field (clustering regime: 3rd stage). Finally, global structures in the system (final attractor regime: 4th stage) appear as a non-equilibrium steady state, or local inelastic collapse occurs, which depends on the restitution coefficient. Since we only focus on the quasi-elastic large system, i.e., r ∼ 1, inelastic collapse do not occur in our present simulation.
In Fig. 1, the typical decaying process of energy (or granular temperature) in the simulation with two packing fractions ν = 0.10 and 0.60 as a functions of t * s = t s (1 − r 2 )/(4t 0 ) are shown, where t 0 is the mean free time in equilibirum for each density. The solid red curve shows the result obtained by the simulations at (N, ν, r) = (1024 2 , 0.10, 0.90). At first, we confirmed that the energy decay is consistent with the Haff's law before clustering regime, in which we confirmed the relationship of t ∼ exp (t s ). After clustering regime, the decay of granular temperature deviates from the Haff's law, at which we found empirical relationship for time t, rescaled time t s and collision number τ are t ∼ t 2 s and τ ∼ t s , respectively. We estimate α ∼ 2 in the form of T * (t s ) ∼ t s −α for asymptotic power decay of temperature, which leads to T * (t) ∼ t −1 . The results are fairly consistent with the prediction of Ref. [10].
In the dense case (ν = 0.60), we performed simulations wtih three sets of parameter by changing system size as (N, r) = (1024 2 , 0.9941), (2048 2 , 0.9964), and (4096 2 , 0.9981), respectively. The largest system consists of up to about N = 16.7 million hard disks (2D). Before clustering regime, the relationship between t and t s are also consistent with the Haff's law as shown in the case of ν = 0.10. After clustering regime, we observe the exponents α ∼ 1, which are obviously smaller than that of the dilute case. Moreover, during the growing process of the clustering after t * s ∼ 10, the exponent seems to have a crossover to the different larger values, i.e., the energy decaying rate is increasing. These peculiar behaviors have not predicted by present theories and suggest unknown physical mechanisms regarding the dynamics of clusters, which are increasing inelastic collisions (i.e. decreasing energy) per unit time after clustering stages.

Collision rate
To investigate the reason for the deviation from the theories, we focus on the collision rate after the clustering regime. In Fig. 2, the time evolution of the collision rate  (2/N)(dτ/dt s ) in terms of t * s is shown, where dτ is the amount of binary collision number occurring during small time interval dt s . We found that the data obtained by the simulation with three different system sizes show almost universal scaling curves and the rapid decreasing collision rate after clustering begins since particles move collectively as a cluster, where velocity correlation near the particle takes positive values along to the flow. However, the collision rate increases again around t * s ∼ 10, which correspond to crossover time that temperature decay shows the different exponents of power decay as shown in Fig. 1.

Volumetric dilatation and shoke waves
To clarify the origin of the rapid increase of collision rate, we focus on the spatial distribution of the rate of volumetric dilatation Q(x), where u(x) is the velocity field at x (i.e., in 2D case, u(x) = (u x (x, y), u y (x, y))). This properties are related to the compressibility of the fluid. In case of the incompressible fluids, the volumetric dilatation rate should be zero, i.e., Q = 0. This have higher values if the strong change is occurred in spatial differentiation of velocity field such as shock wave. In Fig. 3, we visualize Q(x) after well developed clustering regime, where the corresponding density field is also shown as the reference (monotone snapshot drawn as upper right of Fig. 3). We can clearly identify the high Q regions (shown as red) in not only the edge of cluster but also inside of the cluster. Four series of snapshots for Q show the growing semicircular shape density shocks, shown inside of the solid open circles in Fig. 3, which are excited from the contact point of giant impact between clusters. (See, higher values of Q(x) drawn as red in Fig. 3.). These shocks propagates inside the clusters which cause the rapid increase of collision rate. These results indicated that the theoretical prediction of the asymptotic behavior of energy decay after the clustering regime seems not to be valid at high density and large-scale system.

Discussion
Although we confirmed that the theories seem to be consistent in relatively dilute and small system, we found the new regime regarding on collisions between "clusters" spontaneously after the clustering regime in dense case. In the new regime, the power decay of energy shows clearly slower than that of the dilute case. Furthermore, the power decay itself is also violated in the long time limit. Since the collision rate in the new regime increases rapidly, the novel two-step energy relaxation as a clustering scenario is expected as follows: At first, several string-like clusters are nucleated with an aggregation of inelastic hard spheres, which are distributed separately in space. Those nucleated clusters move actively, and deformed with complex behaviors, and gradually aggregating or coarsening. Next, a certain amount of cluster begins to collide each other in the time regimes of t * s ∼ 10. This giant impact between clusters preceded correctively and made round shape density wave propagation (shocks) on the cluster. This shock propagation results in the increase of collision rate, which causes a rapid growth of dissipation. At that time, the volumetric dilatation pattern of semicircular shape originated from shock density propagation in the clusters after cluster impact can be detected. This regime called as "clustering impact regime" can only be identified in more than a few million particles system. In our analysis of volumetric dilatation pattern of semicircular shape and collision rate clarify that the cause of giant impact between cluster is a completely novel physical mechanism in the field of freely cooling granular gas. This has not known yet, since "clustering impact regime" can be detected only in the systems more than several million particles systems. These findings might be relevant in the field of the late stage of planetary formation. With these results, the asymptotic predictions of energy (granular temperature) in the long time limit should be modified in dense and large systems. Theoretical breakthrough in correspondence with simulation results might be demanded. We will address the results for the other packing fraction and 3D case in the near future. Figure 3. Four time series of snapshots for the evolving process of the spatial distribution of the rate of volumetric dilatation Q(x) after clustering impact are shown, in which the propagation of semicircular shape density shock wave in the clusters appears. The corresponding density field is also shown as the reference at upper right.