Towards Non-Invertible Anomalies from Generalized Ising Models

We present a general approach to the bulk-boundary correspondence of noninvertible topological phases, including both topological and fracton orders. This is achieved by a novel bulk construction protocol where solvable $(d+1)$-dimensional bulk models with noninvertible topology are constructed from the so-called generalized Ising (GI) models in $d$ dimensions. The GI models can then terminate on the boundaries of the bulk models. The construction generates abundant examples, including not only prototype ones such as $Z_2$ toric code models in any dimensions no less than two, and the X-cube fracton model, but also more diverse ones such as the $Z_2\times Z_2$ topological order, the 4d $Z_2$ topological order with pure-loop excitations, etc. The boundary of the solvable model is potentially anomalous and corresponds to precisely only sectors of the GI model that host certain total symmetry charges and/or satisfy certain boundary conditions. We derive a concrete condition for such bulk-boundary correspondence. The condition is violated only when the bulk model is either trivial or fracton ordered. A generalized notion of Kramers-Wannier duality plays an important role in the construction. Also, utilizing the duality, we find an example where a single anomalous theory can be realized on the boundaries of two distinct bulk fracton models, a phenomenon not expected in the case of topological orders. More generally, topological orders may also be generated starting with lattice models beyond the GI models, such as those with symmetry protected topological orders, through a variant bulk construction, which we provide in an appendix.


Introduction
Bulk-boundary correspondence has been a key concept for understanding topological phases of matter since the discovery of quantum Hall effect. Transport responses in quantum Hall bars are fundamentally contributed by the chiral gapless edge modes on the boundaries. [1][2][3][4][5][6] More generally, the nontrivial boundary properties of various topological phases in d + 1 spacial dimensions can be understood as a consequence of anomalies of the boundary theory in d dimensions [7][8][9][10][11][12][13][14] -obstructions in realizing a theory in a local lattice model in the d spatial dimensions with a tensor product of local Hilbert spaces, a local Hamiltonian, and an onsite symmetry action if any. The edge theory of an integer quantum Hall insulator has an invertible gravitational anomaly characterized by the imbalanced left and right moving modes. The boundary low energy effective theories for onsite-symmetry protected topological phases have 't Hooft anomalies, which prevent an onsite symmetry realization on a lattice without a topological bulk, as well as the theory to be coupled to gauge fields. Each of the above anomalies is invertible, as the anomaly can be matched by an invertible phase in one dimension formalism. Each term is a product of local terms of the GI model and its dual. By virtue of the properties of the GI models, we show the bulk model has several nice properties: • Any ground state is robust against local perturbations.
• When the GI model has non-local X -type symmetries, or when the dual model has non-local Z-type symmetries, the bulk model is either topologically ordered or fractonordered.
• When the bulk model has a pure topological order, its boundary has non-invertible anomaly. The symmetric sector of the GI model that satisfies certain (generalized) boundary conditions is a boundary termination for the bulk model.
• When the bulk model has fracton orders, it can have an anomaly-free boundary, such that when discrete global symmetries appear on the boundary, the boundary Hilbert space includes all charged sectors, rather than only the symmetric sector.
The rest of this paper is organized as follows. In Section 2, we define GI models and give a few examples. In Section 3, we introduce the generalized Kramers-Wannier duality which plays an important role in our construction of bulk models. In Section 4, we construct the bulk model, and together describe a few prototype examples. Then we prove that it has a stable spectral gap and a robust GSD on a topologically nontrivial space manifold. Hence it has a topological or fracton order. The bulk-boundary correspondence is analyzed in Section 5. In Section 6, we study a collection of examples of topological and fracton orders generated from the construction. Particularly, an interesting example demonstrates that two distinct fracton models can host the same anomalous boundary theory. In the end, we summarize and discuss future questions. In Appendix G, we also give a variant bulk construction that generates topological and/or fracton orders from qubit lattice models beyond the GI model and is applicable to some models with symmetry protected topological orders. Further technical details are also summarized in appendixes.

Generalized Ising models
A GI model is referred to a model on a lattice of qubits in arbitrary spatial dimensions, whose Hamiltonian and 2 symmetries has the following properties. The Hamiltonian consists of two types of terms: GI terms and generalized transverse field terms. A generalized Ising (transverse field) term is a product of Pauli-Z (Pauli-X ) operators acting on a local subset of qubits, and is denoted by O Z α (O X i ) with some index α (i) referring the subset. Generically, α and i are from different index sets. Written explicitly, the Hamiltonian is then where J α and h i are real coefficients. We suppose the model lives on a d-dimensional parallelogram with either periodic or open boundary condition along any direction. We will impose some additional assumptions on the model later in this section. The model may have many 2 symmetries. For our purpose, we only consider one group of 2 symmetries: all Z-type symmetries generated by products of Pauli-Z operators (minus sign factor excluded), and all those X -type symmetries generated by Pauli-X operators (minus sign factor excluded) that commute with all Z-type symmetries. We refer this selected group of symmetries the compatible symmetries in the GI model. Compatibility is to emphasize that the generator of each X -type 2 symmetry commutes with not only the Hamiltonian but also all the Z-type symmetry generators. In fact, this implies that each X -type symmetry generator is a product of several O X i operators, analogous to the standard transverse-field Ising model. For a proof, see Corollary 5 in Appendix B. From now on, X -type 2 symmetries in the GI model only refer to those compatible ones.
The many 2 symmetries may either be local or nonlocal, and it is useful to distinguish them for our purpose. Let {G Z r } with some index r be a complete but not necessarily independent set of generators of the local Z-type symmetries. We say there are n Z number of nonlocal Z-type symmetries if one can find a maximal set of symmetry generators {U Z 1 , U Z 2 , · · · , U Z n Z } satisfying that each U Z k is not a product of the remaining ones and G Z r . Formally, this n Z 2 group is nothing but the quotient of the full Z-type symmetry group over its local symmetry subgroup. In practice, the set {U Z k } can be obtained by repeatedly adding new U Z k that is independent from G Z r and the existing U Z 1 , · · · , U Z k−1 until the list is maximal. Similarly, for the X -type symmetries, we have the local generators {G X s } with some index s which belongs to a generically different index set from that for r, and the independent nonlocal generators {U X 1 , U X 2 , · · · , U X n X } with some n X .

Assumptions
Before stating our assumptions, let us introduce some useful terminologies. Consider a set of commuting local operators {M i } where each M i is a tensor product of I, X , Y, Z. We say that a local operator A is locally generated by {M i } if A is generated by a few M i 's in a neighborhood of A's support, such that the linear size of this neighborhood exceeds that of A by an O(1) constant. We say that {M i } is a complete set of local observables (CSLO) if any local operator A that is a tensor product of I, X , Y, Z and commutes with all M i can be locally generated by {M i }. In fact, one can show that if {M i } forms a CSLO, then any local operator A, not necessarily a product of Paulis, that commutes with all M i can be locally generated. We assume the GI models to have the following properties.
• {O X i } ∪ {G Z r } is a CSLO.
• Any local X -type symmetry generator is locally generated by {G X s }.
The first assumption physically means that when J α = 0, h i ̸ = 0, after restricting to the gauge invariant sector G Z = G X = 1, the system has a spectral gap stable to local perturbations, together with either a unique ground state or a robust GSD [72]. This is analogous to the Ising disordered phase. Independent of these assumptions, the bulk model to be constructed has a Hamiltonian whose local terms all commute. These two assumptions ensure that the ground states of the bulk model to be constructed are robust against local perturbations. The above assumptions may seem technical, but in many cases, are not hard to verify, as we will see. Also note that the choices of G Z and G X operators are not unique. It suffices to make one choice that satisfies the assumptions.

Examples
Let us now introduce some examples. Periodic boundary condition will be taken for convenience. A particularly simple situation is when O X i are just the traditional transverse field terms X i with i labeling the qubits on the lattice. In this case, there is no Z-type symmetry at all, and our first assumption is trivially satisfied. The simplest example of this class is of course the standard one-dimensional Ising model: H = −J i Z i Z i+1 − h i X i , which has an X -type references therein), and the quantum Newman-Moore model [76][77][78][79][80] are other examples of this class.
The two-dimensional plaquette Ising model has the Hamiltonian defined for qubits on the vertices of a 2D square lattice. Here and throughout, we sometimes suppress the summation for simplicity when there is little confusion. The product of Pauli X operators along each row and column generates a 2 symmetry of the model, known as a subsystem symmetry. Point excitations of this model in the ordered phase (J > h) has restricted mobility. The quantum Newman-Moore model has subsystem 2 symmetries acting on Sierpinski triangles, but is otherwise similar to the two-dimensional plaquette Ising model, so let us not write it down explicitly. Our next example contains nontrivial O X terms in its Hamiltonian. Consider the following model whose qubits live on the links of a 2D square lattice, Here, the nearest neighboring two-body Z-type terms along the vertical links are not included as they are not independent: they are equivalent to the two-body Z-type terms along the horizontal directions up to a local Z-type symmetry of the Hamiltonian. The local Z-type symmetries of this model are generated by the product of four Z operators around each vertex. Take the local Z-type symmetries as gauge constraints, the model is a quantum 2 gauge theory, and found to arise in the system of Josephson arrays of superconductor and ferromagnet when deposited on top of a quantum spin Hall insulator [81,82].
There are two nonlocal Z-type symmetries, and we may take U Z 1 (U Z 2 ) to be the product of all vertical-link (horizontal-link) Z operators along some horizontal (vertical) line, see Fig. 1a. There are no local X -type symmetries. Each X -type non-local symmetry generator of the model is a product of all vertical-link X operators along even number of vertical lines, and all horizontal-link X operators along even number of horizontal lines. See Fig. 1b for an example.
Later, we will construct a bulk theory for this model, which has the X-cube fracton order [59].

Generalized Kramers-Wannier duality
In this section, we define generalized Kramers-Wannier dual theories for each GI model. Such dual theories play an important role in our construction of the bulk models. A dual theory lives on a generically different lattice which we dub the dual lattice. The operator map of the duality can be written as where ∆ Z α (∆ X i ) is a local product of Pauli Z (Pauli X ) operators on the dual lattice, such that the commuting or anticommuting relations between the operators are preserved. Moreover, the above operator map should be local: if we place the original and the (generically different) dual lattices together, then each O Z α (O X i ) operator should be closed to the corresponding ∆ Z α (∆ X i ) operator. The dual model Hamiltonian then reads Such a duality exists for any GI model, because we can always let the dual lattice consist of qubits labeled by α, and then let i . This is dubbed the standard dual theory. We may treat the dual theory as a GI model as well, but with the roles of X and Z exchanged, which means we first include all X -type 2 symmetries, and then include all compatible Z-type symmetries. Similarly as in the GI model, here, all Z-type symmetries are generated by products of Hamiltonian local terms ∆ Z i . We denote the local symmetry generators in the dual theory by {Γ X ρ } and {Γ Z σ }. The independent nonlocal symmetry generators are denoted as In Appendix B, we prove that if we restrict to the symmetric sectors on both sides of the duality, then the operator map (4) follows from a Hilbert space isomorphism, i.e. an exact duality.
Similar to the original theory, we make the following assumptions for the dual theory: • Any local Z-type symmetry generator is locally generated by {Γ Z σ }.
In addition, we assume that • n X + m Z ≥ 1.
In other words, either there exist compatible nonlocal X -type symmetries in the original model, i.e. n X ≥ 1, or there exist compatible nonlocal Z-type symmetries in the dual model, i.e. m Z ≥ 1. This will help ensure our bulk model to have a topological and/or fracton order. Later in the Section 5, we discuss the further conditions on the GI model (and its dual) so that the GI model has non-invertible anomaly that can be matched with the bulk model to be constructed. For example, the standard dual theory of the standard one-dimensional Ising model is where we place the dual lattice qubits in between the original ones, reflected by the 1/2 shifts in the indices. Another example is that the standard dual theory of both the two-dimensional plaquette Ising model in Eq. 2 and the model in Eq. 3 is which is nothing but the two-dimensional plaquette Ising model with the substitutions X ↔ Z and J ↔ h. This dual theory has no local symmetry. We emphasize that a single GI model may have multiple dual models. Just from the example above, another possible dual theory of the two-dimensional plaquette Ising model is Eq. 3 with the substitutions X ↔ Z and J ↔ h. As a consequence, multiple bulk models may be constructed from a single GI model, as we will see.

Bulk theory
Given some GI model in d spatial dimensions and a dual model of it, we will now construct a bulk theory in one higher dimensions such that certain charge and boundary condition sector(s) of the GI model can live on its boundary. We will explain later what this precisely means.

Construction and prototype examples
The lattice on which the bulk theory lives is an alternating stack of the original and dual ddimensional lattices; see Fig. 2. As an example, we also show the bulk lattice thus constructed from the standard one-dimensional Ising model and its standard dual model in Fig. 3. Here and throughout, we often use empty circles (solid dots) to represent qubits in layers of the original (dual) lattice. We label the original and dual lattice layers by odd and even indices, then our bulk theory is defined by the following Hamiltonian; see Table 1 for a recap of the many notations. . . . . . . where the second subscript of each operator is the layer index. These various operators are schematically plotted in Fig. 2 terms will be called the Zsuspension (X -suspension) terms; this name is from the special case where ∆ Z = Z (O X = X ). The G and Γ terms will be called the gauge symmetry terms. By construction, all local operators in H bulk commute with each other, thus the model is exactly solvable. In other words, H bulk is a stabilizer Hamiltonian.
The simplest example comes out starting from the standard one-dimensional Ising model and its standard Kramers-Wannier dual. We obtain the following bulk Hamiltonian, which lives on the lattice shown in Fig. 3. There are no gauge symmetry terms in 8. The Hamiltonian represents nothing but the 2 toric code model [45]; it can be cast to the standard form by replacing empty circles and solid dots by vertical and horizontal links, respectively. Similarly, the three-dimensional toric code model can be generated from the standard twodimensional Ising model and its standard dual, the 2 lattice gauge model without matter.
From the two-dimensional plaquette Ising model (2) and its standard dual (6), we obtain where the 3D Cartesian frame is indicated in the bracket with z being the out-of-layer direction. The model is also obtainable via other constructions [83,84]. Point excitations in this bulk model are free to move along z direction, but have restricted mobility along x and y directions like the original two-dimensional model. In other words, the model is fractonic along x and y, but behaves like a topologically ordered system along z. Another example, from the two-dimensional model in Eq. 3 and its standard dual theory in Eq. 6, we obtain where we have associated qubits in the dual lattice layers with z-direction links, and the Cartesian frame is the same as that in Eq. 9. This three-dimensional model is topologically equivalent 5 to the X-cube fracton model [59]. Recall that the two-dimensional plaquette Ising model has an alternative dual theory: Eq. 3 with the substitutions X ↔ Z and J ↔ h. The bulk theory constructed from this pair of models is the same as Eq. 10 but with X and Z exchanged. We have thus found that multiple bulk models may be constructed from a single GI model by choosing different dual models.

Robust ground state degeneracy
We now show that the general bulk theory H bulk has a stable spectral gap, and any possible GSD of it is robust. Therefore any degenerate ground states, say on the lattice in d ≥ 2 dimensions with periodic boundary condition, would imply topological and/or fracton orders. Then we compute the GSD.
Regarding H bulk as the negative sum over a set of stabilizers, then in the ground subspace, all these stabilizers equal to +1, i.e. there is no frustration. 6 According to Ref. [72], the following lemma implies that the model has a spectral gap stable to local perturbations, together with either a unique ground state or a robust GSD. 7 In particular, the lemma shows that the stabilizer Hamiltonian constructed is a quantum code with macroscopic code distance. The logical operators, if any, are all non-local operators that commutes with the stabilizers.

Theorem 1. The stabilizers in H bulk form a CSLO. In other words, any local operator
A that is a tensor product of I, X , Y, Z and commutes with all the stabilizers in H bulk can be locally generated by those stabilizers.
Proof. Up to an unimportant phase factor, we can write A = A Z A X where A Z (A X ) is a product of Pauli Z (Pauli X ) operators on different sites. A Z and A X must themselves be local and commute with all the stabilizers in H bulk , because all the stabilizers in H bulk are either X -type or Z-type. We will show that A Z and A X are both local products of the operators appearing in H bulk .
Each Z operator in A Z has some integer layer index l. Let the maximal and minimal ones of those layer indices be l max and l min , respectively. Suppose l max is odd, i.e. coinciding with an original lattice layer. In order to commute with all the X -suspension terms that span the three layers l max , l max+1 and l max+2 , and also by our assumptions on the original d-dimensional model, the top layer of A Z must be a local product of some G Z terms. Thus we may multiply A Z by those G Z terms and reduce l max by at least 1. Now suppose that l max is even. The top layer of A Z coincides with a dual lattice layer, and commutes with all the Γ X operators on that layer. Given our assumptions on the dual d-dimensional model, it follows that the top layer of A Z must be a local product of some ∆ Z operators. Let h = l max − l min be the height of the A Z operator. Whenever h ≥ 2, we may multiply A Z with some Z-suspension operators and reduce its height by at least 1. Repeat the above two operations to decrease l max , and the similar operations to increase l min . Eventually, the reduced A Z operator acts on a single dual lattice layer (even layer index), if it is not yet fully reduced to the identity. This single-layer operator is a product of some ∆ Z operators, and must commute with all the ∆ X operators acting on that layer (due to the X -suspension operators), thus it is a local Z-type symmetry generator of the dual theory and is a local product of some Γ Z operators by our assumptions on the dual model. Given the locality of the O X , O Z , ∆ Z , and ∆ X operators, as well as the locality of the generalized Kramers-Wannier operator map, the above reduction procedure implies that A Z is locally generated by the stabilizers in H bulk . The claim for A X can be proved similarly.
Next, we examine in what cases the model has GSD. It turns out the GSD, for example, on a (d + 1)-dimensional torus, only depends on properties of non-local symmetries in both the original (generalized Ising) model and its dual. Specifically, we take periodic boundary condition along the out-of-layer direction, i.e. identify the first and the L-th layers for some odd L. We obtain the GSD through a thorough counting, which we elaborate in Appendix D.
Here instead, let us prove the existence of a degeneracy by finding a pair of anticommuting operators that act on the ground subspace. As a general mathematical fact, given a set of independent X -type operators A X 1 , A X 2 , · · · , A X n acting on an arbitrary multiple-qubit Hilbert space, there is always a set of Z-type operators B Z 1 , · · · , B Z n such that B Z k anticommutes with A X k but commutes with all the other A X operators. 8 This means that we can always find some operators V Z k acting on the Hilbert space of our original d-dimensional lattice, such that V Z k anticommutes with U X k but commutes with all the other X -type symmetry generators (local or nonlocal). Let us then consider the operator U X k,2l 0 −1 for some k and l 0 , which is U X k acting on an original lattice layer 2l 0 − 1, and the operator where V Z k,2l−1 is V Z k acting on the (2l −1)-th layer. The two operators U X k,2l 0 −1 and W k both commute with H bulk , thus acting within the ground state subspace, and the two operators anticommute with each other. It follows that the ground state subspace can not be one-dimensional. Indeed, let |ψ〉 be a ground state that is also an eigenstate of U X k,2l 0 −1 with some eigenvalue λ = ±1, then W k |ψ〉 is another ground state with eigenvalue −λ under the U X k,2l 0 −1 operator. This analysis actually implies that the U X k,2l 0 −1 operators for some fixed l 0 and all k can take independent eigenvalues ±1 within the ground subspace. Similarly, the Ω Z k,2l 0 operators can also take independent eigenvalues within the ground subspace. Hence, the GSD of the system is at least 2 n X +m Z . Our previous assumption n X + m Z ≥ 1 guarantees a nontrivial degeneracy.
If n X increases with the system size, which, for example, happens in the plaquette Ising model and the model in Eq. 3, the bulk theory will have a fracton order, with a GSD increases with the system size. The scenario for m Z is similar.
A few simple cases are illuminating, following the formula of GSD given in Theorem 3, which we summarize in the following claims. Corollary 1. When the GI model has n X ≥ 1 number of (compatible) X -type non-local symmetries, the bulk model has degenerate ground states stable against local perturbations, and Corollary 2. When the dual of the GI model has m Z ≥ 1 number of (compatible) Z-type non-local symmetries, the bulk model has degenerate ground states stable against local perturbations, and

Bulk-boundary correspondence
Now let us analyze the boundary physics of the bulk model constructed above, and show that the GI model we start with can terminate on its boundary. This more precisely means the following: When the bulk model is placed on a certain space with boundary, in the absence of bulk excitation, its low-energy physics below the bulk excitation gap is described by the GI model subject to certain global constraints. The bulk-boundary correspondence manifests locally, through the matching of Hamiltonian local terms, as well as globally. Especially, we focus on matching the following two types of global constraints on the GI model to the boundary of the bulk model.
• Global symmetry charge projection: for a set of indices N ⊂ {1, 2, · · · , n X }, • Generalized boundary condition: for a set of indices S such that α∈S O Z α = 1 modulo To understand why we call the latter condition a generalized boundary condition, let us apply the condition to one-dimensional tranverse Ising model on a ring. In this model, {O Z α } can be identified with {Z 1 Z 2 , Z 2 Z 3 , · · · , Z N Z 1 }, such that α O Z α = 1, and one way to change the boundary condition is to flip the sign of the coefficient of the Z N Z 1 term in the Hamiltonian; and correspondingly, the sign of η = α sign(J α ) is changed.
There are two main results. Consider a bulk model with a finite number of layers, and with the top and bottom being odd layers. The first is that the boundary Hilbert space L bdry of the bulk model with topological and/or fracton orders is isomorphic to that of two copies of GI models subject to global symmetry constraints. Particularly, under a condition to be specified below, L bdry is isomorphic to the sector labeled by +1 eigenvalues of at least one of the following two types of non-local operators. One is dubbed global symmetry charges, and the other is dubbed generalized boundary conditions. The notation here deserves some explanation. U X N for N ⊂ {1, 2, · · · , n X } are global symmetry operators of the GI model, and hence the operator in Eq. 16 does divide the Hilbert space of two GI models into different eigenvalue sectors. On the other hand, Ω Z M for M ⊂ {1, 2, · · · , m Z } are symmetry operators of the dual model, then how does the operator in Eq. 17 acts on two copies of the original model? In fact, we will show that the eigenvalues of certain Ω Z M in the dual model imply generalized boundary conditions given by (15), through the duality map. Therefore, Ω Z M ⊗ Ω Z M actually acts on two copies of the original model as η S ⊗ η S for some S.
Furthermore, we find that the boundary Hilbert space L bdry can be divided into many sectors labeled by the eigenvalues of some nonlocal operators, which are all X -type or Z-type symmetry operators of the original or dual d-dimensional models acting on certain layers. The sectors are all isomorphic, and the boundary Hamiltonian H bdry is block-diagonal with respect to these sectors. In different sector, the charge projections (14) and generalized boundary conditions (15) either on the top layer or on the bottom layer may be different. Nevertheless, the combinations of their values on both the top and the bottom layer need to be consistent with that U X N ⊗ U X N = 1 and Ω Z M ⊗ Ω Z M = 1. In this way, the boundary model of a non-invertible phase with long-range orders is described by the GI model with global constraints. When this happens, we dub the constrained GI model to have (weak) non-invertible anomaly.
The second main result is a concrete condition on the non-local symmetries U X N , Ω Z N of the original and dual d-dimensional models that leads to a bulk model whose boundary theory matches with the GI model with constraints (14) and/or (15). We summarize it in the following claim.

Claim 1. (Necessary and sufficient condition for an anomalous boundary)
The boundary theory is anomalous if and only if either of the following two conditions is satisfied: 1. For some nonempty subset N ⊂ {1, 2, · · · , n X }, can be written as a product of O X operators such that its dual -a product of ∆ X operators -equals to the identity modulo local symmetry operators Γ X .
2. For some nonempty subset M ⊂ {1, 2, · · · , m Z }, can be written as a product of ∆ Z operators such that its dual -a product of O Z operators -equals to the identity modulo local symmetry operators G Z .
Since that a product of a few O Z equals the identity modulo G Z 's is a generalized boundary condition (15) in the GI model, and there is an analogy for a product of ∆ X . Thus, colloquially speaking, the conditions say that only for those non-local symmetry operators in either the GI or the dual model, which is dual to a generalized boundary conditions, the projections of them lead to the (weak) non-invertible anomaly.
In the special case where O X = X and ∆ Z = Z, the condition is simple, that is n X + m Z ≥ 1. The reason is that in this case, any non-local symmetry U X or Ω Z is dual to the identity, because the original (dual) theory does not have any Z-type (X -type) symmetry.
The following subsections are devoted to analyze the boundary theory from the simplest to the most general case, which leads to the results above. The examples of the boundary theory of toric code model and the X-cube model are presented. A complete and detailed treatment is given in Appendix E.
To begin with, let us define the boundary of our bulk model. We will take an odd number of layers, with layer indices from 1 to L ∈ 2 + 1, and take open boundary condition along the out-of-layer direction, so the 1-st and the L-th layers are the boundary layers. The two boundary layers both have the original (instead of the dual) lattice structure on which our GI model is defined, cf. Fig. 2 and 3. How about the boundary condition along the intra-layer directions? Previously, we have assumed periodic boundary condition when discussing concrete examples, but our construction does not really demand any particular boundary condition. In the following, we just require that the boundary condition for each original (dual) lattice layer be the same as the original (dual) d-dimensional model, but is otherwise arbitrary. However, we emphasize that if one changes the boundary condition for a GI model, its symmetry, the dual model, and the validity of our assumptions should all be reexamined. An example will be given below.
We define the bulk Hamiltonian H bulk to be of the same form as (7), including all the terms that are completely inside the system. This Hamiltonian determines a degenerate ground state subspace, which we consider as the boundary Hilbert space L bdry . All additional local operators that commute with local terms in H bulk , and thus act within the boundary Hilbert space L bdry , are allowed terms in the most general boundary Hamiltonian H bdry . L bdry together with H bdry is the boundary theory that we are going to determine. Note that there is not a unique choice of H bdry , since the product of any two boundary terms is another allowed boundary term. Instead, we will focus on a canonical choice of H bdry . We prove in Appendix E.5 that the boundary terms given in the canonical choice together with the stabilizers in the bulk Hamiltonian are sufficient to generate any allowed boundary local term. Thus the canonical H bdry we consider is a quite general one.

Simplest situation: O X = X and ∆ Z = Z
Let us start with the simplest situation where O X i = X i and ∆ Z α = Z α , with i and α labeling qubits in the original and dual lattices, respectively. In this case, the original model (the dual model) has no Z-type (X -type) symmetry at all. With periodic boundary condition along the out-of-layer direction, the GSD is determined by the number of non-local X -type symmetries in the original model and the number of non-local Z-type symmetries in the dual model, independent of the number of layers, log 2 GSD = n X + m Z .
Note that one obvious type of operators that commute with H bulk is the nonlocal Z-type symmetry operator Ω Z k,2l 0 in any even layer. Thus, we can divide L bdry into several sectors labeled by the eigenvalues of the nonlocal operators Ω Z k,2l 0 (k = 1, · · · , m Z ) for some fixed l 0 .
Notice that under the Kramers-Wannier duality, Ω Z k from the dual theory corresponds to the identity operator of the original theory. Hence, Ω Z k,2l 0 is related to Ω Z k,2l for any other l by the multiplications of several Z-suspension operators. More explicitly, . This is the reason that we only need to consider Ω Z k,2l 0 operators acting on a single layer. Let L bdry,0 ⊂ L bdry be the particular sector where Ω Z k,2l 0 = 1 for all k. Denote by L the Hilbert space for our original d-dimensional lattice and by L G the gauge invariant subspace of it (a.k.a., the symmetric subspace for all local symmetries).
We claim that L bdry,0 is isomorphic to the following fictitious space, where the two copies of L G represent the two boundary layers of our physical system. That is, L bdry,0 is the Furthermore, L bdry,0 is an invariant subspace of H bdry whose action in this sector, when represented in L fic , can take the form where H I GI and H II GI act on the two copies of L G in Eq. 20, respectively. Let us understand the result of the boundary Hilbert space first. Note that the Pauli-X operator acting on any qubit in the two boundary layers commute with the bulk Hamiltonian. States in L bdry,0 can be labeled by the eigenvalues of these Pauli-X operators, subject to the following two constraints. First, each local X -type symmetry generator equals to 1, since the generator is a local term in the bulk Hamiltonian. This constraint gives rise to the gauge invariance requirement in Eq. 20. Second, since U X k is dual to the identity operator under the Kramers-Wanner duality, 9 U X k,1 U X k,L is equal to the product of several X -suspension operators, and thus, is equal to 1. This constraint leads to the n X 2 symmetry projection. Now we consider the boundary Hamiltonian local terms. The Pauli-X operators on the two boundary layers are allowed, since, as just mentioned, they commute with the bulk Hamiltonian. Furthermore, these operators commute with Ω Z k,2l 0 , and thus act within each sector of the boundary Hilbert space. Under the isomorphism from L bdry,0 to L fic , these operators take the same form, Another set of operators that can be added to They all commute with the bulk Hamiltonian and with Ω Z k,2l 0 as well. Their image in L fic is, This map may seem obvious since it preserves the commuting/anticommuting relations with the boundary X operators, but a careful proof is actually necessary. For example, extra minus sign factors also seem allowed and it is not immediately clear whether they can be gauged away. Our proof for this result is given in Appendix E.2. A crucial ingredient in the proof is that Ω Z k,2l 0 = 1 for all k; one can see this by noticing that each Ω Z k,2 is equal to the product of several O Z α,1 Z α,2 . With the above analysis, we conclude that the L bdry,0 block of H bdry , when represented in L fic , may take the form (21), where H I GI and H II GI act on the two copies of L G in Eq. 20, 9 U X k is a product of several X i operators. One can obtain a dual operator by applying the Kramers-Wannier operator map (4). Such a dual operator has to commute with all the Z α operators, so it must be the identity.
respectively. We will refer to Eq. 21 as the effective boundary Hamiltonian in the L bdry,0 sector, where the "effectiveness" is in the sense that the Hamiltonian acts on the fictitious space L fic .
Other sectors with different eigenvalues of Ω Z k,2l 0 can be analyzed by considering unitary operators that map them to L bdry,0 . Denote by L ′ the Hilbert space for the d-dimensional dual lattice. We can find some X -type operators Θ X k acting on L ′ such that Θ X k anticommutes with Ω Z k but commutes with all the other Z-type symmetry generators (local or nonlocal); this is always possible as we mentioned earlier. It follows that is an operator that commutes with the bulk Hamiltonian and can flip the eigenvalue of Ω Z k,2l 0 . Hence, we have found that each of the m Z sectors of L bdry is isomorphic to L fic defined in Eq. 20. The above unitary operator that can alter the sign of Ω Z k,2l 0 commutes with all the O X = X operators on the two physical boundaries, but necessarily anticommute with some In consequence, the effective boundary Hamiltonian in each of the other sectors still takes the form of Eq. 21, but the signs of some GI terms in both H I GI and H II GI are flipped compared to those in the L bdry,0 sector.
Crucially, such sign changes cannot be canceled by any unitary rotation in L fic . To see this, we write Ω Z k = α∈A Z α for some subset A. Then under the generalized Kramers-Wannier duality map, 10 It leads to that in an arbitrary sector of L bdry , In L fic , we have Suppose there is an isomorphism from this sector of L bdry to L fic , such that then we necessarily have It means that as we go from one sector to another with a different Ω Z k,2l 0 , some of the η α must change signs! A similar statement holds for the Z α,L−1 O Z α,L operators. The Hamiltonian in this sector, is isomorphic to the following one in L fic , We have seen that H bdry has a block-diagonal action on ⊕ a L bdry,a , where a is the sector index. In fact, there are no local operators (but only non-local ones) that can map between sectors, as well as commuting with all bulk Hamiltonian local terms. This is because any local operator commuting with the bulk Hamiltonian local terms can be generated by the set of boundary local terms considered above, as we mentioned previously and proved in Appendix E.5.
Let us now discuss some examples. Consider the two-dimensional toric code model constructed in Eq. 8 as a bulk theory for the standard one-dimensional Ising model with periodic boundary condition. Following our prescription, we shall take open boundary condition along the vertical direction, while keep periodic boundary condition along the horizontal direction; see Fig. 3. The Ising model has only one U X operator, given by U X = i X i , and similarly, its standard dual model has only one Ω Z operator given by Ω Z = i Z i+1/2 . Thus, from our discussion above, L bdry can be divided into two sectors with Ω Z 2l 0 = ±1. Each of the two sectors can be regarded as two spin chains, subject to the symmetry condition U X ⊗ U X = 1. The boundary Hamiltonian in one of the two sectors may be the sum of two Ising model Hamiltonians acting on the two fictitious spin chains, respectively, and with periodic boundary condition. Then the boundary Hamiltonian in the other sector will again be the sum of two Ising model Hamiltonians, but now with antiperiodic boundary condition, i.e. one Ising term on each of the two spin chains changes its sign.
To give a complementary perspective, we may alternatively start from the standard onedimensional Ising model with open boundary condition, namely defined on a chain of N spins labeled by 1, 2, · · · , N . Again, the model has one 2 symmetry defined on a chain of N − 1 spins labeled by 3/2, 5/2, · · · , N − 1/2. This dual model has no 2 symmetry at all! One can construct the bulk theory accordingly, which now lives on a lattice with left and right boundaries. We can take open boundary condition along the vertical direction as well, and analyze its boundary theory with the result established above. We see that the boundary Hilbert space contains only one sector, and can be regarded as two disconnected open spin chains under a 2 symmetry projection U X ⊗ U X = 1. The boundary Hamiltonian may take the form of an Ising Hamiltonian on each of the two chains. The two effective open spin chains are not connected because our formalism does not allow any boundary terms on the left and right boundaries. This is actually just a matter of choice. We may redefine the bulk Hamiltonian by removing certain terms near the left and right boundaries. This will enlarge the boundary Hilbert space a bit, and allow boundary terms acting on the left and right boundaries. One can check that, with a rectangular geometry, the low-energy physics of the toric code model can be a one-dimensional Ising model defined on a closed chain with periodic boundary condition and the 2 even projection, as discussed in Refs. [36,37,47]. We have seen that when n X ≥ 1, there are the symmetry projections U X k ⊗ U X k = 1, k = 1, · · · , n X . When m Z ≥ 1, the boundary conditions for H I GI and H II GI in the effective boundary Hamiltonian will simultaneously change as we alter the values of the Ω Z k,2l 0 operators. Either phenomenon implies the boundary theory to be anomalous. Conversely, when n X + m Z = 0, the whole boundary Hilbert space L bdry is simply isomorphic to L G ⊗ L G , on which the effective boundary Hamiltonian takes the form of Eq. 21, or more generally consists of local terms generated by those in Eq. 21. This is a nonanomalous theory. Therefore, we reach the following conclusion.
Necessary and sufficient condition for an anomalous boundary. In the special case where O X = X and ∆ Z = Z, non-invertible anomaly on the boundary exists if and only if or equivalently, the GSD of the bulk model with periodic boundary condition along the out-oflayer direction satisfies GSD > 1 .
That is, we can always build a bulk model on a lattice with odd layers, such that its boundary has non-invertible anomaly, that can be matched by a GI model with distinct sectors of Hilbert space. The equivalent condition (32) follows from that in the case O X = X and ∆ Z = Z, GSD = 2 n X +m Z .

Less simple situation: ∆ Z = Z
Next, we consider the less simple situation where O X i are general but ∆ Z α = Z α . That is, we adopt the standard dual theory. This includes the model in Eq. 10 with the X-cube fracton order that we constructed as a bulk theory for Eq. 3.
The analysis of the boundary theory is similar to the simplest case. We track how non-local operators Ω Z m , U Z n and U X n ′ manifest on the boundary L bdry . Again, we divide L bdry into different sectors. These sectors are now labeled by the eigenvalues of not only Ω Z k,2l 0 The L bdry,0 sector, defined by Ω Z k,2l 0 = 1 and U Z k,2l−1 = 1 for all the internal layers, is isomorphic to a fictitious space of the same form as Eq. 20. Now, states in L bdry,0 are labeled by the eigenvalues of all the O X and U Z operators acting on the two boundary layers, subject to the constraints G X s,1 = G X s,L = 1 and U X k,1 U X k,L = 1. As in the previous situation, U X k,1 U X k,L is generated by the X -suspension operators. The operator map from L bdry,0 to L fic is again thus the L bdry,0 block of H bdry takes the same form as Eq. 21.
Other sectors of L bdry can again be analyzed by establishing isomorphisms to L bdry,0 , and thus to L fic . The eigenvalue of each Ω Z k,2l 0 can be adjusted without affecting any U Z k,2l−1 by the operator k,2l whose definition is the same as that in Section 5.1. The eigenvalues of the internal-layer U Z operators can be altered with some X -type operators that commute with not only the bulk Hamiltonian, but also with O Z α,1 Z α,2 and Z α,L−1 O Z α,L , and thus act trivially on the effective boundary Hamiltonian. That is, eigenvalues of U Z labels extra degeneracy of the boundary model that are unrelated to symmetry charge projections or boundary conditions. More explicit description of such operators is given in Appendix E.2.
Suppose for a nonempty subset M ⊂ {1, 2, · · · , m Z }, Ω Z M ≡ k∈M Ω Z k is dual to the identity modulo the G Z operators. One can show that altering the eigenvalue of Ω Z M,2l 0 will necessarily flip the signs of some O Z terms in both H I GI and H II GI , with a proof essentially the same as that in Section 5.1. We also prove in Appendix E.2 that, if such M does not exist, and at the same time n X = 0, the boundary theory is nonanomalous. That is to say that the boundary theory in this case is a direct sum of identical sectors. The Hilbert space of each sector is isomorphic to L G ⊗ L G with the operator mapping rule in Eq. 33. The effective boundary Hamiltonian of each sector takes the form of Eq. 21, or more generally consists of local terms generated by those in Eq. 21.
Necessary and sufficient condition for an anomalous boundary. We conclude that, in the special case where ∆ Z = Z, the boundary theory is anomalous if and only if either of the following two conditions is satisfied: 1. n X ≥ 1, so that there are the symmetry charge constraints U X k ⊗ U X k = 1. Figure 4: Examples of (a) a U Z operator, (b) a U X operator, and (c) an Ω Z operator of the bulk model in Eq. 10, all viewed from z direction.

For some nonempty subset
The boundary of the X-cube model. Now we illustrate with the example of the model in Eq. 10, which describes the X-cube fracton order, and is constructed from Eq. 3 and its standard dual in Eq. 6. We put the model on a 3D cubic lattice with L x × L y × L z number of vertices, with periodic boundary condition along x and y, and open boundary condition along z. The two boundary surfaces are "smooth", and L is related to L z by L = 2L z − 1 (the height of the system is L z −1 number of lattice constants). We label the lattice vertices by integer coordinates r = (x, y, z) ∈ 3 such that x ∼ x + L x , y ∼ y + L y , and 1 ≤ z ≤ L z . We denote by Z(r ; ∆r ) the Pauli Z operator acting on the link connecting the two neighboring vertices r and r + ∆r with ∆r =x,ŷ,ẑ; similar for Pauli X operators. Symmetries of the two-dimensional models have been described in words previously. A careful analysis of degeneracy relations shows that n Z = 2, n X = L x + L y − 2, m X = 0, and m Z = L x + L y − 1. More explicitly, the U Z operators acting on the (2z − 1)-th layer can be chosen as and L y y=1 Z(x 0 , y, z;x) for some x 0 .
See Fig. 4a for an example. The U X operators acting on the (2z − 1)-th layer can be chosen as and where we have excluded y = L y in (36) and x = L x in (37) because they are not independent. See Fig. 4b for an example.
The Ω Z operators acting on the 2z 0 -th layer can be chosen as and where we have excluded x = L x in (39) because it is not independent. See Fig. 4c for an example. According to our general theory, L bdry can be divided into 2 m Z +n Z (L−3)/2 number of sectors. The boundary theory in L bdry,0 may be two copies of the model in Eq. 3, subject to the n X number of symmetry projections U X k ⊗ U X k = 1. Other sectors are all isomorphic to L bdry,0 , and the isomorphisms may flip the signs of some GI terms in both H I GI and H II GI in the effective boundary Hamiltonian. More explicitly, the eigenvalue of each U Z operator acting on the (2z − 1)-th layer with 3 ≤ 2z − 1 ≤ L − 3 can be independently flipped by the operators and which correspond to the two operators in (34) and (35), respectively. It is not generically true that the eigenvalues of the internal-layer U Z operators can be adjusted with single-layer operators as in this example, but it is indeed true that these adjustment operators can be chosen to commute with H bdry . Given some fixed even layer 2z 0 , one can independently flip the signs of the Ω Z operators acting on this layer by the string operators and which are in one-to-one correspondence with the operators in (38) and (39). Each of the above anticommutes with some terms in H bdry . For example, the string opera- . One can verify that the two conditions for anomaly are both satisfied, thus the boundary theory of this model is indeed anomalous.

The most general situation
Now we briefly discuss the most general situation: no further assumption on either O X or ∆ Z .
We may again divide L bdry into several sectors, which are now labeled by the eigenvalues , Ω X k,2l for all l, and Ω Z k,2l 0 for some fixed layer 2l 0 , as they are non-local operators commuting with the bulk Hamiltonian. A caveat is that the eigenvalues of either the U Z operators or the Ω X operators may not be totally independent. It may happen that a product of several Z-suspension operators centered on some internal odd layer 2l −1 equals to a nonlocal Z-type symmetry generator (independent of G Z 's) acting on that layer. This will induce some relation between the U Z operators on the same internal layer. In terms of the d-dimensional GI model, this means that a product of several O Z operators equals to a nonlocal Z-type symmetry generator (independent of G Z 's) and is dual to the identity. A similar possibility exists for the Ω X operators. These possible degeneracy relations reduce the apparent number of sectors in L bdry . Without loss of generality, we may assume there are some integers ν and µ, such that the many sectors of L bdry are labeled by the independent eigenvalues of U Z k>ν,2l−1 for all internal layers, Ω X k>µ,2l for all l, and Ω Z k,2l 0 .
As before, we define L bdry,0 to be the sector where the U Z , Ω X and Ω Z operators just mentioned all equal to 1. One can show that L bdry,0 is again isomorphic to the fictitious space in Eq. 20 with a similar operator mapping: The n X 2 symmetry projection exists because U X k,1 U X k,L can be generated by the X -suspension operators, the Γ X operators, and the Ω X operators.
Other sectors are all isomorphic to L bdry,0 , and thus to L fic . The discussions for the Ω Z and internal-layer U Z operators turn out to be very similar to the ∆ Z = Z case, and will not be repeated. The eigenvalue of Ω X k>µ,2l can be adjusted with a Z-type operator that may anticommute with some O X i,1 operators, and thus may flip the signs of some O X terms in H I GI while leaving H II GI invariant. We refer the readers to Appendix E.3 for a more explicit description of those isomorphisms.
The change of eigenvalues of Ω X k>µ,2l may conflict with the global conditions U X k ′ ⊗ U X k ′ = 1 in the following sense. Changing the signs of certain O X terms in a GI model is equivalent to altering some X -type symmetry charges, since U X k ′ for all 1 ≤ k ′ ≤ n X is a product of several O X operators. 11 In order to have an anomalous boundary, we expect some of the symmetry charge projections, such as U X k ′ ⊗ U X k ′ = 1 for some k ′ , should hold in all sectors of the boundary theory.
Fortunately, we prove in Appendix E.3 the following result: If for a nonempty subset N ⊂ {1, 2, · · · , n X }, U X N ≡ p∈N U X p can be written as a product of O X operators such that the product is dual to the identity modulo the Γ X operators, then altering the eigenvalue of Ω X k>µ,2l does not affect the corresponding symmetry charge projection U X N ⊗ U X N = 1. A clue for the claim is that U X N ,1 U X N ,L is a product of the X -suspension and Γ X operators, and thus does not depend on the value of Ω X k>µ,2l . The necessary and sufficient condition for an anomalous boundary in the most general situation is given in Claim 1. The first sufficient condition about U X operators is discussed above. The second sufficient condition about Ω Z operators is a direct generalization of the one in the previous subsection and can be proved analogously. When both sufficient conditions are violated, we prove in Appendix E.3 that the boundary theory is nonanomalous. More precisely, the boundary theory in this case is a direct sum of identical sectors. 12 The Hilbert space of each sector is isomorphic to L G ⊗ L G with the operator mapping rule in Eq. 44. The effective boundary Hamiltonian of each sector takes the form of Eq. 21, or more generally consists of local terms generated by those in Eq. 21.
We would like to also remind the readers that each U X (Ω Z ) operator can always be expanded as a product of some O X (∆ Z ) operators, but there may be multiple ways of doing it. The first (second) condition in Claim 1 holds as long as there is one expansion of U X N (Ω Z M ) in terms of the O X (∆ Z ) operators such that the requirement is satisfied.

Examples
We have shown that our basic construction can generate a bulk model with prototypes of topological orders, such as 2 topological order in any dimensions greater than 2 and the X-cube model. Now let us explore further examples. They exploit the capacity of our construction: (1) the construction can produce lattice gauge theories whose gauge group is beyond 2 , (2) the construction can provide a bulk topological order from a GI model describing a symmetry protected topological (SPT) phase, (3) the construction can provide bulk topological orders with only quasi-loop excitations, (4) the same GI model can be matched with more than one 11 Think about flipping the sign of a single transverse-field term in the one-dimensional transverse field Ising model. The sign-flip can be undone via a basis rotation. Nevertheless, the rotation anti-commutes with the 2 symmetry generator. 12 Note that each (new) sector here may be the sum of several (old) sectors discussed above with different values of U X k ⊗ U X k .

22
SciPost Phys. 15, 150 (2023) bulk model with distinct fracton orders. In other words, when fracton order is present, the boundary cannot determine a unique bulk.

Bulk 2 × 2 topological order in two dimensions
The following example shows that the bulk construction can produce bulk topological orders whose gauge group is beyond 2 . The GI model we start with is a one-dimensional model with two global X -type symmetries.
The 2 × 2 symmetry is generated by i X 3i X 3i+1 and i X 3i+1 X 3i+2 . The standard dual model obtained from the generalized Kramers Wannier duality is the following, The bulk model generated through our construction is the following.
We prove in the appendix that this model is topologically ordered. The GSD on a torus is 2 4 . Alternatively, we may obtain the GSD from observing that the Hamiltonian local terms satisfy in total two constrains, In other words, the anyon theory of this stabilizer code has total quantum dimension D = 2 2 . In fact, the anyon theory of the bulk model is the 2 × 2 topological order, equivalent to that of the stack of two toric code models. This is due to a proof showing that the topological phase of any 2d stabilizer code on qubits with translation symmetry is uniquely determined by its total quantum dimension D = 2 n . Its anyon theory is the same as that of n copies of toric code. [85] In our case, n = 2.
The phase diagram of the model (45) is simple. When J ≥ h, the ground state breaks both 2 symmetries spontaneously, and when 0 < J ≤ h, the ground state is the trivial 2 × 2 symmetric state. Comparing (45) and (46) we observe that the model (45) enjoys a self-duality at J = h. In fact, the ground state at J = h is the critical state described by the 4-state Potts model, whose central charge equals to 1. [86][87][88][89] Due to the bulk-boundary correspondence, all these phases appear as well on the boundary of 2 × 2 topological order.

Bulk topological order from a symmetry protected topological model
The one-dimensional spin systems with 2 × 2 symmetry have more phases than those described by the model (45). Particularly, there is a SPT phase which is usually described by the following stabilizer models, In this convention, the symmetry is generated by We ask if we can obtain a solvable model with a topological order adapting the bulk construction to SPT models. Indeed, we find that given a minimal variation to the bulk construction, we obtain a two-dimensional 2 × 2 topological order from a GI model, in whose phase diagram, the SPT is one gapped phase.
To show this, we begin with the GI model. It is the Hamiltonian (50) with additional transverse field terms.
The next step is to obtain the dual model, whose Hamiltonian is the following, Obviously, the dual model is the same as two copies of Ising models in transverse fields. Variant construction Now to construct the bulk model, note that (52) does not satisfy an assumption on the GI models -the local operators {O X i } and {G Z r } = here do not form a CSLO. The price is that the GSD in the bulk model we would obtain from the basic construction is not robust. Part of the degeneracy originates from symmetry breaking orders.
Nevertheless, the violation is modest, and the construction, with a slight variation, can still generate a topological ordered bulk model. Let us give a minimal variation of the basic construction, essentially we modify the rule how we assign the terms across three layers to be centered on odd or even layers, based on the commutation relations of Hamiltonian local operators O Z α 's and O X i 's. The variation is based on the observation that in (52) {X 2i } and {Z 2i+1 } actually form a CSLO. The spirit is that we now assign the three-layer Hamiltonian local terms built with these terms in a CSLO to be centered in the same (odd) layers. We elaborate on the prescription in Appendix G. In particular, we give a sufficient condition when the bulk model from the variant construction has robust ground state subspace. Applying to the current example with SPT order, the modification is enough to provide us a pure topologically ordered bulk model.
The bulk Hamiltonian is This model has a topological order, and its GSD on a torus is 2 4 , when there are in total even number of layers. This can be derived through the standard polynomial representation, as we show in appendix F. As the model is a stabilizer code with translation symmetry, we can conclude that its ground state is the 2 × 2 topologically ordered state. [85].
We show the phase diagram of the model (52) with 2 × 2 symmetry in Fig.5. The phase diagram as well as the critical phases is most easily determined from the dual model (53).

Pure-loop toric code in four dimensions
There is a toric code model in four dimensions that potentially can be used as a thermally stable quantum memory [90,91]. This comes from that the model has no point-like topological excitations, but only loop-like ones. Correspondingly, the effective topological field theory for the model, given by the action S = 2 2π ad b in 5D Euclidean spacetime, where a, b are twoforms, has two 2-form 2 symmetries generated by membrane operators e a and e b . In this subsection, we show how this model comes out of our construction. Consider the GI model on a 3d cubic lattice, described by the following Hamiltonian where qubits live on the links of a 3D cubic lattice with periodic boundary condition, and p, l label plaquettes and links, respectively. The model is homogeneous in the three directions x, y, z. There is obviously no Z-type symmetry. To find the X -type symmetries, it is helpful to note that the four-Z terms in H also appear in the Hamiltonian of the three-dimensional toric code model. We can thus take the local X -type symmetry generators to be star operators of the following form, namely products of six Pauli X operators sharing a single vertex. The model has n X = 3. Take three planes that cut through links and are perpendicular to the x, y, z directions, respectively. The three nonlocal symmetry generators can be taken as the products of Pauli X operators on the links cut by these three planes, respectively. It is illuminating to view the model as having a 1-form symmetry: given each closed membrane that intersect with links, the product of Pauli X operators on the intersecting links commutes with the Hamiltonian. We consider the standard dual theory of the above model. It is convenient to place the dual model on the same cubic lattice, but with qubits living on plaquettes. Each four-Z term in H is dual to a single-Z term associated with the corresponding plaquette. Each single-X term in H is dual to the product of four Pauli X operators on the plaquettes sharing the corresponding link. The dual model is actually equivalent to the original one up to the substitutions J ↔ h and Z ↔ X .
We can now construct the bulk model according to our prescription. To this end, it is convenient to imagine expanding each link in the original model (55) along the 4th spatial direction w, such that the links become plaquettes parallel to the w axis. In this way, the constructed bulk model can be naturally placed on a 4-dimensional cubic lattice with qubits living on plaquettes. We can write the bulk Hamiltonian as where l, p, and c label links, plaquettes (squares), and cubes, respectively. Each term in the first sum is the product of six Pauli-X operators on plaquettes common to a link, and each term in the second sum is the product of six Pauli-Z operators on plaquettes forming a cube.

Same anomalous boundary for two distinct bulks
We have seen that the two-dimensional plaquette Ising model (2) has two different dual models: one is Eq. 6, and the other is Eq. 3 with the substitutions X ↔ Z and J ↔ h. Periodic boundary condition has been assumed for these models in our previous discussion. As a result, two distinct bulk models can be constructed: one with the Hamiltonian Eq. 9 hosts point-like quasiparticle restricted to move along inter-layer direction; and the other, whose Hamiltonian is Eq. 10 with the substitution X ↔ Z, has X-cube fracton order. This example suggests that two different fracton models can share the same boundary theory. Indeed, let us construct one boundary theory that can terminate the above two distinct bulk models. The literal boundary theory with periodic boundary condition along the intralayer directions does not work straightforwardly. This is because the boundary Hilbert space L bdry for the two bulk models have different numbers of sectors. To get around, we take the open boundary conditions along intra-layer directions. In result, L bdry has a single sector for both models. Now we describe the construction. We place the two-dimensional plaquette Ising model on a 2D open square lattice with L x × L y sites, with the same Hamiltonian (2) (only include terms that are completely inside the system). We take the first dual model to be its standard dual, defined on an open square lattice with (L x − 1) × (L y − 1) sites. One can verify that this dual model has no 2 symmetry at all. We define the second dual model on an open square lattice of the same size (L x × L y number of vertices), but now with qubits living on the links. As in the periodic case, we take the dual of each GI term O Z α of the plaquette Ising model to be the product of four Pauli-Z operators around a plaquette. The X -type gauge symmetries of the dual model are generated by the product of all Pauli-X operators around each vertex. The dual of each transverse-field term X i of the plaquette Ising model is uniquely determined modulo the gauge symmetry terms. One can verify that this second dual model has no Z-type symmetry at all. We can then construct two bulk models according to the two dual theories, and take open boundary condition along the out-of-layer direction as well. The boundary Hilbert space L bdry for both models, following the treatment in the Section 5, only has one sector now.
We note that although open boundary condition is taken along the intra-layer directions, our construction actually forbids boundary terms on the four surface planes perpendicular to the layers, according to our result in Appendix E.5. The X -suspension, Z-suspension, and gauge symmetry operators ("bulk" operators) near these surface areas have coefficients as large as those deep inside the bulk, and therefore fix all local degrees of freedom there. With this somewhat artificial setup, the boundary theory only contains local degrees of freedom near the 1st and the L-th layers.
Consequently, we can have exactly the same anomalous boundary theory for these two models: two copies of the two-dimensional plaquette Ising model (possibly with different J and h coefficients), subject to the symmetry projections U X k ⊗U X k = 1 for k = 1, 2, · · · , n X where n X = L x + L y − 1.

Conclusions and Discussions
In this paper, we systematically construct d + 1-dimensional commuting projector models with long-range orders from d-dimensional GI models -qubit lattice models with non-commuting Hamiltonian local terms. The simplicity of the construction allows us to analyze the precise correspondence between the boundary of the long-range ordered state and the GI model. Under a certain condition, the boundary model is subject to global constraints, which are either symmetry charge projections, and/or boundary conditions, implying that the boundary theory to be anomalous. Furthermore, the anomalous boundary model is isomorphic to two copies of the GI models from which the bulk model is constructed, also subject to global constraints. The condition for the anomaly is a surprising requirement on the non-local symmetries in either the GI model or its dual model. This is in contrary to the most common intuition that for all non-local symmetries appearing on the boundary of a long-range ordered state, up to those in the system trivially stacked onto the boundary, only the charge-neutral states under the symmetry contribute to the boundary Hilbert space.
Many open questions following this work worth future exploration. For example, we have not discussed the consequences of bulk anyon fluxes terminating on the boundary theory. Such results will be part of the properties of non-invertible anomaly in higher dimensions. Our construction also has the potential leading to many tangible and interesting generalizations. In particular, qubit stabilizer models only describe a limited class of topological phases [92]. An immediate generalization to n qudit lattice models is worthwhile. With it, we postulate that the bulk models with topological/fracton orders can be constructed staring with a qudit lattice models with any discrete finite Abelian symmetries G in transverse fields, since G can always be written as G = m i=1 n i , with a positive integer m, and integers n i ≥ 2. One question with further stretch is whether the generalization of our models to topological models beyond stabilizer models can generate new topological lattice models, especially in three dimensions or higher. More particularly, the lattice models of topological phases related by Morita equivalence are related by generalized Kramers-Wannier duality. [93][94][95] It is interesting to adapt our alternating layer construction to Morita equivalent topological lattice models and study the topological properties of the resulting bulk model. The Haah's code model [96] does not seem to fit into our construction. Nevertheless, the stabilizers in its Hamiltonian do have a clear trilayer structure, viewed from the (1, 1, 1)-direction. It is interesting to figure out whether some generalized construction works for this representative type-II fracton model. The low energy effective field description for our alternating layer in generating topological theories in one dimensional higher is also in demand. Our construction might remind the readers of the coupled-layer construction in generating topologically ordered models, fracton models, or their hybrids [83,[97][98][99][100][101][102][103][104], in which the model in each layer to begin with is the same, and inter-layer coupling terms are introduced. We leave it for future works to unravel whether there is a relation between our approach and the conventional coupled-layer constructions. Last but not least, our constructed bulk model is straightforward so that their boundary phase diagrams with various Abelian global symmetries can in principle be accessed via numerics. This suggests the possibility of the numerical verification for the conjecture that gapped phases on the boundary of (d + 1)-dimensional topological order with a discrete gauge group G have one-to-one correspondence with the gapped phases on the d-dimensional system with the global G symmetry.
Note added. We came to know of Ref. [105] which also considers correspondence between bulk and boundary of fracton orders. In particular, it also points out that the same boundary theory may be realized for two distinct bulk fracton orders.

A Elementary results in the stabilizer formalism
In this section, we introduce some elementary results in the stabilizer formalism, and set up a few terminologies that will be used in the rest of the appendix. Basics of the stabilizer formalism can also be found in [106,107].
A set of stabilizers is a set of operators satisfying the following properties: • They are tensor products of I, X , Y, Z multiplied by ± sign factors.
• The group generated by them does not contain −1.
The last property guarantees that all the stabilizers are able to simultaneously take the eigenvalue +1. A set of stabilizers is called independent if any product of a nonempty subset of those operators is not proportional to the identity. Specifying the eigenvalue of each independent stabilizer will reduce the Hilbert space dimension by a half. To see this, suppose we have chosen the eigenvalues λ 1 , λ 2 , · · · , λ k for k independent stabilizers U 1 , U 2 , · · · , U k , respectively, with λ i = ±1 for all i. One can see that the (k + 1)-th independent stabilizer U k+1 has zero trace in this eigen-subspace, i.e.
Tr U k+1 due to the zero trace of Pauli operators and the requirement of independence. Thus, U k+1 has equal number of +1 and −1 eigenvalues in this eigen-subspace, and specifying its eigenvalue will further cut down the Hilbert space dimension by a half. As a consequence, it is impossible to find more than N number of independent stabilizers where N is the total number of qubits. A set of stabilizers is called complete if it includes N number of independent stabilizers. Each possible list of eigenvalues of a complete set of stabilizers uniquely determines a state in the Hilbert space.
The following two lemmas about stabilizer extension are useful.

Lemma 1. Any set of independent stabilizers can be extended to a complete and independent one.
Proof. Suppose we have N − k number of independent stabilizers with k > 0. Any set of eigenvalues for those operators determines 2 k number of degenerate eigenstates. It is therefore always possible to find some operator A, not necessarily a tensor product of I, X , Y, Z, that is independent from and commutes with the existing stabilizers. We can expand A as a linear combination of the 4 N number of tensor product operators of I, X , Y, Z which are linearly independent. Let U be any existing stabilizer, UAU −1 = A implies that any tensor product operator entering the expansion of A with a nonzero coefficient must also commute with U. Therefore, we can always find a tensor product of I, X , Y, Z that is independent from and commutes with the existing stabilizers. We can add this new operator to the list of stabilizers and start over again, until the list is complete.

Lemma 2. Suppose we are given a set of independent stabilizers that are all products of Pauli X operators (and the identity operators on the other sites; same below). We can extend the set of stabilizers to a complete one by adding operators that are products of Pauli Z operators.
Proof. Let N be the number of qubits in the Hilbert space and let M be the number of independent stabilizers given. We can represent those stabilizers by an N × M matrix in 2 , the field of integers modulo 2, such that the (i, j) entry of the matrix is 1 when the j-th stabilizer contains an X operator acting on the i-th qubit, and is 0 otherwise. Each column of the matrix corresponds to a stabilizer, and each row corresponds to a qubit in the Hilbert space. Multiplying one stabilizer onto another one will lead to an equivalent set of independent stabilizers, which amounts to adding one column of the matrix onto another one. We are also free to permute the rows of the matrix, which amounts to reordering the qubits. Using these two types of operations, this stabilizer matrix can be cast to the following canonical form One can check that the N operators represented by those two matrices indeed mutually commute and are independent. This completes the proof.

B Hilbert space isomorphism in the generalized Kramers-Wannier duality
In the main text, we introduced the generalized Kramers-Wannier duality by the operator map in Eq. 4. In this section, we will prove that this operator map follows from an isomorphism of symmetric subspaces.
The following result that follows from our definition of the symmetry groups in the main text and Lemma 2 will be useful.

{O Z α } ({∆ X i }) together with all the symmetry generators of the original (dual) theory form a complete set of stabilizers in L (L ′ ).
It is also useful to establish the following result.
, U X n } be a set of stabilizers acting on some multiplequbit Hilbert space L, such that each U Z p (U X q ) is a product of several Pauli Z (Pauli X ) operators acting on different qubits. Let L 0 be the subspace where U Z p = U X q = 1 for all p and q. There exists an orthonormal basis of L 0 such that

Any operator O Z that is a product of several Pauli Z operators and commutes with all the
stabilizers is diagonal in this basis.

Any operator O X that is a product of several Pauli X operators and commutes with all the stabilizers has matrix elements in this basis equal to either 0 or 1.
Proof. Let |{Z i }〉 be the standard eigenbasis for all Pauli Z operators in L, i.e. tensor products of the states (1, 0) T and (0, 1) T . Consider the following list of states, Notice that these states generate all states in L 0 . Moreover, the above states have the following simple properties: (1) |{Z i }〉 = 0 if and only if the spin configuration {Z i } violates the U Z constraints.
(2) Given any two nonzero states |{Z i }〉 and |{Z ′ i }〉, they are equal if |{Z i }〉 can be mapped to |{Z ′ i }〉 by the action of several U X operators, otherwise they are orthogonal. We can restrict the above list of states to a nonzero and nonequal subset. After normalization, this subset of states form an orthogonal basis of L 0 , and it has the desired properties claimed in the statement of the lemma.
The key result for this section is the following.

Theorem 2. The operator map in the generalized Kramers-Wannier duality follows from an isomorphism of Hilbert spaces when we restrict to the symmetric sectors on both sides.
Proof. We start by defining some notations. Let L and L ′ be the full Hilbert spaces of the original theory and the dual theory, respectively. We denote by L Z ⊂ L the symmetric subspace for all Z-type symmetries, and by L 0 ⊂ L the symmetric subspace for both Z-type and X -type symmetries. Thus L 0 ⊂ L Z ⊂ L. Similarly, we define L ′ X and L ′ 0 such that L ′ 0 ⊂ L ′ X ⊂ L ′ . Each possible list of eigenvalues of a complete set of stabilizers uniquely determines a state in the Hilbert space (see Appendix A). Therefore, according to Corollary 5, the eigenvalues of all O Z α uniquely determine a state in L 0 , and the eigenvalues of all ∆ Z α uniquely determine a state in L ′ X . Denote by |{O Z α }〉 ∈ L 0 the simultaneous eigenstates of all O Z α (here we use the same notation for an operator and its eigenvalue). By Lemma 3 just proved above, we can always choose the phase factors of those states such that the matrix elements of each O X i in this basis are either 0 or 1. Similarly, we denote by |{∆ Z α }〉 ∈ L ′ X the simultaneous eigenstates of all ∆ Z α , such that the matrix elements of ∆ X i are either 0 or 1. We define a homomorphism f : As long as f is well-defined which is to be proved below, the operator O X i is mapped to ∆ X i under this map, i.e. f O X i = ∆ X i f . This is because the commuting or anticommuting relations between {O Z α } and {O X i } are the same as those between {∆ Z α } and {∆ X i }. We will prove that f is actually an isomorphism from L 0 to the subspace Since the ∆ Z α operators are not necessarily independent, we need to first confirm that f is well-defined, i.e. the desired image state always exists in L ′ X . Although ∆ Z α may not be independent, they are generated from an independent subset of operators which are able to independently take eigenvalues ±1. Therefore, to prove that f is well-defined, it suffices to show that whenever α∈S ∆ Z α = 1 in L ′ for some set S, α∈S O Z α = 1 is also true in L 0 . This is not hard; since α∈S ∆ Z α commutes with all ∆ X i , α∈S O Z α must commute with all O X i by the duality, and is therefore either the identity or a Z-type symmetry of the original theory which is set to 1 in L 0 .
Next, we would like to show Im f ⊂ L ′ 0 . Corollary 5 says that {∆ Z α } and all the X -type symmetry generators of the dual theory form a complete set of stabilizers in L ′ . It follows that each Z-type symmetry generator of the dual theory is a product of some ∆ Z α operators, which is then dual to either the identity or a Z-type symmetry generator of the original theory. Given that all symmetry generators are set to 1 in L 0 , the image of f must be contained in L ′ 0 . From the definition, we see that f maps a basis of L 0 to a set of linearly independent states in L ′ 0 . It follows that f is injective and dim L 0 = dim Im f ≤ dim L ′ 0 . Now, we can construct another map g : L ′ 0 → L Z similar to f (the roles of X and Z need to be exchanged), and eventually show that dim L ′ 0 ≤ dim L 0 . Therefore, we must have dim L ′ 0 = dim L 0 = dim Im f , which implies that the restricted map f : L 0 → L ′ 0 is one-to-one. This is the isomorphism claimed by the theorem.

C A toy example for sanity checks
In this subsection, we introduce a toy example that will be used later for testing our results. Consider the model and its dual model both defined on a 2D square lattice of the size L x × L y , with periodic boundary condition, and with qubits living on the links. One may shift the dual lattice relative to the original one along y by half the lattice constant, such that the centers of each dual pair of operators coincide. We emphasize that there is no vertical-link single Z (single X ) term in H (H ′ ). The Z-type symmetries of H are the same as those for Eq. 3. Each X -type symmetry generator of H is a product of all vertical-link X operators along even number of vertical lines, which is nonlocal. The dual model H ′ is equivalent to H under the substitutions J ↔ h and Z ↔ X . One may check that n Z = m X = 2 and n X = m Z = L x − 1. A three-dimensional bulk model can thus be constructed according to our prescription.

D Ground state degeneracy of the bulk model
In this section, we compute the GSD of the bulk model with periodic boundary condition along the out-of-layer direction. We denote by L ∈ 2 the total number of layers. We need to introduce some new parameters that enter our final result. Let G Z = U Z 1 ,· · · ,U Z n Z be the group generated by the U Z operators of the original theory. We define two subgroups of it: Obviously, G Z,0 ⊂ G Z,1 . We can always redefine the U Z operators such that, for some integers ν andν, withν ≥ ν G Z,0 = U Z 1 , · · · , U Z ν and G Z,1 = U Z 1 , · · · U Z ν . Similarly, let G ′ X = Ω X 1 , · · · , Ω X m X be the group generated by the Ω X operators of the dual theory. We define the following two subgroups: ρ∈R Γ X ρ for some sets I and R}.
Our main result for this section is the following.

Theorem 3. Let D be the GSD of the bulk model, then
Examples. As the first example, consider the X-cube order described by the Hamiltonian in Eq. 10 that is constructed from Eq. 3 and Eq. 6. We put the model on a 3D cubic lattice with L x × L y × L z number of vertices, with periodic boundary condition along all three directions. The number of layers L is given by L = 2L z . The model has n Z = 2, n X = L x + L y − 2, m X = 0, m Z = L x + L y − 1, andν = ν =μ = µ = 0. It follows that log 2 D = 2(L x + L y + L z ) − 3, which is a well-known result [59].

32
SciPost Phys. 15, 150 (2023) As the second example, consider the model in Appendix C. We take periodic boundary condition along the out-of-layer direction, and the number of layers L has to be an even integer. The model has n Z = m X = 2, n X = m Z = L x − 1,ν =μ = 1, and It follows that which we have verified numerically. The rest of this section is to prove the theorem. Our strategy for counting the GSD is as follows. We extend the set of stabilizers in H bulk to a complete one by including r number of additional stabilizers C 1 , C 2 , · · · , C r which are independent modulo the stabilizers in H bulk , i.e. each C k can not be generated by the other C k ′ operators and the stabilizers in H bulk . Then the GSD is nothing but 2 r . The first useful result is the following.

Lemma 4. The stabilizers in H bulk together with all the nonlocal symmetry operators, namely
and Ω Z k,2l , form a complete set of stabilizers. Proof. As an equivalent statement, any operator A that is a product of Pauli operators and commutes with all these stabilizers can be generated by them (see Corollary 4). As in Theorem 1, we write A = A Z A X . If suffices to show that both A Z and A X are generated by the stabilizers in the theorem statement.
We will focus on A X , and A Z is similar. Since A X commutes with all the G Z and U Z operators, each odd layer of A X is a product of several O X operators. Therefore, using the Xsuspension operators, we are able to remove the X operators in A X on all but one odd layers. Say this remaining odd layer is just the 1st layer. We may thus assume A X only contains X operators on the 1st layer or with even layer indices. To commute with the Z-suspension terms, different even layers of A X must differ from each other only by some Γ X and Ω X operators. Thus, using the Γ X and Ω X operators, A X can be cast to the form A 1 X p∈Λ L/2 l=1 X p,2l for some set Λ, where A 1 X is an operator that only acts on the 1st layer. This A X has to commute with the Ω Z k,2l and Γ Z σ,2l operators, which are the Z-type symmetry generators for the dual theory. It follows from our Lemma 2 that p∈Λ X p must be a product of ∆ X i , Γ X ρ and Ω X k . Therefore, using the X -suspension operators, the Γ X and the Ω X operators, we are able to reduce A X to an operator that acts on the 1st layer. In order to commute with the G Z , U Z and Z-suspension operators, this single-layer operator is an X -type symmetry generator of the original theory, and thus a product of the G X and U X operators. This completes the proof.
Next, we shall restrict the nonlocal symmetry operators to a subset such that, modulo the stabilizers in H bulk , the subset is independent and equivalent to the original set. First notice that, starting from U X k,2l 0 −1 for some l 0 , using the X -suspension, Γ X and Ω X operators, we are able to generate U X k,2l−1 for all l. This is because each U X k operator acting on the original lattice can be written as a product of several O X operators, and is then dual to an X -type symmetry generator of the dual theory. It thus suffices to retain U X k,2l 0 −1 while dropping the U X operators acting on all the other layers. Similarly, it suffices to retain Ω Z k,2l 0 while dropping the Ω Z operators acting on all the other layers. Secondly, according to the definition of G Z,0 , each U Z k,2l−1 with k ≤ ν is a product of some Z-suspension and G Z operators, and thus can be dropped. Similar for the Ω X k,2l operators with k ≤ µ. Thirdly, according to the definition of G Z,1 , L/2 l=1 U Z k,2l−1 for ν < k ≤ν can be generated by the Z-suspension and G Z operators, thus we can drop the U Z k,1 operators acting on the 1st layer with ν < k ≤ν. Similarly, we can drop the Ω X k,2 operators acting on the 2nd layer with µ < k ≤μ. After all these steps, we are left with the following reduced set of nonlocal symmetry operators: • U X k,2l 0 −1 for all k and some l 0 ; Ω Z k,2l 0 for all k and some l 0 .
• U Z k,2l−1 forν < k ≤ n Z and all l; Ω X k,2l forμ < k ≤ m X and all l.
We complete the proof for Theorem 3 by the following claim.
Lemma 5. The above stabilizers are independent modulo the stabilizers in H bulk . Equivalently, they can independently take the eigenvalues ±1 in the ground state subspace.
Proof. It suffices to find operators that can independently flip the signs of the above stabilizers without affecting the stabilizers in H bulk . In the main text, we have discussed how to flip the signs of U X k,2l 0 −1 or Ω Z k,2l 0 , and let us not repeat it here. According to the definition of G Z,1 andν, each U Z k>ν is independent from (not a product of) the O Z and G Z operators. As a result, there exist some X -type operators V X k>ν acting on the Hilbert space of the original model such that each V X k>ν anticommutes with U Z k>ν but commutes with all the other U Z k ′ >ν , O Z and G Z operators. We can then use the operator V X k>ν,2l−1 to flip the sign of U Z k>ν,2l−1 , without affecting the other stabilizers listed above and the stabilizers in H bulk . The scenario for the Ω X k>μ,2l operators is similar. Flipping the sign of a U Z k,2l−1 operator with ν < k ≤ν and l ̸ = 1 is somewhat complicated, consisting of several steps: (1) Let V X k≤ν be X -type operators acting on the Hilbert space of the original theory such that V X k for each k ≤ν anticommutes with U Z k but commutes with all the other local or nonlocal Z-type symmetry generators of the original theory (G Z r and U Z k ′ for all the other k ′ ). Apply V X k,2l−1 with some k ∈ {ν + 1, ν + 2, · · · ,ν} to flip the sign of U Z k,2l−1 . This will necessarily flip the signs of some Z-suspension operators centering on the (2l − 1)-th layer as well, which we need to fix. (2) We can restrict the ∆ Z operators of the dual theory to an independent and equivalent subset: , · · · . It suffices to fix the Z-suspension operators centering on the (2l − 1)-th layer with these α-indices (α 1 , α 2 , · · · ) because they can generate all the other Z-suspension operators centering on the (2l − 1)-th layer with the help of the G Z r,2l−1 and U Z k≤ν,2l−1 operators. There exists a set of operators E X while commutes with all the other independent ∆ Z operators. We can use ) when l 0 ≥ l (l 0 < l) to fix the eigenvalues of the Z-suspension operators centered on the (2l − 1)-th layer without affecting Ω Z k,2l 0 . This step will flip the signs of some Z-suspension operators centering on the 1st layer. (3) Apply V X k,1 with the same k as in the first step. Notice that the Γ Z operators will not be affected by the above steps, because Γ Z σ,2l 0 are unaffected, and the other Γ Z operators can all be generated from the Γ Z σ,2l 0 , Z-suspension, and G Z operators.

E Further details on the bulk-boundary correspondence
This section is devoted to understanding the boundary theory of our bulk model: the boundary Hilbert space L bdry and the boundary Hamiltonian H bdry which are defined in Section 5. The bulk model is defined in d + 1 spatial dimensions, and we expect the boundary theory to have an effective d-dimensional description that will enable us to examine whether it is anomalous or not.
In Appendix E.1, we derive a complete and independent set of stabilizers characterizing the full Hilbert space of our system, which turns out to be useful for understanding L bdry . We work out a d-dimensional effective description of the boundary theory as well as a condition for anomaly in Appendix E.2 for a special case and in Appendix E.3 for the general situation. The analysis uses a classification result about possible local terms in H bdry , which is proved in Appendix E.5. In Appendix E.4, we discuss in what cases the anomaly condition is violated.

E.1 Stabilizers characterizing the full Hilbert space
Our strategy for analyzing L bdry is as follows. We first find a complete and independent set of stabilizers that characterizes the full Hilbert space, done in this subsection. Then we divide this set into two subsets, such that one subset of stabilizers is equivalent to the stabilizers appearing in the bulk Hamiltonian. It follows that the boundary Hilbert space is characterized by the other subset of stabilizers. Once this is done, we will find a natural d-dimensional description of L bdry .
The first useful result is the following.

Lemma 6.
The following operators form a complete (but not independent) set of stabilizers.
1. All the O X operators on the two boundary layers.

All the X -suspension and Z-suspension terms in the bulk Hamiltonian.
3. All G Z r,2l−1 and U Z k,2l−1 , namely the Z-type symmetry generators of the original theory acting on all odd layers.

All Γ X ρ,2l
and Ω X k,2l , namely the X -type symmetry generators of the dual theory acting on all even layers.
and Ω Z k,2l 0 for all σ, k and some fixed l 0 .
Proof. Our strategy is similar to that for Theorem 1 in the main text. Suppose A = A Z A X commutes with all the above operators, where A Z (A X ) is a product of Pauli Z (Pauli X ) operators. We will show that A is a product of the above operators. Once this is proved, it follows from Corollary 4 that those stabilizers are complete.
Let us start with A Z . Similar to the case of Theorem 1, because of the existence of the X -suspension operators and the O X operators on the boundary layers, we can reduce A Z to an operator that acts on a single layer with an even layer index, by repeatedly multiplying it with the Z-suspension operators, the G Z operators and the U Z operators. This single-layer operator must be a symmetry (local or nonlocal) of the dual theory, namely a product of the Ω Z k,2l and Γ Z σ,2l operators for some l. Note that using the Z-suspension operators, the G Z operators and the U Z operators, we are able to generate Ω Z k,2l and Γ Z σ,2l for all l starting from Ω Z k,2l 0 and Γ Z σ,2l 0 . We have thus proved that A Z is a product of the operators in the statement of the lemma.
Next, we consider A X . Since A X commutes with all the G Z and U Z operators, each odd layer of A X is a product of several O X operators. Therefore, using the O X operators on the boundary layers and the X -suspension operators, we are able to remove all the odd-layer X operators in A X . We may thus assume A X only contains X operators with even layer indices. Then, to commute with all the Z-suspension terms, different even layers of A X must differ from each other only by some Γ X and Ω X operators. Thus, using the Γ X and Ω X operators, A X can be cast to the form p∈Λ (L−1)/2 l=1 X p,2l for some set Λ. This A X has to commute with the Ω Z k,2l 0 and Γ Z σ,2l 0 operators, which are the Z-type symmetry generators for the dual theory.
It follows from our Corollary 5 that p∈Λ X p must be a product of ∆ X i , Γ X ρ and Ω X k . One can

35
SciPost Phys. 15, 150 (2023) then see that A X is a product of the X -suspension operators, the boundary O X operators, the Γ X and the Ω X operators. This completes the proof.
The next step is to restrict the set of stabilizers in the above lemma to an independent (and still complete) subset. This is helpful because only independent stabilizers can independently take eigenvalues ±1. Let us start by defining some notations. The local Z-type symmetry generators G Z r of the original theory may not be independent, and we can restrict it to an independent subset withñ Z number of elements. This process is equivalent to choosing a basis for a vector space over 2 . Similarly, we denote byñ X the number of independent local X -type symmetry generators of the original theory, and denote bym X (m Z ) the number of independent local X -type (Z-type) symmetry generators of the dual theory. Recall that we denote the full Hilbert spaces of the original and the dual GI models by L and L ′ , respectively.
First consider the O X operators on layer 1, one of the two boundary layers. As mentioned above, only N −n Z −ñ Z number of them are independent. For our convenience later, we replace these N − n Z −ñ Z number of independent O X i,1 operators by the following equivalent set of stabilizers: n X number of U X k,1 ,ñ X number of G X s,1 , and additional N −n Z −ñ Z −n X −ñ X number of O X i,1 operators. The same applies to the O X operators on layer L, the other boundary layer. It turns out that U X k,1 and U X k,L are not independent from each other. Using the X -suspension operators, and the Γ X , Ω X operators acting on all even layers, one can generate the operators U X k,1 U X k,L , because each X -type symmetry generator of the original theory can be written as a product of a few O X operators, and the product is dual to the identity or an X -type symmetry generator in the dual theory. Similarly, G X s,1 G X s,L can also be generated. Therefore, we can remove U X k,L and G X s,L as redundancies from our list of stabilizers. The X -suspension operators O X ∆ X O X are in general not independent. In particular, given any relation between the O X i operators, say i∈Λ O X i = 1 for some set Λ, the operator is an X -type symmetry generator acting on the lattice layer 2l, and is therefore a product of the Γ X and Ω X operators. This means that for each l, it suffices to retain N − n Z −ñ Z number of the X -suspension operators Things become even simpler when the standard dual theory is used, i.e. when ∆ Z α = Z α . In this situation, due to the absence of Γ X and Ω X operators, each relation between the O X i operators directly translates to a relation between the X -suspension operators, , which will be useful later. Moreover, all the Z-suspension operators are independent in this special case.
With all the above procedures of removing redundancies, let us count how many stabilizers we are left with. We have 1. n X number of U X k,1 ,ñ X number of G X s,1 , and additional 2(N − n Z −ñ Z − n X −ñ X ) number of O X operators acting on the two boundary layers. Theorem 2 about the generalized Kramers-Wannier duality implies that With this relation, one can check that the total number of stabilizers listed above is N ′ (L − 1)/2 + N (L + 1)/2, same as the total number of qubits. Therefore, the above stabilizers must be independent. A complete basis of states in the Hilbert space are labeled by the independent eigenvalues of those operators.

E.2 Simple situation: using the standard dual theory
In this subsection, we restrict to the special situation where the standard dual theory is used, i.e. ∆ Z α = Z α . This means that the dual theory has no X -type symmetry at all. All the Γ X and Ω X operators disappear, and m X =m X = 0. Now, we divide the complete and independent set of stabilizers found above into two subsets: • Subset A: the Z-suspension and X -suspension operators, G Z r,2l−1 for all l, G X s,1 acting on the first layer, and Γ Z σ,2l 0 acting on the layer 2l 0 .
• Subset B: the 2(N − n Z −ñ Z − n X −ñ X ) number of O X operators acting on the two boundary layers, the n X number of U X k,1 , the m Z number of Ω Z k,2l 0 , and the n Z (L + 1)/2 Subset A is actually equivalent to the set of stabilizers appearing in the bulk Hamiltonian: Firstly, recall that Subset A already contains all the Z-suspension operators, and the Xsuspension operators contained in Subset A are able to generate all the other X -suspension operators, thanks to the absence of Ω X and Γ X . Moreover, it is obvious that all the G Z operators can be generated by the independent ones we retained here. Finally, all the G X operators can be generated using the independent set of G X s,1 acting on the first layer and the X -suspension operators. Similarly, all the Γ Z operators are generated by the independent set of Γ Z σ,2l 0 acting on the layer 2l 0 , the Z-suspension operators and the G Z operators. This result implies that L bdry , the ground subspace of the bulk Hamiltonian, is completely characterized by the stabilizers in Subset B. Example. Consider the bulk model in Eq. 10 that is constructed from Eq. 3 and Eq. 6. We put the model on a 3D cubic lattice with L x × L y × L z number of vertices, with periodic boundary condition along x and y, and open boundary condition along z. The two boundary surfaces are "smooth", and L is related to L z by L = 2L z −1 (the height of the system is L z −1 number of lattice constants). We have n Z = 2,ñ Z = L x L y −1, n X = L x + L y −2,ñ X = 0, m Z = L x + L y −1, N = 2L x L y , and N ′ = L x L y . It follows that log 2 (dimL bdry ) = 2L x L y + L which we have verified numerically.
It is convenient to divide L bdry into several sectors labeled by the eigenvalues of Ω Z k,2l 0 and U Z k,2l−1 . Each of these sectors is characterized by a set of stabilizers that are all products of Pauli X operators. Let us first focus on the sector where Ω Z k,2l 0 = 1 and U Z k,2l−1 = 1, denoted byL bdry,0 (the reason for a tilde symbol will be clear below). We denote by L the Hilbert space of the original d-dimensional lattice, and by L G the gauge-invariant subspace of it (the symmetric sector for all local symmetries). We see thatL bdry,0 is isomorphic to the following fictitious space,L fic := where the 2 symmetries are generated by U X k ⊗ U X k , U Z k ⊗ 1 and 1 ⊗ U Z k . We would like to choose the isomorphism fromL bdry,0 to the above fictitious space such that the independent stabilizers characterizingL bdry,0 have the most natural correspondence, i.e.
Thus, we shall identity the simultaneous eigenstates of the corresponding stabilizers. One simple consequence is that, all the O X operators on the two boundary layers, not just the independent ones in Subset B, will satisfy the above mapping rule: However, we have not uniquely determined the isomorphism yet, since the eigenstates have arbitrary phase factors. Thanks to the fact that all the stabilizers are products of Pauli Xoperators, we can choose orthonormal bases for bothL bdry,0 and the fictitious space according to Lemma 3, with the roles of X and Z exchanged. Vectors in these orthonormal bases are automatically eigenstates of the stabilizers and can be identified according to the eigenvalues. As one important consequence of this special choice of bases, the operators O Z α,1 Z α,2 and Z α,L−1 O Z α,L , which all preserveL bdry,0 and can be added to the boundary Hamiltonian, have the simple mapping rule: We may take the boundary Hamiltonian H bdry to be a linear combination of the operators From the above discussion, we see that theL bdry,0 block of H bdry , when represented in the fictitious spaceL fic , takes the form where H I GI and H II GI act on the two copies of L G in Eq. E.1, respectively. We prove in Appendix E.5 that any allowed local term in H bdry can be generated by the four types of terms considered here and the stabilizers in H bulk . Therefore, our canonical choice of H bdry is a quite general one.
Other sectors with different eigenvalues of Ω Z k,2l 0 or U Z k,2l−1 can be analyzed by looking for unitary operators that map them toL bdry,0 . Given a set of independent Z-type operators A Z 1 , A Z 2 , · · · , A Z n acting on an arbitrary multiple-qubit Hilbert space, there is always a set of Xtype operators B X 1 , · · · , B X n such that B X k anticommutes with A Z k but commutes with all the other A Z operators. It follows that we can always find some X -type operators which commute with the bulk Hamiltonian but can independently change the eigenvalues of Ω Z k,2l 0 and U Z k,2l−1 . They generate isomorphisms from all the other sectors of L bdry toL bdry,0 which is equivalent toL fic as we just shown. These isomorphisms commute with O X i, 1 and O X i,L , but may anticommute with some O Z α,1 Z α,2 and Z α,L−1 O Z α,L operators. Therefore, the effective boundary Hamiltonian in each of the other sectors takes the same form as Eq. E.4, but the signs of some GI terms may be flipped.
Let us try to write down the isomorphisms between different sectors more explicitly, which turn out to be illuminating. We start from the U Z operators acting on the two boundary layers, namely U Z k,1 and U Z k,L . Let V X 1 , V X 2 , · · · , V X n Z be the X -type operators acting on L such that V X k anticommutes with U Z k but commutes with all the other Z-type symmetry generators (local or nonlocal) of the original theory. V X k,1 and V X k,L commute with the bulk Hamiltonian, and thus can be used to adjust the eigenvalues of U Z k,1 and U Z k,L within L bdry , respectively. Each V X k,1 may flip the signs of some O Z terms in the H I GI part of the effective boundary Hamiltonian. We can, however, apply the unitary operator V X k ⊗ 1 toL fic to compensate these sign changes, at the cost of flipping the sign of U Z k ⊗ 1. 13 A similar statement holds for V X k,L . This observation inspires us that there is actually a nicer way of viewing L bdry : We can alternatively divide L bdry into sectors labeled by Ω Z k,2l 0 and U Z k,2l−1 for all internal layers (3 ≤ 2l − 1 ≤ L − 2). Each new 13 V X k ⊗ 1 is not an automorphism ofL fic , so strictly speaking we shall first embedL fic to L G ⊗ L G , and then apply V X k ⊗ 1.
sector now includes several old sectors that differ from each other only by the eigenvalues of U Z k,1 and/or U Z k,L . Let L bdry,0 be the new sector where Ω Z k,2l 0 = 1 and U Z k,2l−1 = 1 for internal layers. L bdry,0 is isomorphic to a new fictitious space where the 2 symmetries are generated by U X k ⊗U X k . The operator mapping rule fromL bdry,0 tõ L fic that we established earlier still works here, now from L bdry to L fic . The effective boundary Hamiltonian takes the form of Eq. E.4 as well. Other new sectors are all isomorphic to L bdry,0 , and thus to L fic : • Denote by L ′ the Hilbert space for the d-dimensional dual lattice. We may find some Xtype operators Θ X k acting on L ′ such that Θ X k anticommutes with Ω Z k but commutes with all the other Z-type symmetry generators (local or nonlocal). It follows that is a multilayer operator that commutes with the bulk Hamiltonian and can flip the eigenvalue of Ω Z k,2l 0 . This multilayer operator may flip the signs of some O Z terms simultaneously in both H I GI and H II GI .
• Adjusting the values of the internal-layer U Z operators is more complicated, consisting of three steps: (1) Flip the sign of U Z k,2l−1 using V X k,2l−1 defined earlier, but this may unexpectedly flip the signs of some Z-suspension operators as well. (2) Fix the Z-suspension operators using string operators of the form operators. Notice that the Γ Z operators will not be affected by the above steps either, because Γ Z σ,2l 0 are unaffected, and the other Γ Z operators can all be generated from the Γ Z σ,2l 0 , Z-suspension, and G Z operators. (3) As an optional step, apply V X k,1 (V X k,L ) when l 0 ≥ l (l 0 < l). With the last step added, the effective boundary Hamiltonian will be invariant under the above operations.
We have seen that altering the eigenvalue of Ω Z k,2l 0 may change the signs of some O Z terms in both the H I GI and H II GI parts of the effective boundary Hamiltonian. Is it possible to cancel this effect by some additional unitary rotation on L fic ? The answer depends on a certain property of the Ω Z operators. Let G ′ Z = Ω Z 1 , · · · , Ω Z m Z be the group generated by Ω Z 's in the dual theory.
We define a subgroup G ′ Z,0 ={g ∈ G ′ Z |g is dual to the identity modulo the G Z operators}. We can always redefine the Ω Z operators such that for some integer m 0 , G ′ Z,0 = Ω Z 1 , · · · , Ω Z m 0 . We claim and will elaborate below that: When k > m 0 (k ≤ m 0 ), it is possible (not possible) to flip the sign of Ω Z k,2l 0 without affecting the effective boundary Hamiltonian.
First consider k ≤ m 0 . We write Ω Z k = α∈A Z α for some subset A, then under the Kramers-Wannier operator map, Ω Z k → α∈A O Z α = r∈R G Z r for some subset R. In an arbitrary sector of L bdry , where we used the fact that Ω Z k≤m 0 ,2l 0 is related to Ω Z k≤m 0 ,2 by the multiplication of several Z-suspension and G Z operators. In L fic , we have Suppose there is an isomorphism from this sector of L bdry to L fic , such that , · · · that is independent modulo is not a product of the remaining ones and G Z 's. Then there exist while commutes with all the others in the independent subset and also commutes with all G Z 's. Now, using the operators Q X β j ⊗ 1 acting on L fic , we can freely adjust the signs of the O Z β j ⊗ 1 terms in the effective boundary Hamiltonian. In other words, we can freely adjust the signs of η β j whose definition is in Eq. E.8 above. In fact, once we correct all the sign changes of η β j due to the sign flip of Ω Z k,2l 0 , the sign changes of all η α are also corrected. To see this, notice that given any relation of the form α∈A O Z α = 1 mod G Z , α∈A Z α must be an element of G ′ Z,0 modulo some Γ Z 's, which follows from the definition of G ′ Z,0 as well as the fact that each Γ Z operator is dual to the product of some G Z 's, and thus α∈A η α is equal to the eigenvalue of this element of G ′ Z,0 acting on the 2l 0 -th layer. Since each O Z α is a product of some O Z β j and G Z operators, each η α is then a product of some η β j and the eigenvalue of some element of G ′ Z,0 acting on the 2l 0 -th layer. Flipping the sign of Ω Z k>m 0 ,2l 0 does not affect any Ω Z k ′ ,2l 0 with k ′ ≤ m 0 , therefore our statement above about η α is indeed true.
With all these discussions, we propose the following necessary and sufficient condition for anomaly.
Claim. When the standard dual theory is used for constructing the bulk model, the boundary theory is anomalous if and only if either of the following two conditions is satisfied: 2. For some nonempty subset K ⊂ {1, 2, · · · , m Z }, k∈K Ω Z k is dual to the identity modulo the G Z operators.
The first condition guarantees a symmetry charge projection U X k ⊗ U X k = 1 acting on the two copies of the GI model in the boundary theory; each copy is allowed to have states with both values of the symmetry charge, but the total charge of the two copies is fixed. When the second condition holds, k∈K Ω Z k is dual to a generalized boundary condition, and a boundary condition projection is applied to the two copies of the GI model in the boundary theory; each copy is allowed to take both values of the boundary condition, but the boundary condition values of the two copies are locked together. When neither condition is satisfied, the boundary theory is a direct sum of identical sectors. The Hilbert space of each sector is isomorphic to L G ⊗ L G with the operator mapping rule in Eq. E.2 and E.3. The effective boundary Hamiltonian of each sector takes the form of Eq. E.4, or more generally consists of local terms generated by those in Eq. E.4.

E.3 General situation
To work with the most general situation, let us recall some notations defined in Appendix D. In general, the nonlocal Z-type symmetry group G Z ∼ = n Z 2 of our GI model contains a subgroup G Z,0 ∼ = ν 2 , whose definition is G Z,0 ={g ∈ G Z | g = α∈A O Z α r∈R G Z r for some sets A and R satisfying α∈A ∆ Z α = 1}. We can always choose our symmetry generators such that G Z,0 is generated by U Z 1 , · · · , U Z ν . Similarly, the nonlocal X -type symmetry group G ′ X ∼ = m X 2 of the dual theory contains a µ 2 subgroup, whose definition is SciPost Phys. 15, 150 (2023) some sets I and R satisfying i∈I O X i = 1}. Without loss of generality, we assume this µ 2 subgroup is generated by Ω X 1 , · · · , Ω X µ . As in the previous case, we divide the complete and independent set of stabilizers for the full Hilbert space into two subsets: • Subset A: the Z-suspension and X -suspension operators, G Z r,2l−1 for all l, G X s,1 acting on the first layer, Γ X ρ,2l for all l, Γ Z σ,2l 0 acting on the layer 2l 0 , U Z 1,2l−1 , · · · , U Z ν,2l−1 for all the internal layers (3 ≤ 2l − 1 ≤ L − 2), and Ω X 1,2l , · · · , Ω X µ,2l for all l.
• Subset B: the 2(N − n Z −ñ Z − n X −ñ X ) number of O X operators acting on the two boundary layers, the n X number of U We take open boundary condition along the out-of-layer direction, and the number of layers L is an odd integer as explained in Section 5. One may check that n Z = m X = 2,ñ Z = L x L y − 1, n X = m Z = L x − 1,ñ X = 0, and It follows that log 2 (dimL bdry ) = 2L x L y + L (L x odd) , 2L x L y + 2L − 2 (L x even) , (E.10) which we have verified numerically. We can divide L bdry into several sectors labeled by the eigenvalues of Ω Z k,2l 0 , U Z k>ν,2l−1 for all internal layers, and Ω X k>µ,2l . Denote by L bdry,0 the sector where these operators all equal to 1. Then L bdry,0 is again isomorphic to the fictitious space in Eq. E.5. The operator mapping rule is also similar: (E.11) By taking H bdry as a linear combination of 2 and ∆ Z α,L−1 O Z α,L , the L bdry,0 block of H bdry may again have the form of Eq. E.4 when represented in L fic . We prove in Appendix E.5 that any allowed local term in H bdry can be generated by the four types of terms considered here and the stabilizers in H bulk .
The next task is to establish isomorphisms from other sectors to L bdry,0 , and therefore to the fictitious space: • As before, the eigenvalue of each Ω Z k,2l 0 can be flipped by the operator without affecting the U Z or Ω X operators. This multilayer operator may flip the signs of some O Z terms simultaneously in both H I GI and H II GI .

E.5 On possible boundary terms
As we defined in the main text, an allowed boundary term is a local operator that commutes with the stabilizers in the bulk Hamiltonian in Section 5 (open boundary condition along the out-of-layer direction). The boundary terms that we have considered so far are This set turns out to be complete in the following sense.

Lemma 7.
Given any local operator A that is a product of Pauli operators and commutes with all the stabilizers appearing in the bulk Hamiltonian in Section 5, A can be locally generated by those bulk stabilizers together with Proof. The proof is not much different from that for Theorem 1. We again write A = A Z A X . Both A Z and A X must commute with all the stabilizers in the bulk Hamiltonian, and we will prove that both of them satisfy the claim of the lemma. Let us start from A Z . Denote by l max and l min the maximal and minimal layer indices of A Z 's support, respectively. Since A Z is supposed to be local, i.e. small compared to the system size, either l min ≫ 1 or l max ≪ L. These two scenarios are nearly identical, so let us just assume l max ≪ L. We can then use the reduction procedure described in the proof of Theorem 1 to reduce l max until l max ≤ 2, if A Z is not yet fully reduced to the identity. Now suppose l max = 2, in order to commute with the Γ X operators on the 2nd layer, the top layer of A Z must be a local product of some ∆ Z operators. Thus, we can further reduce A Z using the O Z α,1 ∆ Z α,2 operators so that it no longer has any support on the 2nd layer. If l max = 1, in order to commute with the X -suspension operators spanning the 1st, 2nd and 3rd layers, this single-layer A Z must be a local product of some G Z operators. This completes the proof for A Z .
Next, we consider A X . We similarly define l max and l min , and assume l max ≪ L without loss of generality. Using the reduction procedure in the proof of Theorem 1, we can reduce l max all the way to 1, if A X is not yet fully reduced to the identity operator. If l max = 1, in order to commute with the G Z operators, this single-layer A X must be a local product of the O X i,1 operators. We have thus completed the proof for A X .
Some minor technical comments: (1) The above result still holds without assuming A to be a product of Pauli operators. We can expand A as a superposition of the linearly independent tensor product operators of I, X , Y, Z. Since Pauli operators either commute or anticommute, SAS −1 = A for a bulk stabilizer S implies that any tensor product operator entering the expansion of A with a nonzero coefficient must also commute with S. (2) If A is a local operator that preserves L bdry but does not commute with the bulk Hamiltonian, its effect on L bdry is equivalent to some local operator that commutes with the bulk stabilizers, so it is unnecessary to consider such an operators as a boundary term candidate. To see this, notice that this operator A overlaps with at most finite number of bulk stabilizers since it is local. We can then repeatedly "symmetrize" the operator by A → (A + SAS −1 )/2 for each overlapping bulk stabilizer S. The resulting new operator has the same action on L bdry , commutes with all the bulk stabilizers, and has a linear size exceeding the old one by at most an O(1) amount.

F The polynomial representation and topological orders
We derive the properties of the bulk model in (54) through the polynomial representation introduced in [58], see also [65] for pedagogical purpose. The Hamiltonian local terms are represented by a stabilizer map, which is a matrix with elements in 2 [x, x −1 , y, y −1 ], the Laurent polynomials whose coefficients are in 2 , wherex ≡ x −1 andȳ ≡ y −1 . The model has a robust GSD. This is determined from that ker ε S = im S, where ε S is the excitation map. [65] The Hamiltonian in the stabilizer formalism with such a property is also called an exact code. We compute the GSD on a square lattice with L x × L y sites and periodic boundary conditions. Treating theŷ as the layer-indexed direction, we also assume L y to be even. Since the number of types of stabilizers in the Hamiltonian is the same as the number of qubits in a unit cell, the GSD can be computed from The block-diagonal structure in S † allows us to reduce the evaluation further [108,109] to where S Z and S X are sub-matrices of the block matrix S = [S Z , 0; 0, S X ], I(σ) is the ideal of σ, and b L is the ideal generated by the polynomials that declare boundary conditions. We find that for the periodic boundary conditions represented by b L = 〈x L x − 1, y L y − 1〉 with even L y , the associated ideals represented in a Groebner basis [65] are The degeneracy is the same as two copies of toric code. Similarly, we can also obtain the GSD for the 2 × 2 model given in (48) from the polynomial representation. The stabilizer map for this model is The model is locally topological ordered, implying the absence of symmetry breaking order. This is determined by that S 2 × 2 satisfy the codimension condition, 15 codim I(S C ) ≥ 2, where S C is the stabilizer map of the classical spin model, which becomes the quantum model after gauging, and I(S C ) is the ideal of S C . In the current case, S C is the map (1 + x + x −1 , y + y −1 ), and codim I(S C ) = 2.

G The variant bulk construction
We explain the variation from the basic bulk construction, which allows us to produce the pure topologically ordered model from the SPT model in (52).
Towards constructing the bulk model, we note that (52) does not satisfy the assumptions in Section 2. The price is that the GSD in the bulk model we would obtain from the basic construction is not robust. Part of the degeneracy originates from symmetry breaking orders. Nevertheless, let us give a minimal variation of the construction. This is enough to provide us a pure topologically ordered bulk.
Firstly, we note that the GI model we begin with has the following property. We group the operators {O Z α } and {O X i } according to their commutation relations.
The two sets have the property that all operators in each set commute with each other; for any Z-type (X -type) operator in one set, there are X -type (Z-type) of operators in the other set that anti-commutes with it. By design, O Z 0 ∈ A 1 . Now we give the modified rule in constructing the bulk model. For local operators in A 1 , their corresponding three-layer local operators are ∆ Z α,l−1 O Z α,l ∆ Z α,l+1 and ∆ X i,l−1 O X i,l ∆ X i,l+1 in the bulk model and center at odd layers i.e.,l ∈ 2 + 1. On the other hand, for local operators in A 2 , their corresponding three-layer local opertors O Z α,l−1 ∆ Z α,l O Z α,l+1 and O X i,l−1 ∆ X i,l O X i,l+1 in the bulk model center at even layers, i.e.,l ∈ 2 . This rule of determining the layer index of local terms supported on three-layers is the only modification in the variant construction.
Explicitly, the bulk Hamiltonian coming from the variant construction is the following, In the example of the GI model capturing the 2 × 2 SPT phase, the two sets of local operators are (G.2) Theorem 5. A sufficient condition for the bulk model to have a robust ground state subspace is that the GI model has the following properties: • A 2 forms a CSLO.
• The dual of A 1 , denoted as A ′ 1 , forms a CSLO.
• Neither the GI model nor the dual model has local symmetries.
We can see that the 2 × 2 model has these properties. Particularly, in this example, the dual of A 1 is A ′ 1 = {Z 2i+1 , X 2i } is a CSLO. Now let us prove the above theorem. To prove that the bulk model has a robust ground state subspace is the same as to prove that the terms in the bulk Hamiltonian form a CSLO.
Suppose there is a local operator that commutes with all terms in the Hamiltonian of the bulk model, we will show it must be either an identity operator or a product of Hamiltonian local terms.
The operator is local in the sense that its support on any layer is finite, independent of total system size along any direction. We begin with considering that the support of the local operator has a single connected component.
First, it cannot be a local operator that is non-trivial only within a single layer. Such an operator, if it existed on an odd layer, would commute with all operators in A 1 and A 2 , thus it would be a local symmetry operator of H GI . However, there is no local symmetry in the GI model, as required in the properties. Similarly, if the operator is within an even layer, it would be a local symmetry of the dual model, and this violates the required properties. In the end, a local symmetry operator within a single layer for the bulk model does not exist.
Second, it cannot be a local operator within two adjacent layers. Suppose there exists such an operator, and let us call it A. Without losing generality, let us suppose its support is on the z-th layer (an odd layer) and the z + 1-th layer (an even layer). A can thus be decomposed as A o A e , with A o (A e ) on the odd (even) layer. A commutes with all terms in the Hamiltonian of the bulk model. In particular, A commutes with all terms O k,z−2 ∆ k,z−1 O k,z , for any operator O k in A 2 . A at most share supports with these operators on the z-th layer. This means A o commutes with A 2 . Through similar steps, one can show that A e commutes with A ′ 1 , the dual of A 1 . Next, A also commutes with O k,z ∆ k,z+1 O k,z+2 ,for any operator O k in A 2 . Because A o on the z-th layer commutes with all O k ∈ A 2 , A e on the z + 1-th layer must commute with all ∆ k ∈ A ′ 2 . Thus, A e commutes with A ′ 1 ∪ A ′ 2 , which is the set of all local operators in the dual of the GI model. Since we have required that the dual model has no local symmetries, A e is an identity operator. Through a similar step, we can show that A o must commute with A 1 . Combined with the derivation several steps above, this means that A o commutes with all Hamiltonian local terms in the GI model. And as we require the GI model has no local symmetries, A o is at most an identity operator. In conclusion, A = A o A e if commutes with all Hamiltonian local terms of the bulk model, must be an identity operator.
As the final case, we consider that the local operator A has a support from the z min -th layer to the z max -th layer. Both z min and z max are finite, and z max − z min ≥ 3. In this case, we can run the same steps as in the proof of Theorem 1. That is, by multiplying Hamiltonian local terms, we can reduce the support of A to be within at most two adjacent layers. At this end, we can use the results above to show that A after the reduction, must be an identity operator. Thus, the operator A we begin with, is a product of Hamiltonian local terms.
Finally, we consider that the operator has multiple disconnected components. In this case, we take the operator as a single component operator, which are identity operators on some layers. Then through the steps above, we can conclude that the operator must be either an identity operator, or a product of Hamiltonian local terms. This completes our proof.