Fractalizing quantum codes

We introduce"fractalization", a procedure by which spin models are extended to higher-dimensional"fractal"spin models. This allows us to interpret type-II fracton phases, fractal symmetry-protected topological phases, and more, in terms of well understood lower-dimensional spin models. Fractalization is also useful for deriving new spin models and quantum codes from known ones. We construct higher dimensional generalizations of fracton models that host extended fractal excitations. Finally, by applying fractalization to a 2D subsystem code, we produce a family of locally generated 3D subsystem codes that are conjectured to saturate a quantum information storage tradeoff bound.


I. INTRODUCTION
Recent years have seen increasing interest in various spin models that are defined on regular lattices and nevertheless exhibit "fractal" properties. These include gapped spin liquid models in which immobile topological excitations are created at the corners of operators with fractal support, or spin models with symmetries that act on a fractal subsystem.
One example of the former is Haah's code 1 , a canonical model of type-II 2 fracton topological order 2-12 . Such 3D phases are characterized by topological quasiparticle excitations that are strictly immobile. As quantum codes, they lack string-like logical operators, and instead have logical operators supported on a fractal subset of sites. The fractal nature of these codes leads to promise as quantum memories [13][14][15] . More generally, fracton phases have received tremendous attention in a wide variety of contexts  .
An example of the latter type of fractal models are the fractal Ising models 50 . These are classical spin models on the square lattice with symmetries that flip the Ising spins on a fractal subset of sites. These have been studied as classical codes [51][52][53] for their information storage capacity and also as translation invariant models of glassiness 54 . In these codes, classical information is stored in the spontaneous symmetry-broken ground states of the model. With the same set of fractal symmetries, there may also be non-trivial phases without spontaneous symmetry breaking, known as fractal symmetry-protected topological (SPT) phases 50,55,56 . An example is the cluster model 57 on the honeycomb lattice, which realizes a non-trivial SPT phase protected by fractal subsystem symmetries 50 . In 2D, such phases have received attention as they have been proven to be useful resources for universal measurement-based quantum computation [58][59][60][61] .
In this paper, we unify all these fractal spin models via a process called "fractalization". Fractalization maps certain operators defined on a D dimensional lattice to one on a D + m dimensional lattice, taking as input a set of m-dimensional linear cellular automaton (LCA) rules. In particular, the various fractal models above may all be understood as fractalized versions of simple lower-dimensional models. A class of type-II fracton topological phases 9 in 3D, including Haah's cubic code, may be understood as fractalized 2D toric codes 62 ; the fractal Ising models in 2D are simply fractalized 1D Ising models; and the 2D cluster models realizing fractal SPT phases are fractalized 1D cluster models. We also present an interpretation of the fractalization map as a three-step process which directly relates many properties of the fractalized model to that of the original. For instance, we show that the ground state manifold of a fractalized commuting Hamiltonian is (non-local) unitarily related to that of stacks of the original model.
We discuss the nature of the excitations in fractalized models in relation to those of the original model, taking as examples the toric code in various dimensions. Pointlike excitations of the original model lead to immobile point-like fracton excitations in the fractalized model, whose immobility is guaranteed by the lack of string-like logical operators. When the original model has looplike excitations, on the other hand, we determine a criterion under which the fractalized model lacks any loop or string-like excitations whatsoever. The excitations of such models are instead "fractalized loop" excitations: extended deformable fractal excitations with energy cost scaling faster than linearly with the linear size.
When applied to quantum codes, fractalization results in a higher dimensional code with information storage capabilities similar to that of decoupled stacks of the original, but with improved code distance. In particular, for a locally generated [n, k, d] code in D dimensions, where n is the number of physical qubits, k is the number of logical qubits, and d is the code distance, the D → D + m fractalized code is a locally generated [L m n, L m k, d ′ ] code in D+m dimensions, where d ′ ≥ d depending on the LCA rules and L is the linear size of the system. Fractalization therefore has a useful application in deriving new fractal codes from known lower-dimensional codes. By fractalizing the 4D toric code, for instance, we construct a 5D model which lacks any membrane logical operators, a direct generalization of Type-II fracton phases in 3D whose hallmark is the absence of string logical opera-tors 2 . We also discuss in detail fractalization applied to the 2D Bacon-Shor subsystem code 63,64 , which results in the 3D fractal Bacon-Shor code. By incorporating ideas from Bravyi 65 and Yoshida 52 , we are able to improve the information storage capabilities of the 3D fractal Bacon-Shor code to [n, k, d] ∼ [L 3 , L 2 , L η ], where we conjecture η → 2 in a limit of large physical qudit dimension. To the best of our knowledge, if our conjecture proves true this would be the first code to saturate the information storage tradeoff bound for locally generated subsystem codes in 3D 65,66 , k √ d ≤ O(n). In Sec II, we review how LCA rules generically lead to fractal structures along with other necessary background. In Sec III, we define and discuss the fractalization procedure and its interpretation in terms of a three-step process. In Sec IV, we go through each of the mentioned examples in detail, and discuss connections to other recent works. In Sec V, we discuss fractalized loop excitations in detail through examples, including higher dimensional toric codes. In Sec VI, we present the fractalized Bacon-Shor code and discuss its memory storage capabilities. Finally, concluding remarks and future directions are laid out in Sec VII.

A. Fractals and linear cellular automata
We begin by reviewing linear cellular automata (LCA) and their connection to fractals. Let c r,t ∈ {0, 1} represent the state of cell r = (r 1 , . . . , r D ) ∈ Z D of a cellular automaton at time t. The LCA update rule determines how the state is determined as a linear function (modulo 2) of the state at previous time. Throughout this paper, we are mostly concerned with first-order, local, and translation invariant LCA update rules. That is, where the binary matrix F specifies the LCA rules, and all arithmetic is implied modulo 2. First-order means that c t only depends on the states c t−1 at time t − 1, locality implies that F is only non-zero for small |r − r ′ |, and translation invariance means that F r,r ′ = f r−r ′ only depends on the difference r − r ′ . Starting at time t = 0 with the state c 0 , the space-time trajectory at all future times is completely determined by The set of cells with c r,t = 1 will generically form a fractal subset of the (D + 1) dimensional space-time lattice. This is most naturally seen by adopting a polynomial representation for LCAs. In the polynomial representation, the state c t is represented by a D-variate polynomial with F 2 coefficients, , and x r ≡ i x ri i is a monomial. Similarly, the update rule is represented by in terms of which Eqs. 1 and 2 take the particularly simple form We will utilize both the vectoral representation and the polynomial representation, depending on whichever one is most convenient.
We always assume f (x) is chosen to have a non-zero constant term and only positive powers of x i , and that f (x) = 1 is non-trivial. To see that this generates fractal structures, notice that for any t = 2 n , we have, due to the properties of Assuming the initial state c 0 (x) only contained finitely high powers of any x i , then at exponentially long times c t (x) = f (x 2 n )c 0 (x) describes copies of the initial state each shifted by 2 n for each term in f (x).
We most often consider LCA defined on finite (L 1 , . . . , L d ) systems with periodic boundary conditions, which is enforced by the identification x Li i = 1. An important property is the reversibility of an LCA rule for a given system size, which corresponds to the existence of f (x) −1 satisfying f (x) −1 f (x) = 1, after taking into account periodic boundary conditions. A simple family of reversible LCA is given by any polynomial f (x) satisfying f (1) = 1 on a system of size L i = 2 n . This follows from the fact that f (x) 2 n = f (x 2 n ) = f (1) = 1, and so the inverse is given by f (x) −1 ≡ f (x) 2 n −1 .
As an example, consider d = 1 spatial dimensions, f (x) = 1 + x (corresponding to F r,r ′ = δ r,r ′ + δ r,r ′ +1 ) and the initial state c 0 (x) = 1 (corresponding to c r,0 = δ r,0 ). Then c t (x) = f (x) t , and listing out terms for the powers of f (x) we find and so on. Zooming out, the fractal generated is the Sierpinski triangle (Pascal's triangle mod 2). We remark that f (x) = 1 + x is irreversible on any system size. A simple example of a reversible LCA is f (x) = 1 + x + x 2 on sizes L = 2 n .

B. Pauli operators
In this paper, we present a procedure by which a set of LCA rules may be used to extend a quantum CSS code into a higher dimensional fractal code. To this end, it is useful to use a vectoral (and polynomial) representation of Pauli operators 9,67,68 . Let us consider qubits on the sites r of a hypercubic lattice. Acting on these qubits, we have Pauli operators X r and Z r . Define a map σ X which maps a binary vector a r to a tensor product of X Pauli operator as and similarly σ Z , The commutation relation between two Pauli operators is straightforward to compute, Using the correspondence between vectors and polynomials, a ↔ r a r x r , σ X/Z may also be interpreted as a map from F 2 polynomials to Pauli operators. Translations are nicely expressed in this language as multiplication by a monomial: σ X [a(x)] translated by r is σ X [x r a(x)]. Note that we use the same symbol (σ X/Z ) in both the vectoral and polynomial representation.
The commutation relation for translations of two operators can be neatly summed up by a single commutation polynomial c(x), where andx ≡ x −1 . In particular, if c(x) = 0, then all translations commute. Throughout this paper, we typically consider more than one physical qubit per lattice site. The generalization to N qubits per site, with operators X (n) r and Z (n) r for n = 1 . . . N , is straightforward: in the vectoral representation, σ X/Z now takes as input N vectors (or a tensor), a r,n , one for each of the qubits. Similarly, in the polynomial representation, σ X/Z takes as input a vector of N polynomials, a(x) = a n (x). The commutation polynomial is then given by c(x) = a T (x)b(x).

A. Definition
We are now ready to discuss fractalization. In the most general case, fractalization maps a D-dimensional X or Z Pauli operator to one in D + m dimensions. We consider a hypercubic lattice with N qubits per site. As input, we take a set of D, m-dimensional LCA rules, Let A be a local X Pauli operator in D dimensions, in the polynomial representation, where x = {x 1 , . . . , x d } and we have chosen an "anchor point" r 0 . We uniquely specify the anchor point for this operator by the requirement that a(x) has a non-zero constant term and only positive powers of x i . Locality means that a(x) contains only finitely high powers of x i . Going to D + m dimension, let s denote the position vector along the final m dimensions, such that the full position of a site is (r, s). Similarly, in the polynomial notation, let y denote the variables corresponding to the final m dimensions. For each s 0 , the fractalized operator A frac s0 is then defined as where Thus, fractalization maps a D-dimensional operator to a set of D +m-dimensional operators, A → {A frac s0 }, one for each shift s 0 in the new dimensions. Most of the examples in this paper will involve only m = 1 new dimension.
Similarly, suppose we start out with the Z Pauli operator, where now we have chosen the anchor r 0 such that b(x) has a non-zero constant term and only positive powers of x ≡ x −1 . Then, the fractalized operator is defined as For completeness, let us also express fractalization using the vectoral representation. The input LCAs are given by matrices we find the anchor point r 0 such that a r0,n is non-zero (for some n), and a r,n = 0 if (r − r 0 ) i < 0 for any i. Then, the D + m dimensional fractalized operator is where the ordering of each F i in the product is arbitrary, as all F i commute due to translation invariance (this is obvious from the polynomial representation, ). Note that since a r,n = 0 if (r−r 0 ) i < 0, F i is always raised to a non-negative power. Similarly, for a Z operator, we must choose the anchor r 0 such that b r0,n is non-zero for some n and b r,n = 0 if (r 0 − r) i < 0. Then, Given a CSS stabilizer group S = {O l } , or Hamiltonian where O l are local X or Z generators of the stabilizer group, we may define the fractalized stabilizer group S frac = {O frac s,l } , or Hamiltonian which, as we will show, inherits many of the properties of H. Note that H frac is always translation invariant along the final m dimensions by construction. Finally, fractalization may also be applied to non-local operators, provided certain conditions are satisfied. If the original system has open boundary conditions, fractalization can be applied straightforwardly. In the case of periodic boundary conditions, however, it may not be possible to choose an anchor point consistently for a non-local operator. Suppose an operator is non-local in the ith direction. If f i (y) Li = 1, then the fractal is commensurate with the system size and the anchor r 0,i can be chosen arbitrarily. If this is not the case, then we must find a set of polynomials Q i = {q(y) : q(y)f i (y) Li = q(y)}, and make a choice of basis polynomials {q j (y)} j for Q i . Then, an operator A or B maps on to a set of fractalized operators, one for each q j (y). Eq 13 and 16 generalize to where the choice of r 0,i can be arbitrary. Thus, each operator only maps on to log 2 |Q i | fractalized operators.
If an operator is non-local in multiple directions i ∈ I, then q j (y) must be chosen from ∩ i∈I Q i such that the fractal is commensurate with all directions.
B. Properties

Locality
We first point out that fractalization is locality preserving. As long as A has support only on a local patch of the d-dimensional lattice, A frac s0 will also have support on a local patch of the D + m-dimensional lattice near s 0 . This follows from the locality of the LCA rules: is only non-zero for small |s − s 0 |.
Non-local operators map on to non-local, potentially fractal, operators. To see what this means, consider the D = 1 chain with length L and open boundary conditions, and the 1-dimensional LCA rule f (y) = 1 + y. The non-local operator maps on to the fractal operators in two dimensions. Recall that f (y) n generates the nth row of the Sierpinski triangle fractal (Eq 6). Hence, S frac s is an operator that acts as X on sites along a Sierpinski triangular fractal subsystem, and s just denotes an overall shift in the y direction. This operator has support on a fraction of sites scaling as L d f , where d f ≈ 1.58 is the Hausdorff dimension of the Sierpinski triangle.

Commutativity
The next property of fractalization is commutativity preservation.
Suppose we have two opera- which has a zero constant term iff [A, B] = 0. The commutation polynomial for A frac s0 and B frac s1 is Therefore, if c(x) has a zero constant term, c frac (x, y) also has a zero constant term regardless of s 0 , s 1 . It is also instructive to work in the vectoral representation. Let AB = (−1) c BA, then, commutativity of is zero. The fractalized operators satisfy which is zero if c = 0. If instead A and B anticommute, then c = 1 and c frac s0,s1 may be zero or non-zero depending on their separation s 1 − s 0 .

C. Three-step process
We now present a three-step process which is equivalent to fractalization, but offers insight into the relation between the original and fractalized operators. The three steps are 1) constructing a layered system, 2) unitary transformation, and 3) choosing a different set of generators. Each of these three steps is explained in detail below. While fractalization may be applied to systems of arbitrary size and LCA rules, the three-step process is only applicable under certain conditions: the boundary conditions of the D original directions must be either open or periodic with F Li i = 1, F i must all be invertible, and the boundary conditions along the m new directions must be periodic.
Let us start with an X or Z operator in the vectoral representation, A = σ X [a r,n ] or B = σ Z [b r,n ], defined in D dimensions, and choose an anchor point r 0 in the same way as before. The first step is to extend these operators to a layered system, with each layer labeled by its coordinate s along the perpendicular direction. The operator A or B now maps on to a set of operators on each layer separately, on the D + m dimensional lattice. The next step involves a non-local (along the s directions) unitary transformation. The unitary transformation involves a series of controlled-X (CX) gates. In particular, for any invertible binary matrix M with M ii = 1, there is a corresponding CX circuit which implements the transformation In Appendix A, we give a general algorithm for determining a CX circuit for any such M. Here, we choose the matrix which are highly non-local operators along the s directions as they involve F i being raised to potentially high powers r i .
The third and final step involves choosing a different set of generators for the groups {Ã s0 } and {B s0 } . For X operators we choose which matches Eq 18. Similarly, for Z operators, results in the fractalized operator. Since F i are all invertible, this just corresponds to a different choice of generators for their generated groups. AlthoughÃ s0 were highly non-local, A frac s0 are all local operators. Thus, locality is restored.
This three-step process has important implications for fractalized quantum codes. Consider the CSS stabilizer Hamiltonian with the 2 k -dimensional ground state manifold obtained as the simultaneous +1 eigenstate of all O l . Let us choose LCA rules and the system sizes such that the three-step process can be applied. After the first step, H → s H s simply describes a layered system with ground state degeneracy 2 Lk . The second step is a unitary transformation, H → U HU † , which does not affect the ground state degeneracy. The ground states themselves are related to the original ground states of the layered system by the non-local unitary transformation |ψ → U |ψ . Finally, since all ground states are +1 eigenstates of each term, choosing a different set of generators for the stabilizer group does not affect the ground state manifold (although higher excited eigenstates are shifted). Thus, starting with a locally generated CSS stabilizer Hamiltonian with ground state degeneracy 2 k , the fractalized Hamiltonian will have a 2 Lk degenerate ground state manifold which is unitarily related to those of the layered system by the non-local unitary transformation U . In particular, this implies that upon compactifying the additional m dimensions the fractalized model is in the same D dimensional gapped quantum phase of matter as the layered system.

Beyond qubits
Everything discussed above straightforwardly generalizes from qubit degrees of freedom to Z p qudits of prime dimension p. The X and Z Pauli operators are replaced by Z p clock operators, satisfying the relation XZ = e 2πi p ZX. Fractalization now takes as input p-state LCA rules, specified by polynomials f (x) with F p coefficients. Such LCA still generate fractal structures due to the property f (x) p n = f (x p n ), where self-similarity occurs at scales p n . Our previous discussions readily applies, with arithmetic now done modulo p.

Higher order LCA
The fractalization procedure generalizes to higher order LCA. In an nth order LCA, the state at time t, c t , may be determined by the states at times back to t − n, where the LCA update rule is now specified by a list of n polynomials {f (j) (y)}. Fractalization from D to D + m dimensions now takes as input D, m-dimensional, n-th order LCA rules, {f (j) i (y)}. One possible generalization for Eq. 18 is and This generalization is chosen such that the fractalized Ising model (see Sec IV A) has as its ground state the set of valid space-time LCA trajectories. However, this higher order generalization no longer has as many of the nice properties as before. For one, [A, B] = 0 is no longer a sufficient criterion for the commutativity of the fractalized operators. Instead, commutativity is only guaranteed if all translations commute. Let A r and B r be translations of A and B, such that their anchor point is at r. Then, [A r0 , B r1 ] = 0 for all r 0/1 , implies that [A frac r0,s0 , B frac r1,s1 ] = 0 for all r 0/1 , s 0/1 . This follows from the fact that, for operators whose translations all commute, Thus, fractalization with higher order LCA should only be applied for translation invariant commuting codes.
We also do not have a generalization of the threestep process for fractalizing higher order LCA, even for translation invariant codes. Thus, the precise connection between the original and fractalized code for a general higher order LCA remains an open question.

IV. EXAMPLES AND DISCUSSION
In this section, we go through a number of examples of fractalization applied to well-understood models, connecting them with more exotic fractal models in the literature. We start with the 1D quantum Ising model, which we show fractalizes to a spin model studied by Newman and Moore 54 as a translation invariant model of glassiness, and also as a classical code [51][52][53] . We then move on to the cluster model 57 and find that it fractalizes into cluster states on higher dimensional lattices, known examples of fractal subsystem SPT phases 50 which have been shown to be useful for universal measurement-based quantum computation 58,59 . We then fractalize Kitaev's 2D toric code model 62 , which results in Yoshida's fractal spin liquids 9 , a family of models of type-II fracton topological order.
A. 1D Ising model The 1D Ising model is defined on a chain with one qubit per site r. The Hamiltonian simply consists of pairwise nearest neighbor Ising interactions (42) which, after fractalization with a 1 dimensional LCA f (y), becomes Using the choice f (y) = 1 + y, we arrive at the Newman-Moore model 54 , where (r, s) labels the x, y coordinates on the square lattice. The global Z 2 symmetry of the Ising model maps on to a set of fractal symmetries, which act on a subsystem of sites corresponding to the Sierpinski triangle fractal, shifted by s along the y direction. The ground states of H frac 1DIsing are spin configurations that are valid space-time LCA trajectories. Let us denote a product state in the Z (r,s) basis by |{z (r,s) } , where Z r |{z (r,s) } = (−1) z (r,s) . Then, z (r,s) forms a valid LCA trajectory, with r playing the role of time, and s playing the role of space. This remains true in higher dimensions, as well as for higher-order LCA. For the D-dimensional Ising model, fractalized with m-dimensional LCA, the resulting D + m-dimensional H frac 1DIsing has as its ground state valid space-time LCA trajectories, where r are time dimensions, and s are space dimensions. Multiple time dimensions, each governed by their own update rule, are possible since translation-invariant LCA rules commute.
The Ising chain may also be viewed as the classical repetition code, whereby the logical bit is represented by the two symmetry-broken ground states. Similarly, the fractalized Ising chain, in 1 + m dimensions, is also useful as a classical code 51-53 , in which logical bits are again represented by the symmetry-broken ground states, of which there are now 2 O(L m ) . The fractal symmetry operators S frac s are logical operators that flip between the symmetry-broken ground states. Upon generalizing to Z p degrees of freedom, this code saturates an information storage tradeoff bound 69 in the limit p → ∞ 52 .
The nature of the excitations of the fractalized Ising model can be deduced from the excitations of the Ising model. The Ising model in 1D has point-like domain-wall excitations, which are created from a ground state |ψ by flipping a segment of spins, making it clear that W only creates two excitations at ℓ 1 and ℓ 2 when acting on the ground state. From Eq. 29, the fractalized versions of these operators satisfy Therefore, B frac s1 creates a single excitation at (ℓ 1 , s 1 ), while creating potentially many excitations along r = ℓ 2 , depending on F ℓ2−ℓ1 . However, if ℓ 2 − ℓ 1 = 2 n , f (y) 2 n = f (y 2 n ) and so (F 2 n ) s1,s is only non-zero for a few s; for f (y) = 1 + y, only two further excitations are created at (ℓ 1 + 2 n , s 1 ) and (ℓ 1 + 2 n , s 1 + 2 n ). The excitations of H frac 1DIsing are therefore point-like, and may be created in special configurations separated by distances 2 n via fractal operators. They are an example of (nontopological) fractons.
This can be generalized to higher dimensions. The 2D Ising model, for example, has loop-like domain wall excitations and is treated in section V.

B. 1D Cluster model
The 1D cluster model, in CSS form, has two qubits per site and is described by the Hamiltonian After fractalization with a 1-dimensional LCA f (y), which describes another cluster state on a 2D bipartite lattice.
With the choice f (y) = 1 + y, for example, H frac clus describes the cluster state on the honeycomb lattice 50 . Higher dimensional models may be created by using m-dimensional LCA rules, for example, using the 2-dimensional LCA rule f (y) = 1 + y 1 + y 2 , one arrives at a 3-dimensional cluster state, with Sierpinski tetrahedron fractal subsystem symmetries.
The cluster model describes a non-trivial symmetryprotected topological (SPT) phase, protected by the Z 2 × Z 2 symmetry that is generated by S 1 = r Z (1) r and S 2 = r X (2) r . The fractalized cluster model also describes an SPT phase, protected by the set of fractal subsystem symmetries generated by S frac α,s , α ∈ {1, 2}. The non-triviality of H frac clus as a fractal subsystem SPT follows from that of H clus . The cluster chain H clus realizes a non-trivial SPT phase, as can be observed from the fact that, at a boundary, the G = Z 2 × Z 2 symmetry acts as a non-trivial projective representation 70 . We may apply fractalization directly to the open cluster chain and its symmetry operators, following the threestep process. In the first step, the system is composed of decoupled cluster chains along x, each with their own symmetry S α,s realizing a projective representation on the left edge. The second step is a unitary transformation, local along x, which does not affect the projective representation. The third step is a different choice of generators for the symmetry group {S α,s } . We can simply choose the origin such that the anchor point for S α,s is r 0 = 0, in which case the change of generators (Eq 34) is trivial. Thus, S frac α,s realizes the same projective representation in H frac clus as S α,s from the layered H clus . If we take periodic boundary conditions along the x direction as well, along with an irreversible LCA, then the three-step process cannot be applied. Following the discussion in Sec. III A, the symmetry group of H frac clus is generated by {S frac α,j }, which is smaller than the symmetry group of the layered system (or may even be trivial in the most extreme case), depending on the number of solutions to q(y)f (y) = q(y). These are pseudo-SPT phases, as defined in Ref. 56.
For a given set of fractal symmetries, there are an infinite number of non-trivial SPT phases 56 . The fractal subsystem SPT phase obtained from fractalizing the cluster state is only the simplest non-trivial phase. The remaining phases are generically not strictly translation invariant (they are translation invariant with a larger period 2 n ), and therefore cannot be realized by fractalization.

C. 2D toric code model
The 2D toric code 62 is defined on the square lattice with a qubit on every bond. Associating with each site the two qubits on bonds straddling it in the +x 1 and +x 2 direction, the Hamiltonian may be written as Applying fractalization with two 1-dimensional LCA rules, f 1 (y), f 2 (y), we get which is Yoshida's fractal spin liquid model 9 . We take H TC on an L 1 × L 2 torus, and H frac TC on an L 1 × L 2 × L 3 3-torus.
The toric code describes a 2D topologically ordered phase. The fractalized toric codes describe a family of (type-II) fracton topological ordered phases. Fracton topological order is characterized by a subextensive ground state degeneracy (scaling with an envelope 2 O(L) ) and quasiparticle excitations with restricted mobility. When the three-step process for fractalization holds, the ground state degeneracy of H frac TC is 2 2L3 , but will be less for generic system sizes where f i (y) Li = 1.
The excitations of the toric code are created at the ends of string operators. Fractalizing these string operators and following a similar analysis as in Sec. IV A, one finds that the excitations of the fractalized toric code are also point-like and must be created in particular configurations. As pointed out in Ref. 9, these quasiparticle excitations are only strictly immobile if f 1 (y) and f 2 (y) are algebraically unrelated. Two polynomials are algebraically related if there is a non-trivial solution to for any finite constants n i ,c, without periodic boundary conditions. Haah's cubic code 1 can also be viewed as a higher order fractal spin liquid 9 . Using the generalization of fractalization to higher order LCA, with f (1) 1 (y) = 1 + y + y 2 f one arrives at a model that is locality preserving unitarily related 9 to Haah's cubic code, up to a redefinition of the lattice vectors.
The toric code may also be regarded as a limit of the Ising gauge theory (IGT), which is obtained by gauging the global Z 2 symmetry of the 2D quantum Ising model H 2DIsing . After gauge fixing, the IGT takes the form of a perturbed toric code, where n = 1, 2 labels the qubits on each site. H IGT is in the topologically ordered phase when J/h, Γ/K ≪ 1. The fractalized IGT, H frac IGT , can either be obtained by fractalizing H IGT or alternatively by applying a generalized gauging procedure 2,11,71 to the fractal subsystem symmetries of the fractalized 2D Ising model. Thus, the generalized gauging procedure naturally arises in this context as the "fractalized" version of the usual gauging procedure.
In Refs. 72 and 73, the effect of dimensional compactification of 3D fracton models was investigated. Compactifying a 3D translation invariant type-II fracton code along one dimension leads to a 2D code which is equivalent to copies of toric code and trivial states by a locality-preserving unitary transformation [74][75][76] . Take the L 1 × L 2 × L 3 fractalized toric code models, and consider compactifying along the third dimension. In cases where the three-step process for fractalization can be applied, we can easily show that the resulting compactified model is local unitarily equivalent to L 3 copies of the toric code: The first step creates the layered toric codes, the second step is a unitary transformation, and the third step is an irrelevant change of generators for the stabilizer group. The unitary transformation, the most important step, is non-local only in the third direction. Thus, once the third dimension is compactified, the second step is a local unitary transformation which transforms the ground state into L 3 copies of the toric code ground state.
In a recent manuscript 77 by the present authors, a coupled cluster state construction was presented for a number of models, including the toric code and many type-II fracton phases. The Hamiltonian takes the form H = H clus + λH coup , where H clus is the cluster Hamiltonian on decoupled chains or layers, and H coup couples them together. By choosing the the cluster states and couplings appropriately, as discussed in Ref. 77, the low energy effective Hamiltonian in the limit λ → ∞ coincides with the desired fracton or toric code model. In particular, a construction of the toric code in terms of coupled wires, and Yoshida's type-II fracton phases 9 in terms of coupled layers was obtained. This coupled layer construction of Yoshida's models follows directly from fractalization applied to the coupled wire construction of the toric code. Specifically, the Hamiltonian H describing the coupled wire construction of the toric code can be written in terms of only X or Z operators. Applying fractalization, H → H frac , the cluster state wires map on to fractal cluster state layers, and in the limit λ → ∞, H frac realizes the fractalized toric code, which is exactly Yoshida's models.

V. FRACTALIZED LOOP EXCITATIONS
In this section we describe several examples that involve fractalizing stabilizer models with loop-like excitations. We find that the excitations in the fractalized models are no longer generically loop or string-like, but form more general extended "fractalized loop" excitations. We determine the criterion for such models to lack any stringlike excitations whatsoever. This implies the lack of any membrane logical operators and can lead to a boost in the energy barrier for applying a logical operator from O(L) in the original model to O(L α ), for some 1 < α ≤ 2 determined by the choice of fractalization, where L is the linear size of the system.
For this section, we make use of a more compact notation for the sum of K Hamiltonian terms, where is an N × K matrix whose columns correspond to each of the Hamiltonian terms (and rows correspond to the physical qubits on each site as before). A similar notation is also used for σ Z [B(x)]. This representation is useful for translation invariant models with multiple types of X or Z terms, such as the higher dimensional toric code models.
A. 2D quantum Ising model As a warmup we fractalize the 2D quantum Ising model which supports looplike domain wall excitations. The Hamiltonian is given by with where, as explained above, the two columns correspond to the Ising interactions along the x 1 and x 2 directions. The logical operators for the code on periodic boundary conditions are given by which represent a global X spin flip and single site Z operator, respectively. Hence, the Ising model is not topologically ordered as the ground space degeneracy can be split locally. A rectangular loop excitation can be created with a truncated membrane operator, where R = R 1 × R 2 is a set of sites in some rectangular The excitations created in response to W acting on the ground state is described by the polynomial vector whose elements denote the locations of excited σ Z terms, and the two rows correspond to the two types of σ Z terms. In this case, this means that a line of σ Z [x r (1+x 1 )] terms are excited along the edges of the rectangle at r 1 = 0, l 1 , and a line of σ Z [x r (1 +x 2 )] excitations are created along the other edges at r 2 = 0, l 2 . This is simply the domain wall excitation of the Ising model.
The 3D fractalized 2D Ising model is described by the Hamiltonian for some choice of F 2 polynomials f 1 , f 2 . For periodic boundary conditions L 1 , L 2 , L 3 satisfying f 1 (y) L1 = f 2 (y) L2 = 1 there are L 3 logical qubits encoded within the ground state manifold with pairs of independent logical operators obtained by fractalizing the logical operators in Eq. (61), We see that the fractalized model is also not topologically ordered as the ground space degeneracy is split by local operators.
Suppose we now act on the ground state with a truncated ℓ X s logical operator, i.e. a fractalized W operator, . This creates excitations located at which we refer to as a fractalized loop excitation. It is clear that such excitations have an energy cost scaling faster than linearly with the size l 1 , l 2 . Thus, these are not loop-like excitations. They appear loop-like when projected on to the (x 1 , x 2 ) plane, but extend nontrivially into fractals in the y direction. A natural question is whether there exists any extended topological excitations which are string-like. By topological, we mean that the excitation is not simply a string of locally-created excitation clusters (recall that H frac 2DIsing is not topologically ordered). To answer this question, observe that the Hamiltonian terms obey a local relation on each cube, where a local relation is defined to be a product of Hamiltonian terms that result in unity. The local relation is represented by a polynomial vector l(x) which satisfies B(x)l(x) = 0. Going back to the original 2D Ising model, the local relation where the two columns correspond to the two Hamiltonian terms, tells us that σ Z [x r B(x)l(x)] = 1: a product of four Ising interaction terms around any square plaquette is trivial. This simply means that excitations (which are domain walls) must come in pairs around every square plaquette, and therefore form closed Z 2 loops. For the fractalized Ising model, we have the fractalized local relation l frac (x,ȳ) = l(f(ȳ) •x) which implies that excitations of Hamiltonian terms are only allowed in certain configurations. Meanwhile, a local cluster of excitations can always be created by acting with a single X on a ground state, described by A topological string-like excitation, E string (x, y), would have to be a valid excitation satisfying l frac (x, y) T E string (x, y) = 0, and also not be made up of a product of E X (x, y). This is equivalent to the problem of finding an X-type string logical operator ℓ X = σ X [E string (x, y)] for the "excitation Hamiltonian" which we define in terms of the local relation and local excitation cluster, (69) and is equivalent to the fractalized toric code H frac T C ( Eq. 51). As shown in Ref. 9, this code lacks string-like logical operators iff f 1 and f 2 are algebraically unrelated, as defined in Eq. 54. Hence, for algebraically unrelated f 1,2 , H frac 2DIsing lacks any topological string or loop-like excitations, otherwise there exists string-like excitations (which will only be allowed to lie along a particular direction).

B. 3D toric code
Next we consider fractalizing the 3D toric code, a topologically ordered model with both loop and particle excitations. It is described by the Hamiltonian where The logical operators for the code on periodic boundary conditions correspond to conjugate string-membrane pairs given by where e 1 = (1, 0, 0) T is a unit vector. ℓ X/Z 2/3 are defined similarly by cyclic permutations of 1, 2, 3. Since this model is topologically ordered, no local operators can split the ground space degeneracy. Truncating the string logical operator ℓ X i leads to topological quasiparticle excitations at its ends, and truncating the membrane logical operator ℓ Z i leads to a topological loop excitation along its boundary, when acting on the ground state manifold.
The 4D fractalized 3D toric code model is described by the Hamiltonian for the A and B matrices introduced in Eq. (71) and a choice of F 2 polynomials f(y) = {f 1 (y), f 2 (y), f 3 (y)}. Excitations of the X term are point-like, while excitations of the Z term are fractalized loops. Similar to before, these fractalized loops resemble loops when projected to (x 1 , x 2 , x 3 ), but extend non-trivially in the y direction. Following Ref. 9, a necessary and sufficient condition for the model to have no string logical operators is that the triple f 1 , f 2 , f 3 is not algebraically related. A triple f 1 , f 2 , f 3 is said to be algebraically related iff an equation of the form holds for some finite constants n i , c, where n i are nonnegative integers and not all n i = 0, or for some other permutation of 1, 2, 3, without periodic boundary conditions. We may also ask whether the excitations of the Z terms can be string-like. Repeating the same analysis as before, we may reduce the problem to that of finding string-like logical operators for an excitation Hamiltonian, on which the arguments of Ref. 9 can again be applied. The condition for the existence of string-like excitations is once again Eq. 74. Hence, algebraic independence of f 1 , f 2 , f 3 imply both that there are no string-like logical operators and also no string-like X excitations.
The lack of string-like X excitations further implies the non-existence of a membrane Z-type logical operator (or some fractal structure than can be embedded onto a membrane), since if one existed then its truncation would produce a string excitation around its perimeter. This model still has X-type logical operators which are fractals that can be embedded on to a single (x i , y) plane, however. By fractalizing the 4D toric code in the next section, we construct a 5D code which lacks any membrane logical operators altogether.
For periodic boundary conditions L 1 , L 2 , L 3 , L 4 satisfying f i (y) Li = 1 there are 3L 4 logical qubits encoded within the ground space and a basis of logical operators is given by ℓ X/Z,frac i,s . To construct an explicit example we take f 1 (y) = 1 + y + y 2 , f 2 (y) = 1 + y + y 3 , f 3 (y) = 1 + y 2 + y 3 , which satisfy the algebraic independence condition, as the chosen polynomials are prime, and the boundary conditions when L i = 2 l .

C. 4D toric code
Next we consider fractalizing the 4D toric code which only supports loop-like excitations. The model is described by the Hamiltonian with where the six rows correspond to qubits living on an (x i , x j )-plaquette, with ij = 12, 13, 14, 23, 24, 34, in that order.
On periodic boundary conditions the ground space encodes six logical qubits with logical operators where e 12 is the unit vector corresponding to the qubit on the (x 1 , x 2 ) plaquette. The logical operators ℓ X/Z ij for the remaining ij are defined similarly by permuting 1, 2, 3, 4. ℓ X ij is a membrane operator acting in the plane orthogonal to ij, while ℓ Z ij acts in the plane ij. This model, in fact, extends to a nontrivial finite temperature quantum phase of matter and forms a self-correcting quantum memory at finite temperature, in part due to the linear scaling of the energy barrier that needs to be overcome to implement any logical operator 78 .
The 5D fractalized 4D toric code is described by the following Hamiltonian for the A and B matrices defined in Eqs. loop excitation which appears loop-like when projected to a particular (x i , x j ) plane but extends non-trivially into the y direction. Following the same argument, we find that the requirement for non-existence of string-like excitations is that f 1...4 satisfy a generalized algebraic relation where either or can be satisfied for a set of finite constants n i , c (where not all n i = 0) or for any permutation of 1, 2, 3, 4, without periodic boundary conditions. Hence, for four algebraically unrelated polynomials f i , H 4DTC lacks any string-like excitations. This also implies the non-existence of any membrane logical operators, including those with dimension less than 2 which can be embedded within a membrane. This is therefore a 5D code which lacks membrane logical operators, a direct generalization of Type-II fracton models in 3D which are characterized by the non-existence of string-like logical operators.
To find an explicit example we can take prime polynomials f i . For instance f 1 (y) = 1 + y + y 2 , f 2 (y) = 1 + y + y 3 , To achieve a more local model we can change our choice by the modification f 4 (y) = 1 + y, at the cost of encoding fewer logical qubits for a given system size.

D. Fractalized membrane excitations and beyond
There is another toric code in 4D given by a direct generalization of the 3D toric code to include a further variable x 4 . This model supports point-like and membranelike excitations. Similar to the above examples it can be fractalized, leading to a model that supports "fractalized membrane" excitations. In 6D there is a toric code supporting only membrane-like excitations, which again can be fractalized. For an algebraically unrelated choice of fractalizing polynomials this produces a 7D model with no 3D (or lower) operators.
This follows from a dimensional reduction scheme generalizing our proof for lack of membrane logical operators. The presence of a 3D logical operator in H implies the existence of a 2D (membrane) excitation, which implies the existence of a 2D logical operator in its excitation Hamiltonian H exc (as defined in Sec V A). This in turn implies the existence of 1D (string-like) excitations in H exc , which implies the existence of a 1D logical operator in its excitation Hamiltonian H exc 2 . At this point, the proof of Ref 9 can be applied to show that string logical operators exist for H exc 2 iff f i are algebraically related. The contrapositive of this chain of implications then means that for an algebraically unrelated choice of f i , H does not have any 3D logical operators.
More generally higher dimensional toric codes exist that can support excitations of any dimensionality which can then be fractalized leading to higher dimensional extended excitation generalizations of fracton particles. In particular in (2n + 1)D there is a fractalized toric code that is a generalized type-II fracton model in the sense that it supports no 2n, or lower, dimensional logical operators.
Fractalizing higher dimensional models requires larger sets of algebraically unrelated polynomials to construct type-II models (where the algebraic independence condition is further generalized following the form of the above generalizations). These are easy to find by simply taking all the f i to be prime polynomials of low degree, as the uniqueness of polynomial factorization over a field guarantees algebraic independence.

A. Fractal Bacon-Shor code
In this section we focus on a subsystem code obtained by fractalizing the Bacon-Shor code 63 and Bravyi's generalization thereof 65 . We investigate the code's quantum information storage capacity which we conjecture asymptotically saturates an information storage tradeoff bound 65,66 .
In a subsystem code 63,64 , the full Hilbert space is structured as H = (H L ⊗ H G ) ⊕ H E , where C = H L ⊗ H G is the codespace, composed of the logical and gauge subsystems, and H E is an error subspace. The information is stored in the state of the logical subsystem H L , while the state of the gauge subsystem can be arbitrary and is not protected.
A subsystem code is completely specified by its gauge group G. From G, the stabilizer group S and logical operators follow. In this section, we explain this structure for the Bacon-Shor (BS) code 63 in parallel with the fractalized Bacon-Shor (FBS) code obtained by applying fractalization to the BS gauge group generators.
We take the 2D BS code defined on an L 1 × L 2 lattice with one qubit per site and open boundary conditions. The gauge group G is generated by products of two adjacent Xs in the same row, or Zs in the same column, The Hamiltonian model based on these generators is known as the quantum compass model 79 . The minus signs, although irrelevant for qubits, are included as we later generalize to Z p qudits.
We apply fractalization to the gauge group generators with a pair of 1-dimensional LCA rules, f 1 (y), f 2 (y). The gauge group of the FBS code is on an L 1 × L 2 × L 3 system, with periodic boundary conditions only along L 3 . In the above, r = (r 1 , r 2 ) where r i ∈ R i ≡ {0 . . . L i − 1}, and s ∈ R 3 with s = L 3 and s = 0 identified. The above local gauge generators define a family of Hamiltonians parameterized by the relative coupling strengths of the X and Z type generators, similar to the quantum compass model. Unlike stabilizer codes, such gauge code Hamiltonians may become gapless for certain values of the relative coupling strengths. We remark that when the three-step process of fractalization holds exactly, then the code generated by G frac is unitarily related to the stack of L 3 BS codes. Code properties that are unaffected by unitary transformations, such as the number of encoded qubits, are therefore equivalent to those of L 3 BS codes. Other properties, such as the code distance which is defined below, are affected by the unitary transformation.
The stabilizer group S is obtained as the center of G, S = C(G) ∩ G, where C(·) is the centralizer in the Pauli group. The codespace C is the simultaneous +1 eigenspace of all operators in S. For the BS code, S is generated by the the product of X along two adjacent columns, or Z along two adjacent rows, where S X r1 is defined for r 1 ∈ R 1 \ {L 1 − 1}, and S Z r2 is defined for r 2 ∈ R 2 \ {0} (due to the boundaries). The generators of the stabilizer group of the FBS code, S frac , is then for all s ∈ R 3 . These operators act on two adjacent rows or columns when projected on to (x 1 , x 2 ), like those of the BS code, but spread out non-locally into a fractal along the y direction. The codespace is further decomposed into the logical and gauge subsystems, C = H L ⊗ H G . The gauge operators act on the codespace as 1 L ⊗ G G . Logical operators are those that act non-trivially in H L , and come in two types: bare and dressed logical operators. Bare logical operators act as ℓ L ⊗ 1 G . Dressed logical operators are products of bare logical operators and gauge operators, and therefore act on the codespace as ℓ L ⊗ G G . The set of bare logical operators is defined as L bare = C(S) \ G.
Bare logical operators of the BS code are products of X along an odd number of columns and/or Z along an odd number of rows. Let us focus on the simplest ones, consisting of products along a single row/column, which we choose to be and ℓ Y = iℓ X ℓ Z , from which all other logical operators can be obtained by multiplying with S. Analogously, the corresponding bare logical operators for the FBS code are ]] = (−1) δs 0 ,s 1 , and thus generate an L 3 -dimensional Pauli algebra on H L . This code therefore protects L 3 logical qubits, as expected from the three-step interpretation of fractalization. Dressed logical operators are then obtained as products of bare logical operators and gauge operators.
The code distance d is defined to be the minimum support of a non-trivial dressed logical operator. Let us take all L i ∼ L. The minimum weight of bare logical operators of the FBS code scale as d ∼ L min(d1,d2) where d 1,2 are the Hausdorff dimensions of the fractal generated by f 1,2 52 . However, estimating the minimum weight of dressed logical operators is a difficult optimization problem in general. This is because it is generally possible to reduce the support of a bare logical operator by multiplying it with gauge operators. Let us assume that d ∼ L η , where η is the scaling dimension of the minimum weight dressed operator. If f 1,2 are algebraically related, then it is straightforward to find a dressed logical operator scaling with η = 1. However, if they are algebraically unrelated, we expect that 1 < η ≤ min(d 1 , d 2 ) ≤ 2. The simplest example of two algebraically unrelated polynomials is f 1 (y) = 1 + y and f 2 (y) = 1 + y + y 2 .
We have taken periodic boundary conditions along the third dimension, which may be difficult to implement in systems with locality constraints. This can be easily circumvented. With open boundaries, we choose to keep all the the local gauge generators A frac r,s and B frac r,s with s ∈ R 3 , truncating them at the edge. With this choice of gauge generators, the logical operators are given by Eq. 90 with s ∈ R 3 , also truncated at the edge. As a result, logical operators near the boundary have reduced weight: ℓ Z,frac 0 = r1 Z (r1,0) is the minimal weight logical operator with support on only L 1 sites, and so d = L 1 . To avoid this, we can increase the system size L 3 → L 3 +2δL with where deg is the degree function, and use only the central L 3 logical qubits, {ℓ X/Z,frac s } δL≤s<L3+δL to store information. The logical operators for qubits within this range are unaffected by the boundaries, and the overall scaling of the code parameters with L remains the same.

B. Information storage capacity and the fractal Bravyi-Bacon-Shor code
We now discuss the information storage capacity of the fractal subsystem code. The 2D BS code with n = L 2 physical qubits can encode k = 1 qubits of quantum information with a code distance d = L. The [n, k, d] parameters of subsystem codes with a locally generated gauge group must satisfy the information storage tradeoff bound 65,66 The 2D BBS code is parameterized by a binary matrix K r , which specifies a set of sites K r = 0 to be removed from the 2D lattice. The gauge generators are modified accordingly to skip over the removed sites, i.e. X r X r+nx , where r,r+nx are two nearest unremoved sites along that row. Although this may induce long range couplings, additional auxiliary qubits may be introduced to ensure all generators are local 65 . The end result is that logical operators from different rows or columns can no longer generically be related to one another by S, and the number of logical qubits is given by k = rank(K). Furthermore, it is possible 65 to find a family of matrices K r such that k ∼ L, while maintaining n ∼ L 2 and d ∼ L. Thus, [n, k, d] ∼ [L 2 , L, L] saturates the 2D bound in Eq. 92, kd ∼ L 2 ≤ O(L 2 ). In 3D, however, the current best known subsystem code does not quite saturate the scaling in the tradeoff bound 80 .
The FBS code has n = L 3 , k = L, and d ∼ L η for some 1 ≤ η ≤ 2 as discussed previously. This is still far from saturating the 3D tradeoff bound, k √ d ∼ L 1+η ≤ O(L 3 ). To further improve the quantum information storage capacity, we may apply fractalization to the BBS code, resulting in the fractal BBS (FBBS) code. Removing a single site K r = 0 from the original model corresponds to removal of a whole line of sites {(r, s)} s from the fractal model. The number of logical qubits is improved to k = L rank(K). Again, there is a choice of binary matrices such that the parameters scale as [n, k, d] ∼ [L 3 , L 2 , L η ]. Thus, the FBBS code nearly saturates the tradeoff bound depending on the value of η, k √ d ∼ L 2+η ≤ O(L 3 ). Next, we may further generalize beyond qubits to qudits of prime dimension p. In the limit of large p, the Hausdorff dimensions of the generated fractals approach d 1,2 → 2 (this was used by Yoshida 52 to asymptotically saturate the information storage bound in a classical spin system). We conjecture that in the limit p → ∞, for f 1 , f 2 algebraically unrelated polynomials, the code distance scaling parameter also approaches η → 2. The FBBS code is therefore conjectured to asymptotically saturate the 3D information storage bound k √ d ∼ n in the limit of large p.
We may also generalize to higher dimension. Suppose we wish to apply fractalization to a code from D → D + m dimensions. Let us start with an [n, k, d] subsystem code on an n = L D system, and apply fractalization via the three-step process as described in Sec III C. The first step is to construct a layered system, in which the code now has parameters [L m n, L m k, d]. The second step is a different choice of gauge generators, which does not affect any parameters of the code. Finally, the third step is a non-local unitary transformation. This does not affect n or k, but does affect the distance, which now scales as d ∼ L η , where 1 ≤ η ≤ m + 1. Thus, starting with the 2D BBS code, the 2 + m dimensional FBBS code has parameters [n, k, d] ∼ [L m+2 , L m+1 , L 1≤η≤m+1 ]. The analogue of our conjecture in higher dimensions would be that, for certain choices of LCA rules and Z p qudits, η → m + 1 in the limit p → ∞. In this case, the tradeoff bound (93) is saturated in the large p limit in any dimension m + 2.

VII. CONCLUSIONS
We have introduced the fractalization procedure, by which certain spin models can be extended into higher dimensional fractal models. Many exotic fractal models, such as type-II fracton models or fractal SPTs, may be understood as fractalized versions of familiar spin models, such as the toric code or the cluster model, as we have discussed in detail with many examples. We provided an interpretation of fractalization in terms of a three-step process, which elucidates how the properties of the fractal model are inherited from those of the original model. Fractalizing models with loop-like or higher dimensional excitations, for example, leads to exotic extended fractal excitations.
We have introduced the fractal Bacon-Shor subsystem code which was, after applying Bravyi's generalization 65 and in the limit of large on-site Hilbert space dimen-sion 52 , conjectured to saturate the information storage tradeoff bound in 3D. A more detailed analysis of the code properties of the fractal code is important, but beyond the scope of this current paper. This would include a numerical estimation of the code distance, especially in the limit of large p, as well as finding and implementing error correction algorithms. Due to the presence of nonlocal stabilizers, similar to the BS code, we expect that the fractalized BS code does not give rise to a positive threshold for local noise, without further modifications.
There are still many open directions regarding the fractalization procedure. Currently, fractalization can only be applied to purely X or Z operators while preserving commutativity; it would be interesting if a generalization beyond such operators is possible. An intriguing potential application would then be the fractalization of Kitaev's honeycomb model 81 in the non-Abelian Ising phase to possibly obtain a non-Abelian type-II fracton model. Other possibilities would be generalizing beyond translation-invariant LCA, perhaps to quantum cellular automata. These would allow fractalization to be applied to more complicated spin models. Unfortunately, the present formalism in this paper does not readily admit such generalizations. Another question is whether there exists an interpretation, similar to the three-step interpretation, of fractalization with higher-order LCA. Such an interpretation would greatly improve our understanding of type-II fracton phases described by higherorder LCA, such as Haah's cubic code.
There are also many more CSS codes to which fractalization can be applied, particularly in higher dimensions. In this work we considered fractalizing 3D and 4D toric codes, the latter lead to 5D models with no membrane operators that are, in some sense, a generalization of the type-II fracton models that exist in 3D. Beyond this, each other D dimensional toric code or type-I fracton model defines a family of D + m dimensional fracton models with fractal logical operators and excitations beyond the examples covered in the current work.

ACKNOWLEDGMENTS
TD acknowledges support from the Charlotte Elizabeth Procter Fellowship at Princeton University. DW acknowledges support from the Simons foundation.