Nonequilibrium steady states, coexistence, and criticality in driven quasi-two-dimensional granular matter

Nonequilibrium steady states of vibrated inelastic frictionless spheres are investigated in quasi-two-dimensional confinement via molecular dynamics simulations. The phase diagram in the density-amplitude plane exhibits a fluidlike disordered and an ordered phase with threefold symmetry, as well as phase coexistence between the two. A dynamical mechanism exists that brings about metastable traveling clusters and at the same time stable clusters with anisotropic shapes at low vibration amplitude. Moreover, there is a square bilayer state which is connected to the fluid by BKTHNY-type two-step melting with an intermediate tetratic phase. The critical behavior of the two continuous transitions is studied in detail. For the fluid-tetratic transition, critical exponents of $\tilde{\gamma }=1.73, {\eta }_4\approx 1/4$, and $z=2.05$ are obtained.

Figure 1

Sketch of geometrical objects defining the 2D local bond-orientational order parameter ${\psi }_4^{\left(n\right)}$. Particle centers are depicted as circles. Facets of the Voronoi tessellation (i.e., the perpendicular bisections of the connections of particle centers) are drawn as dashed, dotted, and dashed-dotted lines. Particles that share a Voronoi facet are considered nearest neighbors. Symbols and coloring of the sketch are described in the text.

Figure 2

(a) Phase diagram in the density-amplitude plane for box height $h=1.83\sigma$ and angular frequency $\omega {\tau }_0=12$. Investigated state points are on a grid with spacings $\mathrm{\Delta }\rho {\sigma }^2=0.02$ and $\mathrm{\Delta }A=0.002\sigma$. The simulation at each state point was carried out with $N=4000$ particles and averaged over a time interval $t=5000{\tau }_0$. The maximum density ${\rho }_{\text{max}}{\sigma }^2\approx 2.31$ is slightly larger than the density of two hexagonal close packed layers, $\rho {\sigma }^2=4/\sqrt{3}$, due to the possibility of buckling. Below the blue squares the fluidized state collapses and all particles drop to the bottom plate. The green crosses indicate evaporation of a threefold cluster upon decreasing $\rho$. Between $0.02\sigma \le A\le 0.04\sigma$ the nucleation density (gray filled triangles) of the threefold cluster upon compression differs from the evaporation density, with a hysteresis region (hatched) between the green crosses and gray filled triangles. The purple open triangles mark the density where the threefold cluster comprises the whole simulation box. The red stars and yellow circles indicate the continuous fluid-tetratic and tetratic-square transitions, respectively. Labels (c)–(j) relate the state points to the snapshots in the following panels. Lines are a guide to the eye. The inset shows a wider view of the phase diagram up to amplitude $A=1\sigma$ (where the dashed line marks the view of the main panel). (b) Sketch of the color coding of the particles in snapshots (c)–(j). The hue is determined by the complex phase $arg\left({\psi }_4^{\left(n\right)}\right)$ of the local order parameter. The color saturation indicates the modulus of ${\psi }_4^{\left(n\right)}$ ($|{\psi }_4^{\left(n\right)}|=0 \to$ gray, $|{\psi }_4^{\left(n\right)}|=1 \to$ fully saturated). (c)–(j) Top view snapshots of the simulation box at state points indicated and labeled in the phase diagram: (c) $\rho {\sigma }^2=1.2, A=0.004\sigma$ inelastic collapse (nonergodic). (d) $\rho {\sigma }^2=1.0, A=0.03\sigma$ fluid phase. (e) $\rho {\sigma }^2=1.2, A=0.1\sigma$ fluid-threefold coexistence (circular cluster shape). (f) $\rho {\sigma }^2=1.46, A=0.03\sigma$ fluid with square bilayer patches. (g) $\rho {\sigma }^2=1.5, A=0.03\sigma$ tetratic phase. The magnifier shows a connected path along grid lines, which follows two lattice sites up, two right, two down, and two left. The fact that the path does not close, as it would in a regular lattice, indicates a dislocation defect. The other dislocations are marked in the same way. (h) $\rho {\sigma }^2=1.7, A=0.03\sigma$ fluid-threefold coexistence (with the cluster wetted by square bilayer phase; triangular cluster shape). (i) $\rho {\sigma }^2=1.7, A=0.05\sigma$ fluid-threefold coexistence (partially wetted). (j) $\rho {\sigma }^2=1.7, A=0.03\sigma$ fluid-threefold coexistence with cluster percolating the box in $y$ direction [metastable; same state point as (h)]. $\phi$ denotes the angle between the $y$ axis and the symmetry axis of the lattice unit cell (see Sec. 5). The big arrow indicates the direction of motion of the cluster.

Figure 3

Order parameter ${\mathrm{\Psi }}_4$ as a function of density $\rho$ for $A=0.03\sigma$ and different box sizes $L$ as indicated. Lines are guides to the eye. ${\rho }_4$ indicates the fluid-tetratic critical density as fitted via Eq. (5). The arrow at $\rho {\sigma }^2=1.47$ marks the density where at $L=160\sigma$ crossover from short-range to quasi-long-range behavior is observed in the fourfold correlation function (cf. Fig. 6).

Figure 4

(a) Static structure factor $S_4$ as a function of wave number $k$ in the critical region exemplarily at densities $\rho$ indicated in the legend obtained at amplitude $A=0.03\sigma$ and box size $L=120\sigma$. Lines are fits to the Ornstein-Zernike form, Eq. (4), only valid while the finite system is fluid for $\rho {\sigma }^2\le 1.464$. (b) Orientational correlation length ${\xi }_4$ as a function of density $\rho$ as extracted from the data in (a). The blue line is a fit to Eq. (5) using data of the fluid regime (red symbols with error bars). Gray symbols without error bars are ordered states. (c) Static susceptibility ${\chi }_4$ as a function of ${\xi }_4$ as extracted from the data in (a). The blue line is a fit to Eq. (6) using data of the fluid regime (red symbols with error bars). Gray symbols without error bars are ordered states.

Figure 5

(a) Intermediate scattering function $F_4$ as a function of time difference $\tau$ at fixed wave number $k\sigma =0.052$ obtained at several densities $\rho$ as indicated, amplitude $A=0.03\sigma$, and box size $L=120\sigma$. Symbols show simulation data, lines show fits to Eq. (7). (b) Correlation time ${\tau }_4$ as a function of orientational correlation length ${\xi }_4$ obtained from the fits of $F_4(k,\tau )$ [see (a) ]. The line is a fit to Eq. (8) using data of the fluid regime (red symbols with error bars). Gray symbols without error bars are ordered states [cf. Fig. 4].

Figure 6

Fourfold correlation function $g_4$ as a function of distance $r$ on a double-logarithmic scale at densities ranging from $\rho {\sigma }^2=1.45$ to 1.5, amplitude $A=0.03\sigma$, and box size $L=160\sigma$. The black straight line corresponding to a power law with exponent –1/4 is a guide to the eye.

Figure 7

(a) Scan through the 2D direct correlation function $h$ along the direction $x$ of the order parameter ${\mathrm{\Psi }}_4$ at some exemplary densities as indicated and box size $L=120\sigma$. Multiple samples have been aligned and averaged over coherently, see main text. The black straight lines show exponential fits of the maxima of $h\left(x\right)$ to Eq. (12). The inset shows the resulting fit values of the positional correlation length ${\xi }_{\mathrm{pos}}$ for all obtained exponentials. (b) Maxima of the data shown in panel (a) on a double-logarithmic scale for a wider range of densities from $\rho {\sigma }^2=1.46$ to 1.70. For $\rho {\sigma }^2=1.54$ two curves are shown to illustrate non-ergodicity: the dashed and dashed-dotted lines were obtained from tetratic and square lattice initial configurations, respectively. Data shown for higher $\rho$ were obtained with square lattice initial configurations. The black straight line (power law $\sim r^{-1/4}$) is a guide to the eye, where the exponent $1/4$ is the upper bound for ${\eta }_{\mathrm{pos}}$ predicted by BKTHNY theory.

Figure 8

(a) Top view sketch of the threefold lattice consisting of two hexagonal layers of particles (not to scale). The rhombus in the center represents a unit cell of the lattice with its long diagonal being one of the three symmetry axes (dashed-dotted lines). The orientation of a cluster is measured by the angle $\phi$ between this diagonal and the $y$ axis of the box. Arrows at the edge show how top layer particles are supported by bottom layer particles from the exterior at the stable and unstable facets, respectively, and hence illustrate the asymmetry between these facets. In the lattice displayed here, $\phi =0deg$. This orientation has two equivalent left and right interfaces in configurations with interfaces parallel to the $y$ direction [as in Fig. 2], whereas a lattice rotated by $\phi =30deg$ would constitute the maximally asymmetric case with a stable facet to the right and an unstable facet to the left. The case $\phi =60deg$ is the same as $\phi =0deg$ mirror-inverted in $y$ direction and therefore again symmetric with respect to the $y$ axis. Therefore, it is sufficient to consider angles between $0deg$ and $30deg$. (b) Sketch of the calculation of the $x$ coordinate of the center of mass $x_{\mathrm{com}}$ with periodic boundary conditions via mapping of the $x$ coordinates of the particles onto the unit circle in the complex plane. Particles are displayed as red (light gray) circles; the average of the mapped coordinate is depicted as blue (dark gray) circle. Symbols are declared in the text.

Figure 9

(a) $x$ coordinate of the center of mass $x_{\mathrm{com}}$ as a function of time $t$ for individual runs at density $\rho {\sigma }^2=1.7$, amplitude $A=0.04\sigma$, and different cluster orientations $\phi$ as indicated. Black straight lines show end-to-end motion at each $\phi$ as used to extract the average velocity. (b) Cluster speed $v_{\mathrm{com}}$ of the center of mass as a function of angle $\phi$ obtained from the end to end distance of motion as displayed in (a) and averaged over four individual runs each. Lines are guides to the eye. By symmetry arguments, it is sufficient to consider the range $0deg\le \phi \le 30deg$, see inset. (c) Same as (a) but for fixed angle $\phi =20deg$ and different $A$ as indicated. (d) Same as main panel of (b) but for the center of mass motions displayed in (c).