Higher-Order Cellular Automata Generated Symmetry-Protected Topological Phases and Detection Through Multi-Point Strange Correlators

In computer and system sciences, higher-order cellular automata (HOCA) are a type of cellular automata that evolve over multiple time steps and generate complex patterns, which have various applications such as secret sharing schemes, data compression, and image encryption. In this paper, we introduce HOCA to quantum many-body physics and construct a series of symmetry-protected topological (SPT) phases of matter, in which symmetries are supported on a great variety of subsystems embbeded in the SPT bulk. We call these phases HOCA-generated SPT (HGSPT) phases. Specifically, we show that HOCA can generate not only well-understood SPTs with symmetries supported on either regular (e.g., line-like subsystems in the 2D cluster model) or fractal subsystems, but also a large class of unexplored SPTs with symmetries supported on more choices of subsystems. One example is \textit{mixed-subsystem SPT} that has either fractal and line-like subsystem symmetries simultaneously or two distinct types of fractal symmetries simultaneously. Another example is \textit{chaotic-subsystem SPT} in which chaotic-looking symmetries are significantly different from and thus cannot reduce to fractal or regular subsystem symmetries. We also introduce a new notation system to characterize HGSPTs. We prove that all possible subsystem symmetries in square lattice can be locally simulated by an HOCA generated symmetry. As the usual two-point strange correlators are trivial in most HGSPTs, we find that the nontrivial SPT orders can be detected by what we call \textit{multi-point strange correlators}. We propose a universal procedure to design the spatial configuration of the multi-point strange correlators for a given HGSPT phase. Specifically, we find deep connections between multi-point strange correlators and the spurious topological entanglement entropy (STEE), both exhibiting long range behavior in SRE states.


I. INTRODUCTION
Cellular automata (CA) are dynamic systems that evolve in discrete time steps, which have been wided used in comptuer and system sciences [1].CA have rather simple evolution rules, but produce rich structures.Because of their ability to model a wide range of phenomena, they have been used to model various real-world systems and can be used for prediction and simulation [2,3].
In the field of condensed matter physics, CA are often adopted to simulate dynamical properties of systems, such as Ref. [7].An example is the quantum cellular automata (QCA), originating from von Neumann and Feynman [8][9][10].QCA consist of arrays of identical finite-dimensional quantum systems that evolve in discrete-time steps by iterating a unitary operator U [11].QCA are useful for simulating quantum systems and processes, such as quantum walks, quantum circuits, and quantum phase transitions [12][13][14][15][16]. Apart from simulation, CA also play a role in the study of symmetry-protected topological (SPT) order.Let us review some basic facts of SPT physics.SPT phases are short-range entangled states that cannot be smoothly deformed into trivial states without breaking some symmetries .SPT phases have been extensively explored through various methods, sparking interests from fields like condensed matter physics, mathemati- * yepeng5@mail.sysu.edu.cncal physics and quantum information.These methods include group cohomology [40], cobordism groups [41,42], nonlinear sigma models (NLSM) [43,44], topological field theories [45][46][47][48], conformal field theories (CFT) [49][50][51], decoration construction [52], topological response theory [53][54][55][56][57][58][59], projective or parton construction [60][61][62][63][64][65], and braiding statistics [66][67][68][69].However, all the above SPTs are limited to the cases where the symmetries are global, meaning that they act uniformly on the whole system.Recently, partially motivated by the field of fracton physics , people have realized that not only symmetries themselves but also where exactly symmetries act on the system matters, catalyzing the research of symmetry protected topological phases protected by symmetries supported on either regular (e.g., line-like, membranelike) [4,5,[107][108][109][110][111][112] or fractal [6,113,114] subsystems.However, systematically designing the Hamiltonian in order to get the model protected by symmetries acting on a specific kind of subsets of the system remains a challenging problem.On one hand, the Hamiltonian of a lattice model usually involves some products of local operators; on the other hand, we need to control the subsets that the symmetries act on (usually the subsets are non-local and spread through the system).
Remarkably, CA coincidentally have simple local rules but exhibit complex behaviors, making them an ideal paradigm to design the Hamiltonian of condensed matter physics.So far, order-1 linear cellular automata have been used to construct FSPT models whose nontrivial edge states are protected by fractal subsystem symmetries [6].Although fractal geometry in physics is also an interesting topic [115][116][117][118][119][120][121][122][123][124][125][126][127][128][129][130][131], order-1 lin-Table I. Representative examples of SPT phases generated by HOCA.Strange correlators of model-IVc and model-Vb can also be calculated by the same procedure given in the paper, but is not explicitly shown in the table as these models are not the main focus of this paper.
Fig. 1.A brief schematic introduction to the types of SPT models produced by (linear) CA of different orders.The order-1 CA can only produce fractal patterns (see Appendix A for detailed argument).It would require HOCA, i.e., CA with order n ≥ 2, to create SPT models like RSPT, CSPT, and MSPT.These are SPT models whose symmetries are supported on regular subsets of lattice, chaotic-looking subsets of lattice, and more than one type of subsets of lattice, respectively.
ear CA cannot produce subsystem symmetries without selfsimilarity (see Appendix A).Therefore, it is natural to ask whether we can include SPT orders protected by symmetries supported on regular subsystems in the framework of CA 1 , and whether there are other possible subsets supporting symmetry action.
In this paper, we go beyond the linear order-1 CA by using linear higher-order cellular automata (HOCA) to generate SPTs with various kinds of subsystem symmetries.We call these phases HOCA-generated SPT (HGSPT) phases.HOCA are cellular automata whose evolution involves multiple time steps [132], and are widely used in computer science [133][134][135][136]. HOCA produce a rich variety of subsystem patterns in the spacetime lattice, including line-like and fractal patterns.HOCA have local update rules that make them useful for constructing Hamiltonians 2 .By using HOCA, we obtain a series of exactly solvable lattice models with various types of subsystem symmetries, as shown in Table I and Fig. 1.These models include not only SPT models with symmetries supported on regular or fractal subsystems, but also more peculiar models.Therefore, in this paper we propose a notation for the types of SPT orders protected by symmetries supported on subsystems, which includes regular(subsystems) SPT (RSPT), fractal(-subsystem) SPT (FSPT), mixed(-subsystem) SPT (MSPT) and chaotic(-subsystem) SPT (SPT) orders.For example, the 2D cluster model with linear subsystem symmetries discussed in Ref. [5] and order-1 CA generated models with fractal subsystem symmetries discussed in Ref. [6] are respectively classified into RSPT 3 and FSPT orders.Besides, we also have MSPT models with both ified.
line-like and fractal-like subsystem symmetries, MSPT models with two distinct types of fractal-like subsystem symmetries and CSPT models with chaotic-looking subsystem symmetries.A more detailed and technical definition of these types of SPT orders is given in Sec.III G.We also introduce a new notation system to characterize these newly discovered subsystem symmetries, claiming that HOCA patterns can be labeled and HOCA rules can be classified by the patterns they produce.This notation plays the role of an attempt towards constructing a universal classification system of all possible configurations of subsystem symmetries, as the spatial form of symmetry elements has not been mathematically labeled as has symmetry group itself.We also discuss how universal the HOCA symmetries are, i.e. if any kind of subsystem symmetries can be understood as an HOCA symmetry.Nevertheless, it is particularly worth noticing that the CSPT is a class of SPT models supported on symmetries with highly exotic spatial distribution, which is very challenging to extract the mathematical properties of these symmetry patterns and labeling them.Because of the application of HOCA in the realm of data processing and encryption, we deem that CSPT models are also applicable in quantum computation, being an resource of quantum encryption algorithm, which is left to future exploration.
To detect the nontriviality of a given SPT ordered ground state, one may use strange correlators [137].By definition, a usual strange correlator is introduced as a two-point correlation function in which bra-and ket-wave functions are a symmetric short-range entangled state and the state to be diagnosed, respectively.If the state to be diagnosed is indeed SPTordered, the strange correlator will either saturate to a constant or decay algebraically at long distances.While nontrivial phenomena of SPT order are fully characterized by the boundary with 't Hooft anomaly, strange correlators enable us to detect the nontrivial SPT order directly from the bulk, which removes potential analytic and numerical complexity induced by intricate boundary conditions.So far, strange correlators have been successfully applied to many SPT phases [4,[137][138][139][140][141][142][143], and have also been applied in intrinsic topological orders and conformal field theories (CFT) [144][145][146][147][148][149][150][151].In particular, in Ref. [4], this tool has been successfully applied to a 2D SPT order protected by line-like subsystem symmetries, i.e., a RSPT order following the nomenclature of the present paper.Therefore, one may wonder how to detect SPT phases with other types of subsystem symmetries, such as FSPT, MSPT, and CSPT discussed in the present paper.It is also interesting to ask whether or not the usual definition of two-point strange correlators is sufficient to detect all HGSPTs.
In this paper, we find that HGSPTs can be efficiently detected by what we call multi-point strange correlators (MPSC).More concretely, in some HGSPT models, 2-point strange correlators are insufficient for probing SPT orders, which leads us to generalize the usual strange correlators to multi-point.This approach reveals the complexity of SPT the previous work corresponds to the special case that we call "RSPT" in this paper, and we will avoid using "SSPT" to prevent confusion.physics induced by HOCA and also expands the research scope of strange correlators.To explore this topic, we design a general procedure to detect the nontrivial SPT orders in the HGSPT models and to determine the class of HOCA update rules for a given HGSPT phase.We have shown that there are models that can only be detected by MPSC with more than 2 points by rigorous mathematical proof.We explicitly present the multi-point strange correlators designed for the models discussed in this paper.By generalizing the strange correlator to multi-point strange correlator, we find that the spatial properties of the symmetry can be reflected by the configuration of multi-point strange correlator, exhibiting the complexity of the HOCA evolution.
It is also worth noticing that we discover MPSC, as the long range behavior in a short range entangled ground state, is inextricably connected to the spurious topological entanglement entropy (STEE) [152], a "spurious" long range behavior in SSPT models devoid of topological order.We show that the nonlocal stabilizers that can run along the boundary in boundary geometries like Levin-Wen prescription [153] or dumbbell-like tripartition in [152], have the exact form of MPSC in these HGSPT models.We demonstrate this claim by many concrete models and mathematical proof in the paper.We discover a large variety of MPSC that can serve as the nonlocal stabilizers giving STEE in multiple HGSPT models, showing that STEE is a common character of most SPT orders with subsystem symmetry, beyond the scope of models protected by line-like symmetries only.To enable spurious contributions of these stabilizers, one must consider more generalized boundary geometries, e.g.staggered boundary and even detached boundary.MPSC and STEE are both long range behaviors in a short range entangled state, between which the relation is a profound topic.
Hopefully, the multi-point strange correlator can potentially serve as a quantum mechanical approach to surpass the computational irreducibility [1,154,155] of HOCA, since we can effectively verify whether an arbitrary step of HOCA evolution is the given configuration by measuring the corresponding multi-point strange correlator (see also Sec.IV B 4).These emergent phenomena naturally urge us to establish a more general and fundamental theory of strange correlators, shedding light on its underlying physics and explaining its efficacy in probing SPT phases, which is left to future exploration.
The rest of this paper is organized as follows.In Sec.II, we provide some basic knowledge of HOCA.In Sec.III, we present the details of HGSPTs (see also Table I), including the edge states, symmetry protection, duality, and concrete examples of models.In Sec.III G, we introduce a notation system to label HOCA rules and patterns, and we give a technical definition of the types of SPT orders protected by subsystems symmetries base on these notations.Some typical examples are summarized in Table II.In Sec.IV, we find that the phases mentioned above can be detected via multi-point strange correlators, and we propose a general procedure to design multipoint strange correlators for a given HGSPT model.We apply the procedure to all models discussed in Sec.III and show the results in Sec.IV B. In Sec.V we show the relation between MPSC and STEE by some concrete examples.In Ap-pendix A, we prove that order-1 CA cannot produce genuine RSPT phases as HOCA can, demonstrating the need for introducing HOCA.In Appendix B we compare the model-IVa and the commonly known 2D cluster model, showing their subtle differences.In Appendix C, we demonstrate that for an FSPT model, strange correlators with 2 onsite operators are all trivial.And as a special case, in Appendix D we prove that there must be at least 3 points in the multi-point strange correlator to detect the nontriviality of an FSPT model.We further add a part comparing different CA approaches in constructing subsystem symmetries in Appendix E, and a part that gives a brief review of calculating several dynamical properties of HOCA in Appendix F. The mathematical discussion of 2 criteria in Sec.IV A on how to detect the class of HGSPT phases using MPSC is given in Appendix G. Finally, a mathematical discussion on the universality of HGSPT phase is given in Appendix H.

II. HIGHER-ORDER CELLULAR AUTOMATA (HOCA)
A. Preliminaries of CA and HOCA Cellular automata (CA), first introduced by von Neumann [8], have been recognized as a good dynamical system for simulating complex physical systems.CA have a simple structure but exhibit a great variety of complex behaviors, and are used to model phenomena with local, uniform, and synchronous processing [156].Formally speaking, a CA consists of an infinite set of identical finite automata placed over a lattice and all taking a state from a finite set called the alphabet of the CA.
There are many possible variants of CA.People have explored CA in higher dimensions to model systems with multiple degrees of freedom, and higher-order CA with memory size n > 1, which is the main focus of this paper.A higherorder cellular automaton (HOCA) is a discrete dynamic system whose evolution involves multiple time steps, first introduced by Toffoli in [132].While ordinary linear cellular automata always generate self-similar patterns (e.g., fractal patterns in the spacetime lattice, which can be proved by the Freshman's Dream theorem of a polynomial over F p [157], see also Appendix A), (linear) HOCA produce many peculiar patterns besides fractal patterns 4 .For example, HOCA can exhibit chaotic behaviors, which are often used in secret sharing schemes [134,135], data compression and image encryption [133].The encryption algorithm based on HOCA can be efficiently implemented in hardware due to the simple structure of CA, and is hard to decipher due to the chaotic behavior of the HOCA.
Despite the above interesting applications in computer science, the understanding of the dynamic behavior of HOCA is still at an early stage, and few results are known for linear 4 Linear CA rules can be written in terms of polynomial representations, making it possible to write Hamiltonians with respect to these rules.
HOCA [156].More properties and applications of HOCA are still to be studied and explored.Furthermore, to the best of our knowledge, HOCA have not been used in the realm of physics so far, and this paper will serve as an attempt to explore the interdisciplinary amalgamation of HOCA and condensed matter physics.Now we introduce some basic notations of HOCA.Consider a set of 1D lattice sites {i}, i ∈ Z with alphabet a i ∈ {0, 1, • • • , p − 1} = F p evolving with time j; the state of any given site at any given time may be expressed as a i (j).We introduce the polynomial representations to simplify our notation.By doing the substitution we express the spacetime configuration of all lattice sites by a polynomial: Also, we define site configuration at time j 0 with respect to x as by picking all terms with y-exponent equals to j 0 .Notice that our model is defined on a semi-infinite plane here, which shows the entire evolution of the HOCA rule.The HOCA model can be also defined on an open slab by truncation, which we will introduce later in subsection III A. We will use r j (x) to denote the configuration at time j from now on.Now we introduce the concept of higher-order cellular automata (HOCA), which are extensions of traditional cellular automata that involve interactions across multiple time steps.In an order-n CA 5 , the state of a site at time j 0 is determined by the states of a neighborhood of sites at times From now on, we focus on HOCA defined on F 2 = {0, 1}.Every HOCA rule mentioned below is defined on F 2 .If a i0j0 can be written as translationally invariant sums of elements in where c pq ∈ F 2 are coefficients, R is radius, a constant describing the maximal range of p, which does not scale with the system size, making the rule local.We concentrate on linear HOCA because update rule of a linear order-n HOCA can be represented by n polynomials, which enables us to construct Hamiltonians with decorated defect construction using these update rules.We demand R < ∞ to make sure the HOCA rule is local, which means that the effect of the HOCA rule (i.e.change of a ij due to the HOCA rule) will not propagate faster than the speed R. We denote an HOCA rule (i.e.update rule) by an n-row vector f , dubbed as the update rule of the HOCA: where the superscript T denotes the transpose of the vector.
When n = 1, the HOCA returns to the normal CA as discussed in [6].The time evolution of local linear HOCA can be denoted by a single formula (here we assume j 0 > n): where r j (x) is defined in Eq. (3).
To ascertain the whole time evolution process of all lattice sites of an order-n CA, one needs to manually specify the configurations of first n time steps r 0 , r 1 , ..., r n−1 , which is called the initial condition of the system.It can also be denoted by an n-row vector q(x): By specifying an HOCA rule f and an initial condition q, the whole spacetime pattern can be uniquely defined.We define and so on, such that r n−1+i (x) = q T (x) • E (i) (f ).E (i) is dubbed as the evolution operator, which can be calculated using Eq. ( 6).We can always write r j (x), j ≥ n as sum of each row in the initial condition multiplied by some update rules.It follows that the whole spacetime pattern can be expressed as where q and F capture the effect of the initial condition and the update rule separately, and the label y p,q is defined as Eq. ( 9) is useful in the calculation of commutation polynomial (Eq.( 15)), which is important to the discussion of the symmetry elements in the HOCA generated SPT phases (to be discussed in section III A).If we treat the time axis as another spatial dimension, we can view the whole time evolution of the given HOCA F (x, y) as a static pattern in a 2D semiinfinite plane.Any given HOCA rule f can generate infinite number of patterns by adjusting initial condition q(x).

B. Spin (qubit) model in terms of polynomial representations
Polynomial representations can also express spin systems by identifying a ij with the state of the spin located at site x i y j .Then the whole HOCA pattern F (x, y) (given in Eq. ( 9)) naturally expresses the spin configuration in the lattice.By introducing polynomial representations, we naturally transplant HOCA into the realm of spin systems.Consider a spin model defined on a d-dimensional square lattice with α sublattices (i.e. each site contains α independent degrees of freedom, which don't have to equal the order of the HOCA).One spin is placed on each site of the sublattice.We introduce the following conventions • The coordinate of site s = (i 1 , i 2 , ..., i d ) are represented by a monomial m with respect to In this paper, we focus on d = 2 case, and we use the notation i 1 ≡ i, i 2 ≡ j, x 1 ≡ x, x 2 ≡ y.
• Previously defined a ij expresses the state of spin located at site x i y j (in a specific sublattice).a ij = 0 represents that the spin at x i y j is at the state |0⟩, and a ij = 1 represents |1⟩.The sublattice that the spin belongs to will be defined in the next point.
where m k denotes the position of site s k in sublattice i that the operator O acts nontrivially on, and P k denotes a polynomial with respect to x.
• We define the coefficients of monomials m (i) k to be in Z 2 , and then we can naturally obtain . . .
III. HOCA GENERATED SYMMETRY-PROTECTED TOPOLOGICAL PHASES A. Lattice models, short-range entanglement, symmetry, and symmetry protected edge states In this section we construct symmetry-protected topological (SPT) phases protected by HOCA generated symmetry.Due to the great variety of the HOCA behavior, we can naturally obtain models with fractal symmetries, line-like symmetries, both of the above, and even chaotic symmetries, which are respectively classified into FSPT, RSPT, MSPT and CSPT orders according to our notation.And as a summary, in Sec.III G, we give the technical definition of all these types of SPT orders based on the concrete examples demonstrated in this section.We'll begin with some basics and notations.
Through HOCA and decorated defect construction, we can obtain symmetry-protected topological (SPT) order with two types of subsystem symmetry, which we refer to as mixedsubsystem SPT (MSPT).Among all MSPT models, there are models with both fractal-like and line-like subsystem symmetries (abbreviated as I-MSPT), and models with 2 different fractal symmetries (abbreviated as II-MSPT).Additionally, there are models with only chaotic-looking symmetries (dubbed as chaotic SPT) and models with line-like and membrane-like symmetries (resembling the previously known SSPT model).These HOCA generated models are all defined on a 2D square lattice with 2 sublattices (a) and (b), and the Hamiltonians can be generally written as where f is an HOCA update rule (see Eq. ( 5)), ȳ is a vector composed of monomials of y (see Eq. ( 10)), and the notation The Hamiltonian in Eq. ( 13) describes an exactly solvable cluster model with a short-range entangled unique ground state on a torus, similar to the usual cluster states.In fact, we can obtain the commonly seen 1D cluster model (which is an SPT phase protected by global Z 2 × Z 2 symmetry) by taking f = 1, and define it on a lattice with L x = 1.The Hamiltonian we get are equivalent to that of the 1D cluster model up to a change of basis (Z ↔ X).The exact solvability of the model can be proved by noting that there are always 0 or 2 overlapping operators between two terms in the Hamiltonian, ensuring that every Hamiltonian term commutes with each other.Two examples are shown in Fig. 2. The property of the ground state can be verified by noting that there are 2L x L y qubits and 2L x L y Hamiltonian terms in a model defined on a L x ×L y torus.Notice that each Hamiltonian term corresponds to a unique onsite Pauli operator, resulting that all Hamiltonian terms are independent with each other.For example, each term Z x i y j (1 + f • ȳ1,n ) x i y j corresponds to operator Z 0 x i y j , and each term X With no other constraints being present, the ground state subspace has dimension 2 2LxLy /2 2LxLy = 2 0 = 1, giving a unique ground state on the torus.Now, let's investigate deeper into the symmetry elements of the model.Suppose these models are defined on an open slab with a size of L x × L y , and all Hamiltonian terms with operators outside of the boundary are excluded.Here, L x is the length in the i direction and L y is the length in the j direction.For a SPT model generated by an order-n CA, the open slab should satisfy L x ≥ p min + p max and L y > n , where p max and −p min are respetively the largest and smallest power of x in f (x) (if p max or p min are less than zero, then it is defined to be zero), to ensure there is at least one valid Hamiltonian term in the model.
We now focus on Hamiltonian terms whose coordinate x i y j is in the slab but contains sites outside of the slab.Assuming the coordinate axis is taken as in Fig. 3, then for sublattice (a) there are n rows of such Hamiltonian terms excluded at the top edge of the system, p max terms at the left edge, and p min terms at the right edge.Each excluded term with coordinate x i y j plays the role of a lost constraint on the ground state manifold, producing a free spin at site x i y j .Similarly, we can obtain extra degrees of freedom at sublattice (b), with everything reversed (top↔down, left↔right, etc.).Suppose there are k such excluded Hamiltonian terms in the model, then the ground state degeneracy of the model will be 2 k .
We are now able to flip free spins at the edge without changing the energy of the system.Flipping spins at the edge will generally affect spins in the bulk following the HOCA update rule, producing symmetry elements in the shape of the HOCA pattern.The operations that flip spins in these HOCA patterns commute with the Hamiltonian of the model, being the symmetries that protect the degenerate edge state.These model has 2 sets of subsystem symmetries S (a) and S (b) , each set per sublattice: where F (x, y) is the truncated HOCA pattern F (x, y) (Eq.( 9)) specified by an HOCA rule f and an initial condition q, and all terms which are not fully in the slab are excluded.The HOCA rule f controls which type of subsystem symmetry can be found in this model, and q controls the specific pattern of the symmetry.We can enumerate these symmetry elements in an HOCA generated SPT model by counting all possible different initial conditions that can be defined in the slab.
These symmetry elements commute with the Hamiltonian terms, which can be verified by examining the commutation polynomial.The commutation polynomial with respect to two polynomials α, β is defined as: P (α, β) ≡ α β.If the coefficient of x 0 y 0 in P (α, β) is zero, X(α) and Z(β) commute with each other [6].Now we verify the commutation relation in sublattice (a).We will do the calculation in the semi-infinite plane (−∞ < i < ∞, j ≥ 0) and do the truncation afterwards.In sublattice (a), the symmetry writes (sublattice (b) are not represented below) X(F (x, y)) and a general Hamiltonian term (we consider the Z-term only) writes Z(x i y j (1 + f • ȳ1,n )).The commutation polynomial is (15) Here, r(x) is defined as follows: (16) where all terms have a power of y that is lower than n.Note that i ∈ Z and j ∈ N. The definition of F can be found in Eq. ( 9).While the exponent of y is given by k − j in Eq. ( 15), and k = 0, 1, ..., n − 1, which means that all terms in the commutation polynomial have a y-power less than n.According to our convention, Hamiltonian terms with j < n are all excluded.Hamiltonian terms with j ≥ n have a null commutation polynomial since the y-power of all terms in the commutation polynomial are less than zero, which are also outside of the slab.Thus, we proved that such symmetry elements indeed commute with the Hamiltonian (Eq.( 13)).
The edge states are protected by the above symmetry elements.These symmetry elements are all in the shape of an HOCA pattern (and their superposition after translation), which can be generated by choosing an initial condition q and truncating the resulting HOCA pattern to fit the open slab.The visual property of the symmetry element is controlled by both initial condition q and HOCA rule f .To explore edge physics, we can define a series of new Pauli operators at the edge.Taking sublattice (a) as an example, the edge Pauli matrices are written as: where these operators are truncated to the slab by default.Here, edge states in sublattice (a) are distributed along the top, left and right edges of the slab.These three matrices all commute with remaining Hamiltonian terms and they form a Pauli algebra.To open the gap of a degenerate edge state, we can add a magnetic field to an edge free spin.This operation must violate approximately 2 k−1 symmetry elements which act nontrivially on this edge spin, reducing the ground state degeneracy by half.Now we give a brief picture of the duality of HGSPT model.On open slab, each ground state of an HGSPT model can be mapped to a symmetry breaking ground state of the dual model.By redefining Pauli operators, an HGSPT Hamiltonian becomes 2 decoupled copies of symmetry breaking orders with the same HOCA generated symmetry, each showing a phase transition at magnetic field h = 1 via Kramers-Wannier duality.So for an HGSPT model, we can add magnetic field in two different directions (h x and h z ), and transition will happen at h x = 1 and h z = 1.If both magnetic fields are smaller than 1, the model remains in the HGSPT order.If one of h x and h z is bigger than 1, the model becomes symmetry breaking phase in one sublattice.If h x , h z are both bigger than 1, the system are in the trivial paramagnetic phase.
We will give several examples of HGSPT models in the next few subsections.
B. Model-I: I-MSPT generated by order-2 CA

An example of I-MSPT Hamiltonian writes
generated by an order-2 CA with update rule The Hamiltonian is pictorially shown in Fig. 3.The order of the HOCA is n = 2, so we need to specify 2 rows of initial conditions by a 2-row vector q.Now we focus on symmetry elements in sublattice (a), and symmetries in sublattice (b) can be obtained similarly (by reversing everything, as mentioned in subsection III A).First, take an initial condition q and plug it into Eq.( 9) to get the symmetry pattern.Then, if the model is defined on an L x × L y open slab, we place the first row of symmetry pattern on top of the slab (row with j = 0) and exclude parts that are not in the slab.Since the HOCA rule is translationally invariant in x-axis, we can move our pattern in the x-direction with terms outside of the slab being excluded.The operation above gives us F , which we plug into Eq.( 14) to get the symmetry element.In the case of sublattice (a), the symmetry elements are made up to Pauli-X operators in the shape of F , just like we defined in Eq. ( 14).Consider symmetry elements generated by the following 4 initial conditions: The overall results are shown in Fig. 4. It can be seen clearly that a Sierpinski triangle (Fig. 4(a)) and a line (Fig. 4(b)) can both be the symmetry element, and there are symmetry elements that look like the attachment of two patterns (Fig. 4(c)).
Although we only consider the initial condition with the absolute values of exponent of x (denoted as P ) in each row less or equal than 1, any q can be chosen in principle if it can fit into the size of the open slab.However, doing so does not bring us extra peculiar phenomenon.So far we have not found guiding principles of ascertaining q for a given type of symmetry element (e.g.line-like or fractal), and 4 initial conditions mentioned above are found by computer enumeration.
The edge states of the model are protected by the symmetries we generated above (and other possible HOCA generated symmetries).In Fig. 5 we explicitly show an example of symmetry protection for model (Eq.( 17)).In the figure we show that to open the gap of an edge spin must violate two symmetries: a line-like symmetry and a fractal-like symmetry (and many other HOCA generated symmetries that act nontrivially on this spin).The only way to modify edge spins while keeping all commutation relation with symmetries is to couple edge spins at different edges, which is either non-local or located at the corner of the system.

C. Model-II: II-MSPT generated by order-3 CA
There are also models with 2 different fractal symmetries.Consider a Hamiltonian generated by an order-3 HOCA rule The Hamiltonian is pictorially shown in Fig. 6.Examine symmetries generated by the following initial conditions: The overall results are pictorially shown in Fig. 7.When there is only one flipped spin in the initial condition, the HOCA rule gives an chaotic pattern (shown in Fig. 7(a)).20)).They are both Sierpinski triangles but with different shapes and orientation.Fig. 7(c) can be viewed as the attachment of two fractal symmetries (with small modifications in the middle).

D. Model-III: CSPT generated by order-3 CA
Chaotic SPT (CSPT) models only contain subsystem symmetries in chaotic patterns.Under various initial conditions, the symmetry elements majorly show chaotic patterns.HOCA rules producing chaotic patterns are often used for encryption algorithm in computer science, as minor change in the initial condition may produce entirely different chaotic pattern.An example writes The Hamiltonian is pictorially shown in Fig. 8. .Two examples of symmetries that protect the edge free spin located at green circle at the top edge of the system.The Hamiltonian is (Eq.( 17)).If we manually break the degenerate edge mode at green circle (e.g. by adding a Zeeman term), such modification of the Hamiltonian will anticommute with two symmetry elements shown in the figure (and many other terms that act nontrivially on this site), showing that the edge mode is indeed protected by our HOCA generated symmetry.Given 4 initial conditions The resulting subsystem symmetries are shown in Fig. 9.
E. Model-IVa: RSPT generated by order-2 CA It has been known that 2D regular SPT (RSPT, previously referred to as SSPT in Ref. [5]) can be generated by a cluster model.Now we want to show that our HOCA framework also includes quantum models with line-like and membranelike symmetry elements just like Z sub 2 strong SSPT discussed in [5].Our model is different from the model discussed in [5], but some of their subsystem symmetries share the same type.A rigorous proof of this statement is shown in Appendix B. In addition, there are also checkerboard-like membrane symmetries (Fig. 11(d)) in our model, which is different from the previously defined SSPT model.A typical update rule of such an RSPT model writes which generates the Hamiltonian Hamiltonian generated by this rule are pictorially shown in Fig. 10.
Consider following initial conditions: There are line-like and membrane-like symmetry elements present in the model, as shown in Fig. 11.
which generates the Hamiltonian The Hamiltonian of Model-IVb is shown pictorially in Fig. 12.Given 4 initial conditions:   29)).Black lattice and blue dashed lattice denote 2 sublattices.
we obtain 4 different subsystem symmetries shown in Fig. 13.Different from Model-IVa which can only generate regu-lar patterns, there are some chaotic-looking symmetries in Model-IVb.At first look some symmetries (e.g.Fig. 13(d)) may seem to possess a fractal structure, but they do not actually have a rigorous self-similarity (see Sec. III G for more detailed discussions).
G. A notation system for labeling the higher-order cellular automata Now, we want to introduce new notations to characterize various subsystem patterns generated by HOCA.As we have shown in the previous sections, HOCA can generate various types of patterns in the spacetime lattice, such as fractal patterns with rigorous self-similarity (e.g., Sierpinski triangle), patterns that consist of some periodic repetition of basic structures (e.g., checkerboard), and even patterns that look like a mixture of the two above.We also found these subsystems have various dimensions (in the sense of Hausdorff dimension).We claim that these properties of an HOCA pattern F (defined on the semi-infinite plane, see Eq. ( 9)) generated by a finite initial condition (i.e.there are finite terms in the initial condition) can be captured by a mathematical object   quantity, M , is dubbed as mix rate, describing how fractal or periodic the pattern is.For an order-n HOCA pattern, M is mathematically defined as where and where A(i) is the number of cells in state 1 of the i-th row.
The definition of M comes from following observations.Now we have observed two possible local behaviors of a HOCA evolution patterns for all HOCA rules with radius r ≤ 3/2 and n ≤ 3: 1. Self-similar fractal structure: Some parts of a HOCA pattern tend to appear recurringly while we increase the time of evolution (i.e.zooming out the pattern).While the whole pattern may not be fully self-similar in general, there are recognizable fractal structure in many patterns.An example is Fig. 7(d).
2. Regular structure: Some parts of a HOCA pattern may appear to be filled by some local repeating structures, like Fig. 11.
While in general a HOCA pattern may not be a fully fractal or regular pattern, a large subset of HOCA patterns can be viewed as some mixture of fractal and regular patterns.This visual observation can be clearly seen for almost all HOCA patterns with r = 3/2 and n ≤ 3, while it is subtle to argue whether HOCA patterns like Fig. 9 can be viewed as this kind of mixture just by watching.To give a more quantitative and rigorous description of this mixing behavior, we observe that fractal patterns and regular patterns are distinguishable in terms of counting cells with state 1 in every row.We denote the number of cells with state 1 in row t to be A(t) as in Eq. 33, 34, then we have following qualitative observations: 1.In a fully fractal pattern, there always exist a infinite sequence {t i } such that A(t i )/t i → 0. This observation is obvious because of the self-similar essence of fractal.Self-similarity means there are infinitely many rows can be represented by simply scaling a single row configurations.An example of this are the sequence {t = 2 n − 1 : n ∈ Z + } of Fig. 7(d).All rows with index in this sequence have 2 cells in state 1, scaled by different proportions.
2. In a fully fractal pattern, there is always a sequence {t i } such that the sequence {A(t i )} grows to infinity.A heuristic explanation of this observation are fractal patterns always tend to grow bigger as t increases.A rigorous proof of this statement would involve the sensitivity of HOCA rule, which will be introduced in Section F 4. In short, this statement is always true for a fractal pattern.
3. In a fully regular pattern, where every row configuration can be viewed as some repetitions of a specific local structure (e.g. each row in Fig. 11(d) can be viewed as the repetitions of a white-black checkerboard structure), the number of sites with state 1 in each row either grows linearly (for a membrane-like pattern) , or remains a constant (for a line-like pattern), or oscillate between two constants (for a zipper-like pattern).In all situations the upper limit and the lower limit of the sequence {A(t)} when t → ∞ will be at the same order, being both O(t) or O(1).
The observation above enables us to define the "mix rate" M of two types of structures in a single HOCA pattern.As seen in the definition of M (Eq.( 32), ( 33),( 34)), for a fully fractal pattern there will be a constant upper limit of {A(t)/t}, and a zero lower limit of {A(t)/t} (all constant lower limits of {A(t)} will be suppressed by the 1/t factor).Thus, M = 1 is always true for the fractal case.As for the case of regular pattern, the upper and lower limit of {A(t)/t} will either be equal to an identical constant (membrane-like pattern) or equal to zero (line-like, zipper-like case).Both cases lead to the result that M = 0, being the character of a regular pattern.Generally, M is in [0, 1] for a general HOCA pattern, and this property is true for any finite HOCA pattern generated by finite initial condition (only sites in a limited area have nonzero values).Specifically, the sum in the numerator of Eq.( 33), (34) is taken from i = k to i = k + n because the fact that there will be at most n − 1 consecutive empty rows in a pattern generated by an order-n HOCA, which can cause some non-fractal pattern to obtain a zero lower limit of {A(t)}, affecting the effectiveness of M .As the research on the dynamical properties of HOCA is still at early stages, there are still no mathematical paper to give rigorous classifcations of these properties.Therefore, we have done an exhaustive and case-by-case verification of the validity of M for all r ≤ 3/2 and n ≤ 3 HOCA rules (512 rules in total), and have found no counterexamples of our classification based on our current observation.While there are certainly mathematical foundations of this quantity, the topic is beyond the scope of this research and is left as a part of future works.
We provide a pictorial description of various HOCA patterns: 1. M = 1: It is a fractal pattern exhibiting a self-similar structure, like Fig. 7(d).Also there are mixed patterns that can be considered as the attachment of (d, 1) fractal pattern and a (1, 0) pattern, as shown in Fig. 4(c).
2. M = 0: (d, 0) pattern can be considered as spatial repeating of some minimal structures or some stable patterns.The overall pattern can extend in the 2D plane (Fig. 11(d)) or propagate along some 1D subsystem of the plane (Fig. 11(b)).
3. 0 < M < 1: Chaotic patterns can show recognizable repeating structures locally or appear to be irregular, but it does not fit into the classification above, as shown in Fig. 9.
More explicit examples of the validity of M in classifying HOCA is shown in Fig. 14, 15.
An HOCA rule can generate infinite patterns by varying the initial condition q(x).However, if we collect all possible patterns generated by an HOCA rule, we find that different rules may produce different types of patterns, dividing HOCA rules into 4 classes.We use X (f ) = (X r , X f ) to denote the classes [f ], where where ⌈ ⌉ and ⌊ ⌋ are ceil and floor functions, respectively.{M } represents the set of all possible M generated by the given HOCA rule.Heuristically, X r and X f describe whether a certain HOCA rule can generate regular pattern or fractal pattern.For example, X r = 1 means that there are at least one regular pattern (e.g.line, membrane, checkerboard, etc.) can be obtained from the HOCA rule by varying the initial condition, and vice versa.Typical examples of patterns above can be found in Fig. 14, 15.Given a specific update rule f , different patterns can emerge when we adjust the initial conditions.Thus, we can classify different update rules by the patterns they can produce, as shown below: • X (f ) = (0, 0): These HOCA rules only produces chaotic patterns like symmetry elements presented in chaotic SPT.Neither fractal nor like-like patterns can be found in this class.
• X (f ) = (0, 1): HOCA rules in this class can produces fractal patterns but not periodic patterns.Sierpinski FSPT [6], and previously mentioned II-MSPT can be generated by CAs in this class.
• X (f ) = (1, 1): HOCA rules in this class produces both fractal and periodic patterns.These rule can generate I-MSPT phases.
To capture finer details of an HOCA rule, we define two sets of characteristic dimension of an HOCA rule f : and where F denotes any possible HOCA pattern generated by the HOCA rule f , and d(F ) denotes the box dimension of tation in hand, we can give a technical definition of RSPT, FSPT, MSPT and CSPT orders as follows: • Regular(-subsystem) SPT (RSPT) is the SPT phases protected by subsystem symmetries that necessarily (i) include regular subsystem symmetries (e.g.line-like symmetry) and (ii) exclude fractal subsystem symmetries.For HGSPT models, it means that RSPT models correspond to HOCA rules f satisfying X (f ) = (1, 0).
• Type-I mixed(-subsystem) SPT (I-MSPT) is the SPT phases protected by subsystem symmetries that necessarily (i) include regular subsystem symmetries (e.g.line-like symmetry) and (ii) include fractal subsystem symmetries.For HGSPT models, it means that I-MSPT models correspond to HOCA rules f satisfying X (f ) = (1, 1).
Given the technical definitions above, there are still points need further clarification.Firstly, we can notice that these definitions are not completely intuitive: for example, when regular symmetries and chaotic-looking symmetries exist simultaneously in one model and fractal symmetries do not, the model is classified as RSPT phase as well.Besides, II-MSPT orders with two different kinds of fractal subsystem symmetries do not have a specific position in this classification, as purely according to the X (f ) of HOCA rules they would be classified into FSPT orders.Furthermore, for a HOCA rule f a rigorous proof between X (f ) = (0, 0) and chaos is still lacked, although in our observation X (f ) = (0, 0) always implies the HOCA rule can generate chaotic-looking patterns.
A finer classification of HGSPT orders naturally depends on a more complete and sophisticated understanding of the dynamics of HOCA, thus it is beyond the scope of this paper, but we expect it to be an important future direction which may lead to further understanding of subsystem symmetries.
Here we list some SPT phases characterized by our notation in Table II

IV. MULTI-POINT STRANGE CORRELATOR DETECTION
Originally proposed in [137], "strange correlator" is a powerful tool to detect nontrivial short-range entangled states.Recently, strange correlators have been used to detect the nontriviality of the RSPT state [4] (see footnote 3).This work shows that RSPT state with line-like subsystem symmetries can be detected by strange correlators with two operators ϕ being in the same straight line corresponding to the anisotropy of the subsystem symmetries.This naturally motivates us to detect the nontriviality of HGSPT phases through strange correlator, with the hope that the configuration of the operators inside the strange correlator will reflect the property of the HOCA generated symmetry of the model.
In the previous sections, we have shown that HOCA are able to successfully generate SPTs protected by various kinds of subsystem symmetries.In the following, given a specific HGSPT, we want to detect its nontriviality and the class that its HOCA rule belongs to.We will show that this task can be completed by what we call "multi-point strange correlator" (MPSC).

A. Definition
The strange correlator is defined as follows in Ref. [137]: where |Ψ⟩ is the short range entangled (SRE) state to be diagnosed, |Ω⟩ is the trivial disordered state in the same Hilbert space as |Ψ⟩, ϕ is some local operator.For nontrivial SRE states, the strange correlator will saturate to a constant or undergo a power law decay for specific ϕ, while that of trivial SRE states will decay exponentially or become null.The strange correlator defined above involves 2 local operators, and the definition can be extended to the case of n local operators, dubbed as multi-point strange correlator: Here, n are dubbed as the correlation number of the multipoint strange correlator.We also introduce the multi-point normal correlator, which can be regarded as the strange correlator of the trivial disordered state, serving as the "background" to be subtracted from the strange correlator: For a given ϕ and a spatial configuration {r i }, we say the strange correlator C({r i }) gives nontrivial result if and only if C({r i }) − N ({r i }) saturate to a constant or decay algebraically.Since if C({r i }) is (quasi-)long range ordered but C({r i }) − N ({r i }) is not, it would mean that we cannot distinguish non-trivial SPT ordered states from the trivial symmetric state with this strange correlator.In this work we demand ϕ to be onsite Pauli operators.By means of multi-point strange correlator, we construct a general procedure that detects the nontriviality of the HGSPT ground state.
It is also worth noticing that using the duality relation between the SPT model and the symmetry breaking model, the MPSC can be mapped to the membrane-like order parameters in the symmetry breaking models, some examples have previously shown in [6,158,159].However, it has not been discussed up to our knowledge that whether the ground state of SPT models with fractal symmetries can only be detected by strange correlators with more than 2 local operators, i.e. if there are "intrinsic" multi-point nature in these models.In the following texts, we explore systematically the behavior of MPSC in various HGSPT models and proved that there is indeed SPT models that can only be detected by strange correlators with more than 2 points, as shown in Appendix D.

B. Detection through multi-point strange correlators
In this subsection we are to probe the nontrivial ground states of HGSPT models.The aim of this subsection is to raise a universal approach that distinguishes HGSPT models in different classes.Now that we have demonstrated that all HGSPT models have degenerate edge states on an open slab in the previous section, it is guaranteed that we can find a specific ϕ and a particular spatial configuration of ϕ to produce nontrivial results.The point here is to find the ϕ and configuration that give nontrivial results while reflecting the symmetry properties of the system.
If we denote the position of operators in the multi-point strange correlator by {r i }, then for a HGSPT generated by an order-n CA, we claim that the nontrivial ground state of the HGSPT model can be detected by the multi-point strange correlator in the following configuration: where and we choose the trivial symmetric state |Ω⟩ = X(a) = Ẑ(b) = 1 .Specifically, m 0 (x) = 0, and m i (x), i > 0 can be calculated by and All calculations of polynomials above can be easily done by computer.
L plays the role of "distance" in this configuration and takes value in N, and is named evolution distance hereafter.Given the HOCA rule of the HGSPT model, we can construct a series of multi-point strange correlator by fixing an initial condition q (as long as the initial condition can be defined in the given bulk) and increase L. Different series of multi-point strange correlators will behave variously depending on the X (f ) and q.Specifically, we define the correlation number n of a strange correlator to be the number of onsite operators ϕ included in the correlator n[D L (q, f ; x, y)].Then, we claim that by observing how n grows with L will be helpful to determine the X (f ) (Sec.III G), the class of the HOCA rule for the given HGSPT phase.
The seemingly complicated configuration can be interpreted below: q(x) can be recognized as the initial condition q(x) that generates the whole configuration.Since generally the strange correlator is acted in the bulk, so the configuration q(x) may possibly violate the HOCA update rule, resulting in trivial result of strange correlator, forcing us to introduce the term m.The term m(x) is added for two reasons: (i) m i (x) is added to each q i (x) to make sure that the product of operators commute with the symmetry, and is determined by q(x); (ii) (a) Model-I (Eq.( 17)), initial condition are given by Eq. (19b), L = 14, correlation number n = 6.
Fig. 16.Examples of q and f that make N (q, f ) satisfies criterion 1 or 2. The configuration Dn,L(q; x, y) generated by the given q and f are shown in the figure above.Detailed results are shown in Table III.
p i (x) has the form of HOCA pattern generated by q(x) in the corresponding rows to meet the commutation relations.Such construction ensures that X 0 D L (q, f ; x, y) (take sublattice (b) as an example) act equivalently as products of Hamiltonian terms in the given sublattice, giving a trivial action on the ground state |Ψ⟩.This guarantees that the multi-point strange correlator C(q, f , L; {r i }) gives nontrivial results since (48) Here H X α,β is the Hamiltonian term defined on the site denoted by α, β and composed of X operators, where α, β sum over all the sites with nontrivial values in D L (q, f ; x, y).Now we give some comments of this construction of configuration.There are two main goals that we consider while designing this specific configuration.First, we want to find out the simplest configuration of ϕ needed to show the nontrivial results in an HGSPT ground state.Second, we hope that the configuration we design will reflect the symmetry property of the given HGSPT model.The second goal is well-achieved in our construction in all models discussed in this paper, while the first goal is not always easy to satisfy.Generally, we conjecture that configuration that meets the first goal are always included in this configuration.This claim is proved rigorously in the case of model-Va (Eq.( 55)), an FSPT model.
For an order-n HOCA generated SPT, we can determine the class [f ] by examining the behavior of correlation number n of the strange correlator in the configuration Eq. ( 42) as L → ∞.
Based on the definition of correlation numbern[D L (q, f ; x, y)], we can further define the following quantity: where and Then the following two criteria holds: • Criterion 1: • Criterion 2: Fig. 16 shows some explicit examples of these two criteria.So far the search for the initial conditions that satisfy criterion 1 or 2 is done by computers, and hopefully an analytical way will be found in future works.The detailed mathematical discussion and concrete examples of these two criteria are shown in Appendix G. Now we apply these criteria to HOCA generated SPT phases that we have discussed before: Model Number Xr X f Criterion 1 Criterion 2 I (Eq.( 17)) 1 1 Yes (Fig. 16(a)) Yes (Fig. 16(b)) II (Eq.( 20)) 0 1 No Yes (Fig. 16(c)) III (Eq.( 23)) 0 0 No No IVa (Eq.( 27)) 1 0 Yes (Fig. 16(d)) No IVb (Eq.( 30)) 1 0 Yes (Fig. 20) No Va (Eq.( 55 It can be seen in Table III that the class [f ] can be detected by multi-point strange correlator.Examples that meet the criteria in the table are shown in Fig. 16.

Detecting FSPT (Model-Va) generated by order-1 CA
In this section, we will probe the nontriviality of the Sierpinski FSPT ground state via the strange correlator.The Sierpinski FSPT model has degenerate edge modes, which are localized states at the boundary of the system that are protected by the fractal symmetry of the model.These edge modes are argued in detail in [6], and they indicate that the ground state of the model must be a nontrivial short-range entangled (SRE) state, which is a quantum state that cannot be transformed into a product state by local unitary operations.We take the Sierpinski model as a model generated by order-1 CA, and explore its nontriviality by means of multi-point strange correlator.The update rule of the Sierpinski FSPT model [6] is and the Hamiltonian is written as Applying our procedure to the model (Eq.( 55)), we find that when where q = 1 and f = 1 + x, then criterion 2 is satisfied: Also, there are no possible configuration that satisfies criterion 1, which can be proved by the self-similarity nature of the order-1 CA, or simply by enumearting all possible initial conditions on an open slab.Thus, we indeed find out that the Sierpinski rule (Eq.( 54)) is in the (0, in the case of the Sierpinski FSPT model, and these 3 Pauli matrices must be placed at the 3 corners of a Sierpinski triangle in the lattice.This claim is proved in appendix D.

Detecting I-MSPT (Model-I) generated by order-2 CA
For model-I (Eq.( 17)), we expect that both criteria above can be satisfied by controlling the initial condition.We observe that the multi-point strange correlator satisfies when the initial condition q is set to be Eq.(19c).The correlation number of the configurations above are of which 3 examples are shown in Fig. 18(b), Fig. 18(c), Fig. 18(d).Also, we find that the multi-point strange correlator Fig. 17.Pictorial illustration of configuration of multi-point strange correlator DL(q, f ; x, y) of update rule Eq. ( 54), which generates an FSPT phase (nontrivial terms are shown in blue cubes in Fig. 17   satisfies when the initial condition q is set to be Eq.(19b).The correlation number of the configurations above are of which an example is shown in Fig. 16(a).

Detecting II-MSPT (Model-II) generated by order-3 CA
For model-II (Eq.( 20)), we expect that criterion 2 can be satisfied by controlling the initial condition.We observe that the multi-point strange correlator satisfies when the initial condition q is set to be Eq.( 22d).The correlation number of the configurations above are of which an example is shown in Fig. 16(c).For model-IV (Eq.( 23)), we expect that no criterion can be satisfied.This claim is confirmed by computational search on initial conditions with size L x ≤ 50, and increasing the size generally do not give any new phenomenon.For all possible configurations we observe the correlation number n generally increase with L. Unlike other models mentioned in this paper, fixing any initial condition q, we will never obtain an infinite sequence of L that makes the multi-point strange correlators share the same correlation number n.Among all strange correlators in this model, the one with minimal correlation number writes where the initial condition is set to be Eq.(25a), of which the figure is shown in Fig. 19.We notice that multi-point strange correlators in CSPT order seem to give a promising procedure to overcome the computation irreducibility of CA.While the computational irreducibility states that we cannot directly calculate an arbitrary step in CA evolution (for CA showing complex behaviors, not CA with regular and predicable patterns, e.g.HOCA rules in CSPT models) without calculating steps before the wanted step, in principle we can efficiently prepare the ground state of this model in an array of qubits and measure the strange correlator by a series of quantum operations in this qubit array.Only the multi-point strange correlator with the correct con- figuration will show nontrivial result.That is to say, we can verify whether the result of an arbitrary step is a given configuration with zero knowledge about the intermediate steps, which can potentially serve as an quantum approach to surpass the well-known principle of computational irreducibility [1,154,155].

Detecting RSPT (Model-IVa) generated by order-2 CA
For model-IVa (Eq.( 27)), we expect that criterion 1 can be satisfied.We observe that the multi-point strange correlator satisfies when the initial condition q is set to be Eq.(28a).The correlation number of the configurations above are of which an example is shown in Fig. 16(d).For model-IVb (Eq.( 29)), we expect that criterion 1 can be satisfied.We observe that the multi-point strange correlator satisfies when the initial condition q is set to be Eq.(31a).The correlation number of the configurations above are of which an example is shown in Fig. 20.

V. MULTI-POINT STRANGE CORRELATOR AND SPURIOUS TOPOLOGICAL ENTANGLEMENT ENTROPY
Firstly studied in [160], it has been pointed out in Ref. [152] that the extraction of topological entanglement entropy (TEE) −γ via S topo (e.g.prescriptions proposed by Levin-Wen [153] and Kitaev-Preskill [161]) can suffer from spurious contributions from the nonlocal string order in the SSPT order, and this  spurious contribution is extensively studied [162][163][164].Authors of Ref. [152] have shown explicitly the string-like nonlocal stabilizer generators that spread across the boundary of subregions can contribute to the S topo in 2D cluster model, and have proposed a quantity S dumb to detect this spurious contri-bution.
In this work, we want to show that MPSC are closely related to the spurious contributions in these models, and the spurious contributions in calculating TEE exist in a large variety of SPT orderes protected by subsystem symmetry.We will show how spurious topological entanglement entropy (STEE) appears in the HGSPT order, and how the configurations of the nonlocal string-like stabilizers that contribute to STEE are exactly mapped to the spatial distributions of local operators in MPSC that can detect the nontrivial SRE ground state of an HGSPT order.
In the next following sections, we will discuss what kind of nonlocal stabilizers can contribute to STEE in Sec.V A, explore the connection between STEE and MPSC in Sec.V B, and show concrete models in Sec.V C.

A. Nonlocal stabilizers contributing to Stopo
In this subsection we will discuss what kind of nonlocal stabilizers can finally contribute to the calculation of S topo .
While extracting the topological entanglement entropy of a given physical model by means of tripartitions (e.g.Kitaev-Preskill and Levin-Wen), we argue that the calculation can be massively simplified by counting only a special type of nonlocal stabilizer generators.First we start with the entanglement entropy of a specific subregion A in the system, which is given by where N A is the number of qubits contained in region A, G A is the stabilizer group whose elements are fully supported in A. In gapped spin systems, the entropy scales with Here c is some non-universal constant, and the −γ term is the so-called topological entanglement entropy (TEE), which is a universal constant.Kitaev and Preskill as well as Levin and Wen proposed two tripartitions and used linear combinations of entanglement entropy of these parts to cancel out the non-universal constant and extract the constant term −γ in the scaling behavior of the entanglement entropy.The corresponding quantity is named topological entropy S topo , defined as It is argued that the combinations in S topo managed to cancel out all the boundary and corner contributions in the entanglement entropy, leaving However it is pointed out in [152] that S topo is not fully topologically invariant, being sensitive to the deformation of the region boundaries in the subsystem symmetry-protected topological (SSPT) phase.In SSPT phase, this extraction process may suffer from unwanted contributions due to long-range string order, giving a nonzero S topo in SSPT phase, being devoid of the topological order, which is different from our expectation that −γ = 0 in the SSPT phase.
Here we want to systematically explore the origin of the spurious contribution ∆S ≡ S topo + γ in the realm of the lattice model.
First we observe that N A∧B is the number of qubits that exist in AB but not in A or B. O A is the number of independent stabilizer generators in A, and O A∧B is the number of independent stabilizer generators that exists in AB but not in A or B. O A∨B is the number of independent stabilizer generators that exists in A or B, but become no longer independent in AB.The following relationship holds: Similar observations can be found for S ABC : Adding up each term in S topo , we can explicitly observe what terms are cancelled out in the linear combination.The calculation process are shown below: To explore what kind of operators can contribute to S topo eventually, we introduce some notation to keep track of terms that a specific stabilizer generator enters.The distribution vector D(g) of a valid generator (to be explained below) g is defined as where T (g) denotes the number of basic partitions (i.e. the unions of partitions are not included here) that g as a total has a support in.P (g) denotes the number of basic partitions that supports at least one local Hamiltonian term making up the generator.From now on, we will denote the support of the operator as total support and the union of supports of Hamiltonian terms that make up the operator as partial support.To specify the contribution of g in S topo , we further define the contribution vector C(g), which writes: where N is the number of partitions and C i (g) denotes how many times does g appear as the minimal generators in i-th order region.An i-th order region is the union of i basic partitions (e.g.ABC is a 3-rd order region), and a minimal generator in i-th order region is a generator that cannot be written as products of generators in j-th (j < i) regions with less area of support.We see explicitly in the definition of C(g) that the contribution of g in S topo is Through enumerations we find that for N = 3 case (including LW,KP prescription): And by calculation we observe that only the stabilizer generators with P (g) = 3 have a nonzero ∆S(g) in N = 3 case (LW, KP prescription).We define the power of set |{g : T (g) = i, P (g) = j}| as Q(i, j).Then we obtain: This quantity shows immediate potential to be generalized to N -partitions.Eq. (75, 76) can be interpreted as following: The stablizers that finally contribute to the calculation of S topo are the ones that have a local support on all 3 subregions, and their contributions depend on how many subregions their global support have.With this in hand, we can massively simplify the calculation of S topo by counting these special nonlocal stabilizers only.

B. Connection to multi-Point strange correlator
In this subsection we want to show that the nonlocal stabilizers will exactly be the operator that gives nontrivial results in the multi-point strange correlator.
In Eq. ( 76) we conclude that the stabilizers which can contribute to the spurious value of TEE have partial support that span across all three partitions.This kind of stabilizer either appears at the triple intersection point of partitions (KP case), or serves as nonlocal stabilizers running along the boundaries (both LW and KP case).Next we will show that the stabilizer in these cases will exactly be the operator that gives nontrivial results in the multi-point strange correlator.We will start with some explicit examples.In the following part of this section, we will show some nonlocal stabilizers that exists in some HGSPT models.These nonlocal stabilizers have nonzero contributions to the value of S topo in certain partition geometries.While it was pointed out in [152] that the string-like nonlocal stabilizers in 2D cluster model can be detected by calculating a tripartite topological entropy S dumb in a dumbbell-like tripartition, we want to show the following facts: 1.In HGSPT models with generally more exotic subsystem symmetries, there will be nonlocal stabilizers that contribute to the spurious values of S topo that cannot be detected by the original definition of S dumb .
2. There will also be string-like nonlocal stabilizers in the HGSPT models devoid of any line-like symmetries, which can also be detected by S dumb .So there is no general correspondence between the presence of linelike symmetries and nonlocal long range string order proposed in [152].
And it is worth noticing that these nonlocal stabilizers happen to be the configuration of operators that can detect the nontrivial SRE ground state of the corresponding HGSPT phase.It is natural to notice that both STEE and MPSC point out the fact that there are hidden long range order in the SPT phases protected by subsystem symmetries: STEE is the unexpected contribution to the TEE in the absence of the topological order, while MPSC is the hidden long range correlation behavior in a short range entangled ground state.By this exact relation we see that two quantities share the same physical origin.
In the following texts, we show some nonlocal stabilizers that can contribute to the tripartite topological entropy where the tripartitions A, B, C are denoted by blue, green, red areas respectively in the figure .As a reminder, we would like to clarify that when we say a nonlocal stabilizer is generated by a certain symmetry, we are actually saying that we pick certain rows from the symmetry pattern (Eq.( 14)) to be the initial condition that generates the MPSC (Eq.( 42)), and the resulting local operator configuration in the MPSC is the nonlocal stabilizer generated by this symmetry.From the discussion above it is clear that a certain symmetry pattern can generate a huge amount of MPSC, but not all MPSC can serve as the nonlocal stabilizers that contribute to S topo calculation.Only the ones that fit the boundary geometry of the tripartition can potentially have a spurious contribution.So we will explicitly draw the boundary geometry that admits spurious contributions of nonlocal stabilizers in the following texts.For the sake of simplifying the picture, we exerted a coarse-graining procedure on the lattice, combining two sublattices.Now the model is defined on a 2D square lattice with 2 qubits per site.A general Pauli operator acting on a site will be represented as O := O 1 O 2 , where O 1 denotes the operator acting on the sublattice (a) and O 2 denotes the operator acting on sublattice (b), respectively.
The reason why there exists such correspondence between the configuration of MPSC and the nonlocal stabilizers with spurious TEE contributions can be explained as follows: 1. Nonlocal stabilizers in the context of STEE are always made up of products of Hamiltonian terms along the boundary, which contain Hamiltonian terms with supports outside of the tripartition.However, the nonlocal stabilizer as a whole does not have support outside of the tripartition, making it an independent stabilizer generator of stabilizer group of area ABC (denoted as G ABC ).Details are discussed in Section V A.
2. Such nonlocal stabilzers naturally act trivially on the SRE ground state of the HGSPT phase as products of Hamiltonian terms, as shown in the second row of Eq. (48).

C. Model study
From the observation above we can see that any product of operators with the form designed in Eq. ( 42) automatically have the form of product of Hamiltonian terms therefore having the potential to contribute to STEE.The only thing we need to do is to analyze the boundary geometry that can admit such nonlocal stabilizers.In the following sections, we will show some concrete examples of nonlocal stabilizers together with the corresponding boundary geometries of tripartition.The possible nonlocal stabilizers that may appear in these models extend beyond the examples we will show below, so we will demonstrate some typical examples only.

Nonlocal Stabilizers in Model-I Model-I (a I-MSPT
model) possesses two types of subsystem symmetries, as mentioned in Fig. 4.Each type of subsystem symmetry corresponds to a class of multi-point strange correlators, giving rise to a class of nonlocal stabilizers with the same geometry.An example of nonlocal stabilizers generated by the fractal-like symmetry (Fig. 4(a)) is shown in Fig. 21(a), which gives a ∆S topo = −1 in the tripartite topological entropy.While the symmetry is fractal-like, the nonlocal stabilizer given by the symmetry has a string-like outlook, and can be prolonged in the i direction by appropriately tuning the initial conditions that generates Fig. 21(a).
There also exists another type of nonlocal stabilizers with different directions, which is generated by the linelike symmetries (Fig. 4(b)) in the model.
The corresponding MPSC configuration D L (q, f ; x, y) is shown in Table IV.The corresponding MPSC configuration D L (q, f ; x, y) is shown in Table IV.The corresponding MPSC configuration D L (q, f ; x, y) is shown in Table IV.The corresponding MPSC configuration D L (q, f ; x, y) is shown in Table IV Model-Vb is named Fibonacci FSPT in [6], with the update rule Despite the lack of line-like symmetries 6 , there are string-like nonlocal stabilizers that can contribute to S topo in this model.An example of nonlocal stabilizers generated by the fractal-like symmetry is shown in Fig. 21(g), which gives a ∆S topo = −1 in the tripartite topological entropy.The string-like stabilizer can be prolonged in the i direction.It is worth noticing that there are stabilizers in this model that can fit into a staggered boundary geometry, as shown in Fig. 22(d).
The corresponding MPSC configuration D L (q, f ; x, y) is shown in Table IV.

VI. SUMMARY AND OUTLOOK
In this work, we find exotic SPT phases protected by a variety of HOCA-generated symmetries.We identify HGSPT models with both fractal and line-like symmetries (e.g., Eq. ( 17)), models with two distinct fractal symmetries (e.g., Eq. ( 20)), and models with chaotic subsystem symmetries (e.g., Eq. ( 23)).These models are derived from the HOCA rule, as explained in Section III.We show that the framework of HOCA naturally encompasses these SPT phases protected by exotic symmetries and previously studied SSPT and FSPT phases, and we introduce labels that classify these SPT phases into different categories.To detect the nontrivial ground HGSPT order, we show the necessity of introducing multi-point strange correlators, which are a generalization of the strange correlator that involves more than two operators.The necessity is demonstrated by proving that all 2-point onsite strange correlators are trivial in the Sierpinski FSPT model, which is a fractal SPT model with a Sierpinski triangle symmetry (see Eq. ( 55)).This model can be recognized as an SPT phase generated by an order-1 CA, which can also be regarded as the simplest case of HOCA.By examining the multi-point strange correlator of the given phase, we can determine the class of the phase.Also, we have found the relation between the multi-point strange correlator and the nonlocal stabilizers resulting in spurious topological entanglement entropy, revealing the connection between these two quantities showing long range behaviors in a short range entangleed state.
There are many interesting topics that remain unsolved.For example, while HOCA can be used to construct SPT models, the symmetries supported on HOCA-generated configurations can also be utilized to build other phases, including symmetry-breaking phases and symmetry-enriched topological (SET) orders, where the order-1 CA case is done in [165].Such future directions may require us to use more than one HOCA rule or use HOCA in higher dimensions, and HOCA patterns in these generalized conditions remain a future direction to explore, where the order-1 CA case is discussed in [166].Moreover, the MPSC may be mapped to the membrane-like order parameters in the symmetry breaking models, of which some examples have previously shown in [6,158,159], showing the potential to probe new kinds of quantum criticality in symmetry breaking models.Another interesting topic is whether all types of subsystem symmetries can be generated by the HOCA framework, and how to develop a unified notation system to label miscellaneous subsystem symmetries.Finally, as chaotic patterns are realized as symmetries in CSPT ordered states, we may expect such quantum states to have diverse applications in computer science involving chaotic systems, such as for encryption [136].Moreover, the multi-point strange correlators in HGSPT models can be studied by the Monte Carlo method, where the strange correlators of RSPT models have been studied in Ref. [4].Entanglements in HGSPT models are also intriguing.For example, whether there are spurious topological entanglement entropy [152] in HGSPT models is an intriguing problem.Whether we can probe the HGSPT order via non-Hermitian perturbation [167] and study the non-Hermitian entanglement in HGSPT order [168,169] are also promising future directions.It will be interesting to probe the HGSPT in the realm of average symmetry-protected topological (ASPT) order [170,171], searching for peculiar behaviors of HGSPT models.Besides, using the relation between MPSC and STEE, we may probe a new way of detecting phases of matter via entanglements with geometric properties, paving the way for a general way to characterize exactly solvable models by analyzing their Hamiltonian terms without actually solving the model.Finally, because of the duality between the symmetry breaking model and the HGSPT model, it will be interesting to map the MPSC back to the symmetry breaking models to detect new kinds of quantum criticality.where |0 • • • 0⟩ denotes all qubits in the system are taken as +1 eigenstate of σ z operator.The trivial state should contain all symmetries in the ground state, so the trivial state is taken as For FSPT states, the first problem is to specify the operator ϕ, we will try a set of different operators as candidates for ϕ.
Overall results are shown in Table V. Candidates of local operators: Now we do a case-by-case calculation for all candidates, and notice that i ̸ = j is always assumed.

ϕ = X(a/b)
ij : There are intuitively 3 cases for the choice above: a) i ′ j ′ : They are both symmetry elements of the model, resulting that Table V. Strange and normal correlator in the FSPT model (Eq.( 55)) with correlation number 2 Operator ϕ(i, j)ϕ(i ′ , j ′ ) X(a) ij X(a) The strange correlator is Only one spin in the (b) sublattice is flipped.Such configuration cannot appear in the ground state, which makes (a) i ′ j ′ : 2 spins in the (a) sublattice are flipped, which is not a configuration in the ground state.
1 spin in the (a) sublattice is flipped, having zero overlap with the ground state.Therefore we have Using the data calculated above, we can easily obtain the result: convenient to consider only one sublattice at a time, and here we will discuss (b) sublattice, in which spins are all |0⟩ in the trivial state |Ω⟩, and operators with nontrivial action on the sublattice is ∆ ij .In the model, each sublattice is isomorphic to a square lattice.For simplicity, we denote the state of spin in (b) sublattice (|0⟩ or |1⟩) and its corresponding lattice site (i, j) by a matrix element s i+1,j+1 = 0 or 1.All s ij form a matrix Spin, shown in Fig. 26.
X(a) ij acts on our state, then we note this by o ij = 1.An example is shown in Fig. 27.
Our claim above is straightforward, which is basically a translation to the matrix language.Finally, we denote all elements in the i-th row of matrix A = {a ij } by a i * , and similarly j-th column by a * j .1. First we truncate Op by deleting all null columns and rows at the edge.We suppose that Op is a a × b matrix after the truncation.Lemma 2: (D4) Proof: 1. We assume that o ab = 0, Lemma 1.3 and the discussion 4(b) in the proof of Lemma 1 we know that N (Spin) ≥ 4, so our assumption is false and o ab = 1.  2. To satisfy the condition above, we observe that condition should be satisfied by all o ij , 1 < i ≤ a, 1 < j < b.
The proof is straightforward, Eq. (D9) is equivalent to s i,j+1 = 0. We also notice that this is intrinsically equivalent to the update rule of Sierpinski fractal (Eq.( 54)).
3. We also observe that otherwise extra flipped spin will be generated, and there are no way to cancel these extra spins.

(D12)
We can see that o ij is symmetric along the counterdiagonal by examining the variable substitution and equation Eq. (D12) is invariant under the substitution:  CA Method Lattice Dimension Lattice Geometry Discovered Subsystem Symmetries CQCA 2D [172,173] Square [172], 11 Archimedean lattices [173] Fractal [172,173], Line-like [172], Ribbon [173], Cone [ mal form, that is to say, it has the following form: where each m i (x) is a polynomial about the formal variable x.Topological conjugacy of LHOCA and LCA: The LHOCA defined by H = ⟨n, Z 2 , r, f ⟩, where f is specified by a j i , j ∈ [1, n], i ∈ [−r, r], can be simulated by an LCA L = ⟨Z n m , r, f ⟩, with f totally determined by f : In Ref. [156], it has been shown the above correspondence is a topological conjugacy that preserves dynamical properties, thus we can investigate the dynamical properties of LHOCA by considering the corresponding LCA.

D-chaos of LCA
In this appendix, we consider the well-accepted notion of chaos of discrete time dynamical system (DTDS) proposed by Devaney (often abbreviated as D-chaos).In general, the criteria of D-chaos is composed of three components: topological transitivity, sensitivity to initial conditions and denseness of periodic orbits [175].In this subsection, we briefly review the definition of these properties for DTDS, and introduce an algorithm that can decide whether an 1D LCA is D-chaotic or not proposed in Ref. [176].
Definition of discrete time dynamical system (DTDS): A discrete dynamic system is a pair (X , F) where X is a space equipped with a metric d, and F is a transformation on X which is continuous with respect to that metric.The dynamical evolution of a DTDS is described by an initial state x (0) ∈ X evolving as x (t) = F t (x (0) ), ∀t ∈ N.
In the realm of 1D CA, the space X is taken as S Z , where S is the alphabet of the CA.Therefore the space X represents the space of configurations at a specific time step.F is naturally the global rule of the CA.Here we take the metric as the standard Cantor distance where n = min{i ≥ 0 : Here c i denotes the state of configuration c at site i, which is an element of the alphabet S.
Definition of topological transitivity of DTDS: A DTDS (X , F) is said to have topological transitivity if for an arbitrary pair of open nonempty subsets of X , ∃n ∈ N, such that F n (U ) ∩ V ̸ = ∅.

Definition of sensitivity to the initial conditions of DTDS:
A DTDS (X , F) is said to be sensitive to the initial conditions if there exists ϵ > 0 such that ∀x ∈ X , δ > 0, then there exists y ∈ X , n ∈ N, such that d(x, y) < δ and d(F n (x), F n (y)) > ϵ.
Definition of denseness of periodic orbits of DTDS: An element x ∈ X is said to be a periodic point if there exists a natural number n > 0 such that F n (x) = x.The denseness of periodic orbits is the denseness of the set composed of all such periodic points.By definition, as a DTDS, an 1D LCA simultaneously satisfying the above three properties is chaotic according to Devaney's notion [175].

Deciding chaos of linear HOCA
In Ref. [176], the authors proposed an efficient method to decide whether a 1D LCA is D-chaotic or not.Obviously, 4. Using the additivity of the HOCA rule we immediately obtain that if a site is at state 0 at time t, then it will remain at state 0 governed by a radius-0 HOCA rule, which contradicts our initial assumption that the rule can produce a self-similar fractal pattern.
5. We conclude that all HOCA rules that can generate a fractal pattern are sensitive to the initial condition.
which suggests S u = S d = Const..This indicates that if we plug the corresponding initial condition into n inf [D L (q, f )] and n sup [D L (q, f )], we will obtain n inf = n sup = Const.So we have Appendix H: Mathematical discussion on the universality of HGSPT phases As HOCA managed to produce a large variety of symmetry patterns, one may wonder if any kind of subsystem symmetry can be generated by HOCA approach mentioned in the main text.In this section we will show the "completeness" of HGSPT model.
Proposition: Given a pattern S(x, y) defined on an open slab with size L x × L y , there will be at least one HOCA rule f and initial condition q, such that the HOCA configuration F (x, y) is identical to S(x, y) in the open slab.We say any finite patterns can always be locally simulated by an HOCA rule.
Proof: The proof of this proposition involves the topological transitivity of the HOCA rule (Appendix F 2): A DTDS (X , F) is said to have topological transitivity if for an arbitrary pair of open nonempty subsets U, V in X , then there exists a positive natural number n, such that F n (U ) ∩ V ̸ = ∅.
We start from a desired HOCA pattern S(x, y), and maps it to a single-row Frobenius LCA configuration α.The validity of this process is guaranteed by the topological conjugacy between an order-n HOCA over Z 2 and an LCA over Z n 2 , which can be considered as a single row of Z 2 vectors with n components.So this is a one-to-one map, without losing or adding any information.Then, we select an open set of configurations V that contains α locally, i.e. each configuration v ∈ V shares the same configuration with α within the domain of α, and can be arbitrarily chosen outside of the domain of α.Then, because of the topological transitivity of the LCA rule, if we pick an LCA rule with topological transitivity, then for any configuration subset U , there exists N ∈ N * such that F N (U )∩V ̸ = ∅, where F is the global rule of the LCA.Without lost of generality, we suppose that F N (u) ∈ V, u ∈ U .Then we can map LCA configuration u back to HOCA configuration u 0 .By selecting u 0 as the initial condition of HOCA rule, we obtain the desired symmetry pattern S(x, y) in an open slab after N steps of evolution governed by HOCA.This finishes our proof of the proposition above. Q.E.D.
The proof indicates that any order-n HOCA rule with topological transitivity have the ability to simulate symmetry patterns within an L x × L y slab with L y ≤ n.A pictorial illustration of the proof above is shown in Fig. 29.

•
If the onsite Pauli operators Pauli X, Pauli Ŷ and Pauli Ẑ operators are represented by O = X, Y, Z respectively, a many-body Pauli operator O can be denoted as

Fig. 2 .
Fig. 2. Two possible overlapping ways of Hamiltonian terms of I-MSPT model (Eq.(17)).Gray circles are the overlapping X and Z Pauli matrices from two terms.Black lattice and blue dashed lattice denote sublattices (a) and (b).

Fig. 7 (
b) and Fig. 7(d) are two fractal subsystem symmetries of the model (Eq.(

Fig. 5
Fig.5.Two examples of symmetries that protect the edge free spin located at green circle at the top edge of the system.The Hamiltonian is (Eq.(17)).If we manually break the degenerate edge mode at green circle (e.g. by adding a Zeeman term), such modification of the Hamiltonian will anticommute with two symmetry elements shown in the figure (and many other terms that act nontrivially on this site), showing that the edge mode is indeed protected by our HOCA generated symmetry.
Here, d is defined as the Hausdorff dimension of the pattern with infinite time evolution steps, which can be approached numerically by box dimension.If we denote the number of evolution time steps by t, and the number of sites with state 1 from time 0 to time t by a(t), then we have d = lim t→∞ ln a(t) ln t .For a Sierpinski triangle with Hausdorff dimension d H = ln 3/ ln 2 ≈ 1.5850, the numerical result with t = 256 gives d = 1.5830, which is quite close to the exact result.Another

Fig. 14 . 54 Fig. 15 .
Fig. 14.Example 1 (2 in total) illustrating the validity of M in classifying HOCA patterns.Fig. 14(a), 14(c), 14(e), 15(a), 15(c), 15(e) show 6 examples of HOCA evolution pattern.Fig. 14(b), 14(d), 14(f), 15(b), 15(d), 15(f) show how i+n k=i A(k) grows with i, where the numerical results of d, M is shown above each figure.Su and S d can be understood as the slope of green and orange straight lines in the subfigures in the right column.For regular patterns two slopes always equal, while for fractal pattern the slope of orange lines are always zero, for chaotic-looking patterns two straight lines have different nonzero slopes.

Fig. 18 .
Fig. 18.Pictorial illustration of configuration of multi-point strange correlator DL(q, f ; x, y) of update rule Eq. (18), which generates an I-MSPT phase (nontrivial terms are shown in blue cubes in Fig. 18(b), Fig. 18(c), Fig. 18(d)).3 configurations above have q(x) = 0 x −1 + 1 + x and n = 10.Fig. 18(a) shows the HOCA pattern generated by the initial condition above.It can be seen in the figure that 3 MPSC with different L share the same correlation number n, showing the scaling invariance of this MPSC.

Fig. 20 .
Fig. 20.Pictorial illustration of configuration of multi-point strange correlator DL(q, f ; x, y) of update rule Eq. (29), which generates an RSPT phase (nontrivial terms are shown in blue cubes).The initial condition is given by Eq. (31a).

2 .
Nonlocal Stabilizers in Model-II Model-II (a II-MSPT model) possesses two types of subsystem symmetries, as mentioned in Fig.7.Despite the lack of line-like symmetries, there are string-like nonlocal stabilizers that can contribute to S topo in this model.An example of nonlocal stabilizers generated by the fractal-like symmetry (Fig.7(d)) is shown in Fig.21(a), which gives a ∆S topo = −1 in the tripartite topological entropy.The string-like stabilizer can be prolonged in the i direction.

3 . 4 .
Nonlocal Stabilizers in Model-III Model-III (a CSPT model) possesses chaotic-looking subsystem symmetries only, as mentioned in Fig.9.So far we have not found any recognizable classes of nonlocal stabilizers that can contribute to the TEE due to the chaotic nature of the symmetry pattern.Nonlocal Stabilizers in Model-IVa Model-IVa (an RSPT model) possesses regular subsystem symmetries, as mentioned in Fig.11.There are string-like nonlocal stabilizers that can run along a smooth boundary (Fig.21(e)) and stabilizers that can fit into more peculiar boundary geometries (Fig.22(a), 22(c)).

5 .
Nonlocal Stabilizers in Model-IVb Model-IVb (an RSPT model) possesses regular subsystem symmetries, as mentioned in Fig.13.There are string-like nonlocal stabilizers that can run along a smooth boundary (Fig.22(b)) and stabilizers that can fit into more peculiar boundary geometries (Fig.22(a), 22(c)).

Fig. 22 .
Fig. 22. Nonlocal stabilizers that can contribute to Stopo in tripartition given in the figure.Red, green, and blue area is respectively the partition A, B, C. Such products of operators all have the same configuration with configuration of local operators in some certain MPSC.Nonlocal stabilizers in this figure can contribute to more exotic boundary geometry (dashed line surrounding green area).

1 Fig. 23 .
Fig. 23. 2 different symmetry elements in 3 different weak RSPT models with q = −1, 0, 1 respectively.It can be observed that the symmetries in weak RSPT are always non-overlapping.White cubes are the spins that the symmetry acts nontrivially on.The figure can be compared with Fig. 11, which shows the symmetries of a genuine RSPT phase.Notice that 3 subfigures above are symmetries from 3 different models while Fig. 11 shows 4 symmetries in the same model.

Fig. 24 .
Fig. 24.Fig. 24(a) shows the Hamiltonian of the 2D cluster model up to some basis transformation.Two sublattices are shown in the Fig. 24(b), where sublattice (1) are drawn in orange and gray dashed lines and sublattice (2) are drawn in green and gray dashed lines.There is 1 qubit living in each site.In Fig. 24(c) we deform 2 sublattices back to the form of a standard 2D square lattice, similar to the HGSPT model discussed in the main text.

:
2 spins in the (b) sublattice are flipped.Such configuration cannot appear in the ground state either, which makes C(i, j; i ′ , j ′ ) = 0.There are similarly 3 cases here:

2. Proof of Theorem 1 Lemma 1 :
If Op is not null, denote the number of nonzero values in matrix A by N (A), then N (Spin) ≥ 3, where Spin = Op ⋄ Ker.(D3) Proof:

2 . 3 . 4 .
Lemma 1.1 (a) Claim: Changing the value of an element o ∈ o a * , o ̸ = o ab from 0 to 1 will increase N (Spin) by 1.(b) Proof: Through direct calculation, o ai (1 ≤ i < b) = s a+1,i+1 .The claim above is obvious.Lemma 1.2 (a) Claim: If o * b is not null, then N (s * ,b+1 ) ≥ 2. (b) Proof: We examine o 1b , o 2b , • • • , o ab in order.Through direct calculation, each nonzero value o ib in the sequence will either increase N (s * ,b+1 ) by 2 (if o i−1,b = 0) or 0 (if o i−1,b ̸ = 0).Since there are at least 1 nonzero elements, so it follows that N (s * ,b+1 ) ≥ 2. It follows that o a * and o * b should contain at least one nonzero value (otherwise it would have been truncated).Note that s a+1, * is completed determined by o a * , and s * ,b+1 by o * b .There are 2 cases here: (a) o ab ̸ = 0: This already meets the condition above.From Lemma 1.2 we know that N (Spin) ≥ N (s * ,b+1 ) ≥ 2. (b) o ab = 0: It means that there are at least 1 nonzero value in {o|o ∈ o a * , o ̸ = o ab } and {o|o ∈ o * b , o ̸ = o ab }.From Lemma 1.1 and Lemma 1.2 we know that N (Spin) ≥ 3.In this case, Lemma 1 is proven.

2 . 3 .
Lemma 2.1 (a) Claim: We call a set of nonzero matrix elements with one common index and one consecutive index a string, and adding any elements to this set should makes it no longer satisfy the definition of a string.For example, {a 11 = 1, a 12 = 1, a 13 = 1} form a string.There are 2 strings in the set {a 11 = 1, a 12 = 1, a 13 = 0, a 14 = 1}, they are separately {a 11 , a 12 } and {a 14 }.We denote the maximum number of strings that can be possibly defined in a set of matrix elements A by S(A).We claim that if N (Spin) = 3, then S(o * b ) = 1.(b) Proof: Through direct calculation, 1 string correspond to 2 endpoints, flipping 2 spins.So we observe that N (s * ,b+1 ) = 2S(o * b ).(D5) To ensure that N (Spin) = 3, we require that N (s * ,b+1 ) = 2.So it follows that S(o * b ) = 1.Lemma 2.2 (a) Claim:

Fig. 29 .
Fig. 29.A pictorial illustration of main idea in proof in Appendix H.

Table II .
: Typical SPT phases denoted by the new notation system.Here, model I and II are respectively I-MSPT and II-MSPT models, model III is a CSPT model, model IVa, IVb, IVc are all RSPT models, model Va and Vb are both FSPT models.Specially, we can notice that though model IVa, IVb and IVc are all classified as 2D RSPT models, their behavior can be very different.
1) class.This configuration reflects the fractal symmetry of the Sierpinski triangle.At the same time, this configuration possesses the minimal correlation number among all possible strange correlators made up of Pauli matrices in this FSPT model.It can be proved that all 2-point strange correlators show trivial results, giving the same result as the normal correlator gives.Detailed calculation can be found in appendix C. It is natural to ask what is the minimal n that gives nontrivial multi-point strange correlators (giving different results from what multi-point normal correlator gives).It is proved that

. 6 .
Nonlocal Stabilizers in Model-Va Model-Va (Eq.(54), an FSPT model) possesses fractal-like symmetries.Despite the lack of line-like symmetries, there are stringlike nonlocal stabilizers that can contribute to S topo in this model.An example of nonlocal stabilizers generated by the fractal-like symmetry is shown in Fig.21(b), which gives a ∆S topo = −1 in the tripartite topological entropy.The string-like stabilizer can be prolonged in the i direction.
Claim: If N (Spin) = 3, then ∀o ∈ o 1 * , o = 1.(b) Proof: Similarly we can prove o 11 = 1 by considering S(o * 1 ) and S(o 1 * ).If there are any null elements in o 1 * , then it follows that S(o 1 * ) ≥ 2 and N (s 1 * − s 1,b+1 ) ≥ 3, contradicting our assumption N (Spin) = 3.Therefore, Lemma 2, we already have N (s 1 * ∪ s * ,b+1 ) = 3, which means that there are already 3 spins flipped by the Hamiltonian terms.If we want to meet the condition of Theorem 1, the rest of the Spin matrix should not occur any nonzero elements.So the rest part of Op should be selected to ensure that there are no other spins flipped by the Hamiltonian term ∆ ij .

Table VI .
Comparison of several CA approaches in constructing subsystem symmetries.