# Lattice Glass Model in Three Spatial Dimensions

Yoshihiko Nishikawa, Koji Hukushima

► To cite this version:

Yoshihiko Nishikawa, Koji Hukushima. Lattice Glass Model in Three Spatial Dimensions. Physical Review Letters, 2020, 125 (6), pp.065501. 10.1103/PhysRevLett.125.065501. hal-03161197

HAL Id: hal-03161197

https://hal.science/hal-03161197

Submitted on 2 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Lattice glass model in three spatial dimensions

Yoshihiko Nishikawa¹,²,* and Koji Hukushima²,³,†

¹Laboratoire Charles Coulomb, UMR 5221 CNRS,
Université de Montpellier, 34095 Montpellier, France

²Department of Basic Science, Graduate School of Arts and Sciences,
The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan

³Komaba Institute for Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan

(Dated: July 14, 2020)

The understanding of thermodynamic glass transition has been hindered by the lack of proper models beyond mean-field theories. Here, we propose a three-dimensional lattice glass model on a simple cubic lattice that exhibits the typical dynamics observed in fragile supercooled liquids such as two-step relaxation, super-Arrhenius growth in the relaxation time, and dynamical heterogeneity. Using advanced Monte Carlo methods, we compute the thermodynamic properties deep inside the glassy temperature regime, well below the onset temperature of the slow dynamics. The specific heat has a finite jump towards the thermodynamic limit with critical exponents close to those expected from the hyperscaling and the random first-order transition theory for the glass transition. We also study an effective free energy of glasses, the Franz-Parisi potential, as a function of the overlap between equilibrium and quenched configurations. The effective free energy indicates the existence of a first-order phase transition, consistent with the random first-order transition theory. These findings strongly suggest that the glassy dynamics of the model has its origin in thermodynamics.

A thermodynamic (or ideal) glass transition in finite low dimensions has been actively discussed both theoretically and experimentally for decades since the seminal works by Kauzmann [1], and Adam and Gibbs [2]. A mean-field theory, the random first-order transition (RFOT) theory of structural glasses, proposed and developed in Refs. [3–6] indeed shows that a thermodynamic glass transition at finite temperature exists with vanishing configuration entropy, or complexity [7]. In finite dimensions, numerical simulation of models of fragile supercooled liquids is a promising way to theoretically explore a possibility of the thermodynamic glass transition. However, the notoriously long relaxation time prevents us to access directly low-temperature thermodynamics of glass-forming supercooled liquids. Although recent progress on particle models simulated using the swap Monte Carlo method [8, 9] allows us to get much more stable glass configurations at lower temperature or higher density, the thermodynamic glass transition is still inaccessible.

In an analogy to phase transitions into long-range ferromagnetic and crystal states, the lower critical dimension for the thermodynamic glass transition in models with discrete symmetry may be lower than that in models with continuous symmetry. Thus, it would be crucial to explore a possibility of a thermodynamic glass transition in finite-dimensional lattice models with discrete symmetry. Although several simple lattice models with the mean-field thermodynamic glass transition have already been proposed [10–12], they are not fully suitable to study the finite-dimensional glass transition: For the models in Refs. [10, 12], their mean-field glass transitions are turned into a crossover in finite dimensions or their low-temperature glassy states are unstable due to crystallization. The other model [11], while having the

typical glassy dynamics, is a monodisperse model, where an efficient algorithm such as the swap method [8, 9] is not known. Another lattice glass model was proposed in Ref. [13], which is shown to have irregular configurations as an ordered state. However, its autocorrelation function decays without any plateau even at low temperature. Thus the model does not have the essential features of the structural glasses.

To find finite-dimensional models that allow us to access equilibrium low-temperature states without crystallization is still actively discussed. In this letter, we propose a simple lattice glass model to study glassy behaviors in three dimensions. Our model, in contrast to several lattice models [10, 12, 13], shows typical two-step relaxation dynamics as observed in fragile supercooled liquids. We study a binary mixture of the model, where the non-local swap dynamics explained below provide the benefits for equilibration. By large-scale Monte Carlo simulations, we equilibrate the system at temperature well below that where two-step relaxation emerges. We study the effect of a coupling field conjugate to the overlap between the system and a quenched configuration using the Wang–Landau algorithm [14, 15], and compute the quenched version of the Franz–Parisi potential [16–18], an effective free energy of the glass transition. Our results show that the system indeed has a thermodynamic behavior consistent with the RFOT theory.

A lattice glass model we study in this letter is a binary mixture of particles defined by the Hamiltonian with a positive constant A

$$H(n, \sigma) = A \sum_{i=1}^{L^3} n_i \left( \sum_{j \in \partial i} n_j - \ell_{\sigma_i} \right)^2, \quad (1)$$

where the first summation runs over all the lattice sites.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Phase diagram of the model in the plane of temperature T and type-2 particle density $\rho_2$ with the total density $\rho = 0.75$. Crystal$_3$ and Crystal$_5$ stand crystalline phases where a crystal structure is mainly composed of particles with $\ell_1 = 3$ and $\ell_2 = 5$, respectively. The gas, liquid, crystalline phases are separated by first-order phase transitions. In the glassy region, the system shows an aging effect and its autocorrelation function relaxes in two steps. We study the system with $\rho_2/\rho = 0.6$ in this work.</figcaption></figure>

The lattice is the simple cubic lattice with linear dimension L. We denote the occupancy of site i with $n_i \in \{0,1\}$, and the type of a particle at site i with $\sigma_i \in \{1,2\}$. The parameter $\ell_{\sigma_i}$ determines the most favorable number of neighboring particles of type $\sigma_i$ particle. Here, we consider a binary mixture of $\ell_1 = 3$ and $\ell_2 = 5$ particles. The boundary conditions in all the directions are periodic. We study this model at finite temperature T by Monte Carlo simulation that preserves the number of each type of particles [19]: Two randomly chosen particles are swapped or a randomly chosen particle is moved to a vacant site of the lattice chosen also randomly with the Metropolis probability. Depending on the density of each type of particle, $\rho_1$ and $\rho_2$, with fixed total density $\rho = \rho_1 + \rho_2 = 0.75$, our model has crystal phases at low temperature, see Fig. 1. When $\rho_1 = 0.3$ and $\rho_2 = 0.45$, we confirmed the absence of a drop in the energy and a large peak in the specific heat in the glassy region, indicating no crystallization.

Our model is somewhat related to a softened model of a lattice glass proposed by Biroli and Mézard (BM) [10, 20]. The Hamiltonian of their model $H_{BM} = A \sum_{i=1} n_i \theta(\sum_{j \in \partial i} n_j - \ell_{\sigma_i}) (\sum_{j \in \partial i} n_j - \ell_{\sigma_i})$ with the Heaviside step function $\theta(x)$. The crucial difference between the two models is the number of neighboring particles that each particle energetically favors: In the soft BM model, any number of neighboring particles lower than $\ell_{\sigma}$ achieves the lowest energy while only $\ell_{\sigma}$ does in our model. This slight change makes the entropy of our model at low temperature smaller than that of the BM model and the dynamics more similar to fragile supercooled liquids.

The dynamics of our model is studied by a local Monte Carlo algorithm [19], which usually gives dynamics qualitatively similar to molecular and Brownian dynamics. In our simulation, particles can move to only neighboring sites. While any physical quantity of our model relaxes rapidly at high temperature, the relaxation gets extremely slow with decreasing temperature. To quantify the slow dynamics of the model, we study the autocorrelation function

$$C(t; t_w) = \left\langle \frac{\frac{1}{N} \sum_i \delta_{\sigma_i(t_w), \sigma_i(t_w+t)} - C_0}{1 - C_0} \right\rangle, \quad (2)$$

where $t_w$ is the waiting time, $C_0 = \rho(\rho_1^2 + \rho_2^2)$, $N = \rho L^3$ the total number of particles, and the summation is taken over sites with $n_i(t_w) = 1$. We also measure the dynamical susceptibility $\chi_4(t; t_w)$ characterizing the dynamical heterogeneity observed in supercooled liquids. In the limit $t_w \to \infty$, the system is in equilibrium and we denote $C(t) = C(t; t_w \to \infty)$ and $\chi_4(t) = \chi_4(t; t_w \to \infty)$. We set $t_w$ to a sufficiently large value to study the equilibrium dynamics of the model so that $C(t; t_w)$ and $C(t; t_w/10)$ agree with each other. Typical values of $t_w$ range from $10^2$ to $10^8$ Monte Carlo sweeps depending on temperature. Whereas the autocorrelation function of the system decays rapidly at high temperature, a two-step relaxation emerges in $C(t)$ at temperature lower than $T/A \simeq 0.6$, see Fig. 2 (see also [21] for a physical interpretation of the plateau in $C(t)$). The dynamical susceptibility $\chi_4(t)$ shows a peak at finite time, indicating the emergence of heterogeneous dynamics of the system [7, 22–25]. The peak value grows with decreasing temperature (see inset of Fig. 2), and it suggests that the dynamics gets more heterogeneous at lower temperature.

The RFOT theory of the glass transition predicts slow dynamics frozen at a dynamical transition temperature without any thermodynamic anomaly. In the vicinity of the dynamical transition temperature, the relaxation time diverges algebraically. However, an exponentially growing relaxation time well described by the Vogel–Fulcher–Tammann (VFT) law has been observed in experiments of fragile super-cooled liquids. Activated process in finite dimensions is supposed to wipe out the mean-field dynamical transition. Here, we study temperature dependence of the relaxation time $\tau_{\alpha}$ measured as a time when $C(t)$ decays to $e^{-1}$. Around temperature where the two-step relaxation emerges, the relaxation time shows a super-Arrhenius growth with decreasing temperature (see Fig. 3). We find that the VFT law fits our data at lower temperatures very well with $T_{VFT}/A = 0.177(6)$ (see also [21] for estimation of the dynamical transition temperature). Note that, as the relaxation time at low temperature increases with the system size, $T_{VFT}$ would shift to higher temperature in larger systems [21].

In the RFOT theory, the overlap between independent replicas has been studied as an order parameter for the
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) The autocorrelation function of the system. Inset shows an enlarged view around the plateau region. (b) The dynamical susceptibility $\chi_4(t)$ of the system. Inset shows temperature dependence of the peak and the plateau values of $\chi_4(t)$ that are denoted as $\chi_4^*$ and $\chi_4^p$, respectively. The system size $N = 6000$ ($L = 20$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. The Arrhenius plot of the relaxation time $\tau_\alpha$ of the system with $N = 6000$ ($L = 20$). The broken and dotted lines are the Arrhenius law with $\Delta E = 4$ and the VFT law with $T_{VFT}/A = 0.177$ respectively for guides to eyes.</figcaption></figure>

thermodynamic glass transition. The overlap successfully detects the glass transition in systems without spatial symmetry [26–29]. The overlap can detect even the slow dynamics at low temperature through its effective free energy [16–18]. In our model, however, the overlap and its distribution function can have no anomaly at any temperature due to its spatial symmetry [30]. The temperature dependence of the long-time limiting value of the dynamical susceptibility $\chi_4(t \to \infty)$, equivalent to the spin glass susceptibility $\chi_{SG}$, is indeed almost independent of temperature (see inset of Fig. 2). This is in contrast to a Potts spin glass model [31] where $\chi_{SG}$ in-

creases with approaching its transition temperature [32].

To study rare fluctuations and the effective free energy of the overlap appropriately, we introduce a field $\epsilon$ in our model: A system at temperature $T$ is coupled by the field to a quenched configuration randomly sampled at the same temperature in equilibrium [16–18]. The field explicitly breaks the invariance associated with the spatial symmetry, and thus the overlap can have nontrivial distributions. In the RFOT theory, the field induces a first-order transition at low temperature and the transition terminates at a critical point belonging to the random-field Ising model universality class [33]. The Franz–Parisìe potential reveals the existence of metastable states at low-temperature and allows us to compute the configurational entropy as a free-energy difference between two minima [34, 35]. Here, we use the Wang–Landau algorithm [14, 15] with the multi-overlap ensemble [36] to compute the density of states $\Omega_T(Q)$ as a function of the overlap $Q$ at a given temperature $T$ with a fixed reference configuration. The Franz–Parisìe potential $V_q(Q)$ is computed from $\Omega_T(Q)$ by averaging over quenched reference configurations. We prepared reference configurations by equilibrating the system at each temperature for $10^8$ Monte Carlo sweeps with non-local swaps of particles. The non-local swap dynamics is faster with a factor of $\sim 10^2$ at low temperature than the local swap dynamics while the factor slightly depends on temperature and the system size. The number of reference configurations is 912 for $T/A = 0.31$ and 0.33, and 48 for higher temperatures. At high temperature, the potential $V_q(Q)$ is convex and the overlap $[\langle Q \rangle]$ as a function of the field $\epsilon$ grows gradually, see Fig. 4 and Fig. 5 (c). Here, the brackets $\langle \cdot \rangle$ and $[\cdot]$ represent the thermal and the reference-configuration averages, respectively. At $T/A = 0.33$ and 0.31, the Franz–Parisìe potential is slightly non-convex (see Fig. 4), and the probability distribution of the overlap $[P_q(Q)]$ with finite $\epsilon$ has two

<figure><img src="image_5.png" /><figcaption>FIG. 4. The Franz–Parisìe potential $V_q(Q)$ of the system with $N = 1296$ ($L = 12$). At high temperature the potential is convex while it is slightly non-convex at low temperature, $T/A = 0.31$. Inset shows an enlarged view. The broken line is a guide to eyes to show the non-convexity.</figcaption></figure>
separated peaks whereas that at temperature higher than $T/A = 0.35$ shows a clear single peak at any $\epsilon$, see Fig. 5 (a) and (b), respectively [37]. We thus conclude that the coupling field induces a first-order transition into our model that terminates at a critical point at finite temperature as in mean-field and particle models in finite dimensions [17, 18, 33, 38, 39] (see Fig. 5 (d) for $T-\epsilon$ phase diagram).

When we assume the existence of the thermodynamic glass transition at $T_K > 0$ for a system with spatial symmetry, we find the system-size dependence of the effective transition point $\epsilon_c(L)$ and the configurational entropy density $s_{conf}(L)$ as follows. In the thermodynamic limit, an infinitesimal coupling field should make them have a finite overlap at $T < T_K$ whereas the overlap is strictly zero when $\epsilon = 0$ due to spatial symmetry [30]. In finite systems, the first-order transition at finite field thus never goes to the expected $T_K$, even in the limit $\epsilon \to +0$, while it may approach zero temperature in the limit. At temperature lower than $T_K$, the effective transition point $\epsilon_c$ should scale as $\epsilon_c = O(L^{-a})$ while at higher temperature $\epsilon_c - \epsilon_c^{(\infty)} = O(L^{-a})$ where $\epsilon_c^{(\infty)} = \lim_{L \to \infty} \epsilon_c(L) > 0$ and $a > 0$, which depends on systems [40, 41]. Regarding that the configurational entropy density $s_{conf}$ measured as the free energy difference in the Franz-Parisi potential is almost equivalent to $\epsilon_c$, we expect $s_{conf}$ of finite sys-

<figure><img src="image_6.png" /><figcaption>FIG. 5. (a, b) The probability distribution of the overlap averaged over quenched reference configurations [$P_q(Q)$] with finite coupling field $\epsilon$ at $T/A = 0.33$ (a) and 0.45 (b), respectively. (c) The overlap [$\langle Q \rangle$] as a function of the field $\epsilon$. With increasing $\epsilon$, [$\langle Q \rangle$] grows smoothly at high temperature while it shows a rapid increase at low temperature. The system size $N = 1296$ ($L = 12$). (d) Temperature-coupling phase diagram of the system. The empty points and their errors are estimated by the peak of the spin glass susceptibility $\chi_{SG}$ from the data for $L = 12$. The filled circle point represents an expected location of the critical point that is not determined precisely here. The triangular point $T_{VFT}$ is from the VFT-law fit of the relaxation time for $L = 20$, see Fig. 3. The diamond point indicates the dynamical transition temperature $T_c \simeq 0.348(4)$ estimated by fitting [21].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Temperature dependence of (a) the specific heat $c$ and (b) its temperature derivative $\partial c/\partial T$ of the system. Inset of (b) shows the peak value of $\partial c/\partial T$ as a function of $L$.</figcaption></figure>

tems to be finite even below expected $T_K$, but decreases with $\sim L^{-a}$, as $\epsilon_c(L)$ does. Studying the finite-size dependence of the first-order transition line and the configurational entropy would be a decisive test for the RFOT theory.

To measure the thermodynamic properties at low temperature, we use the non-local swap dynamics of randomly chosen pairs of particles and the exchange Monte Carlo (or parallel tempering) method [42] to further enhance equilibration. We also utilize the multiple-temperature reweighting technique [43–45], which significantly improves the accuracy of Monte Carlo results. Typical number of Monte Carlo sweeps for equilibration range from $10^8$ to $10^{11}$ per site depending on the system size.

At temperature $T/A \simeq 0.6$ where the two-step relaxation emerges, the specific heat is rather smooth, and it starts to drop without a divergent behavior at lower temperature (see Fig. 6 (a)). To study further the system-size dependence of the specific heat, we compute the temperature-derivative of the specific heat $\partial c/\partial T$. We find that $\partial c/\partial T$ grows with increasing the system size, indicating that the drop gets steeper with the system size. Similar size-dependence has been observed in a mean-field model with RFOT [46] and a three-dimensional Potts glass model that has a RFOT-like spin-glass transition [31].

Assuming the hyperscaling relation of the critical exponents $d\nu = 2 - \alpha$ with $d$ the spatial dimension, we find that the peak value of $\partial c/\partial T$ follows the finite-size scaling relation

$$ \left( \frac{\partial c}{\partial T} \right)^* \propto L^{\theta} \qquad (3) $$

with $\theta = (\alpha + 1)/\nu$. We thus can evaluate the two conventional critical exponents from $\theta$ as $\nu = 3/(\theta + d)$ and $\alpha = 2 - d\nu = 2 - 3d/(\theta + d)$. If the specific heat does not diverge but has a finite jump at a thermodynamic glass transition as in the RFOT theory, the critical exponent $\alpha = 0$, suggesting $\nu = 2/d$. The peak value $(\partial c/\partial T)^*$ thus grows algebraically with exponent $\theta = 1/\nu = d/2$.
This provides a useful way for studying thermodynamic
anomaly with a finite-size scaling analysis even when an
order parameter is unknown. Here, in our model, the
exponent $\theta = 1.44(2)$ (see inset of Fig. 6 (b)), implying
$\nu = 0.68(1)$ and $\alpha = -0.03(1)$, marginally compatible
with the RFOT theory.

In summary, we proposed a lattice glass model that
is very stable against crystallization and shows the typ-
ical dynamics observed in fragile supercooled liquids at
low temperature. Our numerical results show that a glass
transition is detected by a singularity of the temperature-
derivative of the specific heat with exponents compatible
with the RFOT theory. We also numerically computed
the quenched version of the Franz-Parisi potential using
the Wang-Landau algorithm. The potential shows very
large overlap fluctuations at low temperature, and a first-
order transition was found in the coupled system that ter-
minates at a critical point. The large overlap fluctuations
strongly suggest that low-temperature glassy dynamics
emerge from thermodynamics of the system, similar to
the mean-field models with the RFOT. We thus conclude
that our model is useful to study and further examine
the mean-field RFOT predictions. Thanks to the lattice
nature of our model, it is easy to study its mean-field
solution using the cavity method and compare it with
our results. To determine precisely the phase diagram
both in mean-field and finite dimensions will constitute
a crucial test of the validity of the theory.

The authors thank A. Ikeda for many useful discus-
sions. Y.N. is grateful to J. Takahashi, T. Takahashi,
S. Takabe, H. Yoshino, M. Ozawa, L. Berthier and
D. Coslovich for critical comments and useful discussions.
This work was supported by JSPS KAKENHI Grant
Numbers 17J10496, 17H02923 and 19H04125. Numeri-
cal simulation in this work has mainly been performed
by using the facility of the Supercomputer Center, Insti-
tute for Solid State Physics, the University of Tokyo.

<footnote>* yoshihiko.nishikawa@umontpellier.fr</footnote>

<footnote>† k-hukushima@g.ecc.u-tokyo.ac.jp</footnote>

<footnote>[1] W. Kauzmann, Chem. Rev. 43, 219 (1948).</footnote>

<footnote>[2] G. Adam and J. H. Gibbs, The Journal of Chemical
Physics 43, 139 (1965).</footnote>

<footnote>[3] T. R. Kirkpatrick and D. Thirumalai, Physical Review B 36, 5388 (1987).</footnote>

<footnote>[4] T. R. Kirkpatrick and P. G. Wolynes, Physical Review B 36, 8552 (1987).</footnote>

<footnote>[5] T. R. Kirkpatrick and D. Thirumalai, Physical Review A 37, 4439 (1988).</footnote>

<footnote>[6] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Physical Review A 40, 1045 (1989).</footnote>

<footnote>[7] L. Berthier and G. Biroli, Reviews of Modern Physics 83, 587 (2011).</footnote>

<footnote>[8] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Physical Review Letters 116, 238002 (2016).</footnote>

<footnote>[9] A. Ninarello, L. Berthier, and D. Coslovich, Physical Review X 7, 021039 (2017).</footnote>

<footnote>[10] G. Biroli and M. Mézard, Physical Review Letters 88, 025501 (2001).</footnote>

<footnote>[11] M. P. Ciamarra, M. Tarzia, A. de Candia, and A. Coniglio, Physical Review E 67, 057105 (2003).</footnote>

<footnote>[12] G. D. McCullagh, D. Cellai, A. Lawlor, and K. A. Dawson, Physical Review E 71, 030102(R) (2005).</footnote>

<footnote>[13] S.-i. Sasa, Physical Review Letters 109, 165702 (2012).</footnote>

<footnote>[14] F. Wang and D. P. Landau, Physical Review Letters 86, 2050 (2001).</footnote>

<footnote>[15] F. Wang and D. P. Landau, Physical Review E 64, 056101 (2001).</footnote>

<footnote>[16] S. Franz and G. Parisi, Journal de Physique I, 21 (1995).</footnote>

<footnote>[17] S. Franz and G. Parisi, Physical Review Letters 79, 2486 (1997).</footnote>

<footnote>[18] S. Franz and G. Parisi, Physica A: Statistical Mechanics and its Applications 261, 317 (1998).</footnote>

<footnote>[19] K. Kawasaki, Phys. Rev. 145, 224 (1966).</footnote>

<footnote>[20] L. Foini, G. Semerjian, and F. Zamponi, Physical Review B 83, 094513 (2011).</footnote>

<footnote>[21] “See Supplemental Material for a further discussion on the dynamics, the finite-size dependence of the relaxation time, and the determination of the dynamical transition temperature $T_d$, which includes Refs. [7, 10, 17, 18, 33, 47–49],”.</footnote>

<footnote>[22] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Physical Review Letters 79, 2827 (1997).</footnote>

<footnote>[23] R. Yamamoto and A. Onuki, Physical Review Letters 81, 4915 (1998).</footnote>

<footnote>[24] M. D. Ediger, Annual Review of Physical Chemistry 51, 99 (2000).</footnote>

<footnote>[25] C. Toninelli, M. Wyart, L. Berthier, G. Biroli, and J.-P. Bouchaud, Physical Review E 71, 041505 (2005).</footnote>

<footnote>[26] C. Cammarota and G. Biroli, Proceedings of the National Academy of Sciences 109, 8850 (2012).</footnote>

<footnote>[27] C. Cammarota and G. Biroli, The Journal of chemical physics 138, 12A547 (2013).</footnote>

<footnote>[28] W. Kob and L. Berthier, Physical Review Letters 110, 245702 (2013).</footnote>

<footnote>[29] M. Ozawa, W. Kob, A. Ikeda, and K. Miyazaki, Proceedings of the National Academy of Sciences 112, 6914 (2015).</footnote>

<footnote>[30] M. Mézard and G. Parisi, “Glasses and replicas,” in *Structural Glasses and Supercooled Liquids* (John Wiley & Sons, Ltd, 2012) Chap. 4, pp. 151–191.</footnote>

<footnote>[31] T. Takahashi and K. Hukushima, Physical Review E 91, 020102(R) (2015).</footnote>

<footnote>[32] T. Takahashi and K. Hukushima, (unpublished).</footnote>

<footnote>[33] S. Franz and G. Parisi, Journal of Statistical Mechanics: Theory and Experiment 2013, P11012 (2013).</footnote>

<footnote>[34] L. Berthier and D. Coslovich, Proceedings of the National Academy of Sciences 111, 11668 (2014).</footnote>

<footnote>[35] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, Proceedings of the National Academy of Sciences 114, 11356 (2017).</footnote>

<footnote>[36] B. A. Berg and W. Janke, Physical Review Letters 80, 4771 (1998).</footnote>

<footnote>[37] Note that the Franz-Parisi potential is convex at any temperature in the thermodynamic limit as the free-energy barrier is sub-extensive, and the non-convexity is seen only in finite systems.</footnote>

<footnote>[38] L. Berthier, Physical Review E 88, 022313 (2013).</footnote>

<footnote>[39] L. Berthier and R. L. Jack, Physical Review Letters 114,</footnote>
205701 (2015).

[40] M. E. Fisher and A. N. Berker, Physical Review B **26**, 2507 (1982).

[41] M. Mueller, W. Janke, and D. A. Johnston, Physical Review Letters **112**, 200601 (2014).

[42] K. Hukushima and K. Nemoto, Journal of the Physical Society of Japan **65**, 1604 (1996).

[43] A. M. Ferrenberg and R. H. Swendsen, Physical Review Letters **61**, 2635 (1988).

[44] A. M. Ferrenberg and R. H. Swendsen, Physical Review Letters **63**, 1195 (1989).

[45] E. P. Münger and M. A. Novotny, Physical Review B **43**, 5773 (1991).

[46] A. Billoire, L. Giomi, and E. Marinari, Europhysics Letters (EPL) **71**, 824 (2005).

[47] L. Berthier, G. Biroli, D. Coslovich, W. Kob, and C. Toninelli, Physical Review E **86**, 031502 (2012).

[48] D. Coslovich, A. Ninarello, and L. Berthier, SciPost Physics **7**, 077 (2019), 1811.03171.

[49] S. Karmakar, C. Dasgupta, and S. Sastry, Proceedings of the National Academy of Sciences of the United States of America **106**, 3675 (2009).
## Supplemental Item 1: Dynamics at low temperature

In off-lattice particle models for supercooled liquids and glasses, the dynamics at low temperature shows two-step relaxation as well as our model. The physical interpretation of the plateau regime in the particle models is well known that particles vibrate locally inside cages effectively formed by surrounding particles [7]. However, since particles of our model defined on a lattice have discrete degrees of freedom, the appropriate physical interpretation analogous to the local vibration in the plateau regime is unclear. In this supplemental item, we show evidence of local vibrations of particles during the plateau regime.

Thanks to the discrete nature of our model, it is easy to find lattice sites that have no contribution to the autocorrelation function Eq. (2), where $\delta_{\sigma_i(t_w), \sigma_i(t_w+t)} = 0$, i.e. the lattice sites that have different types of particles at $t=0$ and $t>0$ or that are vacant at $t>0$ but occupied at $t=0$. Roughly speaking, a cluster represents a movable region during time interval $t$. We identify connected clusters of those lattice sites at time $t$, and compute the normalized distribution of the cluster size $n_c$. Fig. S1 shows the distribution of the connected-cluster size depending on the time $t$.

At temperature $T/A = 0.35$, the autocorrelation function has a plateau from time $t \approx 10$ to $t \approx 10^3$, see Fig. 2. In this time regime, the size of connected clusters $n_c \lesssim 10$, and the distribution is almost independent of time, see Fig. S1. Although the size distribution does not depend on time very much, the positions of the clusters change with time during the plateau, and some of them disappear with time (see Fig. S2 for real-space visualization of the changed lattice sites). We thus conclude, at short time scale $t \lesssim 10^3$, the system has local vibrations that involve only $O(1)$ number of particles. These local vibrations produce the plateau in the autocorrelation function, even in the lattice model.

As time passes, the distribution of $n_c$ has a broad tail towards larger $n_c$, and eventually becomes double-peaked at large time $t \approx \tau_\alpha \approx 10^7$. The emergence of two peaks in the distribution with small clusters of $n_c = O(1)$ and giant clusters of $n_c = O(N)$ indicates the heterogeneous dynamics. At time $t = 10^7$, the largest cluster spreads over the system while there are still some regions where only small clusters occupy, see Fig. S2.

<figure><img src="image_8.png" /><figcaption>FIG. S1. The normalized distribution function of the size of connected clusters at $T/A = 0.35$ for (a) $10^1 \le t \le 5 \times 10^5$ and (b) $3 \times 10^6 \le t \le 8 \times 10^7$. The system size $L = 20$ ($N = 6000$).</figcaption></figure>

## Supplemental Item 2: Estimation of the dynamical transition temperature

In the RFOT theory, the relaxation time diverges algebraically at temperature $T_d > T_K$ rather than exponentially as observed in experiments. Whereas the dynamical transition turns into a crossover in finite dimensions, a well-defined localization transition of the potential energy landscape was shown to control the dynamical cross-over [48]. However, their method to identify the localization transition is not available in our model. We thus estimate the dynamical transition temperature $T_d$ by a more traditional method, that is, power-law fitting of the relaxation time as done in Ref. [10].

We show in Fig. S3 the best fit of the inverse relaxation time to a power-law behavior $\sim (T-T_d)^\gamma$, with $T_d = 0.348(4)$ and $\gamma = 4.18(12)$. In mean-field models, the dynamical transition temperature is close to the critical temperature of the $\varepsilon$-coupled system [17, 18, 33]. Our estimation of $T_d$ is compatible with the temperature-coupling phase diagram Fig. 6.
<figure><img src="image_9.png" /><figcaption>FIG. S2. Real-space visualization of typical configurations at time $t = 10^1, 10^2, 5 \times 10^2, 10^3$, and $10^7$, from left top to right bottom. The temperature $T/A = 0.35$. The lattice sites which have different particles in the configuration at finite time $t > 0$ from the original configuration at $t = 0$ are filled with cubes. Here, at this temperature, the autocorrelation function has a plateau when $10^1 \lesssim t \lesssim 10^3$. For the configuration at $t = 10^7$, the lattice sites belonging to the largest connected cluster are filled with blue cubes, and those for the other clusters are filled with red cubes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. S3. The inverse relaxation time $\tau_{\alpha}^{-1}$ as a function of $(T - T_d)/T_d$ and a power-law curve $\sim (T - T_d)^{\gamma}$. The dynamical transition temperature $T_c$ and the exponent $\gamma$ are estimated by fitting using the data at $0.4 \le T/A \le 0.6$, as $T_d = 0.348(4)$ and $\gamma = 4.18(12)$, respectively.</figcaption></figure>

## Supplemental Item 3: Finite-size effect in the relaxation time

The relaxation time of glass forming liquids has a decreasing behavior at moderately low temperature [47, 49]. The magnitude of the decrease is especially remarkable just below the dynamical transition temperature (or the mode-coupling temperature) [49]. In our model, with the local dynamics, the relaxation time does not decrease with the system size even around the dynamical transition temperature $T_d = 0.348(3)$, see Fig. S4. At lower temperature, the relaxation time rather increases with the system size. The finite-size effect is more prominent in the relaxation time of the non-local dynamics we used for equilibrium computation, see Fig. S4. Intuitively, the finite-size effect is clearer because of the following reason. The relaxation process at moderately low temperature is mixed with many modes, with each contribution. The slowest mode of the relaxation process, which diverges towards the transition
temperature and has strong finite size effects, is buried in many other modes due to its small contribution in the dynamics of inefficient algorithms, including simple local dynamics. The relaxation time shows a finite-size effect only when the slowest mode contribution becomes very large, which is very low temperature. Efficient algorithms, on the other hand, can reduce contributions of modes faster than the slowest mode but still very slow, and thus a finite-size effect that the slowest mode has is seen more clearly. This is actually quite general in other models, even in simple ferromagnetic spin models, by comparing the dynamics of the simple Metropolis algorithm and more efficient algorithms such as the over-relaxation, for instance.

<figure><img src="image_11.png" /><figcaption>FIG. S4. The relaxation time as a function of inverse temperature for $L = 5, 10$, and 20 for the local and non-local dynamics. The inverse dynamical transition temperature $A/T_d$ is estimated in Supplemental Item 2.</figcaption></figure>