Boundary conditions dependence of the phase transition in the quantum Newman-Moore model

We study the triangular plaquette model (TPM, also known as the Newman-Moore model) in the presence of a transverse magnetic field on a lattice with periodic boundaries in both spatial dimensions. We consider specifically the approach to the ground state phase transition of this quantum TPM (QTPM, or quantum Newman-Moore model) as a function of the system size and type of boundary conditions. Using cellular automata methods, we obtain a full characterization of the minimum energy configurations of the TPM for arbitrary tori sizes. For the QTPM, we use these cycle patterns to obtain the symmetries of the model, which we argue determine its quantum phase transition: we find it to be a first-order phase transition, with the addition of spontaneous symmetry breaking for system sizes which have degenerate classical ground states. For sizes accessible to numerics, we also find that this classification is consistent with exact diagonalization, Matrix Product States and Quantum Monte Carlo simulations.


I. INTRODUCTION
In this paper, we study the ground state phase transition of the quantum Newman-Moore model, or quantum triangular plaquette model.The classical triangular plaquette model (TPM), introduced by Newman and Moore [1], is a model of Ising spins interacting in triplets in (half of) the plaquettes of a triangular lattice.Despite the absence of quenched disorder and its trivial static properties, the model has rich glassy dynamics [1][2][3].The TPM is an important model as it realises in an interacting system the paradigm of slow (super-Arrhenius) relaxation at low temperatures due to effective kinetic constraints.This phenomenon is central to the dynamic facilitation picture of the glass transition [4][5][6].The physics of the TPM can also be generalised to three dimensions, for example in the five-spin interaction square-pyramid model [7], or maintaining the triangular interactions in the models of Ref. [8].A three-dimensional generalisation of the TPM with non-commuting terms [9] actually started what is now the field of fractons [10][11][12][13].
The simplest way to transform the TPM into a quantum model is by adding a transverse field term to the classical Hamiltonian.Such quantum TPM (QTPM, or quantum Newman-Moore model) was considered in the context of fractons in Refs.[14,15].Numerics in Ref. [14] suggested that the ground state undergoes a first-order transition.A related work studying the large deviations of plaquette observables in the stochastic dynamics of independent spins [16] also found numerical evidence for a first-order transition at the self-dual point of the model.
In contrast, the results from Ref. [17] indicated that the transition is continuous, with a particular form of fractal symmetry breaking.The classical TPM and its connection with fractals and topological order was also drawn in * ksfairopoulos@gmail.com Ref. [18].Here we aim to resolve these discrepancies by exploiting a general connection between D-dimensional cellular automata (CA) [19] and the ground states of (D + 1)-dimensional classical spin models [10,20].By using this method in the specific case of the QTPM with periodic boundary conditions, we are able to characterize the approach to its quantum phase transition in the large size limit.Our key observation is that the nature of the transition depends on the specific lattice dimensions, and this is manifested in the finite size scaling.
For the TPM, the relevant CA is Rule 60 [19] and not Rule 90 that might be assumed from comment [18] in Ref. [1].For system sizes where one dimension is a power of two, Rule 60 has a single fixed point [21], implying a single energy minimum for the classical TPM.In such cases, we verify that the quantum phase transition in the QTPM is first-order (that is, a sequence of such system sizes tends to a first-order transition in the large size limit).This also holds for other sizes for which Rule 60 has no non-trivial attractors.However, for certain sizes there can be periodic orbits on top of the fixed point for the CA, giving rise to classical ground state degeneracies in the TPM.For the quantum model, this translates into a mixed order quantum phase transition.We provide evidence for this scenario by means of numerical simulations, namely, for small sizes using exact diagonalisation, and for large sizes using both Matrix Product State approximations of the ground state, and continuous-time Quantum Monte Carlo [22][23][24][25].
The paper is organised as follows.In Sec.II, we review the classical and quantum TPM.In Sec.III, we provide the necessary background on CA and the connection to the ground states of the classical TPM.In Sec.IV, we discuss the ground state phase transition of the QTPM in terms of the symmetries that follow from the properties of the associated CA, and support our predictions with numerical simulations.In Sec.V we give our conclusions.In the Appendix we provide further details, including the case of the QTPM with open boundaries.
The triangular plaquette model (or TPM or Newman-Moore model) [1][2][3] is a model of Ising spins s i = ±1 on the sites i = 1, . . ., N of a triangular lattice, with cubic interactions between the spins on the corners of downward-pointing triangles of the lattice, Fig. 1.The Hamiltonian of this classical model reads In what follows, it will be convenient to consider the equivalent model on a square lattice of size N = L × M , with classical Hamiltonian, where we assume periodic boundary conditions in both directions by identifying x + L = x mod L y + L = y mod L.
In Eq.( 2), we label the spins by s x,y at site with coordinates (x, y) in the square lattice.
The classical TPM has been predominantly studied in the context of the glass transition.For lattices with at least one dimension being a power of two, the energy Eq.( 1) reduces to that of the non-interacting plaquette variables [1][2][3], where d = s i s j s k with i, j, k ∈ for every downwardpointing triangle in the lattice.When at least one dimension is a power of two, the relation between plaquettes and spin variables is one-to-one, exactly proving the above.The thermodynamics of the TPM is therefore one of free binary excitations and, as such, it is essentially trivial.
In contrast to the statics, the single spin-flip dynamics of the TPM is highly non-trivial, as flipping one spin changes three adjacent plaquettes.This implies that at low temperatures, where excited plaquettes are suppressed, cf.Eq.( 3), the dynamics has effective kinetic constraints [2].These dynamical constraints lead to an activated relaxation similar to that of the East model [26], with relaxation times growing as the exponential of the inverse temperature squared [2] (a super-Arrhenius form known as the "parabolic law" [27]).Similar glassy behaviour is seen in generalisations of the TPM with odd plaquette interactions [8,28].The TPM has also been considered in the presence of a (longitudinal) magnetic               field [29] and in the related case of coupled replicas [30], and the TPM with open boundary conditions was studied in Ref. [31] through partial trace methods.The classical TPM was studied in the context of topological order in fractal models in Ref. [18].Autoregressive neural networks were applied to the classical TPM with limited success in [32].
B. Quantum: TPM in a transverse field Taking Eq.( 1) and adding a transverse field, we obtain the Hamiltonian of the quantum TPM (QTPM), where Z i and X i represent Pauli operators acting nontrivially on site i.This quantum model was studied in Refs.[14,15] and its connections to models of fractons were investigated.
The Hamiltonian ( 4) is expected to have a quantum phase transition at the self-dual point J = h [16].Numerical results from [14,16] suggested the transition to be first-order.Ref. [16] used trajectory sampling in systems with linear size a power of two and periodic boundaries.In contrast, Ref. [17] found evidence for a continuous phase transition with fractal symmetry breaking using stochastic series expansion methods, also with periodic boundary conditions but not restricted to power of two sizes.The study of QTPM was connected to Rydberg atoms in Ref. [33].Further studies on the QTPM and its generalisations from the viewpoint of fracton field theory were presented in Refs.[34,35].
Cellular automata (CA) consist of a D-dimensional array of sites evolving under discrete-time and synchronous dynamics given by a CA rule.Under this evolution, the state is updated by a deterministic rule which is local in time (although generalisations exist) [19,[36][37][38].As a result, CA dynamics gives rise to, often rich, (D + 1)dimensional structures.In what follows, we consider D = 1, linear CA with deterministic local transition rules with sites taking values from the finite field F 2 .
Discrete cellular automata are fully specified by an initial configuration of L sites and a local transition rule.The local update rule determines the configuration of every site at each timestep using the local neighbourhood of size r.For D = 1, the simplest CA rules are those defined only by r = 1.For a neighbourhood of three sites, r = 1, there are 8 possible configurations giving rise to 256 possible choices for update rules.These elementary CA were classified by Wolfram [19,36].Figure 2(a) shows Rule 60, which will be the relevant one for the TPM: if we identify the empty and occupied sites of the CA with the up and down spins of the TPM, then Rule 60 is the same as the condition that the product of spins in Eq.( 2) is one, thus maximising the local energy.Figure 2(b) also shows the closely related Rule 90.CA like these are called elementary [36].
Figure 3(a,b) shows the patterns generated for Rules 60 and 90 starting from an initial single seed.Note that these depend on the boundary conditions (for a generic analysis on cellular automata with periodic boundaries, see Ref. [39]).For example, Fig. 3(b) shows Rule 90 for L = 64 and periodic boundaries: the timestep after the last one shown will take the CA to the trivial, empty configuration; in contrast for a length L = 63 the Sierpinski fractal [40] will continue.Similarly, time evolution will bring in one timestep the pattern of Fig. 3(a) to the trivial one for periodic boundaries with L = 32, but the fractal shape would continue to be reproduced for open boundaries.Focusing on Rules 60 and 90 for concreteness, CA evolution can be written as polynomials (or generating functions, see Ref. [36]), as with the evolution of the whole lattice obeying [36] Given a configuration at timestep t for the evolution of any given rule, the configuration at timestep t + 1 will be its successor and the one at t − 1 its predecessor.For Rules 60 and 90, following Ref.[36], we have: • There are no predecessors for configurations with an odd number of sites being equal to 1 for both  rules.For Rule 90, 2 L−1 configurations for a system of size L cannot be reached if L is odd and 3 × 2 L−2 if L is even.For Rule 60, 2 L−1 configurations cannot be reached for any system size.This means that rows with an odd number of ones can only occur as initial states for any given time evolution for Rule 60.
• Similarly, there may exist configurations that can be reached in one timestep by more than one predecessors.More specifically, given a configuration with at least one predecessor, there are 2 predecessors if the length L is odd and 4 if it is even for Rule 90, while there are always 2 for Rule 60.
Regarding the length of cycles for Rule 90, more results can be found in Ref. [36].For Rule 60, which is most relevant for the TPM, we discuss the cycle lengths and their multiplicities next.

B. Periodic orbits of Rule 60
The iterated map for D n : Z n → Z n and x = (x 1 , x 2 , ..., x L ) is known as the Ducci map (also known as rule 102 in Wolfram's notation [19,36], the mirror image of Rule 60, e.g.Ref. [41]).Ducci's map (and therefore Rule 60) exhibits periodic orbits or evolves to a fixed point depending on L being a power of two or not.For L = 2 k , any initial state evolves to the trivial state of zeros [42], while for L = 2 k it evolves to a richer attractor structure [43].
In F 2 and by using periodic boundary conditions, the Ducci map can be brought into matrix form as follows, The matrix D can be expressed as with S L the left shift map [21,43].It is easy to check that for n = 2 k , with k ∈ Z and, thus, a system that has L a power of 2 ends up in the trivial configuration.For Rule 60 the corresponding matrix is D .In order to obtain the cycles of Rule 60, we need the following definitions and results for Rule 102 from Refs.[21,37,38,[43][44][45][46]: • For any given CA, each array of sites can be thought of as a vector, v.For any such vector v in F 2 , its order is defined through the monic polynomial, which satisfies µ v (D)v = 0.
• The order of the minimal annihilating polynomial, ord µ v (λ), is equal to the smallest natural number n, such that • Assuming that µ v (λ) = λ k μv (λ) with k ≥ 0, then the k-th successor of v belongs to a cycle of length c = ord µ v .This applies to a vector v of any positive integer L and any linear map (or cellular automaton rule).
The minimal annihilating polynomial for the Ducci map was calculated in Refs.[21,45], based on the characteristic polynomial of the matrix D. It was found that We are now in a position to obtain the attractor structure of Rule 60 by following Ref.[38] (specifically "Principle C").We decompose the minimal annihilating polynomial into the product of its irreducible polynomials, π i (λ).We call their polynomial powers b i .Thus, and where r is the least common multiple of ord π i and t the smallest integer satisfying p t ≥ max (b 1 , b 2 , . . ., b m ). Figure 4 illustrates the different scenarios for cycles as a function of L. We show three different sizes, L = 5, 6, 8, small enough to be able to visualise the network of states.Configurations of the CA are identified by blue circles,  The second column indicates the cycle lengths and the third the total number of periods for Rule 60 (each period of length C is counted C-times), assuming that the least common multiple of the periods for each given L divides M , lcm C|M .The last column is constructed based on the number of predecessors for all configurations for Rule 60, given a length L.
and arrows indicate to which configurations they evolve to under the CA dynamics.Figure 4(a) shows L = 5: here there is one fixed point to which one other state evolves to (shown at the centre of the figure), and one limit cycle of period C = 15, to which all other configurations flow.Figure 4(b) shows a more general situation of multiple distinct cycles for the case of L = 6: there are two cycles of period C = 6, one cycle of period C = 3 and one fixed point.For the case of L a power of 2, as shown in Fig. 4(c) for L = 8, there is a unique fixed point and all states evolve towards it.
C. Classical ground states of the TPM From the analysis above for Rule 60, we can enumerate all the minimum energy configurations of the TPM, cf.Eq.( 2).The classical ground states for a system of size N = L × M are: (i) the state with all spins up, corresponding to the fixed point of Rule 60, for any value of L and M ; (ii) two-dimensional spin configurations that correspond to periodic trajectories of Rule 60, that is, CA trajectories starting from any of the states of a limit cycle, for all limit cycles whose period is contained an integer number of times in M -this occurs only for certain combinations of L and M (never if L is a power of two).
We show the relevant Rule 60 information in Table I for up to L = 40.Under the column C we give the distinct periods of the limit cycles.In the column labelled by M we give the corresponding degeneracy of classical ground states of the TPM, apart from the uniform spinup state, given that lcm C|M .For example, for L = 15, there is one cycle of length 3, three cycles of length 5, and 1091 cycles of length 15, which means 1 + M = 1+3+3×5+1091×15 = 16384 distinct two-dimensional spin configurations that minimise the energy in a TPM of size N = 15 × M as long as 15|M .In contrast, for L = 15, if 5|M but 3 M and 15 M , then there are 1 + 15 different ground states.The symmetries of a TPM of size N = L × M can be similarly constructed from the number of ground states and their multiplicities, as determined by Rule 60.

IV. PHASE TRANSITION IN THE QUANTUM TPM
As we will now show, the set of minimum energy configurations and the associated symmetries of the classical TPM, as obtained from the Rule 60 CA, determine the properties of the ground state phase transition in the quantum TPM, Eq.(4).
A. Symmetries of the QTPM Like its classical counterpart, the QTPM with periodic boundaries has full translational invariance.The symmetries of the QTPM will then be deduced by the results of Sec.III.A similar procedure has been followed by Ref. [10,15].
For systems with dimensions N = L × M , which can accommodate non-trivial cycles of Rule 60, the symmetries of the corresponding QTPM easily follow.Consider as a simple example the case of N = 3 × 3 with periodic boundaries in both dimensions.From Table I, we see that L = 3 has one cycle of period 3.This means that for M = 3 there are three non-trivial symmetries, given by the operators see Fig. 5(a).Note that translational invariance plays a crucial role in the determination of the symmetries: in the example above x , where T x is the translation operator in the x direction, and these symmetries alongside the identity form the Klein group K 4 , which, in turn, is isomorphic to This approach generalises to other system sizes, and the symmetries are products of K 4 .For example, the symmetries of the N = 5 × 15 and the N = 6 × 6 systems form the group K 4 ⊗ K 4 .Since the symmetry operators are products of the X Pauli matrices, cf.Eqs.(16)(17)(18), they commute with the transverse field term in Eq.( 4), and, therefore, are symmetries of the QTPM for all values of J and h.

B. The character of the quantum phase transition
The character of the phase transition for the QTPM is now easy to predict based on the information of its symmetries for a given system size.By this we mean: given a sequence of increasing system sizes with periodic boundaries, all having the same number of symmetries, the progressively sharper finite size crossovers will be indicative of an eventual phase transition in the large size limit whose character-first-order or continuous-will be determined by the underlying symmetries of the given system sizes in the sequence.As these symmetries in a system of size N = L × M depend on the precise values of both L and M , this analysis has to be done carefully.
In the next section we show numerical evidence for the general considerations we give here.
First consider the system sizes N = L × M such that the underlying Rule 60 has only one fixed point, as described previously.Then, the QTPM has no non-trivial symmetries.In the limit J h, there is a single ground state corresponding to the all-up state, and a vanishing number of excited plaquettes, as the classical TPM has a single minimum.In the opposite limit, J h, the ground state is paramagnetic, with spins aligned in the x direction, with an explicit Z 2 symmetry for its ground state, corresponding to a high density of excited plaquettes.In the limit of large N , we expect the quantum phase transition at the self-dual point J = h to be firstorder.
A second scenario is that of system sizes where the underlying Rule 60 cycles give rise to non-trivial symmetries in the QTPM.In this case, in the limit J h, the ground state is the same paramagnetic one as before, invariant under a global Z 2 symmetry.However, for J h there exists spontaneous breaking of the symmetries emerging from Rule 60.This case, therefore, has characteristics of a first-order phase transition (the explicit symmetry breaking of the global symmetry of the paramagnetic ground state) with those of a continuous transition (spontaneous symmetry breaking).This is similar to the first-order phase transitions observed in kinetically constrained models [47] and suggests that, for local observables, the phase transition will appear firstorder; for example, in a discontinuous jump in the excited plaquette density, Appropriate operators will quantify the symmetry breaking of the degenerate classical ground states.The choice of these operators depends on the size and specific lattice dimensions of the system.Consider, for example, the case of L = 3 and M = 3k with k ∈ N, where we know from Rule 60 that there is a single fixed point and the three non-trivial ground states.For the trivial ground state this operator is just the magnetisation M z = 1 N N i Z i .To detect the symmetry breaking into the other three states, we can define the three operators M m z = 1 N N i (−1) ni Z i , where n i = 1 if the spin i is flipped for a state of the cycle m of Rule 60, and n i = 0 otherwise.For example, for the state associated to Eq.( 16), we have Note that when there are multiple non-trivial ground states connected by translations, there will be no single operator taking the form of a sum of local terms for which it will be possible to discern between them.For example, in the N = 3k case, M z will only distinguish between the trivial ground state and the 3-fold degenerate ones, while the operators M m z will be able to distinguish between one of the non-trivial ground states.

C. Numerics
We now provide evidence for the general observations above from numerical simulations.For small systems we use exact diagonalization (ED) [48][49][50], allowing the study of system sizes up to 28 sites.For larger systems we estimate the properties of the ground states using two different approaches.The first of these is Matrix Product States (MPS) [51][52][53][54], which we "snake" around the 2D lattice, and optimize with the 2D Density-Matrix Renormalization Group (DMRG) [55,56].By employing a bond dimension up to D = 1000, we are able to reliably estimate the ground state properties for system sizes on the square lattice for up to N = 16 × 16.Time and memory constraints hinder progressively the convergence in the paramagnetic phase, where J ∼ h.As discussed below, we are also able to apply MPS to cylindrical systems, which can be considered to be quasi-1D, allowing us to reach much larger sizes than in the case of square geometries.
To confirm the results of 2D DMRG, we also employ Quantum Monte Carlo (QMC) methods.In particular, we use the continuous-time expansion (ctQMC) [22], with local spin updates which re-draw the entire trajectory of a single spin, subject to a time-dependent environment, where the trajectories of unmodified spins are considered to act as a "heat bath", e.g., see Refs.[23,24].We run our simulations with an inverse temperature of β = 128, which we find to be large enough to converge to the ground state [25].

First-order transition
As explained above, when the underlying Rule 60 has a single fixed point, and the classical TPM a single energy minimum with all spins up, we thus expect the phase transition of the QTPM to be first-order.19), in the QTPM as a function of J for fixed h = 1, with N = L × L with L a power of two.We compare results from MPS and ctQMC obtained here with results from Ref. [16].The numerical data indicates a first-order transition at J = h.
Figure 6 shows results for N = L × L with L = 4, 8, 16.We show that our MPS and ctQMC results coincide with the large deviations results of Ref. [16], obtained via transition path sampling (TPS).What is plotted is the plaquette density M zzz as a function of the coupling J for fixed transverse field h = 1.The data indicates a firstorder transition at J = h, as expected.For J > 1.0, we see small deviations in the TPS results, due to the extra field used for the acquisition of this data in Ref. [16].Similar issues are observed for ctQMC close to the J = 1.0 point, where the single spin updates do not allow for the collective effects necessary to move between phases.
In Fig. 7, we show the results for several system sizes in a square geometry, N = L × L. Figure 7(a) shows the ground state energy as a function of J (at h = 1) for L = 5 to 16, obtained from MPS numerics.For the smallest size we also show ED results, which coincide with the MPS ones.The kink near J = 1 indicates a quantum phase transition.Note that this behaviour is similar in systems with a single classical ground state (L = 5, 8, 16) or multiple ones (L = 6, 7), cf.Table I.In Figs.7(b,c) we show the average transverse magnetisation, M x = 1 N i X i , and M zzz , respectively, for systems with L a power of two.We get exactly the same results for different system sizes too.Both MPS and ctQMC show clear indications of a first-order transition at J = 1 in both observables.
Figure 8 shows similar results in a rectangular geometry, N = 3×M .For such thin stripe systems we can perform MPS more efficiently for larger system sizes than for square geometries.Once again, MPS and ctQMC results coincide, and indicate a first-order transition at J = 1 (although weaker than in the square lattice case, in the sense that the discontinuity in the local operators shown is smaller).Note that these results include not only values of M which are multiples of three, for which there ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ Average three-spin interaction as a function of J. are multiple classical ground state cycles, but also values of M for which a single ground state is found.What we see in this case is that the observables M x and M zzz are unable to detect changes related to any given classical ground states.

Symmetry breaking
In Figs. 7 and 8, we show the two terms that compete in the Hamiltonian, M x and M zzz .For system sizes where there is one classical ground state and no nontrivial symmetries, the total longitudinal magnetisation M z can also serve as an order parameter, as it picks up the orientation of the ground state.Figure 9(a) shows that the transition is also clear for this observable for square lattices.
For system sizes where degeneracies are expected for J ≥ h, however, M z is unable to detect the symmetry breaking related to the extra symmetries.For these cases, we need the staggered magnetisations, M m z , such as that for N = 3 × 3 in Eq. (20). Figure 9(b) shows that such operators are able to detect the spontaneous breaking of symmetry for these lattices.Note that Fig. 9(b) was obtained through the use of a small symmetry breaking field.This is a standard method for the detection of the symmetry breaking in the ground state of a degenerate quantum model [57].As a result, the calculations were performed through the use of a modified Hamilto-nian H = H QTPM − p Mz , where p is chosen to be small.The detection of the spontaneous symmetry breaking can be similarly preformed for any of the classical ground states of the given lattice size with the appropriate operator Mz .
In order to more clearly understand the mechanism of the phase transition, in Figs. 10 and 11 we plot the lowlying spectrum of the QTPM from ED as a function of h for fixed J.These results support our above observations: for system sizes where only a first-order phase transition is expected, there is an avoided crossing between the ground state and the first excited state; for system sizes with extra symmetries from the cycles of Rule 60, we see both an avoided crossing (indicative of first-order transitions) and a merging of eigenstates indicative of spontaneous symmetry breaking.As seen in Fig. 11 for the case of N = 3 × M , the avoided crossing becomes apparent only with increasing system size.
We now comment on how our results compare to those in Ref. [17].For the numerics, Ref. [17] used a stochastic series expansion (SSE) approach.We in turn use MPS and ctQMC.Both SSE and ctQMC are Quantum Monte Carlo based methods, which indicates that they, in principle, should be able to roughly access system sizes of the same order of magnitude.
Furthermore, while Ref. [17] also considered periodic boundaries, there was no specific restriction on system size, and therefore no distinction between sizes for which there is a single classical minimum and sizes where there are multiple ones, with the implications for symmetries of the corresponding QTPM.Ref. [17] also used a non-local order parameter, compared to our local ones (the staggered magnetisations) that do reflect the minima of the underlying TPM.In [17], the existence of a phase transition at J = h was confirmed through the study of the Binder cumulant; this was done, however, with limited accuracy on the location of the phase transition point.It is important to note that some of the local observables we calculate here are also studied for specific system sizes in the Appendix of Ref. [17].Since the temperature used for those calculations varied for different system sizes, it is possible that the smoothness observed in Ref. [17] is a consequence of thermal effects.We instead used a fixed inverse temperature β = 128 which we verified is sufficient to make thermal effects negligible.
D. Nature of the phase transition in the thermodynamic limit The discussion above and the numerical results indicate the existence of a quantum phase transition in the thermodynamic limit, N → ∞, at the self-dual point, J = h, of the QTPM.However, the approach to the thermodynamic limit is different across different system size geometries.
There are three different limits to thermodynamics: (i) across one of the two dimensions while the other one remains constant (that is, infinite stripes), (ii) across both dimensions, and (iii) on making the spins continuous.We briefly discuss the differences between these limits and the complications that might arise.
In the case (i), if the limit is taken for fixed L and with M such that lcm C|M (e.g.M = 3k, with k ∈ N), the number of classical ground states remains the same.In our numerics we are restricted to narrow stripes to allow convergence of the MPS algorithm.Fig. 8 suggests that in such quasi-1D systems the transition will eventually be slighly weaker than for square system sizes.
Case (ii) can be more involved.The simplest situation is that of square lattices N = L × L with L a power of two, where it is guaranteed that for all sizes there will be a single classical ground state, and therefore the transition is certainly first-order.For other size sequences, the number of relevant Rule 60 cycles, and therefore symmetries of the QTPM, may grow or decline with system size.For some cases this growth is monotonic (as for example for N = 3 k × 3 k with k → ∞), while in others it is not (as for example when N = 3k × 3k with k → ∞), see Table I.
In case (iii) the nature of the underlying CA is altered [58,59].In this limit, Rule 60 becomes f (p, q, r) = p + q − 2pq, (21) where p, q and r indicate the state of the three sites in the neighbourhood determining the local evolution of the CA, see Fig. 2. Basic arguments [58] indicate a single fixed point in the evolution of this fuzzy CA.We speculate that the same behaviour will be observed in the quantum field theory limit for QTPM; a single ground state across different regions of the whole J − h space and thus a first-order phase transition.However, a field theoretic description of the QTPM might not be as obvious and straightforward to get for the above elementary argument to hold.

V. CONCLUSIONS
In this work, we used the cycles of the cellular automaton Rule 60 to describe the symmetries of the quantum triangular plaquette model.We found that the attractor structure of Rule 60 plays an important role in the characterization of the degeneracies of the ground states of the classical TPM, allowing in turn to construct the symmetry operators of the QTPM.In this way, the existence or absence of stable cycles in Rule 60 imply whether it is possible or not for the QTPM to display spontaneous symmetry breaking of the corresponding symmetries, which in turn impacts on the nature of the quantum phase transition at the self-dual point.These general observations are also consistent with the finite size trends from our numerical simulations.A full description of the QTPM phase transition would require a field theoretical description and a renormalization group treatment; we leave these tasks for future works.
FIG. 11: Same as Fig. 10 but for systems with multiple symmetries.(a The merging of the ground state (black squares) with three degenerate excited states (purple rectangles) is indicative of spontaneous symmetry breaking.
In the main text, we show that ground state properties of the TPM and, consequently, the quantum phase transition of the QTPM depends on the boundary conditions, but we only focused on systems with periodic boundaries.Here, we consider the cases of periodic boundaries in only the x-dimension (PBCx) and of open boundaries (OBC), using the same Rule 60 CA considerations as for the periodic case.
For PBCx, we use the update rule for Rule 60 as before, with the only difference that we do not need to explicitly check for the periodicity across the y-direction.As a result, for an initial array of L sites, there will be 2 L configurations and, thus, 2 L ground states for the classical TPM.In this case, only the number of sites in the x-direction matters for the number of classical ground states.For example, a lattice with size N = 3 × 3 and one with N = 3 × 80 will have the same number, 8, of ground states.The identification of the classical ground states can be worked out from the Rule 60 evolution, as before.
For OBC, the update rule for Rule 60 is modified for the first cell of an L-length array so that it is not updated.This freedom on choosing two of the boundaries of the lattice gives an increased number of ground states for the classical TPM.Specifically, given a lattice of N spins, N = L × M , the number of the classical ground states is 2 L+M −1 .
We now perform a similar numerical analysis as in Sec.IV C, but only using MPS methods.Data is normalised with the system size of the given lattice.The system sizes accessible do not give a clear indication of a well-formed phase transition, but only signatures of it.The first-order transition found is weaker than in the case of the fully PBC, which we attribute to the high number of ground states for the classical TPM, given the system sizes.We note that all these states for h = 0 constitute low-lying excited states which affect the convergence of the MPS algorithm.
As seen from Figs. 12 and 13 for PBCx, the difference between the square lattice size scaling and the quasi-1D rectangular stripes is more pronounced, when compared to the finite-size scaling for PBC.Extra calculations on wider rectangular stripes verify that this difference is only a feature of the quasi-1D geometry of the lattice and not an inherent property of the system.Accuracy is lost with increasing size and the MPS results for the sizes studied are not reflective of the true thermodynamic limit.
The above considerations for PBCx are even more noticeable for the case of OBC, compare Figs. 14 and 15.For the quasi-1D stripes we see smooth behaviour for all values of J, while they seem to have converged to their "thermodynamic" behaviour.However, as seen from the square system sizes, the behaviour of the model remains the same regardless of the boundary conditions.It becomes apparent though that bigger system sizes soon become computationally inaccessible due to the exponential number of classical ground states.This behaviour shows an obvious discrepancy with standard MPS methods; normally, fully periodic system sizes are computationally harder to access.Here, since the model can have an exponential number of low-lying states(ground states for h = 0) for OBC, the convergence of the algorithm is hindered and significantly increases the lattice size where the "thermodynamic limit" has been reached.Therefore, only for PBC, the thermodynamic limit becomes apparent for the sizes we can access.
The significance of these arguments is further evident from Figs. 16 and 17.For the case of PBCx, all degenerate ground states for the J h region are easily found from exact diagonalization calculations and classically excited states are easily tractable too.However, the same is not true for OBC.The number of classically degenerate ground states increases exponentially and this is the reason why it would be pointless to show more ground states.

Appendix B: Gap Scaling Analysis for PBC
In this section we present a restricted and with limited accuracy analysis on the energy difference between the ground state and the first excited state.This analysis was conducted based on ED and MPS methods, which limits the validity of the conclusions that can be reached: it becomes quickly obvious that MPS methods are not powerful enough for the detection of the actual gap, especially in regions of the parameter space with high entanglement or with a high number of low-lying excited states, where MPS often converge to excited states above the lowestlying ones.However, the analysis below still provides an indication of the behaviour of the gap with system size when comparing systems with different symmetries.
This limited accuracy when measuring the first excited state energy is evident in Fig. 18(a).Data is calculated for the J = h = 1.0 point.The gap seems to decrease with increasing the system size, but at the same time, the power of MPS to detect it is significantly reduced.The situation seems clearer for Fig. 18(b).However, it is equally problematic despite the monotonically decreasing gap.The only significance of these results are an upper bounds of the actual gap.For the first case, the gap seems to decrease algebraically to zero, while for the case of multiple classical ground states, it seems to decrease  exponentially.This underlines the different behaviour depending on the existence or not of multiple classical ground states.
The same problems are encountered close to the phase transition from the MPS results in Fig. 19(a).Both plots are normalised by the maximum value of the gap encountered in the region of J values studied.For 0.0 < J < 2.0, the gap appears always to be maximum at J = 2.0 for Fig. 19(a) and at J = 0 for Fig. 19(b).In Fig. 19(b), for J > h = 1.0 the gap approaches zero, as expected from the existence of degenerate ground states.▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲

FIG. 3 :
FIG. 3: (a) Evolution from a single site for Rule 60.(b) Same for Rule 90.(c) A stable cycle generated by Rule 60.

FIG. 4 :
FIG. 4: The fixed point and limit cycles of Rule 60.Filled circles indicate states of the CA and the arrows the flow under the dynamics.(a) For L = 5 there is a cycle of length C = 15 and the fixed point (which flows into itself).(b) For L = 6 there are two cycles of period C = 6, a cycle of period C = 3 and the fixed point.(c) For L = 8, there is the fixed point and no cycles.

FIG. 6 :
FIG.6: Normalized three-spin correlator, Eq.(19), in the QTPM as a function of J for fixed h = 1, with N = L × L with L a power of two.We compare results from MPS and ctQMC obtained here with results fromRef.[16].The numerical data indicates a first-order transition at J = h.

FIG. 7 :
FIG. 7: First-order transition in the QTPM for systems of size N = L × L. (a) The normalised by the system size ground state energy as a function of J at fixed h = 1.Open symbols are results from ED, filled symbols from numerical MPS.(b) Transverse magnetisation as a function of J.In this case the open symbols are from ctQMC.(c)Average three-spin interaction as a function of J.

FIG. 9 :
FIG. 9: (a) Longitudinal magnetisation for systems with no symmetries.(b) Staggered magnetisation for detecting symmetry breaking in systems with multiple symmetries.

FIG. 12 :FIG. 13 :
FIG. 12: (a) Normalised energy, (b) M x and (c) M zzz of ground states from MPS and ED for square lattices with PBCx.Data from ED are denoted as empty squares and empty triangles.

FIG. 14 :FIG. 15 :
FIG. 14: (a) Normalised energy, (b) M x and (c) M zzz of ground states from MPS and ED for square lattices with OBC.Data from ED are denoted as empty squares and empty triangles.

30 FIG. 18 :
FIG.18:The scaling of the gap, g, for different lattice sizes without (a) and with (b) symmetries for the QTPM for the J = h = 1/0 point.Both ED and MPS methods are used (where appropriate) for the calculation of the given gaps.Square and rectangular sizes are equally used.

FIG. 19 :
FIG.19:The gap, g, normalised with the maximum gap, g max , in the given domain for different lattice sizes from MPS without (a) and with (b) symmetries for the QTPM.For (a) g max ≈ 3.43 − 3.50 and for (b) g max = 2.0.

TABLE I :
Cycle lengths for Rule 60.