# Phase transition in inelastic disks

## Abstract

This letter investigates the molecular dynamics of inelastic disks without external forcing. By introducing a new observation frame with a rescaled time, we observe the virtual steady states converted from asymptotic energy dissipation processes. System behavior in the thermodynamic limit is carefully investigated. It is found that a phase transition with symmetry breaking occurs when the magnitude of dissipation is greater than a critical value.

In energy conservative systems, it is well known that macroscopic properties at an equilibrium state can be described by a few state variables . This raise a question as to whether, in dissipative systems, similar state variables might be defined for macroscopic properties at steady states maintained by an external energy source. When a dissipation rate is sufficiently small, the system might be expected to behave as a conservative system; but is this actually the case? To examine this and other related questions, it is interesting to investigate models that can connect a conservative and a dissipative system by varying the system parameters.

An ensemble of elastic hard disks is one of the simplest models for describing the fluid state in conservative systems. Computer simulations played an important role in discovering the existence of fluid-solid transition in this model . We will here consider an ensemble of inelastic hard disks, which has previously been investigated as a model of granular materials. By varying the inelasticity, we can continuously change the system from a conservative to a dissipative system. This ability will be useful for our present purposes of investigating a thermodynamic properties in dissipative systems.

In order to attain a steady state in a dissipative system, an external energy source that compensates for energy dissipation due to inelastic collisions is indispensable. In granular physics, the vibrating bed  is a typical energy source. Nevertheless, attaching such an energy source will break the isotropy of the system from the outset. In this letter, we simulate the dissipative system without energy sources and investigate the virtual steady states using a new observation frame that will be introduced later.

The inelastic disk system without energy input , known as a cooling gas or dissipative gas system, has been investigated using molecular dynamics (MD) and hydrodynamic models. These studies reveal that the homogeneous state is unstable for sufficiently high inelasticity even without attaching any energy source. In this state, a cluster of particles  and collective mean flow (shear mode) are formed. We study how the instability is characterized by means of MD simulation of an ensemble of inelastic hard disks, with a particular focus on the asymptotic collective phenomena and its system size dependence.

The system consists of an ensemble of $`N`$ inelastic hard disks in two dimensional space. For simplicity, each particle is assumed to have a unit mass and its rotational degrees of freedom are ignored. Collisions between circular particles are inelastic; the inelastic collision is implemented in the following manner. At a collision of two particles, the tangential velocities to the collision plane are preserved, while the normal component of the relative velocity $`\mathrm{\Delta }v_n`$ changes to $`\mathrm{\Delta }v_n^{}`$, where

$$\mathrm{\Delta }v_n^{}=e_n\mathrm{\Delta }v_n.$$
(1)   
Here, the coefficient of restitution $`e_n(0e_n1)`$ is constant for all collisions. In the following, we parameterize the inelasticity by $`ϵ(1e_n)`$. In the case $`ϵ=0`$, the system becomes conservative with the elastic collisions; in the case $`ϵ>0`$, energy dissipation occurs due to the inelastic collisions.

The selected boundary conditions consist of a square box enclosed by elastic rigid walls. Thus the energy of the system is not dissipated by the bouncing of particles against the wall. For initial conditions, we adopt spatially homogeneous states equilibrated by setting $`ϵ=0`$. The time evolution of the system is calculated by the event-driven method.

In this letter, we focus on asymptotic behavior of the system after a sufficiently long transient time. We note that our model is almost identical with the model in  except for boundary conditions.

Under our boundary conditions, there are three important system parameters: the number of particles $`N`$, the inelasticity $`ϵ`$ and the area fraction of particles $`\varphi `$. Here $`\varphi (0\varphi <1)`$ is the ratio of the total area covered by particles to that of the box $`S`$.

In the conservative case $`ϵ=0`$, it is known that this system, consisting of an ensemble of elastic hard disks, has two phases: a fluid and a solid phase. In this letter, we investigate the system in the range of parameter $`\varphi `$ corresponding to the fluid phase when $`ϵ=0`$. For this parameter range, our results do not depend on the dispersity of particle radii $`a`$. In the following, we show only the results for monodisperse case($`a=0.5`$); however, the results of our simulations for the polydisperse case (uniform distribution in the range $`0.4a0.5`$) are almost the same. On the other hand, in the range of parameter $`\varphi `$ corresponding to the solid phase when $`ϵ=0`$, the results depend on the dispersity of particle radii. The results for this parameter range will be reported in a future publication.

For sufficiently large values of $`ϵ`$, it is known that the time development by the event driven method is ill-defined because infinite collisions occur in a finite time interval (inelastic collapse) . Our investigation is limited to sufficiently small values of $`ϵ`$ in order to avoid the inelastic collapse.

We will now describe how to observe the system. Since no energy source is attached to the system, the total energy of the system monotonically decreases in time. Nevertheless, by the rescaling of time, it is possible to make the energy of the system virtually conserved. We introduce the rescaled time $`\stackrel{~}{t}`$ and the rescaled velocity of $`i`$-th particle $`\stackrel{~}{v}_i`$:

$$\stackrel{~}{t}=_0^t\sqrt{T(t)}𝑑t,\stackrel{~}{v}_i=v_i(t)/\sqrt{T(t)},$$
(2)   
where $`T(t)`$ is the averaged kinetic energy per one particle at a time $`t`$, $`T(t)=_iv_i^2(t)/2N`$; and $`v_i(t)`$ is the original velocity of $`i`$-th particle. We call the rescaled system “R-system” and the original system “O-system”. All of the present observations are carried out for the “R-system”. Under the translation by Eqs.(2), the averaged kinetic energy per one particle in the “R-system” is normalized to the unity, $`\stackrel{~}{T}(t)_i\stackrel{~}{v}_i^2(t)/2N=1`$, i.e., the total kinetic energy is conserved for any $`ϵ`$ in the “R-system”. Further, asymptotic energy dissipation processes in the “O-system” are translated to steady states in the “R-system”.

It must be noted here that the translation by Eqs.(2) is just the rescaling of time and that the trajectories of the particles in space are not influenced by this rescaling. We calculate the time evolution of the “O-system” by the event driven method, and observe the “R-system” with a special focus on its steady states after a sufficiently long transient time. In the rest of this letter, all variables shown refer to those in the “R-system”.

The dependencies of pressure on $`\varphi `$ for several values of $`ϵ`$ are shown in Fig.1. Pressure $`P`$ is defined as the time averaged sum of the impulses at the bouncing of particles on the walls per unit length per unit time in the “R-system”. The vertical axis in the figure is $`N\stackrel{~}{T}/PS`$, which is unity for the ideal gas limit ($`\varphi 0`$). Here $`S`$ is the area of the box.

In the conservative case $`ϵ=0`$, only the fluid phase exists in the system for $`\varphi <\varphi _c`$. The value of $`\varphi _c`$ is reported to be $`\varphi _c0.7`$ . Consider the dissipative cases $`ϵ>0`$ in Fig.1(a) paying attention to the difference from the case $`ϵ=0`$. When $`ϵ=0.02`$, there is no remarkable difference. For $`ϵ=0.04`$, a decrease in pressure is observed in the intermediate range of $`\varphi `$. For $`ϵ=0.08`$, a similar difference appears in the wider range of $`\varphi `$.

Comparing Fig.1(a) and (b), it is found that similar dependencies of pressure on $`\varphi `$ are observed for the same values of $`Nϵ`$. Thus $`Nϵ`$ is an important parameter of the system and it would characterize the distance from the conservative system.

For sufficiently large $`Nϵ`$, the behavior of the system is quite different from the case $`ϵ=0`$. A snapshot of the system for $`(N,ϵ)=(256,0.08)`$ is shown in Fig.2(b) compared to that for $`ϵ=0`$ (Fig.2(a)).

From Fig.2(b), it is found that a correlation of particle velocity exists and that a mean flow circulating anti-clockwise in the box is formed. As a result of the emergence of the circulating mean flow, the impulse of particles to the walls at bounce decreases compared with the case $`ϵ=0`$. Then the decrease of pressure, i.e., the increase of $`N\stackrel{~}{T}/PS`$, occurs. For larger values of $`ϵ`$, the circulation is formed in a wider range of $`\varphi `$, as seen in Fig.1.

Since the circulation is formed above a certain threshold value of $`ϵ`$, normalized total angular momentum $`M`$ of the system is defined as an order parameter of the system.

$$M=\frac{1}{N\sqrt{S}}\underset{i=1}{\overset{N}{}}(r_iR_0)\stackrel{~}{v}_i.$$
(3)   
Here $`r_i`$ and $`\stackrel{~}{v}_i`$ are coordinates and velocities of $`i`$-th particle, $`R_0`$ are the coordinates of the center of the box, “$``$” refers to outer product, and the notation $`\stackrel{~}{v}`$ is used to remind us that we are observing the “R-system”. $`M`$ is normalized by the box length $`\sqrt{S}`$ to eliminate the dependency on the box size.

Time developments of $`M`$ for some values of $`ϵ`$ are shown in Fig.3, where $`N=256`$ and $`\varphi =0.5`$. For $`ϵ=0`$, the value of $`M`$ fluctuates around $`M=0`$. As $`ϵ`$ is increased, the fluctuations of $`M`$ increase. For $`ϵ=0.036,0.038`$, the circulating mean flow is formed. The direction of the circulation sometimes turns over. As $`ϵ`$ increases further, the events of turning over become rare. In order to investigate the formation of the circulation in detail, the distributions of $`M`$ are shown in Fig.4; because these distributions are symmetric around $`M=0`$, only the region $`M0`$ is shown.

As denoted by Fig.4, it is found that the peak position of $`M`$, $`M_{\mathrm{peak}}`$ becomes nonzero for $`ϵ`$ above a certain threshold value and increases continuously from zero with increasing $`ϵ`$. We also observe a similar behavior when we vary the value of $`\varphi `$ while fixing the value of $`ϵ`$. Thus the circulation appears continuously for any direction in the parameter space($`ϵ,\varphi `$). In Fig.4, the distributions of $`M`$ are broad in shape, which indicates that the value of $`M`$ fluctuates around the value at the peak $`M_{\mathrm{peak}}`$. If the formation of the circulation is a phase transition accompanied by symmetry breaking in the thermodynamic limit ($`N\mathrm{}`$), the width of the distribution will converge to zero in this limit. In order to confirm this scenario, the distributions of $`M`$ are examined in the dependence on $`N`$ in Fig.5, where $`N`$ is varied while $`Nϵ`$ is kept constant.

In Fig.5(a), it is found that $`M_{\mathrm{peak}}`$ will converge in the limit, $`N\mathrm{}`$. This convergence shows that $`M_{\mathrm{peak}}`$ is a function of $`Nϵ`$ for sufficiently large $`N`$, because the value of $`Nϵ`$ is fixed in Fig.5(a).

From Fig.5(b), it is found that the width of the fluctuation of $`M`$ decreases in the manner of $`1/\sqrt{N}`$ as $`N`$ increases. Thus, we conclude that the circulation appears as a phase transition with symmetry breaking. The turning over events of the circulation observed in Fig.3 are finite size effects.

In Fig.6, $`M_{\mathrm{peak}}`$ around the critical point is shown. The horizontal axis is $`Nϵ(13/\sqrt{N})`$, where the second term $`O(1/\sqrt{N})`$ comes from finite size effects. This figure clearly shows that $`M_{\mathrm{peak}}`$ is a function of $`Nϵ`$ for sufficiently large $`N`$. Further calculations are needed for the precise determination of the critical exponents.

In this letter, we simulate inelastic disks in a square box and investigate the steady states of the system using a new observation frame wherein the energy of the system is virtually conserved by the rescaling of time. The important parameter of the system is $`Nϵ`$, where $`N`$ is the number of particles and $`ϵ`$ is inelasticity. At $`ϵ=0`$, the steady states are homogeneous equilibrium states. At sufficiently small $`Nϵ`$, the steady states are also homogeneous. For $`Nϵ`$ above the critical points, the homogeneous steady states are no longer stable . Then the phase transition with symmetry breaking occurs where the circulation appears continuously.

The results shown in this letter are independent of the characteristics of the model. As noted before, the results are independent of the dispersity of the particle radii. The hard-core potential is not essential to the results. We have confirmed that the similar results are obtained in the case of soft-core potential: a simulation of the models using the discrete element method with Nose-Hoover thermostat without the rescaling of time. We have also observed similar circulating mean flow in the system with inelastic hard spheres. Thus similar results would also be obtained in three dimensions.

A number of problems remain to be investigated. For example, for larger $`ϵ`$ up to $`ϵ=1`$, are there one or more additional phases? Is hydrodynamic description by $`N\mathrm{}`$ with fixing $`Nϵ`$ really possible for any boundary conditions? Do fluid- and solid phases exist in $`ϵ0`$? How are the phenomena in the system investigated here related to those in the system with a real energy source, such as vibrating beds?

It is noteworthy that the homogeneous states in the system with finite $`ϵ`$ are always unstable for sufficiently large $`N`$, because the phase transition occurs at finite $`Nϵ`$. Thus the conservative system is a singular-limit system when we take the limit $`N\mathrm{}`$ first.

The author thanks S. Sasa and N. Nakagawa for their useful input and careful reading of the manuscript, N. Ito and Y.-h. Taguchi for their helpful advice, and Dept. of Math. Sci. at Ibaraki Univ. for their hospitality. The author also acknowledges the support from JSPS.