{'fraction_non_alphanumeric': 0.033907788719785135, 'fraction_numerical': 0.00739984333034915, 'mean_word_length': 4.207532051282051, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topological Entanglement Entropy of Fracton Stabilizer Codes

Entanglement entropy provides a powerful characterization of two-dimensional gapped topological phases of quantum matter, intimately tied to their description by topological quantum field theories (TQFTs). Fracton topological orders are three-dimensional gapped topologically ordered states of matter, but the existence of a TQFT description for these phases remains an open question. We show that three-dimensional fracton phases are nevertheless characterized, at least partially, by universal structure in the entanglement entropy of their ground state wave functions. We explicitly compute the entanglement entropy for two archetypal fracton models --- the `X-cube model' and `Haah's code' --- and demonstrate the existence of a topological contribution that scales linearly in subsystem size. We show via Schrieffer-Wolff transformations that the topological entanglement of fracton models is robust against arbitrary local perturbations of the Hamiltonian. Finally, we argue that these results may be extended to characterize localization-protected fracton topological order in excited states of disordered fracton models.


I. INTRODUCTION
The study of topological order is a major theme of modern condensed matter physics. Gapped topologically ordered states of matter are characterized by properties 1-4 such as ground state degeneracy on manifolds of non-trivial topology, the inability to distinguish distinct such ground states via local measurements, and the existence of local excitations that cannot be created by purely local operators, leading to fractionalization of statistics. In two spatial dimensions there exists a fairly complete understanding of topological order, described via topological quantum field theory (TQFT) and related ideas. 5 However, our understanding of topological order in three spatial dimensions remains incomplete -a lacuna brought into sharp relief by the development of the fracton models [6][7][8][9][10][11] . These models exhibit the characteristic properties of topological order, but lack a description in terms of TQFTs. As such, they represent a new chapter in the story of topological order, and have begun to draw intensive interest [12][13][14][15][16][17][18][19][20][21] .
In a parallel line of development, ideas from quantum information have been transplanted to the field of correlated systems, providing potent methods with which to characterize and study complex quantum many-body systems. Foremost among these new tools is the use of quantum entanglement -specifically, its characterization in terms of entanglement entropy -in developing the understanding of topological order (see e.g. Ref. 22 for a review). For two-dimensional gapped phases, the entanglement entropy contains a universal subleading 'constant' term 23,24 . This is intimately related to the TQFT description of the topological phase and provides a partial characterization of the nature of topological order. It is thus natural to ask, what insights are afforded by studying entanglement in fracton phases of matter? A previous work 25 has introduced certain bounds on the non-local part of the entanglement entropy in fracton models.
In this work, we discuss the entanglement entropy in fracton phases. We begin by reviewing some basic facts about entanglement entropy in lattice models, and the formalism developed to study the entanglement structure of stabilizer codes. We then explicitly calculate the entanglement entropy for several d = 3 models, including two paradigmatic fracton phases: the 'X-cube model' and 'Haah's code'. Both fracton phases are found to exhibit a term of non-local origin in the entanglement entropy that scales linearly with subsystem size, with a coefficient that we calculate. We refer to this term as topological entanglement entropy.
We also point out that the existence of topological entanglement entropy indicates that on the boundary of the subregion, certain non-local constraints act on the ground state wavefunction. We explicitly identify these constraints for the X cube model. This provides an argument that the topological entanglement entropy is robust under arbitrary local perturbations of the Hamiltonian. Consider a perturbation H = λV , where λ is a control parameter, and V is of unit norm and contains perturbations local in real space which do not commute with H(λ = 0). When λ = 0, the stabilizers are no longer eigenoperators of the groundstate. However, for λ sufficiently small that it does not induce a phase transition as increased from zero, the new ground state is related to the λ = 0 ground state by a unitary transformation where U is a local unitary operator which can in principle be constructed by adapting the method of Schrieffer-Wolff transformations 26 (see Appendix C for details). Non-local constraints on the wavefunction cannot be altered by local unitary transformations; thus the topological entanglement entropy (being of non-local origin) must be invariant under local perturbations. Finally we generalize this argument to excited states of disordered fracton models, which we argue can dis-play localization protected fracton order, characterized by topological entanglement entropy.

II. REVIEW OF TOPOLOGICAL ENTANGLEMENT ENTROPY
We begin by reviewing some essential facts about entanglement entropy, focusing in particular on topological contributions and how to extract them. We will consider a system in a topological phase, with |Ψ one of its degenerate ground states. We divide the system into a small subregion A and 'everything else' (B), and construct the reduced density matrix ρ A = Tr B |Ψ Ψ|. The entanglement entropy is then defined to be the von Neumann entropy of ρ A , (Here and throughout, we measure logarithms in base 2 to remove unwieldy factors of ln 2.) In a gapped phase -such as the fracton phases considered in this paperthe entanglement entropy of a ground state is expected to follow an 'area law' i.e. to be proportional at leading order to the surface area of the subregion A. If A has linear size R, then we may expand in powers of R, viz.
where d is the spatial dimension. The leading term is non-universal and dominated by short-distance physics. The 'topological' information is contained in the subleading corrections; the challenge is to extract it. For concreteness, let us review how this works in d = 2, where we have and the topological contribution is the constant piece, S topo = −cγ, with c the number of connected components of the boundary of A. In principle it appears as though extracting γ should be a straightforward exercise: we should simply compute S, and extract its constant contribution and identify it as topological using the scaling form Eq. (4). In practice, however, most systems of interest and all those we are concerned with in this paper are defined on a lattice. Then subregions often have sharp corners that can lead to non-universal constant contributions in Eq. (4), assuming the entanglement cut is a sequence of edges forming a continuous path. This complicates a direct identification of S topo from the scaling of the bipartite entanglement entropy in Eq. (4). However, there are prescriptions 23,24 to extract topological entanglement by suitably combining the results for a variety of bipartitions. Two such prescriptions in two dimensions are illustrated in Fig. 1. We refer to the type of prescription illustrated in Fig. 1(a) as an ABC prescription, where A, B and C are three disjoint regions, and the topological entanglement entropy is given by An alternate prescription is shown in Fig. 1(b) and referred to as a PQWT prescription. In this case, the regions P , Q, W and T have the properties that P = Q ∪ W and T = Q ∩ W , and the topological entanglement entropy is given by In both cases, these linear combinations of entropies are chosen to ensure that the dependence on local contributions from boundaries, including corner contributions, cancels out. These topological entanglement entropies are related to Eq. (4) by S ABC topo = −γ and S P QW T topo = −2γ, which can be understood by counting the number of connected components in the boundaries of the various regions involved.
Both ABC type and PQWT type prescriptions have been generalized 27,28 to d = 3, and we will make use of these generalizations in this paper. While different geometries and topologies of the regions are possible in these generalizations, in ABC type prescriptions we require the regions A, B and C be disjoint. In contrast, in the PQWT type prescriptions we employ, the set theoretic properties P = Q ∪ W and T = Q ∩ W will always be satisfied.
Much progress has been made in understanding topological order and topological entanglement in d = 2 by linking these ideas to TQFTs. For instance, we may understand the dependence of the topological contribution on the number of connected components of the boundary by recognizing that S topo reflects the additional information obtained by counting field lines of an gauge field constrained by a lattice analog of Gauss's law. Entanglement entropy in turn can provide an important tool for extracting TQFT data 22 , such as the braiding and statistics of the fractionalized excitations 29 .
Far less is known about three-dimensional topological phases, particularly in situations where there is no TQFT description. In this paper, we study the topological entanglement entropy of fracton models that are paradigmatic examples of these unexplored topological orders. To do so, we leverage an approach complementary to TQFT, namely we study zero-correlation-length Hamiltonians that are sums of commuting projectors, which describe special exactly-solvable points within some topological phases. More specifically, we consider 'stabilizer Hamiltonians,' where the commuting terms in the Hamiltonian are products of Pauli operators. Our strategy will be to first compute entanglement entropy at solvable points, and then argue that the entanglement entropy has a topological contribution that persists under perturbation, using various three-dimensional generalizations of the ABC and PQWT prescriptions for obtaining topological entanglement entropy.

A. Entanglement entropy for stabilizer codes
We begin our analysis by reviewing the method developed in Ref. [30][31][32] for calculating the entanglement entropy for the ground states of stabilizer codes, and discussing the straightforward extension to arbitrary eigenstates. This extension will be used in Sec. VI when we discuss localization protected fracton topological order in excited states. We will illustrate the use of the stabilizer formalism on the toric code before applying these techniques to fracton models. Readers familiar with the stabilizer formalism for ground-state entanglement computations may wish to skip ahead to Sec. V B.
Throughout, we will consider systems of spin-1/2 degrees of freedom ('qubits') that reside on either vertices or links of d-dimensional hypercubic lattices. The Hamiltonians is built solely of terms containing either X type stabilizer operators that are products of Pauli matrices X i , or Z type stabilizer operators that are products of Pauli matrices Z i . The stabilizer operators all mutually commute, so that the Hamiltonian is exactly solvable. By construction, both X-and Z-type stabilizers square to the identity and hence have eigenvalues ±1.
Let us take S = {O s } to be a set of mutually commuting stabilizers. Then, a state |ψ is stabilized by S if it satisfies for all O s ∈ S. The set of states stabilized by S is called a stabilizer code and is the ground state manifold of the Hamiltonian, H stab = − s J s O s with J s > 0. Here, we are assuming that the O s operators are local, but, below, we relax this assumption. It is useful to consider the Abelian group G that is multiplicatively generated by the stabilizers in S. We assume that all the stabilizers in S are independent and collectively act on exactly N spins-1/2. For open boundary conditions, this is generally expected to be the case for a maximal set of independent local stabilizers. For periodic boundary conditions, we usually need to include some non-local stabilizers that, for instance, wrap noncontractible loops. Elements g ∈ G may be labeled by a binary vector n = {n 1 , n 2 , . . . , n |S| }, with n i ∈ {0, 1} via where |S| is the number of stabilizers in S. Note that S is a generating set. Because G is Abelian, we may label states in the Hilbert space by their eigenvalues of the group elements g ∈ G.
We note that G can be viewed as a vector space over the two-element field F 2 , a fact that will be useful in our approach to numerical calculation of entanglement entropy. This statement holds for any group, like G, that is isomorphic to a product of Z 2 factors. Vector addition corresponds to multiplication of stabilizers. The zero vector corresponds to the group identity. Scalar multiplication is trivial; multiplication by 1 ∈ F 2 is the identity operation, and multiplication by 0 ∈ F 2 sends any element of g ∈ G to the identity of G. We will use both vector space and group language to describe operations in G, and will do this without comment when the meaning is clear from the context. We note that, in vector space language, the set S is a basis.
It remains to determine the size of the group G, that we denote |G|. We further suppose that there is a unique eigenstate |ψ of eigenvalue +1 for all g ∈ G; if this is not the case, then we may add elements to S until this is so. Then, G cannot contain any pure scalar element η = 1 since such an element must have an eigenvalue that is not 1 (a scalar element is an element proportional to the identity). We may then write the projector onto |ψ as (To see this, observe using group properties that the RHS squares to itself, i.e. is a projector, and acts as the identity on |ψ , whence the result Eq. (9) follows by the uniqueness of the ground state.) Taking the trace on both sides, we find where we use the uniqueness of the identity I ∈ G and the fact that any non-identity element in G is traceless; the final step follows simply from the fact that the Hilbert space is the tensor product of N spins-1/2. Thus, we see that the size of G is the full dimension of the Hilbert space, |G| = 2 N . From this and Eq. (8) we conclude that there are N independent stabilizers in S. i.e. |S| = N . Furthermore, we see that the group G must be isomorphic to Z ⊗N 2 , the group of spin-flips on N spins (intuitively, we may think of the 'spin' as the eigenvalue of the stabilizer O s .) Therefore, we may completely label eigenstates of H stab in terms of irreducible representations of the spinflip group 31 . The irreps are one-dimensional and are labeled by a binary string k of length N , and defined by the map ρ k [g( n)] = (−1) k· n ; in other words, the eigenvalue of a group element g in state k is given by (−1) k· n . It then follows from standard orthogonality relations that the density matrix of state | k can be written as In this notation, the ground state |ψ corresponds to the vector k = (0, 0, . . .).
We will now demonstrate how to compute the entanglement entropy of a subregion for any such eigenstate | k of H stab . Consider any bipartition (A, B) as described above. Then, we see that the reduced density matrix of subregion A in state | k is Any element g that is not equal to the identity on B must contain at least one X or Z operator acting in B and consequently will have vanishing trace in B. Thus, the only non-zero contributions to the sum in Eq. (12) are from operators supported only on A, i.e. equal to the identity on B. Operators in the sum supported only on A form a subgroup G A , with elements g( n A ), and irreps labeled by k A , where we have made the obvious generalizations of notation. Since the identity on B has trace TrI B = 2 N B where N B is the number of spins in B, we have (using |G| = 2 N = 2 N A +N B ) that Note that the group G A only includes complete stabilizers in A, since stabilizers that 'dangle' across the entanglement cut are not equal to the identity on B. Since Eq. (13) expresses ρ A as a projector, its entanglement entropy follows straightforwardly: Thus, the computation of entanglement entropy of a subregion A reduces to that of determining the size of the subgroup G A ⊂ G that consists of stabilizers contained entirely within A. In vector space language, log 2 |G A | = dim G A , so Viewed as an F 2 vector space, A and G X A are vector spaces of Z and X stabilizers, respectively. Therefore so that Z and X stabilizers can be treated separately in computing S A .

B. Local and non-local stabilizers
Here, we briefly discuss some further properties of the stabilizer group G A in connection with entanglement entropy and topological entanglement entropy. These properties are then used in the numerical procedure for computing entanglement entropy described in Sec. III C. Moreover, we show that S P QW T topo is entirely determined by counting non-local stabilizers, while in many cases S ABC topo can be determined by counting local stabilizers.
Given a stabilizer Hamiltonian, we obtain a set S Z loc of local Z stabilizers, which are just the terms of the Hamiltonian. (For simplicity of discussion we focus for the moment on Z stabilizers; identical statements hold for X stabilizers.) Local Z-stabilizers supported entirely in A generate a subgroup G Z A,loc ⊂ G Z A . If G Z A = G Z A,loc , we say that G Z A is locally generated. In this case, dim G Z A can be obtained by simple counting, accounting for possible constraint equations satisfied by the local stabilizers. (A constraint is a product of local stabilizers which evaluates to the identity, i.e.
A,loc , we call g a non-local stabilizer. In addition to dim G Z A , we will be interested in the number of independent non-local stabilizers, which is defined by By taking the quotient, we are counting non-local stabilizers up to multiplication by local stabilizers. That is, two non-local stabilizers in G Z A that are related by a product of local stabilizers in G Z A,loc are not considered independent in this counting. We have Now we consider topological entanglement entropy S ABC topo , obtained in an ABC type prescription as discussed in Sec. II. We assume that the stabilizer groups G A , G B , G C , G AB , and so on, are all locally generated. Moreover, we assume the generators of these groups are drawn from a set of local stabilizers that do not obey any local constraints. This holds trivially for the d = 2 toric code and Haah's code, where the stabilizers do not obey any local constraints. For the X-cube model, the vertex stabilizers do obey local constraints, but for suitable regions A, B, C it is possible to use only xy plane and xz plane vertex stabilizers, which do not obey local constraints, as shown in Appendix A 1. Under these assumptions, −S topo is simply the number of independent local stabilizers that have non-trivial support in each of the disjoint regions A, B and C. The contributions of other local stabilizers to S ABC topo cancel out. For instance, suppose some local stabilizer is supported entirely in A. Then it is also contained in AB, AC and ABC, and it contributes +1 to each of dim G A , dim G AB , dim G AC and dim G ABC . These contributions cancel in S topo . Similarly, if a local stabilizer has non-trivial support in both A and B, but not C, then it is contained in both AB and ABC, and its contribution to S topo cancels. A more careful argument for this result is given in Appendix B.
In contrast, S P QW T topo , the topological entanglement entropy obtained via a PQWT type prescription is determined in many cases entirely by non-local stabilizers.
Here we give a rough argument that all contributions of local stabilizers cancel out; a more complete treatment is given in Appendix B. First recall that P = Q ∪ W = Q ∪ W ∪ T , so it is enough to consider different types of local stabilizers supported entirely on P . Suppose a local stabilizer is supported entirely in T , then it is contained in all four regions, and its contribution to S P QW T topo cancels. Now suppose a stabilizer is contained in Q but not in T , then it is also contained in P but not in W , so its contribution cancels. This covers all the possibilities for local stabilizers. Because N P + N T = N Q + N W , we have the result (19) where Ω P = Ω Z P + Ω X P is the total number of nonlocal stabilizers in P , and similarly for the other regions. ∆ P QW T is a correction associated with non-local constraints among local stabilizers. In Appendix B, we derive Eq. (19), and show that ∆ P QW T = 0 for all the models discussed in this paper except the d = 3 toric code.

C. Calculating the entanglement entropy
In many cases, it is possible to determine |G A | and hence the entanglement entropy S A analytically. However, this is not always straightforward, and numerical calculation is useful as a check on other methods and sometimes as a primary means of determining the entanglement entropy. Here, we describe a numerical procedure to determine dim G A = log 2 |G A |. We are always interested in the case where A is a subset of a thermodynamically large region. Under a suitable assumption discussed below, which can be verified for particular models of interest, dim G A does not depend on global properties of the large region containing A.
To proceed, we choose a finite enclosing region B with We make the assumption that it is always possible to choose B so that G Z B is locally generated. This assumption implies G Z A ⊂ G Z B,loc ; that is, all stabilizers in A are products of local stabilizers in B. We show in Appendix A that this assumption holds for the stabilizer codes studied in this paper.
Now we introduce F Z B,loc , the group of formal products of local Z stabilizers supported on B. We let M be the number of local Z stabilizers supported entirely in B, and denote these operators by O s (s = 1, . . . , M ). A general product of these stabilizers is can correspond to the same operator, because there can be constraints among the local stabilizers. There is a linear map where P Z B is the group (or F 2 vector space) generated by all Z Pauli operators supported on B. The map φ Z is defined by replacing the O s in a formal product with their expressions in terms of Pauli operators. G Z B,loc is a subspace of P Z B , and moreover G Z B,loc = Imφ Z . If f ∈ Kerφ Z , then f is a formal product of stabilizers that evaluates to the identity operator. This happens when the stabilizers obey some constraint equations. Indeed, the number of independent such constraints is precisely dim Kerφ Z .
We note that We thus have the projection map defined by Elements in the kernel of π B−A are products of Z Pauli operators supported entirely on A. The last step is to consider the composition This means that f corresponds to a non-trivial element of G Z A . We are finally able to express the number of stabilizers in A as 33 Here, the last term is subtracted to avoid incorrectly counting constraints among local stabilizers as non-trivial elements of G Z A . The linear maps φ Z and π B−A can be constructed explicitly as matrices, and numerical linear algebra methods can then be used to compute the dimensions of the kernels. Our calculations were done using routines for linear algebra over F 2 in Mathematica.
There is a minor modification of the above approach that significantly reduces the computational effort required. We consider the subspace F Z A,loc ⊂ F Z B,loc , which consists of formal products of stabilizers supported entirely in A. Then, we have F Z B,loc = F Z A,loc ⊕ F where F consists of formal products of stabilizers in B that either lie completely outside A or are not fully contained in A.
The number of such stabilizers is It can happen that φ Z (f ) is a local stabilizer in G Z A ; this can happen when the local stabilizers obey some local constraints, so that some local stabilizers in G Z A can be written as a product of local stabilizers not supported entirely in A. We let K Z A,loc ⊂ G Z A,loc be the subspace generated by such local stabilizers in A. If the local stabilizers generating G Z B,loc obey no local constraints, theñ which determines dim G Z A via Eq. (18). We use this method to compute topological entanglement entropy for Haah's code as discussed in Sec. V C. In addition, we employ the same method to check our results for the X-cube model in Sec. V B. To illustrate the use of the stabilizer formalism, we now use it to compute the entanglement entropy of the d = 2 toric code model 34 . One qubit resides on each link of the square lattice, and the Hamiltonian is where B p is the product of the four Z operators surrounding the plaquette p, and A v is the product of the four X operators connected to vertex v. We compute the entanglement entropy of a subsystem A of size R × R, shown in Fig. 2 for R = 2. This region contains N A = 2(R + 1)R spins. Additionally, there are R 2 plaquette terms and (R − 1) 2 vertex terms confined entirely within subregion A; these are local stabilizers and they are all independent. Using the fact that A is simply connected, it can be seen easily that G A is generated by the plaquette and vertex stabilizers supported entirely in A. Moreover, these stabilizers obey no local constraints. Therefore, |G A | = 2 R 2 +(R−1) 2 .
Using Eq. (14), we find that the entanglement entropy is Since the boundary has length |∂A| = 4R measured in lattice edges, we have S A = |∂A| + S topo , where S topo = −1 is the well-known topological entanglement entropy of the two dimensional toric code. This is something of an accident; this region has sharp corners, but the corner contributions happen to vanish. To actually obtain S topo , we can break region A into three subregions A, B, C as in Fig. 1a, and apply the ABC prescription to obtain S ABC topo as in Eq. (5).
To illustrate features not arising in the above simple example, we calculate S P for an annulus-shaped region P shown in Fig. (1b). There, due to the non-trivial topology, G P is no longer generated by local stabilizers supported entirely on P . Using the terminology of Sec. III B, G P is generated by the local stabilizers in P , together with two non-local stabilizers. These are products of plaquette and vertex stabilizers, respectively, over the hole in the annulus. Taking the linear size of the hole to be R and that of the exterior edges to be 3R, we obtain S P = 16R − 6 = |∂P | − 4 − 2, where the constant −4 is non-universal and contributed by vertex terms at four concave corners of the subsystem while the −2 results from the number of non-local plaquette and vertex stabilizers acting on this subsystem. Although there are non-universal parts in this entropy due to the detailed geometry, the topological part of the entropy can be extracted by canceling all of those boundary contributions out via the PQWT prescription in Eq. (6).
A different perspective on the topological entanglement entropy, that further clarifies its robustness, is afforded by an understanding of the ground state wave function of the toric code as a condensate of closed loops of Z 2 'electric' field lines. Each connected component of the entanglement surface intersects any loop an even number of times; this topological fact provides exactly one bit of extra information about the ground state, thereby reducing the entanglement entropy by a universal correction of −1 for each connected component of the entanglement surface. Similar ideas can be used to clarify the topological entanglement entropy of other topological orders in two dimensions by constructing their ground states as string-net condensates 24,35 .
As we have discussed, extracting the topological contribution to the entanglement entropy generally requires employing a prescription designed to cancel local contributions. This becomes even more important in d = 3, where Ref. 28 argued that even for a subregion A with a smooth boundary, in the absence of parity and continuous rotation symmetry, every term in the expansion of S A in powers of the inverse linear size R −1 acquires a local, non-universal contribution. Continuous rotation symmetry is absent in fracton models, and, more seriously, the dynamics of fracton excitations leads to an expectation that fracton topological orders lack a continuum description with continuous rotation symmetry. Discrete symmetries like parity may or may not be present in a given solvable fracton model, but certainly need not be present upon perturbing such a model to make it generic. Therefore, we rely on d = 3 generalizations of the ABC and PQWT prescriptions discussed in Sec. II. We use two different PQWT type prescriptions 27 , one is illustrated in Fig. 4, the other in Fig. 8. A naïve extension of the d = 2 ABC prescription ( Fig. 1(a)) is shown in Fig. 3. As noted in Ref. 28, this prescription fails to cancel local contributions from the two points where regions A, B and C all meet at the top and bottom boundaries. This implies that S ABC topo is contaminated by non-universal contributions that are constant in R. However, we will still employ this prescription, because in fracton models we will find a contribution to S topo proportional to R, which is unaffected by the uncanceled constant local contributions. While we do not use them in this paper, we note that Ref. 28 introduced different d = 3 ABC prescriptions that do not suffer from this issue.
To illustrate our approach in d = 3, we briefly discuss the calculation of topological entanglement entropy for the d = 3 toric code, using the PQWT prescription illustrated in Fig. 4. The Hamiltonian is given by Eq. (26) but on the cubic lattice, so that the vertex terms involve six spins, and there are plaquette terms for each face of a cubic unit cell. As found in Eq. 19, S P QW T topo is given entirely by counting non-local stabilizers in each region. Region P has a single non-local stabilizer, which can be constructed by taking a product of local plaquette stabilizers over an xy plane surface that cuts through the inner cube (i.e. the "hole" at the center of P ), as illustrated in Fig. 5. Regions Q, W and T have no non-local stabilizers, so S P QW T topo = −Ω P = −1. This is consistent with the fact that the ground state wave function is a loop condensate, with loops cutting the entanglement surface an even number of times. Ref. 28 computed the topological entanglement entropy of the d = 3 toric code using a different prescription, which produces the same result.

V. TOPOLOGICAL ENTANGLEMENT ENTROPY OF FRACTON MODELS
In this section, we compute the topological entanglement entropy of fracton models. We begin (Sec. V A) with a simple warm-up example, a stack of decoupled two-dimensional toric code layers, which is related to the X-cube model by a coupled-layer construction. 15,16 Then we proceed to consider the X-cube model (Sec. V B) and Haah's code (Sec. V C).

A. Decoupled layers of d = 2 toric codes
We consider a stack of d = 2 toric codes, with layers spaced uniformly and arranged normal to the z axis. We consider the "square torus" PQWT prescription shown in Fig. 4, with topological entanglement entropy defined by Eq. (6). The intersection of each layer with P is the annulus-shaped region discussed in Sec. III D, which has two non-local stabilizers. The regions Q, W and T do not support non-local stabilizers, so the topological entanglement entropy is two bits per layer. Therefore, choosing the unit of length to be the layer spacing, and choosing the regions P , Q, W , T to contain precisely R toric code layers, we obtain In this example of decoupled layers, it is no surprise that we obtain a R-linear term in S topo . This term is of interest because it also appears in fracton models. We remark that slightly changing the detailed specificiation of the regions can alter the constant term in S P QW T topo . For instance, instead of stating the regions contain precisely R layers, we could define them to extend a distance R along the z-axis between toric code layers on the "bottom" and "top" surfaces, and to contain these surface layers. Then each region intersects R + 1 layers, and we obtain S P QW T topo = −2R−2. Note that this changes the constant term by an even integer; it is possible that this term does have a robust meaning modulo two. In a generic system in the same phase as the decoupled stack we are considering, we would not have precise control over the number of layers intersected by these regions, so that only the R-linear term is clearly meaningful. This illustrates a general point that, in d = 3, one should be cautious in ascribing any meaning to the constant term in S topo when an R-linear term is present.
We also consider the topological entanglement entropy using the d = 3 ABC prescription shown in Fig. 3, with the preferred axis along the z-axis. Considering a single layer, this reduces to the d = 2 ABC prescription of Fig. 1, giving a topological entanglement entropy of one bit per layer. Again choosing the regions to intersect precisely R layers, we have The example of decoupled layers of d = 2 toric codes is also instructive in that it illustrates the answer can depend on the orientation of the entanglement cut. Indeed, different results would be obtained for either the ABC or PQWT prescriptions, if we choose the preferred axis to lie in an arbitrary direction, since the number of intersecting layers would be different. Although this somewhat complicates the interpretation of S topo , it has the advantage of providing a means to identify the "natural" axes for entanglement in this system, by rotating the orientation of the preferred axis so as to obtain a maximal answer. In fracton models, this could potentially help to discover new coupled-layer constructions along the lines of Refs. 15 and 16.
We also consider a stacking of decoupled toric code layers along the x, y and z directions simultaneously, as in the coupled-layer construction of the X-cube model. 15,16 In this case, we obtain the same results for topological entropy using the two prescriptions employed above, because these prescriptions do not capture the topological entanglement of the layers normal to the x and y axes. We are not aware of a single ABC or PQWT type prescription that captures all of the topological entanglement in this system in one shot. Instead, it seems to be necessary to compute S topo for different sets of regions to obtain a full picture of the non-local entanglement. A notion of 'recoverable information' that can be defined for stabilizer codes offers a complementary perspective to the one in the present paper and may achieve this goal.

B. X-cube model
We now apply the stabilizer formalism to compute our first new result: namely, the entanglement entropy of the X-cube model, an archetypal example of a 'type I' fracton phase 11 . The model is defined on a cubic lattice with a spin 1/2 variable on each link, with Hamiltonian where the A-type stabilizers involve a product of four Z i operators that surround a vertex in one of three orthogonal planes, and the B-type stabilizers involve a product of twelve X i operators around a elementary cube.
To compute the topological entanglement entropy, we employ the construction Eq. (6), using the regions shown in Fig. 4, and taking the preferred axis to be the z-axis. As discussed in Sec. III B, S P QW T topo is given by counting non-local stabilizers via Eq. (19). It follows from Appendix A 1 that regions Q, W and T have only local stabilizers (i.e. their stabilizer groups are locally generated), so we have S P QW T topo = −Ω P , with Ω P the number of non-local stabilizers in P .
To determine Ω P = Ω X P + Ω Z P , we first note that a product of cube terms over an xy plane layer cutting through the "hole" in the center of P produces a non- 6. (a)A non-local Z stabilizer (red) of the X-cube model in region P is constructed by taking a product of vertex terms over the "hole" in the center of P . (b) A non-local X stabilizer (red) of the X-cube model in region P is constructed by taking a product of cube terms over the "hole" in the center of P . local X stabilizer, as shown in Fig. 6. There are R different layers, and we expect that the resulting R nonlocal stabilizers are independent in the sense that they cannot be deformed into one another by taking products with local stabilizers, so that Ω X P = R. The reason for this expectation is that each of these non-local stabilizers is a closed-loop string operator for a distinct non-trivial quasi-particle excitation confined to move in the corresponding xy plane. 11,16 It should not be possible to change the particle type of a string operator by multiplying it with local operators. We also verified that Ω X P = R using the numerical method of Sec. III C, for R = 2, . . . , 10. Similarly, taking a product of xy-plane "vertex type" stabilizers over the hole in P also gives a non-local stabilizer. There are R+1 such xy-plane layers, giving Ω Z P = R + 1, which we again verified numerically for R = 2, . . . , 8. Therefore we find Ω P = 2R + 1 and We note that the R-linear term in this result is identical to that obtained in a stack of decoupled d = 2 toric codes along x, y and z axes; as discussed in Sec. V A, only the layers normal to the z-axis contribute to the topological entanglement entropy for this choice of the regions P , Q, W , T . As in the case of the toric code, the linear term of this result can be understood within a loop condensate picture. In the Z basis, configurations satisfying the vertex terms of the Hamiltonian can be viewed in terms of strings of flipped links with Z = −1, where in every {100} plane the strings form closed loops. Each {100} plane thus gives a non-local contribution of −c to S A , where c is the number of connected components in the intersection of the plane with the boundary of A. Applying this simple rule to S P QW T topo , we find that each xy plane contributes −2, while yz and xz planes do not contribute; this reproduces the −2R term obtained above. We note that this is the same loop condensate picture as for a stack of decoupled toric codes. In the X-cube model, the layer-by-layer loop constraints in the ground state wave function are not truly independent; a more detailed analysis taking this into account would presumably also reproduce the constant term in S P QW T topo .
It is evident that, as for the case of decoupled layers of d = 2 toric codes, the number of topological constraints (and hence the topological entanglement entropy) will depend on the orientation of the entanglement surface. Our discussion here is for an entanglement surface aligned with the symmetry axes of the problem. The entanglement entropy for arbitrary orientations could be evaluated using analogous methods, but we do not discuss it further here.
We also compute the topological entanglement entropy using the ABC prescription, with regions shown in Fig. 3. Because regions A, B, C, AB and so on all have locally generated stabilizer groups generated by cube stabilizers and xy and xz plane vertex stabilizers (Appendix A 1), following Sec. III B and Appendix B, S ABC topo is determined by counting local stabilizers whose support is split among all of A, B and C. The detailed geometry can be chosen so that this only occurs for cube stabilizers, and the number of these cube stabilizers is R. Therefore we find where the vanishing of the constant term is presumably unimportant, because this prescription does not cancel all constant local contributions. We observe that the coefficients of the topological entanglement entropy using both ABC and PQWT prescriptions in the X-cube model are exactly as in the corresponding system of decoupled toric codes. Evidently, the linear term in the topological entanglement entropy is insensitive to the m-string condensation that occurs going from the decoupled toric codes to X-cube phase. 16

C. Haah's code
We now turn to Haah's code, the archetype of a 'type-II' fracton model 11 . This model is also defined on a cubic lattice, but now with two spin-1/2 variables on every vertex. The Z and X type stabilizer operators now consist of products of Z and X Pauli matrices around an elementary cube with a Hamiltonian of the form where A c and B c denote the products of Z and X Pauli operators around the vertices of a cube specified in Fig. 7. To compute the topological entanglement entropy, we first employ the ABC prescription with regions as shown in Fig. 3. It is shown in Appendix A 2 that the stabilizer groups for regions A, B, C, AB, and so on are all locally generated. Therefore, following the discussion of Sec. III B, −S ABC topo is the number of local stabilizers contained in ABC that have support split among all three regions A,B, C. These stabilizers reside on cubes along the axis where the three regions meet, and there are 2(R − 1) of them, where the factor of 2 accounts for counting both X and Z stabilizers. Therefore, where only the R-linear term is expected to have any universal meaning. We also consider the topological entanglement entropy captured by two different PQWT prescriptions, determined using the numerical method of Sec. III C. We first simplify the problem using the spatial inversion symmetry of Haah's code, which acts non-trivially on the spins, sending X → Z and Z → −X, and also exchanging the two qubits on each site. If a region A is inversionsymmetric, then Ω Z A = Ω X A . On the other hand, if two regions A and B are related to one another by inversion, then Ω Z A = Ω X B and Ω X A = Ω Z B . In the PQWT prescription of Fig. 4, regions P and T are inversion-symmetric, while inversion exchanges Q and W . This implies Our numerical calculations lead to the conclusion that S P QW T topo = 0 for the regions of Fig. 4. (In more detail, we show in Appendix A 2 that the stabilizer groups for regions Q, W and T are all locally generated, i.e. Ω Z Q = Ω Z W = Ω Z T = 0. Numerically, we find that Ω Z P = 0 for R = 4, . . . , 11. While Ω Z P = 2 for R = 2 and Ω Z P = 1 for R = 3, this seems to be a finite-size effect.) This result is strikingly different from S ABC topo , while the linear term in these two entropies only differed by a factor of two for the X-cube model. The contrast with the X-cube model suggests that Haah's code may not have a coupled-layer description where the layers lie in {100} planes.
To find a different PQWT prescription that does capture some of the non-local entanglement in Haah's code, we note that we should not expect the non-local stabilizers of Haah's code to be one-dimensional objects, as they are in the X-cube model. This expectation is based on the fact that none of the topologically charged excitations in Haah's code can be transported by string operators, so we should expect that any non-local stabilizers are higher-dimensional objects. Moreover, this expectation is further substantiated by the fact that Ω Z P = 0 for the solid torus region of Fig. 4. This motivates us to employ the PQWT prescription of Ref. 27, with regions shown in Fig. 8. Here, the region P is more isotropic, allowing for non-local stabilizers wrapping entirely around the interior cube.
The regions of Fig. 8 indeed give a non-zero result for S P QW T topo . Our numerical results are summarized in Table I, and we find based on numerical calculations up through R = 11. This topological entanglement entropy also has a R-linear term. It is interesting to remark that, while Haah's code has a well-known intricate dependence of the ground state degeneracy on system size 36 , the behavior of the topological entanglement entropy is much simpler.  Fig. 8. The functional forms shown are exact from R = 4 up through R = 11 (the largest value of R for which calculations were done). These results determine S P QW T topo via Eq. (35).
At this point it is natural to ask if the topological entanglement entropy in Haah's code may also be given a geometrical interpretation in terms of constraints on the structure of the ground state wavefunction, like the loop condensate picture for toric code and X-cube models. Because isolated fractons in Haah's code are created at corners of fractal operators, it seems likely that the ground state can be viewed as some kind of condensate of fractal objects. At present, to our knowledge there is no clearer meaning that can be given to this picture, or whether it can shed light on non-local entanglement. itself when the Hamiltonian is perturbed away from the stabilizer form. The entanglement structure of ground states is then 'protected' from perturbations by the energy gap (as discussed above), whereas the excited states lack such protection. Upon perturbation, the excited states of translationally invariant fracton models are thus expected to thermalize 17 to volume law entanglement, in the process losing their topological order. Indeed, if one tries to construct the unitary U of Eq. (1) using Schrieffer-Wolff perturbation theory outside the ground space, one finds that the perturbation theory diverges, suggesting that U may not be a local unitary. Since our discussion thus far relied on the existence of a local unitary U , it has nothing further to say about the translationally invariant case.
However, as noted by Ref. 37, this scenario could change dramatically once we break translational invariance by introducing quenched disorder: the topological order can be protected even in excited states by localization, in a manner that we now sketch. Consider a disordered fracton Hamiltonian of the form where the J A i and J B j are random numbers drawn from some distribution of width W . While the excitations are frozen (non-propagating) for any W at the stabilizer point, once we add small perturbations they will be able to propagate for small W , but for sufficiently large W the system can enter a many body localized phase 38 , where excitations cannot propagate freely. In this many-body localized phase, topological order (including fracton order) can persist at non-zero energy densities -i.e., even in highly excited states. The challenge is how to detect this topological order. 'Excited state degeneracy' cannot serve as a diagnostic, since the many body level spacing in the middle of the spectrum is exponentially small in the volume of the system and thus there is no longer a distinction between topological and other degeneracies in the thermodynamic limit. In Ref. 37, non-local correlation functions -related to the 'Fredenhagen-Marcu' order parameters familiar to lattice gauge theoristswere argued to be good diagnostics of topological order. These diagnostics were generalized to certain fractonic models in Ref. 39, but nevertheless such non-local correlation functions can be challenging to compute. Here, we will demonstrate that topological entanglement can diagnose fracton topological order in excited states.
The fact that excited states in the localized regime can support fracton topological order follows straightforwardly from our preceding discussions. First note that at λ = 0, excited states share the entanglement properties of the ground state (including topological entanglement), since the excited states are also eigenstates of the stabilizer operators. Now note that, in the localized phase, the unitary transformation U is local, with at most exponential tails. This follows because of the 'mobility gap' in many body localized systems 40 i.e. the Schrieffer-Wolff perturbation theory has matrix elements in the numerator, and matrix elements vanish between near degenerate eigenstates. Again, the topological entanglement entropy is non-local, and is expected to be unaffected by a local unitary transformation, subject to the same caveats discussed in the previous section.
We also note that the dressed stabilizers are simply the local integrals of motion or 'lbits' of the localized system 41,42 , and these must be localized by postulate. For a detailed construction of dressed integrals of motion via Schrieffer-Wolff perturbation theory, showing that these are local in the MBL regime, see Ref. 43. (For the cognoscenti, we note that our argument here parallels more closely the construction of l-bits via Wegner-Wilson flow in Ref. 44).
We note that thus far we have assumed that a many body localized phase can exist in three dimensional lattice models. There is some debate about whether many body localization can arise in spatial dimensions d > 1 with random short range correlated disorder 45 , because of 'thermalizing avalanches' triggered by rare regions. For truly random short range correlated disorder, our discussion applies to systems that are small enough to lack the relevant rare regions, and perhaps also in the thermodynamic limit, if the argument from Ref. 45-47 can be somehow circumvented. However, the problem may also be sidestepped by making the disorder long range correlated or quasiperiodic, such that the 'rare region obstruction' identified in Ref. 45 does not apply.
We therefore conclude that in disordered fracton models, fracton topological order can arise even in highly ex-cited states, where it may be diagnosed through a 'topological entanglement entropy' linear in the size of the subsystem.

VII. CONCLUSIONS
We have explicitly computed the entanglement entropy of two archetypal fracton models -the X-cube model and Haah's code -and have demonstrated the existence of a topological contribution to the entanglement entropy that is linear in the size of the subsystem. At a minimum, this provides a coarse characterization of fracton topological phases, in that for a given system, two states with distinct topological entanglement must be in distinct phases. There is also an obvious extension of this diagnostic to anisotropic models, wherein one separately considers the scaling of entanglement entropy with the size of subregion A in the x, y and z directions respectively, thereby characterizing the topological entanglement with three indices. (More carefully, one would separately consider the topological entanglement entropy for a torus of thickness R oriented in three orthogonal planes). In general the topological entanglement can depend not just on the size of the region, but also on its orientation, which may provide a useful means of diagnosing the symmetry axes of a phase by rotating the entanglement surface to obtain a maximal answer.
What more information could be extracted from a study of entanglement? In two-dimensional topologically ordered phases, a careful analysis of the action of symmetries such as rotation and reflection within the manifold of degenerate ground states can provide insights into the fractionalized statistics of quasiparticle excitations in the phase; whether such manipulations can shed additional light on the properties of fracton excitations (that do not admit a quasiparticle description) remains an open question. It would also be interesting to study the entanglement spectrum 48 , as this may contain more information than is encapsulated in entanglement entropy. Finally, the dynamics of entanglement has in other contexts (see, e.g. Refs. 49 and 50) provided much insight into the nature of thermalization and the approach to or avoidance of equilibrium. We leave investigation of these issues to future work.

X-cube model
In the X-cube model, we show that G A is locally generated for regions A obtained by taking a simply connected region in the xy plane, and stacking this region along the z-axis. (The choice of plane and normal stacking direction is of course arbitrary, due to cubic symmetry.) We assume the boundary of each xy plane slice is a sequence of edges that are contained in A and form a path in the cubic lattice. Two examples of such regions for the d = 2 square lattice are shown in Fig. 9. We further assume that acute corners of the slice (see Fig. 9) are sufficiently isolated from other points on the boundary to carry out the cleaning procedure for Z stabilizers discussed below. A more precise statement of this assumption is given below; essentially, we are assuming the boundary of the slice is not too rough. This condition allows the slice shown in Fig 9a, but rules out that in Fig. 9b. We refer to cubic lattice links oriented along the x, y and z axes as x, y and z links, respectively. In addition to spins residing on x and y links contained in each slice, the region A contains all z-links joining two adjacent slices.
In fact, under the assumptions given, we show an even stronger property: G A is locally generated using only xy-plane and xz-plane vertex stabilizers, and cube stabilizers. (Alternatively, we can use xy-plane and yz-plane vertex stabilizers.) Because there are no local constraints among xy and xz plane vertex stabilizers, using this as a generating set allows us to simply establish Eq. (5) in Appendix B.
We first consider some operator X 0 supported on A, which we take to be an X-stabilizer. By definition, this means that X 0 commutes with all local Z-stabilizers, including those not entirely supported on A. We would like to show that X 0 is a product of local X-stabilizers supported on A, which can be accomplished by a cleaning procedure, where we successively multiply X 0 by such stabilizers until we obtain the identity operator. We denote the X-stabilizer obtained from X 0 at the current stage of the cleaning procedure by X . FIG. 9. Two square lattice slices that can be stacked to obtain a region A of the cubic lattice. Such slices have two kinds of corners, labeled acute and obsute in (a). Our cleaning procedure goes through as long as acute corners are sufficiently isolated; region (a) satisfies this condition, while the circled corners in (b) cause the cleaning procedure for Z stabilizers to fail for the corresponding region A.
We begin with the bottom layer (i.e. smallest z coordinate) of X 0 . In this slice, commutation between xy plane vertex stabilizers is exactly as in the d = 2 toric code. Using the fact that the stacked region is simply connected, this implies that X restricted to this layer is a product of plaquette operators ∈p X , where p is a square plaquette in the xy plane. We can thus clean this layer by multiplying X 0 with a suitable product of local X stabilizers whose cube centers lie just above the slice; this works because the restriction of these stabilizers to the slice are plaquette operators.
It appears this cleaning step may leave dangling z-links lying just above the bottom slice, where X | = X on such links. The notation X | means the restriction of X to the link , i.e. X is a product over links of Pauli operators, and X | is the Pauli operator (1 or X) at the link appearing in the product. By considering commutation of X with xz-and yz-plane vertex stabilizers whose centers lie in the original bottom layer, we see that X | = 1 for all these dangling links.
These steps reduce the height of the stack by one unit cell, resulting in a new region of the with the same properties as the one we started with. Therefore, we can continue this procedure until X is only supported on a single slice. We can show X = 1 by considering commutation with xz and yz plane vertex stabilizers whose centers lie in the same slice. Consider x-links in the slice of interest with some fixed y coordinate. Starting at large negative values of x, and moving in the positive x direction, find the first link with X | = X. Then X anticommutes with the xz vertex stabilizer whose center lies adjacent to this link in the negative x direction. This is a contradiction, and implies X | = 1 for all x-links, and similarly for all y-links. Therefore X = 1 and we have reached the end of the cleaning process. Now we consider a Z stabilizer Z 0 supported on the same region A, and clean it via multiplication with vertex stabilizers. Again we start with the bottom slice, where commutation between Z 0 and cube stabilizers whose centers lie just below the slice is the same as in the d = 2 toric code. Therefore, Z 0 restricted to this slice is a product of xy plane vertex stabilizers, which we clean off. This may leave dangling z-links just above the bottom slice, where Z| = Z. First consider a z-link just above an acute corner. We assume there is a cube stabilizer that contains this link, but contains no other links of A; this is a more precise statement of our assumption that the acute corners are sufficiently isolated. This cube stabilizer anticommutes with Z unless Z| = 1. All dangling z-links not above an acute corner can be cleaned by multiplying with a xz plane vertex stabilizer, whose vertex lies just above the link to be cleaned.
This cleaning process can be continued until Z is supported on a bilayer consisting of the top slice, the slice just below it, and z-links joining these slices. Both the top and bottom slice can be cleaned off as above, leaving only the z-links. Considering the set of z links with Z| = Z, we go to the leftmost column of this layer (i.e. smallest x), and find the z link in this column with smallest y coordinate. We see that Z anticommutes with the cube stabilizer whose center lies in the same layer, and lies diagonally adjacent to the selected link in the negative x, negative y direction. This is a contradiction, so we must have Z = 1, and the cleaning procedure is complete.

Haah's code
In Haah's code, we first show that G A is locally generated for A a rectangular prism, i.e. a region containing all sites with x coordinate satisfying x min ≤ x ≤ x max , and similarly for y and z. Axes are chosen as in Fig. 7. Next, we generalize this to show that G A is locally generated for certain "L-shaped" regions shown in Fig. 11. In all cases, it is enough to concentrate on Z stabilizers: The rectangular prism regions are inversion-symmetric, so the corresponding result for X stabilizers follow from inversion symmetry. The L-shaped regions are not inversionsymmetric, but each region of one type as shown in Fig. 11 is related to a region of a different type under inversion, so that if all four types of regions have locally generated G Z A , the corresponding statement for X stabilizers follows.
Let A be a rectangular prism region, and choose a Z stabilizer Z supported in A. We write where Z r,i is the Z Pauli operator for the ith qubit at position r = (x, y, z), and where the choice of operator is specified by the binary numbers n( r, i) = 0, 1. The restriction Z r is given by Z| r = (Z r,0 ) n( r,0) (Z r,1 ) n( r,1) .
We often suppress some of the indices when using this notation, so for instance if n( r, 0) = n( r, 1) = 1, we write In order for Z to be a stabilizer, it must commute with all the local X stabilizers. We label local stabilizers A c and B c by the position of the corner with smallest values of the coordinates x, y, z, thus writing A c = A r and B c = B r . Then the condition ZB r = B r Z can be written n( r, 0) + n( r, 1) + n( r +x, 1) + n( r +ŷ, 1) + n( r +ẑ, 1) + n( r +x +ŷ, 0) + n( r +x +ẑ, 0) + n( r +ŷ +ẑ, 0) = 0 mod 2. (A4) As discussed above for the X-cube model, we will carry out a cleaning process where we successively multiply Z by local Z stabilizers contained in A until Z = 1. Equation (A4) looks complicated, but it simplifies if we choose r → r 0 to be the corner of A with largest x, y and z coordinates (see Fig. 10), where it reduces to n( r, 0) + n( r, 1) = 0 mod 2. (A5) This implies that either Z| r0 = 1 or Z| r0 = ZZ. Supposing the latter case, we multiply Z by A r0−(x+ŷ+ẑ) , to obtain a new Z with Z| r0 = 1. We repeat this procedure, moving down in the zdirection, to clean operators in the column below r 0 (red dots in the figure) until r 1 is reached. Z| r1 cannot be cleaned in the same way, because the Z stabilizer we would need to act with lies outside A. However, again we have Z| r1 = 1 or Z| r1 = ZZ. Considering commutation of Z with B r1−ẑ , we find n( r 1 , 1) = 0 mod 2. Therefore we have Z| r1 = 1, and no cleaning is needed.
So far we have cleaned the vertical column below r 0 . To proceed, we move to r 2 , and repeat the same procedure to clean the vertical column below r 2 . Proceeding in this way we can clean until we reach the position r 3 . By the same reasoning as before, either Z| r3 = 1 or Z| r3 = ZZ. We cannot apply the same cleaning procedure for r 3 , because the local Z stabilizer we would need to multiply lies outside of A. However, commutation with B r3−x implies either Z| r3 = 1 or Z| r3 = ZI. Therefore Z| r3 = 1 and no cleaning is needed. The same argument applies to all the sites vertically below r 3 , and therefore we have cleaned off the entire +y face of Z. We can repeat the same steps to clean almost all of Z until position r 4 is reached. Again we have either Z| r4 = 1 or Z| r4 = ZZ. Commutation with B r4−ŷ implies either Z| r4 = 1 or Z| r4 = ZI. Therefore Z| r4 = 1 and no cleaning is needed. We proceed vertically below r 4 , and then on to position r 5 , and so on, to see that Z| r = 1 everywhere on the remaining −y face, and no further cleaning is needed. This completes the argument, and we have shown G A is locally generated for A a rectangular prism. Now we consider L-shaped regions of four types, as shown in Fig. 11, and show that G A is locally generated for such regions. Regions of this geometry appear in the computation of S topo via the ABC prescription (see Sec. V C). Each region is a stack along the z-axis for z min ≤ z ≤ z max of an xy plane slice. The figure shows the top such slice (z = z max ) for each region. Similar regions that are stacks along x and y axes can be obtained from these by three-fold rotational symmetry about the [111] axis, and do not need to be considered separately.
For each type of region, we again consider a stabilizer Z supported in the region. In each case, the entire subregion including and below the solid circled regions in the figure can be cleaned off by following the cleaning procedure described above for a rectangular prism region. For type 1 and 2 regions, this results in a stabilizer Z supported within a rectangular prism, which we have already shown is a product of local Z stabilizers. For type 3 and 4 regions, it remains to consider the subregion including and below the dashed ovals. Starting with the type 3 region, we consider the site r 0 in the top layer (see figure). Considering commutation of Z with B r0−x and B r0−x−ŷ implies Z| r0 = 1. The same reasoning shows that Z| r = 1 for all sites in the column below r 0 , and in the whole subregion containing and below the dashed oval. This again reduces the problem to the already solved rectangular prism case. The argument for type 4 regions proceeds essentially the same way, except that we consider commutation of Z with B r0 and B r0−ŷ .
Similarly, we can show that in the C-shaped regions involved in PQWT constructions like Fig. (4), G A is also locally generated. As shown in Fig. (12), there are four possible C-shaped regions. For each type of region, we still consider a stabilizer Z supported in that region. In each case, the entire subregion including and below the solid circled regions in the figure can be cleaned off by following the cleaning procedure described above for a rectangular prism region. For type 1 region, this results in a stabilizer Z supported within a rectangular prism, which we have already shown is a product of local Z stabilizers. For type 2 and 4 regions, the solid and dashed oval regions can be cleaned off in the same way as we clean the L-shaped regions. Then what is left are type 2 and 1 L-shaped regions respectively, shown in Fig. 11, which can be cleaned using the same procedure. For the type 3 region, it is more convenient to start from the lower left corner, i.e. r 0 . The commutation of Z with B r0−ŷ leads to either Z r0 = 1 or Z r0 = ZI. Suppose the latter one is true, then we can multiply by A r0−x−ẑ to clean this site. All sites below r 0 (with the same x and y coordinates) can be cleaned by a similar procedure, until one arrives at the bottom site r 1 , which has z = z min . The commutation of Z with both B r1−ŷ and B r1−ŷ−ẑ requires Z r1 = 1. The next column beginning with r 2 can be cleaned off in the same way. Next, we can move to r 3 . Notice the commutation of Z with both B r3−ŷ and B r3−x−ŷ implies that Z r3 = 1. In a similar fashion we can proceed to clean the whole subregion in the dashed circle, resulting in a type 4 L-shaped region which we have already be able to clean as discussed above. We derive the formula Eq. (19) for S P QW T topo in terms of the number of non-local stabilizers in each region. Treating X and Z stabilizers together, we pick a basis B T,loc for G T,loc . Because T is a subset of all the other regions, we can extend this to bases B Q,loc , B W,loc and B P,loc for G Q,loc , G W,loc and G P,loc . Each of these bases can then be extended to a basis for the full stabilizer group in the corresponding region, including non-local stabilizers. The topological entanglement entropy is To simplify this expression, we first note that This holds because the local stabilizers added to extend B T,loc to B Q,loc are contained only in Q and not in W , and vice versa. Next we consider the union B Q,loc ∪ B W,loc . This set spans G P,loc , because there is a basis for G P,loc where the basis elements are local operators that are thus fully contained either in Q or W . However it can happen that this set is not linearly independent. Therefore we have |B P,loc | = |B Q,loc ∪ B W,loc | − ∆ P QW T , where ∆ P QW T can be thought of as the number of constraint equations obeyed by elements of B Q,loc ∪ B W,loc . These constraints must be non-local, involving stabilizers contained in Q but not in W , and vice versa. It follows that S P QW T topo = −Ω P + Ω Q + Ω W − Ω T + ∆ P QW T . (B4) We now discuss ∆ P QW T in more detail; among the models we consider, it can be non-zero only for the d = 3 toric code. Each constraint equation contributing to ∆ P QW T is a product of local stabilizers in P that equals the identity operator. Such a constraint can always be obtained as a product of local constraints, involving stabilizers that may lie in a larger region containing P . This can occur for the plaquette stabilizers d = 3 toric code. For example, using the regions P, Q, W, T shown in Fig. 8, a product of plaquette stabilizers in P over a surface S enclosing the inner cube is the identity operator. The plaquette stabilizers satisfy the local constraint that a product of stabilizers over the faces of an elementary cube is unity, and the non-local constraint in P can be obtained by from these local constraints by taking a product over all cubes inside S.
However, there are no such non-local constraints in the other models we consider, where ∆ P QW T = 0. In the d = 2 toric code and Haah's code, the stabilizers obey no local constraints. This is also true for the cube stabilizers of the X-cube model. The vertex stabilizers of the X-cube model do obey local constraints, but because each vertex stabilizer participates in exactly one local constraint, it is not possible to obtain non-local constraints contributing to ∆ P QW T by taking products of local ones. Now we give a similar discussion of S ABC topo , to obtain the result that it is determined by stabilizers whose support is split among all of A, B and C. Here, we consider only models and regions where A, B and C are locally generated by a basis of local stabilizers that obey no local constraints. This property is established for the X-cube model in Appendix A 1 using a basis of cube stabilizers, and xy and xz plane vertex stabilizers. For Haah's code and the d = 2 toric code, the stabilizers obey no local constraints. We choose bases B A , B B and B C for the stabilizer groups B A , B B and B C , respectively.
To obtain a basis for the pairwise unions AB and so on, we consider the set b AB = B A ∪ B B . This set is clearly linearly independent, except possibly for the vertex stabilizers of the X-cube model. Any linear relation would have to involve basis stabilizers in both A and B. It cannot be a non-local constraint as discussed above. It also cannot be a local constraint near the boundary of A and B, because the stabilizers in B A and B B are drawn from a subset of local stabilizers that do not obey any local constraints. Therefore b AB is linearly independent. In order to obtain a basis B AB for G AB , we extend b AB by adding stabilizers δ AB whose support is split between A and B.
Finally we consider the union ABC. To find a basis, we first consider the set b ABC = B A ∪ B B ∪ B C ∪ δ AB ∪ δ BC ∪δ AC . Again, this set is linearly independent, and we extend it to a basis B ABC by adding a set of stabilizers δ ABC , whose support is split among the three regions.
We then obtain for the topological entanglement entropy the desired result.
Appendix C: Schrieffer-Wolff perturbation theory In this section we discuss the construction of dressed stabilizers upon perturbation of the stabilizer Hamiltonian. The derivation is a variation on the standard method of Schrieffer-Wolff transformations, most closely related to the method of Wegner-Wilson flow discussed in Ref. 44. We outline it here mainly in the interests of completeness.
Let the eigenvectors of the unperturbed Hamiltonian (H 0 ) be {|n } with eigenvalues {E n }. Let the new Hamiltonian be written as H = H 0 +λV , where V is some Hermitan operator and λ is any real number. Also let the eigenvectors of this Hamiltonian be {|n } with eigenvalues {E n }. The idea is to find the unitary operator U (λ) such that |n = U (λ) |n . (C1) If we assume that U (λ) is an analytic function, then there exists a Hermitian operator-valued function F (λ) such that U (λ) = exp (iF (λ)) .
We may then expand F (λ) as a power series in λ. To find a relation for U mn = m|n , note that m|H 0 + λV |n = E m m|n +λ m|V |n = E n m|n .
(C3) We may also expand the V term in the original eigenbasis, to obtain the self-consistent equation We now make the standard assumption that V does not act within degenerate subspaces of H 0 (any portion of V that does so act should be absorbed into our definition of H 0 ), so we only have to worry about 'off diagonal' matrix elements of V . If we define we can use (C4) to express the expansion coefficients as assuming the derivatives of A are well behaved at λ = 0. At first order, we have where L can be expressed in terms of its matrix elements in the old eigenbasis, or [H 0 , L] = iV . One can continue the expansion to any desired order. Note that if one is concerned only with the ground state, and the system is gapped within a topological sector, then the denominator has a non-zero lower bound. It may be shown using standard techniques that if the original Hamiltonian has a conserved quantity, ie. [H 0 , S] = 0, then one can find a conserved quantity for the perturbed Hamiltonian,S by solving The expert reader will recognize this as the equation of motion for Wegner-Wilson flow 44 . It is well known that this sort of flow equation preserves locality of the integrals of motion both for gapped systems and for localized systems (see e.g. Refs. 43 and 44 for recent discussions), where 'local' means 'local up to an exponentially decaying tail.'