Non-equilibrium response and slow equilibration in hard disk systems

The relaxation from a non-equilibrium state to the equilibrium depends on the methodologies and initial conditions. To investigate the microscopic mechanisms of equilibration systematically, we focus on the non-equilibrium response during the equilibration process induced by a disturbance of the homogeneous expansion of the simple hard disk systems. Large scale simulations by event-driven molecular dynamics revealed that an anomalous slow equilibration toward the liquid states emerges when starting from the co-existence phase. The origin of the slow decay mechanism is investigated using the probability distribution of local density and orientational order parameter.Their inhomogeneities seem to cause the anomalous slow equilibration.


Introduction
In the molecular simulations, the production runs for obtaining physical properties are generally performed after the "equilibration" (relaxation toward equilibrium states) from initial non-equilibrium states. Any arbitrary states in the many-body systems will relax toward the statistical equilibrium state if you perform standard molecular simulations, which are the molecular dynamics using Newton's equation of motion and sampling using Markov chain Monte Carlo [1]. The "equilibration" is just preliminary tasks on the research on our specific topics; however, it often requires a long computational cost. It sometimes exceeds our lifetime, especially in high-density systems and large scale systems. Recent issues such as crystallization, granular gas, and glass/jamming transition have been actively investigated in these decades, requiring large-scale simulation in dense molecular systems. The investigation of microscopic mechanisms for equilibration is crucial, including developing optimal algorithms to shortcut unnecessary tasks.
In three dimensional (3D) hard-sphere systems, the celebrated paper of computer simulation for both physical new insights and methodologies was published in 1957 on the solid-fluid phase transition [2] by Event-Driven Molecular Dynamics (EDMD) method [3,4], which is the called "Alder transition" [5]. In the melting problem of two-dimensional (2D) hard disk systems [6], equilibration at large-scale hard disk systems was essential to figure out the type of phase transition. Insufficient equilibration caused different conclusions and controversy for many decades. The invention of Event-Chain Monte Carlo (ECMC) [7] was the most crucial contribution to tackle the 2D melting problem [8,9]. Then its variants have been actively developed such as for all-atom model with generalized potential [10], hybrid scheme with EDMD [11] * e-mail: cmm12069@ict.nitech.ac.jp * * e-mail: isobe@nitech.ac.jp and Newtonian Event-Chain MC [12]. In general, particle positions' equilibration is much difficult, especially in highly dense systems rather than velocities due to the excluded volume effect being dominant. In EDMD, even though the particle and cluster composed of a dozen particles get trapped in the metastable states, the velocities can take the Maxwell-Boltzmann distribution after just a few collisions. Due to the equivalence of the equilibrium states generated by Monte Carlo and molecular dynamics, even when its paths to the equilibrium state are completely different, it is reasonable to use ECMC for positional relaxation at first. Because ECMC is much faster than EDMD for equilibration, which is the central idea of a hybrid scheme [11]. There are a few studies that have been reported previously on comparisons of equilibration efficiency between methodologies, on melting and crystallization [9,11,13]. The ECMC needs one parameter chain-length, which depends on the density and probably other factors. We have to optimize it empirically in the given systems. The guideline of surveying the optimal parameter is still elusive.
Our purpose of the present work is to systematically figure out the microscopic mechanisms of equilibration between situations and methodologies, in addition to the acceleration mechanism of algorithms. As a first step, we focus on the non-equilibrium response around the Alder transition in the large-scale 2D hard disk systems and investigate the non-equilibrium response to calculate each physical property against disturbance, i.e., homogeneous expansion as the simplest case, to the system. In Sec, 2, our model and numerical settings, including a definition of relaxation time and physical properties, are described. Sec. 3 show the results on the non-equilibrium response of physical properties, in which unexpected anomalous slow decay of equilibration is observed due to inhomogeneous spatial density around co-existing phase. Finally, concluding remarks are described in Sec. 4.

Model and simulation methods
We focus on the two dimensional (2D) mono-disperse systems consisting of N = 512 2 hard disks without any external force, friction and dissipation between disks. The disks with radius σ are placed in a L x × L y (= A) rectangular box (i.e. L y /L x = √ 3/2) with periodic boundary conditions. The basic units of the system are set by the mass of one particle m, the diameter 2σ, and energy 1/β, from which we derive the unit of time as 2σ √ βm, where β is 1/k B T with Boltzmann's constant k B and the temperature T . The hard sphere system in both 2D and 3D is governed by solid-fluid phase transition (so-called "Alder transition" [2,5,6]) changing the packing fraction ν = Nπσ 2 /A, which is the primal control parameter.
(In the recent large scale hard disk simulation [8,9], the phase changes from the liquid (ν < 0.700), co-existence (0.700 < ν < 0.716), hexatic (0.716 < ν < 0.720) and solid (crystal) (ν > 0.720), respectively.) As the initial states, the systems for each packing fraction above the transition ν = 0.700 ∼ 0.722 are carefully prepared after long-time equilibration by efficient ECMC [7] up to 10 13 collisions. After equilibration, the whole system is expanded homogeneously so as that the packing fraction is set at ν target = 0.698 (corresponding to the liquid phase). The production runs of equilibration are done using EDMD simulations [4,11]. We investigate the non-equilibrium response during equilibration focusing on the relaxation times of various physical properties O(t), such as pressure P * , local and global orientational order parameters Φ L 6 and Φ G 6 , and the number of nearest neighbor disks N i . Since O(t) changes from O(0) to O(∞) = O eq in which O eq is the equilibrium values at the target packing fraction (i.e., ν target = 0.698), we defined a function of relaxation 1 ≥ F ≥ 0) and estimate the relaxation time τ where F (τ) = 0.01. To investigate the origin of the physical mechanism in detail, the system is divided into small grids (N g ) to calculate the local packing fraction ν l (r g ) at the location of each grid r g . We studied the relationship between τ and variance of ν l (r g ) as 2 at t = 0 and the evolution of the probability density distribution P(ν l ) as a function of ν l . The following formula calculates dimensionless pressure P * , where Λ is the collision number per unit time.
The orientation order parameter is defined by equation (2) for each particle i, where i indicates imaginary units, and θ j i is the angle between the position vector from the nearest neighbor j to i and arbitrary axis (e.g. x-axis). The local and global orientation order parameters Φ L 6 and Φ G 6 are defined by eqs. (3) and (5), respectively, The number of nearest neighbors N i is estimated by a fixed-distance cut-off method where the cut-off length is set at the first minimum of the radial distribution function.
We also consider the number of disks N D where N i 6.
To avoid a large scattering of local properties in each grid due to a few disks located, we adopt a Gaussian filter to analyze and visualize the spatial distribution of density and orientational order parameters. We carefully optimize the filter parameters for smoothing not to affect the central claim of the results.  Figure 1 shows the relaxation times τ for various physical properties in terms of initial packing fraction ν. We observed all physical properties studied here, showing a long relaxation time around ν = 0.708 (corresponding to the co-existence phase). In general, it would be naturally expected that the larger difference between the states described by the control parameter (in the present case, the packing fraction ∆ν = ν − ν target ) results in the positive contribution of the relaxation time of physical properties; that is, ∆ν increases, then, τ also increases. To clarify the origin of these anomalous slow relaxation, the local density variance of the initial configurations σ 2 d is calculated under the assumption that the large spatial fluctuation of the local density in the co-existence phase might cause the non-trivial relaxation time. We found a positive correlation between the relaxation time τ and σ 2 d in Fig. 1 (shown as a blue dotted line), although the peak position might have a slight discrepancy with each other. All properties are averaged over six independent ensembles.

Typical equilibration process
In Fig. 2, the typical equilibration process in time at ν = 0.708 (co-existing phase) and 0.720 (solid phase) are shown, respectively. In the early stage of equilibration, both cases' tendency is similar that physical properties relax rapidly within a short time scale. Note that the peculiar slow decay of the global orientational order parameter Φ G 6 is observed compared with other local properties. We observe that the distinct difference at the intermediate stage of equilibration emerges between ν = 0.708 and 0.720. All properties relax immediately to the equilibrium systems' values with an approximately single decaying rate in the case of ν = 0.720. However, the multiple relaxations can be observed in the case of ν = 0.708, where there might be a couple of stages with a few different (slow) decaying rates caused by the unknown physical factor(s).
To explore the unknown physical factor(s), the evolution process of the spatial distributions of both the local packing fraction and the orientational order parameter is visualized in Fig. 3. Note that for illustrating spatial distribution for orientational order parameter in Fig. 3, cos α i for each disk i is considered where α i is the difference of argument between ϕ i 6 and ϕ G 6 as a reference, i.e., α i = arg (ϕ i 6 ) − arg (ϕ G 6 ). In both density and orientational order parameters, we confirm that the system is composing of two phases; the solid (blue area) and the liquid (red area) phases at ν = 0.708. In the evolving process, these two phases persist at least until t = 4, 038, and partially it continues for a more long time. On the contrary, at ν = 0.720, the entire system looks homogeneously ordered phase initially, which changes drastically until t = 2, 089 homogeneously disordered phase until t = 2, 089. The whole system has already been completed for equilibration toward the liquid phase. In other words, it is not necessary to relax regarding the spatial distributions of the local density at ν = 0.720 since the initial and equilibrated states are closely similar distributions with each other.

Time evolution of the probability density distribution of the local density
To confirm the origin of slow relaxation based on local density, we calculated the time evolution of the probability density distribution P(ν l ) of the local density. In Fig. 4, P(ν l ) at ν = 0.708 has a broad distribution with a plateau of ν l in the regions from the liquid phases to the solid phases initially and evolves toward a single peak distribution in the liquid phase at ν target = 0.698. On the contrary, P(ν l ) at ν = 0.720 has the single peak (corresponding to the solid phase) initially and evolves toward the target distribution keeping the shape of the distribution as Gaussian-like. The result showed a positive correlation between local density fluctuations and relaxation time, as described in Sec. 3.1. The difference in the shape of the distribution comes from the spatial inhomogeneity in density, which might be one of the potential candidates of the origin of the anomalous slow relaxation.

Concluding remarks
In the present work, to elucidate the microscopic origin of equilibration in a hard disk system, we investigate the relaxing process toward the liquid states as a nonequilibrium response induced by the simple disturbance.
After the homogeneous expansion of the equilibrated systems around the Alder transition is induced, we estimate the relaxation time of four physical properties. As we do not expect, anomalous slow relaxation can be observed around ν = 0.708 in the co-existing phase for all properties. To identify this physical mechanism, we assumed that the spatial inhomogeneity of the initial equilibrated phases would contribute to the relaxation time. By visualization of spatial distribution and the probability distribution of the local density, it is confirmed that the two wellseparated co-existed phases cause the broad distribution with a plateau in density, which needs time to change the shape of distribution to a Gaussian-like single one. The inhomogeneity of the system might be considered as one of the physical factor(s) of anomalous slow relaxation in hard disk systems. Finally, we remark the relationship of inhomogeneity induced slow relaxation with the dynamical heterogeneity observed in super-cooled liquids near the glass transition [14,15], slower decay of granular gas with cluster impact in dense case [16] and the Mupemba paradox of fast freezing on warming (including inverse case) [17]. The simple model results would provide a useful example to consider the fundamental issues on slow decay observed in the various physical systems. . The spatial distributions' evolving process for both the local packing fraction and the orientational order parameter at ν = 0.708 (upper two lines) and ν = 0.720 (lower two lines) are shown, respectively. At ν = 0.720, the entire system is a homogeneous (crystal) phase like the liquid phase initially where the location of disorder (melting) patches seems to emerge homogeneously. On the contrary, at ν = 0.708, a partial domain of the crystallized phase exists at first and persists for a long time. These inhomogeneities due to the co-existing of two phases might cause a dominant contribution to the relaxation time of physical properties. Figure 4. The time evolution of the local density distribution at ν = 0.700, 0.708, and 0.720, respectively. The difference of shape between initial and final distribution might cause the potential substantial factor of slow decay.