Solvable Models of Supercooled Liquids in Three Dimensions

TL;DR

This paper presents a three-dimensional kinetically-constrained model for supercooled liquids, providing parameter-free predictions that match simulations well, especially in low-temperature regimes where traditional theories fail.

Contribution

It introduces a solvable 3D supercooled liquid model with quantitative predictions aligning with simulations, extending understanding beyond Mode-Coupling Theory.

Findings

Excellent agreement with numerical simulations in low-temperature regimes.

Model's potential applicability to real supercooled liquids.

Challenges traditional Mode-Coupling Theory predictions.

Abstract

We introduce a supercooled liquid model and obtain parameter-free quantitative predictions that are in excellent agreement with numerical simulations, notably in the hard low-temperature region characterized by strong deviations from Mode-Coupling-Theory behavior. The model is the Fredrickson-Andersen Kinetically-Constrained-Model on the three-dimensional $M$-layer lattice. The agreement has implications beyond the specific model considered because the theory is potentially valid for many more systems, including realistic models and actual supercooled liquids.

Equations

$$\varphi (t)-{\varphi }_{plat}\approx \frac{1}{(t/t_0{)}^a},t\gg 1,$$

$$g(x,t)\equiv \left(\frac{1}{M}\sum _{\alpha =1}^M{\varphi }_i^{\alpha }(t)\right)-{\varphi }_{plat},$$

$$⟨g(x,t)g(y,t)⟩-⟨g(x,t)⟩⟨g(y,t)⟩\stackrel{\text{MF}}{\approx }\frac{t^{a(2-d/2)}}{M}f\left(\frac{x-y}{\xi (t)}\right).$$

$$\frac{t_G^{a(2-d/2)}}{M}\approx \frac{1}{t_G^{2a}}\to t_G\approx M^{\frac{1}{a(4-d/2)}}.$$

$$⟨g(x_1,t_1)\dots g(x_K,t_K)⟩\approx [\hat{g}(x_1,t_1)\dots \hat{g}(x_K,t_K)].$$

$$\sigma +s(x)=-\alpha {\nabla }^2\hat{g}(x,t)-\lambda {\hat{g}}^2(x,t)+\frac{d}{dt}{\int }_0^t\hat{g}(x,t-s)\hat{g}(x,s)ds.$$

$$[s(x)]=0,[s(x)s(y)]=\Delta {\sigma }^2\delta (x-y).$$

$$\Delta {\sigma }^2=\frac{0.285}{M},\alpha =0.411,\sigma =0.222\times (T_c-T).$$

$$g(t)\equiv \frac{1}{N_{tot}}\sum _{i,\alpha }{g}_i^{\alpha }(t),{\chi }_4(t)\equiv N_{tot}(⟨g^2(t)⟩-⟨g(t){⟩}^2).$$

$$g(t)=\frac{1}{(t/t_0{)}^a}+\frac{{\delta }_1}{(t/t_0{)}^{2a}}+\dots$$

$$\frac{\Delta \ln g}{\Delta \ln t}=-a-a{\delta }_{1}g+\dots$$

$$a\approx 0.352,t_0\approx 2.30.$$

$$k=c-m+1$$

$$p=\frac{e^{\beta }}{1+e^{\beta }}.$$

$$P^{(i)}(\{s_j:j\in \partial i\})=\prod _{j\in \partial i}{P}_B(s_j),$$

$$P_B(s)\equiv {\delta }_{s,-1}P+{\delta }_{s,1}(1-P).$$

$$P=p{P}_{k-1,c-1},P_{site}=p{P}_{k,c}.$$

$$c_{L}(p)\approx{p_{1}^{2}\,\Delta\over\mu}\,L\,\mu^{L}+O(\lambda^{L})\$$

$$\mu =p\left(\genfrac{}{}{0ex}{}{c-2}{k-2}\right)P^{k-2}(1-P{)}^{c-k}$$

$$\Delta\equiv{1-\mu\over\mu}P(1-P)\$$

$$p_{1}\equiv p{c-1\choose k-1}P^{k-1}(1-P)^{c-k}\$$

$$P=\frac{3}{4},\mu =\frac{1}{3},P_{site}=\frac{21}{32}\mathrm{for}p=p_c=\frac{8}{9}\to T_c=\frac{1}{\ln 8}\approx .480898$$

$$⟨g(x,t)g(y,t)⟩-⟨g(x,t)⟩⟨g(y,t)⟩=\frac{1}{M}\sum _{L=0}^{\infty }{N}_L(x,y)c_L^{Bethe}(t)$$

$$c_L^{Bethe}(t)\equiv ⟨g_i(t)g_j(t)⟩-⟨g_i(t)⟩⟨g_j(t)⟩$$

$$N_L(x,y)\approx \frac{c(c-1{)}^{L-1}}{\rho }G(x-y).$$

$$N_L(k,k^{\prime })\approx \delta (k+k^{\prime })\frac{c}{\rho (c-1)}(c-1{)}^L\exp [-L{D}_{NBW}{k}^2]$$

$$D_{NBW}\equiv \underset{t\to \infty }{\lim }\frac{⟨||{\mathbf{x}}_t|{|}^2⟩}{2dt}.$$

$$\underset{t\to \infty }{\lim }⟨g(x,t)⟩=P_{site}(p)-P_{site}(p_c)\approx \frac{27}{16\sqrt{2}}\delta p^{1/2}$$

$$\underset{t\to \infty }{\lim }{c}_L^{Bethe}(t)=c_L^{Bethe}(p)\approx c_{1}L{\mu }^L(p)$$

$$\mu (p)\approx \frac{1}{c-1}(1-c_2\delta p^{1/2})$$

Full text

The ubiquity of glass in nature and technology has driven research in this area for decades but there is still no agreement on the basic mechanism by which a supercooled liquid forms a glass Biroli and Garrahan (2013); Wolynes and Lubchenko (2012); Götze (2008). An essential difficulty for deciding over competing theories is the fact that often they make only qualitative statements. For instance, the classic thermodynamics vs. dynamics controversy resolves around the putative divergences of susceptibilities and correlation lengths but, while the debate has produced many conceptual developments Biroli and Garrahan (2013); Wolynes and Lubchenko (2012), there are no precise quantitative predictions to be be matched with experiments and numerical simulations. In other words, models of supercooled liquids in three dimensions, either realistic or on lattice, have been studied previously only by means of numerical simulations. In the following we present a model in three dimensions in which the gap between theory and experiments can instead be filled. Furthermore the agreement has implications beyond the specific model considered because the theory is potentially valid for many more physical systems, including the most realistic models and actual supercooled liquids.

Mode-Coupling-Theory (MCT) Götze (2008) is widely popular in the experimental literature because it captures many qualitative features of the physics of liquids upon supercooling, notably two-step relaxation and stretched exponential decay. Furthermore it agrees quantitatively with numerical simulations although one has to replace the values of some MCT parameters with values extracted from data fits Nauroth and Kob (1997); Kob (1999); Sciortino and Kob (2001); Weysser et al. (2010). The essential problem is that it predicts a dynamical arrest transition at a temperature where actual systems are still in the liquid phase. In spite of this serious drawback many believe that the MCT transition is still relevant and determines a crossover from power-law to exponential increase of the relaxation time that is widely observed. Further support to this scenario comes from the fact that simple liquid models in the limit of infinite physical dimension $d$ display a sharp transition qualitatively similar to the one of MCT, although MCT itself is quantitatively wrong in this limit Charbonneau et al. (2017). In that case the sharp transition is clearly a mean-field artifact due to the $d\rightarrow\infty$ limit and it should become a crossover as soon as the dimension is finite. Similarly, the success of kinetically constrained models (KCM) in reproducing the physics of supercooled liquids is often attributed to the presence of an avoided MCT-like transition that becomes sharp when one switches from lattices in finite dimensions to the Bethe lattice where mean-field theory is correct Fredrickson and Andersen (1984, 1985); Sellitto et al. (2005); Sellitto (2015); De Candia et al. (2016); Franz and Sellitto (2013); Ikeda et al. (2017); Sausset et al. (2010).

Building on the analogy between MCT and spin-glass models with one step of Parisi’s Replica-Symmetry-Breaking (1RSB) discovered more that thirty years ago Kirkpatrick and Thirumalai (1987), it has been recently proposed that to fix MCT one has to replace it with a set of stochastic dynamical equations called Stochastic-Beta-Relaxation (SBR) Rizzo (2014, 2016). SBR describes the $\beta$ regime, i.e. the time-scale when dynamic correlations stay close to a plateau value, and it has a simple and intuitive interpretation: it is basically MCT supplemented with random spatial fluctuations of the temperature that are effectively quenched on the $\beta$ time-scale. SBR is promising as it seems to cure the drawbacks of MCT, displaying in particular the crossover and also dynamical heterogeneities, without spoiling its successes, in particular two-step relaxation and stretched exponential decay Rizzo and Voigtmann (2015a, b). SBR is a universal theory potentially valid for many different microscopic models, the specific model determining the values of its few (five) quantitative parameters. In the following we show that SBR holds for the models we study. To do so we have first computed the five SBR parameters, then solved numerically the dynamical stochastic equations and finally performed Monte-Carlo simulations to be compared with the theoretical predictions. It turned out that SBR provides an accurate quantitative parameter-free description of the dynamics, i.e. the models are solvable beyond mean-field theory.

The model considered is a particular realization of the classic Fredrickson-Andersen (FA) KCM Fredrickson and Andersen (1984, 1985). The FA model is made by Ising spins on the sites of a lattice that are independent, the Hamiltonian being $H=\sum_{i}s_{i}$, but obey a kinetically constrained dynamics: a spin can flip only if it has at least $m$ of its $c$ nearest neighbors in the excited (up) state. An equilibrium configuration is thus easily generated numerically at time zero and one typically measures the persistence. More precisely, we define the local persistence $\phi_{i}(t)$ as equal to one if $s_{i}(t^{\prime})=-1$ for all $0\leq t^{\prime}\leq t$ and zero otherwise, thus the averaged persistence is the number of negative sites that have never flipped at time $t$ divided by the total number of spins 111Note that, at variance with ours, the definition of persistence often used in the literature takes into account all sites that did not flip, both negative and positive. The quantitative agreement between numerical data and theory is expected to hold also using the conventional definition of the persistence or other proxies of the correlation, e.g. the overlap. This can be understood within the Renormalization Group framework to be discussed later in the text as all observables have the same behavior provided they have a non-zero projection on the critical mode. Quantitatively, different (shifted) observables are just related by a rescaling factor. The rationale for our choice is that analytic computations of the quantitative SBR parameters, see discussion later and details in the supplemental material, are simpler because our definition is as close as possible to bootstrap percolation.. The FA model on the Bethe lattice is known to exhibit a dynamical arrest transition of the MCT type Sellitto et al. (2005); Sellitto (2015); De Candia et al. (2016); Franz and Sellitto (2013); Ikeda et al. (2017); Sausset et al. (2010): at the critical temperature $T_{c}$ the persistence remains blocked to a plateau value $\phi_{plat}$ that is approached in a power-law fashion. The FA dynamical transition is intimately related to bootstrap percolation (BP) and both $T_{c}$ and $\phi_{plat}$ can be computed from its solution on the Bethe lattice as discussed in the supplemental material. In particular for connectivity $c=4$ and $m=2$ the average persistence $\phi(t)$ obeys at $T_{c}=0.480898$:

[TABLE]

where $\phi_{plat}=21/32$. At present, analytic expressions of $t_{0}$ and $a$ are not available but they can be estimated from numerical simulations as $a\approx 0.352$ and $t_{0}\approx 2.30$. It is well-known that the sharp transition is instead avoided when the FA model is studied on regular lattices in two and three dimensions but no first-principle theoretical description of the dynamics can be obtained in those cases. As we will show in the following, such a description can instead be obtained on the finite dimensional lattice we consider here.

We studied the FA model (with $m=2$) on the (random) lattice in three dimensions yielded by the application of $M$-layer construction of Altieri et al. (2017) to the diamond cubic lattice (that has connectivity $c=4$, see fig. 1). The $M$-layer construction can be applied to any lattice: to obtain a random instance one considers $M$ copies of the original lattice, say the square lattice in $d=2$ as in the figure, rewires through a random permutation the $M$ links corresponding to a given link on the original $(M=1)$ lattice and finally repeats the procedure for each link of the original lattice as shown in Fig. (1). It can be easily seen that short loops in the lattice are rare for large values of $M$ and the lattice is locally tree-like. Loops are nevertheless present at large distances, thus at any finite $M$ the lattice is finite dimensional although it looks like a Bethe lattice at short distances. Given that for each site $i=1,\dots,N$ of the original lattice ($M=1$) there are $M$ spins $s_{i}^{\alpha}$, $\alpha=1,\dots,M$ the total number of sites is $N_{tot}=M\times N$ and the natural local order parameter is the average over the layers of the local persistence minus the plateau value

[TABLE]

where $x$ is the spatial coordinate of site $i$. For $M$ large but finite it is natural to expect that any observable takes the same value it has on the Bethe lattice with small $O(1/M)$ corrections: the model should be solvable by the Bethe approximation. This is indeed true except at the Bethe lattice critical temperature: while the Bethe approximation predicts that the averaged order parameter $g(t)$ never reaches zero the three-dimensional nature of the $M$-layer lattice implies that this cannot be true and leads to a dramatic deviation from mean-field behavior.

This is indeed confirmed by numerical data. In fig. (2) we plot the time decay of the shifted persistence $\langle g(x,t)\rangle$ where here and in the following the angle brackets mean average with respect to: i) different instances of the random lattice generated by the $M$-layer construction, ii) different initial equilibrium configurations and iii) different thermal histories. Monte Carlo simulations were carried at the critical temperature of the Bethe lattice $T_{c}=0.480898$ for systems with $M=3000$ (with size of the original diamond cubic lattice $L=8$), $M=50000$, $M=100000$, $M=200000$ ($L=4$). At initial times the data follow the mean-field (MF) curve corresponding to $M=\infty$ (obtained from numerical simulations on the Bethe lattice (see the SI) and deviate from it at larger times (increasing with $M$) reaching the plateau value ($g=0$) in finite time.

The difficult problem is to obtain theoretical predictions in the region where the data deviate from the mean-field curve and to proceed it is useful to examine the role of fluctuations. Within the $M$-layer a mean-field approximation to fluctuations can be obtained as a sum over non-backtracking walks of Bethe lattice fluctuations Altieri et al. (2017). At the critical temperature this yields 222See details in the Supplemental Material, which includes Refs. Parisi and Rizzo (2013); Schwarz et al. (2006); Franz et al. (2011); Fitzner and van der Hofstad (2013):

[TABLE]

where $f(x)$ is a scaling function, the correlation length diverges with time as $\xi(t)\propto t^{a/2}$ ansd $d$ is the space dimension. We can now invoke the Ginzburg criterion and argue that the MF approximation $\langle g(x,t)\rangle\stackrel{{\scriptstyle\text{MF}}}{{\approx}}1/(t/t_{0})^{a}$ is accurate as long as fluctuations around the mean are small. Generically they are small due to the $1/M$ prefactor but we see that there is a time-scale $t_{G}$ when, due to the $t^{a(2-d/2)}$ prefactor they become comparable with the (squared) order parameter $\langle g(x,t)\rangle\stackrel{{\scriptstyle\text{MF}}}{{\approx}}1/(t/t_{0})^{a}$:

[TABLE]

Thus on this time-scale MF theory must fail and most notably the spurious transition will be avoided, for instance the dressed propagator on the LHS of (3) will deviate from the bare expression on the RHS and the actual correlation length will cease to grow. We note that the Ginzburg time grows with $M$ and thus deviations from MF occur at later time for increasing values of $M$ in agreement with the data of fig. (2). Most importantly since $t_{G}\approx M^{1\over a(4-d/2)}$ is large for large $M$, the order parameter is $O(1/t_{G}^{a})$ small and the correlation length is $O(t_{G}^{a/2})$ large: this grants that deviations from mean-field theory are described by an effective Landau theory, because one can retain only the lowest orders in the Taylor expansion of the order parameter and its space and time derivatives. Following the arguments and computations of Rizzo (2014, 2016) we argue that the effective theory is SBR, meaning that the generic $K$-point average obeys for $1\ll M<\infty$:

[TABLE]

where $\hat{g}(x,t)$ in the RHS is the solution of the SBR equations:

[TABLE]

The separation parameter $\sigma$ measures the distance from the critical point and vanishes at $T=T_{c}$. The square brackets mean average with respect to the field $s(x)$ that is a quenched random fluctuation of $\sigma$, Gaussian and delta-correlated in space:

[TABLE]

the SBR equations have to be solved with the small-time condition $\lim_{t\rightarrow 0}\hat{g}(x,t)(t/t_{0})^{a}=1$ where $\lambda$ and $a$ are related by the MCT relationship $\lambda={\Gamma^{2}(1-a)\over\Gamma(1-2a)}$. In practice for times smaller than $t_{G}$ the observables on the LHS of eq. (5) can be accurately approximated with the values they have on the Bethe lattice while on times of order $t_{G}$ they are described by the RHS. This explains the peculiar initial conditions of the SBR equations: the short-time behavior on times $O(t_{G})$ matches long-time behavior for times $1\ll t\ll t_{G}$, i.e. the mean-field result given by eq. 1.

Equation (5) embodies the power of the effective theory approach: on the LHS we have a model with a complex microscopic dynamics for which no analytic treatment of dynamics is available (not even on the Bethe lattice), on the RHS we have a (numerically) solvable set of equations that were derived in Rizzo (2014, 2016) starting from symmetry considerations (essentially the detailed balance property of the dynamics) but without reference to any specific microscopic model. The microscopic details determine the actual values of the five SBR parameters $a$, $t_{0}$, $\alpha$, $\Delta\sigma$ and $\sigma$ that are needed to get quantitative predictions. Using recent developments on bootstrap percolation on the Bethe lattice Rizzo (2019) and some lattice-dependent geometrical constants, we obtain (details in the SI):

[TABLE]

Within SBR, mean-field theory is recovered setting $\Delta\sigma^{2}=0$, in this case $\hat{g}(x,t)$ is constant in space, the gradient term plays no role and one recovers the critical MCT equation Götze (2008), in particular for $\sigma\geq 0$ ($T<T_{c}$) $\phi(t)$ never goes below the plateau value. The $M$-layer construction allows to have a finite but small $\Delta\sigma$ so that the MCT transition is avoided and $\phi(t)$ crosses the plateau at a finite time for all values of $\sigma$. The SBR predictions corresponding to the data shown in fig. (2) were obtained solving numerically (by space-time discretization) eq. (6) for many instances of the $s(x)$ in a box of size $L$. From the figure we note that the quality of the SBR predictions increases with $M$ and is excellent for $M=200000$, especially considering that there is no single fitting parameter.

SBR is a powerful theory that provides not only the average dynamical order parameter but, according to eq. (5), also all possible fluctuations. To demonstrate this, in the inset of fig. (2), we plot parametrically the $\chi_{4}(t)$ function that yields the fluctuations of the persistence density:

[TABLE]

According to the MF expression (3) $\chi_{4}(t)$ should diverge with time as $t^{2a}$, (leading to a MF asymptote $\chi_{4}(g)\propto g^{-2}$), instead on the Ginzburg time scale $t_{G}$ over which $g(t)$ deviates from MF and reaches zero $\chi_{4}(g)$ deviates from the MF law and remains finite. Note that the agreement between numerical data and the numerical solution of the SBR equations is even better in the parametric representation.

SBR can be applied in other dimensions as well, we have considered (SI) the $d=0$ case that corresponds to finite-size effects in mean-field models on fully-connected or sparse random graphs with 1RSB at the so-called dynamical temperature $T_{d}$ Mezard and Montanari (2009). The interplay between the parameters $M$ and $L$ can be also be clarified in terms of the SBR equations (SI).

To discuss the results in a broader context we note that the $M$-layer construction can be applied virtually to all supercooled liquids lattice models, including different KCM’s Ritort and Sollich (2003) and plaquette models Jack and Garrahan (2016); Biroli et al. (2016), leading to analogous solvable non-MF models described by SBR. Furthermore SBR can provide quantitative theoretical predictions for generic tunable models Franz and Toninelli (2004); Franz and Montanari (2007); Caltagirone et al. (2011); Mari et al. (2009); Mari and Kurchan (2011); Charbonneau et al. (2014); Berthier et al. (2012) that in earlier studies could only be studied by means of numerical simulations. On the other hand, the lattice for $M\gg 1$ is rather different from the original $M=1$ lattice one is ideally interested in: the latter has many short loops while the former has very few. Thus the condition $M\gg 1$ alters artificially the three-dimensional geometry at the microscopic scale and one may ask if this hampers the applicability of SBR to realistic models and actual supercooled liquids. To clarify this point we stress that $M\gg 1$ is a sufficient but not necessary condition. A necessary condition in a generic system is that the dynamical correlation length is large enough to justify the use of a coarse-grained description: numerical simulations do indeed report correlation lengths significantly larger than the microscopic scale in supercooled liquids Lačević et al. (2003); Berthier et al. (2007); Flenner and Szamel (2013); Karmakar et al. (2014) while unfortunately they cannot be measured in current experimental settings. The natural framework to discuss coarse-grained observables is Wilson’s renormalization-group (RG) theory where each system corresponds to a particular point in the space of RG Hamiltonians that display all possible powers of the order parameter and its spatial and time derivatives and thus depends on a infinite number of coupling constants. In practice these additional terms lead to higher powers of $\hat{g}(x,t)$ and higher orders spatial and time derivatives in eqs. (6). SBR assumes that these additional terms can be neglected and this can be motivated by the following RG argument. The absence of a sharp dynamical transition in finite dimension implies the absence of a stable fixed point (FP), as a consequence all Hamiltonians flow under RG towards the high-temperature FP, howhever, if the correlation length is large it will take many RG steps for it to decrease to one. Since standard dimensional analysis implies that the coupling constants of the additional terms decrease close to the Gaussian fixed point (they are irrelevant operators in RG jargon) it is possible that on the scale of the correlation length the coarse-grained theory is driven near SBR by the RG flow. This explain why many different systems, including experimental ones, are potentially described by SBR and thus share the same qualitative features that does not depend on the actual values of the SBR parameters: notably power-law to exponential increase of the $\beta$ time and dynamical heterogeneities below the avoided transition.

SBR is thus potentially valid for $M=1$ as well, because long-range correlations may develop also in presence of short-range interactions. The deviations at smaller $M$ in the plots are indeed expected because the SBR parameters, being model-dependent, change with $M$. The values we computed in eqs. (8) have actually $1/M$ corrections that can be also computed systematically through a feasible but tedious power expansion. The only special feature of the large-$M$ regime is that the SBR parameters can be computed exactly from the Bethe lattice while the computation in the $M=1$ case is less straightforward (one should take into account the presence of small loops) but it is still feasible in principle.

While it is satisfying to compute the SBR parameter independently as we have done here, one could also extract some or all of them from fits. This means that SBR can be a useful tool to rationalize experimental data in the region where the widely used ideal MCT scalings fail. The outcome would still be highly non trivial because the SBR eqs. (6) yields predictions for many more quantities than those needed to determine the SBR parameters through fits. In particular here we considered only one temperature, but one may consider a whole range of temperatures Rizzo and Voigtmann (2015a) (corresponding to different values of $\sigma$ in eqs. (6)) and also study spatial correlations Rizzo and Voigtmann (2015a) and finite-size effects.

I Numerical Simulations

All simulations were performed at the critical temperature of the FA Model with $m=2$ on the Bethe lattice with connectivity $c=4$. The corresponding value can be obtained from the solution of Bootstrap percolation on the Bethe lattice given below. We considered the random lattice obtained by the application of the $M$-layer construction to the diamond cubic lattice with periodic boundary conditions. The diamond lattice is generated by repeating in the three directions a basic unit cell of length $L=4$. In each unit cell there are eight lattice points that can be divided into two groups: blue lattice points have coordinates $(0,0,0)$, $(2,0,2)$, $(0,2,2)$ and $(2,2,0)$, red lattice points have coordinates $(3,3,3)$, $(3,1,1)$, $(1,3,1)$ and $(1,1,3)$. Each red (blue) lattice point is connected with its four blue (red) nearest neighbors. We have used a Chessboard/Metropolis setting: all red spins are updated sequentially and then all the blue spins are updated sequentially. The single spin update is made with a Metropolis move: a negative mobile spin is flipped with probability $e^{-\beta}$ and a positive mobile spin is flipped with probability one. The chessboard setting is more convenient than the standard Random-site/Metropolis setting typically considered in the FA model as will be discussed below.

II Mean-Field Behavior

In figure (3) we plot $g(t)=\phi(t)-\phi_{plat}$ for a system with $L=4$ and $M=1.6\times 10^{6}$. In the range of times shown these data have converged on the Bethe lattice solution corresponding to $M\rightarrow\infty$.

In fig. (4) we plot parametrically the (discrete) logarithmic derivative of $g(t)$ vs. $g(t)$. For a $g(t)$ decaying at large times as $1/t^{a}$ the logarithmic derivative should converge to $-a$ in the limit of $g(t)\rightarrow 0$. The data display linear behavior in the small $g$ region corresponding to a correction of order $1/t^{2a}$. This is precisely the small-time correction that one would expect within MCT: the leading term is given by a quadratic equation while various subleading cubic terms induce a $1/t^{2a}$ correction. In full generality we can write

[TABLE]

where all possible constants in front of the leading term are absorbed in the definition of $t_{0}$. The above expression leads to:

[TABLE]

and thus from a linear fit we can extract $a$ and $\delta_{1}$. The constant $t_{0}$ can then be estimated fitting the numerical data with the asymptotic form (10), in fig. (3) we show both the $1/(t/t_{0})^{a}$ term and the corrected expression. We thus obtain the values quoted in the main text:

[TABLE]

From this we also have $\lambda=0.634$. We stress that in the range of times accessed numerically the $1/t^{2a}$ correction is significant and a linear fit of the data in fig. (3) would give an incorrect smaller exponent $a$. Instead the analysis of the second derivative is much safer and inconsistent with the value $a=.28-.3$ reported in earlier studies Sellitto et al. (2005); Sellitto (2015); De Candia et al. (2016). We also note that within the replica treatment the parameter exponent (and thus the exponent $a$) is available from the statics Caltagirone et al. (2012); Parisi and Rizzo (2013) and only the coefficient $t_{0}$ must be extracted numerically. Unfortunately a static replica treatment of the FA model is not available at present.

II.1 Comparison between different dynamics

The choice of the Chessboard setting is more convenient than the Random-Site/Metropolis dynamics typically used in the literature because each Monte-Carlo step (MCS) requires less CPU time, besides in MCS unit the relaxation is faster as can be seen in fig. (5) where we display the relaxation for a system with $L=4$ and $M=4\times 10^{5}$ with different dynamics, including Chessboard with heat-bath (Glauber) update. The key point is that at large enough times the different curves differ only by a constant shift in time: this can be seen more clearly considering the parametric plot of the logarithmic derivative of $g$ vs. $g$.

In fig. (6) we see that at small values of $g$ the parametric curves collapse onto a single curve independently of the dynamics. Note that the collapsing curve deviates from the linear fit corresponding to the asymptotic behavior implying that finite-size critical deviations are also independent on the dynamics. This is consistent with the fact that as we saw already these deviations are described by SBR in which the different microscopic dynamics enter only through the parameter $t_{0}$. We also mention that the Chessboard dynamic satisfies detailed balance but is not invertible, furthermore irreducibility of the Markov chain is not granted in general in sequential update with Metropolis moves, but these features, as usual, appear to be harmless in the interesting region of long-time critical behavior.

III Bootstrap Percolation on the Bethe lattice

The dynamics transition of FA models on the Bethe lattice of connectivity $c$ is intimately related with $k$-core, or bootstrap, percolation Sellitto et al. (2005). Let us recall the definition of BP: first the sites of a given lattice are populated with probability $p$, then each site with less than $k$ neighbors is removed and the culling process is repeated until each site has at least $k$ neighbors. When the process is completed the remaining occupied sites, if any, form the so-called the $k$-core.

One can argue that the cluster of blocked negative sites on the FA model is exactly the $k$-core with $k$ given by

[TABLE]

therefore $k=3$ for the $(c=4,m=2)$ FA models considered here. The relationship between temperature in the FA model and $p$ in BP is easily found:

[TABLE]

The outcome of the culling process leading to the $k$-core does not depend on the sequence of which sites are culled, in particular one can consider the modified problem in which a given site $i$ is occupied with probability one and no culling is applied to it and then apply an extraction and culling (eac) move to obtain its probability $P_{site}$ to be in the $k$-core. $P_{site}$ is given by $p$ times the sum over all configuration of its neighbors $\{s_{j},j\in\partial i\}$ satisfying the constraint of their occupancy probability $P^{(i)}(s_{j},j\in\partial i)$ before eac on site $i$. On the Bethe lattice with connectivity $c$ one can argue that the probability $P^{(i)}(s_{j},j\in\partial i)$ is factorized:

[TABLE]

where $s=-1$ if the site is occupied and $s=1$ otherwise and

[TABLE]

With the definition $P_{s,t}\equiv\sum_{i=s}^{t}{t\choose i}P^{i}(1-P)^{t-i}$ the Bethe solution is determined by the following equations:

[TABLE]

For all $k>2$ a solution with non-zero $P$ is found at large values of $p$; the solution disappears at a critical value $p=p_{c}$ with a square-root singularity thus exhibiting the celebrated mixed first-order/second-order character.

As shown in the Supplemental Material for Rizzo (2019) the connected probability that two points at large distances $L$ on the Bethe lattice are both on the $k$-core is given:

[TABLE]

where

[TABLE]

Note that in Rizzo (2019) the parameter $\mu$ is called $\lambda$ but we changed notation to avoid confusion with the MCT parameter exponent.

[TABLE]

[TABLE]

In the case $(c=4,k=3)$ we have

[TABLE]

From which the expressions for $c_{1}$ and $c_{2}$ given in section below can be obtained.

IV Computing the parameters of SBR

In this section we will compute the parameters $\alpha$, $\Delta\sigma$ and $\sigma$ of SBR for FA models on $M$-layer lattices by comparing the equal-time fluctuations of the density $g(x,t)$ in the Mean-Field regime. According to the analysis of Altieri et al. (2017) on the $M$-layer lattice the leading (mean-field) expression of the correlation of the order parameter at two points of the lattice is given by the correlation on all non-backtracking paths connecting the two sites on the original lattice divided by $M$:

[TABLE]

where $c_{L}^{Bethe}(t)$ is the correlation of the shifted persistence between two sites $i$ and $j$ at distance $L$ on the Bethe lattice:

[TABLE]

and $N_{L}(x,y)$ is the number of non-backtracking walks between point $x$ and point $y$ on the original lattice (corresponding to $M=1$). As we have discussed in the main text, for $M$ finite but large, in the regime where SBR provides a quantitatively accurate description, the length-scale of the fluctuations of $g(x,t)$ is large and thus we are interested in the regime where $x-y$ (and thus $L$) in $N_{L}(x,y)$ is also large. In this regime $N_{L}(x,y)$ tends to a Gaussian with a $O(L)$ variance:

[TABLE]

The above relationship is valid only if $x$ and $y$ corresponds to coordinates of points of the lattice and is zero otherwise. This explains the prefactor, indeed since as we will see the variance of the Gaussian is much larger than the lattice spacing we can replace the sum over lattice points as an integral on the continuum $\sum_{j}\rightarrow\int\rho\,d^{d}x)$. The integral on the other hand must be equal to the total number of paths of length $L$ originating from a point, i.e. $c(c-1)^{L-1}$. In terms of the unitary Fourier transform we can then write

[TABLE]

Where $D_{NBW}$ is by definition the diffusion coefficient of non-backtracking random walks $\{\mathbf{x}_{1},\dots,\mathbf{x}_{t}\}$ on the original lattice with $M=1$:

[TABLE]

In section IX we will provide a simple expression for $D_{NBW}$ in terms of the dimension of the lattice and of its connectivity valid for a huge class of lattices.

If we compute quantities at the level of the MF approximation we will get the spurious dynamical arrest transition. In particular, the FA model on the Bethe lattice below the critical temperature (corresponding to $p<p_{c}=8/9$) displays a glassy phase. In the glassy phase it is convenient to study the persistence and its fluctuations in the infinite-time limit. Furthermore in this limit the critical properties of the blocked (negative) sites are exactly the same of the $k$-core of bootstrap percolation. In particular we have

[TABLE]

and

[TABLE]

where $P_{site}(p)$ is the $k$-core density, $\delta p\equiv p-p_{c}$ and $c_{L}^{Bethe}(p)$ is the two-point correlation of the $k$-core. The above formula for $c_{L}^{Bethe}(p)$ (first appeared in Schwarz et al. (2006)) is derived in the supplemental material of Ref. Rizzo (2019) (see also section III above) where it is also shown that $\mu(p)$ tends to the critical value $(c-1)^{-1}$ at large times with a correction of order $\sqrt{\delta p}$

[TABLE]

and the numerical constants $c_{1}$ and $c_{2}$ read for $c=4$ and $m=2$ Rizzo (2019):

[TABLE]

In section III we wrote the general formulas from which $c_{1}$ and $c_{2}$ were obtained. Putting eqs. (23), (26) and (29) together and performing the summation over $L$ we can write for the unitary Fourier transform of the fluctuations:

[TABLE]

[TABLE]

We stress the difference from a simple Ornestein-Zernicke form due to the square that appears because of the $O(L)$ prefactor in (29). The above expression has to be compared with the long-time limit of the mean-field approximation to SBR in the glassy phase $\sigma>0$. The MF approximation to SBR corresponds to assume $\Delta\sigma^{2}$ is negligible and leads to:

[TABLE]

The MF approximation to fluctuations can be computed treating the fields $h(x)$ as a small perturbation and is given by:

[TABLE]

[TABLE]

By equating expressions (33) and (28) we obtain:

[TABLE]

By equating expressions (34) and (32) we obtain:

[TABLE]

[TABLE]

Thus the microscopic properties of the original lattice ($M=1$) enter through the density of lattice points $\rho\equiv N/V$ and through the diffusion coefficient on non-backtracking walks $D_{NBW}$. For the diamond lattice we have $\rho=1/8$ and $D_{NBW}=1$ as can be obtained from a general formula derived in section IX. Replacing $\rho=1/8$, $D_{NBW}=1$ and $\lambda=.634$ (from the numerical solution on the Bethe lattice discussed before) we have obtained the values quoted in the main text.

Note that although we have considered the three-dimensional diamond lattice the previous relationships are valid for a generic lattice of connectivity four and can be used to repeat the numerical analysis on many other lattices, including notably the regular two-dimensional lattice where $\rho=1$ and $D_{NBW}=1/2$. Formulas for lattices with generic connectivity and topology can be obtained from section II and IX.

Note that the MF approximations we have discussed are utterly wrong because neither the FA model on the $M$-layer nor the SBR equation actually display a dynamic arrest transition and there is no static limit. Nevertheless we are making in both cases the same (wrong) approximation and it is thus correct to compare the outcome.

V Dynamical Mean-Field Fluctuations

In the previous section we have computed the mean-field expression (23) of the fluctuations in the glassy phase and in the infinite time limit. In the following we will obtain its expression at large times at the critical point ($p=p_{c}=8/9$). In this case one expects (and may confirm numerically) that:

[TABLE]

where $\mu(t)$ tends to the critical value $(c-1)^{-1}$ at large times with a correction of order $t^{-a}$

[TABLE]

The constants $c_{1}$ is the same obtained above in the glassy phase, while $\tilde{c}_{2}$ is different from $c_{2}$ and at present can only be extracted from the numerics. Repeating the same steps above we then obtain:

[TABLE]

this expression has to be compared with the mean-field small-time approximation of the SBR equations. This is obtained solving the equations perturbatively around the small-time limit $\hat{g}(x,t)\approx(t/t_{0})^{-a}$. In Fourier space one computes the correction due to the $s(x)$, then after averaging one easily obtains:

[TABLE]

where we have used Gotze’s definition (eq. 6.63a in Götze (2008))

[TABLE]

Note the above expression will be modified at larger times and should be replaced with the full SBR average. Equating (40) and (41) provides an alternative way to determine the SBR parameters $\Delta\sigma$ and $\alpha$ from the dynamics of FA models on the Bethe lattice (through the constants $c_{1}$ and $\tilde{c}_{2}$).

Setting for simplicity $\Delta\sigma^{2}=1/M$, $\alpha=t_{0}=1$ we obtain

[TABLE]

Thus at the MF level the total susceptibility diverges with time as $t^{2a}$ and the correlation length diverges as $t^{a/2}$. Both these behaviors will change when the system starts to deviate from MF: the correlation length will not increase indefinitely and $g(x,t)$ will become negative at a finite time. In this regime corrections of all orders in $1/M$ become equally relevant and one must abandon MF theory, technically using SBR in place of the MF espressions eq. (43) amounts to include corrections at all orders in powers of $1/M$ Rizzo (2016). We note en passant that MF fluctuations do not have the Ornstein-Zernicke form but rather the square of it. Thus it should be noted that, while the use of the OZ form to fit numerical data of dynamical fluctuations is widespread in the literature Lačević et al. (2003); Berthier et al. (2007); Flenner and Szamel (2013); Karmakar et al. (2014), it has no theoretical justification. Both MF theory and SBR certainly do not support OZ in the $\beta$ regime, while the $\alpha$ regime predictions are at present not available, not even in MF theory. In real space the above expression leads to:

[TABLE]

that has been used in the main text to discuss the time-dependent Ginzburg criterion. From the above expressions we also obtain the MF behavior of the $\chi_{4}(t)$ susceptibility at large times is

[TABLE]

that gives the coefficient of the MF asymptote discussed in the main text.

VI The SBR equations in universal form

The equations of SBR are Rizzo (2014):

[TABLE]

where the field $s(x)$ is a time-independent random fluctuation of the separation parameter, Gaussian and delta-correlated in space:

[TABLE]

and they have to be solved with the condition

[TABLE]

Thus the equations depend on five coupling constants: $\lambda$, $t_{0}$, $\alpha$, $\sigma$ and $\Delta\sigma$. However one can fix $t_{0}=\alpha=\Delta\sigma=1$ and then the general solution can be obtained by rescalings. This means that in practice the SBR equations at fixed $\lambda$ need only to be solved for varying values of $\sigma$. More precisely one can easily verify that for a generic $K$-point function we have

[TABLE]

where the notation $[\dots]_{\alpha,\Delta\sigma,t_{0},\sigma,\lambda}$ means that the SBR equations above are to be solved with the corresponding values of the five parameters. The rescaling parameters read:

[TABLE]

VII The interplay between $M$ and the system size $L$

In the main text we saw that in the MF regime valid at initial times the correlation length grows with time as $t^{a/2}$ while the time where SBR replaces MF theory grows with $M$ as $O(M^{{1\over a(4-d/2)}})$. Thus if the correlation length becomes comparable with the system size before deviations from MF occur the SBR equations do not have significant spatial variation and can be replaced by the SBR equations without the gradient term that will be written below. This situation will always occur if we increase $M$ keeping the size $L$ of the original ($M=1$) lattice fixed. This can be seen starting from the fact (see section VI) that the solution of SBR corresponding to generic values of $\alpha$, $\Delta\sigma$ and $t_{0}$ can be expressed in terms of the solution with $\alpha=\Delta\sigma=t_{0}=1$ through appropriate rescalings of the correlators ($g\rightarrow b_{\phi}\,g$), of distances ($x\rightarrow b_{x}\,x$) and of times ($t\rightarrow b_{t}\,t$).

From the corresponding formulas we get that the $M$-dependence of the rescaling factors is the following:

[TABLE]

for generic dimension $d$ smaller than the critical dimension $d_{c}=8$. The above relationships allow to understand the phenomenology on the $M$-layer lattice at the avoided MCT transition. They tell us that at $T_{c}$ deviations from the Bethe lattice behavior occur at times that increase with $M$ as $O(M^{{1\over a(4-d/2)}})$. Furthermore at times $O(M^{{1\over a(4-d/2)}})$ the order parameter, i.e. the deviations of the persistence from the plateau value, is small $O(M^{-{1\over 4-d/2}})$ and fluctuates over a large length-scale $O(M^{{1\over 8-d}})$, thus justifying a posteriori the validity of an effective theory description and the expectation that SBR predictions becomes increasingly accurate for larger values of $M$. In the main text we have given instead an a priori justification, starting from the mean-field expression of fluctuations in the $M$-layer and applying a time-dependent Ginzburg criterion that also lead to the above scaling eq. (54).

According to the previous formulas the system size $L$ corresponds to an effective system size $L/b_{x}=L\,M^{-{1\over 8-d}}$ in the universal theory. Thus if we if we go from $M$ to a larger value $M^{\prime}$ while keeping the system size $L$ fixed we are systematically reducing the effective system size until we have a system so small that spatial fluctuations of $g(x,t)$ are negligible because the gradient term suppresses fluctuations on small distances. In this regime the only source of fluctuations is the fluctuations of $\sigma$. The system becomes effectively a zero-dimensional system in the sense that the predictions of the finite-dimensional SBR equations are indistinguishable from those of the the zero-dimensional SBR equations discussed below, eqs. (57).

To complete the discussion on the interplay between $M$ and $L$ it is interesting to consider what happens while changing $L$ at fixed large $M$. When $L=\infty$ we are in the thermodynamic limit and this is so also in terms of the effective system size $L/b_{x}$. For $L<\infty$ finite-size effects will be significant as soon as $L$ is of the order of the correlation length. In this regime the solution of the SBR equations depends on the size of the system and it is important to solve them in a box of the correct size to get accurate results. If we decrease the system size further we will eventually reach the zero-dimensional limit: the spatial fluctuations of the solution are so small that the finite-dimensional SBR equation can be replaced by the zero-dimensional equation discussed below for all practical purposes.

The data for $L=4$ and values of $M=50000$, $M=100000$, $M=200000$ are actually in the zero-dimensional regime and the corresponding SBR predictions plotted of fig. 1 in the main text were obtained indeed from the zero-dimensional SBR equation. The data for $M=3000$ were obtained with $L=8$ an are instead strictly three-dimensional because, although the system size is finite, we had to consider the full three-dimensional SBR equation that depends on the value of $D_{NBW}$ on the lattice and on the system size.

VIII Finite-Size Corrections to Mean-Field Models (SBR in Zero Dimension)

The zero-dimensional equations can be written solely in terms of the total shifted persistence

[TABLE]

Then we have for a generic cumulant of order $K$:

[TABLE]

where in the RHS $g(t)$ is the solution of of the zero-dimensional SBR equations Rizzo (2014, 2016); Rizzo and Voigtmann (2015a):

[TABLE]

The square brackets mean average with respect to the field $s$ that is a time-independent Gaussian random fluctuation of the separation parameter:

[TABLE]

The equations have to be solved with the small-time condition

[TABLE]

Note that the gradient term is not present and indeed the actual geometrical structure of the lattice is irrelevant, the only relevant control parameter for SBR is the total number $N_{tot}=M\,L^{3}\,\rho$ of sites of the lattice that controls $\Delta\sigma^{2}$ through:

[TABLE]

that follows from expression (36). The above expression does not depend on how the lattice is actually generated and in particular it is also correct for a random-regular-graph of size $N_{tot}$ Franz and Sellitto (2013). In the present setting $N_{tot}=M\,L^{3}\,\rho=8\,M$.

As we said mentioned in previous section the data for $L=4$ and values of $M=50000$, $M=100000$, $M=200000$ are in the zero dimensional regime. In fig. (7) we rescaled them horizontally and vertically with the appropriate powers of $M$ that lead to collapse on a unique SBR curve. The rescaling factors are $M^{1/4}$ for the shifted persistence and $M^{-1/(4a)}$ for time and are obtained setting $d=0$ in eqs. (54), which is the technical reason why finite-size corrections correspond to zero dimension.

Note that the agreement between data and SBR in the scaling plot of the figure is not yet perfect and this can be tracked to the fact that, at the values of $M$ considered, the curves start to deviate from the mean-field curve at times where the $M=\infty$ mean-field curve itself has still relevant small-time corrections of order $(t/t_{0})^{-2a}$ to the leading asymptotic behavior $(t/t_{0})^{-a}$, as discussed in section II. As a consequence the data for small rescaled times $tM^{1/(4a)}$ of order $5-20$ display an approximate $20\%$ deviation from the asymptotic expression $(t/t_{0})^{-a}$ that should describe the data for small rescaled times at large values of $M$ but the agreement improves considering larger values of $M$. On the other hand small-time corrections gets smaller at larger rescaled times $\approx 50$ and correspondingly the data agree much better with the SBR curve. We note that the rescaling factor $M^{1/4}$ that leads to the collapse of the curves was obtained for the first time in Franz et al. (2011) from a static computation, but the resulting scaling function, apart from being independent of time, is ill-defined and one needs a full-fledged dynamical treatment to compute the time-dependent scaling curve.

IX Diffusion Coefficient of Non-Backtracking random walks on generic lattices

We consider lattices with connectivity $c$ such that we can classify sites in two classes, say red and blue, such that a blue site is connected to the $c$ red sites located in the $c$ directions ${\mathbf{v}}_{\mu}$. Similarly a red site is connected to the $c$ blue sites located in the $c$ directions $-{\mathbf{v}}_{\mu}$. We will not require the condition $\{{\mathbf{v}}_{\mu}\}=-\{{\mathbf{v}}_{\mu}\}$, meaning that the set of directions need not to be invariant under inversion of the coordinate axes. Thus we extend previous results Altieri et al. (2017); Fitzner and van der Hofstad (2013) to include e.g. the honeycomb lattice in $d=2$ and the diamond lattice in $d=3$. We assume instead that

[TABLE]

A non-backtracking random walk can be described as a sequence of steps at consecutive times $s$ in the directions $\mu(s)$:

[TABLE]

where the minus sign comes from the fact that the set of possible directions changes as the walker moves from a blue and to a red site. The average can then be written in terms of the joint probability $P^{(s,s^{\prime})}_{\mu,\mu^{\prime}}$

[TABLE]

note that we use Einstein’s convention of implicit summation over repeated indexes $\mu$. Since the directions are uniformly distributed at zero time we have exactly

[TABLE]

and we can write for large $t$

[TABLE]

The probability can be computed recursively

[TABLE]

where $T_{\mu,\mu^{\prime\prime}}$ is the $c\times c$ matrix with zero diagonal elements and off-diagonal elements equal to $1/(c-1)$. We can thus write

[TABLE]

where ${\mathcal{P}}_{s}$ is the projector on the vector with all equal components and ${\mathcal{Q}}_{s}=I-{\mathcal{P}}_{s}$ is the orthogonal projector. It follows that

[TABLE]

since ${\mathcal{P}}_{s}$ is a matrix with all elements equal it gives zero contribution when summed over the directions because of the condition $\sum_{\mu}{\bf\hat{v}}_{\mu}=0$ and we have

[TABLE]

replacing the above expression in the large time expression eq. (65) we obtain:

[TABLE]

and therefore (taking into account that each vector has $d$ components)

[TABLE]

On regular lattice $c=2d$ and we recover the result $D_{NBW}=1/(2\,d-2)$. On the diamond cubic lattice studied in the paper we have $c=4$ and $v^{2}=3$ leading to $D_{NBW}=1$. The honeycomb lattice can be realized repeating a $3\times\sqrt{3}$ unit cell with two blue points at coordinates $(0,0)$,$(\sqrt{3}/2,3/2)$ and two red points with coordinates $(0,1)$, $(\sqrt{3}/2,5/2)$. This leads to $\rho=4\,3^{-3/2}$ and $D_{NBW}=3/4$.