Using matrix product states to study the dynamical large deviations of kinetically constrained models
UUsing matrix product states to study the dynamical large deviations of kineticallyconstrained models
Mari Carmen Ba˜nuls
1, 2 and Juan P. Garrahan
3, 4 Max-Planck-Institut f¨ur Quantenoptik, Hans-Kopfermann-Str. 1, D-85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 M¨unchen School of Physics and Astronomy, University of Nottingham, Nottingham, NG7 2RD, UK Centre for the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems,University of Nottingham, Nottingham, NG7 2RD, UK (Dated: April 19, 2019)Here we demonstrate that tensor network techniques — originally devised for the analysis ofquantum many-body problems — are well suited for the detailed study of rare event statisticsin kinetically constrained models (KCMs). As concrete examples we consider the Fredrickson-Andersen and East models, two paradigmatic KCMs relevant to the modelling of glasses. We showhow variational matrix product states allow to numerically approximate — systematically and withhigh accuracy — the leading eigenstates of the tilted dynamical generators which encode the largedeviation statistics of the dynamics. Via this approach we can study system sizes beyond whatis possible with other methods, allowing us to characterise in detail the finite size scaling of thetrajectory-space phase transition of these models, the behaviour of spectral gaps, and the spatialstructure and “entanglement” properties of dynamical phases. We discuss the broader implicationsof our results.
Introduction.–
Dynamics equipped with local kineticconstraints provides a general mechanism for slow coop-erative relaxation [1–4]. Kinetically constrained models(KCMs) — of which the Fredrickson-Andersen (FA) [2]and East [3] facilitated spin models are the simplest ex-ponents — give many insights into the nature of glassforming systems, in particular by showing that systemswith simple thermodynamics can have rich, spatially fluc-tuating and slow dynamics [5]. (For reviews on the glasstransition see [6–8], and on KCMs see [9–11].) Beyondglasses, classical KCMs (and related deterministic models[12–16]) are relevant to the problem of operator spread-ing in quantum systems [17–24] and to non-equilibriumdynamics of ensembles of Rydberg atoms [25–27], whilequantum KCMs provide a template for complex non-equilibrium dynamics under unitary evolution in the ab-sence of disorder [28–31].To characterise dynamics it is natural to study ensem-bles of stochastic trajectories, just like one does in equi-librium statistical mechanics with ensembles of configura-tions. For long-times one can then apply the methods ofdynamical large deviations (LDs) [32] to compute quan-tities that play the role of thermodynamic potentials forthe dynamics. For the case of KCMs this “thermody-namics of trajectories” approach reveals the existence ofa first-order phase transition in the space of trajectoriesbetween active and inactive dynamical phases, indicativeof the singular change when fluctuating away from typ-ical behaviour [33, 34]. Many other systems have beenalso shown to have similar LD transitions, see e.g. [35–44]. Understanding the phase structure of the dynamicsis clearly as important in dynamical problems as it is instatic ones.The standard way of accessing LD statistics of a dy- namical observable is by computing its scaled cumulantgenerating function (SCGF) — see below for definitions— from the largest eigenvalue of an appropriate deforma-tion, or tilting , of the generator of the dynamics [11, 32].Except for the handful of non-trivial cases in which itcan be calculated exactly [16, 36], obtaining the SCGFby diagonalising the tilted generator is only possible forsmall system sizes. To access the LD behaviour for largersizes one has to resort to numerical methods for samplingrare trajectories based on splitting/cloning, importancesampling or optimal control [37, 45–51].By exploiting the similarity between tilted generatorsand quantum Hamiltonians, here we show how to usevariational matrix product states (MPS) to compute nu-merically with high accuracy (and precise control on er-rors) leading eigenvalues and eigenstates of the tiltedgenerator for system sizes way beyond those accessiblethrough other methods. We study in detail the FA andEast models, focusing on the finite size scaling of theiractive-inactive phase transitions and the spatial struc-ture that emerges in the dynamical phases. While in cer-tain special cases MPS can be used to obtain exact LDstatistics, such as in simple exclusion processes [52–56],hard core brownian particles [57], and certain cellularautomata [16], the systematic application of numericalMPS methods to stochastic lattice systems has been lim-ited [58]. Our results for KCMs — together with thevery recent ones [59] for simple exclusion processes —show the potential of numerical tensor network meth-ods for the detailed study of dynamical fluctuations instochastic dynamics.
FA and East models.–
The FA [2] and East [3] modelsare defined in terms of binary variables, { n i = 0 , } Ni =1 ,on the sites of a one dimensional lattice of size N , with a r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r FIG. 1:
Finite size scaling of trajectory transition. (a) SCGF θ ( s ) /N as a function of s for the East model( c = 0 .
2) for system sizes N = 20 to N = 200. The critical s c ( N ) can be obtained from the (extrapolated) crossingof the first two energy levels. The dot-dashed lines correspond to the asymptotic values θ ( s → ∞ ) = − c . (b) Thecorresponding dynamical susceptibilities, χ ( s ) = θ (cid:48)(cid:48) ( s ), exhibit a peak at s c ( N ) that gets narrower and higher as N increases. For s > s c ( N ) we find an almost universal behavior χ ∝ s − γ with γ ≈ .
4. (c) s c ( N ) as a function of N for N ∈ [20 , c in the East model (top) and FA model (bottom). Asexpected the data is compatible with lim N →∞ s c ( N ) →
0, but s c appears to scale as s c ( N ) ∝ N − α with α > a/N + b/N , dashed). (d) The scaling exponents α (bluediamonds) and fitting parameters b/a (red squares) as a function of c (top, East model; bottom, FA model). Thedeparture from 1 /N scaling (dotted-dashed) appears to be more pronounced the lower the c is. (e) Rate functions ϕ ( k ) for N ∈ [20 , c = 0 . c = 0 .
05 (bottom). Thedashed lines correspond to Poisson distributions with average (cid:104) k (cid:105) = − θ (cid:48) (0) /N .single-spin flip dynamics subject to a kinetic constraintsuch that a spin can flip up (with rate c ) or down (withrate 1 − c ) only if either nearest neighbour is in the upstate (FA model) or only if the leftmost nearest neighbouris in the up state (East model). The generators for thecorresponding continuous time Markov chains are [9–11] W FA = (cid:88) i ( n i − + n i +1 ) (cid:2) cσ + i + (1 − c ) σ − i − c (1 − n i ) − (1 − c ) n i ] , (1) W East = (cid:88) i n i − (cid:2) cσ + i + (1 − c ) σ − i − c (1 − n i ) − (1 − c ) n i ] , (2)where σ ± i flips the site i up/down, and the factor in frontof the square brackets is the kinetic constraint. In thisformulation the master equation is ∂ t | P (cid:105) = W | P (cid:105) , where | P (cid:105) is the probability vector over configurations.We consider open boundary conditions which formallycorresponds to setting n = n N +1 = 0 in Eqs. (1) and (2).This is the best setup for the MPS method we use below.Due to the kinetic constraints configuration space can bedisconnected, and we consider the dynamics within thelargest ergodic component: the set of all configurationswith at least one up site for the FA model, and all theconfigurations with fixed n = 1 for the East model.The dynamics has as stationary distribution | P eq (cid:105) given by a projection of the product state | c (cid:105) ⊗ N , where | c (cid:105) = (1 − c ) | (cid:105) + c | (cid:105) , into the relevant ergodic compo-nent, giving | P FAeq (cid:105) = [ | c (cid:105) ⊗ N − (1 − c ) N | (cid:105) ⊗ N ] / [1 − (1 − c ) N ] , (3) | P Easteq (cid:105) = | (cid:105) ⊗ | c (cid:105) ⊗ N − . (4)These are the equilibrium distributions with energy E = (cid:80) i n i at inverse temperature ln(1 − c ) /c in the corre-sponding ergodic components. Dynamical LDs and tilted generators.–
As trajec-tory observable we will consider the dynamical activity [33, 35, 60], given by the total number of configurationchanges K ( ω t ) (i.e., number of spin flips) in a trajectory ω t of time extent t . For large t the distribution of K obeysa LD principle, P t ( K ) = (cid:104) δ [ K ( ω t ) − K ] (cid:105) ≈ e − tϕ ( K/t ) ,where ϕ ( x ) is the LD rate function [32]. The correspond-ing moment generating function Z T ( s ) = (cid:104) e − sK ( ω t ) (cid:105) alsoobeys a LD principle, Z T ( s ) ≈ e tθ ( s ) , where θ ( s ) isthe scaled cumulant generating function (SCGF), whosederivatives at s = 0 give the cumulants of K (scaled by t ) [32]. The LD functions are connected by a Legendretransform, θ ( s ) = − min k [ sk + ϕ ( k )] [32] and play therole of thermodynamic potentials for trajectories.The SCGF can be obtained from the largest eigen-value of a tilted generator, W s [32]. For the case of thedynamical activity, the tilt corresponds to multiplyingthe off-diagonal terms of W by a factor e − s [33, 35].Since the dynamics obeys detailed balance, the gener-FIG. 2: Structure of active phase. (a) Mean density (cid:104) n (cid:105) s in the active phase, s <
0, in the East model for thecase c = 0 .
05 (shown as function of − ν = e s − c the plateau structure of the density is evident (ascompared to c = 0 . H s at ν = 0 .
081 ( s = − . c = 0 .
05 for sizes N = 20 ,
100 (top andmiddle panels) and density profiles across the active phase for N = 20 (bottom panel). (d) Same for the FA modelat c = 0 .
05. For the East model the state has pronounced anticorrelations which are absent in the FA case model.(e) Extreme limit of the active phase, s → −∞ , for the East and FA models (top and bottom, respectively). In thepanels on the right the symbols show the rescaled ˜ θ ( s = −∞ ) /N := e s θ ( s = −∞ ) / [ N (cid:112) c (1 − c )] (black circles), (cid:104) n (cid:105) s = −∞ (blue squares) and (cid:104) n x (cid:105) s = −∞ (green diamonds) for N ∈ [20 , a/N + b to extractthe values in the thermodynamic limit: lim N →∞ θ ( s = −∞ ) /N, (cid:104) n x (cid:105) s = −∞ , (cid:104) n (cid:105) s = −∞ = 0 . , . , .
75 (East) and1 . , . , .
75 (FA). The right panels show (for N = 20) that the density profiles at s = −∞ are uniform, up toboundaries, in both models.ators can be made hermitian by a similarity transforma-tion which is independent of s [34]. That is, if we define H s = − Q − W s Q , where Q is a diagonal matrix withelements (cid:104) n | Q | n (cid:105) = (1 − c ) N/ [ c/ (1 − c )] (cid:80) i n i / in theconfiguration basis {| n (cid:105)} , we get H FA s = − (cid:88) i ( n i − + n i +1 ) (5) × (cid:104) e − s (cid:112) c (1 − c ) σ xi − c (1 − n i ) − (1 − c ) n i (cid:105) ,H East s = − (cid:88) i n i − (cid:104) e − s (cid:112) c (1 − c ) σ xi (6) − c (1 − n i ) − (1 − c ) n i (cid:105) , The SCGF therefore corresponds to (minus) the groundstate energy of H s , θ ( s ) = − E GS ( s ) . (7)The relation between the ground state | Φ s (cid:105) of thetilted Hamiltonian, H s | Φ s (cid:105) = E GS ( s ) | Φ s (cid:105) , and the left (cid:104) L s | and right | R s (cid:105) leading eigenvectors of the tilted gen-erator, W s | R s (cid:105) = θ ( s ) | R s (cid:105) , (cid:104) L s | W s = (cid:104) L s | θ ( s ), is | Φ s (cid:105) = (cid:88) n (cid:112) l n ( s ) r n ( s ) | n (cid:105) (8)where l n ( s ) = (cid:104) L s | n (cid:105) and r n ( s ) = (cid:104) n | R s (cid:105) . The aim nowis to compute E GS ( s ) and | Φ s (cid:105) for Eqs. (5) and (6). Variational MPS method.–
For a lattice of
N d -dimensional quantum systems, a MPS [61] is a vector | Ψ (cid:105) = (cid:80) di ,...i N =1 tr (cid:0) A i A i . . . A i N N (cid:1) | i i . . . i N (cid:105) , where i k labels a local basis of the k − th subsystem, and each A k is a rank-3 tensor of dimensions d × D × D [62].Such a state is described by O ( dN D ) parameters. The bond dimension D limits the entanglement of the state.More precisely, in an MPS of bond dimension D , forany subchain A the entanglement entropy (defined as S E = − Tr A ρ A log ρ A , where ρ A = Tr N \ A | Ψ (cid:105)(cid:104) Ψ | [63])is upper-bounded by S E ≤ D , independent of thesubchain length. Namely, MPS satisfy an entanglement area law [64], and conform a hierarchy of increasingly en-tangled states, with D = d N/ sufficing to describe thewhole Hilbert space.Conversely, MPS can efficiently approximate statesthat satisfy an area law [65], such as ground states ofgapped local Hamiltonians. They thus are the basisfor numerical methods like the celebrated density matrixrenormalization group (DMRG) algorithm [66] which canbe understood as a variational minimization of energyover MPS [67–71], by sequientially optimizing a singletensor, while keeping the rest constant, and iterativelysweeping over the chain until convergence [72]. Formu-lated in terms of tensor networks this algorithm allowsa number of extensions, including simulating dynamics,and the calculation of a few excited states above theground state.FIG. 3: Entanglement. (a) Half-chain entanglemententropy S E of the ground state of H s as a function of s for c = 0 . , . , .
05 in the East model at N = 200. (b) S E for s < c = 0 . N . The peak iscorrelated with the change in shape of the spectral gap∆ E of H s shown in (c).We apply this strategy to find MPS approximations tothe ground state and first excitations of the Hamiltonians(5) and (6). In this case, d = 2 and the basis is {| n (cid:105)} . Aswe show below, MPS with D (cid:28) N provide accurate ap-proximations for systems sizes at an order of magnitudelarger than those accessible by other methods [73]. Results. Finite size scaling of active-inactive tra-jectory transition.–
The key property of KCMs like theFA and East is their first-order phase transition betweenan active phase for s < inactive dynamical phaseat s > N → ∞ . Like for all phase tran-sitions, to characterise the transition and its associatedfluctuations, it is necessary to understand how the singu-larity is approached as the system size increases. Theo-retical and numerical considerations [74–76] suggest thatfor finite N the (rounded) transition occurs at s c ( N ) > s = 0, is perturbatively connectedto the active phase), and s c ( N ) → + as 1 /N . Thesepredictions can be tested with our MPS method.Figure 1(a) shows (minus) the energy density − E GS ( s ) /N = θ ( s ) /N of the MPS solution as a func-tion of s . The transition at s c ( N ) occurs where thetwo branches cross. The leftmost branch is linear in s and proportional to N , corresponding to the linear re-sponse for s (cid:38) s (cid:38) θ ( ∞ ) = − c . The corresponding suscep-tibility χ s = θ (cid:48)(cid:48) ( s ) /N shows a diverging peak at s c ( N ),see Fig. 1(b) [73].We can estimate the location of s c ( N ) from the sus-ceptibility peak. For both models we find a departurefrom the expected 1 /N scaling. Figure 1(c) shows that s c ( N ) can be fit to a power law, s c ( N ) ∝ N − α with α > /N , see Fig. 1(c)(dashed lines), and Fig. 1(d) for the dependence of the scaling parameters with c . We also show in Fig. 1(e) thebroadening with N of the LD rate function, indicative ofthe first order transition [33, 34]. For more details on thefinite size scaling analysis including comparison with thepredictions of the Ref. [74] see [77]. Structure of active phase.–
While both models havesimilar active-inactive transitions, their active phases dif-fer. Figures 2(a,b) show the average density of excita-tions, (cid:104) n (cid:105) s = N − (cid:80) Ni =1 (cid:104) Φ s | n i | Φ s (cid:105) , in the MPS thatapproximates the ground state of H s for s <
0. In theEast model and for small c , (cid:104) n (cid:105) s shows a series of plateausas s becomes more negative, as predicted in Ref. [78].These plateaus are absent in the FA model at the same c , Fig. 2(b), and also when the equilibrium concentration c is high, see insets to Figs. 2(a,b).Figures 2(c,d) show the difference in spatial structureof the active phases. The top two panels in Fig. 2(c,d)give the density profile at s = − . ν = 0 . (cid:104) n (cid:105) s ≈ /
3. For the East model, Fig. 2(c, top two panels),the state is anticorrelated in space, with an occupied sitefollowed by two nearly empty ones. This is evident in the N = 20 case, shown in the figure, while for N = 100 wealso observe a longer ranged modulation of this pattern[73]. In contrast, in the FA model the density is essen-tially uniform, Fig. 2(d, top two panels). This differencein structure is present throughout the s < s → −∞ . We find that a MPS of D ∼ O (10)is enough to obtain a very precise approximation tothe ground state over the whole range of sizes com-puted, N ∈ [20 , N → ∞ . We obtain, Fig. 2(e), for the limiting SCGFslim N →∞ lim s →−∞ e s θ E ( s ) / [ N (cid:112) c (1 − c )] ≈ . .
337 for the FA model, while the densitiesare the same in both models, namely lim N →∞ (cid:104) n (cid:105) −∞ ≈ .
754 and lim N →∞ (cid:104) n x (cid:105) −∞ ≈ .
824 (where n x is up toconstants the “transverse” magnetisation, 2 n x = 1 − N − (cid:80) Ni =1 σ xi ). The panels on the right of Fig. 2(e) showthat the corresponding density profiles are essentially flatin this limit [79]. Entanglement.–
The states at s (cid:54) = 0 have spatial cor-relations absent in equilibrium ( s = 0) and which varieswith s . This can be quantified via their entanglemententropy, which together with other quantum informationmeasures can capture changes in dynamical behaviourthat might escape classical order parameters [80]. Theentanglement entropy is easily computed for a state inMPS form. Figure 3(a) shows the half-chain S E of thestate | Φ s (cid:105) as a function of s in the East model at size N = 200. It is zero in the equilibrium state, cf. Eq. (4),and very small in the inactive phase, where the leadingeigenvector is close to a product state of all sites emptyin the bulk. For s < S E does notseem to scale with system size. Thus, in the languageof quantum many-body systems, the ground state ful-fils an area law. This is also the case for other entropicquantities [73], which justifies the accuracy of the MPSapproximation.The peak in S E nevertheless is sensitive to changes inthe structure of the active phase. Fig. 3(c) shows thecorresponding gap between E GS ( s ) and the eigenvalueof the first excited state: its s dependence changes at avalue of s located by the peak in S E . (Note also that thegap is has no significant N dependence.) The maximumof the entropy depends on the value of c , and we finda larger peak for smaller values, corresponding to richerstructure in the active phase, see Fig. 3(a) and [73].Even if the entanglement is low throughout the phasediagram, cf. Fig. 3(a), this does not guarantee that thevariational method will easily find an MPS approxima-tion. In fact, we find that both for the region close to thephase transition at s = 0 and for the values of s where S E shows a peak, cf. Fig. 3(a,b), the numerical convergenceis slower than would have been expected. We believe thisis a consequence of how the spectrum of the Hamiltonianchanges when approaching these regimes [73]. Conclusions.–
As we have shown here, the MPS meth-ods often employed in quantum many-body problems[71], are also well suited for the study of the dynami-cal generators of classical stochastic systems [12–16, 52–59]. We focused on the LD statistics of KCMs such asthe FA and East models, and showed how variationalMPS approximations allow to efficiently access systemsizes which are larger by an order of magnitude com-pared to previous studies, thus providing detailed infor-mation about the properties of the transitions in thesemodels and the nature of the dynamical phases. We fore-see many other applications of tensor networks in clas-sical stochastic dynamics, including when the dynamicaltransition is continuous rather than first-order, and inthe study of systems in dimension larger than one. Morebroadly, the crossover of ideas and techniques betweenquantum many-body and classical stochastics remains afruitful area of investigation.
Acknowledgements –
This work was supported by theDeutsche Forschungsgemeinschaft (DFG, German Re-search Foundation) under Germany’s Excellence Strat-egy – EXC-2111 – 390814868, by EPSRC Grant No.EP/R04421X/1 and by the Leverhulme Trust Grant No.RPG-2018-181. We acknowledge the hospitality of theKavli Institute for Theoretical Physics at the Univer-sity of California, Santa Barbara, where this work wasstarted, and support from the National Science Founda-tion under Grant No. NSF PHY-1748958. [1] R. G. Palmer, D. L. Stein, E. Abrahams, and P. W.Anderson, Phys. Rev. Lett. , 958 (1984).[2] G. H. Fredrickson and H. C. Andersen, Phys. Rev. Lett. , 1244 (1984).[3] J. J¨ackle and S. Eisinger, Z. fur Phys. B , 115 (1991).[4] W. Kob and H. C. Andersen, Phys. Rev. E , 4364(1993).[5] J. P. Garrahan and D. Chandler, Phys. Rev. Lett. (2002).[6] K. Binder and W. Kob, Glassy materials and disor-dered solids: An introduction to their statistical mechan-ics (World Scientific, 2011).[7] L. Berthier and G. Biroli, Rev. Mod. Phys. , 587(2011).[8] G. Biroli and J. P. Garrahan, J. Chem. Phys. ,12A301 (2013).[9] F. Ritort and P. Sollich, Adv. Phys. , 219 (2003).[10] J. P. Garrahan, P. Sollich, and C. Toninelli, in Dynam-ical Heterogeneities in Glasses, Colloids, and GranularMedia , International Series of Monographs on Physics,edited by L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipel-letti, and W. van Saarloos (Oxford University Press,Oxford, UK, 2011).[11] J. P. Garrahan, Physica A , 130 (2018).[12] T. Prosen and C. Mej´ıa-Monasterio, J. Phys. A ,185003 (2016).[13] A. Inoue and S. Takesue, J. Phys. A , 425001 (2018).[14] T. Prosen and B. Buˇca, J. Phys. A , 395002 (2017).[15] K. Klobas, M. Medenjak, T. Prosen, and M. Vanicat,arXiv:1807.05000 (2018).[16] B. Buˇca, J. P. Garrahan, T. Prosen, and M. Vanicat,arXiv:1901.00845 (2019).[17] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys.Rev. X , 031016 (2017).[18] D. A. Rowlands and A. Lamacraft, Phys. Rev. B ,195125 (2018).[19] X. Chen and T. Zhou, arXiv:1808.09812 (2018).[20] S. Gopalakrishnan, Phys. Rev. B , 060302 (2018).[21] M. Knap, Phys. Rev. B , 184416 (2018).[22] M. C. Tran, A. Y. Guo, Y. Su, J. R. Garrison, Z. El-dredge, M. Foss-Feig, A. M. Childs, and A. V. Gorshkov,arXiv:1808.05225 (2018).[23] S. Gopalakrishnan, D. A. Huse, V. Khemani, andR. Vasseur, Phys. Rev. B , 220303 (2018).[24] V. Alba, J. Dubail, and M. Medenjak, arXiv:1901.04521(2019).[25] I. Lesanovsky and J. P. Garrahan, Phys. Rev. Lett. ,215305 (2013).[26] A. Urvoy, F. Ripka, I. Lesanovsky, D. Booth, J. P. Shaf-fer, T. Pfau, and R. L¨ow, Phys. Rev. Lett. , 203002(2015).[27] M. M. Valado, C. Simonelli, M. D. Hoogerland,I. Lesanovsky, J. P. Garrahan, E. Arimondo,D. Ciampini, and O. Morsch, Phys. Rev. A ,040701 (2016).[28] M. van Horssen, E. Levi, and J. P. Garrahan, Phys. Rev.B , 100305 (2015).[29] A. Smith, J. Knolle, D. L. Kovrizhin, and R. Moessner,Phys. Rev. Lett. , 266601 (2017).[30] Z. Lan, M. van Horssen, S. Powell, and J. P. Garrahan,Phys. Rev. Lett. , 040603 (2018). [31] C. Turner, A. Michailidis, D. Abanin, M. Serbyn, andZ. Papi´c, Nature Phys. , 745 (2018).[32] H. Touchette, Phys. Rep. , 1 (2009).[33] J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. vanDuijvendijk, and F. van Wijland, Phys. Rev. Lett. ,195702 (2007).[34] J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. vanDuijvendijk, and F. van Wijland, J. Phys. A , 075007(2009).[35] V. Lecomte, C. Appert-Rolland, and F. van Wijland, J.Stat. Phys. , 51 (2007).[36] C. Appert-Rolland, B. Derrida, V. Lecomte, and F. vanWijland, Phys. Rev. E , 021122 (2008).[37] L. O. Hedges, R. L. Jack, J. P. Garrahan, and D. Chan-dler, Science , 1309 (2009).[38] T. Speck, A. Malins, and C. P. Royall, Phys. Rev. Lett. , 195703 (2012).[39] J. K. Weber, R. L. Jack, and V. S. Pande, J. Am. Chem.Soc. , 5501 (2013).[40] C. P. Espigares, P. L. Garrido, and P. I. Hurtado, Phys.Rev. E , 032115 (2013).[41] R. L. Jack, I. R. Thompson, and P. Sollich, Phys. Rev.Lett. , 060601 (2015).[42] D. Karevski and G. M. Sch¨utz, Phys. Rev. Lett. ,030601 (2017).[43] Y. Baek, Y. Kafri, and V. Lecomte, Phys. Rev. Lett. , 030604 (2017).[44] T. Oakes, S. Powell, C. Castelnovo, A. Lamacraft, andJ. P. Garrahan, Phys. Rev. B , 064302 (2018).[45] C. Giardina, J. Kurchan, and L. Peliti, Phys. Rev. Lett. , 120603 (2006).[46] F. C´erou and A. Guyader, Stoch. Anal. Appl. , 417(2007).[47] V. Lecomte and J. Tailleur, J. Stat. Mech. , P03004(2007).[48] T. Nemoto, F. Bouchet, R. L. Jack, and V. Lecomte,Phys. Rev. E , 062123 (2016).[49] U. Ray, G. K.-L. Chan, and D. T. Limmer, J. Chem.Phys. , 124120 (2018).[50] K. Klymko, P. L. Geissler, J. P. Garrahan, and S. White-lam, Phys. Rev. E , 032123 (2018).[51] G. Ferr´e and H. Touchette, J. Stat. Phys. , 1525(2018).[52] B. Derrida and J. L. Lebowitz, Phys. Rev. Lett. , 209(1998).[53] J. de Gier and F. H. L. Essler, Phys. Rev. Lett. ,010602 (2011).[54] A. Lazarescu and K. Mallick, Journal of Physics A: Math-ematical and Theoretical , 315001 (2011).[55] M. Gorissen, A. Lazarescu, K. Mallick, and C. Van-derzande, Phys. Rev. Lett. , 170601 (2012).[56] N. Cramp´e, E. Ragoucy, V. Rittenberg, and M. Vanicat,Phys. Rev. E , 032102 (2016).[57] A. Lapolla and A. Godec, New J. Phys. , 113021(2018).[58] M. Gorissen, J. Hooyberghs, and C. Vanderzande, Phys.Rev. E , 020101 (2009).[59] P. Helms, U. Ray, and G. K.-L. Chan, arxiv:1904.07336(2019).[60] M. Baiesi, C. Maes, and B. Wynants, Phys. Rev. Lett. , 010602 (2009).[61] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I.Cirac, Quantum Inf. Comput. , 401 (2007).[62] In the case of open boundary conditions, as used in this work, the first and last tensors reduce to rank-2 tensorsof dimensions d × D .[63] M. A. Nielsen and I. L. Chuang, Quantum Computationand Quantum Information: 10th Anniversary Edition ,10th ed. (Cambridge University Press, New York, NY,USA, 2011).[64] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. , 277 (2010).[65] Strictly speaking, the statement holds for states whichfulfill an area law in Renyi entropies S α = log(tr ρ α ) / (1 − α ) with 0 < α < , 2863 (1992).[67] G. Vidal, Phys. Rev. Lett. , 147902 (2003).[68] F. Verstraete, D. Porras, and J. I. Cirac, Phys. Rev.Lett. , 227205 (2004).[69] I. P. McCulloch, J. Stat. Mech. , P10014 (2007).[70] F. Verstraete, V. Murg, and J. Cirac, Adv. Phys. ,143 (2008).[71] U. Schollw¨ock, Ann. Phys. , 96 (2011).[72] Notice that it is also possible to define MPS directly inthe thermodynamic limit, and optimize them numericallywith appropriate methods [70, 71].[73] For details on the MPS numerics, their convergence, andfor the comprehensive set of results for both the FA andEast models, see Supplemental Material.[74] T. Bodineau, V. Lecomte, and C. Toninelli, J. Stat.Phys. , 1 (2012).[75] T. Bodineau and C. Toninelli, Commun. Math. Phys. , 357 (2012).[76] T. Nemoto, R. L. Jack, and V. Lecomte, Phys. Rev.Lett. , 115702 (2017).[77] Supplemental Material.[78] R. L. Jack and P. Sollich, J. Phys. A , 015003 (2013).[79] The GS in the limit s → −∞ seems to be gapped andwith low entanglement ( D ∼
10 provides a very goodapproximation [73]). The GS energy of the FA in thislimit is almost exactly twice the one for East, and theoverlap of their states is very high, suggesting they havesimilar GS, or rather the FA one is the superposition ofthat of the East and the reflected “West” model.[80] C. Castelnovo, C. Chamon, and D. Sherrington, Phys.Rev. B , 184303 (2010).[81] B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, NewJournal of Physics , 025012 (2010), 0804.3976.[82] See e.g. [71] for more details on initialization, convergencecriteria, etc.[83] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac,Phys. Rev. Lett. , 030504 (2008). SUPPLEMENTAL MATERIALNUMERICAL METHOD
The main MPS algorithm employed for this work is the variational optimization of a MPS with open boundaryconditions, in order to solve the minimization | Ψ (cid:105) = argmin (cid:104) Ψ | H | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) , (S1)over the set of MPS with fixed bond dimension D . The solution is the MPS approximation to the ground stateof the Hamiltonian H . There are many reviews in the literature describing the development, technical details andapplications of tensor network algorithms like this one, as well as their extensions to infinite systems, finite temperatureand dynamics, and possible extensions to higher dimensions [70, 71].It is convenient to express the algorithm fully in terms of tensor networks, by writing the Hamiltonian as a matrixproduct operator (MPO) [69, 81], i.e. a MPS vector in the tensor product basis of operators (that is, as a linearcombination of products of Pauli matrices). Local Hamiltonians as the ones considered in this work have an exactMPO expression with small constant bond dimension D H that does not depend on the system size. Evaluating itsexpectation value in a MPS of bond dimension D , which is the fundamental ingredient for the variational minimizationof the energy, has then a cost that scales as O ( dD H D ) in terms of the tensor dimensions, and linearly with the systemsize. This is crucial for the efficiency of the variational algorithm. In particular, we can write the East Hamiltonianmodel for open boundary conditions with bond dimension D H = 3 (or 4 for periodic chains) and the FA Hamiltonianwith D H = 4 (or 6 for periodic boundary conditions). (a) Pictorial representation of aMPS. (b) A local Hamiltonian has anexact MPO form.(c) The norm as contraction ofthe state with its adjoint overthe physical indices. (d) The energy is computed as theexpectation value of an MPO in aMPS state.(e) The effective norm at onesite N eff is obtained leaving outfrom the norm the tensorscorresponding to the site. (f) Effective Hamiltonian H eff asa TN. FIG. S1: Tensor networks and their contractions can be represented in a convenient pictorial language, whichsimplifies the description of algorithms and operations. A solid geometrical form (e.g. circles or squares above)represents a tensor, with as many indices as depicted legs. A contraction of two tensors over a certain index isrepresented as a connecting line. The pictures show the graphical representation of MPS, MPO and theircontractions, as they appear in the variational algorithm.The variational optimization then proceeds by fixing all tensors of the ansatz | Ψ (cid:105) = d (cid:88) i ,...i N =1 tr (cid:0) A i A i . . . A i N N (cid:1) | i i . . . i N (cid:105) (S2)but the one for site k , A k , and rewriting the optimization (S1) as a local problem in terms of the single variable tensor.The local problem boils down to a generalized eigenvalue problem for the vectorized tensor, H eff A k = λ min N eff A k [70, 71], where H eff ( N eff ) is an effective Hamiltonian (norm matrix) of dimension dD × dD , obtained by contractingall tensors in (cid:104) Ψ | H | Ψ (cid:105) and in (cid:104) Ψ | Ψ (cid:105) , except for A k ; see a pictorial representation in Fig. S1. This problem can thenbe solved with a standard eigensolver from a linear algebra numerical package, and the minimum eigenvalue λ min corresponds to the estimate of the ground state energy. Using a sparse eigensolver allows to keep the cost scaling as D and to deal with very large values of the bond dimension. The k -th tensor is updated with the solution of thelocal optimization, and then the procedure is repeated for all the tensors in the chain, sweeping back and forth until acertain convergence criterion (typically on the energy) is met. The same algorithm can be used to find higher excitedstates by imposing that the solution is orthogonal to already found levels. This can be imposed at the level of thelocal problem, without changing the scaling of the leading cost, which is always O ( D ) (and grows polynomially withthe number of computed levels).For a run with fixed bond dimension, the algorithm is guaranteed to converge, because it can only decrease theenergy in every step, although it may do so to a local minimum. To improve the precision, one increases the bonddimension of the ansatz, typically using the previous solution with smaller D as initial guess. In a typical application,the algorithm is repeatedly run with increasing bond dimension, until the energy of the state is converged to thedesired precision [82]. A notorious case in which the algorithm is slow to converge is that of critical systems, wherethe ground state requires a bond dimension that grows polynomially with the system size in order to achieve a fixedprecision, a situation that is well understood by DMRG practitioners. But on the other hand, having a state thatcan be well approximated by a MPS does not guarantee convergence of the algorithm. A large density of states alsohinders convergence, as happens for instance when trying to approximate excited states in the middle of the spectrum. Convergence -5 (a) East s = − − -6 (b) East s = 3 · − -5 (c) FA s = − − -5 (d) FA s = 3 · − FIG. S2: Convergence of the energy as a function of the bond dimension in some of the most difficult cases ( c = 0 . s ). -10 -10 -2 -10 -4 -10 -5 (a) c = 0 . s < -4 -10 -5 (b) c = 0 . s > -10 -10 -2 -10 -4 -10 -6 -10 -5 (c) c = 0 . s < -8 -4 -10 -5 (d) c = 0 . s > FIG. S3: Energy standard deviation (square root of variance) ∆ H = (cid:112) (cid:104) H (cid:105) − (cid:104) H (cid:105) in the MPS approximation tothe ground state for the East model. The plots show, for system sizes N = 20 (blue) and 400 (red) the systematiclowering of the variance as the bond dimension is varied from D = 2 (squares), to 10 (diamonds) and 20 (triangles).We show the detail of the most difficult regions, namely the region of the plateaus for small s < s > D = 100, the variance of the peaks is reduced to ∆ H (cid:46) − .We find the ground states over the largest part of the parameter space to be very well approximated by MPS withsmall bond dimension. The quality of the MPS approximation can be gauged from the convergence of observables -10 -10 -1 -10 -2 -10 -3 -10 -5 (a) c = 0 . s < D = 20 -5 -3 -1 -10 -5 (b) c = 0 . s > D = 20 -10 -10 -2 -10 -4 -10 -5 (c) c = 0 . s < D = 20 -6 -4 -2 -10 -5 (d) c = 0 . s > D = 20 FIG. S4: Energy standard deviation (square root of variance) ∆ H = (cid:112) (cid:104) H (cid:105) − (cid:104) H (cid:105) in the MPS approximation tothe ground state for the FA model. Qualitatively, we observe similar effects as for the East model, described infigure S3, As in figure S3, we show the systematic lowering of the variance as the bond dimension is varied from D = 2 (squares), to 10 (diamonds) and 20 (triangles) for system sizes N = 20 (blue) and 400 (red). Again, bonddimension D <
100 is enough to ensure very small variance over the most challenging range of parameters.as the bond dimension is increased. We find this to be in general very fast, even for system sizes of several hundredsites. We let the algorithm use bond dimensions as large as D = 100, but in most of the cases analyzed, we find thata bond dimension D = 20 is enough for the energy to be sufficiently converged. As illustrated in figure S2, only in afew cases, mostly for small values of c and around the phase transition, we find that a larger bond dimension allowsus to reach a lower energy. We also find that convergence becomes difficult for large systems when we try to explorethe region of the phase transition at small positive s in both models. Since we have demonstrated that the statesdo not develop a large entropy, even in this region, we attribute this behaviour to the density of states at the lowestenergy becoming larger for increasing system size.A more accurate measure of how close the approximation is to an actual eigenstate is however provided by the energyvariance, ∆ H = (cid:104) H (cid:105) − (cid:104) H (cid:105) , which can be computed efficiently for any MPS. In the cases studied in the paper, wefind that a very small bond dimension, D = 20, is already enough to obtain a very small variance ∆ H (cid:46) O (10 − )for a wide range of values of s and all system sizes up to N = 400 (see the upper row of figures S3 and S4). Theexceptions are the region of the phase transition at small s > c decreases (see figuresS3b, S3d, S4b and S4b). DETAILED NUMERICAL RESULTSFinite size scaling of the active-inactive phase transition -6 -4 -2 -10 -3 -10 -4 -10 -5 -10 -6 -10 -7 (a) (Minus) energy density as afunction of s -6 -5 -4 (b) Activity as a function of s -6 -5 (c) Susceptibility as a function of s FIG. S5: Scaling of the phase transition location for the East model with c = 0 . N ∈ [20 , c ∈ [0 . , . c = 0 . -4 -2 -10 -2 -10 -3 -10 -4 -10 -5 -10 -6 (a) (Minus) energy density as afunction of s -5 -4 -3 (b) Activity as a function of s -4 (c) Susceptibility as a function of s FIG. S6: Scaling of the phase transition location for the FA model with c = 0 . -4 -2 -10 -2 -10 -3 -10 -4 -10 -5 -10 -6 (a) (Minus) energy density as afunction of s -5 -4 -3 (b) Activity as a function of s -4 -3 (c) Susceptibility as a function of s FIG. S7: Scaling of the phase transition location for the FA model with c = 0 . s , in agreement with the perturbative calculation around s = 0. The grey dashed line in the figures shows the linear response prediction in the thermodynamic limit, namely θ/N = − c (1 − c ) s for the East and θ/N = − c (1 − c ) s for the FA model. The intersection of this line with theasymptotic value for s → ∞ , shown as dashed coloured lines in the plots, scales as 1 /N . However, in the figures it isevident that the value of s c at which the actual crossing occurs can be orders of magnitude away from this prediction.The activity k = − θ (cid:48) ( s ) N = 1 N dE GS ( s ) ds (S3)can also be directly evaluated from the MPS ansatz, since dE GS ds = (cid:104) Ψ | dHds | Ψ (cid:105) , (S4)only requires computing the expectation value of local and two-body observables. The results are shown in figuresS5b, S6b and S7b. The numerical derivative of the activity yields the susceptibility χ = dθ ds , shown in Figs. S5c, S6cand S7c. The location s c of the phase transition for each system size is most precisely determined from the positionof the peak in χ .An alternative FSS analysis of the transition can be made by following the approach of Bodineau, Lecomte andToninelli in Ref. [75] (hereafter BLT) which considers in detail the FA model. BLT use the fact that within a regionof size 1 /N around the transition, that is for s of order 1 /N or equivalently for λ = sN with λ = O (1), the SCGF interms of λ , φ ( λ ) = θ ( λ/N ) is of O (1) [rather than of O ( N ) as for s finite] and should progressively interpolate betweentwo behaviours at large N (see also [74]). The two behaviours are that of linear response, φ = −(cid:104) k (cid:105) λ , for λ < λ c , anda regime of constant φ = − Σ for λ ≥ λ c . Here Σ is a “surface tension” related to the creation of an interface betweenthe active and inactive phases, while (cid:104) k (cid:105) is the equilibrium activity per unit time. Figure S8 presents the SCGF in1this representation for both models, cf. Fig. 2 of Ref. [75]. The crossover between these two regimes is apparent bothin the FA, as described by BLT, and also in the East model. Note that accessing the constant regime on the right ismore difficult for lower c .The prediction of BLT is that φ ( λ ) + Σ should behave as − Aλ ν for λ > λ c . Figure S9 shows such a scaling for boththe FA model and the East model, cf. Fig. 3 of [75]. We have estimated the exponent ν and the constant Σ in thefollowing way. Since the activity is (minus) the first derivative of the SCGF, in the relevant region it should scale as λ ν − [cf. the scaling of the susceptibility in Fig. 1(b) of the main text]. This allows to obtain ν without the need tosimultaneously fit Σ. As discussed above, the activity can be calculated efficiently as it corresponds to the contractionof an MPO with the MPS. Figure S10 shows the activity for both models. The smaller c is the larger the system size N is required to be to accurately extract ν . With the exponent ν in hand we can then estimate Σ by subtracting the λ dependence from φ . This is illustrated in Fig. S11(a,b) for c = 0 . ≈ .
077 for the surface tension at c = 0 .
5. In our case, for the FA model wefind Σ ≈ . / c = 0 .
5, see Fig. S11(a), the factor of a half coming from the fact that we use open boundaryconditions (which allows a single interface to be created, in contrast to the periodic boundary condition case). Thisresult seems to confirm the BLT prediction. For the East model we find a slightly lower value at c = 0 .
5, Σ ≈ . c , with Σ seemingly going asΣ ∝ c ξ with ξ ≈ . c for the East model). For the exponent,BLT predicted ν = 2 /
3. While for large c this is compatible with our findings, see Fig. S10, we seem to find that ν increases with decreasing c . Nevertheless, this discrepancy might be due to the fact that for smaller values of c extracting both Σ and ν gets progressively more difficult. (a) FA c = 0 . -3 (b) FA c = 0 . (c) East c = 0 . -4 (d) East c = 0 . FIG. S8: SCGF φ ( λ ) with λ = sN in transition region, cf. Fig. 2 of Ref. [75], for various sizes. The dashed line isthe linear response behaviour. Structure of active phase
As discussed in the main text, the active phase of both models exhibits very different features, which we canexplore with our results. In the East model, for small values of c , we recover the hierarchy of plateaus with welldefined average density predicted in [78], extending between values of ν = 1 − e s equal to integer powers of c . Thisis clearly appreciated in Fig. S12a for c = 0 .
02, where the plateau at (cid:104) n (cid:105) = 1 / ν = 0 .
02 is alreadyconverged in system size. For c = 0 . ν ∼ .
14, as shown in Fig. S12b, while for yet larger values of c , the curve is smooth (see for instance the insetof figure 2(a) in the main text). For the FA model, on the other hand, there are no similar features in the averagedensity, as shown explicitly by Fig. S12c and figure 2(b) in the main text.To explore in more detail the structure of the active phase in both models, we have computed the spatial dependenceof the density across the same range of values of s spanned by Fig. S12. We show the results for a system size N = 20in figures S13 (for the East) and S14 (for the FA model). In the case of the East model, the figure shows how the fixedaverage density of the plateaus is achieved by means of a regular modulation of the local density. For the (cid:104) n (cid:105) s = 1 / s as c increases,and they have disappeared completely at c = 0 . c , as explicitly shown by Fig. S14 for c ∈ [0 . , . (a) FA c = 0 . (b) FA c = 0 . (c) FA c = 0 . (d) FA c = 0 . (e) East c = 0 . (f) East c = 0 . (g) East c = 0 . -3 (h) East c = 0 . FIG. S9: Collapse of SCGF φ ( λ ) for λ > λ c , cf. Fig. 3 of Ref. [75]. -4 -2 (a) FA c = 0 . -4 -2 (b) FA c = 0 . -4 -2 (c) East c = 0 . -5 -1 (d) East c = 0 . FIG. S10: Activity k ( s ). From the region where the curves collapse we can extract the exponent ν used in Fig. S9.We only show two values of c for comparison, but this procedure was used to extract ν for all other c .The period three modulation of density shown in Fig. S13 for the case of N = 20 is present also in larger systems,but the detailed structure depends on the congruence modulo three of the system size. This is shown explicitly inFig. S15, for the particular case c = 0 . s = − . Entanglement
The states we are looking for contain very little entanglement. We have shown explicitly in the main text that theentanglement entropy with respect to a division of the chain in two S E remains bounded by a constant even in theregion of the phase transition for the East model [see Fig. 3 in the main text]. The same is true in the case of the FAmodel, as shown in Fig. S16a for a chain of length N = 200. As one can expect from the previous discussion, in thiscase, no peaks of the entropy occur within the active phase, and the entropy is smooth for all values of c ; see zoom inFig. S16b, to be compared to Fig. 3(b) in the main text. The maximum of the entropy occurs, instead, around thephase transition, for small positive values of s , but the magnitude of the peak shows only a mild dependence on thesystem size, similar to what we observed in the East model.There are however differences between both models in the behaviour of the entropy at very small s , as shown inFig. S17. In the East model, the entropy at s = 0 is strictly zero, corresponding to a product ground state [see Eq. (4)of the main text and Fig. S17a], and builds up to a peak around the transition. In the case of the FA model, the3 -1 -0.038-0.036-0.034-0.032 (a) FA c = 0 . -0.032-0.031-0.03-0.029 (b) East c = 0 . -2 -10 -4 -10 -6 (c) East c = 0 . FIG. S11: (a) SCGF after subtracting the term proportional to λ ν in order to estimate Σ for the FA model. (b)Same for East model. (c) Surface tension Σ as a function of c for both models extracted by the procedure of panels(b,c) for all available c . For the FA model we get Σ ∝ c ξ with ξ ≈ . -10 -10 -1 -10 -2 -10 -3 (a) East model, c = 0 . -10 -10 -1 -10 -2 (b) East model, c = 0 . -10 -10 -1 -10 -2 -10 -3 -10 -4 (c) FA model FIG. S12: Average occupation in the active phase as a function of − ν = e s −
1. For the East model (two leftmostplots), plateaus appear for small values of the equilibrium concentration, up to c = 0 . ν ∼ .
14 ( s ∼ − . c , the curve is smooth (see the inset of Fig. 2(a) in themain text), the same as for the FA model over the whole range of values of c .ground state at s = 0, in the subspace orthogonal to the state with no excitations, is not a product state [see Eq. (3)of the main text], and thus the entropy is not exactly zero, but has a finite value at s = 0, namely S FAE ( s = 0) = H (cid:115) − (cid:18)
21 + (1 − c ) − N/ (cid:19) , (S5)where H( x ) = − x log x − (1 − x ) log (1 − x ) is the Shannon entropy. The value of S FAE ( s = 0) decreases fast withthe system size N (Fig. S17b). For s > S E grows also in the case of the FA model, towards a value of order one.Afterwards we find an almost vanishing energy. Notice that to the right of the transition, in the limit s → ∞ , theground state is doubly degenerate, and the component that is symmetric under parity, thus in the same subspace asthe ground state at s = 0, would have entanglement S E = 1. However, at sufficiently large s , the algorithm prefers tobreak the symmetry to find a solution with smaller bond dimension.The bounded entropy alone is not enough to guarantee the approximability of the ground state by a MPS [83].To gather more compelling evidence we can also study the scaling of Renyi entropies, defined as S α = Tr ρ α / (1 − α )for α >
0, and which in the limit α → S α with α < S α can be computed efficiently. We show the values of thevon Neumann and Renyi entropies for the active phase of the East model with c = 0 .
05 in the region of the plateaus(since it is the region with the largest entropy we found), in Fig. S18. Again we see that, although the magnitude ofthe peaks is not fully converged, the dependence on the system size is very mild.4 -10 -1 -10 -2 -10 -3 (a) c = 0 . -10 -1 -10 -2 -10 -3 (b) c = 0 . -10 -1 -10 -2 -10 -3 (c) c = 0 . -10 -1 -10 -2 -10 -3 (d) c = 0 . FIG. S13: Spatial distribution of density n i = (1 − (cid:104) σ zi (cid:105) ) / s for a chain of length N = 20 and increasing value of c . -10 -1 -10 -2 -10 -3 (a) c = 0 . -10 -1 -10 -2 -10 -3 (b) c = 0 . -10 -1 -10 -2 -10 -3 (c) c = 0 . -10 -1 -10 -2 -10 -3 (d) c = 0 . FIG. S14: Spatial distribution of density n i = (1 − (cid:104) σ zi (cid:105) ) / s fora chain of length N = 20 and increasing value of c . (a) N = 18 (b) N = 39 (c) N = 90 (d) N = 19 (e) N = 40 (f) N = 100 (g) N = 20 (h) N = 38 (i) N = 80 FIG. S15: Distribution of density along the chain for c = 0 .
02 and s = − . (cid:104) n (cid:105) s = 1 / -2 -1 0 100.511.5 (a) N = 200 -10 -10 -1 -10 -2 -10 -3 -10 -4 (b) Active phase of N = 200 FIG. S16: Entropy in the ground state of the FA model. The left plot shows the overall behaviour with respect to s for a chain of size N = 200 and different values of c , and the right one shows the detail of the active phase. -5 (a) East model c = 0 . -1 -0.5 0 0.5 110 -3 (b) FA model c = 0 . FIG. S17: Entropy around s = 0 in the ground state of the East and FA model at c = 0 .
05 for different system sizes. -10 -10 -1 -10 -2 -10 -3 (a) S α , α = 0 . -10 -10 -1 -10 -2 -10 -3 (b) S α → = S E -10 -10 -1 -10 -2 -10 -3 (c) S α , α = 2 FIG. S18: Entanglement in the active phase of the East model, as measured by the von Neumann and Renyientropies S α of the half-chain for c = 0 ..