Entropic long-range ordering in an adsorption-desorption model

Entropic long-range ordering in an adsorption-desorption model Adam Lipowski1 and Dorota Lipowska2 arXiv:1905.04802v2 [cond-mat.stat-mech] 11 Jun 2019 2 1 Faculty of Physics, Adam Mickiewicz University, Poznań, Poland Faculty of Modern Languages and Literature, Adam Mickiewicz University, Poznań, Poland We examine a two-dimensional nonequilibrium lattice model where particles adsorb at empty sites and desorb when the number of neighbouring particles is greater than a given threshold. In a certain range of parameters the model exhibits entropic ordering similar to some hard-core systems. However, contrary to hard-core systems, upon incresing the density of particles the ordering is destroyed. In the heterogenous version of our model, a regime with slow dynamics appears, that might indicate formation of some kind of glassy structures. I. INTRODUCTION Models with hard-core interactions serve as an idealization of a number of important physical systems. Indeed, various aspects of liquids [1], glasses [2], liquid crystals [3], or certain adsorbates [4] were successfully examined using hard-core models mainly by means of numerical simulations. Studies of particular importance are those related to the emergence of long-range ordering such as, for example, freezing of hard spheres [5] or of hard disks [6], which proceeeds via an intermediate hexatic phase. Let us emphasize that hard-core interactions render the temperature irrelevant and ordering in such systems is of purely entropic origin [7] with coverage (or pressure) as a control parameter. In the computationally less demanding lattice hard-core systems, some more detailed insight into the ordering process is available. For example, on a square lattice, and when hard-core exclusion prevents nearest neighbors of a given particle from being occupied, the ordering transition turns out to belong to the Ising model universality class [8, 9], which is related to the double degeneracy of the ordered phase. When nearest- and next-nearest-neighbor repulsions are present, a four-fold degenerate columnar order is formed. Although it is more difficult to establish the nature of the ordering transition in this case [10], most works suggest the Ashkin-Teller universality class [11]. Having in mind formation of some adsorbate structures such as, for example, He on graphite [12] or on graphyne [13], or H on W [14, 15], or O on Pt [16], we should take into account that equilibrium hard-core models provide only a very approximate description of these complex physical phenomena [17]. An important process, which often accompanies adsorption and affects, for example, a surface diffusion [18] or various surface chemical reactions, is desorption [19, 20]. In certain statistical mechanics studies, the role of desorption in some equilibrium as well as nonequilibrium hard-core systems has already been examined [9]. In the present paper, we describe a nonequilibrium model where the desorption rather than the hard-core exclusion plays the primary role in the formation of an entropy stablized long-range order. What is, in our opinion, interesting is that the resulting ordered structures, and perhaps accompanying phase transitions, are the same as those in the hard-core systems but the nature of the ordering process is much different: the ordered structures are destroyed when the density of particles increases (not decreases, as in hard-core systems). Our work thus suggests that an alternative mechanism may play the role in the formation of entropy stabilized long-range ordering. II. MODEL We examine a collection of particles, which adsorb at a two-dimensional surface, but when a particle gets surrounded by too many neigbouring particles, it desorbes. Thus, in a statistical mechanics fashion, in our model we have N particles distributed (without overlaps) over sites of a square lattice of linear size L with periodic boundary conditions. In an elementary step of the dynamics of our model, one selects randomly a particle and if it is unstable, it is relocated to one of the randomly selected empty sites. A particle is considered unstable if the number of particles on neighboring sites is greater than a given value k. Some of the model characteristics are time dependent and we define the unit of time (1 MC step) as N elementary steps (one step per particle). Let us emphasize that dynamics of our model shifts a constant number of particles (N ) and desorption is always followed by adsorption which are thus not independent processes. Such an approach bears some resemblance to the method of constant coverage ensemble used in the context of some surface-reaction models [21]. What is more important, is the lack of detailed balance in our model since a stable particle has a zero probability of desorption. Dynamics of our model might thus get trapped in an absorbing state where each particle is stable. In general, it is impossible to describe the stationary state of such models in terms of equilibrium Gibbs distributions and they belong to the realm of nonequibrium statistical physics. Models of this kind include some versions of the contact process [22–25], but might describe also some adsorption-desorption systems [26, 27]. The density of particles ρ = N/L2 and the parameter k thus control the behavior of the model. Of course, when the density ρ is sufficiently small, after a short transient each particle finds a stable position surrounded by at most k neighbors. Such a state is an absorbing 2 0 -0.5 -1 -1.5 log10(ρa) state of the dynamics. Analogously, when ρ is large, relocated particles are unlikely to find stable positions that, in addition, do not destabilize their neighbors, and consequently, a fraction of particles are constantly reshuffled. More interesting, and less obvious, is the behavior in an intermediate density range. To examine it in more details, we carried out Monte Carlo simulations. We used two types of neighborhoods: (i) nearest neighbors (4 sites) and (ii) nearest and next-nearest neighbors (8 sites), and the results we obtained are presented in the next sections. -2 0.2 0.22 0.221 0.222 0.223 0.224 0.225 0.23 0.3 0.45 0.55 -2.5 -3 -3.5 -4 -4.5 0 III. NEAREST-NEIGHBOR INTERACTIONS In this case each site has four neighbors, which implies that 0 ≤ k ≤ 4. To introduce the methodology, we first examine k = 0. After generating a random initial configuration with N = ρL2 particles, we redistribute them using the model dynamics. We calculated the average density of active (i.e., unstable) particles ρa as a function of time t, and the results are presented in Fig. 1. As expected, for small ρ (ρ ≤ 0.221), we observe a fast decay of ρa and eventually the model reaches an absorbing state of the dynamics, where nearest neighbors of each particle are empty. On the other hand, for ρ ≥ 0.223, the model remains in a state with a finite fraction of active particles. Let us notice that various periodic structures would satisfy the condition k = 0 and the densest of them (with ρ = 0.5) is the checkerboard ordering, analogous to the hard-core model with the nearest-neighbor exclusion [8, 9]. We do not observe formation of any of such global periodic structures and their dynamical creation is apparently unlikely. However, ordering might appear on a small scale. For example, for ρ = 0.225 (Fig. 2), the snapshot configuration shows various clusters with a checkerboard pattern. Such clusters, however, are very unstable. If one of the empty sites is chosen and gets occupied, both this site as well as its four neighbors turn into unstable sites (Fig. 3). As a result, large clusters of this kind do not form. For k = 2, a much different scenario takes place. In this case, the transition between absorbing and active regimes of the model takes place around ρ = 0.5035(10) (Fig. 4). While the absorbing state, similarly to the case k = 0, is disordered, the active regime is different. Indeed, for ρ = 0.51 formation of long-range ordered structures is clearly seen in Fig. 5. Let us notice that a checkerboard structure satisfies actually the stronger limit k = 0, since in this case the number of occupied nearest neighbors is zero. It implies a stronger stability of such structures: an empty site that gets occupied does not destabilize its neighbors (Fig. 6). One can easily find higher-density periodic structures that satisfy the limit k = 2 (Fig. 7), but they are not dynamically stable (as is the checkerboard structure in the k = 0 case) and we did not observe their formation during the evolution of our model. Upon increasing the density ρ, the number of unsta- 0.5 1 1.5 2 log10(t) 2.5 3 3.5 4 FIG. 1. Time dependence of the density of active particles ρa for the nearest-neighbor model with k = 0 and ρ = 0.2 (bottom curve), 0.22,..., 0.55 (top curve). Simulations were run for L = 103 and the results are averaged over 100 independent runs. The statistical error is of the order of noise seen in the plotted curves. ble particles ρa in the steady state also increases, wich gradually destroys the long-range ordering. To examine the process in more detail, we carried out simulations for ρ > 0.5, which started from the predefined checkerboard ordering. We divided the lattice into two sublattices, A and B, and placed L2 /2 particles on the sublattice A, while the remaining ones were randomly distributed on the sublattice B. Running the model dynamics, we relaxed the system until it reached the steady state and then we measured the order parameter m defined as ! X 1 X m= 2 (1) ni , ni − L i∈A i∈B where ni = 0 or 1 for a site i being empty or occupied by a particle, respectively. The results (Fig. 8) show that m decays to 0 at ρ = ρc ≈ 0.556. Assuming the power-law decay (m ∼ (ρc − ρ)β ), we estimate β ≈ 0.17(5), and the fit is based on data close to the critical point (0.555 < ρ < 0.556). Taking into account that the checkerboard structure is double-degenerate, one might expect that the transition at ρc belongs to the Ising model universality class and the obtained estimate of β is marginally consistent with the Ising model value 0.125. We also measured the variance χm of the order parameter m in the ρ > ρc regime (simulations started from a random initial configuration). The results (inset in Fig. 8) show that χ has a power-law divergence χm ∼ (ρ − ρc )−γ with γ ≈ 1.75, which is in a very good agreement with the Ising model value. Presented numerical results are obtained for the system size L = 103 . Except the very vicinity of the critical point (ρ = ρc ) the examined systems seem to be sufficiently large and the finite size effects are negligable. More precise esti- 3 0 100 -0.5 -1 log10(ρa) -1.5 50 -2 -2.5 0.49 0.495 0.5 0.503 0.504 0.505 0.51 0.55 0.7 -3 -3.5 -4 -4.5 -5 0 0.5 1 1.5 2 2.5 log10(t) 3 3.5 4 4.5 0 0 50 100 FIG. 2. The distribution of particles in the stationary state (after relaxation of a random initial configuration for 104 MC steps) for the nearest-neighbor model with k = 0 and ρ = 0.225. FIG. 4. Time dependence of the density of active particles ρa for the nearest-neighbor model with k = 2 and ρ = 0.49 (bottom curve), 0.495,..., 0.7 (top curve). Simulations were run for L = 103 and the results are averaged over 100 independent runs. 100 50 0 FIG. 3. In the nearest-neighbor model with k = 0, when an empty site in a cluster with a checkerboard structure gets occupied, it makes unstable all its four neighbors. A subsequent move is likely to erode the cluster ordering. 0 50 100 FIG. 5. The distribution of particles in the stationary state (after relaxation of a random initial configuration for 104 MC steps) for the nearest-neighbor model with k = 2 and ρ = 0.51. mations of the critical behaviour would certainly require more systematic analysis of finite size effects. We carried out some simulations for k = 1 and k = 3, and we did not observe formation of a long-range ordering. The behavior of the model in these cases seems to be similar to the k = 0 case. As a final remark in this section, let us notice that the decay of ordering via the Ising-like phase transition takes place upon the density increase. This is opposite to the behavior of hard-core systems, where the high-density phase is long-range ordered. IV. NEXT-NEAREST-NEIGHBOR INTERACTIONS We also carried out simulations for the model with the nearest- and next-nearest-neighbor interactions, where each site has 8 such neighbors. The formation of longrange ordering was observed only for k = 4 and densities 4 100 FIG. 6. In the nearest-neighbor model with k = 2, the empty site in the checkerboard structure that gets occupied is the only unstable site. The inital configuration is likely to be restored unless some other nearby empty site gets occupied. 50 0 0 FIG. 7. Periodic structure with the density ρ = 7/12, where particles satisfy the limit k = 2 in the nearest-neighbor version, and k = 4 in the next-nearest-neighbor version. greater, but not much, than 0.5. The snapshot configuration for ρ = 0.51 (Fig. 9) clearly shows the formation of the columnar ordering. To examine such structures in more detail, we ran simulations with initial configurations with a predefined columnar ordering, similarly to the nearest-neighbor k = 2 version. We measured the columnar order parameter l, which basically counts the 0.5 0.45 0.4 0.35 0.25 0.2 log10(χ) m 0.3 0.15 0.1 0.05 3 2.7 2.4 2.1 1.8 1.5 1.2 0.9 0.6 -2.7 γ=1.75 -2.4 -2.1 -1.8 log10(ρ-ρc) -1.5 0 0.5 0.51 0.52 0.53 0.54 0.55 0.56 ρ FIG. 8. The order parameter m as a function of the density ρ for the nearest-neighbor model with k = 2 (L = 103 ). The inset suggests that the variance of the order parameter χm diverges at the critical point ρ = 0.556 with the exponent γ = 1.75. 50 100 FIG. 9. The distribution of particles in the stationary state (after relaxation of a random initial configuration for 104 MC steps) for the nearest- and next-nearest-neighbor model with k = 4 and ρ = 0.51. number of sites with horizontally versus vertically placed neighbors: l= 1 X li , L2 i (2) where li = 1 (or −1) for the site i, which has its two horizontal (or vertical) neighbors occupied (otherwise li = 0). Our numerical results show (Fig. 10) that similarly to the nearest-neigbour version, l takes the maximum value 1/2 (perfect columnar ordering) at ρ = 0.5. When ρ increases, the number of unstable particles also increases, which gradually destroys an ordering. The least-square fitting to the numerical data close to the transition point gives ρc = 0.5207(5) and β = 0.25(5), and the fit was made using data for 0.52 < ρ < 0.5207. Moreover, from the behavior of the variance χl of the order parameter (inset in Fig. 10), we estimate γ = 1.07(3). The four-fold degeneracy of the columnar ordering suggests that, similarly to some hard-core systems with a columnar ordering [11, 28], the critical behavior of our model may belong to the Ashkin-Teller universality class. In such a case, one expects β = 1/12 and γ = 7/6 [29], and our estimate of γ is very close to the expected value. The deviation of β might be related to strong finite-size effects or the fact that the true asymptotic regime was not yet reached in our simulations. More detailed analysis would be clearly desirable. 5 0.5 -0.5 0.45 -1 0.4 2 0.2 0.15 γ=1.07 log10(ρa) log10(χ) l 0.3 0.25 -1.5 2.4 0.35 1.6 1.2 0.32 0.34 0.35 0.355 0.36 0.37 0.38 0.4 0.45 0.48 -3.5 0 -3.6 0.05 -2.5 -3 0.8 0.4 0.1 -2 -3.2 -2.8 -2.4 log10(ρ-ρc) -2 -4 -1.6 -4.5 0 0 0.5 0.505 0.51 0.515 0.52 1 2 0.525 ρ FIG. 10. The columnar order parameter l as a function of the density ρ for the next-nearest-neighbor model with k = 4 (L = 103 ). The inset shows that the variance of the order parameter χl diverges at the critical point ρ = 0.5207 with the exponent γ = 1.07. V. NEAREST-NEIGHBOR INTERACTIONS WITH HETEROGENEITIES Heterogeneity of size, mass, or shape of particles is known to play an important role in hard-core systems. For example, it might lead to the phase separation of different particles [30] or to the formation of multiple glassy phases [31, 32]. Studying analogous phenomena in lattice models, which are usually computationally more tractable, might provide a valuable insight into the role of space dimension, range of interactions or symmetries. Despite decades of intensive research, the formation of a glassy state is a particularly challenging problem. While its existence is well documented in the three-dimensional systems [2, 33], the status of a two-dimensional glass is not certain. Although some works report certain dynamical glassy features in two-dimensional systems [34], some other question the existence of a glassy transition in such systems [35, 36]. It is not our objective to address these important general questions but rather to show that a heterogeneous version of our model, which may mimick the bi- or polydisperse hard-core systems, develops some slowly evolving characteristics which could suggest some relations with glassy systems. In particular, we examine a heterogeneous version of our nearest-neighbor model, where a fraction p of particles obeys the dynamical rule with k = 2, and the remaining fraction (1 − p) with k = 0. Simulations for p = 0.9 show that k = 0 particles hinder reaching the absorbing state and an active regime extends up to ρ ∼ 0.4 (Fig. 11). Moreover, the absorbing regime seems to be separated into two sub-regimes. For lower densities (ρ = 0.32, 0.34), an ordinary, fast (presumably exponential) decay of the density of active sites ρa can be seen. However, for larger densities (ρ = 0.35 ∼ 0.37), a much 3 log10(t) 4 5 6 FIG. 11. Time dependence of the density of active particles ρa for the heterogeneous nearest-neighbor model with k = 2(90%), and k = 0(10%) and for ρ = 0.32 (bottom curve), 0.34,..., 0.48 (top curve). Simulations were run for L = 300 and the results are averaged over 100 independent runs. slower decay of ρa can be clearly seen. For example, for ρ = 0.355 from the estimation of the asymptotic slope of the numerical data, we obtain ρa ∼ t−0.25 . To examine in more detail the structure of the model, we calculated the time-dependent variance χm of the order parameter (1) (Fig. 12). For a moderately large density of particles (ρ = 0.45 ∼ 0.48), the variance χm rapidly increases in time, which indicates formation of long-range ordered checkerboard structures. For larger density (ρ = 0.51), the density of active particles ρa is too large, which destroys a long-range order and χm saturates at a finite value. In the homogeneous case (p = 1), we did not calculate the time dependent χm but an analogous behavior would be observed. Also similarly to the homogeneous case, at a low density of particles (ρ = 0.34), the variance χm saturates at a small value, which indicates an absence of long-range ordered structures. Less evident is the behavior for intermediate densities (ρ = 0.35 ∼ 0.37), where a noticeable but slow increase of χm can be seen. It may indicate a very slow growth of domains, thus providing a further evidence that in this regime the model exhibits some glassy characteristics. VI. CONCLUSIONS In the present paper, we introduced a simple adsorption-desorption model that may generate an entropic long-range ordering. The structures formed (checkerboard, columnar) bear some similarity to the entropic order in hard-core systems, but the mechanism that generates ordering in our model is much different. In particular, the ordered phase exists for a certain intermediate particle density and gets destroyed upon a density increase—not upon decrease as in ordinary hard-core systems. The order-disorder transitions are likely to be- 6 Simulations show that in the nearest-neighbor version, where each site has four neighbors, the checkerboard ordering appears for k = 2. In the next-nearest-neighbor version with 8 neighbors, we found the columnar ordering for k = 4. One of the questions is why an ordering appears only for k equal to the half of the number of neighbors. It might be related to the fact that in both cases the ordered phase originates at (or very close to) ρ = 0.5, but let us notice that for both types of the ordering, even the stronger stability criterion is satisfied and particles have less than k occupied neighbors. One might thus imagine that, in principle, the dynamics with smaller k could reproduce such an ordering as well, but apparently this is not the case. Most likely, for smaller k, ordered structures are not dynamical attractors of the model, however, more convincing evidence of such a scenario would be desirable. Having in mind some adsorbing systems, it would be certainly more realistic to move an unstable particle with a diffusive dynamics rather than place it on a randomly selected empty site. Preliminary calculations (results will be presented elsewhere) show, however, that a similar long-range ordering should form also in such a version. The present dynamics, where a desorbed particle hopes to randomly chosen empty site, might be more suitable in the context of some evaporation-recondensation systems. In such a case, however, temperature dependent effects should be taken into account. Another issue, which in our opinion is worth further studies, is a slow dynamics in a heterogeneous version of our model with particles having different values of k. Simulations show that in such a case a new regime appears, where the evolution toward the absorbing state is very slow (but power-law). One might hope that a mixture of particles with different values of k or with different ranges of interactions will lead to even slower dynamics, which would be more relevant in the context of glassy systems. A glassy state often appears in various hard-core systems and a heterogeneity (e.g., polydispersity) is known to enhance it. The dynamics of our model in some cases exhibits a considerable slowdown and its further studies may contribute to a better understanding of a somewhat unclear status of glassy dynamics in two-dimensional systems. [1] W. G. Hoover and Francis H. Ree, Melting transition and communal entropy for hard spheres, J. Chem. Phys. 49, 3609 (1968). [2] L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, Rev. Mod. Phys. 83.2, 587 (2011). [3] D. Frenkel, Structure of hard-core models for liquid crystals, J. Phys. Chem. 92.11, 3280-3284 (1988). [4] I. N. Yakovkin, Monte Carlo simulation of order-disorder transitions in linear structures of alkali and alkaline earth adsorbates, Surf. Sci. 282.1-2, 195-201 (1993). [5] B. J. Alder and T. E. Wainwright, Phase transition for a hard sphere system, J. Chem. Phys. 27, 1208 (1957). [6] A. Jaster, The hexatic phase of the two-dimensional hard disk system, Phys. Lett. A 330.1-2, 120-125 (2004). [7] R. Dickman, Ising meets Ornstein and Zernike, Debye and Hückel, Widom and Rowlinson, and others, Braz. J. Phys. 30.4, 711-719 (2000). [8] G. Kamieniarz and H. W. J. Blöte, The non-interacting hard-square lattice gas: Ising universality, J. Phys. A 26, 6679 (1993). [9] D.-J. Liu and J. W. Evans, Ordering and percolation transitions for hard squares: Equilibrium versus nonequilibrium models for adsorbed layers with c(2x2) superlattice ordering, Phys. Rev. B 62, 2134 (2000). [10] H. C. M. Fernandes, J. J. Arenzon, and Y. Levin, Monte Carlo simulations of two-dimensional hard core lattice gases, J. Chem. Phys. 126, 114508 (2007). [11] K. Ramola, K. Damle, and D. Dhar, Columnar order and Ashkin-Teller criticality in mixtures of hard squares and 4 0.34 0.35 0.355 0.36 0.37 0.38 0.4 0.45 0.48 0.51 3.5 3 log10(χ) 2.5 2 1.5 1 0.5 0 -0.5 -1 0 1 2 3 log10(t) 4 5 6 FIG. 12. Time dependence of the variance of the order parameter χm (t) for the heterogeneous nearest-neighbor model with k = 2(90%), and k = 0(10%) and for ρ = 0.34 (bottom curve), 0.35,..., 0.51 (top curve, that bends horizontally around t ∼ 103 ). Simulations were run for L = 300 and the results are averaged over 100 independent runs. . long to the expected universality classes. For the double degenerate checkerboard ordering, the transition belongs to the Ising model universality class, and for the columnar ordering, our data are marginally consistent with the Ashkin-Teller universality class. Precise estimation of critical exponents in our model would certainly require more careful analysis of statistical errors and of finite size effects. The main objective of the present paper is, however, to present a new mechanism of entropic ordering in lattice models and its further analysis is postponed for the future. 7 dimers, Phys. Rev. Lett. 114, 190601 (2014). [12] W. E. Carlos and W. C. Milton, Interaction between a He atom and a graphite surface. Surf. Sci. 91.1, 339-357 (1980). [13] J. Ahn et al., Symmetry-changing commensurateincommensurate solid transition in the He 4 monolayer on 6, 6, 12-graphyne, Phys. Rev. B 99.2, 024113 (2019). [14] G. Doyen, G. Ertl, and M. Plancher, Order–disorder phenomena in adsorbed layers described by a lattice gas model, J. Chem. Phys. 62, 2957 (1975). [15] N. V. Petrova, I. N. Yakovkin, and Yu G. Ptushinskii, Monte Carlo simulations of hydrogen adsorption on the W (110) and Mo (110) surfaces. Eur. Phys. J. B 38.3, 525-531 (2004). [16] A. P. J. Jansen, Island formation without attractive interactions, Phys. Rev. B 77.7, 073408 (2008). [17] G. D. Garcı́a et al., Monte Carlo study of binary mixtures adsorbed on square lattices. Surf. Sci. 606.1-2, 8390 (2012). [18] J.-S. Wang, P. Nielaba, and V. Privman, Locally frozen defects in random sequential adsorption with diffusional relaxation, Physica A 199.3-4, 527-538 (1993). [19] D. G. Vlachos, L. D. Schmidt, and R. Aris, The effects of phase transitions, surface diffusion, and defects on surface catalyzed reactions: Fluctuations and oscillations. J. Chem. Phys. 93.11, 8306-8313 (1990). [20] G. Ertl, Non-linear dynamics: Oscillatory kinetics and spatio-temporal pattern formation, in Handbook of Heterogeneous Catalysis, edited by G. Ertl et al. (Wiley-VCH Verlag GmbH & Co. KGaA., 2008), 1492-1516. [21] R.M.Ziff and B. J. Brosilow, Investigation of the firstorder phase transition in the A-B 2 reaction model using a constant-coverage kinetic ensemble. Physical Review A 46 4630 (1992). [22] T. Tomé, and R. Dickman, Ziff-Gulari-Barshad model with CO desorption: An Ising-like nonequilibrium critical point. Physical Review E 47, 948 (1993). [23] R. Durrett, Stochastic spatial models. SIAM review bf 41, 677 (1999). [24] D.-J.Liu, C.-J. Wang, and J. W. Evans, Discontinuous [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] Phase Transitions in Nonlocal Schloegl Models for Autocatalysis: Loss and Reemergence of a Nonequilibrium Gibbs Phase Rule, Phys. Rev. Lett. 121, 120603 (2018). J. Marro and R. Dickman. ”Nonequilibrium phase transitions in lattice models” Cambridge University Press (2005). D.-J. Liu, X. Guo and J. W. Evans, Quadratic contact process: Phase separation with interface-orientationdependent equistability. Physical Review Letters 98, 050601 (2007). D.-J. Liu and J. W. Evans, Symmetry-breaking and percolation transitions in a surface reaction model with superlattice ordering. Physical Review Letters 84, 955 (2000). X. Feng, H. W. J. Blöte, and B. Nienhuis, Lattice gas with nearest-and next-to-nearest-neighbor exclusion, Phys. Rev. E 83.6, 061153 (2011). R. J. Creswick and S-Y. Kim, Critical exponents of the four-state Potts model, J. Phys. A 30.24, 8785 (1997). M. Dijkstra, R. van Roij, and R. Evans, Phase diagram of highly asymmetric binary hard-sphere mixtures, Phys. Rev. E 59.5, 5744 (1999). Th. Voigtmann, Multiple glasses in asymmetric binary hard spheres, EPL 96, 36006 (2011). E. Zaccarelli, S. M. Liddle, and W. C. K. Poon, On polydispersity and the hard sphere glass transition, Soft Matter 11.2, 324-330 (2015). P.N. Pusey and W. Van Megen, Phase behavior of concentrated suspensions of nearly hard colloidal sphere, Nature 320.6060, 340 (1986). B. Doliwa and A. Heuer, Cooperativity and spatial correlations near the glass transition: Computer simulation results for hard spheres and disks, Phys. Rev. E 61, 6898 (2000). A. Donev, F. H. Stillinger, and S. Torquato, Do binary hard disks exhibit an ideal glass transition?, Phys. Rev. Lett. 96.22, 225502 (2006). L. Santen and W. Krauth, Absence of thermodynamic phase transition in a model glass former, Nature 405, 550 (2000).