Bifurcating entanglement-renormalization group flows of fracton stabilizer models

We investigate the entanglement-renormalization group flows of translation-invariant topological stabilizer models in three dimensions. Fracton models are observed to bifurcate under entanglement renormalization, generically returning at least one copy of the original model. Based on this behavior we formulate the notion of bifurcated equivalence for fracton phases, generalizing foliated fracton equivalence. The notion of quotient superselection sectors is also generalized accordingly. We calculate bifurcating entanglement-renormalization group flows for a wide range of examples and, based on those results, propose conjectures regarding the classification of translation-invariant topological stabilizer models in three dimensions.

The renormalization group (RG) is a ubiquitous concept throughout theoretical physics.Intuitively, realspace renormalization [1] involves the rescaling of a system such that the short-distance correlations are integrated out.At the scale-invariant fixed points of an RG transformation, the system has either zero or infinite correlation length, corresponding to gapped or critical phases respectively.As such, RG serves as an important tool for the analysis of long-wavelength physics in a given theory.
Since rescaling transformations increase the local Hilbert space dimension, conventional RG approaches have relied on truncating the local Hilbert space in a manner that does not affect long-range correlations.To further ensure that the long-range entanglement structure is preserved, which is particularly crucial for the classification and characterization of quantum matter, a more careful approach must be taken.A real-space renormalization procedure that exclusively removes shortrange entanglement consists of applying local unitary circuits and projecting out degrees of freedom that have been completely disentangled into trivial states only [2].Along with the coarse-graining of degrees of freedom to rescale the system, such a procedure is referred to as Entanglement Renormalization (ER) [3].Stable fixed points under ER are identified as representatives of quantum phases of matter.For many exactly solvable models with zero correlation length, ER can be implemented directly at the level of the Hamiltonian [4,5].
In this paper we apply ER to study the long-range entanglement structures of topological stabilizer models.In two dimensions, every translation-invariant topological stabilizer model is equivalent to copies of the 2D toric code [6][7][8] and hence their classification is complete.The 2D toric code is a fixed point under ER.Hence, under ER, any 2D translation-invariant topological stabilizer model flows towards copies of the 2D toric code.In contrast, translation-invariant topological stabilizer models in three dimensions exhibit a rich variety of quantum phases due to the existence of fracton topological order .While no systematic classification theorem or procedure has yet been established for these models, they can be organized into broad classes [30] based on the properties of their excitations and compactifications [31].
In Ref. 2 it was found that Haah's cubic code, the canonical example of a type-II fracton topological order with no string operators, bifurcates under ER.More precisely, under ER the cubic code splits into two decoupled models on a coarse-grained lattice: the original cubic code and cubic code B. The possibility of bifurcating ER was previously envisioned in Refs.[32][33][34] with the goal of describing critical states that violate the area law in two or more dimensions.Our goal is to extend the classification of topological phases in terms of ER fixed points to models that may bifurcate under ER.This approach faces an immediate conceptual hurdle: bifurcating models are not scale-invariant.In fact, under the conventional notion of quantum phase [35], a bifurcating model defined on a lattice of spacing a and the same model defined on a lattice of spacing 2a belong to different phases of matter.Therefore one needs to revisit the definition of thermodynamic quantum phases under such circumstances.
Fracton models are examples of bifurcating models.Under ER they may self-bifurcate or bifurcate into distinct models.Since self-bifurcating models give rise to an arbitrarily large number of copies of themselves under repeated ER, it is reasonable to consider them as free resources in the infrared limit.One can then formulate a generalized notion of fixed point where such free resources, i.e. the self-bifurcating models, are quotiented out.This generalizes the usual disentangling and projection steps in conventional ER where trivial product state degrees of freedom, which are the simplest example of a self-bifurcating state, serve as the free resource.The classification of bifurcating models via ER is then divided into two steps: first the classification of self-bifurcating fixed points and then the classification of quotient fixed points.The classification of self-bifurcating fixed points was previously studied from a resource oriented renormalization group perspective in Ref. 36  A similar point of view on fixed points of bifurcating ER has already played a key role in the understanding of foliated fracton models [38][39][40][41][42].The notion of foliated equivalence, local unitary equivalence up to adding stacks of 2D toric code, is rooted in the ER flow of foliated fracton models which produce stacks of 2D toric codes when they bifurcate.A stack of 2D toric codes is the simplest nontrivial self-bifurcating topological state in 3D.It is easy to see that it self-bifurcates under coarse-graining in the direction orthogonal to the toric code planes.The X-cube model is the canonical example with nontrivial foliated fracton order: under ER that coarse-grains by a factor of two along any axis it bifurcates, returning a copy of itself and a stack of 2D toric codes orthogonal to the axis.This reveals the X-cube model's foliation structure and that it is a foliated fixed point, as it is foliated equivalent to itself after coarse-graining.
In this work, we propose to generalize the notion of foliated fracton equivalence to bifurcated equivalence which allows any self-bifurcating state as a resource, thus providing an equivalence relation that is relevant for all fracton models.We study the ER of a large range of fracton models, including 17 of Haah's cubic codes [43] and all of Yoshida's first-order fractal spin liquids [44], and find that they are all bifurcating fixed points.That is, under ER they either self-bifurcate, producing several identical copies, or they bifurcate into a copy of themselves along with some distinct models, referred to as B models.We demonstrate that the B models are self-bifurcating for all the examples considered.Our ER results are presented in table I.Moreover, we find that the form of the B models is constrained by the mobilities of the topological quasi-particles in the original models.For instance, the presence of a planon in the original model causes a stack of 2D toric codes to appear amongst the B models.Such constraints on the B models motivate us to put forward several conjectures concerning the structure of 3D topological stabilizer models that have implications for their classification.
The paper is laid out as follows: in Sec.I B, we introduce ER for gapped quantum phases of matter, discuss conventional and bifurcating ER fixed points, and define bifurcated equivalence and quotient superselection sectors.In Sec.II, we review Haah's polynomial framework.In Sec.III, we explicitly identify bifurcating behavior in the ERG flows of a large range of models including 17 cubic codes [43] and all first-order fractal spin liquids [44].In Sec.IV, we discuss possible ERG flows for the distinct classes of topological order [30].In Sec.V, we study quotient superselection sectors in several examples.In the appendix, we provide numerical results for the number of encoded qubits for all models discussed, complimented by derivations of analytical expressions for a select few.We also present the ER for the X-cube model explicitly.The explicit ER process for all other models as listed is shown in the MATHEMATICA file SMERG.nbprovided as supplementary material.
In this section we introduce the notions of quantum phase of matter and entanglement renormalization.We discuss how the definition of phase equivalence is informed by ERG fixed points.We then introduce the more general notion of bifurcated equivalence based upon bifurcating ERG fixed points.We also discuss how the notion of superselection sector is generalized to quotient superselection sector for bifurcating ERG fixed point models.

A. Gapped quantum phases of matter
Throughout the paper we consider lattice Hamiltonians in 3D with short-range interactions.We only consider translation-invariant Hamiltonians defined on cubic lattices.
A model is in a zero temperature gapped quantum phase if, in the energy spectrum, there is a finite gap between a nearly degenerate ground state subspace and the first excited state.Technically the energy gap must be uniformly lower bounded by a positive constant for a sequence of increasing system sizes approaching the thermodynamic limit [45].The energy splitting for the nearly degenerate ground state subspace must also vanish super-polynomially in the thermodynamic limit.The thermodynamic limit is approached as the ratio of system size L to the lattice constant, or the short-distance cutoff length scale, a goes to infinity L/a → ∞.Two Hamiltonians are in the same phase if they are connected by a path of uniformly gapped Hamiltonians, possibly after adding in additional spins governed by the trivial paramagnetic Hamiltonian whose ground state is a product state.The addition of trivial degrees of freedom serves to stabilize the equivalence relation and allows for the comparison of models with a different number of spins per unit cell.In particular, models with the same L but different a can be meaningfully compared.Conversely, a quantum phase transition between gapped phases necessarily involves the gap closing.
The physical characteristics of zero temperature gapped phases can be studied via their ground states.Such ground states are in the same phase if and only if they are related by a quasi-adiabatic evolution [46], possibly after stabilization by adding spins in a trivial product state.Constant-depth local unitary circuits are often used as a proxy for phase equivalence [35], although strictly speaking they provide only a sufficient condition [47].To compare bulk phase equivalence more generally one should consider stabilized (approximate) locality-preserving unitary maps, or quantum cellular automata.For dispersionless commuting projector Hamiltonians it appears sufficient to consider stabilized exact locality-preserving unitary maps, up to a change in the choice of local Hamiltonian terms that preserves the ground space and gap but may shift higher energy levels.
Throughout this work, topological orders are stable gapped quantum phases defined by a topological degeneracy on torus and characterized by the existence of nontrivial quasiparticles.This includes: • Topological quantum liquid phases, whose ground state degeneracies on the torus have a constant upper bound.It is widely believed that these phases can be described by topological quantum field theories (TQFT) at low energy.We will thus refer to topological orders belonging to class as TQFT phases.
• Fracton phases whose ground state degeneracies do not have a constant upper bound.The low-energy behaviors of these phases are not described by any conventional TQFTs.

B. Entanglement renormalization transformations
We now describe ER transformations following the pioneering work [2].Given a gapped Hamiltonian, an ER transformation consists of the following steps: 1. Coarse-grain the system by enlarging the unit cell by a factor c > 1.
2. Apply local unitary transformations to the Hamiltonian to remove short-range entanglement.
3. Project out local degrees of freedom that are completely disentangled into a trivial state.
The result is a new Hamiltonian in the same phase of matter, defined on the coarse-grained lattice.
In order to maintain a well-defined notion of phase throughout a renormalization procedure, we fix the ratio L/a to be infinite by taking an infinite system size L while successively coarse-graining a finite lattice constant a by a factor c > 1.For a system of finite size, coarsegraining amounts to reducing the number of unit cells in the system.In the thermodynamic limit this remains true, as the ratio L1/a1 L2/a2 of the number of unit cells before and after coarse graining is equal to c even though both numerator and denominator diverge.
Conventionally it is expected that ERG flows within a gapped phase carry models towards a representative fixed point model that is invariant under the ERG flow.We depict such a situation schematically in Fig. 1, which shows Hamiltonians flowing towards RG fixed points in two gapped phases.The set of stable ERG fixed points then suffice to classify the gapped phases.
We remark that under ER copies of the trivial Hamiltonian are produced at an exponential rate, as the trivial Hamiltonian in D spatial dimensions splits into c D copies of itself under coarse-graining by a factor c. Hence to find any fixed points we must mod out by this selfreplicating trivial model.Furthermore, for an ERG fixed point model to lie in the same phase after a change of scale, the trivial Hamiltonian must also be modded out in the definition of phase.This is clearly a necessary condition for a definition of phase that does not rely on a choice of lattice scale, as is commonly desired.Fixed point Hamiltonians under the form of ER we consider here must satisfy for some finite-depth quantum circuit U , where c is the coarse-graining factor.The Hamiltonian H(ca) has the same local terms as the original Hamiltonian H(a), but the degrees of freedom sit on the sites of a coarse-grained lattice with spacing ca.The equivalence ≡ denotes equality up to the addition of disentangled spins governed by the trivial Hamiltonian to either side, and changing the choice of local Hamiltonian terms in a way that preserves the ground space.In particular, this implies that H(a) and H(ca) are in the same phase.This relation holds for commuting projector Hamiltonians that describe TQFT phases [48][49][50][51][52], and hence they are fixed points under ER [4,5].This allows the lattice scale to be ignored in the definition of TQFT phases.

C. Bifurcating entanglement renormalization and bifurcated equivalence
Models that are governed by conventional ERG fixedpoint Hamiltonians with topological order fall into the category of TQFT phases.All known two-dimensional fixed-point models are of this type.However, in three dimensions, fracton models with topological order have been discovered that bifurcate under ER.This means that performing one step of ER on a model produces multiple nontrivial decoupled models.That is, for some finite-depth quantum circuit U we have where c is again the coarse-graining factor and b is the number of nontrivial decoupled models, or branches.We FIG.
2. An illustration of bifurcating ER for a 1D model.We start with a translation-invariant spin chain, perform coarsegraining by grouping pairs of spins into new sites and then apply a local unitary consisting of two-spin gates that is translation-invariant on the coarse-grained lattice.This disentangles two decoupled copies of the original spin chain on the coarse-grained lattice.
remark that this decomposition into decoupled models may not be unique.An illustration of a bifurcating ER transformation is shown in Fig. 2.
A model H(a) is a bifurcating fixed point if any of the resulting models on the right hand side of Eq. ( 2) are equivalent to it, i.e.H 1 (a) = H(a) without loss of generality.In particular, all of the examples that appear in this work are bifurcating fixed points.Bifurcating fixed points can be either self-bifurcating fixed points or quotient fixed points which are defined as follows: is given by where H SB1 and H SB2 both satisfy Eq.(3).
In Fig. 3, we represent the above self-bifurcating and quotient ERG fixed point examples using a diagrammatic notation.In our bifurcating ERG fixed point diagrams the models resulting from one step of ER, corresponding to the right hand side of Eq. ( 2), are found by following all arrows leaving a model, corresponding to the left hand side of Eq. ( 2).Such a diagram represents a generalized fixed point when all arrows leaving models in the diagram return to models within the diagram.This captures conventional fixed points, bifurcating fixed points and limit cycles, which can be removed by increasing the amount of coarse-graining performed during one step of ER.Further examples of bifurcating ERG fixed point models are presented in section III, with the results summarized in table I.We remark that since Eq. ( 2) is not unique for a given Hamiltonian H(a) it is possible for a model to be a self-bifurcating fixed point under one ERG flow, and a quotient fixed point under another, see section III C for such an example.
In contrast to conventional ER fixed points, bifurcating ER fixed point Hamiltonians on different lattices are not in the same phase, i.e.H(a) is not phase equivalent to H(ca).This is evident from the presence of nontrivial models H i with i ≥ 2 on the right hand side of Eq. ( 2).For a self-bifurcating model, H SB (a), that satisfies Eq.( 3) with c = 2 the reason for this inequivalence is especially 4.Under one step of ER indicated by an arrow, a nontrivial self-bifurcating model denoted HSB with a lattice constant a becomes two copies of the model on the coarse-grained lattice HSB(2a).Since HSB(2a) ⊗2 is not in the same phase as HSB(2a), HSB(a) is not in the same phase as HSB(2a).A conventional phase boundary is indicated by a dotted line.
clear: two copies of H SB (2a) cannot be equivalent to a single copy unless H SB is in the trivial phase.This is depicted in Fig. 4.
Self-bifurcating fixed point models produce copies of themselves at an exponential rate under ERG flow.This rate is given by the branching number, which obviously satisfies b ≤ c D , but can be further shown to satisfy b < c D−1 for states that satisfy an area law [2], which are the relevant ones in the study of gapped phases.This is analogous to the exponential splitting of trivial models in the ERG flow of a conventional fixed point model.Inspired by the modding out of trivial models in the definition of gapped phase equivalence, here we introduce the notion of bifurcated equivalence where all self-bifurcating models are modded out.More specifically, we write H 1 (a ).See Fig. 10 (b) for a nontrivial example involving both foliated fracton equivalence and the more general bifurcated equivalence.
The bifurcated equivalence relation serves to essentially remove dependence on the lattice scale from the equivalence class of a quotient fixed point Hamiltonian since H(a) ∼ H(ca) in that case.We remark that models may be bifurcated equivalent, even when they are not in the same conventional gapped phase.In particular any self-bifurcating model is in the trivial bifurcated equivalence class, even though the model may be in a nontrivial conventional phase.In Fig. 4 the dotted line denotes a conventional phase boundary, whereas the whole diagram lies in the same trivial bifurcated equivalence class.It is also useful to generalize the equivalence relation ≡ accordingly by allowing for stacking with arbitrary selfbifurcating models, which we denote ∼ =.With this definition in hand the condition for a quotient fixed point resembles the conventional fixed point condition U H(a)U † ∼ = H(ca) . (5)

D. (Quotient) superselection sectors
A nontrivial superselection sector on some region R is an excitation that can be supported on R but not created by any operator within a neighborhood of R. Superselection sectors are equivalent if they are related by the application of an operator within a neighborhood of R, or equivalently fusion with a local excitation in R.Under ER excitations may split due to a change in the choice of local Hamiltonian terms allowed in the ≡ relation.This may create local excitations in the trivial Hamiltonians that are modded out by the ≡ relation.Hence for superselection sectors to be invariant under ER, local excitations must be modded out in their definition.
For the same reason, to define quotient superselection sectors (QSS) that are invariant under quotient ER we must mod out any excitations that can flow into a selfbifurcating fixed point model, as these models are modded out by the ∼ = relation.Representatives of potential QSS are then given by the fixed point excitations under quotient ER.This captures the notion of QSS for foliated fracton models as a special case when the self-bifurcating model is taken to be a stack of 2D topological orders [42].

II. ENTANGLEMENT RENORMALIZATION IN THE POLYNOMIAL FRAMEWORK
In this section we introduce translation-invariant stabilizer models, Clifford ER transformations and phase equivalence relations, along with their descriptions in the language of polynomial rings from commutative algebra.

A. Translation-invariant stabilizer models
Translation-invariant stabilizer Hamiltonians are specified by a choice of mutually commuting local Pauli stabilizer generators h (i) .The generators become the interaction terms in a Hamiltonian, where v are lattice vectors.In the above equation, h indicates a local generator h (i) after translation by a lattice vector v.The local generator h (i) is a tensor product of local Pauli operators acting on a set of qubits or qudits.Without loss of generality, we consider stabilizer models on a cubic lattice.

Clifford phase equivalence and entanglement renormalization
Maps between translation-invariant stabilizer Hamiltonians are given by locality-preserving Clifford operations, which map local Pauli operators to local Pauli operators.These Clifford operations include local Clifford circuits which are generated by CNOT, Phase and Hadamard gates, and nontrivial Clifford cellular automata, which are required to disentangle certain invertible phases [47].In addition, locality-preserving automorphisms of the lattice, such as the redefinition of coordinates via modular transformations and coarse graining, are also included.

Phase equivalences of translation-invariant stabilizer
Hamiltonians are given by locality-preserving Clifford operations up to stacking with trivial models.
For our Clifford ER transformations, we restrict to local Clifford circuits, coarse-graining, and discarding trivial models, as the modular transformations and other nontrivial locality-preserving operations can be moved to a single final step when comparing models.A change in the choice of local stabilizer generators that preserves the stabilizer group is also allowed when comparing two models, as in the ≡ relation above.As CNOT, Phase and Hadamard gates generate the Clifford group, they are sufficient to implement ER of Pauli stabilizer codes.
In 2D all translation-invariant topological stabilizer models were classified and shown to be equivalent, under locality-preserving Clifford operations, to copies of the 2D toric code.This implies that all translationinvariant topological stabilizer models in 2D flow to ER fixed points.Conversely in 3D examples of translationinvariant topological stabilizer models that bifurcate under ER are known, and the classification problem remains completely open, due to the existence of fracton models.Our goal is to study the bifurcating ERG flows of known fracton stabilizer models to gain clues about the 3D classification problem.
The examples considered in this work are all in CSS form [53,54], which should be preserved under ER, hence we have found it sufficient to consider Clifford circuits that consist of CNOT gates alone.In particular, the unitaries used in the ER of our examples are given by Clifford circuits U = U 1 U 2 ...U N where N is finite and each layer of gates U i is a translation-invariant tensor product of CNOT gates.

B. The polynomial framework
Translation-invariant stabilizer Hamiltonians can be conveniently expressed in terms of polynomials.The use of a polynomial description in a similar context dates back to work on classical cyclic codes [55][56][57].The polynomial approach for quantum codes on a lattice was primarily developed by Haah.Interested readers are directed to Ref. 58 for further details.We proceed by introducing several definitions from the polynomial language that are necessary and sufficient to describe ER.These definitions are demonstrated via examples.

The stabilizer map
For a stabilizer model on a cubic lattice, the stabilizer generators supported on a cubic unit cell can be expressed in terms of the position labels of the vertices on the cube as shown below The canonical example of a type-II model is Haah's code, or cubic code 1, which has the following stabilizer generators In the X-stabilizer generator, the sites on which the first qubit is acted upon by the Pauli X operator are at positions 1, xy, xz and yz of the unit cell.We take this set of positions (1, xy, xz, yz) and write a polynomial corresponding to this set 1 + xy + xz + yz.Similarly, the polynomial corresponding to the action of the Pauli X operator on the second qubit on vertices in the unit cell is 1 + x + y + z.The polynomials corresponding to the action of the Z operator on the first and second qubits on each site involved in the Z-stabilizer generator are given by xy + xz + yz + xyz and x + y + z + xyz.These polynomials refer to exponents of Pauli operators, and hence the addition of polynomials corresponds to the multiplication of operators.One can consider this polynomial representation to be a map from the position labels to the set of Pauli operators.This map is called the stabilizer map which can be written as a 2q × t matrix where q is the number of qubits on each vertex in the unit cell and t is the number of stabilizer generators per unit cell in the translation-invariant Hamiltonian.In such a matrix, each column represents a term in the Hamiltonian and all translates of these terms can be generated by acting on the column by multiplication with monomials of translation variables.We illustrate the polynomial representation for a non-CSS model using the example of Wen's plaquette model in Fig. 5.
In this paper, we focus on CSS models for which the stabilizer map takes the form where σ X (Z) is the map for the X (Z)-sector.Each column, labeled C i , specifies a stabilizer generator.The rows of the σ X (Z) block are labeled R X i (R Z i ), where the row index i is the index of the qubit in the unit cell.For example, the stabilizer map for cubic code 1 is We remark that all stabilizers of the model can be generated by the action of the stabilizer map on a column of translation variables.The translation action can be expressed in terms of the position variables x, y, z.For example, the X-stabilizer on the unit cell at position x relative to the origin is denoted Using this, or any other translation of the X-stabilizer generator in Eq (10) lead to the same Hamiltonian.Hence, multiplying columns by monomials results in an equivalent stabilizer map.For example, we could divide the second column of the cubic code stabilizer map by xyz, corresponding to a unit translation in the negative direction along each axis.The resulting equivalent stabilizer map can be written as where we have introduced the inverse variables x, y and z which satisfy xx = 1, yy = 1 and zz = 1.In the terminology of commutative algebra, introducing negative powers for translation variables involves going from a polynomial ring to a Laurent polynomial ring.

Entanglement renormalization in the polynomial language
As discussed above, in the Clifford ER procedure after coarse-graining certain operations are allowed.These operations consist of acting on the stabilizer generators with Clifford gates, changing the choice of generators for the same stabilizer group and shifting the lattice sites.
In the polynomial language, these operations are represented by matrices acting on the stabilizer map, from the left for Clifford gates and qubit shifts, and from the right for the shifting and redefinition of stabilizer generators.These operations were described in Ref. 2. As our focus is on CSS models, we restrict our discussion to Clifford circuits made up of CNOT gates.In which case the action of ER operations on the stabilizer map are as follows: • Row operations -Elementary row operations: An elementary row operation on a stabilizer map with rows R X(Z) i is specified by two row indices a = b and a monomial f and acts as follows in the X-sector, This operation corresponds to a translationinvariant implementation of CNOT gates between the target qubits specified by a and f with the control qubits specified by b.The corresponding action in the Z-sector is given by -Row multiplication by a monomial: Multiplying any of the rows in the stabilizer map by a monomial corresponds to shifting those qubits in some direction.For the polynomial entries α ab , in a row specified by constant a, the transformation α ab → x i y j z k α ab is allowed, for any finite integers i, j, k.
• Column operations -Elementary column operations: An elementary column operation on a stabilizer map with columns C i is specified by two column indices a = b and a monomial f and acts as follows: This changes the choice of stabilizer generators by replacing those corresponding to C a by products of themselves with translates of the generator corresponding to C b .Such a change of choice of generators results in a phase equivalent Hamiltonian with the same ground space and shifted excitation energy levels.
-Column multiplication by a monomial: This corresponds to changing the choice of a stabilizer generator, translating it by a monomial.This has no effect on the Hamiltonian described by the stabilizer map, as it already includes all translations of the generators.For example, we used such a transformation above to modify the polynomial representation for the Z-stabilizer term of the cubic code to go from Eq. ( 10) to Eq. ( 12).

Modular transformations
The Clifford ER process only involves equivalences between models generated by coarse-graining, local Clifford gates, shifting the lattice sites and changing the choice of column generators.More generally modular transformations, which are locality-preserving automorphisms of the cubic lattice including shear transformations, preserve the bulk properties of a model.Hence, when deciding phase equivalence, modular transformations need to be taken into account.Formally, they correspond to the redefinition of translation variables (x, y, z) to (f 1 (x, y, z), f 2 (x, y, z), f 3 (x, y, z)), where f i (x, y, z) for i = 1, 2, 3 are monomials in the translation variables that induce a bijection of the lattice sites.
Modular transformations are not used during the ER process.However, they are used when checking the equivalence of models that result from ER.We have also used them to find equivalences between some of the cubic codes.For example, cubic code 5 is related to cubic code 9 (and cubic code 15 is related to cubic code 16) via a modular transformation and hence they are in the same phase [30].

The excitation map
The excitations created by the action of a Pauli operator can be found by applying the excitation map to the polynomial representation of that Pauli operator.The map is expressed in terms of the stabilizer map σ via is the symplectic matrix and 1 denotes the q × q identity matrix.The rows in the excitation map correspond to excitations of different stabilizer generators.The operators that do not excite translations of a particular generator correspond to columns that give 0 when acted upon by the corresponding row of .For example, one can see from λ that no X stabilizer generators are excited by the action of a Pauli X operator, as expected.More generally, the kernel of contains the operators that commute with all local Hamiltonian terms.In particular, the condition that the Hamiltonian terms themselves commute can be recast as σ = 0, i.e.Im σ ⊆ ker .This containment is saturated, Im σ = ker , for infinite boundary conditions if and only if the Hamiltonian described by σ is topologically ordered [58,59].Since ker consists of operators of finite extent in the polynomial formalism, Im σ = ker implies that any operator that commutes with the Hamiltonian is in the span of the stabilizer group and so can be written as a sum of products of generators.This implies that any local operator must act as the identity on the degenerate ground space of the Hamiltonian with periodic boundary conditions, up to a proportionality constant that may be 0. In Fig. 6, we depict how the excitation map gives the positions of excited stabilizers due to the action of a local operator Z on a qubit at position xy in the Wen-plaquette model.Another phase equivalent example is the Z 2 toric code whose stabilizer map is given by and whose excitation map is given by Considering the action of a local Z operator on the first qubit at position x.The action of this local operator is represented by the following column and the action of the excitation map on this operator gives This implies that the Z stabilizers at positions x and xy are excited due to the action of ô on Z 2 toric code.Similarly as for the polynomial description of Pauli operators, addition of excitation polynomials corresponds to fusion of excitations.
The image of the excitation map im( ) contains topologically trivial configurations of excitations.For cubic code stabilizer maps, which take the following form The excitation map is given by For any CSS model, the X and the Z generator excitation sectors are decoupled, and in the case of the cubic codes they are related by a spatial inversion transformation.
For cubic codes, the polynomials in the image of the excitation map for Z generators belong to the ideal a generated by f (x, y, z) and g(x, y, z) i.e. p(x, y, z)f +q(x, y, z)g where p(x, y, z) and q(x, y, z) are arbitrary polynomials with Z 2 coefficients.We refer to this ideal, f, g , as the stabilizer ideal [58].

Coarse-graining the stabilizer and excitation maps
Coarse-graining by a factor of 2 enlarges the unit cell by the same factor.Hence, after coarse-graining, the original translation variables x 2 , y 2 and z 2 are transformed to x , y and z on the new lattice.Suppose coarse-graining by a factor of 2 in the x-direction is performed, i.e. x 2 → x .The coarse-grained unit cell has double the number of qubits and stabilizer generators.This coarse-graining transformation is implemented by the transformation of the original translation variables where x = x 2 is a new translation variable.For example, this coarse-graining sends the stabilizer map of the 2D toric code from Eq. ( 17) to It is shown below that after the application of local CNOT gates, column operations and the removal of qubits in the trivial state, the original stabilizer map is recovered.
The coarse-grained excitation map is defined in terms of the coarse-grained stabilizer map via = σ † λ.

Coarse-graining factor and trivial charge configurations
For the qubit-stabilizer models studied in this paper, we consider coarse-graining by factors of 2, i.e. c = 2.It was shown in Ref. 2 and 59 for cubic code 1 that the set of trivial charge configurations referred to as the annihilator of the charge module [59], denoted A, shows self-reproducing behavior under coarse-graining by a factor of 2. The annihilator of the cubic code 1, A, is given by the ideal 1 + x + y + z, 1 + xy + yz + xz .
a An ideal I of a polynomial ring R contains elements r I such that r I r ∈ I for all r ∈ R and all r I ∈ I.
Here, for example, 1 + x + y + z specifies a trivial charge configuration with the stabilizers excited at positions 1, x, y and z under the action of a local operator.After coarse-graining, the annihilator A c is given by 1 + x + y + z , 1 + x y + y z + x z in terms of the coarse-grained variables x = x 2 , y = y 2 , z = z 2 and hence has the same form as A. This suggests that cubic code 1 renormalizes into a model similar to itself.In fact, the annihilators of the charge modules for the two codes that are extracted after ER of cubic code 1 i.e. itself and cubic code 1B, are exactly the same.We find this self-reproducing behavior for all the cubic codes and the corresponding B models; the form of the annihilator after coarse-graining by a factor of 2, A c retains the original form as A just like in the case of cubic code 1.
Conversely, under coarse-graining by a factor of 3, the annihilator does not retain the original form.For any selfbifurcating fixed point model, it obviously follows that the models extracted after ER retain the same annihilator.For the bifurcating quotient fixed point models, which split into a copy of themselves and some B models under ER, having the same form of annihilator after coarse-graining implies that the annihilators of the B models contain the original annihilator A.

III. EXAMPLES OF BIFURCATING ENTANGLEMENT RENORMALIZATION
In this section, we find bifurcating ERG flows for explicit examples of fracton models.Some of the fracton models treated are found to be self-bifurcating fixed points while others are quotient bifurcating fixed points.Moreover, we find that some models may admit several qualitatively different ERG flows.We give such an example that is either a self-bifurcating fixed point or a quotient bifurcating fixed point, depending on the ER transformation chosen.
The models we consider include several examples with known ERG flows: the 3D toric code (3DTC), a stack of 2D toric codes along an axis î (STC î), the X-cube model (X-cube), and Haah's cubic codes 1 (CC 1 ) and 1B (CCB 1 ).Beyond these known examples we also find the ERG flows of the 16 other CSS cubic codes [43,59] 2-17 (CC 2-17 ) and the B codes thus produced (CCB 11-17 ), as well as all first-order fractal spin liquids [44] (FSL), a simple example of which is the Sierpinski fractal spin liquid [12] (SFSL).In the next section we put the ERG flows found for these examples into context by organizing them according to the type of 3D topological order each model displays [30].
We have followed a simple heuristic to find ER transformations for the example models listed in table I.We outline this process in appendix B, and apply it to the X-cube model as a demonstrative example.

A. Self-bifurcating fixed points
The coarse-graining step of the ER procedure can involve one, two or all three lattice directions.Selfbifurcating models can be divided into three categories according to the minimal number of directions one must coarse-grain for the model to bifurcate.Sorting the selfbifurcating models in this way is convenient when it comes to organizing them according to the type of topological order they exhibit, as discussed in next section.
An important necessary condition for self-bifurcation with c = b = 2 is that the number of encoded qubits under periodic boundary conditions doubles when the number of sites L/a along each axis is doubled.This is because a self-bifurcating model with c = b = 2 on a lattice with 2L/a sites along each axis splits into two copies of itself when coarse grained by a factor of 2, each of which lives on a decoupled lattice with L/a sites.
We discuss examples of stabilizer models that we have found to self-bifurcate under ER below.Examples of self-  bifurcating models and their ERG flows are depicted in Fig. 7 using our diagrammatic notation for bifurcating ERG fixed points.

A stack of 2D toric codes
The stabilizer map for a stack of 2D toric codes along the ẑ direction, parallel to the xy plane, appears identical to that in Eq. ( 17).After coarse-graining in the x direction the stabilizer map is identical to that in Eq. (24).Applying a row operation to the coarse-grained stabilizer matrix from Eq. ( 24) results in After performing additional row operations CNOT(3, 4, x ), CNOT(1, 2, x ), CNOT(1, 3, 1 + y) and column operations Col(1, 2, 1) and Col(4, 3, 1), the model becomes which is nothing but a stack of 2D toric codes and two qubits per site in a product state.Hence, the stack of toric codes along the ẑ axis is a fixed point under ER in x.
The stack of 2D toric codes stabilizer map in Eq. ( 17) has x ↔ y symmetry up to relabeling of the qubits.Hence, it is also a fixed point under ER in y.Furthermore, for the stack of 2D toric codes along the ẑ direction, we can trivially coarse-grain the stabilizer map along ẑ by taking as z does not enter the stabilizer map.This simply results in two decoupled copies of the stack of 2D toric codes Hamiltonian.Hence, a stack of 2D toric codes is self-bifurcating under ER.This provides an example of a model that self-bifurcates after ER in only one direction and is a fixed point under ER along either of the orthogonal directions.

Yoshida's fractal spin liquids
We now move on to show that a far more interesting class of examples, the first-order fractal spin liquids of Yoshida [44], are all self-bifurcating under ER.The gen-eral form of the stabilizer map for these models is where f and g are polynomials in the single translation variable x.Such a model is type-II if and only if f and g are not algebraically related [44].This class of models also contains fractal type-I lineon models for f = 1, g = 0, 1, stacks of 2D toric code for f = g = 1, stacks of 2D fractal subsystem symmetry-protected models [60][61][62][63][64] for f = 0, g = 0, 1, and decoupled 1D cluster states for f = 0, g = 1.We remark that f and g can be exchanged in the above models up to a redefinition of the lattice and an on-site qubit swap.Under coarse-graining along ŷ: where y = y 2 is the translation variable on the coarsegrained lattice, the stabilizer map becomes Applying row operations: CNOT(2, 1, f ), CNOT(3, 1, 1 + gz), CNOT(3, 4, f y) and column operations Col(2, 1, f (x)y), Col(4, 3, f (x)) leads to the stabilizer map The first and third qubits are disentangled in the above stabilizer map and hence can be removed b .The resulting stabilizer map is b Up until this step, the same ER process as shown works if 1 + g(x)z is generalized to 1 + g 1 (x)z + g 2 (x)z 2 + • • • , along with a similar generalization for 1 + g(x)z.
which is also a first-order FSL, where f (x) has been replaced by f 2 (x).We now notice that due to the symmetry between f and g in a first-order FSL one can apply essentially the same ER transformation, this time coarsegraining z, so that g is replaced by g 2 .For a polynomial over F 2 , a useful property is that f 2 (x) ≡ f (x 2 ).This leaves only functions of x 2 in the coarse-grained stabilizer map, so we coarse-grain again, along x this time, sending This results in the stabilizer map splitting into two copies of the original model, and hence all qubit first-order FSL models are self-bifurcating fixed points with b = 2.This includes the stack of 2D toric codes for f = g = 1 and the trivial model for f = g = 0. Furthermore, if f = 0, 1 then f 2 = f and the model is a conventional fixed point under ER that coarse-grains along y only, while being a self-bifurcating fixed point under ER that coarse-grains both x and z, and similarly if f and g are swapped.
A particular example, the SFSL model, is obtained when f (x) = 1 + x, g(x) = 1, this model has a fractal logical operator in the xy plane and a lineon operator along ẑ [44].Our results indicate that the model is invariant under ER along the direction of the string operator and self-bifurcates under ER in the plane of the fractal operator.
The above ER transformations were essentially based on the observation that the 2D first-order fractal subsystem symmetry breaking (classical) spin models, from which the first-order FSLs are built, are self-bifurcating under ER.This can easily be seen by following the ER transformation of the FSL before Eq. ( 33) with the second qubit and generator, as well as the z coordinate, dropped from the stabilizer map.This connection generalizes straightforwardly to provide ER transformations for 3D FSLs based upon self-bifurcating 2D fractal subsystem symmetry breaking spin models that may not be first-order.We remark that higher-order FSLs need not be self-bifurcating, as cubic code 1, which is not selfbifurcating, is equivalent to a second-order FSL [44].FSL forms for this, and some other cubic codes are presented in appendix C. The explicit mapping transformations can be found in the supplementary MATHEMAT-ICA file SMERG.nb.An interesting open problem is the classification and characterization of higher-order selfbifurcating 2D fractal subsystem symmetry breaking spin models and the 3D FSLs built from them [65].
We remark that a form of real-space RG was considered for the FSLs in Ref. 44, however it does not conform to the strict definition of ER used here, where the phase of matter cannot change.Instead, degrees of freedom that were not fully disentangled were projected out, which is capable of changing the phase by projecting out an arbitrary nontrivial decoupled model.Our ER results are consistent with the RG results in Ref. 44.
For the details of the ER procedures for the other examples discussed below, we refer the reader to the MATH- EMATICA file SMERG.nb in the Supplementary Material.
3. Cubic codes 5, 6 and 9 These cubic codes are fixed points under ER along one of the lattice directions, while they self-bifurcate under ER along the orthogonal plane, see Fig 7 (b).This is consistent with these models being fractal type-I lineon models [30].These cubic codes are not fixed points under ER along any single lattice directions, and self-bifurcate only after doing ER along all three lattice directions together, see Fig 7 (c).This is consistent with them being either fractal type-I, or type-II, fracton models [30].

B.
Quotient bifurcating fixed points A quotient bifurcating fixed point is defined by its branching under ER into a copy of itself and an inequivalent self-bifurcating model, which may consist of several decoupled self-bifurcating models.Verifying the inequivalence of models can be quite subtle, as it requires one to consider arbitrary locality-preserving unitaries, including coarse graining and modular transformations.Ref. 2 introduced an approach to proving the inequivalence of models based upon the behavior of their ground space degeneracies.In particular, the following sufficient condition was proposed to confirm that a model could not be a self-bifurcating fixed point for integers α > 1, β > 0.Where k(L) = log 2 GSD(L) is the number of encoded qubits in the ground space on an L×L×L system with periodic boundary conditions.This is in contrast with the necessary condition that β = 0 in Eq. ( 35) for any self-bifurcating model.hence is a quotient fixed point, see Fig. 8.To prove this inequivalence the scaling of the number of encoded qubits in the ground space, k(2L) = 2k(L) + 2, was used.This scaling can be directly inferred from the following formula for CC 1 [2,58,59] Thus, any model equivalent to CC 1 cannot self-bifurcate.
Since CCB 1 was demonstrated to be a self-bifurcating model, CC 1 and CCB 1 cannot be in the same phase.

Cubic codes 11-17
We have found that these cubic codes are quotient fixed points under ER.CC i , for i = 11, . . ., 17, branches into a copy of CC i , a stack of 2D toric codes along one direction and a self-bifurcating model CCB i , see Fig. 9.This is consistent with these models being fractal type-I topological orders that support planons.See Appendix D of Ref. 30 for a description of these planons.For CC 13 and CC 17 the quotient fixed point behavior occurs after an initial step of coarse-graining, see Fig. 10.
C. Quotient bifurcating and self-bifurcating fixed point behavior for cubic code 14 We have found two distinct ERG flows for CC 14 : under one ER procedure CC 14 is self-bifurcating with b = 2.Under another ER procedure CC 14 appears to be a quotient fixed point as it bifurcates into a copy of itself, a stack of 2D toric codes and another self-bifurcating model CCB 14 , see Fig. 10.Due to the absence of any correction factor in Eq. ( 35) for CC 14 , as it is a self-bifurcating fixed point, the standard method to argue for inequivalence to CCB 14 and a stack of 2D toric codes does not apply.On the other hand, we have been unable to find a direct equivalence between CC 14 and CCB 14 plus a stack of 2D toric codes.It is clear from the ERG flows that two copies of CC 14 is phase equivalent to CC 14 , CCB 14 and a stack of toric codes.Thus in the case that the ERG flows are inequivalent, this would provide an interesting example of a phase equivalence that is catalyzed by the addition of a copy of CC 14 .In any case, all of the aforementioned models are in the same trivial bifurcated equivalence class.This example raises the interesting question of whether all ERG flows for which a model is a bifurcating fixed point are equivalent.We remark that this is trivially true for conventional fixed points but requires a careful definition of equivalence for more general bifurcating ERG flows.

IV. ENTANGLEMENT RENORMALIZATION IN DIFFERENT TYPES OF 3D TOPOLOGICAL PHASES
In the recent fracton literature [30] 3D topological stabilizer models have been coarsely organized into four qualitatively distinct classes: TQFT, foliated and fractal type-I, and type-II topological orders.In this section we describe the salient features of these different classes, which are determined by the mobilities of their topological quasiparticles, and discuss the influence they have on possible ERG flows.

A. TQFT topological order
TQFT topological order is characterized by a constant topological ground space degeneracy as the system size increases and deformable logical operators.In two dimensions, there is essentially only one type of translation-invariant topological stabilizer model, the 2D toric code which is described by a TQFT at low energy.This is due to a structure theorem [47,66] stating that any translation-invariant topological stabilizer model is equivalent under locality-preserving unitary to copies of the toric code and some disentangled trivial qubits.
Hence all TQFT stabilizer models in 2D are fixed points under ER, equivalent to a number of copies of 2D toric code.We expect similar behavior in 3D, that all TQFT stabilizer models are fixed points under ER equivalent to a number of copies of 3D toric code and hence are described by a TQFT at low energies, although this remains to be shown.
In 3D the logical operator pairs of the toric code are composed of deformable string and membrane operators.More general 3D translation-invariant topological stabilizer models that have a constant ground space degeneracy as system size increases, and so are TQFT topological orders, have been shown to resemble this behavior as their logical operators also come in string-membrane pairs [67].We conjecture that a stronger structure theorem holds for such models in three dimensions: i.e. a TQFT stabilizer model in 3D is locality-preserving unitary equivalent to copies of the 3D toric code (possibly with fermionic point particle [51,68]), and disentangled trivial qubits.

B. Type-I fracton topological order
Type-I fracton topological order is characterized by a sub-extensive ground space degeneracy and rigid string operators, which correspond to excitations with subdimensional mobilities.This type of topological order is divided into two broad categories which we treat separately below.

Foliated type-I topological order
Foliated topological order is characterized via a foliation structure [39,40,69].Such models can be grown by adding layers of a 2D topological order, such as the 2D toric code, according to the foliation structure.The most studied example of this type is the X-cube model which can be grown by adding layers of 2D toric code as shown in Fig. 12.More formally, two Hamiltonians are foliated equivalent [41] if they are connected by adiabatic evolution and addition of layers of two dimensional gapped Hamiltonians.When translation invariance is enforced, foliated equivalence is the same as bifurcated equivalence with respect to stacks of 2D topological orders.Foliated stabilizer models can have topological quasiparticles with a hierarchy of mobilities such as fractons, lineons and planons.Given their foliation structure in terms of 2D toric codes, they must always support planons.In fact, the X-cube model has planons in all the three lattice directions due to its foliation structure that allows layering with 2D toric code in the three lattice directions.Due to the underlying foliation structure of this class of models, the ground space degeneracy scales exponentially with the size of the system.
The canonical self-bifurcating example within this class of models is the stack of 2D toric codes.This serves as the B model for X-cube, which is a quotient fixed point that bifurcates into a copy of itself and a stack of 2D toric codes along each axis.This is a consequence of the fact that X-cube supports planons in all the three lattice planes.The ER process for X-cube in the polynomial language is presented in appendix B.
Since all 2D topological stabilizer model are essentially copies of 2D toric code, we expect that stacks of 2D toric code are the only self-bifurcating stabilizer models in the foliated class.Furthermore, we expect that all nontrivial foliated stabilizer models flow to quotient fixed points with stacks of 2D toric code serving as the B models.An interesting open question is whether all such quotient fixed point foliated stabilizer models are foliated equivalent to a suitably generalized X-cube model which is allowed to take on different foliation structures [40][41][42]70].

Fractal type-I topological order
Fractal type-I topological order captures type-I models that support fractal logical operators and symmetries, due to which no foliation structure in terms of 2D Hamiltonians is known or likely possible.In fact, these models need not support planons at all.Due to the underlying fractal symmetry of this class of models, the ground space degeneracy fluctuates with the system size.The simplest example of this type of model is the Sierpinski fractal spin liquid model [12,44] which supports two dimensional fractal operators in the xy plane along with rigid string operators in the ẑ direction as all particles in the model are lineons.This model was mentioned above as a special case of Yoshida's first-order FSLs which were shown to be self-bifurcating.As expected, this model is invariant under ER along the ẑ direction alone but selfbifurcates under ER in the xy plane or for all the three directions at once.
There are also fractal type-I models amongst the cubic codes, some of which support fractons along with lineons and planons.We discuss the ER of two such examples below.The first is CC 6 which is self-bifurcating and the other one is CC 11 which is a quotient bifurcating fixed point: CC 6 is a self-bifurcating fixed point under ER that coarse-grains all three lattice directions together.One can perform a modular transformation on the stabilizer map of CC 6 to bring the lineon string operator along a lattice direction.The transformed model is invariant under ER that coarse grains in the direction of the string operator, while it self-bifurcates under ER that coarse grains the directions orthogonal to it.The number of encoded qubits of CC 6 is given by k(L) = 2 l+1 deg z gcd((z + 1) L + 1, (z 2 + z + 1) L + 1) .

This obeys
which implies that the number of encoded qubits always doubles upon doubling the system size.This is indeed consistent with the self-bifurcating behavior.
CC 11 is an important fractal type-I example as it supports topological quasiparticles of sub-dimensional mobilities 0, 1 and 2. Performing ER that coarse-grains all the three directions together causes CC 11 to bifurcate into a copy of itself, a lineon model CCB 11 , and a stack of 2D toric codes along x.We conjecture that this bifurcating behavior is directly related to the mobilities of topological quasiparticles in CC 11 .The presence of planons in the yz planes leads to the extraction of a stack of toric codes parallel to yz planes.Similarly, a subset of lineons from the original model CC 11 flow into the lineon model CCB 11 under ER.In fact, for each of the cubic codes 11-17, except cubic code 14, there is an ER procedure that leads to a self-bifurcating lineon B model.For CC 14 , we find a self-bifurcating fracton B model CCB 14 which supports a composite lineon as mentioned in Table III of the appendix.For CCB 11 , the scaling of the number of encoded qubits with the system size, derived in the appendix, is found to be consistent with the self-bifurcating behavior.
On the other hand, we have found that the number of encoded qubits for CC 11 is given by Using this result, we notice that there exists an infinite family of system sizes L for which the scaling of number of encoded qubits has a correction factor: for 3|L, we have The argument from Ref. 2 can then be used to show that the correction factor of 4 cannot be eliminated for any choice of coarse-graining of the original model.Hence, CC 11 cannot self-bifurcate and is therefore inequivalent to CCB 11 .Numerical results for the number of encoded qubits in the ground space of CC 12-17 indicate that such an inconsistency with self-bifurcating behavior may hold for all codes except CC 14 for which the numerics are consistent with k(2L) = 2k(L).In fact, as mentioned before, we show explicitly that there are two ERG flows for CC 14 , one in which it self-bifurcates and the other in which it splits into a copy of itself, a B model and a stack of 2D toric codes.The ER of CC 14 is shown in Fig. 10 (c).

C. Type-II topological order
Type-II topological order is defined by the absence of any logical string operators and characterized by a subextensive ground space degeneracy that fluctuates with the system size.For the models in this class, none of the elementary topological excitations or their nontrivial composites are mobile.In the previous section we discussed the canonical example from this class, cubic code 1, which is known to be a quotient fixed point.CC 1 supports a fractal operator that moves excitations apart in three dimensions.We show in the MATHEMATICA file SMERG.nb that the other type-II cubic code models c : CC 2-4 , CC 7 , CC 8 , CC 10 and a type-II first-order quantum FSL [44] are self-bifurcating fixed points.

V. QUOTIENT SUPERSELECTION SECTORS
In this section we explore the flow of excitations under ER.We find the set of quotient superselection sectors for several quotient fixed point models by calculating their fixed point excitations under quotient ERG flow.
To find fixed point excitations under a quotient ERG flow it is important to understand how excitations split during each step of ER.The local unitaries performed as part of the ER process alter the form of stabilizer generators but do not cause any splitting.It is the redefinition of the stabilizer generators allowed by the equivalence relation ≡ during ER that determines how an excitation splits.This is illustrated in Fig. 13.In the polynomial description of ER, this redefinition is captured by column operations on the stabilizer map.A representation ExS of the excitation splitting map can be found by taking the transpose of the matrix representation of the composition of all column operations involved in an ER transformation.For CSS models ExS is block-diagonal and can be decomposed into separate maps ExS X and ExS Z for the X and Z sectors, respectively.
Quotient superselection sectors are defined to be equivalence classes of superselection sectors modulo any sectors that can flow into a self-bifurcating model.For stabilizer models the QSS form an abelian group under fusion.Nontrivial QSS only arise in nontrivial quotient fixed point models.While self-bifurcating fixed point models may support highly nontrivial superselection sectors, these sectors all collapse into the trivial QSS.A nontrivial QSS must maintain support on a quotient fixed point model along its ERG flow.Conversely, any excitation that flows fully into self-bifurcating B models lies in the trivial QSS.Hence, representatives of potentially nontrivial QSS are given by fixed point excitations under quotient ER which mods out any B models.
c CC 2-4 were rigorously proven to be type-II in Ref. [43] while CC 7 , CC 8 and CC 10 were found to be type-II according to the methods in Ref. [30].
According to the above ordering, the position of each stabilizer generator in the coarse-grained unit cell can be mapped to unit vectors e i .This implies, for example, that the generator originally at position x becomes the 5 th generator in the coarse-grained unit cell, corresponding to the unit vector e 5 = 0 0 0 0 1 0 0 0 , and similarly the generator at position xyz becomes the last, corresponding to e 8 = 0 0 0 0 0 0 0 1 .As mentioned above, for all examples studied we have found that a copy of the original model appears after one step of ER.To find the QSS by looking for fixed point excitations under quotient ERG flow using the ExS map, we need to consider only those rows that correspond to the copy of the original model on the coarse-grained lattice.For a CSS model with one type of X stabilizer generator, we need to consider only one row of the ExS X map, denoted ExSR X .This map, ExSR X , specifies how X type of excitations flow under quotient ER.The quotienting out of excitations that flow into self-bifurcating B models is achieved by ignoring the corresponding rows in ExS X , which specify how excitations flow to those B models.The fixed points of ExSR X capture the potentially nontrivial QSS.
To find the fixed points of ExSR X we take the infrared limit, corresponding to many steps of ER, after which an arbitrary excitation will be supported on the coarse-grained unit cell.This allows us to restrict our attention to columns that contain only Z 2 entries with no polynomial variables.We then utilize the mapping from polynomials to unit vectors described below Eq. ( 39) to translate the monomials in ExSR X into columns, resulting in a square matrix with Z 2 entries, which we also denote ExSR X .With this replacement, the flow of excitations within the unit cell of a quotient fixed point model under quotient ER has been specified.The fixed points of ExSR X , after the replacement, provide representatives for the QSS.Several examples are worked out below.

X-cube model
Under ER, the X-cube model is mapped to itself, stacks of 2D toric code along each axis, and decoupled trivial qubits.In the Quotient Superselection Sectors section of the supplementary MATHEMATICA file SMERG.nb,column 8 contains the X-stabilizer term of X-cube model while Columns 11 and 19 contain the Z stabilizer terms.The map ExS can be found in terms of the column operations used during ER, as described above.The row corresponding to the X-stabilizer term of X-cube in the map ExSX , i.e.ExSRX , is given by x + y + x y x + y + x y 1 1 1 1 1 1 , where x , y and z are again translation variables on the coarse-grained lattice.Using the fact that the charge config-FIG.14.The flow of X-sector excitations in cubic code 1 under the ER transformation depicted in Fig. 8.The coarsegrained unit cell of CC1 is shown as a cube containing circles that represent X stabilizers, which are filled in black to represent excitations.Under ER, the X excitation at position z in the coarse-grained unit cell, denoted e2, is mapped to excitations of the X stabilizer of the coarse grained copy of CC1 represented by a white box, both X stabilizers of CCB1 represented by the purple box and a local excitation in the trivial sector represented by the gray box.uration given by 1 + x + y + x y is trivial in the X-cube model, we rewrite this row as To capture the flow of excitations from the coarse-grained unit cell to the new copy of X-cube produced by ER, we write this row in terms of unit vectors as described above in Eq. ( 39).The 1 entries are replaced by column vectors e1 e1 e1 e1 e1 e1 e1 e1 e1 , where e1 = 1 0 0 0 0 0 0 0 .Here, e1 specifies that all X stabilizer excitations within the unit cell flow to the X stabilizer excitation at position 1 of the new copy of X-cube on the coarse-grained lattice under quotient ER.Accordingly, the fixed point of this matrix is simply e1, corresponding to the Z2 fracton excitation from the X-sector of the X-cube model.
On the other hand, the Z-sector splits into 16 different stabilizer terms during ER.Hence there are 16 columns in the map ExSZ .The rows corresponding to the two Z-stabilizer terms in the copy of X-cube on the coarse-grained lattice are Using the triviality of charge configurations 1 + x and 1 + z in the X-cube model, we rewrite these rows as 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 .

Cubic code 1
After ER, cubic code 1 bifurcates into a copy of itself, cubic code B and some disentangled trivial qubits.The matrix ExSRX is given by Using the triviality of the charge configuration 1 + x + y + z in cubic code 1, we rewrite this row as which specifies the flow of excitations into the CC1 model on the coarse-grained lattice.Converting the position labels into unit vectors, we can write this as a matrix e1 e1 e1 e5 e1 e3 e2 e1 . ( This means, for example, that the X-excitation at position yz in the unit cell shown flows to an excitation of the X-stabilizer at position x in the coarse-grained lattice.Considering the next-coarse-grained unit cell, this position is equivalently expressed by the unit vector e5.The fixed point of the matrix in Eq. ( 40) is e1 which corresponds to the Z2 fracton in X-sector of cubic code 1.We remark that the Z-sector behaves similarly due to the relation between sectors in the cubic codes and hence the QSS group is given by Z2 ⊕ Z2.
In Fig. 14, we depict how the excited X stabilizer at position z splits into different excitations in the new models after ER.One of the split excitations is supported in the copy of the original model CC1 and hence the excitation at position z on the original lattice flows to the nontrivial fixed point under quotient ER and therefore is in the nontrivial QSS.In addition, both X stabilizers of CCB1 are excited and the local stabilizer on a decoupled qubit is excited.The excitation in the coarse-grained lattice of CC1 is found at position 1 in the next-coarse-grained unit cell.Upon further quotient ER, the excitation remains at position 1 as it is a fixed point.

Cubic code 11
Cubic code 11 supports fractons, lineons and planons, where the lineons and planons are composites of fractons.Under ER, cubic code 11 splits into itself, a self-bifurcating lineon model and a stack of 2D toric codes.In this case, the ExSRX map is given by Using the triviality of charge configurations 1 + y + y 2 + z + y z + y 2 z and 1 + x + y + y z in cubic code 11, we rewrite ExSRX as Converting the position labels to unit vectors, we get find the matrix e3 e3 e1 e1 e3 e1 + e3 + e2 e1 e1 + e5 + e2 , whose only nonzero fixed point is given by e1 + e3.Since the Z-sector behaves similarly we find that the QSS group is Z2 ⊕ Z2.

VI. DISCUSSIONS AND CONCLUSIONS
In this work we have systematically studied unconventional bifurcating ERG flows of gapped stabilizer Hamiltonian models.We introduced the notions of self-bifurcating and quotient bifurcating fixed points to organize the possible behavior of bifurcating fixed point models.This inspired us to define the notion of bifurcated equivalence class, generalizing the conventional notion of gapped phase and foliated fracton equivalence.This also led to a natural definition of quotient superselection sectors.These ideas were then brought to bear on a large range of stabilizer model examples, including 17 of Haah's cubic codes and Yoshida's fractal spin liquids.All models were found to be bifurcating ERG fixed points.Furthermore, many cubic codes and all first-order fractal spin liquids were found to be self-bifurcating.These results are summarized in table I.
We found that the long-range entanglement features of stabilizer models provide insights into their structure.The mobilities of particles in each model was found to constrain the nature of the self-bifurcating models produced by ER.For models that support planon particles, we found that a stack of 2D toric codes could be extracted during ER.For example, CC11 supports planons with mobility in the lattice planes perpendicular to x and consequently a stack of 2D toric codes is extracted during ER.This leads us to conjecture [71] that the existence of a planon in a 3D translation-invariant topological stabilizer model implies that a copy of the 2D toric code can be extracted via a local unitary transformation.Extending this in a translation-invariant way leads to the observed extraction of a stack of 2D toric codes.We also conjecture that if the original model has a particle with three dimensional mobility, a 3D toric code can be extracted via a local unitary circuit.
Models in which all the elementary excitations can move along a certain direction î are observed to be invariant under ER that coarse-grains along î.This is consistent with the fact that the number of encoded qubits k does not change as the system size grows along this particular direction i.e. k(2Li) = k(Li).Conversely, it was shown in Ref. 67 that all particles must have three-dimensional mobility in any 3D topological stabilizer model that has a constant ground space degeneracy as the system size increases along all three axes, and hence such models are in TQFT-type phases.Inspired by this result, we conjecture that this principle can be extended to cover models with a ground space degeneracy that is constant as the system size grows along two axes, in which case we posit that all particles in the theory must be at least mobile in a 2D plane spanned by the two axes, and consequently the model must be a stack of planon and TQFT models.Similarly we conjecture for models with a ground space degeneracy that is constant as the system size grows along one axis, that all particles must be at least mobile along that axis, and consequently such a model is a stack of lineon, planon and TQFT models.Combining this with conjectures posed previously in Ref. 30: that all TQFT stabilizer models are equivalent to copies of the 3D toric code (possibly with fermionic point particle) and planon stabilizer models are equivalent to a stack of 2D toric codes, and the fact that lineon models can be compactified [31] along the lineon direction to produce 2D subsystem symmetry breaking models, points the way towards general classification results for 3D translationinvariant topological stabilizer models, inspired by the ERG flows we have found in examples.We leave the proofs of these conjectures to a future work.
For bifurcating models, when there is an additive correction to the exponential scaling of the number of encoded qubits k as the system size increases by a factor c, i.e. k(cL) = αk(L) + β, the model cannot self-bifurcate under coarse graining by c.For models that do self-bifurcate, the additive correction vanishes i.e. k(cL) = αk(L).We have calculated the number of encoded qubits for many examples over a range of system sizes, see table IV in the appendix, and confirmed that they are consistent with the ERG flows found in table I.We remark that a correction to the linear scaling of k(L) was the first indication that X-cube is a nontrivial foliated model [14].It would be interesting to search for invariant quantities that characterize nontrivial bifurcated models, such as the correction β or similarly inspired corrections to entanglement entropy, generalizing ideas that arose in the study of foliated fracton models [41].
Our examples have focused on coarse-graining by a factor of two, which appeared natural for qubit systems.In particular each quotient fixed point example is bifurcated equivalent to itself on all lattices with spacings that are related by a multiple of two.It would be interesting to extend this to other primes, in an attempt to remove the lattice scale entirely from the bifurcated equivalence class of these models.
We also note the appearance of the 2-adic norm of L (for lattice spacing a = 1) in the ratio of the number of quotient fixed point branches to all branches, the vast majority of which are generated by self-bifurcating B models, when one starts with an initial system on the L × L × L 3D torus and applies ERG until all factors of 2 have been removed from L. Similarly for L2 = 2 n L1 the ratio GSD(L1)/GSD(L2) should be given by the 2-adic norm of L2/L1.The further appearance of p-adic norms in fracton models remains an intriguing prospect.
The bifurcating ER concepts and methods employed here apply equally well to subsystem symmetry breaking models.It would be interesting if they could shed light onto the classification of 2D translation-invariant spontaneous symmetry breaking stabilizer models, which has been accomplished in 1D [59] but remains open in 2D due to the presence of complicated fractal like symmetries [15,44].This classification problem appears to be contained within the classification of 3D translation-invariant topological stabilizer models due to the existence of lineon models that can be compactified along the lineon direction to produce any known 2D subsystem symmetry breaking model.The bifurcating ER methods can also be extend directly to subsystem symmetry-protected topological phases (SSPT), by imposing a subsystem symmetry constraint on the individual gates in each ER step.This should shed light on the question: what is the appropriate definition of SSPT equivalence class?[63,72,73] It would also be interesting to extend the ERG approach applied in this work to fractonic U(1) tensor gauge theories [74,75].It has been found that upon Higgsing such theories may transition to either fractonic or conventional topological orders [76][77][78][79], so it raises the question of how the "parent" fractonic U(1) gauge theories behave under ERG.
Finally, exact ERG flows have been found for TQFT Hamiltonians beyond stabilizer models [4].It would be interesting to extend these ERG flows to bifurcating ERG flows for nonabelian fracton models [22,80].It is currently unclear if foliated equivalence extends straightforwardly to nonabelian fracton models as their ground state degeneracies behave somewhat differently to the abelian case [24].Note added : During the writing of this paper we learnt of the work by Shirley, Slagle and Chen on closely related topics [65].
The stabilizer map for the X-cube model is given by After coarse-graining in x, the X and Z sectors of the stabilizer map are respectively given by where x is now the translation variable in the coarse-grained unit cell.The goal is to decouple different models or stabilizer terms such that they are supported on non-overlapping sets of qubits.We first do this for the X-sector.In doing so, the Z-sector will also get modified.In the two columns for the X-sector, the 4 th and 6 th row elements in both columns are the same.Thus, applying the column operation Col(2, 1, 1) simplifies the sectors of the stabilizer map Now, since there are same polynomials in certain row elements of the first column, one can get rid of them by applying CNOT(3, 4, 1) and CNOT(5, 6, 1) such that the second column is not affected due to the zero entries in it.
Finally, applying CNOT(1, 6, 1 + z) and CNOT(1, 2, 1) decouples the two stabilizer terms in the X-sector as follows, It turns out that once the X-sector is decoupled, doing certain column operations in the Z-sector is enough to decouple the whole stabilizer map into two different models.This is also observed for other models.Applying Col(5, 6, 1), Col(3, 4, 1) and Col (6,4,1) gives Removing the disentangled qubit (the 1st qubit), we have found that the model splits into two, one is given by which is just a stack of 2D toric codes parallel to the yz plane, and the other is a coarse-grained X-cube model.
Appendix C: Cubic codes as fractal spin liquids Certain cubic codes could be mapped to the following fractal spin liquid (FSL) [44] form of the stabilizer map, which represents a first-order FSL for g2(x) = 0 and a second-order FSL otherwise.Here xi ≡ x −1 i where xi denotes x, y or z.In table III, we write down the polynomials that appear in the FSL forms of various cubic codes.The explicit mapping from the cubic codes to FSLs is contained in the supplementary MATHEMATICA file SMERG.nb.
a CC 9 is equivalent to CC 5

Appendix D: Numerical results on number of encoded qubits
In table IV, we list the number of encoded qubits w.r.t.system size L in the cubic codes and the Sierpinski fractal spin liquid (SFSL) on an L × L × L lattice.The results are calculated in the supplementary MATHEMATICA file encodedqubits.nb.Results for some of these models were recovered analytically in the next section using the polynomial framework.TABLE IV.Number of encoded qubits for different stabilizer models on a system of size L × L × L where L is the number of stabilizers in each lattice direction.CCi's refer to the cubic codes, CCB11 to the cubic code 11B and SFSL to the Sierpinski Fractal spin liquid.For models that self-bifurcate under coarse-graining by a factor of 2, the number of encoded qubits also doubles as the system size is increased by a factor of 2. For CCB11, one needs a minimum cubic system size of L = 3 since the stabilizer generators are supported on 2 × 2 × 3 unit cells.The results in the table are calculated in the supplementary MATHEMATICA file encodedqubits.nbChecking conditions of Lemma 1, we have the isomorphism (F[x, y, z]/Ia)m ∼ = F[x, y, z]/Ja as rings and hence also as F-modules which are of course simply vector spaces over F. Thus we need to calculate dim F (F[x, y, z]/Ja).Note that if l = 0, then (F[x, y, z]/Ia)m ∼ = F[x, y, z]/Ja ∼ = F is the residue field at m and hence we simply get dim F (F[x, y, z]/Ja) = 1.We assume l ∈ Z>0 henceforth.
Most of the time we can eliminate one of the three variables, say z, and use the substitutions x → x+1 and y → y +1 to show F[x, y, z]/Ja ∼ = F[x, y]/J z a , as rings and hence as vector spaces over F, where J z a = (h, x 2 l + a, y 2 l + b), a = a0 + 1, b = b0 + 1 and h is the polynomial derived from f and g by eliminating the variable z.This further simplifies the problem to the calculation of dim F (F[x, y]/J z a ).In the following sections, we calculate this for various cubic codes by first finding appropriate Gröbner bases with the help of the computer algebra system SageMath.SageMath cannot calculate Gröbner bases for a general l ∈ Z>0 but we test various values of l to accurately predict the Gröbner bases in terms of l and then we prove that it is indeed so.Then we use the Gröbner bases to calculate ka,m.To calculate this, we use [82, Proposition 2.1.6]which we restate in Theorem 2.
Theorem 2. Let n ∈ Z>0 and R = F[x1, x2, . . ., xn] be a polynomial ring endowed with a monomial ordering.Let I ⊂ R be an ideal and G ⊂ R be a Gröbner basis for I. Then the set {m + I ∈ R/I : m ∈ R is a monomial which is reduced with respect to G} is a basis for the vector space R/I over F.
Finally, we use ka = m ka,m, where the sum is over all maximal ideals m ⊂ F[x, y, z]/Ia, which is feasible to compute explicitly whenever we can explicitly count the number of points in certain appropriate subvarieties of V (Ia).
Remark.All Gröbner bases in the rest of the document are using the lexicographic ordering of monomials with x < y < z.Polynomial divisions with respect to a Gröbner basis are also using the same ordering of monomials.
Remark.We often need the roots of the polynomial x 2 + x + 1 ∈ F[x] which are the primitive cube roots of unity.We denote by ζn ∈ F any choice of a primitive n th root of unity, for all n ∈ Z>0.

Cubic code 1
The stabilizer generators of cubic code 1 (CC1) are given by Hence, the stabilizer ideal that defines CC1 is a = (f, g) = (x + y + z + 1, xy + yz + zx + 1).We eliminate the variable z and use the substitutions x → x + 1 and y → y + 1 to obtain Recalling definitions, we know that a = x0 Proof.It is a straight forward calculation using the S-polynomials of Buchberger's algorithm to verify that G is a Gröbner basis for the ideal (G) for all the cases.It is also a straight forward calculation by polynomial divisions by elements in G to verify that J z a ⊂ (G) for all the cases.It remains to show that (G) ⊂ J z a for all the cases.This is shown below assuming the same hypotheses and using the same set G as in the lemma for the corresponding cases.
Case 1.The only nontrivial containment we need to show is yx 2 l −1 ∈ J z a , i.e., yx 2 l −1 ≡ 0 (mod J z a ).This follows if we show our claim nx 2 n + x 2 n −1 y + y 2 n ≡ 0 (mod J z a ) for all integers 0 ≤ n ≤ l by taking the n = l case because x 2 l ≡ y 2 l ≡ 0 (mod J z a ).We show this by induction.The base case n = 0 is trivial.Suppose that the claim holds for some integer 0 ≤ n − 1 ≤ l − 1, i.e., (n − 1)x 2 n−1 + x 2 n−1 −1 y + y 2 n−1 ≡ 0 (mod J z a ).Squaring this and using x 2 + xy + y 2 ≡ 0 (mod J z a ), we get where we use the fact that n 2 ≡ n (mod 2) by Fermat's little theorem or by directly checking both cases n ≡ 0 (mod 2) and n ≡ 1 (mod 2).
Case 2. The only nontrivial containment we need to show is (ζ3 + l)x + y ∈ J z a , i.e., (ζ3 + l)x + y ≡ 0 (mod J z a ).The equation nx 2 n + x 2 n −1 y + y 2 n ≡ 0 (mod J z a ) holds for all integers 0 ≤ n ≤ l in this case as well.The proof is exactly the same as above.By hypothesis, a, b = 0. Using the case n = l and multiplying by x a , we get x a (lx Now we wish to calculate ka = m ka,m where the sum is over all maximal ideals m ⊂ F[x, y, z]/Ia.For the calculation, we first need to explicitly calculate the points in the variety V (Ia).To solve h = 0, we divide by x 2 to get y x 2 + y x + 1 = 0 and hence the solutions are {(x 0 , ζ3x 0 ) ∈ A 2 F : x 0 ∈ F} ∪ {(x 0 , ζ3 2 x 0 ) ∈ A 2 F : x 0 ∈ F}.Thus we calculate that the solutions of f = g = 0 is the set of points V (a) = V1(a) ∪ V2(a) where we put everything together using the formula ka Now we simply read off the dimension using the Gröbner basis and obtain the formula We wish to calculate ka = m ka,m where the sum is over all maximal ideals m ⊂ F[x, y, z]/Ia.For the calculation, we first need to explicitly calculate the points in the variety V (Ia).By first finding solutions of h = 0, we easily calculate that the solutions of f = g = 0 is the set of points Since (z 2 + 1) L + 1 = ((z + 1) L + 1) 2 , we can in fact simplify the equation to #V (Ia) = deg gcd((z + 1) L + 1, (z 2 + z + 1) L + 1) .
We put everything together using the formula In general, it is difficult to obtain a more explicit formula.We note however that calculating this for a single value of L can already provide a formula for an infinite family of values of L, namely, for all L ∈ {2 l L ∈ Z>0 : l ∈ Z ≥0 }.It is not difficult to do this this explicitly by hand for small values of L .We do this now for L = 1 and L = 3.

Moreover, we have the explicit formulas
Remark.Theorem 4 implies that ka as a function of L, obeys the scaling relation for all r ∈ Z ≥0 .
Case 1.(b).The only nontrivial containment we need to show is yz + (b + 1)z ∈ J x a , i.e., yz + (b + 1)z ≡ 0 (mod J x a ).By a similar calculation as in the proof of Case 1.(a), we have This time, the hypothesis 3 | 2 l + 1 implies b 2 l +1 + 1 = 0. Hence (b + 1) Then the above implies ((b + 1)y − 1)P ≡ 1 (mod J x a ).Again using z(y + (b + 1))(y + b) = z(y 2 + y + 1), we multiply by (b + 1)P to get Case 2. The only nontrivial containment we need to show is y + (b + l) ∈ J x a , i.e., y + (b + l) ≡ 0 (mod J x a ).This follows if we show our claim y 2 l−n + (b + n) ≡ 0 (mod J x a ) for all integers 0 ≤ n ≤ l by taking the n = l case.We show this by induction.The base case n = 0 is trivial.Suppose that the claim holds for some integer 0 ≤ n − 1 ≤ l − 1, i.e., y 2 l−(n−1) + (b + (n − 1)) ≡ 0 (mod J x a ).Then The hypothesis c = 0 implies y 2 l−n + (b + n) ≡ 0 (mod J x a ) as desired.Case 3. The only nontrivial containment we need to show is z ∈ J x a , i.e., z ≡ 0 (mod J x a ).We have Now we simply read off the dimensions using the Gröbner bases and obtain the formula Recalling the definitions, the above gives the formula Note that this holds if l = 0 as well.Now we wish to calculate ka = m ka,m where the sum is over all maximal ideals m ⊂ F[x, y, z]/Ia.For the calculation, we first need to explicitly calculate the points in the variety V (Ia).By first finding solutions of h = 0, we easily calculate that the solutions of f = g = 0 is the set of points V (a) = V1(a) ∪ V2(a) ∪ V3(a) where Note that the intersection is Note that #(V2(Ia) \ V∩(Ia)) and #(V3(Ia) \ V∩(Ia)) are in fact equal in all cases using any extension σ ∈ Gal(F/F2) of σ ∈ Gal(F 2 2 /F2) specified by σ(ζ3) = ζ3 2 , noting that F2(ζ3) ∼ = F 2 2 .We put everything together using the formula ka = (2 l+2 − 2)#V∩(Ia) + 2 l+1 #(V1(Ia) \ V∩(Ia)) + 2 l+1 #(V2(Ia) \ V∩(Ia)) In general, it is difficult to obtain a more explicit formula for the case 3 | L .However, it is possible for a special case.Suppose 3 | L = 4 n − 1 for some n ∈ Z>0.We calculate that (z + 1) Thus we have the factorization We apply the same equation for (ζ3z + 1) 4 n −1 + 1 and use ζ3 3 = 1 to get the factorization (ζ3z + 1) It is easy to see that z(z+ζ3)(z 4 n −4 +z 4 n −7 +• • •+z 3 +1) is a greatest common divisor in Eq. (E7) whose degree is 4 n −2 = L −1.
Recalling definitions, we know that b0 = y0 2 l , c0 = z0 2 l where y0 L = z0 L = 1 such that b0 2 + b0 + 1 = 0. Notice that we have not made the substitutions y → y + 1, z → z + 1 in this example as the ideal is simple enough as it is.
Proof.It is a straight forward calculation using the S-polynomials of Buchberger's algorithm to verify that G is a Gröbner basis for the ideal (G) for all the cases.It is also a straight forward calculation by polynomial divisions by elements in G to verify that J x a ⊂ (G) for all the cases.It remains to show that (G) ⊂ J x a for all the cases.This is shown below assuming the same hypotheses and using the same set G as in the lemma for the corresponding cases.
We summarize the results in the following theorem.for all r ∈ Z ≥0 .This is consistent with the self-bifurcating behavior of CCB11.

Appendix F: Charge annihilators
In our argument for choosing the coarse-graining factor to be 2, we argued how the charge annihilator or the set of trivial charges shows self-reproducing behavior under coarse-graining by a factor of 2. To be precise, the charge annihilator given by the stabilizer ideal f, g changes to f 2 , g 2 after coarse-graining.In the MATHEMATICA file SMERG.nb,we show the calculation for all the cubic codes.In fact, for all of them the annihilator after coarse-graining is given by f 2 , g 2 where f and g are polynomials that give the stabilizer ideal as well as the charge annihilator of the original model.We now explain how we use elimination theory with Gröbner bases to compute the annihilator of the charge module after coarse-graining.The charge annihilator after coarse-graining is defined as the intersection of the original annihilator and the coarse-grained Laurent polynomial ring, i.e., f, g ∩ F2[x ±2 , y ±2 , z ±2 ].
Let R = F2[x, y, z] and R = F2[x 2 , y 2 , z 2 ] ⊂ R be a subring.Consider an ideal I = f, g ⊂ R for some f, g ∈ R. We wish to compute the annihilator Ann R (R/I) of the left R -module R/I.Since Ann R (R/I) = I ∩ R , we develop a method to compute I ∩ R .
Let R = F2[x, y, z, a, b, c] and R = F2[a, b, c] ⊂ R be a subring.Consider the ideal Ĩ = I + x 2 − a, y 2 − b, z 2 − c ⊂ R. Then elimination theory directly provides an algorithm to compute Ĩ ∩ R using a Gröbner basis according to [82,Theorem 2.3.4].We now focus on arguing that computing Ĩ ∩ R is essentially the same as computing I ∩ R .
Define the surjective homomorphism φ : R → R uniquely determined by the mappings Note that φ( R ) = R and φ( Ĩ) = I.A simple set theoretic calculation gives φ( Ĩ ∩ R ) ⊂ φ( Ĩ) ∩ φ( R ) = I ∩ R .However, the reverse containment does not hold in general and this is the nontriviality which is to be shown.

FIG. 5 .
FIG.5.The stabilizer map of Wen's plaquette model.The qubit positions are expressed in terms of the translation variables x, y, z in blue.The polynomial representation for each plaquette stabilizer generator is written in purple, where the first (second) row contains the positions acted upon by the X (Z) operator.These representations can be obtained by applying the stabilizer map σ to monomials that specify the positions of the respective stabilizer generators.

FIG. 6 .
FIG.6.Excitation map for Wen's plaquette model.The positions of stabilizer generators are written on the dual lattice.All but one of the labels for qubits and Pauli operators have been omitted for simplicity.A local operator Z acting at position xy and its polynomial column representation are depicted.The first (second) row of the column contains the position where X (Z) acts.Applying the excitation map , also depicted, to the polynomial column returns the positions of excited stabilizers 1 + xy (indicated by shaded circles).

FSLFIG. 7 .
FIG. 7.Self-bifurcating ERG flow diagrams.(a) ER of a stack of 2D toric codes along the ẑ direction, parallel to the xy plane, denoted by STC ẑ .The arrows labeled x and ŷ indicate directions in which the model is invariant under coarse-graining.Conversely, under coarse-graining along ẑ, STC ẑ self-bifurcates into two copies, which is indicated by a pair of arrows labeled ẑ.An arrow with no label indicates coarse-graining in all three lattice directions, in this case that produces the same result as coarse-graining in ẑ.(b) ER of a lineon model, denoted LM x, where the topological excitations of the model can all move along x.Cubic codes 5, 6, 7 and 9 are examples of similar lineon models.(c) ER of selfbifurcating cubic codes CCi for i in the range 2-10 or 14.(d) ER of any first-order fractal spin liquid[44].

14 FIG. 10 .
FIG. 10.Examples of bifurcating ERG flows.a) ER of CC13.The xy label on the arrows indicates the directions under which the pre-coarse-graining is done.After an initial coarsegraining step, CC13 splits into two copies of CC 13 .CC 13 is a quotient fixed point of the type shown in Fig. 9 (b).b) CC17 splits into two copies of CC 17 and two stacks of 2D toric codes.CC 17 is another quotient fixed point of the type shown in Fig. 9 (b).c) CC14 appears to be a quotient fixed point of the type shown in Fig. 9 (b) under one ERG flow, depicted with solid lines, but is a self-bifurcating fixed point under another ERG flow, depicted with dotted lines.

FIG. 13 .
FIG.13.An illustration of the flow of excitations under ER for a quotient fixed point model.On the left hand side, the coarse-grained unit cell of a 2D quotient fixed point model with one type of excitation per site is depicted.The presence of an excitation is indicated by a filled circle.Under ER, this excitation can split into multiple excitations that may be supported on a coarse grained copy of the original model (white box) and the B models (blue and grey boxes).The excitation in (a) is a fixed point under quotient ER and hence represents a QSS.On the other hand the excitation in (b) corresponds to a trivial QSS as it flows to the trivial fixed point under quotient ER.

2 l + 1 ,Lemma 2 . 1 .
b = y0 2 l + 1 where x0 L = y0 L = 1 such that a 2 + ab + b 2 = 0.The last equation is satisfied if and only if a = b = 0 or b a 2 + b a + 1 = 0.In the latter case b a is a primitive cube root of unity and we assume the choice ζ3 = b a .For the ideal J z a = (x 2 + xy + y 2 , x 2 l + a, y 2 l + b), we have the following Gröbner bases.Suppose a = b = 0. Then G = {x 2 + xy + y 2 , yx 2 l −1 , x 2 l } is a Gröbner basis.2. Suppose a 2 + ab + b 2 = 0 with a, b = 0. Then G = {(ζ3 + l)x + y, x 2 l + a} is a Gröbner basis.
arXiv:1909.12304v1 [cond-mat.str-el]26 Sep 2019 FIG.1.Conventional RG flows within gapped phases.Hi indicate different models that flow towards fixed points models HF P i .Here, H1-4 are in one phase while H5-9 are in another.These phases are represented by regions in parameter space, separated by a dotted phase transition line.

1 ) ∼ H B H 2 (a 2 ) if H 1 (a 1 ) stacked with some number of copies of H B lies in the same con- ventional phase as H 2 (a 2 ) stacked with some number of copies of H B . This is bifurcated equivalence with respect to a self-bifurcating fixed point model
H B .When H B is a trivial Hamiltonian we recover the conventional phase equivalence, when H B is a stack of 2D topological orders we recover foliated fracton equivalence.Similarly we write H 1 (a 1 ) H B 1 ∼ H B 2 H 2 (a 2 ) if H 1 (a 1 ) stacked with copies of H B1 is equivalent to H 2 (a 2 ) stacked with copies of H B2 .More generally we denote bifurcated equivalence by H 1 (a 1 ) ∼ H 2 (a 2 ) whenever there exist self-bifurcating models H B1 , H B2 such that H 1 (a 1 ) H B 1 ∼ H B 2 H 2 (a 2 The ER results shown require an initial coarse-graining step for CC 13 and CC 17 .b This QSS was calculated for CC 11 in Sec.V 3. The QSS of CC 12-17 can be calculated similarly.
a c CC 14 shows both self-bifurcating and bifurcating ERG behavior

TABLE III .
Polynomials in the fractal spin liquid form (C1) of cubic codes.