Fractal Symmetric Phases of Matter

We study spin systems which exhibit symmetries that act on a fractal subset of sites, with fractal structures generated by linear cellular automata. In addition to the trivial symmetric paramagnet and spontaneously symmetry broken phases, we construct additional fractal symmetry protected topological (FSPT) phases via a decorated defect approach. Such phases have edges along which fractal symmetries are realized projectively, leading to a symmetry protected degeneracy along the edge. Isolated excitations above the ground state are symmetry protected fractons, which cannot be moved without breaking the symmetry. In 3D, our construction leads additionally to FSPT phases protected by higher form fractal symmetries and fracton topologically ordered phases enriched by the additional fractal symmetries.


Introduction
Symmetries are indispensable in the characterization and classification of phases of matter. In many cases, knowledge of the systems symmetries and how they are respected or spontaneously broken provide a complete description of a phase. Beyond the usual picture of spontaneously broken symmetries, it has been recently appreciated that multiple phases with the same unbroken symmetry can also exist, known as symmetry protected topological (SPT) phases, which has been the subject of great interest in a variety of systems . These phases cannot be connected adiabatically while maintaining the symmetry, but can be so connected if the symmetry is allowed to be broken. The vast majority of these cases deal with global symmetries, symmetries whose operation acts on an extensive volume of the system. On the opposite end of the spectrum, one also has systems with emergent local gauge symmetries [24], which act on a strictly local finite part of the system; such symmetries may lead to topologically ordered phases [25,26]. In between these two, one also has symmetries which act on some subextensive d-dimensional (integer d) subsystem of a system, such as along planes or lines in a 3D system; these are intermediate between gauge and global symmetries, and have as such also been called gauge-like symmetries. Systems with such symmetries display interesting properties [27][28][29], and are intimately related to models of fracton topological order [30][31][32][33][34][35] through a generalized gauging procedure [36,37] (Type-I fracton order in the classification of Ref 37). Fracton topological order is a novel type of topological order, characterized by subextensive topologydependent ground state degeneracy and immobile quasiparticle excitations, and has inspired much recent research . In a recent work by the present authors, such subsystem symmetries were shown to also lead to new phases of matter protected by the collection of such symmetries [63].
The subject of this paper is yet another type of symmetry, which may be thought of as being "in between" two of the aforementioned subsystem symmetries: fractal symmetries. These act on a subset of sites whose volume scales with linear size L as L d with some fractal dimension d that is in general not integer. Note that these models have symmetries which act on a fractal subset of a regular lattice, and should be distinguished from models (with possibly global symmetries) on fractal lattices [64][65][66][67][68][69]. Systems with such symmetries appear most notably in the context of glassiness in translationally invariant systems [34], such as the Newman-Moore model [70][71][72][73][74][75][76][77][78]. Via the gauge duality [36,37], systems with such symmetries may describe theories with (Type-II) fracton topological order [32,33]; these have ground state degeneracies on a 3-torus that are complicated functions of the system size, and immobile fracton excitations which appear at corners of fractal operators. Indeed, the recent excitement in the study of fracton phases arose from the discovery of the Type-II fracton phase exemplified by Haah's cubic code [33].
We focus on a class of fractal structures on the lattice that are generated by cellular automata (CA), from which many rich structures emerge [79][80][81][82][83]. In particular, we will focus on CA with linear update rules, from which self-similar fractal structures are guaranteed to emerge, following Ref 32. We construct a number of spin models which are symmetric under operations that involve flipping spins along these fractal structures. Unlike a global symmetry, the order of the total symmetry group may scale exponentially with system size, and therefore their case is more like that of subsystem symmetries.
We first present in Sec 2 a brief introduction to CA, and how fractal structures emerge naturally from them. When dealing with such fractals, a polynomial representation makes dealing with seemingly complicated fractal structures effortlessly tractable (see Ref [32]), and we encourage readers to become familiar with the notation. In Sec 3, we take these fractal structures to define symmetries on a lattice in 2D. These symmetries are most naturally defined on a semi-infinite lattice; here, symmetries flip spins along fractal structures (e.g. translations of the Sierpinski triangle). We describe in detail how such symmetries should be defined on various other lattice topologies, including the infinite plane. Simple Ising models obeying these symmetries are constructed in Sec 4, which demonstrate a spontaneously fractal symmetry broken phase at zero temperature, and undergoes a quantum phase transition to a trivial paramagnetic phase.
In Sec 5 we use a decorated defect approach to construct fractal SPT (FSPT) phases. The nontriviality of these phases are probed by symmetry twisting experiments and the existence of symmetry protected ungappable degeneracies along the edge, due to a locally projective representation of the symmetries. Such phases have symmetry protected fracton excitations that are immobile and cannot be moved without breaking the symmetries or creating additional excitations.
Finally, we discuss 3D extensions in Sec 6, these include models similar to the 2D models discussed earlier, but also novel FSPT phases protected by a combination of regular fractal symmetries and a set of symmetries which are analogous to higher form fractal symmetries. These FSPT models with higher form fractal symmetries, in one limit, transition into a fracton topologically ordered phase while still maintaining the fractal symmetry. Such a phase describes a topologically ordered phase enriched by the fractal symmetry, thus resulting in a fractal symmetry enriched (fracton) topologically ordered (fractal SET [84][85][86][87][88][89][90][91][92], or FSET) phase.

Cellular Automata Generate Fractals
We first set the stage with a brief introduction to a class of one-dimensional CA, from which it is well known that a wide variety of self-similar fractal structures emerge. In latter sections, these fractal structures will define symmetries which we will demand of Hamiltonians.
Consider sites along a one-dimensional chain or ring, each site i associated with a p-state variable a i ∈ {0, 1, . . . , p − 1} taken to be elements of the finite field F p . We define the state of the CA at time t as the set of a (t) i . We will typically take p = 2, although our discussion may be easily generalized to higher primes. We consider CA defined by a set of translationally-invariant local linear update rules which determine the state {a (t+1) i } given the state at the previous time {a (t) i }. Linearity here means that the future state of the ith cell, a (t+1) i , may be written as a linear sum of a (t) j for j within some small local neighborhood of i. Throughout this paper, all such arithmetic is integer arithmetic modulo p, following the algebraic structure of F p . Figure 1 shows two sets of linear rules which we will often refer to: 1. The Sierpinski rule, given by a with p = 2, so called because starting from the state a (0) i = δ i,0 one obtains Pascal's triangle modulo 2, who's nonzero elements generate the Sierpinski triangle with fractal Hausdorff dimension d = ln 3/ ln 2 ≈ 1.58. In the polynomial representation (to be introduced shortly), this rule is given by f (x) = 1 + x.
Fractal dimensions for CA with linear update rules may be computed efficiently [95].
To see why such linear update rules always generate self-similar structures, it is convenient to pass to a polynomial representation. We may represent the state a In the polynomial representation, the row t is given by Notice that self-similarity at every row t = 2 l (here, we show evolution up to t = 40).
for an infinite chain. Alternatively, periodic boundary conditions may be enforced by setting x L = 1.
In this language, these update rules take the form for some polynomial f (x) containing only small finite powers (both positive or negative) of x. For the Sierpinski rule we have f (x) = 1 + x, and for the Fibonacci rule we have f (x) = x −1 + 1 + x. Then, given an initial state s 0 (x), we have that A neat fact about polynomials in F p is that they obey what is known as the Freshman's Dream, whenever t is a power of p. This can be shown straightforwardly by noting that the binomial coefficient p k n is always divisible by p unless n = 0 or n = p k . It thus follows that such CA generate fractal structures. Let us illustrate for the Fibonacci rule starting from the initial configuration s 0 (x) = 1, i.e. the state where all a i = 0 except for a 0 = 1. Looking at time t = 2 l , the state is s t (x) = x −2 l + 1 + x 2 l . In the following evolution, each of the non-zero cells a −2 l = a 0 = a 2 l = 1 each look locally like the initial configuration s 0 , and thus the consequent evolution results in three shifted structures identical to the initial evolution of s 0 (up until they interfere), as can be seen in Figure 1. At time t = 2 k+1 , this process repeats but at a larger scale. Thus, we can see that any linear update rule of this kind will result in self-similar fractal structures when given the initial state s 0 (x) = 1. As the rules are linear, all valid configurations correspond to superpositions of this shifted fractal.
The entire time evolution of the CA may be described at once by a single polynomial F (x, y) over two variables x and y, and we have that the coefficient of The two-dimensional fractal structures in Figure 1 generated by these CA emerge naturally due to a set of simple local constraints given by the update rules. In the next section, we will describe 2D classical spin Hamiltonians which energetically enforce these local constraints. The ground state manifold of these classical models is described exactly by a valid CA evolution, which we will then take to define symmetries.

Fractal Symmetries
To discuss physical spin Hamiltonians and symmetries, it is useful to also use a polynomial representation of operators. Such polynomial representations are commonly used in classical coding theory [96], and refined in the context of translationally invariant commuting projector Hamiltonians by Haah [97]. We will utilize only the basic tools (following much of Ref [32]), and specialize to Pauli operators (p = 2 from the previous discussion), although a generalization to p-state Potts spins is straightforward.
Let us consider in 2D a square lattice with one qubit (spin-1/2) degree of freedom per unit cell. Acting on the qubit at site (i, j) ∈ Z 2 , we have the three anticommuting Pauli matricesẐ ij ,X ij , and Y ij . We define the function Z(·) from polynomials in x and y over F 2 to products of Pauli operators, such that acting on an arbitrary polynomial we have and similarly for X(·) and Y (·). For example, we have Z(1 + x + xy) = Z 0,0 Z 1,0 Z 1,1 . Some useful properties are that the product of two operators is given by the sum of the two polynomials, Z(α)Z(β) = Z(α + β), and a translation of Z(α) by (i, j) is given by Z(x i y j α).
Perhaps the most useful property of this notation is that two operators Z(α) and X(β) anticommute if and only if [αβ] x 0 y 0 = 1, where [·] x i y j denotes the coefficient of x i y j in the polynomial, and we have introduced the dual, which may be thought of as the spatial inversion about the point (0, 0). We will also often usex to represent x −1 for convenience. More usefully, we may express the commutation relation between Z(α) and translations of X(β) (given by X(x i y j β)) as where d ij may be computed directly from the commutation polynomial of α and β, which may easily be computed directly given α and β. In particular, P = 0 would imply that every possible translations of the two operators commute.

Semi-infinite plane
We may now transfer our discussion of the previous section here. Let us consider a semi-infinite plane, such that we only have sites (i, j) with x i y j≥0 . We may then interpret the jth row as the state of a CA at time j, starting from some initial state at row j = 0. Consider the linear CA with update rule given by the polynomial f (x), as defined in Eq 2. The classical Hamiltonian which energetically enforces the CA's update rules is given by where we have excluded terms that aren't fully inside the system. As an example, consider the Sierpinski rule f = 1 + x (f will always refer to a polynomial in only x). Equation 10 for this rule gives, (11) which is exactly the Newman-Moore (NM) model originally of interest due to being an exactly-solvable translationally invariant model with glassy relaxation dynamics [70]. The NM model was originally described in a more natural way on the triangular lattice as the sum of three-body interactions on all downwards facing triangles, H NM = − ▽ ZZZ. This model does not exhibit a thermodynamic phase transition (similar to the 1D Ising chain). Fractal codes based on higher-spin generalizations of this model have also been shown to saturate the theoretical information storage limit asymptotically [98].
We will be interested in the symmetries of such a model that involve flipping subsets of spins. Due to the deterministic nature of the CA, such operation must involve flipping some subset of spins on the first row, along with an appropriate set of spins on other rows such that the total configuration remains a valid CA evolution. Operationally, symmetry operations are given by various combinations of F (x, y) (Eq 5). That is, for any polynomial q(x), we have a symmetry element Here, q(x) has the interpretation of being an initial state s 0 , and S(q(x)) flips spins on all the sites corresponding to the time evolution of s 0 . As the update rules are linear, this operation always flips between valid CA evolutions. For example, S(1) will correspond to flipping spins along the fractals shown in Fig 1. To confirm that this symmetry indeed commutes with the Hamiltonian, we may use the previously discussed technology (Eq 8 and 9) to compute the commutation polynomial between S(q(x)) and translations of the Hamiltonian term Z(1 +fȳ), Since terms which have shift y 0 are not included in the Hamiltonian (Eq 10), this operator therefore fully commutes with the Hamiltonian. We may pick as a basis set of independent symmetry elements, S(q α ), for α ∈ Z with q α (x) = x α . These operators correspond to flipping spins corresponding to the colored pixels in Fig 1, and horizontal shifts thereof. Each of these symmetries act on a fractal subset of sites, with volume scaling as the Hausdorff dimension of the resulting fractal.

Cylinder
Rather than a semi-infinite plane, let's consider making the x direction periodic with period L, such that x L = 1, while the y direction is either semi-infinite or finite. In this case, there are a few interesting possibilities.

Reversible case
In the case that there exists some ℓ such that f ℓ = 1, then the CA is reversible. That is, for each state s t , there exists a unique state s t−1 such that s t = f s t−1 , given by s t−1 = f ℓ−1 s t . A proof of this is straightforward, suppose there exists two distinct previous states s t−1 , s ′ t−1 , such that f s t−1 = f s ′ t−1 = s t . As they are distinct, s t−1 + s ′ t−1 = 0. However, there is a contradiction. Hence, the state s t−1 must be unique. The inverse statement, that a reversible CA must have some ℓ such that f ℓ = 1, is also true.
In this case, all non-trivial symmetries extend throughout the cylinder, and their patterns are periodic in space with period dividing ℓ. An example of this is the Fibonacci rule with L = 2 m , for which f L/2 = 1. There are L independent symmetries on either the infinite or semi-infinite cylinder, and the total symmetry group is simply Z L 2 . The symmetries on an infinite cylinder are given by

Trivial case
If there exists ℓ such that f ℓ = 0, then the model is effectively trivial. All initial states s 0 will eventually flow to the trivial state s ℓ = 0. On a semi-infinite cylinder, possible "symmetries" will involve sites at the edge of the cylinder, but will not extend past ℓ into the bulk of the cylinder. On an infinite cylinder, there are no symmetries at all. An example of this is the Sierpinski rule with L = 2 m , for which f L = 0.

Neither reversible nor trivial
If the CA on a cylinder is neither reversible nor trivial, then every initial state s 0 must eventually evolve into some periodic pattern, such that s t = s t+T for some period T at large enough t (this follows from the fact that there are only finitely many states). Thus, there will be symmetry elements whose action extends throughout the cylinder, like in the reversible case. Interestingly, however, irreversibility also implies the existence of symmetry elements whose action is restricted only to the edge of the cylinder, much like the trivial case. Let us take two distinct initial states s 0 , s ′ 0 that eventually converge on to the same state at time ℓ. Then, lets 0 = s 0 + s ′ 0 = 0 be another starting state. After time ℓ,s ℓ = s ℓ + s ′ ℓ = 0, this state will have converged on to the trivial state. Thus, the symmetry element corresponding to the starting states 0 will be restricted only to within a distance ℓ of the edge on a semi-infinite cylinder.
On an infinite cylinder, only the purely periodic symmetries will be allowed, so the total number of independent symmetry generators is reduced to between 0 and L.

On a torus
Let us next consider the case of an L x × L y torus. Symmetries on a torus must take the form of valid CA cycles on a ring of length L x with period L y . The order of the total symmetry group is the total number of distinct cycles commensurate with the torus size, which in general does not admit a nice closed-form solution, but has been studied in Ref 94. Equivalently, there is a one-to-one correspondence between elements in the symmetry group and solutions to the equation with x Lx = 1. This may be expressed as a system of linear equations over F 2 , and can be solved efficiently using Gaussian elimination. For each solution q(x) of the above equation, the action of the corresponding symmetry element is given by As an example, consider the Sierpinski model on an L × L torus. Let k(L) = log 2 (N sym (L)) be the number of independent symmetry generators, where N sym (L) is the order of the symmetry group. We are free to pick some set of k(L) independent symmetry operators as a basis set (there is no most natural choice for basis), which we label by q α (x) with 0 ≤ α < k. To illustrate that k(L) is in general a complicated function of L, we show in Table 2 k(L) and a choice of q where the number of cycles can be solved for exactly. An interesting point is that for the Sierpinski rule, f (x) 2 l = 0, thus for L = 2 l , there are no non-trivial solutions to Eq 15 and so k(2 l ) = 0. To contrast, the Fibonacci rule has f (x) 2 l = 1, and so k(L = 2 l ) = L.

Infinite plane
Now, let us consider defining such symmetries directly on an infinite plane, where we allow all x i y j .
In the CA language, we are still free to pick the CA state at time, say t = 0, s 0 (x), which completely determines the CA states at times t > 0. However, we run into the issue of reversibility -how do we determine the history of the CA for times t < 0 which lead up to s 0 ? For general CA, there may be zero or multiple states s −1 which lead to the same final state s 0 . For a linear CA on an infinite plane, however, there is always at least one s −1 . We give an algorithm for picking out a particular history for s 0 , and discuss the sense in which it is a complete description of all possible symmetries, despite this reversibility issue. For this discussion it is convenient to, without loss of generality, assume f contains at least a positive power of x (we may always perform a coordinate transformation to get f into such a form). The basic idea is as follows: we have written the Hamiltonian (Eq 10) in a form that explicitly picks out a direction (y) to be interpreted as the time direction of the CA. However, we may always write the same term as a higher-order linear CA that propagates in thex direction, where a > 0 is the highest power of x in f , n max is finite, and g n (y) is a polynomial containing only non-negative powers of y. This describes an n max -order linear CA. For the Sierpinski rule, we have only g 1 (y) = 1 + y, and for the Fibonacci rule we have both g 1 (y) = 1 + y and g 2 (y) = 1. We then further define g(x, y) for convenience, which only contains non-negative powers of x and y. Now, consider the fractal pattern generated byx which describes a higher-order CA evolving in thex direction. Note that powers ofḡ no longer have the nice interpretation of representing an equal time state in terms of this CA, due to it containing both powers ofȳ as well asx (but evaluating the series up to theḡ nxn does give the correct configuration up tox n ). Asḡ contains only negative powers of y, this fractal pattern is restricted only to the half-plane with y j<0 . It thus lives entirely in the "past", t < 0, of our initial CA. The full fractal given by unambiguously describes a history of the CA with the t = 0 state s 0 = 1. This is shown in Figure 3 for the Fibonacci model, with the forward propagation of f in red and the propagation ofḡ in orange. Going back to operator language, it can be shown straightforwardly that the symmetry action for arbitrary q(x) commutes with the Hamiltonian (Eq 10 but with all ij included in the sum) everywhere. The only term with y 0 in F is 1, so this operator only flips the spins q(x) on row y 0 . Furthermore, the choice of choosing the y 0 row for defining this symmetry does not affect which operators can be generated, as it is easy to show that f (x)F (x, y) =ȳF (x, y), so that S(q(x)f (x)) flips any set of spins q(x) on the rowȳ instead. Simple counting would then suggest that the total number of symmetry generators thus scales linearly with the size of the system, like on the semi-infinite cylinder. This result seems to contradict the irreversibility of the CA. It would suggest that one can fully determine s t at time t < 0 by choosing the state s 0 appropriately, which would seemingly imply that the evolution is always reversible. The resolution to this paradox lies in the fact that we are on an infinite lattice, and in this procedure we have chosen the particular f −1 such that it only contains finitely positive powers of x (there are in general multiple inverses f −1 ). Defining h(x) = [g(x, y)] y 0 such that f = x a (1 +hx), then we are choosing the inverse from which it can be readily verified that f −1 f = 1. In this language, F (x, y) looks like which obviously commutes with the Hamiltonian. As an example, with the Sierpinski rule, the two possible histories for the state s 0 = 1 are s By this inverse, we would only get s (−) −1 . However, if we wanted to generate the state with history s (+) −1 , we would instead find that the t = 0 state should be the limit s 0 = 1 + x ∞ . If we were just interested in any finite portion of the infinite lattice, for example, we may get any history by simply pushing this x ∞ beyond the boundaries.

Open slab
Finally, consider the system on an open slab with dimensions L x ×L y . Elements of the symmetry group are in correspondence with valid CA configurations on this geometry. The state at time t = 0 may be chosen arbitrarily, giving us L x degrees of freedom. Furthermore, at each time step the state of the cells near the edge may not be fully specified by the CA rules. Hence, each of these adds an additional degree of freedom. Let x −pmin , x pmax , be the smallest and largest powers of x in f (if p min/max would be negative, then set set it to 0). Then, we are free to choose the cell states in a band p max × L y along the left (x i=0 ) edge, and p min × L y along the right edge as well. Thus, the total number of choices is N sym = 2 Lx+(pmin+pmax)(Ly−1) , and there are log 2 N sym independent symmetry generators. Note that some of these symmetries may be localized to the corners.
One may be tempted to pick a certain boundary condition for the CA, for example, by taking the state of cells outside to be 0, which eliminates the freedom to choose spin states along the edge and reduces the order of the symmetry group down to simply 2 Lx . What will happen in this case is that there will be symmetry elements from the full infinite lattice symmetry group which, when restricted to an L x × L y slab, will not look like any of these 2 Lx symmetries. With the first choice, we are guaranteed that any symmetry of the infinite lattice, restricted to this slab, will look like one of our N sym symmetries. This is a far more natural definition, and will be important in our future discussion of edge modes in Sec 5.3.

Spontaneous fractal symmetry breaking
At T = 0, the ground state of H classical is 2 k -degenerate and spontaneously breaks the fractal symmetries, where k is the number of independent symmetry generators (which will depend on system size and choice of boundary conditions). Note that k will scale at most linearly with system size, so it represents a subextensive contribution of the thermodynamic entropy at T = 0. As a diagnosis for long range order, one has the many-body correlation function C(ℓ) given by which has C(ℓ) = 1 in the ground states of H classical as can be seen by the fact that Eq 23 is a product of terms in the Hamiltonian. If M is the number of terms in f , then this becomes an M + 1-body correlation function when ℓ = 2 l is a power of 2. Long range order is diagnosed by lim ℓ→∞ C(ℓ) = const. At any finite temperature, however, these models are disordered and have C(ℓ) vanishing asymptotically as C(ℓ) ∼ p −ℓ d , where d is the Hausdorff dimension of the generated fractal, and p = 1/(1 + e −2β ). This can be seen by mapping to the dual (defect) variables in which the Hamiltonian takes the form of a simple non-interacting paramagnet [70], and the correlation function C(ℓ) maps on to a O(ℓ d )-body correlation function. Thus, there is no thermodynamic phase transition in any of these models, although the correlation length defined through C(ℓ) diverges as T → 0. Even without a thermodynamic phase transition, much like in the standard Ising chain, there is the possibility of a quantum phase transition at T = 0. We may include quantum fluctuations via the addition of a transverse field h, One can confirm that a small h will indeed correspond to a finite correction lim l→∞ C(2 l ) = 1 − const(h), and so does not destroy long range order. This model now exhibits a zero-temperature quantum phase transition at h = 1, which is exactly pinpointed by a Kramers-Wannier type selfduality transformation which exchanges the strong and weak-coupling limits. This self-duality is readily apparent by examining the model in terms of defect variables, which interchanges the role of the coupling and field terms. This should be viewed in exact analogy with the 1D Ising chain, which similarly exhibits a T = 0 quantum phase transition but fails to have a thermodynamic phase transition. The transition at h = 1 is a spontaneous symmetry breaking transition in which all 2 k fractal symmetries are spontaneously broken at once (although under general perturbations they do not have to all be broken at the same time). Numerical evidence [69] suggests a first order transition. If one were to allow explicitly fractal symmetry breaking terms in the Hamiltonian (Z-fields, for example) then it is possible to go between these two phases adiabatically. Thus, as long as the fractal symmetries are not explicitly broken in the Hamiltonian, these two phases are properly distinct in the usual picture of spontaneously broken symmetries. In the following, we will only be discussing ground state (T = 0) physics.

Fractal symmetry protected topological phases
Rather than the trivial paramagnet and spontaneously symmetry broken phases, we may also generate cluster states [99] which are symmetric yet distinct from the trivial paramagnetic phase. These cluster states have the interpretation of being "decorated defect" states, in the spirit of Ref 100, as we will demonstrate. These fractal symmetry protected topological phases (FSPT) are similar to recently introduced subsystem SPTs [63], and were hinted at in Ref 36. In contrast to the subsystem SPTs, however, there is nothing here analogous to a "global" symmetry -the fractal symmetries are the only ones present!

Decorated defect construction
To describe these cluster Hamiltonians, we require a two-site unit cell, which we will refer to as sublattice a and b. For the unit cell (i, j) we have two sets of Pauli operatorsẐ Our previous polynomial representation is extended as and similarly for X(·) and Y (·). This notation is easily generalized to n spins per unit cell, represented by n component vectors. Our cluster FSPT Hamiltonian is then given by which consists of commuting terms and is exactly solvable at h = h x = h z = 0, which we will assume for now. There is a unique ground state on a torus (regardless of the symmetries). The ground state is short range entangled, and may be completely disentangled by applications of controlled-X (CX) gates at every bond between two different-sublattice sites that share an interaction, as per the usual cluster states -however, this transformation does not respect the fractal symmetries of this model. These fractal symmetries come in two flavors, one for each sublattice: where we have assumed an infinite plane with F (x, y) as in Eq 22, and q(x) may be any polynomial. The picture of the ground state is as follows. Working in theẐ (a) ,Ẑ (b) basis, notice that ifẐ type symmetry (Eq 27). In (c), we perform a symmetry twist on the Sierpinski FSPT on a 7 × 7 torus. The chosen symmetries g 1 (g 2 ) corresponds to operations on all spins highlighted by a large blue (red) circle. The green triangles correspond to terms in the twisted Hamiltonian H twist (g 1 ) that have flipped sign. The charge response T (g 1 , g 2 ) = −1 is given by the parity of red circles that also lie in the green triangles, and is independent of where we make the cut j 0 .
CA rule is followed. The second term in the Hamiltonian transitions between states with different configurations of such defects. The ground state is therefore an equal amplitude superposition of all possible configurations. The same picture can also be obtained from theX (a) ,X (b) basis, in terms of the CA rules acting on theX (b) ij spins.

Sierpinski FSPT
As a particularly illustrative example, let us consider the FSPT generated from the Sierpinski rule. The resulting model is the "decorated defect" NM paramagnet, which we refer to as the Sierpinski FSPT. The Hamiltonian is given by It is particularly enlightening to place this model on a honeycomb lattice, as shown in Fig 4a. Fig 4b  shows the action of two symmetries as an example.
We may then redefineẐ ij , after which the Hamiltonian takes the particularly simple form of a cluster model where s = (i, j, a/b) labels a site on the honeycomb lattice and Γ(s) is the set of its nearest neighbors. However, we will generally not use such a representation. Note that this model is isomorphic to the 2D fractal SPT obtained in Ref [107].

Fibonacci FSPT
Our other example is the Fibonacci FSPT. The Hamiltonian takes the form which we illustrate in Fig 5a. Unlike with the Sierpinski FSPT, this model does not have as nice of an interpretation of being a cluster model with interactions among sets of nearest neighbors on some simple 2D lattice.

Symmetry Twisting
To probe the nontriviality of the FSPT symmetric ground state, we may place it on a torus and apply a symmetry twist to the Hamiltonian, and observe the effect in the charge of another symmetry [102][103][104][105]. To be concrete, let H twist (g 1 ) be the g 1 symmetry twisted Hamiltonian. The g 2 charge of the ground state of H twist (g 1 ) relative to its original value tells us about the nontriviality of the phase under these symmetries. That is, let with Z the partition function, then, we define the charge response where g 2 1 is simply the g 2 charge of the ground state of the untwisted Hamiltonian. On a torus, we may twist along either the horizontal or vertical direction -here we first consider twisting along the vertical direction. Let us be more concrete. Take the FSPT Hamiltonian (Eq 26) on an L x × L y torus, and let k be the number of independent symmetry generators of the type Z 2 ). We assume L x , L y have been chosen such that k > 0. The total symmetry group of our Hamiltonian is therefore Z where 0 ≤ α < k and q To apply a g-twist, we first express the Hamiltonian as a sum of local terms H FSPT = ij H ij . We then pick a horizontal cut j = j 0 , dividing the system between j < j 0 and j ≥ j 0 . For each term that crosses the cut, we conjugate H ij → g < H ij g −1 < , where g < is the symmetry action of g restricted to j < j 0 . For an Ising system, this will simply have the effect of flipping the sign of some terms in the Hamiltonian. The resulting Hamiltonian is H twist (g).
To understand which terms in the Hamiltonian change sign under conjugation, consider the choice of symmetry g 1 in Fig 4c, which consists of flipping all spins in the large blue (dark and transparent) circles. Restricting g 1 to j < j 0 leaves g 1,< , flipping only spins in the dark circles. Conjugating by g 1,< results in the terms in the green triangles appearing in H twist (g 1 ) with a relative minus sign.
Doing this explicitly for a symmetry element S (a) α , we find that the incomplete symmetry restricted to j < j 0 is given by The terms in the Hamiltonian that pick up a minus sign when conjugated with S (a) α,< are exactly translations of the first term in H FSPT (Eq 26) given by the non-zero coefficients of the commutation polynomial along j 0 : P = q (a) α (x)(f y) j0 . However, the same twisted Hamiltonian may also be obtained by conjugating the entire H FSPT by (36) such that H twist (S where we have used y Ly = 1 and the definition of a symmetry on the torus, Eq 15. As expected, the result is independent of our choice of j 0 , and it is also apparent that T (g 1 , g 2 ) = T (g 2 , g 1 ) for any g 1 ,g 2 . If we choose the same symmetry basis for both sublattices, q α (x), then we additionally get that T αβ = T βα . Figure 4c is an illustration of this twisting calculation for the Sierpinski FSPT on a 7 × 7 torus. Letting x 0 y 0 label the unit cell in the top left of the figure, g 1 is an (a) type symmetry with q (a) (x) = x 3 + x 4 and g 2 is a (b) type symmetry with q (b) (x) = x 4 + x 5 . Then, Eq 37 gives T (g 1 , g 2 ) = −1, which can be confirmed by eye in the figure.
The exact same procedure may also be applied for twists across the horizontal direction, which will provide yet another set of independent relations between the symmetries (but will not have as nice of a form).

Degenerate edge modes
Upon opening boundaries, the ground state manifold becomes massively degenerate. Away from a corner, we will show that these degeneracies cannot be broken by local perturbations as long as the fractal symmetries are all respected, much like in the case of SPTs with one-dimensional subsystem symmetries [63].
Let us review the open slab geometry from Sec 3.5 for the FSPT. We take the system to be a rectangle with L x × L y unit cells, such that we are restricted to x 0≤i<Lx y 0≤j<Ly . as before, let x −pmin , x pmax , be the smallest and largest powers of x in f (and let p min/max = 0 if they would be negative).
The total symmetry group is Z and we assume L x > R (otherwise there are no allowed terms in the Hamiltonian at all). A Z (a) 2 type symmetry acts as X (a) ij on a subset of unit cells, and a particular symmetry is fully specified by how it acts on the top row x i y 0 , the band x i<pmax y j (on the left side), and the band x i≥Lx−pmin y j (on the right side). A Z (b) 2 type symmetry acts as Ẑ (b) ij and a particular one is fully specified in a similar manner, but spatially inverted (top↔bottom, left↔right). Alternatively, we may simply think of the symmetries as those of the infinite plane, but truncated to the L x × L y slab.
On the open slab, we take our Hamiltonian (Eq 26) with h = 0 on the infinite plane and simply exclude terms that contain sites outside of the sample. For each term with shift x i y j that are excluded, but for which the unit cell x i y j is still in the system, we lose a constraint on the ground state manifold and hence gain a two-fold degeneracy. The number of terms excluded is given by exactly the same counting as before. Along the top (bottom) edge, there is one excluded Z (X) term per unit cell. Along the left edge, there are p max Z terms excluded and p min X terms, for a total of R excluded terms per unit cell, and similarly for the right edge. Hence, there are a total of 2 2k ground states, coming from a 2 R -fold degeneracy per unit cell along the left/right edges, and 2-fold degeneracy per unit cell along the top/bottom (with some correction for overcounting).
To describe the edge physics on the slab geometry, let us introduce some additional notation. Let us denote the truncation of an arbitrary polynomial p(x, y) to the slab as [p(x, y)] slab , where only the terms with x i y j with (i, j) ∈ slab are kept, where is the set of sites (i, j) which exist on the L x × L y slab, and [a, b] = {a, a + 1, . . . b}. We may further make the distinction between those on the edge or the bulk of the slab. Let us denote two types of bulks, which we denote by bulk a and bulk b , such that the Hamiltonian on the slab is given by Finally, we denote the edge simply as those sites in the slab that are not in the bulk, For each excluded Z term in H slab , i.e. each (i, j) ∈ edge a , we may define a set of three Pauli operators,X and similarly, for each excluded X term at (i, j) ∈ edge b , we may definê where the [·] slab truncation ensures that only those sites physically in the slab are involved. We will call such operators "edge" Pauli operators. There are 2k such sets of edge Pauli operators, one for each excluded term. It may readily be verified thatX satisfy the Pauli algebra while being independent of and commuting with every term in the Hamiltonian and each other at different sites. They therefore form a Pauli basis for operators which act purely within the 2 2k dimensional ground state manifold.
In principle, any local perturbation, projected on to the ground state manifold, will have the form of being some local effective Hamiltonian in terms of these edge Pauli operators, and may break the exact degeneracy. However, we wish to consider only perturbations commuting with all fractal symmetries. To deduce what type of edge Hamiltonian is allowed, we must find out how our many symmetry elements act in terms of these edge operators.
Consider the action of a Z (a) 2 symmetry on the slab, which is written as the truncation of a symmetry on an infinite plane to a slab, as discussed in Sec 3.5. By construction, we have that F (x, y)(1 + f y) = 0, so we may also write Let us denote for convenience γ ≡ q(x)F (x, y), which can be decomposed into three parts: slab c , the bulk b part, the edge b part, and the parts external to the slab (denoted by the complement slab c ). Then, we may be decompose the symmetry action as The first factor, the bulk action, is made out of products of terms in H slab (i.e. is an element of the stabilizer group) and therefore acts trivially on the ground state manifold. The second factor acts only on edge b , and operates within the ground state manifold as a product ofX In a similar fashion, we may show that a Z We claim that it is always possible to find a particular symmetry element that acts locally on one edge in any way (but it will generally extend non-trivially into the bulk and act in complicated way on the other boundaries). For example, for any (i 0 , j 0 ) on the left edge, there exists a Z There is no non-trivial operator acting on a single edge that commutes with bothX andẐ, and therefore we are prohibited from adding anything non-trivial to the effective Hamiltonian on this edge which therefore guarantees that no degeneracy can be broken while respecting all fractal symmetries. Note that we don't even have the possibility of spontaneous symmetry breaking at the surface -even simpleẐẐ couplings along the edge violate the symmetries. The only way the ground state degeneracy may be broken without breaking the symmetry is by terms which couple edge Paulis along different edges; these terms are either non-local, or located at a corner of the system. Suppose we have found a particular Z 2 with generators g 1 , g 2 , would have (g 1 g 2 ) 2 = 1. However, if we look at the action on this particular edge, then we have that (g edge 1 g edge 2 ) 2 = (XẐ) 2 = −1. Since we know that as a whole g 1 and g 2 must commute, the action of g 1 and g 2 on the other edges must again anticommute (to cancel out the −1 from this edge). Small manipulations of the edges (such as adding or removing sites) or local unitary transformations respecting the symmetry cannot change the fact that the actions of g 1 and g 2 are realized projectively on this edge.
Near particular corners, some symmetry elements may act essentially locally. As a symmetry element (as a whole) must commute with all others, nothing prevents the addition of the full symmetry action itself as a term in the effective Hamiltonian when it is local. For example, when h = 0 there will be terms appearing in the effective Hamiltonian at finite order in perturbation theory near such corners, which commute with all symmetries. The magnitude of such terms will decay exponentially away from a corner, however, and therefore we still have an effective degeneracy per unit length along the boundaries.

Local action of symmetries on edges
To prove our claim that there is always a symmetry element which acts locally along an edge, let us first consider finding a particular Z 2 symmetry element acting locally as well then follows. Such a symmetry will act locally in some way on the edge, but extend into the bulk in a non-trivial way. Note that there is no "most natural basis" for these symmetries, unlike in the case of integer d subsystem symmetries [63].
Let us take a general Z (a) 2 symmetry element defined according to Eq 47 in terms of a single polynomial q(x). However, multiple q(x) may lead to the same symmetry element on the slab. Recall that in the CA picture, q(x) corresponds to the CA state at time 0, and the L x ×L y slab is a space-time trajectory for the CA. There is a strictly defined "light-cone" determined by the CA rules for which cells at time 0 can affect a future cell in our L x × L y slab. It is easy to verify that only the coefficients in q(x) of x i for −p max (L y − 1) ≤ i < L x + p min (L y − 1) can affect the way the symmetry acts within the slab. Let us therefore take q(x) to only contain powers of x within this range. Furthermore, we see that there are 2 pmax(Ly−1)+Lx+pmin(Ly−1) = 2 Lx+R(Ly−1) = 2 k possible q(x)s, which is also the order of the (Z (a) 2 ) k symmetry group, which means that there is a one-to-one correspondence between q(x)s and different elements of the symmetry group.
Top edge Finding a symmetry element that acts locally on the top edge is simple. The only possibilities on the top edge are for it to act asX (a) i,0 operators. For example, we may simply choose any q(x) such that [q(x)] slab = x i0 , and the corresponding symmetry element will act locally as onlyX There is always such a q(x) that does this, as we showed for the infinite plane (Sec 3.4) that one can always find a history for any CA state).
Left/right edge Along the left/right edges, things are slightly trickier. Let us look at only the left edge for now. A symmetry element may act asX (a) ij for 0 ≤ j < p max , or asX ij for 0 ≤ j < p min . Per unit cell along the left edge, there are 2 pmin+pmax = 2 R possible ways to act. From Eq 50, we see that the non-zero coefficients of q(x)F (x, y) in the columns 0 ≤ j < p min of the slab directly correspond to how the symmetry element acts asX (a) ij on the left edge. Once these have been fixed, the coefficients on the p min ≤ j < R columns must be chosen to specify how the symmetry element acts (asX ij ) on the left edge. Thus, to find a particular symmetry element that will act in a particular way on the left edge, we must specify the leftmost R columns of q(x)F (x, y). By a similar lightcone argument as before, these R columns are affected by coefficients of x i in q(x) with −p max (L y − 1) ≤ i < R + p min (L y − 1). As there are a total of 2 RLy possible histories, and also 2 RLy cells within the leftmost R columns, we may fully specify the action of the symmetry within these leftmost R columns by an appropriate choice of q(x). The remaining degrees of freedom in q(x) means that there are a total of 2 k−RLy = 2 Lx−R symmetry elements that act in the same way on the left edge. ij on the left edge, for the Fibonacci FSPT (Eq 31), whose terms are shown in Fig 5(left). The freedom to choose how the symmetry acts on the right edge (question marks) exactly corresponds to the 2 Lx−R symmetry elements with the specified action on the left edge. These form a Z Lx−R 2 subgroup of the total symmetry group. We choose to show the Fibonacci FSPT here rather than the Sierpinski FSPT, as the latter has R = 1 and is straightforward.

Excitations
On the infinite plane, the lowest lying excitations are strictly immobile. They are therefore fractons protected by the total fractal symmetry group.
Take h = 0, the lowest lying excited states consist of excitations of a single term in the Hamiltonian, say the Z term at site x 0 y 0 . This excited state can be obtained by acting on the ground state witĥ X (b) 0,0 . One may alternatively think in terms of symmetries. Take an independent set of symmetry generators g (a/b) α of the form Eq 27 with the basis choice q (a/b) α = x α . We find that this excited state is uncharged, g (a/b) α = 1, with respect to all symmetry elements except g In fact, the only state with a single excitation with g (b) 0 = −1 is this one with the excitation at the origin.
Let us consider the block of the Hamiltonian with symmetry charges g (b) α = (−1) dα . The blocks containing states with single fractons will have for which the excitation is strictly localized at site x i y j . The excitation may move away from x i y j , but at the cost of creating additional excitations as well, such that the charge of all symmetries are unchanged. If one allows breaking of the fractal symmetries, then these charges are no longer conserved and nothing prevents the excitation from moving to a different site.
On lattices with different topology, these fractons may not be strictly immobile. For example, on a torus, depending on the symmetries, a fracton may be able to move to some subset of other sites (or all other sites, if the symmetry group is trivial). However, such hopping terms are exponentially suppressed with system size. In fact, for the Sierpinski FSPT on a torus with no symmetries, it is actually easier perturbatively to hop a fracton a large power of 2 away than it is to hop a short distance (mimicking some form of p-adic geometry with p = 2).
On an open slab, the ground state manifold is degenerate and all charge assignments are possible in the ground state, protected by the symmetry. Therefore, a fracton may be created, or moved, in any way. However, the amplitude for doing so will decay exponentially away from the edges, and certain processes may only be possible near certain types of edges or corners. The possibilities will depend on the details of the model.

Duality
Here we outline a duality that exist generally for these models, which maps the FSPT phase to two copies of the spontaneous symmetry broken phase of the quantum Hamiltonian in Sec 4. This duality involves non-local transformations and maps the 2 2k ground states of the FSPT on the open slab to the 2 2k symmetry breaking ground states of the dual model. This duality is most naturally described on an L x × L y cylinder (with x Lx = 1) or slab. Let us define new Pauli operatorsZ(·) andX(·) as and translations thereof. It can be readily verified that the latter two commute, and as a whole the set of these operators satisfy the correct Pauli algebra. The fractal symmetries only involve operators in line 52, and so are unchanged. The interaction terms are modified however: in terms of these operators, we haveZ and so the Hamiltonian H FSPT (Eq 26) becomes two decoupled copies of H Quantum (Eq 24) with their own set of symmetries. From this, it follows that the order parameter measuring long-range order in H Quantum , C(ℓ) (Eq 23), maps on to a fractal order parameter in our original basis which is pictorially shown for the Sierpinski FSPT in Figure 6, and approaches a constant in the FSPT phase, or zero in the trivial paramagnet, as ℓ = 2 l → ∞. By the self-duality of H Quantum , we also know the FSPT to trivial transition happens at exactly h = 1.
Finally, this duality allows us to determine the full phase diagram even as h x = h z . Keeping h x small and making h z large, one of the H Quantum is driven into its paramagnetic phase where spins are polarized asẐ The Hamiltonian H FSPT then looks like a single H Quantum , and therefore has spontaneously symmetry broken ground states. By the duality transformation, we know this transition happens at exactly h z = 1. The phase diagram is summarized in Fig 7(left).

Three dimensions
Here, we briefly examine the possible physics available in higher dimension. We consider our symmetrydefining CA in 3D in two ways: via one 2D CA, or two 1D CA. The first will have similar properties to our earlier models, while the latter in certain limits also lead to exotic fractal spin liquids introduced by Yoshida [32] and Haah [33], and may be thought of as (Type-II [37]) symmetry-enriched fracton topologically ordered (FSET) phases.

One 2D cellular automaton
A 2D CA has a two-dimensional state space, combined with one time direction. The state of such a CA may be straightforwardly represented by a polynomial in two variables, s t (x, z), where the state of the (i, k)th cell is given by the coefficient of x i z k . The update rule is given as a two variable polynomial f (x, z), such that s t+1 = f s t as before. Two dimensional CA also result in a rich variety of fractal structures [101]. The classical Hamiltonian takes the form with symmetries on the semi-infinite system (with y j≥0 ) given by which commutes with H 1CA everywhere. On an infinite system, an inverse evolution f −1 may be defined analogous to Eq 21 and the symmetry action takes the form S(q(x, z)) = X(q(x, z)F (x, y, z)) with The discussion of Sec 4 and 5 may then be generalized in a straightforward manner. The phase diagram is exactly the same as in 2D, given by Fig 7(left). As an example model, consider the Sierpinski Tetrahedron model, given by the update rule f (x, z) = 1 + x + z. The Hamiltonian is given by The fractal structure of the symmetries for this model are Sierpinski Tetrahedra, with Hausdorff dimension d = 2. The quantum model may be constructed which exhibit the same properties: selfduality about h = 1, spontaneous fractal symmetry breaking, and instability to non-zero temperatures. A cluster FSPT version may also be constructed, with the Hamiltonian This cluster FSPT also has the nice interpretation of being the cluster model (Eq 29) on the diamond lattice. In the presence of an edge, terms in the Hamiltonian must be excluded leading to degeneracies, and in exactly the same way as in 2D one finds these degeneracies along a surface cannot be gapped, thus leading to a 2 O(L 2 ) overall symmetry protected degeneracy for an open system. with s t1+1,t2 = f 1 (x)s t1,t2 and s t1,t2+1 = f 2 (x)s t1,t2 . Interpreting the y, z, directions as the t 1 , t 2 , directions, the classical 3D Hamiltonian takes the form

Two 1D cellular automata
where α = 1 + f 1 y and β = 1 + f 2 z are defined, and in the second line for notational convenience we have suppressed the x i y j z k factor, when summation over translations is apparent (and we will continue to do so). The fractal symmetries on a semi-infinite system (with x i y j≥0 z k≥0 are of the form) which can be readily verified to commute with everything in the Hamiltonian. On an infinite system some inverse may again be defined and the symmetry takes the form with F 1/2 each defined as in Eq 22 with f 1/2 . The decorated defect construction starting from H 2CA results in the following Hamiltonian, with three spins per unit cell, on which we have operatorsẐ which is illustrated in Fig 8,  : but now the remaining independent symmetry elements are more complicated, which arises because there is a further local operator that commutes with the first three terms in H FSPT , given bŷ Due to the existence ofB ijk , given any symmetry operation S,B ijk S is also a valid symmetry. Thus, these should be thought of as higher form fractal symmetries [106]. Consider the analogy with, say, a 1-form symmetries in 3D: these are symmetries which act along a 2 dimensional manifold which may be deformed by local operations. Here, we have the symmetry operations acting on only b or only c sublattice sites which may be made to live on a single plane, but we are also free to deform such symmetries using products ofB ijk . Such higher form fractal symmetries are an interesting subject by themselves, and we leave a more thorough investigation as a topic for future study. One may confirm that when h x = h z = 0, all these symmetries are products of terms in the Hamiltonian, and therefore must have expectation value 1 in the ground state. As every term is independent, and there are three terms that must be satisfied per unit cell of three sites, the ground state is unique. This model in fact describes an FSPT protected by the combination of the "global" fractal symmetries Z symmetries may be chosen to act as Z on a single site. Thus, the only Hamiltonian we can write down on the surface must be composed of Z (to commute with a local Z) and must commute with the fractal symmetry. The only possibility is therefore the classical Hamiltonian (as in Eq 10), which exhibits spontaneous fractal symmetry breaking in the ground state. Thus, the surface is non-trivial and must either be gapless or spontaneous symmetry breaking. Figure 7(right) shows a sketch the phase diagram for this model. Increasing h x/z drives this model out of the FSPT phase. If we increase only h z while keeping h x small, we arrive at the spontaneously fractal symmetry broken phase like in the 2D FSPT. Increasing both h x and h z too large will result in the trivial paramagnetic phase. However, if we only increase h x while keeping h z small, the system enters into a symmetric fracton topologically ordered phase, which is the subject of the following discussion.

Connection to fracton topological order
The decorated defect approach of the previous sections may be thought of alternatively as the following process: 1. Start with a classical Hamiltonian and some symmetries involving flipping some spins 2. Introduce additional degrees of freedom at each site and couple them to the interaction terms via a cluster-like interaction (this is exactly what one would get following the gauging procedure of Refs [36,37], and adding the gauge constraint as a term in the Hamiltonian).
3. The resulting theory still has the original symmetries, along with some additional symmetry which we may define acting on the new spins, which we take to be the defining symmetries our model. 4. Perturbations respecting these symmetries may then be added to the Hamiltonian (note these may break the gauge constraint from earlier: we are now interpreting both matter and gauge fields as physical).
Most of our models, except the preceding one, were special under this gauging procedure as they allowed for no local gauge fluctuations terms and exhibited a self-duality between the topological and trivial phases. As we will show, in 3D with symmetries defined by two 1D CA, gauge fluctuations are allowed (these are theB ijk operators we found in Eq 67) and there is a phase in which these models exhibit fracton topological order. They may be thought of as the simplest fractal symmetry enriched topological (FSET) phases (this possibility was already hinted at in Ref 36). The phenomenology of the resulting topological orders are the same as those of the Yoshida fractal codes [32]. The Z (a) 2 symmetry will serve the purpose of the enriching symmetry, while the other symmetries will have the interpretation of being logical operators for the underlying Yoshida code.
To avoid complications, let specialize to an L×L×L 3-torus with f L 1 = f L 2 = 1 (x L = y L = z L = 1). The symmetries in this case are given by Eq 66 and 68, but with F 1 = L l=0 (f 1 y) l and F 2 = L l=0 (f 2 z) l instead of F 1 , F 2 , with q still arbitrary. There are L independent Z symmetries. An independent basis for these symmetries are, for α = 0 . . . L − 1, given by and All remaining symmetry elements may be written as products of these andB ijk (as S (b/c) are higherform fractal symmetries). The fracton topologically ordered phase corresponds to the limit in which we take h x in Eq 65 to be large. Expanding about this limit, the Hamiltonian looks like where we have now specified an energy scale G for the second term, the third term is the leading order perturbative correction to the Hamiltonian, and we neglect all the other perturbations. Fixing allX (a) ij = 1 results in exactly the Yoshida code which exhibits a ground state degeneracy (with our geometry and choice of f 1/2 ) of 2 k with k = 2L. From the perspective of the original FSPT, one finds that the charge of all the S (a/b) α (Eq 70) in the ground state of this phase no longer has to be +1, but instead may be ±1. These are exactly the logical operators of the Yoshida fractal code [32]. This transition may also be thought of as some kind of non-local spontaneous symmetry breaking of the higher form fractal symmetries Z The ground state must still be uncharged under the Z (a) 2 . We define the fracton excitation as an excitation of only the first term in the H FSET (these are the relevant charge excitations when G is large). Such an excitation may be created in multiplets by (for example) an operator of the form which creates only excitations of the first term at locations given by the non-zero coefficients in 1 + (f 1ȳ ) ℓ , and is a few-body creation operator whenever ℓ = 2 l . A single such excitation clearly carries charge −1 under some Z symmetries. This Hamiltonian therefore describes a fracton topologically ordered phase, enriched by an additional Z (a) 2 symmetry, and is a genuine FSET. As a helpful analogy, in Appendix A we show how, in exactly the same way, a relaxed Ising gauge theory may be interpreted as an SPT protected by a global Z 2 and 1-form Z 2 symmetry, and in a certain limit describe an SET phase enriched by a global Z 2 .
A single charge is immobile, as discussed in Sec 5.4, provided that f 1 and f 2 are not algebraically related, the same condition which implies the lack of a string-like logical operator in the Yoshida code [32]. Finally, we note that Haah's cubic code [33] is isomorphic to this type of model, but with a second-order CA along one time direction [32].

Conclusion
We have constructed and characterized a family of Ising Hamiltonians that are symmetric under symmetry operations which involve acting on a fractal subset of spins. Fractal structures on a lattice are taken to be those defined by cellular automata with linear update rules. We discuss some possible phases in systems with such symmetries.
These include the trivial symmetric and spontaneously symmetry broken phases which are symmetric under the fractal symmetry. These fractal symmetries form the total symmetry group (Z 2 ) k , where k will depend strongly on system size and topology. We then construct non-trivial symmetric phases, FSPT phases, via a decorated defect approach. For fractal symmetry groups generated by a single CA, the decorated defect construction leads to a family of cluster type Hamiltonians which have a non-trivial gapped ground state under the symmetry group (Z 2 × Z 2 ) k of fractal symmetries. We characterize such a phase by means of symmetry-twisting, ungappable edge modes, and immobile excitations protected by the set of all fractal symmetries.
In three dimensions, our construction leads to an FSPT protected by a combination of the usual fractal symmetry along with a higher form fractal symmetry. Aside from the FSPT phase one also has the possibility of fracton topological order, enriched by the fractal Z 2 symmetries. The topological order in these models are those of the Yoshida fractal codes [32]. While maintaining our fractal symmetries, these topologically ordered phases may be thought of as simple fractal symmetry enriched topological phases (FSET), in which an elementary excitation is charged under the fractal symmetries.
This construction is may also be generalized to higher D-dimensional systems, where one may consider fractal symmetries generated by n d-dimensional CA, with D = n + d. The D = 3 dimensional models examined in this paper have (n, d) = (1, 2) and (2,1). This suggests an avenue towards constructing higher dimensional fractal SPT or topological phases. Finally, we note that a generalization to p-state Potts variables, rather than Ising, is also possible.
Recent work in Ref [107] develops a gauging/ungauging procedure for quantum error-correcting codes. They provide a prescription for obtaining a D dimensional SPT from a gapped domain wall of a D + 1 dimensional quantum code. The example provided of a 2D FSPT obtained from this method is exactly isomorphic to our Sierpinski FSPT.