Bifurcating entanglement-renormalization group flows of fracton stabilizer models
Arpit Dua Pratyush Sarkar Dominic J. Williamson Meng Cheng






Department of Physics, Yale University, New Haven, CT 06520-8120, USA
Yale Quantum Institute, Yale University, New Haven, CT 06520, USA


Department of Mathematics, Yale University, New Haven, CT 06511, USA 


Current Address: Stanford Institute for Theoretical Physics, Stanford University, Stanford, CA 94305, USA
Department of Physics, Yale University, New Haven, CT 06520-8120, USA 



Department of Physics, Yale University, New Haven, CT 06520-8120, USA
Yale Quantum Institute, Yale University, New Haven, CT 06520, USA


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.  

  
Bifurcating entanglement-renormalization group flows of fracton stabilizer models
    Meng Cheng
    December 30, 2023
=================================================================================



The renormalization group (RG) is a ubiquitous concept throughout theoretical physics. 
Intuitively, real-space renormalization <cit.> 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 short-range entanglement consists of applying local unitary circuits and projecting out degrees of freedom that have been completely disentangled into trivial states only <cit.>. 
Along with the coarse-graining of degrees of freedom to rescale the system, such a procedure is referred to as Entanglement Renormalization (ER) <cit.>. 
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 <cit.>. 

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 <cit.> 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 <cit.>. 
While no systematic classification theorem or procedure has yet been established for these models, they can be organized into broad classes <cit.> based on the properties of their excitations and compactifications <cit.>. 

In haah2014bifurcation 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. Evenbly2014real,Evenbly2014class,Evenbly2014scaling 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 <cit.>, 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. swingle_greevy1,swingle_greevy2. 

A similar point of view on fixed points of bifurcating ER has already played a key role in the understanding of foliated fracton models <cit.>. 
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 <cit.> and all of Yoshida's first-order fractal spin liquids <cit.>, 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 <ref>. 
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. <ref>, 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. <ref>, we review Haah's polynomial framework. 
In Sec. <ref>, we explicitly identify bifurcating behavior in the ERG flows of a large range of models including 17 cubic codes <cit.> and all first-order fractal spin liquids <cit.>. 
In Sec. <ref>, we discuss possible ERG flows for the distinct classes of topological order <cit.>.
In Sec. <ref>, 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.nb provided as supplementary material. 



§ ENTANGLEMENT RENORMALIZATION AND PHASE EQUIVALENCE

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.  




 §.§ 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 <cit.>. 
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 <cit.>, 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 <cit.>, although strictly speaking they provide only a sufficient condition <cit.>. 
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. 




 §.§ Entanglement renormalization transformations


We now describe ER transformations following the pioneering work <cit.>. Given a gapped Hamiltonian, an ER transformation consists of the following steps:

    
  * Coarse-grain the system by enlarging the unit cell by a factor c>1.
    
  * Apply local unitary transformations to the Hamiltonian to remove short-range entanglement.
    
  * 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, coarse-graining amounts to reducing the number of unit cells in the system. 
In the thermodynamic limit this remains true, as the ratio L_1/a_1/L_2/a_2 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. <ref>, 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 self-replicating 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 

    UH(a)U^†≡ H(ca)    ,
 
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 <cit.>, and hence they are fixed points under ER <cit.>. This allows the lattice scale to be ignored in the definition of TQFT phases. 



 §.§ Bifurcating entanglement renormalization and bifurcated equivalence

Models that are governed by conventional ERG fixed-point 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

    UH(a)U^†≡ H_1(c a)+ H_2(c a)+...+H_b(c a)   ,

where c is again the coarse-graining factor and b is the number of nontrivial decoupled models, or branches. We remark that this decomposition into decoupled models may not be unique. 
An illustration of a bifurcating ER transformation is shown in Fig. <ref>. 

















































A model H(a) is a bifurcating fixed point if any of the resulting models on the right hand side of Eq. (<ref>) 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:


    
  * A model H(a) is a self-bifurcating fixed point with branching number b if all of the resulting models  on the right hand side of Eq. (<ref>) are equivalent to it, i.e. H_i(a)=H(a) for i=1,…,b. For example, a self-bifurcating fixed point with branching number b=2 satisfies
    
    UH_SB(a)U^†≡ H_SB(c a)+ H_SB(c a)   .

    
  * A bifurcating fixed point model H(a) is a quotient fixed point if the models H_i(a) for i≥ 2 are self-bifurcating fixed points that are not equivalent to H_1(a). More specifically we may refer to such an H(a) as a quotient fixed point with respect to the self-bifurcating Hamiltonian ∑_i≥ 2H_i(a). An example of a quotient fixed point model, with respect to two decoupled b=2 self-bifurcating fixed points, is given by 
    
    UH_B(a)U^†≡ H_B(c a)+ H_SB1(c a)+ H_SB2(c a)   ,

where H_SB1 and H_SB2 both satisfy Eq.(<ref>).
  
























In Fig. <ref>, 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. (<ref>), are found by following all arrows leaving a model, corresponding to the left hand side of Eq. (<ref>). 
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 <ref>, with the results summarized in table <ref>. We remark that since Eq. (<ref>) 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 <ref> 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(c a). This is evident from the presence of nontrivial models H_i with i≥ 2 on the right hand side of Eq. (<ref>). 
For a self-bifurcating model, H_SB(a), that satisfies Eq.(<ref>) with c=2 the reason for this inequivalence is especially 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. <ref>. 

















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 <cit.>, 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_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 conventional 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_B_1 is equivalent to H_2(a_2) stacked with copies of H_B_2. 
More generally we denote bifurcated equivalence by H_1(a_1) ∼ H_2(a_2) whenever there exist self-bifurcating models H_B_1,H_B_2 such that H_1(a_1) H_B_1∼_H_B_2 H_2(a_2). 
See  Fig. <ref> (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. <ref> 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 self-bifurcating models, which we denote ≅. 
With this definition in hand the condition for a quotient fixed point resembles the conventional fixed point condition 

    UH(a)U^†≅ H(ca)   .






 §.§ (Quotient) superselection sectors


A nontrivial superselection sector on some region ℛ is an excitation that can be supported on ℛ but not created by any operator within a neighborhood of ℛ. Superselection sectors are equivalent if they are related by the application of an operator within a neighborhood of ℛ, or equivalently fusion with a local excitation in ℛ. 
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 self-bifurcating 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 <cit.>. 



§ 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. 







 §.§ 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, 

    H = ∑_v⃗ ( -  h^(i)_v⃗)    ,

where v⃗ are lattice vectors. In the above equation, h^(i)_v⃗ 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 <cit.>.  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 translation-invariant topological stabilizer models in 2D flow to ER fixed points. Conversely in 3D examples of translation-invariant 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 <cit.>, 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.
 


 §.§ 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 <cit.>. 
The polynomial approach for quantum codes on a lattice was primarily developed by Haah. Interested readers are directed to Ref. haah2013commuting 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

    [ xzxyzxxyy1yzz ]  .



The canonical example of a type-II model is Haah's code, or cubic code 1, which has the following stabilizer generators 

    [ XIIIIXXIIXXXXIIX  ZIZZIZZIIZIIZIIZ ]  .

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. <ref>. 


In this paper, we focus on CSS models for which the stabilizer map takes the form  

    σ = ([ σ_X   0;   0 σ_Z ])
          ,

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 

    σ = ([   1+xy+xz+yz            0;      1+x+y+z            0;            0 xy+xz+yz+xyz;            0    x+y+z+xyz ])
         .

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 

    ([ (1+xy+xz+yz)x;    (1+x+y+z)x;             0;             0 ])
          .

Using this, or any other translation of the X-stabilizer generator in Eq (<ref>) 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

    σ = ([ 1+xy+xz+yz          0;    1+x+y+z          0;          0    1+x+y+z;          0 1+xy+xz+yz ])

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. haah2014bifurcation. 
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, 
        
    —CNOT(a,b,f):R^X_a↦ R^X_a +f(x,y,z) R^X_b  .

        This operation corresponds to a translation-invariant 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 
        
    —CNOT(a,b,f): R^Z_b ↦ R^Z_b+f(x,y,z) R^Z_a   .
 
       
    
  * 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: 
        
    Col(a,b,f): C_a↦ C_a +f(x,y,z) C_b  .
 
        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. (<ref>) to Eq. (<ref>). 
    




  §.§.§ 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 <cit.>. 




  §.§.§ 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

    ϵ:= σ^†λ,
          where      λ=(
    [  0  1; -1  0 ])

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. σ⊆ ϵ. 
This containment is saturated, σ =  ϵ,  for infinite boundary conditions if and only if the Hamiltonian described by σ is topologically ordered <cit.>. Since ϵ consists of operators of finite extent in the polynomial formalism, σ =  ϵ 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. <ref>, 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 ℤ_2 toric code whose stabilizer map is given by

    σ=([ 1+y   0; 1+x   0;   0 1+x;   0 1+y ])   ,

and whose excitation map is given by 

    ϵ= ([   0   0 1+y 1+x; 1+x 1+y   0   0;     ])   .

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

    ô=([ 0; 0; x; 0 ])   ,
 
and the action of the excitation map on this operator gives

    ϵô= ([    0; x+xy ])   .

This implies that the Z stabilizers at positions x and xy are excited due to the action of ô on ℤ_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  

    σ=([ f 0; g 0; 0 g; 0 f ]).

The excitation map is given by 

    ϵ=  ([ 0 0 f g; g f 0 0;   ])  .
 
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[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.] 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 ℤ_2 coefficients. We refer to this ideal, f,g, as the stabilizer ideal <cit.>. 





  §.§.§ 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 

    x↦([   0 x^';   1   0 ])   ,
          
    y↦(
    [ y 0; 0 y ])   , 
           z↦(
    [ z 0; 0 z ])   ,

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. (<ref>) to 

    σ^'=  ([ 1+y   0        ;   0 1+y    ;   1 x^'        ;   1   1        ;           1   1;         x^'   1;         1+y   0;           0 1+y ])   .

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. haah2013,haah2014bifurcation for cubic code 1 that the set of trivial charge configurations referred to as the annihilator of the charge module <cit.>, denoted 𝐀, shows self-reproducing behavior under coarse-graining by a factor of 2. 
The annihilator of the cubic code 1, 𝐀, is given by the ideal 1+x+y+z,1+xy+yz+xz. 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 𝐀_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 𝐀. 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, 𝐀_c retains the original form as 𝐀 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 self-bifurcating 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 𝐀.







§ 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 <cit.> 2-17 (CC_2–17) and the B codes thus produced (CCB_11–17), as well as all first-order fractal spin liquids <cit.> (FSL), a simple example of which is the Sierpinski fractal spin liquid <cit.> (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 <cit.>. 

We have followed a simple heuristic to find ER transformations for the example models listed in table <ref>. We outline this process in appendix <ref>, and apply it to the  X-cube model as a demonstrative example. 




 §.§ Self-bifurcating fixed points






The coarse-graining step of the ER procedure can involve one, two or all three lattice directions. Self-bifurcating 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 self-bifurcating 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. <ref> 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. (<ref>). 
After coarse-graining in the x̂ direction the stabilizer map is identical to that in Eq. (<ref>). 
Applying a row operation 

    CNOT(2,4,1+y)  = ([    1    0    0    0               ;    0    1    0  1+y               ;    0    0    1    0                    ;    0    0    0    1               ;                        1    0    0    0;                        0    1    0    0;                        0    0    1    0;                        0 1+y̅    0    1 ])

to the coarse-grained stabilizer matrix from Eq. (<ref>) results in

    CNOT(2,4,1+y)
       ([  1+y    0          ;    0  1+y     ;    1  x^'          ;    1    1          ;              1    1;           x̅^'    1;           1+y̅    0;              0 1+y̅ ])
       =
        ([        1+y          0                      ;        1+y          0           ;          1        x^'                      ;          1          1                      ;                                1          1;                             x̅^'          1;                             1+y̅          0;                       x̅^'(1+y̅)          0 ])
          .

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

    σ^'=  ([     0     0            ;   1+y     0      ; 1+x^'     0            ;     0     1            ;                 1     0;                 0 1+x^';                 0   1+y;                 0     0 ])

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. (<ref>) 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 

    x↦(
    [ x 0; 0 x ])   ,        y↦(
    [ y 0; 0 y ])   ,

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 <cit.>, are all self-bifurcating under ER. 
The general form of the stabilizer map for these models is 

    [ 1+f(x) y        0; 1+g(x) z        0;        0  1+g(x)z;        0  1+f(x)y ],

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 <cit.>. 
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 <cit.> 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 ŷ: 

    x↦(
    [ x 0; 0 x ])   ,        
    y↦([   0 y^';   1   0 ]) 
           z↦(
    [ z 0; 0 z ])   ,

where y^'=y^2 is the translation variable on the coarse-grained lattice, the stabilizer map becomes

    [     1  fy^'            ;     f     1            ;  1+gz     0            ;     0  1+gz            ;             1+g z     0;                 0 1+g z;                 1     f;             f y^'     1 ]  .

Applying row operations: CNOT(2,1,f), CNOT(3,1,1+gz), CNOT(3,4,fy) and column operations Col(2,1,f(x)y), Col(4,3,f(x)) leads to the stabilizer map

    [           1           0                        ;           0 1+f^2(x)y^'                        ;           0           0                        ;           0     1+g(x)z                        ;                                   0           0;                                   0     1+g(x)z;                                   1           0;                                   0 1+f^2(x)y^' ].

The first and third qubits are disentangled in the above stabilizer map and hence can be removed[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.]. 
The resulting stabilizer map is 

    [ 1+f^2y^'        0;     1+gz        0;        0     1+gz;        0 1+f^2y^' ]  ,

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 coarse-graining  z, so that g is replaced by g^2. For a polynomial over 𝔽_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

    x^2 ↦[ x^'   0;   0 x^' ]  .
 
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 ẑ <cit.>.  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. (<ref>) 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 self-bifurcating, is equivalent to a second-order FSL <cit.>. 
FSL forms for this, and some other cubic codes are presented in appendix <ref>. The explicit mapping transformations can be found in the supplementary MATHEMATICA file SMERG.nb. 
An interesting open problem is the classification and characterization of higher-order self-bifurcating 2D fractal subsystem symmetry breaking spin models and the 3D FSLs built from them <cit.>. 

We remark that a form of real-space RG was considered for the FSLs in Ref. yoshida2013exotic, 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. yoshida2013exotic. 

For the details of the ER procedures for the other examples discussed below, we refer the reader to the MATHEMATICA file SMERG.nb in the Supplementary Material.  



  §.§.§ 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 <ref> (b).  
This is consistent with these models being fractal type-I lineon models <cit.>. 



  §.§.§ Cubic codes 2–4, 7, 8 and 10

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 <ref> (c).  
This is consistent with them being either fractal type-I, or type-II, fracton models <cit.>.  















 §.§  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. haah2014bifurcation 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

    k(c L)=α k(L) + β  ,

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. (<ref>)
for any self-bifurcating model. 



  §.§.§ Cubic code 1

It was shown in Ref. haah2014bifurcation that CC_1 bifurcates into a copy of itself and an inequivalent model CCB_1 under ER and hence is a quotient fixed point, see Fig. <ref>. 
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 <cit.> 

    k(L) = 2^l + 2_x(((   x + 1 )^L' + 1, (ζ_3 x + 1)^L' + 1, 
       (ζ_3^2 x + 1)^L' + 1))- 2
      .

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. <ref>. 
This is consistent with these models being fractal type-I topological orders that support planons. See Appendix D of Ref. Dua_Classification_2019 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. <ref>. 



















































 §.§ 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. <ref>. 
Due to the absence of any correction factor in Eq. (<ref>) 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. 



§ ENTANGLEMENT RENORMALIZATION IN DIFFERENT TYPES OF 3D TOPOLOGICAL PHASES
 


In the recent fracton literature <cit.> 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. 



 §.§ 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 <cit.> 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 <cit.>. 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 <cit.>), and disentangled trivial qubits. 



 §.§ 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 sub-dimensional 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 <cit.>. 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. <ref>. More formally, two Hamiltonians are foliated equivalent <cit.> 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 <ref>. 

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 <cit.>. 



  §.§.§ 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 <cit.> 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 self-bifurcates 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_z(((z + 1)^L' + 1, (z^2 + z + 1)^L' + 1)).

This obeys

    k(2^r L) = 2^r k(L)

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 <ref> 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

    2L,      3 ∤ L 
    
    2L - 4  + 2^l + 2_z(((z + 1)^L' + 1,    
    2L - 4  + 2^l + 2_z ( gcd (ζ_3 z + 1)^L' + 1)),     3 | L   .

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

    k(2L)= 2k(L)+4   .

The argument from Ref. haah2014bifurcation 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. <ref> (c). 



 §.§ Type-II topological order

Type-II topological order is defined by the absence of any logical string operators and characterized by a sub-extensive 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[CC_2–4 were rigorously  proven to be type-II in Ref. <cit.> while CC_7, CC_8 and CC_10 were found to be type-II according to the methods in Ref. <cit.>.]: CC_2–4, CC_7, CC_8, CC_10 and a type-II first-order quantum FSL <cit.> are self-bifurcating fixed points. 




§ 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. <ref>. 
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. 
































































Each row of the map ExS  corresponds to a different type of stabilizer generator in the coarse-grained model. 
The entries in a row are the coefficients for expressing a new stabilizer generator as a linear combination of the original stabilizer generators within the coarse-grained unit cell, 

    [ xzxyzxxyy1yzz ].

Here, x, y and z are the translation variables of the original lattice. For the choice of coarse-graining map described in Sec. <ref>, the original stabilizer generators within the coarse-grained unit cell are ordered as follows 

    1-z-y-yz-x-xz-xy-xyz   .

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 ℤ_2 entries with no polynomial variables. We then utilize the mapping from polynomials to unit vectors described below Eq. (<ref>) to translate the monomials in ExSR_X into columns, resulting in a square matrix with ℤ_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 ExS_X, i.e. ExSR_X, 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 configuration given by 1+x^'+y^'+x^' y^' is trivial in the X-cube model, we rewrite this row as 

    ([ 1 1 1 1 1 1 1 1 ])
      .

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. (<ref>). The 1 entries are replaced by column vectors e_1 

    ([ e_1 e_1 e_1 e_1 e_1 e_1 e_1 e_1 ])   ,

where e_1=([ 1 0 0 0 0 0 0 0 ]). 
Here, e_1 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 e_1, corresponding to the ℤ_2 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 ExS_Z. The rows corresponding to the two Z-stabilizer terms in the copy of X-cube on the coarse-grained lattice are 


    ([          1        z^'          1          1          1          z          1          1;      1+x^' z^'+x^'z^'          0          0          0          0          0          0 ][               1+z^'                   0                   0                   0                   0          x^'+x^'z^'                   0                   0;                 x^'        1+z^'+x^'z^'                   1                   1                 x^' x^'+x^'^2z^'+x^'z^'                 x^'                   1 ])
      .


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 ])
      .

Again, converting the position labels to unit vectors, we obtain the map for the flow of excitations 

    ([ e_1 e_1 e_1 e_1 e_1 e_1 e_1 e_1; e_0 e_0 e_0 e_0 e_0 e_0 e_0 e_0 ][ e_0 e_0 e_0 e_0 e_0 e_0 e_0 e_0; e_1 e_1 e_1 e_1 e_1 e_1 e_1 e_1 ])   ,

where e_0 =([ 0 0 0 0 0 0 0 0 ]).  
The fixed points are

    ([ e_1; e_0 ])       and      ([ e_0; e_1 ])

which correspond to the two generating ℤ_2 lineons in the Z-sector of the X-cube model. 
Hence we find that the QSS group of X-cube is given by ℤ_2 ⊕ℤ_2^2, which agrees with the result in Ref. shirley2018Fractional.  

















































  §.§.§ Cubic code 1



After ER, cubic code 1 bifurcates into a copy of itself, cubic code B and some disentangled trivial qubits. The matrix ExSR_X is given by

    ([           1           1           1         x^'           1   1+x^'+z^'   1+x^'+y^' x^'+y^'+z^' ]).

Using the triviality of the charge configuration 1+x^'+y^'+z^' in cubic code 1, we rewrite this row as 

    ([   1   1   1 x^'   1 y^' z^'   1 ])   ,

which specifies the flow of excitations into the CC_1 model on the coarse-grained lattice. Converting the position labels into unit vectors, we can write this as a matrix

    ([ e_1 e_1 e_1 e_5 e_1 e_3 e_2 e_1 ])   .

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 e_5. 
The fixed point of the matrix in Eq. (<ref>) is e_1 which corresponds to the ℤ_2 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 ℤ_2 ⊕ℤ_2. 

In Fig. <ref>, 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 CC_1 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 CCB_1 are excited and the local stabilizer on a decoupled qubit is excited. The excitation in the coarse-grained lattice of CC_1 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 ExSR_X map is given by

    ([            y^'            y^'              1              1            y^'      1+y^'+z^'              1 1+y^'^2+y^'^2z ])   .

Using the triviality of charge configurations  1+y^'+y^'^2+z^'+y^'z^'+y^'^2z^' and 1+x^'+y^'+y^'z^' in cubic code 11, we rewrite ExSR_X  as 

    ([       y^'       y^'         1         1       y^' 1+y^'+z^'         1 1+x^'+z^' ]).

Converting the position labels to unit vectors, we get find the matrix 

    ([         e_3         e_3         e_1         e_1         e_3 e_1+e_3+e_2         e_1 e_1+e_5+e_2 ])   ,

whose only nonzero fixed point is given by e_1+e_3. Since the Z-sector behaves similarly we find that the QSS group is ℤ_2 ⊕ℤ_2. 



§ 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 <ref>. 


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, CC_11 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 <cit.> 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(2L_i)=k(L_i). 

Conversely, it was shown in Ref. yoshida2011classification 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. Dua_Classification_2019: 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 <cit.> along the lineon direction to produce 2D subsystem symmetry breaking models,  
points the way towards general classification results for 3D translation-invariant 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(c L)=α 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(c L)=α k(L). 
We have calculated the number of encoded qubits for many examples over a range of system sizes, see table <ref> in the appendix, and confirmed that they are consistent with the ERG flows found in table <ref>. 
We remark that a correction to the linear scaling of k(L) was the first indication that X-cube is a nontrivial foliated model <cit.>. 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 <cit.>. 

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 L_2 = 2^n L_1 the ratio GSD(L_1)/GSD(L_2) should be given by the 2-adic norm of L_2/L_1. 
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 <cit.> but remains open in 2D due to the presence of complicated fractal like symmetries <cit.>. 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? <cit.> 

It would also be interesting to extend the ERG approach applied in this work to fractonic U(1) tensor gauge theories <cit.>. It has been found that upon Higgsing such theories may transition to either fractonic or conventional topological orders <cit.>, 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 <cit.>. It would be interesting to extend these ERG flows to bifurcating ERG flows for nonabelian fracton models <cit.>. 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 <cit.>. 



Note added: During the writing of this paper we learnt of the work by Shirley, Slagle and Chen on closely related topics <cit.>. 

AD and DW thank Jeongwan Haah for useful discussions. DW also thanks Yichen Hu and Abhinav Prem. PS thanks Jayameenakshi Venkatraman for introducing him to quantum error correction. This work is supported by start-up funds at Yale University (DW and MC), NSF under award number DMR-1846109 and the Alfred P. Sloan foundation (MC). 
 
apsrev_1




§ LIST OF CODES AND THEIR POLYNOMIAL REPRESENTATIONS


In Table <ref>, we list the polynomials f and g for various cubic codes, denoted CC_i where i runs from 1 to 17, and their B-models if any, denoted CCB_i. The stabilizer map for a cubic code in terms of f and g is given by 

    σ=([ f 0; g 0; 0 g; 0 f ])   ,

and the stabilizer ideal is defined as the ideal f,g, which contains the polynomials in the image of the excitation map, see Eq. (<ref>). 
The third column of the table below contains polynomial entries, each of which specifies the positions of two non-trivial clusters of charges that are created at the ends of string operators. For example, y+xz^2 implies that the small segment of string operator in CC_5 creates elementary excitations at y and xz^2. This means that the direction of each string operator can be read off the binomial factor in its description. 





§ DEMONSTRATING OUR HEURISTIC PROCEDURE FOR ENTANGLEMENT RENORMALIZATION ON THE X-CUBE MODEL


In this appendix, we outline our heuristic procedure for entanglement renormalization on CSS codes and demonstrate it, using the X-cube model as an example. 
The heuristic is quite simple, essentially we coarse grain until a monomial appears in a column. We next use that entry to set the rest of the column to 0 via CNOT gates. The monomial entry can then be used to set its row to 0 via column operations. 
After this the qubit has been disentangled into the trivial state and its row and column can be removed. 
This is repeated until there are no monomials left. At which point the model has either bifurcated or we perform more steps of coarse graining. 
There is an important subtlety to the above recipe, as we do not necessarily need a monomial to set the rest of a column and row to 0. In some cases, such as the X-cube example below, a non-monomial entry, pair of entries, or several entries, in a column can be used to set the rest of the column to zero. This more general step can be used in place of the monomial step mentioned above. 

The stabilizer map for the X-cube model is given by 

    [ (1+y)(1+z)                      ; (1+x)(1+z)           ; (1+x)(1+y)           ;                   1+x        1+x;                     0        1+y;                   1+z          0 ]  .
 
After coarse-graining in x, the X and Z sectors of the stabilizer map are respectively given by

    [ (1+y)(1+z)          0;          0 (1+y)(1+z);        1+z     x(1+z);        1+z        1+z;        1+y     x(1+y);        1+y        1+y ]  ,
              [   1   1   1   1;   x   1   x   1;   0   0 1+y   0;   0   0   0 1+y; 1+z   0   0   0;   0 1+z   0   0 ]  ,

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

    [ (1+y)(1+z) (1+y)(1+z);          0 (1+y)(1+z);        1+z (1+x)(1+z);        1+z          0;        1+y (1+x)(1+y);        1+y          0 ]  , 
              [   1   1   1   1;   x   1   x   1;   0   0 1+y   0;   0   0   0 1+y; 1+z   0   0   0;   0 1+z   0   0 ]  .

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. 

    [ (1+y)(1+z) (1+y)(1+z);          0 (1+y)(1+z);          0 (1+x)(1+z);        1+z          0;          0 (1+x)(1+y);        1+y          0 ]  ,
              [   1   1   1   1;   x   1   x   1;   0   0 1+y   0;   0   0 1+y 1+y; 1+z   0   0   0; 1+z 1+z   0   0 ]  .

Finally, applying CNOT(1,6,1+z) and CNOT(1,2,1) decouples the two stabilizer terms in the X-sector as follows, 

    [          0          0;          0 (1+y)(1+z);          0 (1+x)(1+z);        1+z          0;          0 (1+x)(1+y);        1+y          0 ]  ,
              [   1   1   1   1; 1+x   0 1+x   0;   0   0 1+y   0;   0   0 1+y 1+y; 1+z   0   0   0;   0   0 1+z 1+z ]  .

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 

    [          0          0;          0 (1+y)(1+z);          0 (1+x)(1+z);        1+z          0;          0 (1+x)(1+y);        1+y          0 ]  ,
              [   0   1   0   0; 1+x   0 1+x   0;   0   0 1+y   0;   0   0   0 1+y; 1+z   0   0   0;   0   0   0 1+z ]  .

Removing the disentangled qubit (the 1st qubit), we have found that the model splits into two, one is given by

    [ 1+z    ; 1+y    ;     1+y;     1+z ],

which is just a stack of 2D toric codes parallel to the yz plane, and the other is a coarse-grained  X-cube model.



§ CUBIC CODES AS FRACTAL SPIN LIQUIDS


Certain cubic codes could be mapped to the following fractal spin liquid (FSL) <cit.> form of the stabilizer map,

    ([            1+f(x) y                   0; 1+g_1(x)z+g_2(x)z^2                   0;                   0 1+g_1(x)z+g_2(x)z^2;                   0             1+f(x)y ])
          ,

which represents a first-order FSL for g_2(x)=0 and a second-order FSL otherwise. Here x_i≡ x_i^-1 where x_i denotes x, y or z. In table <ref>, 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.






§ NUMERICAL RESULTS ON NUMBER OF ENCODED QUBITS

In table <ref>, 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. 




§ NUMBER OF ENCODED QUBITS


In Ref. haah2013commuting, Haah used techniques from commutative algebra to derive a formula for the number of encoded qubits for the cubic code CC_1 which we recount in <ref>. We use the general strategy used by Haah to derive similar formulas for the number of encoded qubits for other cubic codes. We consider CC_6 and CC_11 as examples. We first explain the general definitions and explain the steps used by Haah in his derivation. 

Suppose 𝔞 = (f, g) ⊂𝔽_2[x, y, z] is an ideal corresponding to a code. Fix some L ∈ℤ_> 0. Imposing the periodicity conditions x^L - 1 = y^L - 1 = z^L - 1 = 0, the number of encoded qubits is given by

    k_𝔞 = 2_𝔽_2(𝔽_2[x, y, z]/I_𝔞)

where _𝔽_2(𝔽_2[x, y, z]/I_𝔞) is the dimension of 𝔽_2[x, y, z]/I_𝔞 as a vector space and I_𝔞 is the ideal defined as

    I_𝔞 = 𝔞 + (x^L - 1, y^L - 1, z^L - 1) = (f, g, x^L - 1, y^L - 1, z^L - 1).

Due to algebrogeometric reasons, it is preferable to work over an algebraically closed field. Hence we take the algebraic closure 𝔽 = 𝔽_2. By extension of scalars of the vector spaces, we have

    k_𝔞 = 2_𝔽(𝔽[x, y, z]/I_𝔞)

where I_𝔞 is now the ideal generated in 𝔽[x, y, z]. In order to calculate this, we use the structure theorem for Artinian rings from <cit.> which we restate in <ref>.



Let A be an Artinian ring. Then A = ∏_j = 1^n A_j for some n ∈ℤ_>0 where A_j is a local Artinian ring for all integers 1 ≤ j ≤ n. Moreover, the decomposition is unique up to isomorphism.



In the above theorem, if A is also a vector space over 𝔽, we have A = ⊕_j = 1^n A_j and in particular _𝔽(A) = ∑_j = 1^n _𝔽(A_j).


We apply <ref> to obtain the formula k_𝔞 = ∑_𝔪 k_𝔞, 𝔪 where

    k_𝔞, 𝔪 = 2_𝔽((𝔽[x, y, z]/I_𝔞)_𝔪)

and the sum is taken over all maximal ideals 𝔪⊂𝔽[x, y, z]/I_𝔞. By the weak form of Hilbert's Nullstellensatz, any maximal ideal 𝔪⊂𝔽[x, y, z]/I_𝔞 is of the form 𝔪 = (x - x_0, y - y_0, z - z_0)/I_𝔞 where (x_0, y_0, z_0) ∈ V(I_𝔞) ⊂𝔸_𝔽^3, i.e., (x_0, y_0, z_0) ∈𝔸_𝔽^3 satisfies

    f(x_0, y_0, z_0) = g(x_0, y_0, z_0) = x_0^L - 1 = y_0^L - 1 = z_0^L - 1 = 0.

In particular x_0, y_0, z_0 ≠ 0. Since x_0^L-1=x_0^2^l L^'-1=(x_0^L^'-1)^2^l=0, we also have x_0^L^'=1 and similarly, y_0^L^'=z_0^L^'=1. We now calculate k_𝔞, 𝔪 for all maximal ideals 𝔪⊂𝔽[x, y, z]/I_𝔞.

Fix a maximal ideal 𝔪 = (x - x_0, y - y_0, z - z_0)/I_𝔞⊂𝔽[x, y, z]/I_𝔞. Now write L = 2^l L' where 2 ∤ L' and l ∈ℤ_≥ 0. Consider the factorization

    x^L - 1    = x^L - x_0^L = (x^L')^2^l - (x_0^L')^2^l = (x^L' - x_0^L')^2^l
       = ((x - x_0)(x^L' - 1 + x^L' - 2x_0 + … + xx_0^L' - 2 + x_0^L' - 1))^2^l
       = (x - x_0)^2^l(x^L' - 1 + x^L' - 2x_0 + … + xx_0^L' - 2 + x_0^L' - 1)^2^l.

The right most factor has no further factors of x - x_0 since putting x = x_0 results in L'x_0^L' - 1≠ 0 since 2 ∤ L' and x_0 ≠ 0. We obtain similar factorization for the variables y and z. We recognize that (x - x_0)^2^l = x^2^l + x_0^2^l and similarly for the variables y and z since (𝔽) = 2, which motivates us to define the ideal

    J_𝔞 = 𝔞 + (x^2^l + a_0, y^2^l + b_0, z^2^l + c_0) = (f, g, x^2^l + a_0, y^2^l + b_0, z^2^l + c_0)

where a_0 = x_0^2^l, b_0 = y_0^2^l and c_0 = z_0^2^l. We have the canonical map 𝔽[x, y, z]/I_𝔞→𝔽[x, y, z]/J_𝔞 since I_𝔞⊂ J_𝔞. We would like to apply <cit.> which we restate in <ref>.



Let A and B be rings and S ⊂ A be a multiplicatively closed set. If g: A → B is a ring homomorphism such that


  * for all s ∈ S, the element g(s) ∈ B is a unit

  * for all a ∈(g), we have as = 0 for some s ∈ S

  * for all b ∈ B, there are a ∈ A and s ∈ S such that b = g(a)g(s)^-1

then there is a unique isomorphism h: S^-1A → B such that g = h ∘ f where f: A → S^-1A is the canonical homomorphism.


Checking conditions of <ref>, we have the isomorphism

    (𝔽[x, y, z]/I_𝔞)_𝔪≅𝔽[x, y, z]/J_𝔞

as rings and hence also as 𝔽-modules which are of course simply vector spaces over 𝔽. Thus we need to calculate _𝔽(𝔽[x, y, z]/J_𝔞). Note that if l = 0, then (𝔽[x, y, z]/I_𝔞)_𝔪≅𝔽[x, y, z]/J_𝔞≅𝔽 is the residue field at 𝔪 and hence we simply get _𝔽(𝔽[x, y, z]/J_𝔞) = 1. We assume l ∈ℤ_>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 𝔽[x, y, z]/J_𝔞≅𝔽[x, y]/J_𝔞^z, as rings and hence as vector spaces over 𝔽, where J_𝔞^z = (h, x^2^l + a, y^2^l + b), a = a_0 + 1, b = b_0 + 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 _𝔽(𝔽[x, y]/J_𝔞^z). 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 ∈ℤ_> 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 k_𝔞, 𝔪. To calculate this, we use <cit.> which we restate in <ref>.



Let n ∈ℤ_>0 and R = 𝔽[x_1, x_2, …, x_n] 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 𝔽.


Finally, we use k_𝔞 = ∑_𝔪 k_𝔞, 𝔪, where the sum is over all maximal ideals 𝔪⊂𝔽[x, y, z]/I_𝔞, which is feasible to compute explicitly whenever we can explicitly count the number of points in certain appropriate subvarieties of V(I_𝔞).

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.


We often need the roots of the polynomial x^2 + x + 1 ∈𝔽[x] which are the primitive cube roots of unity. We denote by ζ_n ∈𝔽 any choice of a primitive nth root of unity, for all n ∈ℤ_>0.




 §.§ Cubic code 1


The stabilizer generators of cubic code 1 (CC_1) are given by 

    [ XIIIIXXIIXXXXIIX  ZIZZIZZIIZIIZIIZ ]

Hence, the stabilizer ideal that defines CC_1 is 𝔞 = (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

    J_𝔞^z = (h, x^2^l + a, y^2^l + b) = (x^2 + xy + y^2, x^2^l + a, y^2^l + b).

Recalling definitions, we know that a = x_0^2^l + 1, b = y_0^2^l + 1 where x_0^L' = y_0^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 = (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.

  * 	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.




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 ⊂ (G) for all the cases. It remains to show that (G) ⊂ J_𝔞^z 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.

<ref>. The only nontrivial containment we need to show is yx^2^l - 1∈ J_𝔞^z, i.e., yx^2^l - 1≡ 0 J_𝔞^z. This follows if we show our claim nx^2^n + x^2^n - 1y + y^2^n≡ 0 J_𝔞^z for all integers 0 ≤ n ≤ l by taking the n = l case because x^2^l≡ y^2^l≡ 0 J_𝔞^z. 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 - 1y + y^2^n - 1≡ 0 J_𝔞^z. Squaring this and using x^2 + xy + y^2 ≡ 0 J_𝔞^z, we get

    (n - 1)^2x^2^n + x^2^n - 2y^2 + y^2^n≡ 0 J_𝔞^z
       (n^2 - 1)x^2^n + x^2^n - 2(x^2 + xy) + y^2^n≡ 0 J_𝔞^z
       nx^2^n + x^2^n - 1y + y^2^n≡ 0 J_𝔞^z

where we use the fact that n^2 ≡ n 2 by Fermat's little theorem or by directly checking both cases n ≡ 0 2 and n ≡ 1 2.

<ref>. The only nontrivial containment we need to show is (ζ_3 + l)x + y ∈ J_𝔞^z, i.e., (ζ_3 + l)x + y ≡ 0 J_𝔞^z. The equation nx^2^n + x^2^n - 1y + y^2^n≡ 0 J_𝔞^z 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^2^l + x^2^l - 1y + y^2^l) ≡ 0 J_𝔞^z
       1/a(lx^2^l + 1 + x^2^ly + bx) ≡ 0 J_𝔞^z
       1/a(lax + ay + bx) ≡ 0 J_𝔞^z
       (b/a + l)x + y ≡ 0 J_𝔞^z
       (ζ_3 + l)x + y ≡ 0 J_𝔞^z.


Now we simply read off the dimension using the Gröbner basis and obtain the formula

    _𝔽(𝔽[x, y]/J_𝔞^z) =
    
    2^l + 1 - 1,     a = b = 0 
    
    2^l,     a^2 + ab + b^2 = 0  with  a, b ≠ 0.

Recalling definitions, the above gives the formula

    k_𝔞, 𝔪 =
    
    2^l + 2 - 2,     x_0 = y_0 = 1 
    
    2^l + 1,     (x_0 + 1)^2 + (x_0 + 1)(y_0 + 1) + (y_0 + 1)^2 = 0  with  x_0, y_0 ≠ 1.

Now we wish to calculate k_𝔞 = ∑_𝔪 k_𝔞, 𝔪 where the sum is over all maximal ideals 𝔪⊂𝔽[x, y, z]/I_𝔞. For the calculation, we first need to explicitly calculate the points in the variety V(I_𝔞). 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', ζ_3 x_0') ∈ A_𝔽^2: x_0' ∈𝔽}∪{(x_0', ζ_3^2 x_0') ∈ A_𝔽^2: x_0' ∈𝔽}. Thus we calculate that the solutions of f = g = 0 is the set of points V(𝔞) = V_1(𝔞) ∪ V_2(𝔞) where

    V_1(𝔞)    = {(x_0' + 1, ζ_3 x_0' + 1, ζ_3^2 x_0' + 1) ∈ A_𝔽^3: x_0' ∈𝔽}
    
    V_2(𝔞)    = {(x_0' + 1, ζ_3^2 x_0' + 1, ζ_3 x_0' + 1) ∈ A_𝔽^3: x_0' ∈𝔽}.

Note that the intersection is V_∩(𝔞) = V_1(𝔞) ∩ V_2(𝔞) = {(1, 1, 1)} which consists of the only point satisfying x_0 = y_0 = 1. Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0, we have the corresponding set of solutions V(I_𝔞) = V_1(I_𝔞) ∪ V_2(I_𝔞) where V_j(I_𝔞) = {(x_0, y_0, z_0) ∈ V_j(𝔞): x_0^L' = y_0^L' = z_0^L' = 1} for all j ∈{1, 2}. The intersection is V_∩(I_𝔞) = {(1, 1, 1)}. Thus, #V_∩(I_𝔞) = 1 and

    #(V_1(I_𝔞) ∖ V_∩(I_𝔞)) = #(V_2(I_𝔞) ∖ V_∩(I_𝔞)) 
    
    =   (((x + 1)^L' + 1, (ζ_3 x + 1)^L' + 1, (ζ_3^2 x + 1)^L' + 1)) - 1.

we put everything together using the formula

    k_𝔞   = (2^l + 2 - 2)#V_∩(I_𝔞) + 2^l + 1#(V_1(I_𝔞) ∖ V_∩(I_𝔞)) + 2^l + 1#(V_2(I_𝔞) ∖ V_∩(I_𝔞)) 
       = 2^l + 2(((x + 1)^L' + 1, (ζ_3 x + 1)^L' + 1, (ζ_3^2 x + 1)^L' + 1)) - 2.

In general, it is difficult to obtain a more explicit formula. However, we obtain a more explicit formula for some special cases.

First suppose L' = 1. Then it is easy to see that x is a greatest common divisor in <ref>. Hence <ref> gives k_𝔞 = 2^l + 2 - 2 = 4L - 2.

Now suppose L' = 2^n + 1 for some n ∈ℤ_≥ 0. Then we have the factorization

    (x + 1)^2^n + 1 + 1    = (x + 1)(x + 1)^2^n + 1 = (x + 1)(x^2^n + 1) + 1 
       = x^2^n + 1 + x^2^n + x = x(x^2^n + x^2^n - 1 + 1).

We apply the same equation for (ζ_3 x + 1)^2^n + 1 + 1 and (ζ_3^2 x + 1)^2^n + 1 + 1 to get the factorization

    (ζ_3 x + 1)^2^n + 1 + 1    = ζ_3x((ζ_3 x)^2^n + (ζ_3 x)^2^n - 1 + 1) 
    
    (ζ_3^2 x + 1)^2^n + 1 + 1    = ζ_3^2 x((ζ_3^2 x)^2^n + (ζ_3^2 x)^2^n - 1 + 1).

Let r be a root of the second factor x^2^n + x^2^n - 1 + 1. Suppose that 2^n ≡ 1 3. In this case we have

    (ζ_3 r)^2^n + (ζ_3 r)^2^n - 1 + 1 = ζ_3 r^2^n + r^2^n - 1 + 1 = ζ_3 r^2^n + r^2^n = r^2^n(ζ_3 + 1)^2^n.

So (ζ_3 r)^2^n + (ζ_3 r)^2^n - 1 + 1 = 0 implies r = 0 which is a contradiction. Suppose that 2^n ≡ 2 3. In this case we have

    (ζ_3 r)^2^n + (ζ_3 r)^2^n - 1 + 1 = ζ_3^2 r^2^n + ζ_3 r^2^n - 1 + 1 = ζ_3^2 r^2^n + ζ_3 r^2^n + ζ_3 + 1.

So (ζ_3 r)^2^n + (ζ_3 r)^2^n - 1 + 1 = 0 implies

    ζ_3 r^2^n(ζ_3 + 1) + (ζ_3 + 1) = 0 
       r^2^n + ζ_3^2^n = 0 
       (r + ζ_3)^2^n = 0.

Thus r = ζ_3. But then

    (ζ_3^2 r)^2^n + (ζ_3^2 r)^2^n - 1 + 1 = (ζ_3^2 ·ζ_3)^2^n + (ζ_3^2 ·ζ_3)^2^n - 1 + 1 = 1 + 1 + 1 = 1 ≠ 0.

Thus, in any case, there is no common root among the second factors. So x is a greatest common divisor in the <ref>. Now <ref> gives k_𝔞 = 2^l + 2 - 2.

Now suppose L' = 4^n - 1 for some n ∈ℤ_>0. We calculate that

    (x + 1)^4^n = x^4^n - 1 = (x - 1)(x^4^n - 1 + x^4^n - 2 + … + x + 1) 
       (x + 1)^4^n - 1 + 1 = x^4^n - 1 + x^4^n - 2 + … + x^2 + x.

Thus we have the factorization

    (x + 1)^4^n - 1 + 1    = x(x^4^n - 2 + x^4^n - 3 + … + x + 1) 
       = x(x^2 + x + 1)(x^4^n - 4 + x^4^n - 7 + … + x^3 + 1) 
       = x(x + ζ_3)(x + ζ_3^2)(x^4^n - 4 + x^4^n - 7 + … + x^3 + 1).

We apply the same equation for (ζ_3 x + 1)^4^n - 1 + 1 and (ζ_3^2 x + 1)^4^n - 1 + 1 to get the factorizations

    (ζ_3 x + 1)^4^n - 1 + 1    = x(x + 1)(x + ζ_3)(x^4^n - 4 + x^4^n - 7 + … + x^3 + 1) 
    
    (ζ_3^2 x + 1)^4^n - 1 + 1    = x(x + 1)(x + ζ_3^2)(x^4^n - 4 + x^4^n - 7 + … + x^3 + 1).

It is easy to see that x(x^4^n - 4 + x^4^n - 7 + … + x^3 + 1) is a greatest common divisor in <ref> whose degree is 4^n - 3 = L' - 2. Using this in <ref> gives k_𝔞 = 4L - 2^l + 3 - 2.

Now suppose L' = 2^2n + 1 - 1 for some n ∈ℤ_≥ 0. Suppose r is a root of (x + 1)^2^2n + 1 - 1 + 1. We calculate that

    (x + 1)^2^2n + 1 = x^2^2n + 1 - 1 = (x - 1)(x^2^2n + 1 - 1 + x^2^2n + 1 - 2 + … + x + 1) 
       (x + 1)^2^2n + 1 - 1 + 1 = x^2^2n + 1 - 1 + x^2^2n + 1 - 2 + … + x^2 + x.

Thus we have the factorization

    (x + 1)^2^2n + 1 - 1 + 1 = x(x^2^2n + 1 - 2 + x^2^2n + 1 - 3 + … + x + 1).

We apply the same equation for (ζ_3 x + 1)^2^2n + 1 - 1 + 1 and (ζ_3^2 x + 1)^2^2n + 1 - 1 + 1 to get the factorizations

    (ζ_3 x + 1)^2^2n + 1 - 1 + 1    = ζ_3 x((ζ_3 x)^2^2n + 1 - 2 + (ζ_3 x)^2^2n + 1 - 3 + … + ζ_3 x + 1) 
    
    (ζ_3^2 x + 1)^2^2n + 1 - 1 + 1    = ζ_3^2 x((ζ_3^2 x)^2^2n + 1 - 2 + (ζ_3^2 x)^2^2n + 1 - 3 + … + ζ_3^2 x + 1).

Suppose r is a root of the second factor x^2^2n + 1 - 2 + x^2^2n + 1 - 3 + … + x + 1. By multiplying by r - 1, we find that r^2^2n + 1 - 1 = 1, i.e., r is a (2^2n + 1 - 1)th root of unity. Similarly, so are ζ_3 r and ζ_3^2 r. But then we have

    (ζ_3 r)^2^2n + 1 - 1 = 1  (ζ_3)^2^2n + 1 - 1 = 1

which is a contradiction because 2^2n + 1 - 1 ≡ 1 3 for all n ∈ℤ_≥ 0. Thus x is a greatest common divisor in <ref>. Now <ref> gives k_𝔞 = 2^l + 2 - 2. We summarize the results in the following theorem.



Consider CC_1 defined by the ideal 𝔞 = (x + y + z + 1, xy + yz + zx + 1). Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0 for some L = 2^l L' where 2 ∤ L' ∈ℤ_>0 and l ∈ℤ_≥ 0, the number of encoded qubits is given by

    k_𝔞 = 2^l + 2(((x + 1)^L' + 1, (ζ_3 x + 1)^L' + 1, (ζ_3^2 x + 1)^L' + 1)) - 2.

Moreover, we have the explicit formulas

    k_𝔞 =
    
    4L - 2,     L = 2^l 
    
    2^l + 2 - 2,     L = 2^l · (2^n + 1) 
    
    4L - 2^l + 3 - 2,     L = 2^l · (4^n - 1) 
    
    2^l + 2 - 2,     L = 2^l · (2^2n + 1 - 1).




<ref> implies that k_𝔞 as a function of L, obeys the scaling relation

    k_𝔞(2^r L) = 2^r k_𝔞(L) + 2(2^r - 1)

for all r ∈ℤ_≥ 0.




 §.§ Cubic code 6


The stabilizer generators of cubic code 6 (CC_6) are given by 

    [ XIIIIXIIXXXXXIIX  ZZZZIZZIIZIIZIII ]

Hence, the stabilizer ideal that defines CC_6 is 𝔞 = (f, g) = (x + y + z + 1, yz + zx + y + 1). We eliminate the variable y to obtain

    J_𝔞^y = (h, x^2^l + a, z^2^l + c) = (z^2 + x, x^2^l + a, z^2^l + c).

Recalling definitions, we know that a = x_0^2^l + 1, c = z_0^2^l + 1 where x_0^L' = z_0^L' = 1 such that c^2 + a = 0.



For the ideal J_𝔞^y = (z^2 + x, x^2^l + a, z^2^l + c), we have that

    G = {z^2 + x, x^2^l - 1 + c}

is a Gröbner basis.



We assume the same hypotheses and use the same set G as in the lemma. 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). It is also a straight forward calculation by polynomial divisions by elements in G to verify that J_𝔞^y ⊂ (G). It remains to show that (G) ⊂ J_𝔞^y.

The only nontrivial containment we need to show is x^2^l - 1 + c ∈ J_𝔞^y, i.e., x^2^l - 1 + c ≡ 0 J_𝔞^y. We calculate that

    x^2^l + a = (x^2^l - 1)^2 + c^2 = (x^2^l - 1 + c)^2.

Hence x^2^l + a ≡ 0 J_𝔞^y implies x^2^l - 1 + c ≡ 0 J_𝔞^y as desired.


Now we simply read off the dimension using the Gröbner basis and obtain the formula

    k_𝔞, 𝔪 = 2_𝔽(𝔽[x, y]/J_𝔞^y) = 2^l + 1.

We wish to calculate k_𝔞 = ∑_𝔪 k_𝔞, 𝔪 where the sum is over all maximal ideals 𝔪⊂𝔽[x, y, z]/I_𝔞. For the calculation, we first need to explicitly calculate the points in the variety V(I_𝔞). By first finding solutions of h = 0, we easily calculate that the solutions of f = g = 0 is the set of points V(𝔞) = {(z_0'^2 + 1, z_0'^2 + z_0' + 1, z_0' + 1) ∈𝔸_𝔽^3: z_0' ∈𝔽}. Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0, we have the corresponding set of solutions V(I_𝔞) = {(x_0, y_0, z_0) ∈ V(𝔞): x_0^L' = y_0^L' = z_0^L' = 1}. Thus

    #V(I_𝔞) = (((z^2 + 1)^L' + 1, (z^2 + z + 1)^L' + 1, (z + 1)^L' + 1)).

Since (z^2 + 1)^L' + 1 = ((z + 1)^L' + 1)^2, we can in fact simplify the equation to

    #V(I_𝔞) = (((z + 1)^L' + 1, (z^2 + z + 1)^L' + 1)).

We put everything together using the formula

    k_𝔞 = 2^l + 1#V(I_𝔞) = 2^l + 1(((z + 1)^L' + 1, (z^2 + z + 1)^L' + 1)).

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^lL' ∈ℤ_>0: l ∈ℤ_≥ 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.

First suppose L' = 1. Then we have the factorizations

    (z + 1) + 1    = z
    
    (z^2 + z + 1) + 1    = z^2 + z = z(z + 1).

So we calculate that

    #V(I_𝔞) = (((z + 1) + 1, (z^2 + z + 1) + 1)) = (z) = 1.

Thus k_𝔞 = 2^l + 1 = 2L.

Now suppose L' = 3. Then we have the factorizations

    (z + 1)^3 + 1 = z^3 + 3z^2 + 3z + 1 + 1 = z^3 + z^2 + z = z(z + ζ_3)(z + ζ_3^2)

and

    (z^2 + z + 1)^3 + 1    = (z^2)^3 + 3(z^2)^2(z + 1) + 3z^2(z + 1)^2 + (z + 1)^3 + 1 
       = z^6 + z^5 + z^3 + z 
       = z(z^5 + z^4 + z^2 + 1).

So ((z + 1)^3 + 1, (z^2 + z + 1)^3 + 1) = z because we check that ζ_3 and ζ_3^2 are not roots of z^5 + z^4 + z^2 + 1 by putting z = ζ_3 and z = ζ_3^2. Again #V(I_𝔞) = 1. Thus k_𝔞 = 2^l + 1 = 2/3L.

We summarize the results in the following theorem.



Consider CC_6 defined by the ideal 𝔞 = (x + y + z + 1, yz + zx + y + 1). Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0 for some L = 2^l L' where 2 ∤ L' ∈ℤ_>0 and l ∈ℤ_≥ 0, the number of encoded qubits is given by

    k_𝔞 = 2^l + 1(((z + 1)^L' + 1, (z^2 + z + 1)^L' + 1)).

Moreover, we have the explicit formulas

    k_𝔞 =
    
    2L,     L = 2^l 
    2/3L,     L = 2^l · 3.




<ref> implies that k_𝔞 as a function of L, obeys the scaling relation

    k_𝔞(2^r L) = 2^r k_𝔞(L)

for all r ∈ℤ_≥ 0.




 §.§ Cubic code 11


The stabilizer generators of cubic code 11 (CC_11) are given by 

    [ IXIIXIXXIXXXIIXI  ZIZZIIIZZIIIIZZZ ]

Hence, the stabilizer ideal that defines CC_11 is 𝔞 = (f, g) = (yz + x + y + 1, xy + x + y + z). We eliminate the variable x and use the substitutions y ↦ y + 1 and z ↦ z + 1 to obtain

    J_𝔞^x = (h, y^2^l + b, z^2^l + c) = (z(y^2 + y + 1), y^2^l + b, z^2^l + c).

Recalling definitions, we know that b = y_0^2^l + 1, c = z_0^2^l + 1 where y_0^L' = z_0^L' = 1 such that c(b^2 + b + 1) = 0. The last equation is satisfied if and only if b^2 + b + 1 = 0 or c = 0.



For the ideal J_𝔞^x = (z(y^2 + y + 1), y^2^l + b, z^2^l + c), we have the following Gröbner bases.


  * Suppose b^2 + b + 1 = 0 and c = 0. Then we have the following cases.


  * 	Suppose 3 ∤ 2^l + 1. Then

    G = {yz + bz, y^2^l + b, z^2^l}

is a Gröbner basis.

  * 	Suppose 3 | 2^l + 1. Then

    G = {yz + (b + 1)z, y^2^l + b, z^2^l}

is a Gröbner basis.


  * 	Suppose b^2 + b + 1 = 0 and c ≠ 0. Then

    G = {y + (b + l), z^2^l + c}

is a Gröbner basis.

  * 	Suppose b^2 + b + 1 ≠ 0 and c = 0. Then

    G = {z, y^2^l + b}

is a Gröbner basis.




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 ⊂ (G) for all the cases. It remains to show that (G) ⊂ J_𝔞^x 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.

<ref>. The only nontrivial containment we need to show is yz + bz ∈ J_𝔞^x, i.e., yz + bz ≡ 0 J_𝔞^x. We calculate that

    (by - 1)((by)^2^l - 1 + (by)^2^l - 2 + … + by + 1) 
    ≡   (by)^2^l - 1 J_𝔞^x
    ≡   b^2^l· y^2^l - 1 J_𝔞^x
    ≡   b^2^l + 1 + 1 J_𝔞^x.

The hypothesis b^2 + b + 1 = 0 implies b ≠ 1 with b^3 = 1. Hence the hypothesis 3 ∤ 2^l + 1 implies b^2^l + 1 + 1 ≠ 0. Let P = (by)^2^l - 1 + (by)^2^l - 2 + … + by + 1. Realizing that z(y + b)(y + (b + 1)) = z(y^2 + y + 1), we multiply by b/b^2^l + 1 + 1P to get

    z(y + b)(y + (b + 1)) ·b/b^2^l + 1 + 1P ≡ 0 J_𝔞^x
       z(y + b) ·(by - 1)P/b^2^l + 1 + 1≡ 0 J_𝔞^x
       yz + bz ≡ 0 J_𝔞^x.


<ref>. The only nontrivial containment we need to show is yz + (b + 1)z ∈ J_𝔞^x, i.e., yz + (b + 1)z ≡ 0 J_𝔞^x. By a similar calculation as in the proof of <ref>, we have

    ((b + 1)y - 1)(((b + 1)y)^2^l - 1 + ((b + 1)y)^2^l - 2 + … + (b + 1)y + 1) 
    ≡   (b + 1)^2^l + 1 + 1 J_𝔞^x.

This time, the hypothesis 3 | 2^l + 1 implies b^2^l + 1 + 1 = 0. Hence

    (b + 1)^2^l + 1 + 1    = (b + 1)(b + 1)^2^l = (b + 1)(b^2^l + 1) = b^2^l + 1 + b^2^l + b + 1 
       = b^2^l + b = b^2^l + 1 + b^2/b = b^2 + 1/b = b/b = 1.

Let P = ((b + 1)y)^2^l - 1 + ((b + 1)y)^2^l - 2 + … + (b + 1)y + 1. Then the above implies ((b + 1)y - 1)P ≡ 1 J_𝔞^x. Again using z(y + (b + 1))(y + b) = z(y^2 + y + 1), we multiply by (b + 1)P to get

    z(y + (b + 1))(y + b) · (b + 1)P ≡ 0 J_𝔞^x
       z(y + (b + 1)) · ((b + 1)y - 1)P ≡ 0 J_𝔞^x
       yz + (b + 1)z ≡ 0 J_𝔞^x.


<ref>. The only nontrivial containment we need to show is y + (b + l) ∈ J_𝔞^x, i.e., y + (b + l) ≡ 0 J_𝔞^x. This follows if we show our claim y^2^l - n + (b + n) ≡ 0 J_𝔞^x 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 J_𝔞^x. Then

    z^2^l - 2^l - n· (z(y^2 + y + 1))^2^l - n + z^2^l(y^2^l - (n - 1) + (b + (n - 1))) ≡ 0 J_𝔞^x
       z^2^l(y^2^l - n + 1 + y^2^l - n + 1) + z^2^l(y^2^l - n + 1 + (b + (n - 1))) ≡ 0 J_𝔞^x
       z^2^l(y^2^l - n + (b + n)) ≡ 0 J_𝔞^x
       c(y^2^l - n + (b + n)) ≡ 0 J_𝔞^x.

The hypothesis c ≠ 0 implies y^2^l - n + (b + n) ≡ 0 J_𝔞^x as desired.

<ref>. 
The only nontrivial containment we need to show is z ∈ J_𝔞^x, i.e., z ≡ 0 J_𝔞^x. We have

    ((y^2 + y) - 1)((y^2 + y)^2^l - 1 + (y^2 + y)^2^l - 2 + … + (y^2 + y) + 1) 
    ≡   (y^2 + y)^2^l - 1 J_𝔞^x
    ≡   (y^2^l)^2 + y^2^l + 1 J_𝔞^x
    ≡   b^2 + b + 1 J_𝔞^x.

Now b^2 + b + 1 ≠ 0 by hypothesis. Let P = (y^2 + y)^2^l - 1 + (y^2 + y)^2^l - 2 + … + (y^2 + y) + 1. Then the above implies (y^2 + y + 1) ·P/b^2 + b + 1≡ 1 J_𝔞^x. Thus

    z ≡ z(y^2 + y + 1) ·P/b^2 + b + 1≡ 0 J_𝔞^x.


Now we simply read off the dimensions using the Gröbner bases and obtain the formula

    _𝔽(𝔽[x, y]/J_𝔞^x)
       =
    
    2^l + 1 - 1,     b^2 + b + 1 = 0  and  c = 0 
    
    2^l,     b^2 + b + 1 = 0  and  c ≠ 0 
    
    2^l,     b^2 + b + 1 ≠ 0  and  c = 0.

Recalling the definitions, the above gives the formula

    k_𝔞, 𝔪   =
    
    2^l + 2 - 2,     y_0 ∈{ζ_3, ζ_3^2} and  z_0 = 1 
    
    2^l + 1,     y_0 ∈{ζ_3, ζ_3^2} and  z_0 ≠ 1 
    
    2^l + 1,     y_0 ∉{ζ_3, ζ_3^2} and  z_0 = 1.

Note that this holds if l = 0 as well. Now we wish to calculate k_𝔞 = ∑_𝔪 k_𝔞, 𝔪 where the sum is over all maximal ideals 𝔪⊂𝔽[x, y, z]/I_𝔞. For the calculation, we first need to explicitly calculate the points in the variety V(I_𝔞). By first finding solutions of h = 0, we easily calculate that the solutions of f = g = 0 is the set of points V(𝔞) = V_1(𝔞) ∪ V_2(𝔞) ∪ V_3(𝔞) where

    V_1(𝔞)    = {(1, y_0, 1) ∈𝔸_𝔽^3: y_0 ∈𝔽}
    
    V_2(𝔞)    = {(ζ_3 z_0' + 1, ζ_3, z_0' + 1) ∈𝔸_𝔽^3: z_0' ∈𝔽}
    
    V_3(𝔞)    = {(ζ_3^2 z_0' + 1, ζ_3^2, z_0' + 1) ∈𝔸_𝔽^3: z_0' ∈𝔽}.

Note that the intersection is

    V_∩(𝔞) = V_1(𝔞) ∩ V_2(𝔞) ∩ V_3(𝔞) = {(1, ζ_3, 1), (1, ζ_3^2, 1)}

which consists of the only points satisfying both y_0 ∈{ζ_3, ζ_3^2} and z_0 = 1. Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0, we have the corresponding set of solutions V(I_𝔞) = V_1(I_𝔞) ∪ V_2(I_𝔞) ∪ V_3(I_𝔞) where V_j(I_𝔞) = {(x_0, y_0, z_0) ∈ V_j(𝔞): x_0^L' = y_0^L' = z_0^L' = 1} for all j ∈{1, 2, 3}. The intersection is

    V_∩(I_𝔞) = V_1(I_𝔞) ∩ V_2(I_𝔞) ∩ V_3(I_𝔞) =
    ∅,     3 ∤ L' 
    {(1, ζ_3, 1), (1, ζ_3^2, 1)},     3 | L'.

Thus we immediately calculate

    #V_∩(I_𝔞)    =
    
    0,     3 ∤ L' 
    
    2,     3 | L'
    
    #(V_1(I_𝔞) ∖ V_∩(I_𝔞))    =
    
    L',     3 ∤ L' 
    
    L' - 2,     3 | L'
    
    #(V_2(I_𝔞) ∖ V_∩(I_𝔞))    =
    
    0,     3 ∤ L' 
    (((z + 1)^L' + 1, (ζ_3 z + 1)^L' + 1)) - 1,     3 | L'
    
    #(V_3(I_𝔞) ∖ V_∩(I_𝔞))    =
    
    0,     3 ∤ L' 
    (((z + 1)^L' + 1, (ζ_3^2 z + 1)^L' + 1)) - 1,     3 | L'
    .

Note that #(V_2(I_𝔞) ∖ V_∩(I_𝔞)) and #(V_3(I_𝔞) ∖ V_∩(I_𝔞)) are in fact equal in all cases using any extension σ̃∈(𝔽/𝔽_2) of σ∈(𝔽_2^2/𝔽_2) specified by σ(ζ_3) = ζ_3^2, noting that 𝔽_2(ζ_3) ≅𝔽_2^2. We put everything together using the formula

    k_𝔞 =   (2^l + 2 - 2)#V_∩(I_𝔞) + 2^l + 1#(V_1(I_𝔞) ∖ V_∩(I_𝔞)) + 2^l + 1#(V_2(I_𝔞) ∖ V_∩(I_𝔞)) 
       + 2^l + 1#(V_3(I_𝔞) ∖ V_∩(I_𝔞))
    
    
    =   
    2L,     3 ∤ L 
    
    2L - 4 + 2^l + 2(((z + 1)^L' + 1, (ζ_3 z + 1)^L' + 1)),     3 | L.

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 ∈ℤ_>0. We calculate that

    (z + 1)^4^n = z^4^n - 1 = (z - 1)(z^4^n - 1 + z^4^n - 2 + … + z + 1) 
       (z + 1)^4^n - 1 + 1 = z^4^n - 1 + z^4^n - 2 + … + z^2 + z.

Thus we have the factorization

    (z + 1)^4^n - 1 + 1    = z(z^4^n - 2 + z^4^n - 3 + … + z + 1) 
       = z(z^2 + z + 1)(z^4^n - 4 + z^4^n - 7 + … + z^3 + 1) 
       = z(z + ζ_3)(z + ζ_3^2)(z^4^n - 4 + z^4^n - 7 + … + z^3 + 1).

We apply the same equation for (ζ_3 z + 1)^4^n - 1 + 1 and use ζ_3^3 = 1 to get the factorization

    (ζ_3 z + 1)^4^n - 1 + 1 = z(z + 1)(z + ζ_3)(z^4^n - 4 + z^4^n - 7 + … + z^3 + 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 <ref> whose degree is 4^n - 2 = L' - 1. Using this in <ref> gives k_𝔞 = 6L - 4(2^l + 1).

We summarize the results in the following theorem.



Consider CC_11 defined by the ideal 𝔞 = (yz + x + y + 1, xy + x + y + z). Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0 for some L = 2^l L' where 2 ∤ L' ∈ℤ_>0 and l ∈ℤ_≥ 0, the number of encoded qubits is given by

    k_𝔞 =
    
    2L,     3 ∤ L 
    
    2L - 4 + 2^l + 2(((z + 1)^L' + 1, (ζ_3 z + 1)^L' + 1)),     3 | L.

Moreover, if 3 | L' = 4^n - 1 for some n ∈ℤ_>0, then we have the explicit formula

    k_𝔞 = 6L - 4(2^l + 1).




<ref> implies that k_𝔞 as a function of L, obeys the scaling relation

    k_𝔞(2^r L) =
    
    2^r k_𝔞(L),     3 ∤ L 
    
    2^r k_𝔞(L) + 4(2^r - 1),     3 | L.

for all r ∈ℤ_≥ 0.




 §.§ Cubic code 11B

The stabilizer generators of cubic code 11B (CCB_11) are given by 

    [ @!0   @-[rr]       @-[rr]        IX @-[dl]
    
    XI @-[ur]        XX @-[ll]@-[ur]@-[rr]        XX
       @.[uu]@.[dl]@.[rr]       @.[rr] @.[uu]       @-[uu]
    @-[uu]@-[rr]        IX @-[uu]@.[ur] @-[rr]       @-[uu] @-[ur]
      @!0   @-[rr]        ZI @-[rr]       @-[dl] 
    @-[ur]       @-[ll]@-[ur]       @-[ll]  
         ZZ@.[uu]@.[dl]@.[rr]        ZZ@.[rr] @.[uu]        IZ @-[uu]  
    
      ZI @-[uu]@-[rr]       @-[uu]@.[ur] @-[rr]       @-[uu] @-[ur] ]

Hence, the stabilizer ideal that defines CCB_11 is given by 𝔞 = (f, g) = (1 + y + y^2,1 + x + y + yz). We eliminate the variable x to obtain

    J_𝔞^x = (h, y^2^l + b_0, z^2^l + c_0) = (y^2 + y + 1, y^2^l + b_0, z^2^l + c_0).

Recalling definitions, we know that b_0 = y_0^2^l, c_0 = z_0^2^l where y_0^L' = z_0^L' = 1 such that b_0^2 + b_0 + 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. 



For the ideal J_𝔞^x = (y^2 + y + 1, y^2^l + b_0, z^2^l + c_0), we have the following Gröbner bases.


  * 	Suppose 3 ∤ 2^l + 1. Then

    G = {y + b_0, z^2^l + c_0}

is a Gröbner basis.

  * 	Suppose 3 | 2^l + 1. Then

    G = {y + b_0 + 1, z^2^l + c_0}

is a Gröbner basis.




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 ⊂ (G) for all the cases. It remains to show that (G) ⊂ J_𝔞^x 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.

<ref>. The only nontrivial containment we need to show is y + b_0 ∈ J_𝔞^x, i.e., y + b_0 ≡ 0 J_𝔞^x. We calculate that

    (b_0y - 1)((b_0y)^2^l - 1 + (b_0y)^2^l - 2 + … + b_0y + 1) 
    ≡   (b_0y)^2^l - 1 J_𝔞^x
    ≡   b_0^2^l· y^2^l - 1 J_𝔞^x
    ≡   b_0^2^l + 1 + 1 J_𝔞^x.

The hypothesis b_0^2 + b_0 + 1 = 0 implies b_0 ≠ 1 with b_0^3 = 1. Hence the hypothesis 3 ∤ 2^l + 1 implies b_0^2^l + 1 + 1 ≠ 0. Let P = (b_0y)^2^l - 1 + (b_0y)^2^l - 2 + … + b_0y + 1. Realizing that (y + b_0)(y + (b_0 + 1)) = (y^2 + y + 1), we multiply by b_0/b_0^2^l + 1 + 1P to get

    (y + b_0)(y + (b_0 + 1)) ·b_0/b_0^2^l + 1 + 1P ≡ 0 J_𝔞^x
       (y + b_0) ·(b_0y - 1)P/b_0^2^l + 1 + 1≡ 0 J_𝔞^x
       y + b_0 ≡ 0 J_𝔞^x.


<ref>. The only nontrivial containment we need to show is y + (b_0 + 1) ∈ J_𝔞^x, i.e., y + (b_0 + 1) ≡ 0 J_𝔞^x. By a similar calculation as in the proof of <ref>, we have

    ((b_0 + 1)y - 1)(((b_0 + 1)y)^2^l - 1 + ((b_0 + 1)y)^2^l - 2 + … + (b_0 + 1)y + 1) 
    ≡   (b_0 + 1)^2^l + 1 + 1 J_𝔞^x.

This time, the hypothesis 3 | 2^l + 1 implies b_0^2^l + 1 + 1 = 0. Hence

    (b_0 + 1)^2^l + 1 + 1    = (b_0 + 1)(b_0 + 1)^2^l = (b_0 + 1)(b_0^2^l + 1) = b_0^2^l + 1 + b_0^2^l + b_0 + 1 
       = b_0^2^l + b_0 = b_0^2^l + 1 + b_0^2/b_0 = b_0^2 + 1/b_0 = b_0/b_0 = 1.

Let P = ((b_0 + 1)y)^2^l - 1 + ((b_0 + 1)y)^2^l - 2 + … + (b_0 + 1)y + 1. Then the above implies ((b_0 + 1)y - 1)P ≡ 1 J_𝔞^x. Again using (y + (b_0 + 1))(y + b_0) = (y^2 + y + 1), we multiply by (b_0 + 1)P to get

    (y + (b_0 + 1))(y + b_0) · (b_0 + 1)P ≡ 0 J_𝔞^x
       (y + (b_0 + 1)) · ((b_0 + 1)y - 1)P ≡ 0 J_𝔞^x
       y + (b_0 + 1) ≡ 0 J_𝔞^x.



Now we simply read off the dimensions using the Gröbner bases and obtain the formula

    _𝔽(𝔽[x, y]/J_𝔞^x)
       =
    2^l

Recalling the definitions, the above gives the formula

    k_𝔞, 𝔪   =
    2^l + 1


Now we wish to calculate k_𝔞 = ∑_𝔪 k_𝔞, 𝔪 where the sum is over all maximal ideals 𝔪⊂𝔽[x, y, z]/I_𝔞. For the calculation, we first need to explicitly calculate the points in the variety V(I_𝔞). By first finding solutions of h = 0, we easily calculate that the solutions of f = g = 0 is the set of points V(𝔞) = V_1(𝔞) ∪ V_2(𝔞) where

    V_1(𝔞)    = {(ζ_3 z_0' + 1, ζ_3, z_0' + 1) ∈𝔸_𝔽^3: z_0' ∈𝔽}
    
    V_2(𝔞)    = {(ζ_3^2 z_0' + 1, ζ_3^2, z_0' + 1) ∈𝔸_𝔽^3: z_0' ∈𝔽}.


Note that the intersection is V_∩(𝔞) = V_1(𝔞) ∩ V_2(𝔞) = ∅. Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0, we have the corresponding set of solutions V(I_𝔞) = V_1(I_𝔞) ∪ V_2(I_𝔞) where V_j(I_𝔞) = {(x_0, y_0, z_0) ∈ V_j(𝔞): x_0^L' = y_0^L' = z_0^L' = 1} for all j ∈{1, 2}. The intersection is V_∩(I_𝔞) = ∅. Thus, #V_∩(I_𝔞) = 0 and

    #(V_1(I_𝔞) ∖ V_∩(I_𝔞)) = #(V_2(I_𝔞) ∖ V_∩(I_𝔞)) = 
    0,     3 ∤ L' 
    (((z + 1)^L' + 1, (ζ_3 z + 1)^L' + 1)),     3 | L'.

Recalling arguments from <ref>, we put everything together using the formula

    k_𝔞   = 2^l + 1#V_∩(I_𝔞) + 2^l + 1#(V_1(I_𝔞) ∖ V_∩(I_𝔞)) + 2^l + 1#(V_2(I_𝔞) ∖ V_∩(I_𝔞)) 
       = 
    
    0,     3 ∤ L' 
    
    2^l + 2(((z + 1)^L' + 1, (ζ_3 z + 1)^L' + 1)),     3 | L'.

We also recall from <ref> that for the special case 3 | L' = 4^n - 1 for some n ∈ℤ_>0, we have the explicit calculation

    (((z + 1)^L' + 1, (ζ_3 z + 1)^L' + 1)) = 4^n - 2 = L' - 1.


We summarize the results in the following theorem.



Consider CCB_11 defined by the ideal 𝔞 = (1 + y + y^2,1 + x + y + yz). Imposing the periodic conditions x^L - 1 = y^L - 1 = z^L - 1 = 0 for some L = 2^l L' where 2 ∤ L' ∈ℤ_>0 and l ∈ℤ_≥ 0, the number of encoded qubits is given by

    k_𝔞 =
    
    0,     3 ∤ L 
    
    2^l + 2(((z + 1)^L' + 1, (ζ_3 z + 1)^L' + 1)),     3 | L.

Moreover, if 3 | L' = 4^n - 1 for some n ∈ℤ_>0, then we have the explicit formula

    k_𝔞 = 4(L - 2^l).




<ref> implies that k_𝔞 as a function of L, obeys the scaling relation

    k_𝔞(2^rL)= 2^rk_𝔞(L)

for all r ∈ℤ_≥ 0. This is consistent with the self-bifurcating behavior of CCB_11.




§ 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∩𝔽_2[x^± 2,y^± 2,z^± 2].

Let R = 𝔽_2[x, y, z] and R' = 𝔽_2[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 _R'(R/I) of the left R'-module R/I. Since _R'(R/I) = I ∩ R', we develop a method to compute I ∩ R'.

Let R̃ = 𝔽_2[x, y, z, a, b, c] and R̃' = 𝔽_2[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 <cit.>. 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

    x    ↦ x     a    ↦ x^2 
    
    y    ↦ y     b    ↦ y^2 
    
    z    ↦ z     c    ↦ z^2.

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.


We have I ∩ R' = ϕ(Ĩ∩R̃'). Hence, viewing I as an ideal generated in 𝔽_2[x^± 1, y^± 1, z^± 1], the charge annihilator after coarse-graining can be calculated by I ∩𝔽_2[x^± 2, y^± 2, z^± 2] = ⟨ϕ(Ĩ∩R̃')⟩⊂𝔽_2[x^± 2, y^± 2, z^± 2].



First we show the useful fact (ϕ) = ⟨ x^2 - a, y^2 - b, z^2 - c ⟩. Define the surjective homomorphism ϕ_a: 𝔽_2[x, y, z, a, b, c] →𝔽_2[x, y, z, b, c] uniquely determined by the mappings

    x    ↦ x     a    ↦ x^2 
    
    y    ↦ y     b    ↦ b 
    
    z    ↦ z     c    ↦ c.

and also define the surjective homomorphisms ϕ_b: 𝔽_2[x, y, z, b, c] →𝔽_2[x, y, z, c] and ϕ_c: 𝔽_2[x, y, z, c] →𝔽_2[x, y, z] in a similar fashion. Then we have the commutative diagram

    𝔽_2[x, y, z, a, b, c] [r, "ϕ_a"] [rrr, bend right = 15, "ϕ"]    𝔽_2[x, y, z, b, c] [r, "ϕ_b"]    𝔽_2[x, y, z, c] [r, "ϕ_c"]    𝔽_2[x, y, z].

In light of the containments

    𝔽_2[x, y, z] ⊂𝔽_2[x, y, z, c] ⊂𝔽_2[x, y, z, b, c] ⊂𝔽_2[x, y, z, a, b, c]

we have (ϕ) = (ϕ_c ∘ϕ_b ∘ϕ_a) = (ϕ_a) + (ϕ_b) + (ϕ_c). Now, making the identification 𝔽_2[x, y, z, c] = 𝔽_2[x, y, z][c] and using the division algorithm with the divisor being the monic polynomial c - z^2 ∈𝔽_2[x, y, z][c], we conclude that (ϕ_c) = ⟨ z^2 - c ⟩⊂𝔽_2[x, y, z, c]. By similar arguments, we also have (ϕ_b) = ⟨ y^2 - b ⟩⊂𝔽_2[x, y, z, b, c] and (ϕ_a) = ⟨ x^2 - a ⟩⊂𝔽_2[x, y, z, a, b, c]. Thus (ϕ) = ⟨ x^2 - a, y^2 - b, z^2 - c ⟩⊂𝔽_2[x, y, z, a, b, c].

We now compute that ϕ^-1(I) = I + (ϕ) = I + ⟨ x^2 - a, y^2 - b, z^2 - c ⟩ = Ĩ. Hence the surjective homomorphism ϕ|_R̃' + Ĩ: R̃' + Ĩ→ R' + I lifts to the isomorphism ϕ|_R̃' + Ĩ: R̃' + Ĩ/Ĩ→R' + I/I. We use this and the second isomorphism theorem to get

    R̃'/Ĩ∩R̃'≅R̃' + Ĩ/Ĩ≅R' + I/I≅R'/I ∩ R'.

More explicitly tracing the maps in the isomorphisms above, we have the commutative diagram

    R̃' + ĨĨ[rrr, "ϕ|_R̃' + Ĩ"]    [-20pt]       [-20pt] R' + II[ddd] 
       r̃' + Ĩ[r, mapsto]    ϕ(r̃') + I [d, mapsto]    
       r̃' + Ĩ∩R̃' [u, mapsto] [r, mapsto, dashed]    ϕ(r̃') + I ∩ R'    
    R̃'Ĩ∩R̃'[uuu] [rrr, dashed]          R'I ∩ R'.
 
We see that the isomorphism R̃'/Ĩ∩R̃'≅R'/I ∩ R' is induced by the surjective homomorphism ϕ|_R̃': R̃' → R'. This implies Ĩ∩R̃' = {r̃' ∈R̃': ϕ(r̃') ∈ I ∩ R'} = ϕ^-1(I ∩ R'). Hence we conclude ϕ(Ĩ∩R̃') = I ∩ R'.

