#  Fractalizing quantum codes

Trithep Devakul  Department of Physics, Princeton University, Princeton, NJ
08540, USA  Department of Physics, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Dominic J. Williamson  Stanford Institute for
Theoretical Physics, Stanford University, Stanford, CA 94305, USA

(April 15, 2021)

###### Abstract

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.

##  1 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 [3, 4, 5, 6, 7, 8, 9, 10, 2, 11, 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 [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47,
48, 49, 50, 51, 52, 53, 54].

An example of the latter type of fractal models are the fractal Ising models
[55]. 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 [56, 57, 58] for their information storage capacity
and also as translation invariant models of glassiness [59]. 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 [55, 60, 61]. An example is the
cluster model [62] on the honeycomb lattice, which realizes a non-trivial SPT
phase protected by fractal subsystem symmetries [55]. In 2D, such phases have
received attention as they have been proven to be useful resources for
universal measurement-based quantum computation [63, 64, 65, 66].

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 [67]; 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. Point-like 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 loop-like 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\rightarrow D+m$ fractalized code is a locally generated
$[L^{m}n,L^{m}k,d^{\prime}]$ code in $D+m$ dimensions, where $d^{\prime}\geq
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 operators [2]. We also discuss in
detail fractalization applied to the 2D Bacon-Shor subsystem code [68, 69],
which results in the 3D fractal Bacon-Shor code. By incorporating ideas from
Bravyi [70] and Yoshida [57], we are able to improve the information storage
capabilities of the 3D fractal Bacon-Shor code to
$[n,k,d]\sim[L^{3},L^{2},L^{\eta}]$, where we conjecture $\eta\rightarrow 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
[70, 71], $k\sqrt{d}\leq\mathcal{O}(n)$.

In Sec 2, we review how LCA rules generically lead to fractal structures along
with other necessary background. In Sec 3, we define and discuss the
fractalization procedure and its interpretation in terms of a three-step
process. In Sec 4, we go through each of the mentioned examples in detail, and
discuss connections to other recent works. In Sec 5, we discuss fractalized
loop excitations in detail through examples, including higher dimensional
toric codes. In Sec 6, we present the fractalized Bacon-Shor code and discuss
its memory storage capabilities. Finally, concluding remarks and future
directions are laid out in Sec 7.

##  2 Preliminaries

###  2.1 Fractals and linear cellular automata

We begin by reviewing linear cellular automata (LCA) and their connection to
fractals. Let $c_{\mathbf{r},t}\in\\{0,1\\}$ represent the state of cell
$\mathbf{r}=(r_{1},\dots,r_{D})\in\mathbb{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,

$c_{\mathbf{r},t}=\sum_{\mathbf{r}^{\prime}}F_{\mathbf{r},\mathbf{r}^{\prime}}c_{\mathbf{r}^{\prime},t-1}\equiv\mathbf{F}\mathbf{c}_{t-1}$
(1)

where the binary matrix $\mathbf{F}$ specifies the LCA rules, and all
arithmetic is implied modulo 2. First-order means that $\mathbf{c}_{t}$ only
depends on the states $\mathbf{c}_{t-1}$ at time $t-1$, locality implies that
$\mathbf{F}$ is only non-zero for small $|\mathbf{r}-\mathbf{r}^{\prime}|$,
and translation invariance means that
$F_{\mathbf{r},\mathbf{r}^{\prime}}=f_{\mathbf{r}-\mathbf{r}^{\prime}}$ only
depends on the difference $\mathbf{r}-\mathbf{r}^{\prime}$.

Starting at time $t=0$ with the state $\mathbf{c}_{0}$, the space-time
trajectory at all future times is completely determined by

$\mathbf{c}_{t}=\mathbf{F}^{t}\mathbf{c}_{0}$ (2)

The set of cells with $c_{\mathbf{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 $\mathbf{c}_{t}$ is represented by
a $D$-variate polynomial with $\mathbb{F}_{2}$ coefficients,

$c_{t}(\mathrm{x})=\sum_{\mathbf{r}}c_{\mathbf{r},t}\mathrm{x}^{\mathbf{r}}\in\mathbb{F}_{2}[\mathrm{x}],$
(3)

where $\mathrm{x}=\\{x_{i}\\}_{i\in[1\dots D]}$, and
$\mathrm{x}^{\mathbf{r}}\equiv\prod_{i}x_{i}^{r_{i}}$ is a monomial.
Similarly, the update rule is represented by

$f(\mathrm{x})=\sum_{\mathbf{r}}f_{\mathbf{r}}\mathrm{x}^{\mathbf{r}}$ (4)

in terms of which Eqs. 1 and 2 take the particularly simple form

$c_{t}(\mathrm{x})=f(\mathrm{x})c_{t-1}(\mathrm{x})=f(\mathrm{x})^{t}c_{0}(\mathrm{x}).$
(5)

We will utilize both the vectoral representation and the polynomial
representation, depending on whichever one is most convenient.

We always assume $f(\mathrm{x})$ is chosen to have a non-zero constant term
and only positive powers of $x_{i}$, and that $f(\mathrm{x})\neq 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 $\mathbb{F}_{2}$ polynomials,
$f(\mathrm{x})^{2^{n}}=f(\mathrm{x}^{2^{n}})=f(\\{x_{i}^{2^{n}}\\})$. Assuming
the initial state $c_{0}(\mathrm{x})$ only contained finitely high powers of
any $x_{i}$, then at exponentially long times
$c_{t}(\mathrm{x})=f(\mathrm{x}^{2^{n}})c_{0}(\mathrm{x})$ describes copies of
the initial state each shifted by $2^{n}$ for each term in $f(\mathrm{x})$.

We most often consider LCA defined on finite $(L_{1},\dots,L_{d})$ systems
with periodic boundary conditions, which is enforced by the identification
$x_{i}^{L_{i}}=1$. An important property is the reversibility of an LCA rule
for a given system size, which corresponds to the existence of
$f(\mathrm{x})^{-1}$ satisfying $f(\mathrm{x})^{-1}f(\mathrm{x})=1$, after
taking into account periodic boundary conditions. A simple family of
reversible LCA is given by any polynomial $f(\mathrm{x})$ satisfying $f(1)=1$
on a system of size $L_{i}=2^{n}$. This follows from the fact that
$f(\mathrm{x})^{2^{n}}=f(\mathrm{x}^{2^{n}})=f(1)=1$, and so the inverse is
given by $f(\mathrm{x})^{-1}\equiv f(\mathrm{x})^{2^{n}-1}$.

As an example, consider $d=1$ spatial dimensions, $f(x)=1+x$ (corresponding to
$F_{\mathbf{r},\mathbf{r}^{\prime}}=\delta_{\mathbf{r},\mathbf{r}^{\prime}}+\delta_{\mathbf{r},\mathbf{r}^{\prime}+1}$)
and the initial state $c_{0}(x)=1$ (corresponding to
$c_{\mathbf{r},0}=\delta_{\mathbf{r},0}$). Then $c_{t}(x)=f(x)^{t}$, and
listing out terms for the powers of $f(x)$ we find

$\begin{matrix}f(x)^{0}&=&1&&&&&\\\ f(x)^{1}&=&1&x&&&&\\\
f(x)^{2}&=&1&&x^{2}&&&\\\ f(x)^{3}&=&1&x&x^{2}&x^{3}&&\\\
f(x)^{4}&=&1&&&&x^{4}&\\\ f(x)^{5}&=&1&x&&&x^{4}&x^{5}\\\
f(x)^{6}&=&1&&x^{2}&&x^{4}&&x^{6}\\\
f(x)^{7}&=&1&x&x^{2}&x^{3}&x^{4}&x^{5}&x^{6}&x^{7}\\\ \end{matrix}$ (6)

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}$.

###  2.2 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 [72, 73, 9]. Let us consider qubits on the sites $\mathbf{r}$ of a
hypercubic lattice. Acting on these qubits, we have Pauli operators
$X_{\mathbf{r}}$ and $Z_{\mathbf{r}}$. Define a map $\sigma_{X}$ which maps a
binary vector $a_{\mathbf{r}}$ to a tensor product of $X$ Pauli operator as

$\sigma_{X}[\mathbf{a}]=\prod_{\mathbf{r}}X_{\mathbf{r}}^{a_{\mathbf{r}}},$
(7)

and similarly $\sigma_{Z}$,

$\sigma_{Z}[\mathbf{b}]=\prod_{\mathbf{r}}Z_{\mathbf{r}}^{b_{\mathbf{r}}}.$
(8)

The commutation relation between two Pauli operators is straightforward to
compute,

$[[\sigma_{X}[\mathbf{a}],\sigma_{Z}[\mathbf{b}]]]=(-1)^{\mathbf{b}^{T}\mathbf{a}}$
(9)

where $[[A,B]]=A^{-1}B^{-1}AB$ is the group commutator.

Using the correspondence between vectors and polynomials,
$\mathbf{a}\leftrightarrow\sum_{\mathbf{r}}a_{\mathbf{r}}\mathrm{x}^{\mathbf{r}}$,
$\sigma_{X/Z}$ may also be interpreted as a map from $\mathbb{F}_{2}$
polynomials to Pauli operators. Translations are nicely expressed in this
language as multiplication by a monomial: $\sigma_{X}[a(\mathrm{x})]$
translated by $\mathbf{r}$ is
$\sigma_{X}[\mathrm{x}^{\mathbf{r}}a(\mathrm{x})]$. Note that we use the same
symbol $(\sigma_{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(\mathrm{x})$,

$[[\sigma_{X}[\mathrm{x}^{\mathbf{r}}a(\mathrm{x})],\sigma_{Z}[\mathrm{x}^{\mathbf{r}^{\prime}}b(\mathrm{x})]]]=(-1)^{c_{\mathbf{r}^{\prime}-\mathbf{r}}}$
(10)

where

$c(\mathrm{x})=\sum_{\mathbf{r}}\mathrm{x}^{\mathbf{r}}c_{\mathbf{r}}=a(\mathrm{x})b(\bar{\mathrm{x}})$
(11)

and $\bar{\mathrm{x}}\equiv\mathrm{x}^{-1}$. In particular, if
$c(\mathrm{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_{\mathbf{r}}^{(n)}$ and $Z_{\mathbf{r}}^{(n)}$ for $n=1\dots N$, is
straightforward: in the vectoral representation, $\sigma_{X/Z}$ now takes as
input $N$ vectors (or a tensor), $a_{\mathbf{r},n}$, one for each of the
qubits. Similarly, in the polynomial representation, $\sigma_{X/Z}$ takes as
input a vector of $N$ polynomials, $\mathbf{a}(\mathrm{x})=a_{n}(\mathrm{x})$.
The commutation polynomial is then given by
$c(\mathrm{x})=\mathbf{a}^{T}(\mathrm{x})\mathbf{b}(\bar{\mathrm{x}})$.

##  3 Fractalization

###  3.1 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,
${\mathrm{f}(\mathrm{y})\equiv\\{f_{i}(\mathrm{y})\\}_{i\in[1\dots d]}}$,
where $\mathrm{y}=\\{y_{j}\\}_{j\in[1\dots m]}$.

Let $A$ be a local $X$ Pauli operator in $D$ dimensions, in the polynomial
representation,

$A=\sigma_{X}\left[\mathrm{x}^{\mathbf{r}_{0}}\mathbf{a}(\mathrm{x})\right],$
(12)

where $\mathrm{x}=\\{x_{1},\dots,x_{d}\\}$ and we have chosen an “anchor
point” $\mathbf{r}_{0}$. We uniquely specify the anchor point for this
operator by maximizing each $(\mathbf{r}_{0})_{i}$ subject to the requirement
that $\mathbf{a}(\mathrm{x})$ has only positive powers of $x_{i}$. Locality
means that $\mathbf{a}(\mathrm{x})$ contains only finitely high powers of
$x_{i}$.

Going to $D+m$ dimension, let $\mathbf{s}$ denote the position vector along
the final $m$ dimensions, such that the full position of a site is
$(\mathbf{r},\mathbf{s})$. Similarly, in the polynomial notation, let
$\mathrm{y}$ denote the variables corresponding to the final $m$ dimensions.
For each $\mathbf{s}_{0}$, the fractalized operator
$A^{\mathrm{frac}}_{\mathbf{s}_{0}}$ is then defined as

$A^{\mathrm{frac}}_{\mathbf{s}_{0}}=\sigma_{X}\left[\mathrm{x}^{\mathbf{r}_{0}}\mathrm{y}^{\mathbf{s}_{0}}\mathbf{a}(\mathrm{f}(\mathrm{y})\circ\mathrm{x})\right]$
(13)

where

$\mathrm{f}(\mathrm{y})\circ\mathrm{x}\equiv\left\\{f_{i}(\mathrm{y})x_{i}\right\\}_{i\in[1\dots
D]}.$ (14)

Thus, fractalization maps a $D$-dimensional operator to a set of
$D+m$-dimensional operators,
$A\rightarrow\\{A^{\mathrm{frac}}_{\mathbf{s}_{0}}\\}$, one for each shift
$\mathbf{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,

$B=\sigma_{Z}\left[\mathrm{x}^{\mathbf{r}_{0}}\mathbf{b}(\bar{\mathrm{x}})\right],$
(15)

where the anchor point $\mathbf{r}_{0}$ is chosen similarly as before but
spatially inverted: each $(\mathbf{r}_{0})_{i}$ is minimized subject to the
constraint that $\mathbf{b}(\bar{\mathrm{x}})$ contains only positive powers
of $\bar{\mathrm{x}}$. Then, the fractalized operator is defined as

$B^{\mathrm{frac}}_{\mathbf{s}_{0}}=\sigma_{Z}\left[\mathrm{x}^{\mathbf{r}_{0}}\mathrm{y}^{\mathbf{s}_{0}}\mathbf{b}(\mathrm{f}(\bar{\mathrm{y}})\circ\bar{\mathrm{x}})\right]$
(16)

For completeness, let us also express fractalization using the vectoral
representation. The input LCAs are given by matrices
$\\{\mathbf{F}_{i}=F_{\mathbf{s},\mathbf{s}^{\prime}}^{(i)}\\}_{i\in[1\dots
D]}$. Starting with an $X$ operator

$A=\sigma_{X}\left[a_{\mathbf{r},n}\right],$ (17)

we find the anchor point $\mathbf{r}_{0}$ by maximizing each
$(\mathbf{r}_{0})_{i}$ subject to the constraint that $a_{\mathbf{r},n}=0$ for
any $(\mathbf{r}-\mathbf{r}_{0})_{i}<0$. Then, the $D+m$ dimensional
fractalized operator is

$\begin{split}A^{\mathrm{frac}}_{\mathbf{s}_{0}}&=\sigma_{X}\left[a^{\mathrm{frac}}_{(\mathbf{r},\mathbf{s}),n}\right]\\\
&=\sigma_{X}\left[\left(\prod_{i}\mathbf{F}^{(\mathbf{r}-\mathbf{r}_{0})_{i}}_{i}\right)_{\mathbf{s},\mathbf{s}_{0}}a_{\mathbf{r},n}\right].\end{split}$
(18)

where the ordering of each $\mathbf{F}_{i}$ in the product is arbitrary, as
all $\mathbf{F}_{i}$ commute due to translation invariance (this is obvious
from the polynomial representation, $f_{i}(x)f_{j}(x)=f_{j}(x)f_{i}(x)$). Note
that since $a_{\mathbf{r},n}=0$ if $(\mathbf{r}-\mathbf{r}_{0})_{i}<0$,
$\mathbf{F}_{i}$ is always raised to a non-negative power.

Similarly, for a $Z$ operator,

$B=\sigma_{Z}\left[b_{\mathbf{r},n}\right],$ (19)

we must choose the anchor $\mathbf{r}_{0}$ such that each
$(\mathbf{r}_{0})_{i}$ is minimized subject to the constraint that
$b_{\mathbf{r},n}=0$ for any $(\mathbf{r}_{0}-\mathbf{r})_{i}<0$. Then,

$B^{\mathrm{frac}}_{\mathbf{s}_{0}}=\sigma_{Z}\left[\left(\prod_{i}\mathbf{F}^{(\mathbf{r}_{0}-\mathbf{r})_{i}}_{i}\right)_{\mathbf{s}_{0},\mathbf{s}}b_{\mathbf{r},n}\right].$
(20)

Given a CSS stabilizer group
$\mathcal{S}=\langle\\{\mathcal{O}_{l}\\}\rangle$, or Hamiltonian

$H=-\sum_{l}\mathcal{O}_{l},$ (21)

where $\mathcal{O}_{l}$ are local $X$ or $Z$ generators of the stabilizer
group, we may define the fractalized stabilizer group
$\mathcal{S}^{\mathrm{frac}}=\langle\\{\mathcal{O}^{\mathrm{frac}}_{\mathbf{s},l}\\}\rangle$,
or Hamiltonian

$H^{\mathrm{frac}}=-\sum_{\mathbf{s},l}\mathcal{O}^{\mathrm{frac}}_{\mathbf{s},l}$
(22)

which, as we will show, inherits many of the properties of $H$. Note that
$H^{\mathrm{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 $i$th direction. If $f_{i}(y)^{L_{i}}=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(\mathrm{y}):q(\mathrm{y})f_{i}(\mathrm{y})^{L_{i}}=q(\mathrm{y})\\}$,
and make a choice of basis polynomials $\\{q_{j}(\mathrm{y})\\}_{j}$ for
$Q_{i}$. Then, an operator $A$ or $B$ maps on to a set of fractalized
operators, one for each $q_{j}(\mathrm{y})$. Eq 13 and 16 generalize to

$\begin{split}A^{\mathrm{frac}}_{j}&=\sigma_{X}\left[\mathrm{x}^{\mathbf{r}_{0}}q_{j}(\mathrm{y})\mathbf{a}(\mathrm{f}(\mathrm{y})\circ\mathrm{x})\right]\\\
B^{\mathrm{frac}}_{j}&=\sigma_{Z}\left[\mathrm{x}^{\mathbf{r}_{0}}q_{j}(\bar{\mathrm{y}})\mathbf{b}(\mathrm{f}(\bar{\mathrm{y}})\circ\bar{\mathrm{x}})\right].\end{split}$
(23)

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\in I$, then $q_{j}(\mathrm{y})$ must be chosen from
$\cap_{i\in I}Q_{i}$ such that the fractal is commensurate with all
directions.

###  3.2 Properties

####  3.2.1 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^{\mathrm{frac}}_{\mathbf{s}_{0}}$ will also have support on a local patch
of the $D+m$-dimensional lattice near $\mathbf{s}_{0}$. This follows from the
locality of the LCA rules: $F^{(i)}_{\mathbf{s},\mathbf{s}_{0}}$ is only non-
zero for small $|\mathbf{s}-\mathbf{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

$S=\prod_{r}X_{r}=\sigma_{X}\left[\sum_{r}x^{r}\right]$ (24)

maps on to the fractal operators

$S^{\mathrm{frac}}_{s}=\sigma_{X}\left[y^{s}\sum_{r}f(y)^{r}x^{r}\right].$
(25)

in two dimensions. Recall that $f(y)^{n}$ generates the $n$th row of the
Sierpinski triangle fractal (Eq 6). Hence, $S^{\mathrm{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}\approx 1.58$ is the Hausdorff dimension of the Sierpinski triangle.

####  3.2.2 Commutativity

The next property of fractalization is commutativity preservation. Suppose we
have two operators
${A=\sigma_{X}[\mathrm{x}^{\mathbf{r}_{0}}\mathbf{a}(\mathrm{x})]}$ and
${B=\sigma_{Z}[\mathrm{x}^{\mathbf{r}_{1}}\mathbf{b}(\bar{\mathrm{x}})]}$
satisfying $[A,B]=0$. Then,
$[A^{\mathrm{frac}}_{\mathbf{s}_{0}},B^{\mathrm{frac}}_{\mathbf{s}_{1}}]=0$
for all $\mathbf{s}_{0},\mathbf{s}_{1}$ as well. This can be seen by computing
the commutation polynomial. The commutation polynomial for $A$ and $B$ is

$c(\mathrm{x})=\mathrm{x}^{\mathbf{r}_{0}-\mathbf{r}_{1}}\mathbf{a}^{T}(\mathrm{x})\mathbf{b}(\mathrm{x}),$
(26)

which has a zero constant term iff $[A,B]=0$. The commutation polynomial for
$A^{\mathrm{frac}}_{\mathbf{s}_{0}}$ and $B^{\mathrm{frac}}_{\mathbf{s}_{1}}$
is

$\begin{split}c^{\mathrm{frac}}(\mathrm{x},\mathrm{y})&=\mathrm{x}^{\mathbf{r}_{0}-\mathbf{r}_{1}}\mathrm{y}^{\mathbf{s}_{0}-\mathbf{s}_{1}}\mathbf{a}^{T}(\mathrm{f}(\mathrm{y})\circ\mathrm{x})\mathbf{b}(\mathrm{f}(\mathrm{y})\circ\mathrm{x})\\\
&=\mathrm{y}^{\mathbf{s}_{0}-\mathbf{s}_{1}}c(\mathrm{f}(\mathrm{y})\circ\mathrm{x})\end{split}$
(27)

Therefore, if $c(x)$ has a zero constant term,
$c^{\mathrm{frac}}(\mathrm{x},\mathrm{y})$ also has a zero constant term
regardless of $\mathbf{s}_{0},\mathbf{s}_{1}$.

It is also instructive to work in the vectoral representation. Let
$AB=(-1)^{c}BA$, then, commutativity of $A=\sigma_{X}[(a_{\mathbf{r},n})]$ and
$B=\sigma_{Z}[(b_{\mathbf{r},n})]$ implies that

$c=\sum_{\mathbf{r},n}b_{\mathbf{r},n}a_{\mathbf{r},n}$ (28)

is zero. The fractalized operators satisfy
$[[A^{\mathrm{frac}}_{\mathbf{s}_{0}},B^{\mathrm{frac}}_{\mathbf{s}_{1}}]]=(-1)^{c_{\mathbf{s}_{0},\mathbf{s}_{1}}^{\mathrm{frac}}}$
where

$c^{\mathrm{frac}}_{\mathbf{s}_{0},\mathbf{s}_{1}}=\left(\prod_{i}\mathbf{F}_{i}^{(\mathbf{r}_{1}-\mathbf{r}_{0})_{i}}\right)_{\mathbf{s}_{1},\mathbf{s}_{0}}c$
(29)

which is zero if $c=0$. If instead $A$ and $B$ anticommute, then $c=1$ and
$c^{\mathrm{frac}}_{\mathbf{s}_{0},\mathbf{s}_{1}}$ may be zero or non-zero
depending on their separation $\mathbf{s}_{1}-\mathbf{s}_{0}$.

###  3.3 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
$\mathbf{F}_{i}^{L_{i}}=\mathbb{1}$, $\mathbf{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=\sigma_{X}[a_{\mathbf{r},n}]$ or $B=\sigma_{Z}[b_{\mathbf{r},n}]$, defined
in $D$ dimensions, and choose an anchor point $\mathbf{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 $\mathbf{s}$ along the perpendicular
direction. The operator $A$ or $B$ now maps on to a set of operators on each
layer separately,

$A_{\mathbf{s}_{0}}=\sigma_{X}[a_{\mathbf{r},n}\delta_{\mathbf{s},\mathbf{s}_{0}}],\;\;\;\;B_{\mathbf{s}_{0}}=\sigma_{Z}[b_{\mathbf{r},n}\delta_{\mathbf{s},\mathbf{s}_{0}}].$
(30)

on the $D+m$ dimensional lattice.

The next step involves a non-local (along the $\mathbf{s}$ directions) unitary
transformation. The unitary transformation involves a series of controlled-$X$
(CX) gates. In particular, for any invertible binary matrix $\mathbf{M}$ with
$M_{ii}=1$, there is a corresponding CX circuit which implements the
transformation

$\sigma_{X}[\mathbf{v}]\rightarrow\sigma_{X}[\mathbf{M}\mathbf{v}],\;\;\;\;\sigma_{Z}[\mathbf{w}]\rightarrow\sigma_{Z}[\mathbf{M}^{-1,T}\mathbf{w}].$
(31)

In Appendix A, we give a general algorithm for determining a CX circuit for
any such $\mathbf{M}$. Here, we choose the matrix

$M_{(\mathbf{r},\mathbf{s},n),(\mathbf{r}^{\prime},\mathbf{s}^{\prime},n^{\prime})}=\delta_{\mathbf{r},\mathbf{r}^{\prime}}\delta_{n,n^{\prime}}\left(\prod_{i}\mathbf{F}_{i}^{r_{i}}\right)_{\mathbf{s},\mathbf{s}^{\prime}}$
(32)

which maps
$A_{\mathbf{s}_{0}}\rightarrow\tilde{A}_{\mathbf{s}_{0}},\leavevmode\nobreak\
B_{\mathbf{s}_{0}}\rightarrow\tilde{B}_{\mathbf{s}_{0}}$,

$\begin{split}\tilde{A}_{\mathbf{s}_{0}}=\sigma_{X}\left[\left(\prod_{i}\mathbf{F}_{i}^{r_{i}}\right)_{\mathbf{s},\mathbf{s}_{0}}a_{\mathbf{r},n}\right],\\\
\tilde{B}_{\mathbf{s}_{0}}=\sigma_{Z}\left[\left(\prod_{i}\mathbf{F}_{i}^{-r_{i}}\right)_{\mathbf{s}_{0},\mathbf{s}}b_{\mathbf{r},n}\right].\end{split}$
(33)

which are highly non-local operators along the $\mathbf{s}$ directions as they
involve $\mathbf{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 $\langle\\{\tilde{A}_{\mathbf{s}_{0}}\\}\rangle$ and
$\langle\\{\tilde{B}_{\mathbf{s}_{0}}\\}\rangle$. For $X$ operators we choose

$\displaystyle A^{\mathrm{frac}}_{\mathbf{s}_{0}}$ $\displaystyle=$
$\displaystyle\prod_{\mathbf{s}^{\prime}}\tilde{A}_{\mathbf{s}^{\prime}}^{\left(\prod_{i}\mathbf{F}^{-r_{0,i}}_{i}\right)_{\mathbf{s}^{\prime},\mathbf{s}_{0}}}$
$\displaystyle=$
$\displaystyle\sigma_{X}\left[\sum_{\mathbf{s}^{\prime}}\left(\prod_{i}\mathbf{F}_{i}^{-r_{0,i}}\right)_{\mathbf{s}^{\prime},\mathbf{s}_{0}}\left(\prod_{i}\mathbf{F}_{i}^{r_{i}}\right)_{\mathbf{s},\mathbf{s}^{\prime}}a_{\mathbf{r},n}\right]$
$\displaystyle=$
$\displaystyle\sigma_{X}\left[\left(\prod_{i}\mathbf{F}_{i}^{(\mathbf{r}-\mathbf{r}_{0})_{i}}\right)_{\mathbf{s},\mathbf{s}_{0}}a_{\mathbf{r},n}\right]$
(35)

which matches Eq 18. Similarly, for $Z$ operators,

$\begin{split}B^{\mathrm{frac}}_{\mathbf{s}_{0}}&=\prod_{\mathbf{s}^{\prime}}\tilde{B}_{\mathbf{s}^{\prime}}^{\left(\prod_{i}\mathbf{F}^{r_{0,i}}_{i}\right)_{\mathbf{s}_{0},\mathbf{s}^{\prime}}}\end{split}$
(36)

results in the fractalized operator. Since $\mathbf{F}_{i}$ are all
invertible, this just corresponds to a different choice of generators for
their generated groups. Although $\tilde{A}_{\mathbf{s}_{0}}$ were highly non-
local, $A^{\mathrm{frac}}_{\mathbf{s}_{0}}$ are all local operators. Thus,
locality is restored.

This three-step process has important implications for fractalized quantum
codes. Consider the CSS stabilizer Hamiltonian

$H=-\sum_{l}\mathcal{O}_{l}$ (37)

with the $2^{k}$-dimensional ground state manifold obtained as the
simultaneous $+1$ eigenstate of all $\mathcal{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\rightarrow\sum_{\mathbf{s}}H_{\mathbf{s}}$ simply describes
a layered system with ground state degeneracy $2^{Lk}$. The second step is a
unitary transformation, $H\rightarrow UHU^{\dagger}$, 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 $\ket{\psi}\rightarrow U\ket{\psi}$. 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.

###  3.4 Generalizations

####  3.4.1 Beyond qubits

Everything discussed above straightforwardly generalizes from qubit degrees of
freedom to $\mathbb{Z}_{p}$ qudits of prime dimension $p$. The $X$ and $Z$
Pauli operators are replaced by $\mathbb{Z}_{p}$ clock operators, satisfying
the relation $XZ=e^{\frac{2\pi i}{p}}ZX$.

Fractalization now takes as input $p$-state LCA rules, specified by
polynomials $f(\mathrm{x})$ with $\mathbb{F}_{p}$ coefficients. Such LCA still
generate fractal structures due to the property
$f(\mathrm{x})^{p^{n}}=f(\mathrm{x}^{p^{n}})$, where self-similarity occurs at
scales $p^{n}$. Our previous discussions readily applies, with arithmetic now
done modulo $p$.

####  3.4.2 Higher order LCA

The fractalization procedure generalizes to higher order LCA. In an $n$th
order LCA, the state at time $t$, $c_{t}$, may be determined by the states at
times back to $t-n$,

$c_{t}(\mathrm{y})=\sum_{j=1\dots n}f^{(j)}(\mathrm{y})c_{t-j}(\mathrm{y})$
(38)

where the LCA update rule is now specified by a list of $n$ polynomials
$\\{f^{(j)}(\mathrm{y})\\}$.

Fractalization from $D$ to $D+m$ dimensions now takes as input $D$,
$m$-dimensional, $n$-th order LCA rules, $\\{f^{(j)}_{i}(\mathrm{y})\\}$. One
possible generalization for Eq. 18 is

$\begin{split}A^{\mathrm{frac}}_{\mathbf{s}_{0}}=&\sigma_{X}\left[\mathrm{x}^{\mathbf{r}_{0}}\mathrm{y}^{\mathbf{s}_{0}}\mathbf{a}\left(\left\\{\sum_{j}{f}^{(j)}_{i}(\mathrm{y})x_{i}^{j}\right\\}_{i}\right)\right]\end{split}$
(39)

and

$\begin{split}B^{\mathrm{frac}}_{\mathbf{s}_{0}}=&\sigma_{Z}\left[\mathrm{x}^{\mathbf{r}_{0}}\mathrm{y}^{\mathbf{s}_{0}}\mathbf{b}\left(\left\\{\sum_{j}{f}^{(j)}_{i}(\bar{\mathrm{y}})\bar{x}_{i}^{j}\right\\}_{i}\right)\right]\end{split}$
(40)

This generalization is chosen such that the fractalized Ising model (see Sec
4.1) 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_{\mathbf{r}}$ and $B_{\mathbf{r}}$ be
translations of $A$ and $B$, such that their anchor point is at $\mathbf{r}$.
Then, $[A_{\mathbf{r}_{0}},B_{\mathbf{r}_{1}}]=0$ for all $\mathbf{r}_{0/1}$,
implies that
$[A^{\mathrm{frac}}_{\mathbf{r}_{0},\mathbf{s}_{0}},B^{\mathrm{frac}}_{\mathbf{r}_{1},\mathbf{s}_{1}}]=0$
for all $\mathbf{r}_{0/1},\mathbf{s}_{0/1}$. This follows from the fact that,
for operators whose translations all commute,

$\mathbf{a}^{T}(\mathrm{x})\mathbf{b}(\mathrm{x})=\mathbf{a}^{T}(\mathrm{f}(\mathrm{y})\circ\mathrm{x})\mathbf{b}(\mathrm{f}(\mathrm{y})\circ\mathrm{x})=0\,.$
(41)

Thus, fractalization with higher order LCA should only be applied for
translation invariant commuting codes.

We also do not have a generalization of the three-step 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.

##  4 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 [59] as a translation
invariant model of glassiness, and also as a classical code [56, 57, 58]. We
then move on to the cluster model [62] and find that it fractalizes into
cluster states on higher dimensional lattices, known examples of fractal
subsystem SPT phases [55] which have been shown to be useful for universal
measurement-based quantum computation [63, 64]. We then fractalize Kitaev’s 2D
toric code model [67], which results in Yoshida’s fractal spin liquids [9], a
family of models of type-II fracton topological order.

###  4.1 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

$H_{\mathrm{1DIsing}}=-\sum_{r}Z_{r-1}Z_{r}=-\sum_{r}\sigma_{Z}[x^{r}(1+\bar{x})]$
(42)

which, after fractalization with a $1$ dimensional LCA $f(y)$, becomes

$H_{\mathrm{1DIsing}}^{\mathrm{frac}}=-\sum_{r,s}\sigma_{Z}[x^{r}y^{s}(1+f(\bar{y})\bar{x})].$
(43)

Using the choice $f(y)=1+y$, we arrive at the Newman-Moore model [59],

$H_{\mathrm{NM}}=-\sum_{r,s}Z_{(r,s)}Z_{(r-1,s)}Z_{(r-1,s-1)}.$ (44)

where $(r,s)$ labels the $x,y$ coordinates on the square lattice. The global
$\mathbb{Z}_{2}$ symmetry of the Ising model maps on to a set of fractal
symmetries,

$S=\sigma_{X}\left[\sum_{r}x^{r}\right]\rightarrow
S^{\mathrm{frac}}_{s}=\sigma_{X}\left[y^{s}\sum_{r}f(y)^{r}x^{r}\right]$ (45)

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_{\mathrm{1DIsing}}^{\mathrm{frac}}$ 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)}\\}\rangle$, where
$Z_{\mathbf{r}}|\\{z_{(r,s)}\\}\rangle=(-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_{\mathrm{1DIsing}}^{\mathrm{frac}}$ has as its ground state valid space-
time LCA trajectories, where $\mathbf{r}$ are time dimensions, and
$\mathbf{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 [57, 56, 58], in which logical bits are again represented by
the symmetry-broken ground states, of which there are now
$2^{\mathcal{O}(L^{m})}$. The fractal symmetry operators
$S^{\mathrm{frac}}_{\mathbf{s}}$ are logical operators that flip between the
symmetry-broken ground states. Upon generalizing to $\mathbb{Z}_{p}$ degrees
of freedom, this code saturates an information storage tradeoff bound [74] in
the limit $p\rightarrow\infty$ [57].

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 $\ket{\psi}$ by
flipping a segment of spins,
$W\ket{\psi}\equiv\prod_{r\in[\ell_{1}\dots\ell_{2}-1]}X_{r}\ket{\psi}$. Let
$B_{r}=Z_{r-1}Z_{r}$ be a term in $H_{\mathrm{1DIsing}}$, which satisfies the
commutation relation

$[[W,B_{r}]]=(-1)^{\delta_{r,\ell_{1}}+\delta_{r,\ell_{2}}}\,,$ (46)

making it clear that $W$ only creates two excitations at $\ell_{1}$ and
$\ell_{2}$ when acting on the ground state. From Eq. 29, the fractalized
versions of these operators satisfy
$[[W^{\mathrm{frac}},B^{\mathrm{frac}}_{r,s}]]=(-1)^{c}$ with

$c=\delta_{r,\ell_{1}}\delta_{s,s_{1}}+\delta_{r,\ell_{2}}(\mathbf{F}^{\ell_{2}-\ell_{1}})_{s,s_{1}}\,.$
(47)

Therefore, $B^{\mathrm{frac}}_{s_{1}}$ creates a single excitation at
$(\ell_{1},s_{1})$, while creating potentially many excitations along
$r=\ell_{2}$, depending on $\mathbf{F}^{\ell_{2}-\ell_{1}}$. However, if
$\ell_{2}-\ell_{1}=2^{n}$, ${f(y)^{2^{n}}=f(y^{2^{n}})}$ and so
$(\mathbf{F}^{2^{n}})_{s_{1},s}$ is only non-zero for a few $s$; for
$f(y)=1+y$, only two further excitations are created at
$(\ell_{1}+2^{n},s_{1})$ and $(\ell_{1}+2^{n},s_{1}+2^{n})$. The excitations
of $H^{\mathrm{frac}}_{\mathrm{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 (non-topological) 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 5.

###  4.2 1D Cluster model

The 1D cluster model, in CSS form, has two qubits per site and is described by
the Hamiltonian

$\begin{split}H_{\mathrm{clus}}&=-\sum_{r}X_{r}^{(1)}X_{r}^{(2)}X_{r+1}^{(1)}-\sum_{r}Z_{r-1}^{(2)}Z_{r}^{(1)}Z_{r}^{(2)}\\\
&=-\sum_{r}\sigma_{X}\left[x^{r}\begin{pmatrix}1+x\\\
1\end{pmatrix}\right]-\sum_{r}\sigma_{Z}\left[x^{r}\begin{pmatrix}1\\\
1+\bar{x}\end{pmatrix}\right].\end{split}$ (48)

After fractalization with a $1$-dimensional LCA $f(y)$,

$\begin{split}H_{\mathrm{clus}}^{\mathrm{frac}}=&-\sum_{r,s}\sigma_{X}\left[x^{r}y^{s}\begin{pmatrix}1+f(y)x\\\
1\end{pmatrix}\right]\\\
&-\sum_{r,s}\sigma_{Z}\left[x^{r}y^{s}\begin{pmatrix}1\\\
1+f(\bar{y})\bar{x}\end{pmatrix}\right],\end{split}$ (49)

which describes another cluster state on a 2D bipartite lattice. With the
choice $f(y)=1+y$, for example, $H^{\mathrm{frac}}_{\mathrm{clus}}$ describes
the cluster state on the honeycomb lattice [55]. Higher dimensional models may
be created by using $m$-dimensional LCA rules, for example, using the
$2$-dimensional LCA rule $f(\mathrm{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 symmetry-protected topological (SPT)
phase, protected by the ${\mathbb{Z}_{2}\times\mathbb{Z}_{2}}$ symmetry that
is generated by $S_{1}=\prod_{r}Z_{r}^{(1)}$ and $S_{2}=\prod_{r}X_{r}^{(2)}$.
The fractalized cluster model also describes an SPT phase, protected by the
set of fractal subsystem symmetries generated by
$S^{\mathrm{frac}}_{\alpha,s}$, $\alpha\in\\{1,2\\}$.

The non-triviality of $H^{\mathrm{frac}}_{\mathrm{clus}}$ as a fractal
subsystem SPT follows from that of $H_{\mathrm{clus}}$. The cluster chain
$H_{\mathrm{clus}}$ realizes a non-trivial SPT phase, as can be observed from
the fact that, at a boundary, the $G=\mathbb{Z}_{2}\times\mathbb{Z}_{2}$
symmetry acts as a non-trivial projective representation [75]. We may apply
fractalization directly to the open cluster chain and its symmetry operators,
following the three-step process. In the first step, the system is composed of
decoupled cluster chains along $x$, each with their own symmetry
$S_{\alpha,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 $\langle\\{S_{\alpha,s}\\}\rangle$. We can
simply choose the origin such that the anchor point for $S_{\alpha,s}$ is
$r_{0}=0$, in which case the change of generators (Eq 3.3) is trivial. Thus,
$S^{\mathrm{frac}}_{\alpha,s}$ realizes the same projective representation in
$H^{\mathrm{frac}}_{\mathrm{clus}}$ as $S_{\alpha,s}$ from the layered
$H_{\mathrm{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. 3.1, the symmetry group of
$H^{\mathrm{frac}}_{\mathrm{clus}}$ is generated by
$\\{S_{\alpha,j}^{\mathrm{frac}}\\}$, 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. [61].

For a given set of fractal symmetries, there are an infinite number of non-
trivial SPT phases [61]. 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.

###  4.3 2D toric code model

The 2D toric code [67] 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
$+\mathbf{x}_{1}$ and $+\mathbf{x}_{2}$ direction, the Hamiltonian may be
written as

$\displaystyle H_{\mathrm{TC}}=$
$\displaystyle-\sum_{\mathbf{r}}\sigma_{X}\left[\mathrm{x}^{\mathbf{r}}\begin{pmatrix}1+x_{2}\\\
1+x_{1}\end{pmatrix}\right]-\sum_{\mathbf{r}}\sigma_{Z}\left[\mathrm{x}^{\mathbf{r}}\begin{pmatrix}1+\bar{x}_{1}\\\
1+\bar{x}_{2}\end{pmatrix}\right]$ $\displaystyle=$
$\displaystyle-\sum_{\mathbf{r}}A_{\mathbf{r}}-\sum_{\mathbf{r}}B_{\mathbf{r}}\,.$
(50)

Applying fractalization with two 1-dimensional LCA rules, $f_{1}(y),f_{2}(y)$,
we get

$\begin{split}H^{\mathrm{frac}}_{\mathrm{TC}}=&-\sum_{\mathbf{r},s}\sigma_{X}\left[\mathrm{x}^{\mathbf{r}}y^{s}\begin{pmatrix}1+f_{2}(y)x_{2}\\\
1+f_{1}(y)x_{1}\end{pmatrix}\right]\\\
&-\sum_{\mathbf{r},s}\sigma_{Z}\left[\mathrm{x}^{\mathbf{r}}y^{s}\begin{pmatrix}1+f_{1}(\bar{y})\bar{x}_{1}\\\
1+f_{2}(\bar{y})\bar{x}_{2}\end{pmatrix}\right].\end{split}$ (51)

which is Yoshida’s fractal spin liquid model [9]. We take $H_{\mathrm{TC}}$ on
an $L_{1}\times L_{2}$ torus, and $H^{\mathrm{frac}}_{\mathrm{TC}}$ on an
${L_{1}\times L_{2}\times L_{3}}$ $3$-torus.

The logical operators of the toric code,

$\begin{split}\ell^{X}_{1}=\sigma_{X}\left[\sum_{r_{1}}\begin{pmatrix}x_{1}^{r_{1}}\\\
0\end{pmatrix}\right]\;\;\;\ell^{Z}_{1}=\sigma_{Z}\left[\sum_{r_{2}}\begin{pmatrix}\bar{x}_{2}^{r_{2}}\\\
0\end{pmatrix}\right]\\\
\ell^{X}_{2}=\sigma_{X}\left[\sum_{r_{2}}\begin{pmatrix}0\\\
x_{2}^{r_{2}}\end{pmatrix}\right]\;\;\;\ell^{Z}_{2}=\sigma_{Z}\left[\sum_{r_{1}}\begin{pmatrix}0\\\
\bar{x}_{1}^{r_{1}}\end{pmatrix}\right]\end{split}$ (52)

generate the $2$ qubit Pauli algebra on the ground state manifold of
$H_{\mathrm{TC}}$. If $f_{1}(y)^{L_{1}}=f_{2}(y)^{L_{2}}=1$, then they can be
straightforwardly fractalized,

$\begin{split}\ell^{X,\mathrm{frac}}_{1,s}=\sigma_{X}\left[y^{s}\sum_{r_{1}}\begin{pmatrix}f_{1}(y)^{r_{1}}x_{1}^{r_{1}}\\\
0\end{pmatrix}\right]\\\
\ell^{Z,\mathrm{frac}}_{1,s}=\sigma_{Z}\left[y^{s}\sum_{r_{2}}\begin{pmatrix}f_{2}(\bar{y})^{r_{2}}\bar{x}_{2}^{r_{2}}\\\
0\end{pmatrix}\right]\end{split}$ (53)

and similarly for $\ell^{X/Z,\mathrm{frac}}_{2,s}$. It can be easily checked
that
$[[\ell^{X,\mathrm{frac}}_{\alpha,s_{0}},\ell^{Z,\mathrm{frac}}_{\beta,s_{1}}]]=(-1)^{\delta_{\alpha\beta}\delta_{s_{0}s_{1}}}$,
thus generating a $2L_{3}$ qubit Pauli algebra on the ground state manifold.
If $f_{i}(y)^{L_{i}}\neq 1$, then some smaller subgroup is realized instead.

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^{\mathcal{O}(L)}$) and quasiparticle
excitations with restricted mobility. When the three-step process for
fractalization holds, the ground state degeneracy of
$H^{\mathrm{frac}}_{\mathrm{TC}}$ is $2^{2L_{3}}$, but will be less for
generic system sizes where $f_{i}(y)^{L_{i}}\neq 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. 4.1, 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

$f_{1}(y)^{n_{1}}=cf_{2}(y)^{n_{2}}$ (54)

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

$\begin{split}f_{1}^{(1)}(y)&=1+y+y^{2}\\\ f_{1}^{(2)}(y)&=0\\\
f_{2}^{(1)}(y)&=1+y\\\ f_{2}^{(2)}(y)&=1+y+y^{2}\end{split}$ (55)

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 $\mathbb{Z}_{2}$ symmetry of
the 2D quantum Ising model $H_{\mathrm{2DIsing}}$. After gauge fixing, the IGT
takes the form of a perturbed toric code,

$H_{\mathrm{IGT}}=-h\sum_{\mathbf{r}}A_{\mathbf{r}}-K\sum_{\mathbf{r}}B_{\mathbf{r}}-J\sum_{\mathbf{r},n}Z_{\mathbf{r}}^{(n)}-\Gamma\sum_{\mathbf{r},n}X_{\mathbf{r}}^{(n)}$
(56)

where $n=1,2$ labels the qubits on each site. $H_{\mathrm{IGT}}$ is in the
topologically ordered phase when $J/h,\Gamma/K\ll 1$. The fractalized IGT,
$H^{\mathrm{frac}}_{\mathrm{IGT}}$, can either be obtained by fractalizing
$H_{\mathrm{IGT}}$ or alternatively by applying a generalized gauging
procedure [11, 2, 76] 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. [77, 78], 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 [79, 80, 81]. Take the $L_{1}\times L_{2}\times 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 [82] 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_{\mathrm{clus}}+\lambda H_{\mathrm{coup}}$, where $H_{\mathrm{clus}}$ is
the cluster Hamiltonian on decoupled chains or layers, and $H_{\mathrm{coup}}$
couples them together. By choosing the the cluster states and couplings
appropriately, as discussed in Ref. [82], the low energy effective Hamiltonian
in the limit $\lambda\rightarrow\infty$ 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\rightarrow H^{\mathrm{frac}}$, the
cluster state wires map on to fractal cluster state layers, and in the limit
$\lambda\rightarrow\infty$, $H^{\mathrm{frac}}$ realizes the fractalized toric
code, which is exactly Yoshida’s models.

##  5 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 string-like 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 $\mathcal{O}(L)$ in
the original model to $\mathcal{O}(L^{\alpha})$, for some $1<\alpha\leq 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,

$\sum_{k=1\dots
K}\sigma_{X}[\mathbf{a}_{k}(\mathrm{x})]\equiv\sigma_{X}[\mathbf{A}(\mathrm{x})]$
(57)

where

$\mathbf{A}(\mathrm{x})=\begin{pmatrix}\mathbf{a}_{1}(\mathrm{x})&\dots&\mathbf{a}_{K}(\mathrm{x})\end{pmatrix}$
(58)

is an $N\times 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 $\sigma_{Z}[\mathbf{B}(\mathrm{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.

###  5.1 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

$\displaystyle
H_{\text{2DIsing}}=-\sum_{\textbf{r}}\sigma_{Z}[\mathrm{x}^{\mathbf{r}}\mathbf{B}(\bar{\mathrm{x}})]$
(59)

with

$\mathbf{B}(\bar{\mathrm{x}})=\begin{pmatrix}1+\bar{x}_{1}&1+\bar{x}_{2}\end{pmatrix}\,,$
(60)

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

$\displaystyle\ell^{X}=\sigma_{X}\left[\sum_{\textbf{r}}\mathrm{x}^{\mathbf{r}}\right]\,,$
$\displaystyle\ell^{Z}=\sigma_{Z}\left[1\right]\,,$ (61)

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,

$W=\sigma_{X}\left[\sum_{\mathbf{r}\in\mathcal{R}}\mathrm{x}^{\mathbf{r}}\right]\equiv\sigma_{X}[\mathbf{w}(\mathrm{x})]\,,$
(62)

where $\mathcal{R}=\mathcal{R}_{1}\times\mathcal{R}_{2}$ is a set of sites in
some rectangular region $\mathcal{R}_{i}=\\{0\dots l_{i}-1\\}$. The
excitations created in response to $W$ acting on the ground state is described
by the polynomial vector

$\mathbf{E}_{W}(\mathrm{x})=\mathbf{B}(\mathrm{x})^{T}\mathbf{w}(\mathrm{x})=\begin{pmatrix}(1+x_{1}^{l_{1}})\sum_{r_{2}\in\mathcal{R}_{2}}x_{2}^{r_{2}}\\\
(1+x_{2}^{l_{2}})\sum_{r_{1}\in\mathcal{R}_{1}}x_{1}^{r_{1}}\end{pmatrix},$
(63)

whose elements denote the locations of excited $\sigma_{Z}$ terms, and the two
rows correspond to the two types of $\sigma_{Z}$ terms. In this case, this
means that a line of $\sigma_{Z}[\mathrm{x}^{\mathbf{r}}(1+\bar{x}_{1})]$
terms are excited along the edges of the rectangle at $r_{1}=0,l_{1}$, and a
line of $\sigma_{Z}[\mathrm{x}^{\mathbf{r}}(1+\bar{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

$\displaystyle
H_{\text{2DIsing}}^{\text{frac}}=-\sum_{\textbf{r},s}\sigma_{Z}[\mathrm{x}^{\mathbf{r}}y^{s}\textbf{B}(\mathrm{f}(\bar{y})\circ\bar{\mathrm{x}})]$
(64)

for some choice of $\mathbb{F}_{2}$ polynomials $f_{1},f_{2}$.

For periodic boundary conditions $L_{1},L_{2},L_{3}$ satisfying
$f_{1}(y)^{L_{1}}=f_{2}(y)^{L_{2}}=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),

$\displaystyle\ell_{s}^{X}=\sigma_{X}\left[y^{s}\sum_{\textbf{r}}f_{1}(\bar{y})^{r_{1}}f_{2}(\bar{y})^{r_{2}}\bar{\mathrm{x}}^{\textbf{r}}\right]\,,$
$\displaystyle\ell_{s}^{Z}=\sigma_{Z}\left[y^{s}\right].$ (65)

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 $\ell^{X}_{s}$ logical
operator, i.e. a fractalized $W$ operator,
$W^{\mathrm{frac}}\equiv\sigma_{X}[\mathbf{w}(\mathrm{f}(\mathrm{y})\circ\mathrm{x})]$.
This creates excitations located at

$\mathbf{E}_{W}^{\mathrm{frac}}(\mathrm{x})=\begin{pmatrix}(1+f_{1}(y)^{l_{1}}x_{1}^{l_{1}})\sum_{r_{2}\in\mathcal{R}_{2}}f_{2}(y)^{r_{2}}x_{2}^{r_{2}}\\\
(1+f_{2}(y)^{l_{2}}x_{2}^{l_{2}})\sum_{r_{1}\in\mathcal{R}_{1}}f_{1}(y)^{r_{1}}x_{1}^{r_{1}}\end{pmatrix},$
(66)

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 non-trivially
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^{\mathrm{frac}}_{\mathrm{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 $\mathbf{l}(\bar{\mathrm{x}})$ which satisfies
$\mathbf{B}(\bar{\mathrm{x}})\mathbf{l}(\bar{\mathrm{x}})=0$. Going back to
the original 2D Ising model, the local relation

$\mathbf{l}(\bar{\mathrm{x}})=\begin{pmatrix}1+\bar{x}_{2}\\\
1+\bar{x}_{1}\end{pmatrix}\,,$ (67)

where the two columns correspond to the two Hamiltonian terms, tells us that
$\sigma_{Z}[\mathrm{x}^{\mathbf{r}}\mathbf{B}(\bar{\mathrm{x}})\mathbf{l}(\bar{\mathrm{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
$\mathbb{Z}_{2}$ loops. For the fractalized Ising model, we have the
fractalized local relation
$\mathbf{l}^{\mathrm{frac}}(\bar{\mathrm{x}},\bar{y})=\mathbf{l}(\mathrm{f}(\bar{y})\circ\bar{\mathrm{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

$\mathbf{E}_{X}(\mathrm{x},y)=\mathbf{B}^{T}(\mathrm{f}(y)\circ\mathrm{x})=\begin{pmatrix}1+f_{1}(y)x_{1}\\\
1+f_{2}(y)x_{2}\end{pmatrix}\,.$ (68)

A topological string-like excitation,
$\mathbf{E}_{\mathrm{string}}(\mathrm{x},y)$, would have to be a valid
excitation satisfying
$\mathbf{l}^{\mathrm{frac}}(\mathrm{x},y)^{T}\mathbf{E}_{\mathrm{string}}(\mathrm{x},y)=0$,
and also not be made up of a product of $\mathbf{E}_{X}(\mathrm{x},y)$. This
is equivalent to the problem of finding an $X$-type string logical operator
$\ell^{X}=\sigma_{X}[\mathbf{E}_{\mathrm{string}}(\mathrm{x},y)]$ for the
“excitation Hamiltonian” which we define in terms of the local relation and
local excitation cluster,

$H_{\mathrm{exc}}=-\sum_{\mathbf{r},s}\sigma_{Z}[\mathrm{x}^{\mathbf{r}}y^{s}\mathbf{l}^{\mathrm{frac}}(\bar{\mathrm{x}},\bar{y})]+\sigma_{X}[\mathrm{x}^{\mathbf{r}}y^{s}\mathbf{E}_{X}(\mathrm{x},y)]\\\
$ (69)

and is equivalent to the fractalized toric code $H^{\mathrm{frac}}_{TC}$ ( 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^{\mathrm{frac}}_{2D\mathrm{Ising}}$
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).

###  5.2 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

$\displaystyle
H_{\text{3DTC}}=-\sum_{\textbf{r}}\sigma_{X}[\mathrm{x}^{\textbf{r}}\textbf{A}(\mathrm{x})]+\sigma_{Z}[\mathrm{x}^{\textbf{r}}\textbf{B}(\bar{\mathrm{x}})]\,,$
(70)

where

$\displaystyle\textbf{A}(\mathrm{x})=\begin{pmatrix}0&1+x_{3}&1+x_{2}\\\
1+x_{3}&0&1+x_{1}\\\ 1+x_{2}&1+x_{1}&0\end{pmatrix},$
$\displaystyle\textbf{B}(\bar{\mathrm{x}})=\begin{pmatrix}1+\bar{x}_{1}\\\
1+\bar{x}_{2}\\\ 1+\bar{x}_{3}\end{pmatrix}.$ (71)

The logical operators for the code on periodic boundary conditions correspond
to conjugate string-membrane pairs given by

$\displaystyle\ell_{1}^{X}=\sigma_{X}\left[\sum_{r_{1}}x_{1}^{r_{1}}\mathbf{e}_{1}\right]$
$\displaystyle\ell_{1}^{Z}=\sigma_{Z}\left[\sum_{r_{2},r_{3}}\bar{x}_{2}^{r_{2}}\bar{x}_{3}^{r_{3}}\mathbf{e}_{1}\right]\,,$
(72)

where $\mathbf{e}_{1}=(1,0,0)^{T}$ is a unit vector. $\ell^{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 $\ell^{X}_{i}$ leads to
topological quasiparticle excitations at its ends, and truncating the membrane
logical operator $\ell_{i}^{Z}$ 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

$\displaystyle
H_{\text{3DTC}}^{\mathrm{frac}}=-\sum_{\textbf{r},s}\sigma_{X}[\mathrm{x}^{\textbf{r}}y^{s}\textbf{A}(\mathrm{f}(y)\circ\mathrm{x})]+\sigma_{Z}[\mathrm{x}^{\textbf{r}}y^{s}\textbf{B}(\mathrm{f}(\bar{y})\circ\bar{\mathrm{x}})]\,,$
(73)

for the A and B matrices introduced in Eq. (71) and a choice of
$\mathbb{F}_{2}$ polynomials $\mathrm{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

$\displaystyle f_{1}(y)^{n_{1}}=c\leavevmode\nobreak\
f_{2}(y)^{n_{2}}f_{3}(y)^{n_{3}}\,,$ (74)

holds for some finite constants $n_{i},c,$ where $n_{i}$ are non-negative
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)^{L_{i}}=1$ there are $3L_{4}$ logical qubits encoded within the
ground space and a basis of logical operators is given by
$\ell^{X/Z,\mathrm{frac}}_{i,s}$. To construct an explicit example we take

$\begin{split}f_{1}(y)&=1+y+y^{2}\,,\\\ f_{2}(y)&=1+y+y^{3}\,,\\\
f_{3}(y)&=1+y^{2}+y^{3}\,,\end{split}$ (75)

which satisfy the algebraic independence condition, as the chosen polynomials
are prime, and the boundary conditions when $L_{i}=2^{l}$.

###  5.3 4D toric code

Next we consider fractalizing the 4D toric code which only supports loop-like
excitations. The model is described by the Hamiltonian

$\displaystyle
H_{\text{4DTC}}=-\sum_{\textbf{r}}\sigma_{X}[\mathrm{x}^{\textbf{r}}\textbf{A}(\mathrm{x})]+\sigma_{Z}[\mathrm{x}^{\textbf{r}}\textbf{B}(\bar{\mathrm{x}})]\,,$
(76)

with

$\displaystyle\textbf{A}(\mathrm{x})$
$\displaystyle=\begin{pmatrix}1+x_{2}&1+x_{1}&0&0\\\ 1+x_{4}&0&0&1+x_{1}\\\
1+x_{3}&0&1+x_{1}&0\\\ 0&1+x_{3}&1+x_{2}&0\\\ 0&1+x_{4}&0&1+x_{2}\\\
0&0&1+x_{4}&1+x_{3}\end{pmatrix}\,,$ (77)
$\displaystyle\textbf{B}(\bar{\mathrm{x}})$
$\displaystyle=\begin{pmatrix}0&0&1+\bar{x}_{4}&1+\bar{x}_{3}\\\
0&1+\bar{x}_{3}&1+\bar{x}_{2}&0\\\ 0&1+\bar{x}_{4}&0&1+\bar{x}_{2}\\\
1+\bar{x}_{4}&0&0&1+\bar{x}_{1}\\\ 1+\bar{x}_{3}&0&1+\bar{x}_{1}&0\\\
1+\bar{x}_{2}&1+\bar{x}_{1}&0&0\end{pmatrix}\,,$ (78)

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

$\displaystyle\ell_{12}^{X}=\sigma_{X}\left[\sum_{r_{3},r_{4}}x_{3}^{r_{3}}x_{4}^{r_{4}}\mathbf{e}_{12}\right],$
$\displaystyle\ell_{12}^{Z}=\sigma_{Z}\left[\sum_{r_{1},r_{2}}\bar{x}_{1}^{r_{1}}\bar{x}_{2}^{r_{2}}\mathbf{e}_{12}\right],$
(79)

where $\mathbf{e}_{12}$ is the unit vector corresponding to the qubit on the
$(x_{1},x_{2})$ plaquette. The logical operators $\ell_{ij}^{X/Z}$ for the
remaining $ij$ are defined similarly by permuting $1,2,3,4$. $\ell^{X}_{ij}$
is a membrane operator acting in the plane orthogonal to $ij$, while
$\ell^{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 [83].

The 5D fractalized 4D toric code is described by the following Hamiltonian

$\displaystyle
H_{\text{4DTC}}^{\mathrm{frac}}=-\sum_{\textbf{r},s}\sigma_{X}[\mathrm{x}^{\textbf{r}}y^{s}\textbf{A}(\mathrm{f}(y)\circ\mathrm{x})]+\sigma_{Z}[\mathrm{x}^{\textbf{r}}y^{s}\textbf{B}(\mathrm{f}(\bar{y})\circ\bar{\mathrm{x}})]\,,$
(80)

for the A and B matrices defined in Eqs. (77) & (78), and some
$\mathbb{F}_{2}$-polynomials
$\mathrm{f}(y)=\\{f_{1}(y),f_{2}(y),f_{3}(y),f_{4}(y)\\}.$

On periodic boundary conditions $L_{i},$ for $i=1\dots 5$, such that
$f_{i}(y)^{L_{i}}=1,$ for $i=1\dots 4$, there are $6L_{5}$ encoded qubits. The
logical operator pairs $\ell^{X/Z,\mathrm{frac}}_{ij,s}$ are given by
fractalizing Eq. (79). Acting on the ground state with a truncated version of
these operators, as in earlier examples, results in the creation of a
fractalized 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\dots 4}$ satisfy a
generalized algebraic relation where either

$f_{1}(y)^{n_{1}}=c\leavevmode\nobreak\
f_{2}(y)^{n_{2}}f_{3}(y)^{n_{3}}f_{4}(y)^{n_{4}}$ (81)

or

$f_{1}(y)^{n_{1}}f_{2}(y)^{n_{2}}=c\leavevmode\nobreak\
f_{3}(y)^{n_{3}}f_{4}(y)^{n_{4}}$ (82)

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_{\mathrm{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. We note that a 5D model
without string-like excitations has appeared before [84].

To find an explicit example we can take prime polynomials $f_{i}$. For
instance

$\displaystyle f_{1}(y)=1+y+y^{2}\,,\ \,$ $\displaystyle
f_{2}(y)=1+y+y^{3}\,,$ (83) $\displaystyle f_{3}(y)=1+y^{2}+y^{3}\,,$
$\displaystyle f_{4}(y)=1+y+y^{4}\,.$ (84)

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.

###  5.4 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 membrane-like 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_{\mathrm{exc}}$ (as defined in Sec 5.1). This in turn implies the existence
of 1D (string-like) excitations in $H_{\mathrm{exc}}$, which implies the
existence of a 1D logical operator in _its_ excitation Hamiltonian
$H_{\mathrm{exc}^{2}}$. At this point, the proof of Ref [9] can be applied to
show that string logical operators exist for $H_{\mathrm{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.

##  6 Application: Fractal Bacon-Shor subsystem code

###  6.1 Fractal Bacon-Shor code

In this section we focus on a subsystem code obtained by fractalizing the
Bacon-Shor code [68] and Bravyi’s generalization thereof [70]. We investigate
the code’s quantum information storage capacity which we conjecture
asymptotically saturates an information storage tradeoff bound [70, 71].

In a subsystem code [69, 68], the full Hilbert space is structured as
$\mathcal{H}=(\mathcal{H}_{L}\otimes\mathcal{H}_{G})\oplus\mathcal{H}_{E}$,
where $\mathcal{C}=\mathcal{H}_{L}\otimes\mathcal{H}_{G}$ is the codespace,
composed of the logical and gauge subsystems, and $\mathcal{H}_{E}$ is an
error subspace. The information is stored in the state of the logical
subsystem $\mathcal{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 $\mathcal{G}$.
From $\mathcal{G}$, the stabilizer group $\mathcal{S}$ and logical operators
follow. In this section, we explain this structure for the Bacon-Shor (BS)
code [68] 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}\times L_{2}$ lattice with one
qubit per site and open boundary conditions. The gauge group $\mathcal{G}$ is
generated by products of two adjacent $X$s in the same row, or $Z$s in the
same column,

$\displaystyle\mathcal{G}=$
$\displaystyle\langle\\{X_{\mathbf{r}}X_{\mathbf{r}+\mathbf{x}},Z_{\mathbf{r}}Z_{\mathbf{r}-\mathbf{y}}\\}\rangle$
(85) $\displaystyle=$
$\displaystyle\left\langle\\{\sigma_{X}[\mathrm{x}^{\mathbf{r}}(1-x_{1})],\sigma_{Z}[\mathrm{x}^{\mathbf{r}}(1-\bar{x}_{2})]\\}\right\rangle\equiv\langle\\{A_{\mathbf{r}},B_{\mathbf{r}}\\}\rangle\,.$

The Hamiltonian model based on these generators is known as the quantum
compass model [85]. The minus signs, although irrelevant for qubits, are
included as we later generalize to $\mathbb{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
$\mathcal{G}^{\mathrm{frac}}=\langle\\{{A}_{\mathbf{r}}^{\mathrm{frac}},{B}_{\mathbf{r}}^{\mathrm{frac}}\\}\rangle$
where

$\begin{split}{A}_{\mathbf{r},s}^{\mathrm{frac}}=&\sigma_{X}\left[\mathrm{x}^{\mathbf{r}}y^{s}(1-f_{1}(y)x_{1})\right]\,,\\\
{B}_{\mathbf{r},s}^{\mathrm{frac}}=&\sigma_{Z}\left[\mathrm{x}^{\mathbf{r}}y^{s}(1-f_{2}(\bar{y})\bar{x}_{2})\right]\,,\end{split}$
(86)

on an $L_{1}\times L_{2}\times L_{3}$ system, with periodic boundary
conditions only along $L_{3}$. In the above, $\mathbf{r}=(r_{1},r_{2})$ where
$r_{i}\in R_{i}\equiv\\{0\dots L_{i}-1\\}$, and $s\in 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 $\mathcal{G}^{\mathrm{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 $\mathcal{S}$ is obtained as the center of $\mathcal{G}$,
$\mathcal{S}=C(\mathcal{G})\cap\mathcal{G}$, where $C(\cdot)$ is the
centralizer in the Pauli group. The codespace $\mathcal{C}$ is the
simultaneous $+1$ eigenspace of all operators in $\mathcal{S}$. For the BS
code, $\mathcal{S}$ is generated by the the product of $X$ along two adjacent
columns, or $Z$ along two adjacent rows,

$\begin{split}S^{X}_{r_{1}}=&\sigma_{X}\left[x_{1}^{r_{1}}(1-x_{1})\sum_{r_{2}}x_{2}^{r_{2}}\right]\\\
S^{Z}_{r_{2}}=&\sigma_{Z}\left[x_{1}^{L_{1}-1}x_{2}^{r_{2}}(1-\bar{x}_{2})\sum_{r_{1}}\bar{x}_{1}^{r_{1}}\right]\end{split}\,,$
(87)

where $S^{X}_{r_{1}}$ is defined for $r_{1}\in R_{1}\setminus\\{L_{1}-1\\}$,
and $S^{Z}_{r_{2}}$ is defined for $r_{2}\in R_{2}\setminus\\{0\\}$ (due to
the boundaries). The generators of the stabilizer group of the FBS code,
$\mathcal{S}_{\mathrm{frac}}$, is then

$\begin{split}S^{X,\mathrm{frac}}_{r_{1},s}=&\sigma_{X}\left[x_{1}^{r_{1}}y^{s}(1-f_{1}(y)x_{1})\sum_{r_{2}}f_{2}(y)^{r_{2}}x_{2}^{r_{2}}\right]\,,\\\
S^{Z,\mathrm{frac}}_{r_{2},s}=&\sigma_{Z}\left[x_{1}^{L_{1}-1}x_{2}^{r_{2}}y^{s}(1-f_{2}(\bar{y})\bar{x_{2}})\sum_{r_{1}}f_{1}(\bar{y})^{r_{1}}{\bar{x}_{1}}^{r_{1}}\right]\,.\end{split}$
(88)

for all $s\in 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,
$\mathcal{C}=\mathcal{H}_{L}\otimes\mathcal{H}_{G}$. The gauge operators act
on the codespace as $\mathbb{1}_{L}\otimes G_{G}$. Logical operators are those
that act non-trivially in $\mathcal{H}_{L}$, and come in two types: bare and
dressed logical operators. Bare logical operators act as
$\ell_{L}\otimes\mathbb{1}_{G}$. Dressed logical operators are products of
bare logical operators and gauge operators, and therefore act on the codespace
as $\ell_{L}\otimes G_{G}$. The set of bare logical operators is defined as
$L_{\mathrm{bare}}={C}(\mathcal{S})\setminus\mathcal{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

$\begin{split}\ell^{X}=\sigma_{X}\left[x_{1}^{L_{1}-1}\sum_{r_{2}}x_{2}^{r_{2}}\right]\,,\,\ell^{Z}=\sigma_{Z}\left[x_{1}^{L_{1}-1}\sum_{r_{1}}{\bar{x}_{1}}^{r_{1}}\right]\,,\end{split}$
(89)

and $\ell^{Y}=i\ell^{X}\ell^{Z}$, from which all other logical operators can
be obtained by multiplying with $\mathcal{S}$. Analogously, the corresponding
bare logical operators for the FBS code are

$\begin{split}\ell^{X,\mathrm{frac}}_{s}=\sigma_{X}\left[x_{1}^{L_{1}-1}y^{s}\sum_{r_{2}}f_{2}(y)^{r_{2}}x_{2}^{r_{2}}\right]\,,\\\
\ell^{Z,\mathrm{frac}}_{s}=\sigma_{Z}\left[x_{1}^{L_{1}-1}y^{s}\sum_{r_{1}}f_{1}(\bar{y})^{r_{1}}{\bar{x}_{1}}^{r_{1}}\right]\,,\end{split}$
(90)

of which there are now $L_{3}$ independent $X$ and $Z$ type bare logical
operators. These satisfy the commutation relation
$[[\ell^{X,\mathrm{frac}}_{s_{0}},\ell^{Z,\mathrm{frac}}_{s_{1}}]]=(-1)^{\delta_{s_{0},s_{1}}}$,
and thus generate an $L_{3}$-dimensional Pauli algebra on $\mathcal{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}\sim L$. The minimum weight of
bare logical operators of the FBS code scale as $d\sim
L^{\mathrm{min}(d_{1},d_{2})}$ where $d_{1,2}$ are the Hausdorff dimensions of
the fractal generated by $f_{1,2}$ [57]. 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\sim L^{\eta}$, where $\eta$ 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 $\eta=1$.
However, if they are algebraically unrelated, we expect that
$1<\eta\leq\mathrm{min}(d_{1},d_{2})\leq 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^{\mathrm{frac}}_{\mathbf{r},s}$ and
$B^{\mathrm{frac}}_{\mathbf{r},s}$ with $s\in R_{3}$, truncating them at the
edge. With this choice of gauge generators, the logical operators are given by
Eq. 90 with $s\in R_{3}$, also truncated at the edge. As a result, logical
operators near the boundary have reduced weight:
$\ell^{Z,\mathrm{frac}}_{0}=\prod_{r_{1}}Z_{(r_{1},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}\rightarrow L_{3}+2\delta L$
with

$\delta L=L_{1}\mathrm{deg}(f_{1})+L_{2}\mathrm{deg}(f_{2})$ (91)

where $\mathrm{deg}$ is the degree function, and use only the central $L_{3}$
logical qubits, $\\{\ell^{X/Z,\mathrm{frac}}_{s}\\}_{\delta L\leq
s<L_{3}+\delta 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.

###  6.2 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 [70, 71]

$kd^{\frac{1}{D-1}}\leq O(n)$ (92)

(Ref. [70] only derived the bound in $D=2$, but the proof can be extended to
higher dimensions as written down in Ref. [71]). Bravyi’s generalization [70]
of the BS code, henceforth referred to as the Bravyi-Bacon-Shor (BBS) code,
saturates this bound in 2D.

The 2D BBS code is parameterized by a binary matrix $K_{\mathbf{r}}$, which
specifies a set of sites $K_{\mathbf{r}}=0$ to be removed from the 2D lattice.
The gauge generators are modified accordingly to skip over the removed sites,
i.e. $X_{\mathbf{r}}X_{\mathbf{r}+n\mathbf{x}}$, where
$\mathbf{r}$,$\mathbf{r}+n\mathbf{x}$ 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 [70]. The end
result is that logical operators from different rows or columns can no longer
generically be related to one another by $\mathcal{S}$, and the number of
logical qubits is given by $k=\mathrm{rank}(K)$. Furthermore, it is possible
[70] to find a family of matrices $K_{\mathbf{r}}$ such that $k\sim L$, while
maintaining $n\sim L^{2}$ and $d\sim L$. Thus, $[n,k,d]\sim[L^{2},L,L]$
saturates the 2D bound in Eq. 92, $kd\sim L^{2}\leq\mathcal{O}(L^{2})$. In 3D,
however, the current best known subsystem code does not quite saturate the
scaling in the tradeoff bound [86].

The FBS code has $n=L^{3}$, $k=L$, and $d\sim L^{\eta}$ for some
$1\leq\eta\leq 2$ as discussed previously. This is still far from saturating
the 3D tradeoff bound, $k\sqrt{d}\sim L^{1+\eta}\leq\mathcal{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_{\mathbf{r}}=0$ from the original model corresponds
to removal of a whole line of sites $\\{(\mathbf{r},s)\\}_{s}$ from the
fractal model. The number of logical qubits is improved to
$k=L\,\mathrm{rank}(K)$. Again, there is a choice of binary matrices such that
the parameters scale as $[n,k,d]\sim[L^{3},L^{2},L^{\eta}]$ (for a potentially
different $\eta$). Thus, the FBBS code nearly saturates the tradeoff bound
depending on the value of $\eta$, $k\sqrt{d}\sim L^{2+\eta}\leq 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}\rightarrow 2$ (this was used by Yoshida [57] to
asymptotically saturate the information storage bound in a classical spin
system). We conjecture that in the limit $p\rightarrow\infty$, for
$f_{1},f_{2}$ algebraically unrelated polynomials, the code distance scaling
parameter also approaches $\eta\rightarrow 2$. The FBBS code is therefore
conjectured to asymptotically saturate the 3D information storage bound
$k\sqrt{d}\sim 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\rightarrow 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 3.3. 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\sim L^{\eta}$, where
$1\leq\eta\leq m+1$. Thus, starting with the 2D BBS code, the $2+m$
dimensional FBBS code has parameters
$[n,k,d]\sim[L^{m+2},L^{m+1},L^{1\leq\eta\leq m+1}]$. The analogue of our
conjecture in higher dimensions would be that, for certain choices of LCA
rules and $\mathbb{Z}_{p}$ qudits, $\eta\rightarrow m+1$ in the limit
$p\rightarrow\infty$. In this case, the tradeoff bound

$kd^{\frac{1}{m+1}}\sim L^{m+1+\eta/(m+1)}\xrightarrow[\eta\rightarrow
m+1]{p\rightarrow\infty}L^{m+2}\leq O(L^{m+2})$ (93)

is saturated in the large $p$ limit in any dimension $m+2$.

##  7 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 [70] and in the limit of large on-site
Hilbert space dimension [57], 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 non-local 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 [87] 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 higher-order
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.

## Acknowledgements

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

##  Appendix A CX circuit for arbitrary M

In this section, we present an algorithm for constructing a CX circuit $U$
that implements the transformation

$U\sigma_{X}[\mathbf{v}]U^{\dagger}=\sigma_{X}[\mathbf{M}\mathbf{v}],\;\;\;\;U\sigma_{Z}[\mathbf{w}]U^{\dagger}=\sigma_{Z}[\mathbf{M}^{-1,T}\mathbf{w}].$
(94)

for an arbitrary invertible binary matrix
$\mathbf{M}\in\mathbb{F}_{2}^{L\times L}$ with $M_{ii}=1$. We will say that
$U$ is the unitary corresponding to the matrix operation $\mathbf{M}$.

A single $\mathrm{CX}_{ij}$ gate ($i\neq j$), defined as

$\mathrm{CX}_{ij}=(Z_{i}-1)/2+(Z_{i}+1)X_{j}/2\,,$ (95)

corresponds to the matrix operation

$\mathbf{T}_{ij}=\mathbb{1}+\mathbf{e}_{i}\mathbf{e}_{j}^{T}$ (96)

where $\mathbf{e}_{i}$ is the $i$th unit vector. Thus, the problem of finding
the CX circuit for $\mathbf{M}$ is equivalent to finding a sequence of pairs
$(i_{n},j_{n})$, $n=1\dots N$, that satisfies

$\mathbf{M}=\mathbf{T}_{i_{1}j_{1}}\dots\mathbf{T}_{i_{N}j_{N}}$ (97)

corresponding to the CX circuit

$U=\mathrm{CX}_{i_{1}j_{1}}\dots\mathrm{CX}_{i_{N}j_{N}}.$ (98)

We can construct $\mathbf{M}$ iteratively by noticing that

$(\mathbf{T}_{ij}\mathbf{M})_{kl}=M_{kl}+\delta_{ik}M_{jl}$ (99)

corresponds to a row action on $\mathbf{M}$ where the $j$th row of
$\mathbf{M}$ is added to the $i$th row. As $\mathbf{M}$ is invertible with
unit diagonal, a series of such row operations can be used to transform
$\mathbf{M}$ into the identity matrix. The sequence of row operations
$(i_{n},j_{n})$ may be generated via the following algorithm,

  * Start with $\mathbf{M}$ and $n=1$.

  * For each $j=1\dots L$,

    * For each $i=1\dots L$, except $i=j$,

      * If $M_{ij}=1$, set $(i_{n},j_{n})=(i,j)$,

$\mathbf{M}\rightarrow\mathbf{M}^{\prime}=\mathbf{T}_{ij}\mathbf{M},$

and then increment $n\rightarrow n+1$.

By the end of this algorithm, $\mathbf{M}$ will have been transformed into the
identity matrix. Thus, we have generate a list of pairs $(i_{n},j_{n})$ such
that

$\mathbf{T}_{i_{N},j_{N}}\dots\mathbf{T}_{i_{1},j_{1}}\mathbf{M}=\mathbb{1}$
(100)

Since $\mathbf{T}_{ij}^{2}=\mathbb{1}$, this can be rewritten

$\mathbf{M}=\mathbf{T}_{i_{1},j_{1}}\dots\mathbf{T}_{i_{N},j_{N}}$ (101)

which matches Eq. 97 and we have therefore found a CX circuit corresponding to
$\mathbf{M}$.

## References

  * Haah [2016a] Jeongwan Haah.  An Invariant of Topologically Ordered States Under Local Unitary Transformations.  _Commun. Math. Phys._ , 342(3):771–801, 2016a.  ISSN 14320916.  doi: 10.1007/s00220-016-2594-y.
  * Vijay et al. [2016] Sagar Vijay, Jeongwan Haah, and Liang Fu.  Fracton topological order, generalized lattice gauge theory, and duality.  _Phys. Rev. B_ , 94:235157, Dec 2016.  doi: 10.1103/PhysRevB.94.235157.  URL https://link.aps.org/doi/10.1103/PhysRevB.94.235157.
  * Nandkishore and Hermele [2019] Rahul M. Nandkishore and Michael Hermele.  Fractons.  _Annual Review of Condensed Matter Physics_ , 10(1):295–313, 2019.  doi: 10.1146/annurev-conmatphys-031218-013604.  URL https://doi.org/10.1146/annurev-conmatphys-031218-013604.
  * Pretko et al. [2020] Michael Pretko, Xie Chen, and Yizhi You.  Fracton phases of matter.  _International Journal of Modern Physics A_ , 35(06):2030003, 2020.  doi: 10.1142/S0217751X20300033.  URL https://doi.org/10.1142/S0217751X20300033.
  * Castelnovo and Chamon [2012] Claudio Castelnovo and Claudio Chamon.  Topological quantum glassiness.  _Philosophical Magazine_ , 92(1-3):304–323, 2012\.  doi: 10.1080/14786435.2011.609152.  URL https://doi.org/10.1080/14786435.2011.609152.
  * Bravyi et al. [2011] Sergey Bravyi, Bernhard Leemhuis, and Barbara M. Terhal.  Topological order in an exactly solvable 3d spin model.  _Annals of Physics_ , 326(4):839 – 866, 2011\.  ISSN 0003-4916.  doi: https://doi.org/10.1016/j.aop.2010.11.002.  URL http://www.sciencedirect.com/science/article/pii/S0003491610001910.
  * Haah [2011] Jeongwan Haah.  Local stabilizer codes in three dimensions without string logical operators.  _Phys. Rev. A_ , 83:042330, Apr 2011.  doi: 10.1103/PhysRevA.83.042330.  URL https://link.aps.org/doi/10.1103/PhysRevA.83.042330.
  * Kim [2012] Isaac H. Kim.  3d local qupit quantum code without string logical operator.  2012\.  URL https://arxiv.org/abs/1202.0052.
  * Yoshida [2013] Beni Yoshida.  Exotic topological order in fractal spin liquids.  _Phys. Rev. B_ , 88:125122, Sep 2013.  doi: 10.1103/PhysRevB.88.125122.  URL https://link.aps.org/doi/10.1103/PhysRevB.88.125122.
  * Vijay et al. [2015] Sagar Vijay, Jeongwan Haah, and Liang Fu.  A new kind of topological quantum order: A dimensional hierarchy of quasiparticles built from stationary excitations.  _Phys. Rev. B_ , 92:235136, Dec 2015.  doi: 10.1103/PhysRevB.92.235136.  URL https://link.aps.org/doi/10.1103/PhysRevB.92.235136.
  * Williamson [2016] Dominic J. Williamson.  Fractal symmetries: Ungauging the cubic code.  _Phys. Rev. B_ , 94(15):155128, oct 2016.  ISSN 24699969.  doi: 10.1103/PhysRevB.94.155128.  URL https://link.aps.org/doi/10.1103/PhysRevB.94.155128.
  * Dua et al. [2019a] Arpit Dua, Isaac H. Kim, Meng Cheng, and Dominic J. Williamson.  Sorting topological stabilizer models in three dimensions.  _Phys. Rev. B_ , 100:155137, Oct 2019a.  doi: 10.1103/PhysRevB.100.155137.  URL https://link.aps.org/doi/10.1103/PhysRevB.100.155137.
  * Bravyi and Haah [2011] Sergey Bravyi and Jeongwan Haah.  Energy landscape of 3D spin hamiltonians with topological order.  _Phys. Rev. Lett._ , 107(15):150504, oct 2011\.  ISSN 00319007.  doi: 10.1103/PhysRevLett.107.150504.  URL https://link.aps.org/doi/10.1103/PhysRevLett.107.150504.
  * Bravyi and Haah [2013] Sergey Bravyi and Jeongwan Haah.  Quantum self-correction in the 3D cubic code model.  _Phys. Rev. Lett._ , 111(20):200501, nov 2013\.  ISSN 00319007.  doi: 10.1103/PhysRevLett.111.200501.  URL https://link.aps.org/doi/10.1103/PhysRevLett.111.200501.
  * Brown et al. [2016] Benjamin J. Brown, Daniel Loss, Jiannis K. Pachos, Chris N. Self, and James R. Wootton.  Quantum memories at finite temperature.  _Rev. Mod. Phys._ , 88(4):45005, nov 2016.  ISSN 15390756.  doi: 10.1103/RevModPhys.88.045005.  URL https://link.aps.org/doi/10.1103/RevModPhys.88.045005.
  * Prem et al. [2017] Abhinav Prem, Jeongwan Haah, and Rahul Nandkishore.  Glassy quantum dynamics in translation invariant fracton models.  _Phys. Rev. B_ , 95(15):155133, apr 2017.  ISSN 24699969.  doi: 10.1103/PhysRevB.95.155133.  URL https://link.aps.org/doi/10.1103/PhysRevB.95.155133.
  * Pai et al. [2019] Shriya Pai, Michael Pretko, and Rahul M. Nandkishore.  Localization in fractonic random circuits.  _Phys. Rev. X_ , 9:021003, Apr 2019.  doi: 10.1103/PhysRevX.9.021003.  URL https://link.aps.org/doi/10.1103/PhysRevX.9.021003.
  * Gromov et al. [2020] Andrey Gromov, Andrew Lucas, and Rahul M. Nandkishore.  Fracton hydrodynamics.  _Phys. Rev. Research_ , 2:033124, Jul 2020.  doi: 10.1103/PhysRevResearch.2.033124.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033124.
  * Pai and Pretko [2019] Shriya Pai and Michael Pretko.  Dynamical scar states in driven fracton systems.  _Phys. Rev. Lett._ , 123:136401, Sep 2019.  doi: 10.1103/PhysRevLett.123.136401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.123.136401.
  * He et al. [2020] Huan He, Yizhi You, and Abhinav Prem.  Lieb-schultz-mattis–type constraints on fractonic matter.  _Phys. Rev. B_ , 101:165145, Apr 2020.  doi: 10.1103/PhysRevB.101.165145.  URL https://link.aps.org/doi/10.1103/PhysRevB.101.165145.
  * Yan [2019] Han Yan.  Hyperbolic fracton model, subsystem symmetry, and holography.  _Phys. Rev. B_ , 99:155126, Apr 2019.  doi: 10.1103/PhysRevB.99.155126.  URL https://link.aps.org/doi/10.1103/PhysRevB.99.155126.
  * Pretko and Radzihovsky [2018] Michael Pretko and Leo Radzihovsky.  Fracton-Elasticity Duality.  _Phys. Rev. Lett._ , 120(19):195301, may 2018\.  ISSN 10797114.  doi: 10.1103/PhysRevLett.120.195301.  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.195301.
  * Gromov [2019a] Andrey Gromov.  Chiral Topological Elasticity and Fracton Order.  _Phys. Rev. Lett._ , 122(7), dec 2019a.  ISSN 10797114.  doi: 10.1103/PhysRevLett.122.076403.  URL http://arxiv.org/abs/1712.06600.
  * Pretko et al. [2019] Michael Pretko, Zhengzheng Zhai, and Leo Radzihovsky.  Crystal-to-fracton tensor gauge theory dualities.  _Phys. Rev. B_ , 100:134113, Oct 2019.  doi: 10.1103/PhysRevB.100.134113.  URL https://link.aps.org/doi/10.1103/PhysRevB.100.134113.
  * Gromov and Surówka [2019] Andrey Gromov and Piotr Surówka.  On duality between Cosserat elasticity and fractons.  _SciPost Phys._ , 8(4), aug 2019.  doi: 10.21468/SciPostPhys.8.4.065.  URL http://dx.doi.org/10.21468/SciPostPhys.8.4.065.
  * Halász et al. [2017] Gábor B. Halász, Timothy H. Hsieh, and Leon Balents.  Fracton Topological Phases from Strongly Coupled Spin Chains.  _Phys. Rev. Lett._ , 119(25):257202, dec 2017\.  ISSN 10797114.  doi: 10.1103/PhysRevLett.119.257202.  URL https://link.aps.org/doi/10.1103/PhysRevLett.119.257202.
  * Prem et al. [2018] Abhinav Prem, Sagar Vijay, Yang-Zhi Chou, Michael Pretko, and Rahul M. Nandkishore.  Pinch point singularities of tensor spin liquids.  _Phys. Rev. B_ , 98:165140, Oct 2018.  doi: 10.1103/PhysRevB.98.165140.  URL https://link.aps.org/doi/10.1103/PhysRevB.98.165140.
  * Doshi and Gromov [2020] Darshil Doshi and Andrey Gromov.  Vortices and Fractons.  may 2020.  URL http://arxiv.org/abs/2005.03015.
  * Pretko [2017] Michael Pretko.  Emergent gravity of fractons: Mach’s principle revisited.  _Phys. Rev. D_ , 96(2):24051, jul 2017.  ISSN 24700029.  doi: 10.1103/PhysRevD.96.024051.  URL https://link.aps.org/doi/10.1103/PhysRevD.96.024051.
  * Brown and Williamson [2020] Benjamin J. Brown and Dominic J. Williamson.  Parallelized quantum error correction with fracton topological codes.  _Phys. Rev. Research_ , 2:013303, Mar 2020.  doi: 10.1103/PhysRevResearch.2.013303.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013303.
  * Vijay and Fu [2017] Sagar Vijay and Liang Fu.  A Generalization of Non-Abelian Anyons in Three Dimensions.  2017\.  URL http://arxiv.org/abs/1706.07070.
  * Williamson and Cheng [2020] Dominic J. Williamson and Meng Cheng.  Designer non-Abelian fractons from topological layers.  apr 2020.  URL http://arxiv.org/abs/2004.07251.
  * Bulmash and Barkeshli [2019] Daniel Bulmash and Maissam Barkeshli.  Gauging fractons: Immobile non-abelian quasiparticles, fractals, and position-dependent degeneracies.  _Phys. Rev. B_ , 100:155146, Oct 2019.  doi: 10.1103/PhysRevB.100.155146.  URL https://link.aps.org/doi/10.1103/PhysRevB.100.155146.
  * Prem and Williamson [2019] Abhinav Prem and Dominic Williamson.  Gauging permutation symmetries as a route to non-Abelian fractons.  _SciPost Phys._ , 7(5):068, nov 2019.  ISSN 2542-4653.  doi: 10.21468/scipostphys.7.5.068.  URL https://scipost.org/10.21468/SciPostPhys.7.5.068.
  * Stephen et al. [2020] David T. Stephen, José Garre-Rubio, Arpit Dua, and Dominic J. Williamson.  Subsystem symmetry enriched topological order in three dimensions.  _Phys. Rev. Research_ , 2:033331, Aug 2020.  doi: 10.1103/PhysRevResearch.2.033331.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033331.
  * Pai and Hermele [2019] Shriya Pai and Michael Hermele.  Fracton fusion and statistics.  _Phys. Rev. B_ , 100:195136, Nov 2019.  doi: 10.1103/PhysRevB.100.195136.  URL https://link.aps.org/doi/10.1103/PhysRevB.100.195136.
  * Haah [2014] Jeongwan Haah.  Bifurcation in entanglement renormalization group flow of a gapped spin model.  _Phys. Rev. B_ , 89:075119, Feb 2014.  doi: 10.1103/PhysRevB.89.075119.  URL https://link.aps.org/doi/10.1103/PhysRevB.89.075119.
  * Shirley et al. [2018] Wilbur Shirley, Kevin Slagle, Zhenghan Wang, and Xie Chen.  Fracton models on general three-dimensional manifolds.  _Phys. Rev. X_ , 8:031051, Aug 2018.  doi: 10.1103/PhysRevX.8.031051.  URL https://link.aps.org/doi/10.1103/PhysRevX.8.031051.
  * Dua et al. [2020] Arpit Dua, Pratyush Sarkar, Dominic J. Williamson, and Meng Cheng.  Bifurcating entanglement-renormalization group flows of fracton stabilizer models.  _Phys. Rev. Research_ , 2:033021, Jul 2020.  doi: 10.1103/PhysRevResearch.2.033021.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033021.
  * Slagle et al. [2019a] Kevin Slagle, Abhinav Prem, and Michael Pretko.  Symmetric tensor gauge theories on curved spaces.  _Annals of Physics_ , 410:167910, 2019a.  ISSN 0003-4916.  doi: https://doi.org/10.1016/j.aop.2019.167910.  URL http://www.sciencedirect.com/science/article/pii/S0003491619301654.
  * Slagle et al. [2019b] Kevin Slagle, David Aasen, and Dominic Williamson.  Foliated Field Theory and String-Membrane-Net Condensation Picture of Fracton Order.  _SciPost Phys._ , 6:43, 2019b.  doi: 10.21468/SciPostPhys.6.4.043.  URL https://scipost.org/10.21468/SciPostPhys.6.4.043.
  * Wen [2020] Xiao-Gang Wen.  Systematic construction of gapped nonliquid states.  _Phys. Rev. Research_ , 2:033300, Aug 2020.  doi: 10.1103/PhysRevResearch.2.033300.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033300.
  * Aasen et al. [2020] David Aasen, Daniel Bulmash, Abhinav Prem, Kevin Slagle, and Dominic J. Williamson.  Topological defect networks for fractons of all types.  _Phys. Rev. Research_ , 2:043165, Oct 2020.  doi: 10.1103/PhysRevResearch.2.043165.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043165.
  * Wang [2020] Juven Wang.  Non-liquid cellular states.  2020\.
  * Bulmash and Barkeshli [2018] Daniel Bulmash and Maissam Barkeshli.  Generalized $U(1)$ Gauge Field Theories and Fractal Dynamics.  _arXiv_ , (1):1–7, 2018.  URL http://arxiv.org/abs/1806.01855.
  * Gromov [2019b] Andrey Gromov.  Towards Classification of Fracton Phases: The Multipole Algebra.  _Phys. Rev. X_ , 9(3), dec 2019b.  ISSN 21603308.  doi: 10.1103/PhysRevX.9.031035.  URL https://arxiv.org/pdf/1812.05104.pdf.
  * Pretko [2018] Michael Pretko.  The fracton gauge principle.  _Phys. Rev. B_ , 98:115134, Sep 2018.  doi: 10.1103/PhysRevB.98.115134.  URL https://link.aps.org/doi/10.1103/PhysRevB.98.115134.
  * Williamson et al. [2019a] Dominic J. Williamson, Zhen Bi, and Meng Cheng.  Fractonic matter in symmetry-enriched $u(1)$ gauge theory.  _Phys. Rev. B_ , 100:125150, Sep 2019a.  doi: 10.1103/PhysRevB.100.125150.  URL https://link.aps.org/doi/10.1103/PhysRevB.100.125150.
  * Seiberg [2020] Nathan Seiberg.  Field Theories With a Vector Global Symmetry.  _SciPost Phys._ , 8:50, 2020.  doi: 10.21468/SciPostPhys.8.4.050.  URL https://scipost.org/10.21468/SciPostPhys.8.4.050.
  * Tantivasadakarn and Vijay [2020] Nathanan Tantivasadakarn and Sagar Vijay.  Searching for fracton orders via symmetry defect condensation.  _Phys. Rev. B_ , 101:165143, Apr 2020.  doi: 10.1103/PhysRevB.101.165143.  URL https://link.aps.org/doi/10.1103/PhysRevB.101.165143.
  * Shirley [2020] Wilbur Shirley.  Fractonic order and emergent fermionic gauge theory.  feb 2020.  URL http://arxiv.org/abs/2002.12026.
  * Tantivasadakarn [2020] Nathanan Tantivasadakarn.  Jordan-wigner dualities for translation-invariant hamiltonians in any dimension: Emergent fermions in fracton topological order.  _Phys. Rev. Research_ , 2:023353, Jun 2020.  doi: 10.1103/PhysRevResearch.2.023353.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.023353.
  * Devakul et al. [2020] Trithep Devakul, Wilbur Shirley, and Juven Wang.  Strong planar subsystem symmetry-protected topological phases and their dual fracton orders.  _Phys. Rev. Research_ , 2:012059, Mar 2020.  doi: 10.1103/PhysRevResearch.2.012059.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.012059.
  * Devakul et al. [2018] Trithep Devakul, S. A. Parameswaran, and S. L. Sondhi.  Correlation function diagnostics for type-i fracton phases.  _Phys. Rev. B_ , 97:041110(R), Jan 2018.  doi: 10.1103/PhysRevB.97.041110.  URL https://link.aps.org/doi/10.1103/PhysRevB.97.041110.
  * Devakul et al. [2019] Trithep Devakul, Yizhi You, F. J. Burnell, and S. L. Sondhi.  Fractal Symmetric Phases of Matter.  _SciPost Phys._ , 6:7, 2019.  doi: 10.21468/SciPostPhys.6.1.007.  URL https://scipost.org/10.21468/SciPostPhys.6.1.007.
  * Chowdhury et al. [1994] D. R. Chowdhury, S. Basu, I. S. Gupta, and P. P. Chaudhuri.  Design of caecc - cellular automata based error correcting code.  _IEEE Transactions on Computers_ , 43(6):759–764, 1994.  doi: 10.1109/12.286310.
  * Yoshida [2011] Beni Yoshida.  Information storage capacity of discrete spin systems.  _Ann. Phys. (N. Y)._ , 338:134–166, nov 2011.  doi: 10.1016/j.aop.2013.07.009.  URL http://arxiv.org/abs/1111.3275http://dx.doi.org/10.1016/j.aop.2013.07.009.
  * Nixon and Brown [2021] G. M. Nixon and B. J. Brown.  Correcting spanning errors with a fractal code.  _IEEE Transactions on Information Theory_ , pages 1–1, 2021.  doi: 10.1109/TIT.2021.3068359.
  * Newman and Moore [1999] M. E. J. Newman and Cristopher Moore.  Glassy dynamics and aging in an exactly solvable spin model.  _Phys. Rev. E_ , 60:5068–5072, Nov 1999.  doi: 10.1103/PhysRevE.60.5068.  URL https://link.aps.org/doi/10.1103/PhysRevE.60.5068.
  * Kubica and Yoshida [2018] Aleksander Kubica and Beni Yoshida.  Ungauging quantum error-correcting codes.  2018\.  URL http://arxiv.org/abs/1805.01836.
  * Devakul [2019] Trithep Devakul.  Classifying local fractal subsystem symmetry-protected topological phases.  _Phys. Rev. B_ , 99:235131, Jun 2019.  doi: 10.1103/PhysRevB.99.235131.  URL https://link.aps.org/doi/10.1103/PhysRevB.99.235131.
  * Briegel and Raussendorf [2001] Hans J. Briegel and Robert Raussendorf.  Persistent entanglement in arrays of interacting particles.  _Phys. Rev. Lett._ , 86(5):910–913, jan 2001\.  ISSN 00319007.  doi: 10.1103/PhysRevLett.86.910.  URL https://link.aps.org/doi/10.1103/PhysRevLett.86.910.
  * Raussendorf et al. [2019] Robert Raussendorf, Cihan Okay, Dong-Sheng Wang, David T. Stephen, and Hendrik Poulsen Nautrup.  Computationally universal phase of quantum matter.  _Phys. Rev. Lett._ , 122:090501, Mar 2019.  doi: 10.1103/PhysRevLett.122.090501.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.090501.
  * Devakul and Williamson [2018] Trithep Devakul and Dominic J. Williamson.  Universal quantum computation using fractal symmetry-protected cluster phases.  _Phys. Rev. A_ , 98(2):022332, aug 2018.  ISSN 24699934.  doi: 10.1103/PhysRevA.98.022332.  URL https://link.aps.org/doi/10.1103/PhysRevA.98.022332http://arxiv.org/abs/1806.04663.
  * Stephen et al. [2019] David T. Stephen, Hendrik Poulsen Nautrup, Juani Bermejo-Vega, Jens Eisert, and Robert Raussendorf.  Subsystem symmetries, quantum cellular automata, and computational phases of quantum matter.  _Quantum_ , 3:142, May 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-05-20-142.  URL https://doi.org/10.22331/q-2019-05-20-142.
  * Daniel et al. [2020] Austin K. Daniel, Rafael N. Alexander, and Akimasa Miyake.  Computational universality of symmetry-protected topologically ordered cluster phases on 2D Archimedean lattices.  _Quantum_ , 4:228, February 2020.  ISSN 2521-327X.  doi: 10.22331/q-2020-02-10-228.  URL https://doi.org/10.22331/q-2020-02-10-228.
  * Kitaev [2003] A. Yu Kitaev.  Fault-tolerant quantum computation by anyons.  _Ann. Phys. (N. Y)._ , 303(1):2–30, jan 2003\.  ISSN 00034916.  doi: 10.1016/S0003-4916(02)00018-0.  URL http://www.sciencedirect.com/science/article/pii/S0003491602000180.
  * Bacon [2006] Dave Bacon.  Operator quantum error-correcting subsystems for self-correcting quantum memories.  _Phys. Rev. A_ , 73:012340, Jan 2006.  doi: 10.1103/PhysRevA.73.012340.  URL https://link.aps.org/doi/10.1103/PhysRevA.73.012340.
  * Poulin [2005] David Poulin.  Stabilizer formalism for operator quantum error correction.  _Phys. Rev. Lett._ , 95:230504, Dec 2005.  doi: 10.1103/PhysRevLett.95.230504.  URL https://link.aps.org/doi/10.1103/PhysRevLett.95.230504.
  * Bravyi [2011] Sergey Bravyi.  Subsystem codes with spatially local generators.  _Phys. Rev. A_ , 83:012320, Jan 2011.  doi: 10.1103/PhysRevA.83.012320.  URL https://link.aps.org/doi/10.1103/PhysRevA.83.012320.
  * Flammia et al. [2016] Steven T. Flammia, Jeongwan Haah, Michael J. Kastoryano, and Isaac H. Kim.  Limits on the storage of quantum information in a volume of space.  _Quantum_ , 1:4, oct 2016.  doi: 10.22331/q-2017-04-25-4.  URL http://dx.doi.org/10.22331/q-2017-04-25-4.
  * MacWilliams and Sloane [1977] F. J. MacWilliams and N. J. A. Sloane.  _The theory of error correcting codes_.  North-Holland Pub. Co. ; sole distributors for the U.S.A. and Canada, Elsevier/North-Holland Amsterdam ; New York : New York, 1977.  ISBN 0444850090 0444850104.
  * Haah [2013] Jeongwan Haah.  Commuting Pauli Hamiltonians as Maps between Free Modules.  _Commun. Math. Phys._ , 324(2):351–399, 2013\.  ISSN 00103616.  doi: 10.1007/s00220-013-1810-2.  URL http://dx.doi.org/10.1007/s00220-013-1810-2.
  * Bravyi et al. [2010] Sergey Bravyi, David Poulin, and Barbara Terhal.  Tradeoffs for reliable quantum information storage in 2d systems.  _Phys. Rev. Lett._ , 104:050503, Feb 2010.  doi: 10.1103/PhysRevLett.104.050503.  URL https://link.aps.org/doi/10.1103/PhysRevLett.104.050503.
  * Chen et al. [2011] Xie Chen, Zheng-Cheng Gu, and Xiao-Gang Wen.  Complete classification of one-dimensional gapped quantum phases in interacting spin systems.  _Phys. Rev. B_ , 84:235128, Dec 2011.  doi: 10.1103/PhysRevB.84.235128.  URL https://link.aps.org/doi/10.1103/PhysRevB.84.235128.
  * Shirley et al. [2019] Wilbur Shirley, Kevin Slagle, and Xie Chen.  Foliated fracton order from gauging subsystem symmetries.  _SciPost Phys._ , 6:41, 2019.  doi: 10.21468/SciPostPhys.6.4.041.  URL https://scipost.org/10.21468/SciPostPhys.6.4.041.
  * Dua et al. [2019b] Arpit Dua, Dominic J. Williamson, Jeongwan Haah, and Meng Cheng.  Compactifying fracton stabilizer models.  _Phys. Rev. B_ , 99:245135, Jun 2019b.  doi: 10.1103/PhysRevB.99.245135.  URL https://link.aps.org/doi/10.1103/PhysRevB.99.245135.
  * Williamson et al. [2019b] Dominic J. Williamson, Arpit Dua, and Meng Cheng.  Spurious topological entanglement entropy from subsystem symmetries.  _Phys. Rev. Lett._ , 122:140506, Apr 2019b.  doi: 10.1103/PhysRevLett.122.140506.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.140506.
  * Haah [2016b] Jeongwan Haah.  Algebraic Methods for Quantum Codes on Lattices.  _Rev. Colomb. Matemáticas_ , 50(2):299–349, 2016b.  ISSN 2357-4100.  doi: 10.15446/recolma.v50n2.62214.  URL http://arxiv.org/abs/1607.01387.
  * Haah [2021] Jeongwan Haah.  Classification of translation invariant topological pauli stabilizer codes for prime dimensional qudits on two-dimensional lattices.  _Journal of Mathematical Physics_ , 62(1):012201, 2021.  doi: 10.1063/5.0021068.  URL https://doi.org/10.1063/5.0021068.
  * Bombín [2014] Héctor Bombín.  Structure of 2D Topological Stabilizer Codes.  _Commun. Math. Phys._ , 327(2):387–432, 2014\.  ISSN 14320916.  doi: 10.1007/s00220-014-1893-4.  URL http://dx.doi.org/10.1007/s00220-014-1893-4.
  * Williamson and Devakul [2020] Dominic J. Williamson and Trithep Devakul.  Type-II fractons from coupled spin chains and layers.  2020\.  URL https://arxiv.org/abs/2007.07894.
  * Dennis et al. [2001] Eric Dennis, Alexei Kitaev, Andrew Landahl, and John Preskill.  Topological quantum memory.  _J. Math. Phys._ , 43(9):4452–4505, oct 2001\.  doi: 10.1063/1.1499754.  URL http://arxiv.org/abs/quant-ph/0110143.
  * Haah [2017] Jeongwan Haah.  Two generalizations of the cubic code model.  Talk at KITP Conference: Frontiers of Quantum Information Physics, 2017\.  URL https://online.kitp.ucsb.edu/online/qinfo_c17/haah/.
  * Dorier et al. [2005] Julien Dorier, Federico Becca, and Frédéric Mila.  Quantum compass model on the square lattice.  _Phys. Rev. B_ , 72:024448, Jul 2005.  doi: 10.1103/PhysRevB.72.024448.  URL https://link.aps.org/doi/10.1103/PhysRevB.72.024448.
  * Bacon et al. [2015] Dave Bacon, Steven T. Flammia, Aram W. Harrow, and Jonathan Shi.  Sparse quantum codes from quantum circuits.  In _Proc. Annu. ACM Symp. Theory Comput._ , volume 14-17-June-2015, pages 327–334, New York, New York, USA, jun 2015. Association for Computing Machinery.  ISBN 9781450335362.  doi: 10.1145/2746539.2746608.  URL http://dl.acm.org/citation.cfm?doid=2746539.2746608.
  * Kitaev Alexei [2006] Kitaev Alexei.  Anyons in an exactly solved model and beyond.  _Ann. Phys. (N. Y)._ , 321(1):2–111, 2006.  ISSN 0003-4916.  doi: 10.1016/j.aop.2005.10.005.  URL http://www.sciencedirect.com/science/article/pii/S0003491605002381.