Subsystem symmetry enriched topological order in three dimensions

We introduce a model of three-dimensional (3D) topological order enriched by planar subsystem symmetries. The model is constructed starting from the 3D toric code, whose ground state can be viewed as an equal-weight superposition of two-dimensional (2D) membrane coverings. We then decorate those membranes with 2D cluster states possessing symmetry-protected topological order under line-like subsystem symmetries. This endows the decorated model with planar subsystem symmetries under which the loop-like excitations of the toric code fractionalize, resulting in an extensive degeneracy per unit length of the excitation. We also show that the value of the topological entanglement entropy is larger than that of the toric code for certain bipartitions due to the subsystem symmetry enrichment. Our model can be obtained by gauging the global symmetry of a short-range entangled model which has symmetry-protected topological order coming from an interplay of global and subsystem symmetries. We study the non-trivial action of the symmetries on boundary of this model, uncovering a mixed boundary anomaly between global and subsystem symmetries. To further study this interplay, we consider gauging several different subgroups of the total symmetry. The resulting network of models, which includes models with fracton topological order, showcases more of the possible types of subsystem symmetry enrichment that can occur in 3D.


I. Introduction
A new paradigm in the classification of gapped phases of matter has recently begun thanks to the discovery of models with novel sub-dimensional physics.This includes the fracton topological phases [1][2][3][4][5][6][7][8][9][10][11][12], in which topological quasi-particles are either immobile or confined to move only within subsystem such as lines or planes, as well as models with subsystem symmetries, which are symmetries that act non-trivially only on rigid subsystems of the entire system.These two types of models are dual: in the same way that topological order can be obtained by gauging a global symmetry, gauging subsystem symmetries results in fracton topological order [7,8,13].Models with subsystem symmetries are also interesting in their own right.For example, one may use them to define symmetry-protected topological (SPT) phases, giving rise to the so-called subsystem SPT (SSPT) phases [14][15][16][17][18][19][20][21][22].SSPT phases display new physics compared to conventional SPT phases such as extensive edge degeneracy and unique entanglement properties.They also serve as resources for quantum computation via the paradigm of measurement-based quantum computation [14,21,[23][24][25][26][27].
In this paper, we build on this paradigm of sub-dimensional physics by considering how topological order may be enriched in the presence of subsystem symmetries.We call such order subsystem symmetry-enriched topological (SSET) order.For global symmetries, the theory of symmetry-enriched topological (SET) phases of matter is well developed, and is characterized by a non-trivial action of the symmetry on the topological excitations [28][29][30][31][32][33].In 2D, these excitations are point-like anyons, and the symmetry can act on the anyons by permuting them.The anyons may also carry a fractional symmetry charge, which means that the action of the symmetry on a state containing several anyons reduces to non-trivial projective actions localized around each anyon.For example, in spin liquids with SO(3) spin rotation symmetry, the whole system forms a spin singlet, but individual anyons, i.e. spinons, may have non-trivial spin, and hence carry a degeneracy that is protected by the symmetry.Crucially, single anyons cannot be created, and any physical state containing some anyons will transform under symmetry overall in a trivial manner, despite the non-trivial action on each anyon.In 3D, the situation is made more complicated by the presence of extended loop-like excitations, but significant progress has nonetheless been made [34][35][36][37][38][39][40].
In the presence of subsystem symmetries, the kind arXiv:2004.04181v1[cond-mat.str-el]8 Apr 2020 of fractionalization on mobile point-like excitations described in the preceding paragraph cannot occur.In anticipation of the main model of this paper, consider a 3D system with planar subsystem symmetries, and suppose there are two point-like excitations living on a single plane.It may seem as though it is possible for each excitation to carry a fractional charge, as before.However, in this case, we can simple move one of the excitations away from this plane, such that the subsystem symmetry now only acts on a single excitation.Since the total action on the system must be trivial, the local action on a single excitation must also be trivial, and there is no fractionalization or permutation.One way around this is to consider excitations with restricted mobility, as in fracton topological order.In this case, it may not be possible to move an excitation away from the symmetry plane, or there may be a conservation law that restricts the number of excitations on each plane to be e.g.even, so fractionalization again becomes possible.Such an option was explored in Ref. [16].
In this paper, we show that topological excitations of higher spatial dimension lead to another opportunity for subsystem symmetry fractionalization.Namely, consider a 1D loop-like excitation in a 3D system.A global symmetry cannot fractionalize on the loop, at least not in the simple way described above, since these is no way to decompose the symmetry action as a product of disjoint local actions.However, when we consider subsystem symmetries, fractionalization on a single loop is possible.Consider a plane which intersects a loop-like excitation.Acting with a symmetry on the plane reduces to a product of local actions at each intersection point.Since the number of intersection points is necessarily even, it is possible for the symmetry to fractionalize at each point.This is precisely what occurs in the model described in this paper.Additionally, we find models containing both point-like excitations with restricted mobility, as well as loop-like excitations, and a non-trivial symmetry action that couples the two types.
Our SSET model can be understood as a decorated 3D toric code model.[41,42].This toric code consists of qubits on the faces of a cubic lattice, and the ground states can be visualized as equal-weight superpositions over all basis states where the faces in the state |1 form unions of closed 2D membranes.To enrich this model with subsystem symmetries, we introduce new qubits on the edges of the lattice and couple them to the faces in such a way that edges lying along membranes form 2D cluster states, and those away from membranes remain in trivial product states.The 2D cluster state [25] has line-like subsystem symmetries and is the prototypical example of SSPT order [14,15], see Fig. 1(a).Constructed in this way, the decorated model has planar subsystem symmetries acting on the decorating qubits.This is due to the fact that, along the intersections of a given plane and the membranes, the symmetry action of this plane reduces to the line-like symmetries of the decorating cluster states, see Fig 1(b).Furthermore, the loop excitations of the toric code, which appear along the boundary of open membranes, now coincide with the boundary of cluster states.Due to the SSPT order of cluster states, these boundaries, and hence the loop excitations, fractionalize under the planar symmetries.The immediate consequence of this fractionalization is that the loop excitations are endowed with a degeneracy per unit length which is protected by the subsystem symmetry.The procedure of decorating topological orders with lower-dimensional SPT orders is a well-established way to create SET orders [43][44][45][46], and our construction here can be seen as a generalization of this procedure to subsystem symmetries, with the notable feature that d-dimensional subsystem symmetries of the decorating SSPT translate into (d + 1)-dimensional subsystem symmetries of the decorated model.
Another effect of the subsystem symmetry fractionalization that we discover is an increased value of the topological entanglement entropy (TEE).The TEE refers to a constant correction to the area law when computing bipartite entanglement entropy, and is a topological invariant, in the sense that it takes a uniform value within a given topological phase of matter [47,48].Recently, it has been understood that there are certain topological states for which the TEE, for certain bipartitions, does not match the expected value [17,19,22,[49][50][51][52][53].In a majority of known cases (see Ref. [53] for a possible counterexample), this is due to the presence of lower SSPT Sec.II

SE Fracton
Sec. IV B 3 Flowchart describing the various models obtained by gauging and ungauging the symmetries of the SSET in different ways.The arrows are labelled by the symmetry which is gauged along the direction they point, where g, s1, and s2 represent the global symmetry, lattice-plane subsystem symmetries, and dual-plane subsystem symmetries, respectively.A model is called "fracton" if all excitations display some mobility restrictions, and "panoptic" if restricted-mobility excitations appear alongside fully mobile and loop-like excitations."SE" denotes models with symmetry enrichment due to global symmetries alone, while "SSE" denotes models enriched by subsystem symmetries, possibly alongside global symmetries.
dimensional SPT order around the boundary of the bipartition, which can in turn be related to the presence of SSPT order [17].In fact, it was shown in Ref. [22] that this non-zero TEE takes a uniform value within an SSPT phase, and can be used to detect SSPT order and phase transitions.To emphasize its relation to the SSPT order, it was dubbed the symmetry-protected entanglement entropy (SPEE).For the case of an SSET, which naturally has a non-zero TEE, one might therefore expect a larger value of the TEE due to the additional presence of the SPEE.This is exactly what we confirm in our SSET model.Our SSET model stems from a unique interplay of global and subsystem symmetries.To better understand this interplay, we consider gauging and ungauging the symmetries of the SSET in various combinations, resulting in a network of eight different models, as pictured in Fig. 2. At the root of this network is a short range entangled model with SSPT order under a combination of global and subsystem symmetries.We calculate the cocycle that encodes the non-trivial action of the symmetries on the boundary of this model, revealing a mixed anomaly between global and subsystem symmetries.Using this cocycle, we calculate the effect of the symmetries on the extrinsic symmetry defects, which allows us to predict the outcome of gauging the symmetries in different combinations.We then discuss the nature of the topo-logical order and the symmetry enrichment of each of the eight models, and in particular show that the model resulting from gauging all symmetries, i.e. the model resulting from gauging the subsystem symmetries of the SSET, displays the "panoptic" order that was recently discovered in Refs.[54,55], and contains non-abelian fractons.
The rest of the paper is structured as follows.In Section II, we introduce a model with both global and subsystem symmetries, and show that it has non-trivial SPT order by analysing its boundary.Then, in Section III, we gauge the global symmetries of this model to obtain our model of SSET order and study its properties.In Section IV, we consider gauging the subsystem symmetries and investigate the network of models in Fig. 2. Finally, in Section V, we discuss some principles that could help guide a general theory of subsystem symmetry enrichment, give an argument against the existence of SSET order in 2D, and discuss possible routes of future work.

II. Symmetry protected topological order with global and subsystem symmetries
To derive our model of SSET order, we begin with a short-ranged entangled model that has SPT order with respect to a combination of global and subsystem symmetries.The SSET is obtained by gauging only the global symmetries of this model.The process of partially gauging symmetries of an SPT order to obtain SET order is well understood in the case of global symmetries [32,40,56,57], and is reviewed in a simple 2D example in Appendix A.
We first define the short-ranged entangled model, and demonstrate its non-trivial SPT order by identifying the non-trivial action of the symmetries on the boundary, as encoded by a certain 3-cocycle.Notably, the non-trivial order arises from an interplay between the global and subsystem symmetries, and the system becomes trivial if only the global symmetry or only the subsystem subsystem symmetry is preserved.We will refer to this type of order as SSPT order, despite the equal importance of both global and subsystem symmetries.
The SSPT model lives on a simple 3D cubic lattice, with qubits in the body centers (C) and on the edges (E).To begin, consider a trivial paramagnetic Hamiltonian acting on this system, whose unique ground state is a product state of |+ = 1 √

(|0 + |1
) on every edge and body qubit.Our model can be defined by acting on this trivial system with a finite depth unitary circuit, where the product runs over all triples of qubits consisting of one body qubit and two of its nearestneighbouring edge qubits, as pictured in Fig. 3(a), and CCZ acts on the three qubits as CCZ|i |j |k = (−1) ijk |i |j |k .This gives the Hamiltonian, where C e is as depicted in Fig. 3(a), and where f ∈ c runs over the six faces of c, U f is a product of four CZ operators in a diamond on face f , and CZ|i |j = (−1) ij |i |j , see Fig. 3(a).
H SSP T has a unique ground state which we denote by |SSP T .We can get some intuition for this ground state using the viewpoint of decorated domain walls (DDW) [43], as shown in Fig. 3 where we have defined the state |G C on the edge qubits as, Since the CZ gates cancel out between neighboring cubes in C, |G C describes a state of 2D cluster states on the domain walls of C, and |+ states away from them.The 2D cluster state is the paradigmatic example of a state with SSPT order under line-like symmetries [14,15].Therefore, |SSP T can be described as an equal weight superposition over all body qubit configurations in which the domain walls of the body qubits are decorated with 2D SSPT states.Let us turn to the symmetries of |SSP T .First, by nature of the DDW structure of |SSP T , we have a global Z 2 symmetry acting on the body qubits X C = c∈C X c , which follows simply from the fact that the domain walls are invariant under flipping all body spins.The edge qubits have planar subsystem symmetries.For any plane moving parallel to one of the coordinate planes of the cubic lattice, we define the subset P ⊂ E as the set of edges that are intersected by this plane.We remark that these planes come in two distinct types, determined by whether they are parallel or perpendicular to the edges they intersect, and we refer to the two types as latticeplanes and dual-planes, respectively.We then define the subsystem symmetry for each plane as X P = e∈P X e .The fact that X P is a symmetry of |SSP T for all planes P can be seen by first noticing that [U c , X P ] = 0 for all c ∈ C and all P. Then it follows that X P is a symmetry of |G C for all C ⊂ C, and therefore is also a symmetry of |SSP T thanks to Eq. ( 5).
There is a more insightful way to understand the presence of the planar symmetries.Observe that the intersection of a plane P with a domain wall configuration forms closed 1D loops, as pictured in Fig. 1(a).Since X P acts trivially away from the domain walls, we can restrict the action of X P onto these closed loops.This is a symmetry, since the 2D cluster states living on the domain walls have line-like subsystem symmetries [58].Therefore, the planar symmetries of the 3D SSPT follow from the linelike symmetries of the 2D SSPT states which we use to decorate domain walls.
We can see explicitly that the SPT order of |SSP T is trivial if only the global or subsystem symmetries is enforced by constructing disentangling circuits which respect one of the symmetries.Namely, if we group all of the CCZ gates from U CCZ that live in a given cube, we obtain a local unitary that respects all of the subsystem symmetries.Applying this unitary to all cubes is therefore a subsystem symmetry-respecting circuit that disentangles |SSP T .A circuit which respects the global symmetries can be obtained by similarly grouping CCZs into octahedrons, one for each face of the lattice.Importantly, neither circuit respects both types of symmetry, and we argue in the next section that such a symmetric circuit does not exist by virtue of the non-trivial SSPT order.

A. Boundary of the SSPT
We now demonstrate the non-trivial nature of our SSPT model by analyzing its boundary.Let us consider the geometry pictured in Fig. 4, where the boundary is a 2D square lattice with qubits on the edges and periodic boundary conditions.In the presence of this boundary, the whole 3D state may no longer be symmetric under the global or subsystem symmetries.In particular, it may be necessary to dress symmetry operators with additional action on the boundary qubits in order to leave the system invariant.For the planar symmetries, this turns out to be unnecessary, and the action of the planar symmetries on the boundary, for planes perpendicular to the boundary, is simply a line of X's.The global symmetry, on the other hand, must be decorated with additional CZ's acting between every nearest neighboring pair of edges on the boundary.The action of the symmetries on the boundary is summarized in Fig. 4.
We see from the above analysis that neighbouring plane symmetries commute on the boundary.This is different from the 2D SSPT order of the 2D cluster state, where neighbouring line symmetries commute in the bulk but anti-commute on the boundary.However, we will see that, in the presence of a global symmetry flux threading the cylinder, which is introduced by adding the domain wall in Fig. 4, neighbouring plane symmetries will anticommute.This highlights the importance of both global and subsystem symmetries in our SSPT model.
To formalize this, we can use the language of group cohomology which classifies SPT order.Imagine "compactifying" our 3D system into a quasi-2D system [59].This is achieved simply by fixing the width R in the y direction, as indicated in Fig. 4. We then get a different quasi-2D system for each compactification radius R. The subsystem symmetries perpendicular to the compactification direction, i.e. xz planes, become standard global symmetries of the compactified system.There are 2R of these symmetries, labelled by indicies for i = 1, . . ., 2R where even (odd) i corresponds to lattice-plane (dualplane) symmetries, as pictured in Fig. 4. Including the global symmetry as well, we can consider the symmetry group, where the generators g and s i correspond to the global and subsystem symmetries, respectively.
For each R, we can study the 2D SPT order of the compactified system under the symmetry group G R .This order is determined by a ) which characterizes the action of G R on the boundary [60,61].We can straightforwardly calculate this cocycle using the dimensional reduction procedure introduced in Ref. [62].The result is, We can determine the effect of inserting a global symmetry flux using the slant product, as described in Refs.[61,63].The slant product corresponding to a group element a is a function When ω is a 3-cocycle, χ a will be a 2-cocycle, corresponding to a cohomology class ).In particular, if two 3-cocycles belong to the same cohomology class, then their slant products do as well.The physical meaning of the slant product is the following.
For any a ∈ G, let V a (b) be the action of b ∈ G on the boundary in the presence of an a-flux.This action may be projective, i.e.V a (b)V a (c) may equal V a (bc) only up to a phase.This phase is precisely given by the slant product, V a (b)V a (c) = χ a (b, c)V a (bc) .In the case that χ a belongs to a non-trivial cohomology class, the a-flux carries a non-trivial 1D SPT order, as characterized by the 2-cocycle χ a .We remark that, for finite abelian G, χ a is in a trivial cohomology class if and only , the cohomology of χ a is trivial if and only if the action of the symmetry on the boundary commutes.
If we compute the slant product with a = (1, 0) corresponding to the global symmetry, we find, In particular, the commutation relation reads, This says precisely that neighbouring plane symmetries anti-commute on the boundary in the presence of a global symmetry flux.We can also compute the slant product for a subsystem symmetry flux, corresponding to a = (0, î), where the vector î contains 1 at position i, and 0 elsewhere.We find, which tells us that, in the presence of a lattice-plane (dual-plane) subsystem symmetry flux, the global symmetry anti-commutes on the boundary with the dualplane (lattice-plane) subsystem symmetries neighbouring the flux.
One might wonder if we are missing out on any important information by only considering subsystem symmetries in one direction (the xz planar symmetries).Indeed, it could be the case that, e.g.perpendicular symmetry planes anti-commute in the presence of a certain symmetry flux.However, we have checked that this is not the case: the cases considered above capture all of the non-trivial fractionalization that occurs in our model.
From the non-trivial slant products computed above, we see that the 3-cocycle ω belongs to a non-trivial cohomology class for all compactification radii R. The general arguments of Ref. [60] then show that the boundary, when considered as a quasi-1D system, cannot be gapped and symmetric; either the symmetry will be spontaneously broken, resulting in a boundary degeneracy, or the boundary will be gapless.However, it may be that, when we enforce a true 2D notion of locality on the boundary, it gains even more non-trivial features that we miss by employing compactification.In Appendix B, we examine some possible Hamiltonians which respect the boundary symmetries.We find that, unlike in conventional SPT phases, neither of the simplest choices of boundary Hamiltonian lead to a gapless boundary, suggesting that it is necessary to look beyond this compactified description to fully understand this boundary.We leave such an analysis to future work.

III. Gauging the global symmetry: Subsystem symmetry enriched topological order
In this section, we gauge the global Z 2 symmetry of |SSP T .The resulting model is a Z 2 gauge theory in which the loop-like topological excitations fractionalize under the subsystem symmetries.We therefore call this model an example of SSET order.We show that this implies an extensive degeneracy of the loop excitations which is protected by the symmetry.We also show that the model has an enlarged value of the topological entanglement entropy, as compared to the underlying topological order, due to the subsystem symmetry enrichment.For simplicity, we assume for the time being that our system is defined on a topologically trivial 3D manifold, such that there will be no topological ground state degeneracy.
Let us briefly describe the gauging procedure we which we employ (more details may be found in Refs.[7,8,13]).The gauging procedure maps body qubits to qubits on the faces (F ) of the lattice, in such a way that the face qubits f ∈ F take the state |1 on the domain walls of the body qubits, and |0 elsewhere.That is, given a state |C on the body qubits, the effect of the gauging map Γ can be written as Γ|C = |∂C , where ∂C ⊂ F is the set of faces on the boundary of C, and |∂C describes a state on the face qubits where all qubits in ∂C are in the state |1 , and the rest are in |0 .We can then extend the map Γ to arbitrary states on the body qubits by linearity.Applying this procedure to |SSP T , we get, This state can again be visualized using Fig. 3(b), where now the colored faces indicate elements of ∂C, which can be viewed as a configuration of closed membranes on the faces of the lattice.Therefore, |SSET is a superposition over all closed membranes configurations on the face qubits, where the membranes are decorated with 2D cluster states on the edge qubits.It is clear that |SSET has the same subsystem symmetries on the edge qubits as |SSP T .
If we were to remove |G C from the above equation, the state described would be exactly the 3D toric code, |ψ T C = C⊂C |∂C , which has topological order.In fact, we can disentangle the edge qubits from the face qubits using a unitary circuit U CCZ which places four CCZ's on each face as in Fig. 5(a).Applying this circuit to |SSET , we get, Therefore, |SSET is related to |ψ T C by a unitary circuit, up to trivial degrees of freedom, so the two states have the same topological order.However, this circuit does not respect the subsystem symmetries.In fact, we will show that, when these symmetries are enforced, |SSET is in a distinct phase from |ψ T C , as indicated by symmetry enrichment of the topological excitations.Therefore, we can say that |SSET has subsystem symmetry enriched topological order.

A. Excitations and symmetry enrichment
The Hamiltonian for |SSET can be obtained by gauging H SSP T , and has the following form, Excitations appear on edges marked by points.(c) A 2D cut of our 3D lattice with a loop excitation corresponding to the half-infinite membrane Σ ∞ , as viewed from above.There is an effective qubit degree of freedom for each edge on the boundary of Σ ∞ , indicated by dots.The horizontal dashed lines represent the intersection of two symmetry planes P with this 2D slice, and the action of XP on the edge degrees of freedom, according to Eq. 27, is shown beside each line. where, where f e denotes all faces incident on edge e. C e is pictured in Fig. 5(a).The new gauge term A e enforces that the faces qubits form closed membranes.The modified terms B c and C e are simply the original terms B c and C e rewritten as functions of the new face qubits, using the rules X c → f ∈c X f and c f Z c → Z f that follow from the gauging procedure.In addition, we project the C e term onto the closed membrane subspace (A e = 1 ∀e ∈ E) in order to ensure that the Hamiltonian respects the subsystem symmetries.We note that all terms in H SSET commute.
Let us now construct the excitations of H SSET , which come in three types.Violations of C e are topologically trivial particles that can be created locally by acting with Z e , so we do not discuss them further.The other two types of excitation are topologically non-trivial, and must be created by extended non-local operators.Violations of B c , called electric excitations, are point-like, and are created in pairs at the end points of string operators, where Λ ⊂ F is a curve on the dual lattice.Violations of A e , called magnetic excitations, are loop-like.They appear along the boundary of open membrane operators, where Σ ⊂ F is an open membrane of faces, as pictured in Fig. 5(b).If we were to braid an electric excitation through the loop of a magnetic excitation, we would pick up a minus sign due to the anti-commutation of X and Z.This is the same as in the toric code.
The difference from the toric code comes from the degeneracy of the excitations.In particular, we can decorate the boundary of the membrane operator S m Σ with Z operators without changing the energy of the resulting excitation: Z ae e , a e = 0, 1, where ∂Σ denotes all edges on the boundary of the membrane Σ.Thus the loop-like excitation is two-fold degenerate per unit length.This degeneracy is protected by the subsystem symmetries.Intuitively, this is because the degeneracy originates from the SSPT order of the 2D cluster states that decorate the membranes, which have an exponential edge degeneracy.The edge degeneracy of the 2D cluster states is protected by the line-like subsystem symmetries.Since the planar symmetries of the 3D lattice act as line-like symmetries on the 2D membranes (See Fig. 1), they in turn protect the degeneracy of the magnetic excitations.More precisely, consider a membrane Σ ∞ which corresponds to a half-infinite plane, as shown in Fig. 5(c).This creates an excitation lying along ∂Σ ∞ which corresponds to a line of N edges denoted e i for i = 1, . . .N .The degenerate subspace associated to this excitation is spanned by the states |a 1 , . . .a n defined as, Let us now determine the action of the subsystem symmetries in this N -qubit space.Consider those symmetry planes that are perpendicular to Σ ∞ and also cross ∂Σ ∞ , as pictured in Fig. 5(c).We can index these planes by P i and P i+1/2 , corresponding to dual-planes intersecting edge e i , or lattice-planes intersecting between edges e i and e i+1 , respectively.Then we can calculate, Therefore, if we let X i , Z i denote the logical Pauli operators in the degenerate edge subspace, we have X Pi ∼ = Z i and X P i+1/2 ∼ = X i X i+1 .We see that neighboring plane symmetries anticommute on the edge of the excitation.This is the same pattern of symmetry fractionalization found on the boundary of the 2D cluster state [15].In that case, the symmetry protects the exponential edge degeneracy of the 2D cluster state.In analogy, the symmetry here protects the exponential degeneracy of the line-like excitation.A similar discussion holds for an arbitrary membrane Σ, although there can be some finite size effects due to corners, as discussed in Ref. [15].Therefore, the loop-like magnetic excitations of |ψ 3D SSET carry a two-fold degeneracy per unit length that is protected by the planar subsystem symmetries.

B. Topological entanglement entropy
In this section, we show that the topological entanglement entropy (TEE) of our 3D SSET state is larger than that of the 3D toric code when the boundaries of bipartitions are aligned with the symmetry planes.This is in analogy to the fact that states with 2D SSPT order, such as the cluster state, have a non-zero TEE when biparitions are aligned with the line-like symmetries, despite the absence of topological order [17,22,50,52].To emphasize that this value comes from the SSPT order, rather than topological order, we call it the symmetryprotected entanglement entropy (SPEE) [22].Therefore, in this case, we say that the 3D SSET has a non-zero SPEE in addition to the TEE inherited from the topological order of the 3D toric code.
The origin of the SPEE can be imagined as follows.Consider one closed membrane in the membrane soup defining |SSET .If we bisect the membrane with a plane aligned with the subsystem symmetries, as in Fig. 1, we can calculate the entanglement between the two halves of the membrane.Since the membranes are decorated by 2D cluster states, and the cut is aligned with the linelike symmetries of the cluster state, we find a non-zero SPEE.Since |SSET is a fluctuating soup of such membranes, one can imagine that it also exhibits a non-zero SPEE.In the rest of this section, we show that this is indeed the case, with details of the calculation presented in Appendix C.
For simplicity we calculate the entropy of a finite section of a 3D torus, such that the boundary between subsystems A and B is two disconnected 2D tori.However, we expect the same results to hold for any geometry that is appropriately aligned with the subsystem symmetries, as in Ref. [52].We will aim to determine the 2-Rényi entropy S (2) A ) where ρ A is the reduced state of subsystem A [65].On a 3D torus, H SSET has eight degenerate ground states.We choose |SSET to be one of the minimal entropy states which have the largest TEE [66].This is done by picking the ground state that is +1 eigenstate of the membrane operators S m Σz where Σ z is the non-contractible membrane , as well as the loop operators S e Λ x/y where Λ x/y are the non-contractible loops, as pictured in Fig. 6(a).In this way, the subsystem A has maximum knowledge of the whole state, and hence minimum entropy.We define G to be the abelian group generated by all Hamiltonian terms A e , B c , C e as well as Σz and S e Λ x/y .For simplicity, we assume that region A contains L × L × L vertices, such that L 2 edges are cut on each boundary.In Appendix C, we show that the entropy can be expressed in the following way for large L, S where |A| = 6L 3 + 3L 2 is the number of qubits in A, G A is the subgroup of G containing those operators that act non-trivially on A only, and F(β) is the (extensive) free energy of a 2D square lattice Ising model at inverse temperature β.By counting independent generators in the same way as for the 3D toric code, [42], we can compute that . Compared to the toric code case, G A here contains an additional non-local operator for each boundary, namely the product of C e for all edges e on the boundary, which is equal to the subsystem symmetry on the corresponding dual plane.
To compute the free energy, we can use Onsager's result [67], which yields F(ln √ 2) = ln 2( 1 2 + ln 2)L 2 for large L (for more details, see Appendix C).Crucially, the temperature β = ln √ 2 lies in the disordered phase of the Ising model, such that the free energy is extensive with no constant term.Putting everything together, the entropy is, S where c = (3 − ln 2) ln 2, γ TEE = γ SPEE = ln 2, and the dots represent terms that go to zero as L goes to infinity.We include the factors of 2 to emphasize the fact that A has two disconnected boundaries.γ TEE comes simply from |G A | in the same way as for the 3D toric code, whereas γ SPEE is due to the subsystem symmetries forming non-local constraints on the boundary of A.
While we have only shown the existence of a non-zero SPEE for this specific bipartition, we expect that it would also be present for other bipartitions whose geometry are aligned with the symmetry planes, as is the case for the SPEE of the 2D cluster state [52].Also, in Ref. [52], the authors constructed a method to extract the SPEE due to line-like subsystem symmetries, such that it can be separated from the TEE.They show that the combination of entropies is equal to the SPEE, where subsystems A, B, C form a dumbbell shape as in Fig. 6(b).This is due to the fact that the line-like symmetries (blue line in Fig. 6(b)) can be terminated by applying local operators in the circled regions, such that we get a constraint that reduces S ABC .Crucially, the red regions have a larger radius than the width of B, such that the same truncated line cannot fit into B, AB, or BC.On the other hand, the TEE is the same for each of the four terms in Eq. ( 30), and hence will cancel out, as will all extensive parts of the entropy.Following the same logic, we propose a set of subsystems A, B, C such that S SPEE will be equal to the SPEE due to planar symmetries, see Fig. 6(c).The subsystem B is a thin slab, while A and C are thickened, giving a geometry similar to that of a picture frame.As in the dumbbell, the planar subsystem symmetries can be truncated by applying operators along the boundary.Such an operator will fit into region ABC of the frame, but none of the other combinations.Therefore, we conjecture that S SPEE should be equal to γ SPEE for the picture frame geometry.If there is also a SPEE due to linear subsystem symmetries, then S SPEE would grow as the perimeter of region B. If we take the alternate picture frame geometry in Fig. 6(d), we would only detect SPEE due to planar symmetries, since any truncated line operator that contributes to S ABC will also contribute to one of S AB or S BC .

IV. Gauging the subsystem symmetries
In this section, we examine the variety of topological phases that arise from gauging some or all of the subsystem symmetries, starting either from |SSP T or |SSET .More precisely, starting from |SSP T , we can independently choose to gauge the global symmetry, lattice-plane symmetry, and dual-plane symmetry, resulting in eight possible models (including |SSP T ).We denote the global, lattice-plane, and dual-plane symmetries as Z glob , respectively.We use this notation with the understanding that there are a subextensive number of subsystem symmetry generators, so Z sub 1/2 2 describes the local action of the symmetry, not the total symmetry group.It is known that gauging subsystem symmetries can result in models with fracton topological order [7,8,13].Here, by "fracton" topological order (or more succinctly fracton order), we mean any model in which all topological excitations have restricted mobility of some sort.We use the usual terminology of fracton, lineon, and planeon to describe particles which are fully immobile, constrained to a 1D line, or constrained to a 2D plane, respectively.Later, we consider models in which such particles coexist with fully-mobile excitations, and we refer to the order in these models as "panoptic" [54].
Throughout this section, we do not explicitly perform the gauging as we did in the previous section (except in one case).Rather, we rely on our understanding of |SSP T and its symmetries to determine both the mobility and type of symmetry enrichment of the resulting topological excitations.This is done by examining the symmetry defects of |SSP T , whose relation to gauging we now describe.

A. Symmetry defects and gauging
To define a symmetry defect, we first define domain wall operators.Consider the operator U R which applies some symmetry operator to a compact region R of the lattice.For a global symmetry, R is some 3D region, whereas for a subsystem symmetry, R would be a 2D region confined to a plane.Then, for Hamiltonians composed of local, symmetric terms, the only non-trivial effect of U R will be near the boundary ∂R of R. Therefore, we can write U R ∼ = V ∂R , where ' ∼ =' indicates that the two operators act in the same way within the ground state subspace, and V ∂R is an operator acting on ∂R which we call the domain wall operator.V ∂R is precisely a membrane of CZ's for Z glob 2 , while for the planar subsystem symmetries it is a 1D loop operator also consisting of CZ's., on the other hand, are point-like and can only move within a single plane without creating additional excitations, i.e. they are planons.This is because the symmetry defects themselves are mobile only within a given plane.However, we find that in some cases these planons can be decomposed into a pair of excitations of lower mobility, as in the X-cube model [7].This is due to the symmetry domain walls decomposing further, such as a planar domain wall decomposing into a product of two cageedge domain walls as in X-cube [10].The properties of the symmetry defects before gauging, Including the action of symmetry on them, determine the braiding and fusion statistics of the gauge fluxes, as well as possible symmetry fractionalization under any ungauged symmetries [32,56,57].
The other type of topological excitations that emerge from gauging are called the gauge charges.The pure gauge charges are gauged versions of symmetry charges, which are objects that locally anti-commute with the symmetry, corresponding in our case to a single Z operator on the lattice.There are also composites of gauge charge and existing topological excitations or fluxes, sometimes referred to as dyons.These arise from the symmetry charges created by the action of symmetry on topological excitations or defects that are then promoted to gauge charges via gauging.Even for an abelian symmetry group this can lead to nonabelian topological excitations after gauging due to permutation or fractionalization [32].As described in Ref. [13], the mobility of gauge charges can be determined directly from the spatial structure of the subsystem symmetries.Namely, if a symmetry charge is acted on by planar symmetries in one, two, or three orthogonal directions, then the corresponding gauge charge will be a planon, lineon, or fracton, respectively.Relations between the symmetry generators that imply parity conservation for the number of subsystem charges on all planes in orthogonal directions can be used to determine whether such lineons and fractons are irreducible, meaning they are not a composite of particles of higher mobility [68][69][70].
The slant products calculated in Section II can be used to determine the relationship between the symmetry generators, symmetry charges, and symmetry fluxes.Indeed, although we initially performed these calculations to understand the action of the symmetry on the boundary in the presence of symmetry flux, a symmetry flux is in fact inserted using a domain wall operator, which terminates on the boundary via a symmetry defect.So we were actually calculating the action of the symmetry on the defect, regardless of whether or not it is pushed to the boundary.The algebra of symmetry actions on defects in |SSP T , as described by the slant products, allows us to calculate the stable topological excitations when a (subsystem) symmetry is gauged, by condensing domain walls and projecting topological excitations and defects onto irreps of the appropriate symmetry actions.We can further calculate the action of the remaining ungauged symmetry on the topological excitations via the algebra formed by the remaining symmetry operators and the fluxes and charges of the gauged symmetries.
More precisely, for the SSPT model, if symmetries b and c anticommute on a a-defect, i.e. χ a (b, c)/χ a (c, b) = −1, then this will manifest either as a fractionalization or permutation of the excitations in the gauged theory, depending on which symmetries are gauged.Specifically, if a is among the gauged symmetries, while b and c are not, the a gauge flux will carry a fractional charge under b and c.If a and b are gauged and c is not, then acting with c symmetry on an a (b) gauge flux attaches a b (a) gauge charge.For example, Eq. 10 says that acting with a Z sub1 2 symmetry generator on a Z glob 2 symmetry flux attaches Z sub2 2 symmetry charges to it on the neighbouring planes.This leads to the fractionalization observed in Section III.For a subsystem symmetry flux, Eq. 12 says that, if we act with a Z glob

B. Gauging only subsystem symmetries:
Fracton order We now turn to gauging the subsystem symmetries of |SSP T .The resulting models all have fracton topological order, characterized by topological excitations with restricted mobility.The ungauged symmetries will either permute the excitations or fractionalize on them, depending on whether we gauge all or only some of the subsystem symmetries.

Gauging lattice-plane symmetries
Gauging Z sub1 2 results in a non-trivial fracton model that is not obviously equivalent to any known models.Symmetry charges before gauging are acted on by two perpendicular subsystem symmetries, so the gauge charges are lineons.These are irreducible as the product of all lattice-plane symmetries is the identity, so there is a charge parity conservation symmetry for all lattice-plane symmetries taken together.Interestingly, the gauge fluxes transform as a non-trivial projective representation under Z glob 2 and adjacent Z sub2 2 symmetry generators.Therefore, this model represents a fracton model where excitations carry a fractional charge under the combination of global and subsystem symmetries.

Gauging dual-plane symmetries
Each edge qubit is acted on by only one generator of Z sub2

2
. Therefore, the gauge charges are planons, and furthermore we can gauge the symmetry in each dual-plane separately, such that we arrive at a stack of decoupled 2D toric codes in all three directions, one for each dual plane.The layers are coupled via the body qubits in such a way that the Z sub2 2 gauge fluxes in a given plane transform as a non-trivial projective representation under the adjacent Z sub1 2 generators and Z glob 2 .

Gauging all subsystem symmetries
Gauging all subsystem symmetries of |SSP T results in a model with fracton excitations, where the global symmetry permutes the fractons.This permutation action can be understood by the fact that the global symmetry, when acting on a lattice-plane (dual-plane) symmetry defect, attaches a pair of symmetry charges to the neighbouring dual-planes (lattice-planes).Hence, in the gauged model, the global symmetry attaches gauge charges to gauge fluxes.As in the general case, the gauge flux is a planon.However, in this case, a planar domain wall can be decomposed into a pair of domain walls and thus the planon can be split into a pair of fractons.We can isolate these two fractons by considering the symmetry defect corresponding to a stack of dual-plane symmetries.Gauging this defect gives two fractons, one at the top of the stack and one at the bottom.The fact that these are fractons follows from relations between the gauge constraints, as we elaborate on in Fig. 7. Applying the global symmetry attaches a gauge charge to each fracton.Since the symmetry charges are acted on by planar symmetries in all three directions, the gauge charges are fractons.In particular, there are relations between all subsystem symmetries in two orthogonal directions that imply the fractons are irreducible.Therefore, the global symmetry action causes fracton permutation.
This permutation action will be useful for understand the model obtained by gauging all symmetries in Section IV C 3, so here we will explicitly write down the gauged Hamiltonian.To perform this gauging, we follow the procedure for gauging subsystem symmetries outlined in Ref. [13].First, we identify the minimal coupling terms which commute with all subsystem symmetries, which in the present case takes the form of four-body interactions near an edge, pictured in Fig. 7.The gauge qubits are placed in the nexus of these interactions as shown in Fig. 7.We therefore have four qubits for each edge e, one for each cube c surrounding e.We may therefore label gauge qubits uniquely by a pair (c, e).We can view the gauge qubits as living on the edges of a smaller cube within each cube of the lattice, and we use the labelling in Fig. 7 to reference each of the twelve gauge qubits in a given cube.Following Ref. [13], we now express existing Hamiltonian terms B c and C e in terms of the new gauge qubits, resulting in the terms B c and C e .Furthermore, we add additional flux terms which enforce a zero-flux constraint on the gauge qubits.There is one flux term A e for each edge, and nine flux terms for each cube c.We define C e and the flux terms in Fig. 7.The precise form of B c is not required for our purposes here, as it is not associated with any topological excitations.Overall, the gauged Hamiltonian is,  Now let us investigate the excitations of H f rac .Due to the sheer number of terms, and the fact that they are not all independent, the full spectrum of excitations is tedious to describe.Instead, let us focus on two types which allow us to see the non-trivial action of the global symmetry.Gauge charges, i.e violations of C e , can be created at the corners of a membrane operator, where R ⊂ C is a rectangle of cubes in the xz plane, and excitations appear on the edges at the four corners of R. Gauge fluxes, corresponding to violations of A e , can be created by the following membrane operator, where gauge fluxes again appear at the four corners of R. The fact that these are fractons follows from the three orthogonal relations among flux terms described in Fig. 7. Now we observe that the global symmetry X C permutes the membrane operators, This shows that acting with the global symmetry attaches gauge charges onto gauge fluxes, thereby permuting the fractons of the gauged theory.We remark that Eq. 34 is implied directly by a similar equation describing the action of the global symmetry on the symmetry defect of a stack of dual-planes, and that similar equations hold for all permutation actions described in the following section.
C. Gauging global and subsystem symmetries: Panoptic order In this section, we consider gauging the global symmetry of |SSP T along with some or all of the subsystem symmetries.Equivalently, we are gauging the subsystem symmetries of |SSET .In each case, we will find that fully mobile point-like and loop-like excitations, coming from the gauged global symmetry, coexist with the restricted-mobility excitations coming from the gauged subsystem symmetries.Such a system was dubbed to have "panoptic" order in Ref. [54].

Gauging lattice-plane and global symmetries
If we gauge Z sub1 gauge fluxes in adjacent planes.Thus, the symmetry enrichment manifests as an interesting permutation involving fully mobile excitations and those of restricted mobility.

Gauging dual-plane and global symmetries
Gauging Z sub1

Gauging all symmetries
To understand the model obtained by gauging all symmetries, it is easiest to start from H f rac in Eq.31, where all subsystem symmetries have been gauged.Then, what remains is to gauge the global symmetry, which we saw enacts a non-trivial permutation on the fractons of the model.For this, we can use the general arguments of Ref. [54] (see also Ref. [55]).Therein, it is argued that gauging a fracton permuting symmetry results in a model with non-abelian fractons.Furthermore, gauge fluxes will be loop-like excitations that braid non-trivially with the excitations of reduced mobility.Since we expect the gauging of different symmetries to commute (as is indeed the case in Appendix A), we can conclude that gauging the subsystem symmetries of the SSET will result in the same panoptic order with non-abelian fractons.
It is interesting to compare to the model obtained by gauging the layer-swap symmetry of the bilayer X-Cube model, as presented in Refs.[54,55], which also has panoptic order with non-abelian fractons.The potential equivalence between this model and our own is suggested by the equivalence between the model obtained by gauging the layer-swap symmetry of the bilayer toric code [54,55], and the model obtained by gauging all symmetries of the non-trivial Z 2 × Z 2 × Z 2 SPT with type-III cocycle [71] (which is somewhat analogous to |SSP T ), as discussed in Appendix A.

V. Discussion & Conclusions
We investigated the phenomenon of subsystem symmetry enrichment in 3D systems.We began with a base model possessing SPT order under a mix of global and subsystem symmetries.By gauging the global symmetries of this model, we obtained a topological model with loop-like excitations that carry fractional charge of the subsystem symmetries, which we called an example of SSET order.We showed that this fractionalization leads to a extensive degeneracy of the excitations, as well as an increased value of the topological entanglement entropy.We then considered also gauging the subsystem symmetries of the base model, resulting in a network of models all related by gauging and ungauging symmetries (Fig. 2).Using the algebra of the symmetry defects of the SPT model, we were able to understand the nature of each gauged model, uncovering several distinct types of subsystem symmetry enrichment.In particular, we found models supporting subsystem symmetry fractionalization and permutation between mobile and restricted mobility excitations, as well as a model with non-abelian fractons.
To conclude, we present the first steps towards a general theory of subsystem symmetry enrichment, which allow us to argue against the existence of nontrivial subsystem symmetry enrichment in 2D systems.We then give an outlook on generalizations and possible applications of our results.

A. Towards a general theory of subsystem symmetry enrichment
As discussed briefly in the introduction, symmetry enrichment is defined by a symmetry action on topological excitations and defects, which can involve permutation and fractionalization.For subsystem symmetries, the same is true.However, as demonstrated in our examples, there is an additional non-trivial interplay between the spatial structure of the subsystem symmetries and the mobility and geometry of the topological excitations.The mobility restrictions on fractons are often formulated in terms of abelian conservation rules for the particles supported on various subsystems, which may be planar or fractal [69,70].Furthermore, such conservation rules on deformable co-dimension k subsystems encode that a local excitation must appear as part of an extended kdimensional excitation.The spatial intersection of these subsystem conservation rules and the subsystem symmetry generators lead to constraints on possible consistent symmetry actions that generalize those of a global symmetry in 2D.
To capture this we introduce the concept of fusion rules restricted to subsystems, assuming no fusion degeneracy for simplicity, we write N c ab | S = 0, 1, to indicate whether topological excitations a and b, supported on a subsystem S, may fuse to c, also supported on S, creating and annihilating no further excitations on S. We remark that this definition allows arbitrary excitations to be created and annihilated outside S, and applies to segments of extended topological excitations that are supported on S. For example, a looplike excitation restricts to pointlike excitations where it intersects a codimension-1 subsystem symmetry.This allows us to generalize symmetry enrichment with global symmetries to include subsystem symmetries, by replacing constraints on the action of a global symmetry from the quasiparticle fusion rules with constraints on the action of subsystem symmetries from the fusion rules restricted to the appropriate subsystems.
This paints a general picture for how subsystem symmetries can act on topological phases, via permutation and generalized projective representations that satisfy the consistency equations coming from all restricted fusion rules, or equivalently conservation laws.This generalizes the familiar classification of 2D SET phases.The projective representation of the full subsystem symmetry group should further respect locality in the sense that the action of subsystem generators that are far separated in space should commute, while those that act on a common excitation in the same spatial region need not commute.
As an example consider loop-like excitations in 3D, which obey the 1-form conservation law that every sphere is pierced by a loop an even number of times.Codimension-1 symmetries may act via nontrivial projective representations with Z 2 fusion (as in Section III) and locally permute the string excitations, which could involve attaching particles (as in Section IV C 1) or a segment of a general 1-dimensional defect [72,73] with Z 2 fusion rules.For a second example, fractons with planar charge conservation rules such as in the X-cube model [6,74] and other foliated fracton phases [9,75] can be acted upon nontrivially by planar subsystem symmetries that are aligned with the planar conservation rules [16].

No nontrivial SSET in 2D
Our discussion of SSETs above leads to the conclusion that there can be no nontrivial intrinsic SSETs in 2D, in other words all such SSETs are equivalent to a stack of a conventional SET and an SSPT.To see this we first point out that there are no fracton topological orders in 2D and hence all topological particles are fully mobile [76].Thus any topological charge can be moved to the complement region not acted upon by any subsystem symmetry generator.This implies the action of this subsystem symmetry on any topological sector must be trivial, involving no permutation nor projective representation.Although the string operators for the anyons may have unavoidable intersections with the subsystem symmetries, the defects created by the subsystem symmetries are in the trivial topological superselection sector, by definition of an SSET.Hence there can be no statistical processes between the anyonic string operators and subsystem symmetries either.
The argument above immediately generalizes to rule out nontrivial subsystem symmetry-enrichment on fully mobile point charges, and nontrivial actions by codimension-k, or above, subsystem symmetries on fully mobile (k − 1)-dimensional, or lower, extended excitations.

B. Outlook
Our SSET model can be obtained by decorating the 3D toric code with 2D cluster states.We can straightforwardly generalize our model by changing both the underlying 3D topological model, as well as the 2D SSPT model used to decorate.This raises the question of classification for SSET phases in 3D, and whether all phases can be captured by such decorated constructions.We remark that there should be some issues of compatibility, in that only certain combinations of topological model and SSPT are allowed.For example, in 2D SET models, the possible kinds of fractionalization that a point-like excitation can carry are restricted by its braiding statistics with other excitations [32].In analogy, we suspect that there is a connection between subsystem symmetry fractionalization in loop-like excitations and their braiding with point-like excitations.As a further generalization, it is also plausible that, by decorating with a 2D SSPT possessing 2D fractal subsystem symmetries [8,18,21,77], one could obtain a 3D SSET enriched by 3D fractal subsystem symmetries.These questions are closely related to the classification of SPTs with mixed global and subsystem symmetry, where the approach of Ref. [78] should be applicable.
The models analyzed in this work inherit their properties from the non-trivial interplay between global and subsystem symmetries.In particular, gauging all symmetries results in panoptic order, rather than "pure" fracton order in which all excitations have mobility constraints.It would be interesting to remove the global symmetries from the equation and consider systems with subsystem symmetries alone.Aside from providing novel mechanisms for subsystem symmetry fractionalization, this provides a route towards a systematic construction of "pure" fracton models with non-abelian fracton excitations [68,79,80], which would arise from gauging a symmetry-enriched fracton model in which fracton excitations carry a fractional subsystem symmetry charge.Additionally, by considering fractal subsystem symmetries, it may be possible to obtain type-II fracton models with non-abelian fractons which, as of yet, have proved elusive, previous attempts having resulted in panoptic models [54,55].
Finally, we address the possible applications of our models to the storage and processing of quantum information.In general, it would be interesting to investigate whether the introduction of subsystem symmetries and subsystem symmetry defects into topological models can augment their computational capabilities, as is the case in 2D topological systems [81][82][83].Regarding the SSET model from Section III, the fact that it combines the universal (measurement-based) quantum computational power of the SSPT order [14,21,26,27] together with the information storage capabilities of the topological order [84][85][86] suggests some potential applications.A relevant phenomenon that may be of use is the emergence of a symmetry protected degeneracy on the loop-like excitations.code Hamiltonian in the presence of this symmetry.We note that this same model appeared in Ref. [46], although it was derived from a different perspective.We recall that the ground states of the toric code can be described as loop condensates, which are equal weight superpositions of configurations in which the edge qubits in state |1 form closed loops.Furthermore, a configuration including an open string of 1's will create anyonic excitations, corresponding to violations of A v , at the endpoints.After placing CCZ e,∂(e) on every edge, these loops become decorated by 1D cluster states.Therefore, the anyons at the endpoints of open strings are accompanied by the edges of 1D SPTs.Since these edges transform non-trivially under the symmetry, the anyons fractionalize.
More precisely, consider the following string operator, where CZ ∂e acts on the two vertices touching e, and Γ is some open string of edges with terminal vertices v i and v f , which we assume to both be yellow without loss of generality.If we apply this operator to a ground state |ψ 2D SET , we get an excited state with excitations at v i and v f , corresponding to A v = −1.Because A v = −1 at these points, the projection 1+Av 2 annihilates the Hamiltonian terms involving C v i/f , so we can dress the endpoints of S Γ with Z's without changing the energy of the resulting excitations.This means we have four different string operators: Thus each anyon carries a two-fold degeneracy.Furthermore, Z 2 × Z 2 symmetry acts projectively on each anyon.
To see this, consider the subspace of degenerate states |a, b := S Γ (a, b)|ψ 2D SET .We compute, Therefore, X B ∼ X ⊗ X and X Y ∼ Z ⊗ Z in this subspace.Since X and Z anticommute, each anyon carries a projective representation of Z 2 × Z 2 , which demonstrates the fractionalization.

Gauging vertex symmetries
Now, we instead gauge the Z 2 ×Z 2 symmetry subgroup generated by X B and X Y .For convenience, denote the subsets of all blue/yellow vertices by V B/Y ⊂ V .Notice that V B/Y each form a rotated square lattice.We can therefore gauge each Z 2 factor individually in the same way that we gauged X F in the previous subsection, such that the vertex qubits are mapped onto a pair of new qubits on each face.This results in three qubits per face which we label f r,b,y with f r labelling the ungauged face qubit.
The gauged Hamiltonian reads, where, if we let F v ⊂ F denote the four faces surrounding vertex v, We have again projected certain terms onto the zeroflux subspace to ensure the Hamiltonian is symmetric.If we apply CCZ fr,f b ,fy to every face, this Hamiltonian is disentangled to two copies of the toric code on the blue and yellow face qubits, and a trivial Hamiltonian on the red face qubits.Once again, this unitary does not respect the residual Z 2 symmetry.This time, the residual symmetry acts in a way that permutes the anyons of the two toric codes, rather than fractionalizing.To see this, consider the following string operator, where Λ is a connected path of faces.If we apply this operator to a ground state |ψ g of H SET , we get an excitated state with a pair of excitations of A b v at the two yellow vertices at the endpoints of the path.Now, if we apply the symmetry to this state, we find, The string operator f ∈Λ Z fy acting on |ψ g create excitations of B y v at the same yellow vertices.Therefore, acting with X F permutes the anyons of the gauge theory.We can repeat the same procedure starting with different string operators.

Gauging all symmetries
We can also consider gauging the whole Z 2 × Z 2 × Z 2 symmetry group.While we will not perform this calculation explicitly, we can use some general results to determine the resulting gauge theory.First, suppose we start from H SET and gauge the fractionalized Z 2 × Z 2 symmetry.In general, when a symmetry acts on the anyons as a higher-dimensional representation, gauging it results in non-abelian anyons [32].The fractionalization pattern observed for H SET corresponds precisely to that shown in Fig. 2(c) of Ref. [90], and it is shown therein that gauging the symmetry results in a model with D 8 topological order, where D 8 is the symmetry group of a square.Now, let us rather start from H SET and gauge the anyon-permuting symmetry.In Ref. [54], it was shown in general that gauging an anyon-permuting symmetry leads to non-abelian anyons.The symmetry considered in Ref. [54] was a layer swap of a bilayer toric code, which appears superficially different from the permutation observed for H SET .However, upon a relabelling of anyons that preserves the braiding and fusion rules, the two symmetries turn out to be the same.Therefore, the resulting gauge theory should be the same as the one found in Ref. [54], namely D 8 topological order.
As the two above cases show, we expect a non-abelian D 8 topological order after gauging all symmetries of |SP T , regardless of the order in which the symmetries are gauged.This reinforces the intuition that the gauging operations for a pair of commuting symmetry subgroups should commute.

Symmetry defect analysis
Now, we show that examining symmetry defects in the ungauged model allows us to draw the same conclusions without explicitly gauging the symmetries.Using the definition of symmetry defects from Section IV A, we find that the symmetry defects associated to X F appear at the ends of 1D lines consisting of CZ's on a path along the edges of the lattice.Then, we find that the endpoint of the X F defect line transforms projectively under the X B/Y symmetries, in the exact same way as the string operator S Γ in Eq. (A10).The defects associated to X B (X Y ) appear at the ends of strings of CZs on a path alternating between face qubits and yellow (blue) vertex qubits.Consider the X B defects, with endpoints on yellow vertex qubits.Then, the X F symmetry transforms the defects by dressing these endpoints with Z operators.These correspond to X Y symmetry charges.If we were to then gauge the Z 2 × Z 2 symmetry generated by X B/Y , these two charges would become a pair of anyons.Therefore, the permutation action of the symmetry in the gauged theory can be seen by how the symmetry decorates the symmetry defects with symmetry charges.

B. Boundary Hamiltonians for |SSP T
In this Appendix, we consider some possible Hamiltonians which respect the boundary symmetries of |SSP T .Recall that, on the boundary, the planar subsystem symmetries act like lines of X operators, while the global symmetry acts like CZ's between neighbouring edges.The two simplest Hamiltonians that respect these symmetries are, X e e ∈n(e) Z e (B1) where ∂ E denotes the set of edges on the boundary, and n(e) contains the set of four edges that are nearestneighbouring edges to e, as measured by distance to the center-points of each edge.In H cSSP T , the global symmetry interchanges the two sums.Note that the edge qubits lie on the vertices of the medial square lattice.H cSSP T is exactly a 2D cluster Hamiltonian in an external field, tuned to its critical point [23,91].This critical point corresponds to a first order phase transition between the 2D SSPT phase of the cluster state and the trivial phase [22].Being first order, the phase transition is caused by a level crossing, such that H cSSP T has two degenerate ground states with a gap above them.In fact, these two grounds states are related by the boundary action of the global symmetry.H SSB , on the other hand, corresponds to two decoupled plaquette Ising models [15], one on the vertical edges, one on the horizontal edges.The ground states of this model spontaneously break the subsystem symmetries, such that there is an extensive number of degenerate ground states, with a gap above them.
We see that both of the above Hamiltonians have a degenerate ground space with a finite gap above.It is interesting to compare this to the analogous situation which arises on the boundary of a 2D SPT order [89].In that case, the boundary system is a 1D chain.The minimal Hamiltonian terms that commute with the boundary symmetry correspond to the 1D cluster Hamiltonian in a magnetic field tuned to its critical point, and two decoupled 1D Ising models.Thus the situation is similar to the current one.A crucial difference, however, lies in the fact that the ground state of the critical 1D cluster Hamiltonian respects the anomalous symmetry, and is therefore gapless [89].Furthermore, adding the Ising interaction on top preserves the criticality in a finite region [92].Therefore, the 2D SPT supports a boundary with symmetry-protected gaplessness.Conversely, it is not immediately clear how to engineer a gapless boundary for |SSP T .This is similar to the situation for 2D SSPT phases, which also only support degenerate boundaries [15].We leave a more detailed analysis of the boundary of 3D SSPT phases to future work.

C. Calculation of topological entanglement entropy
In this Appendix, we compute S (2) A for the bipartition depicted in Fig. 6(a).The subsystem A is defined by two planar cuts that run parallel to subsystem symmetry planes and intersect edges/faces of the lattice, such that the qubits on intersected edges and faces lie within subsystem A (ie.A has rough boundaries rather than smooth ones).The boundary between the two subsystems is therefore two disconnected tori.In a slight misuse of our notation, let ∂A E ⊂ E denote the edges intersected by the bipartitioning planes on each end of A. Likewise, let ∂A F ⊂ F denote the intersected faces, and let ∂A = ∂A E ∪∂A F .Since the boundary is made of two disconnected pieces, we can break each set into left and right parts, i.e.
Let G be the abelian stabilizer group defined in the main text.Then, since |SSET is the unique state satisfying g|ψ = |ψ ∀g ∈ G, we can write [50,93], In this case, the trace over B traces over one of the qubits involved in half of the CZ operators in C e .This does not give 0, rather we have Tr a CZ ab = 2P b where P = 1+Z 2 .Using this, we find, where we have defined the projector P E as, with, where C e is a unitary operator obtained from C e by removing half of the CZ's, see Fig. 9.The subgroup G A is defined to be generated by all elements of G which act non-trivially only on A, except for the two operators e∈∂A L/R E C e , which we exclude for notational convenience.
To compute ρ 2 A , we observe that [C E , g] = [P E , g] = 0 for all g ∈ G A and E ⊂ ∂A E .Then we get, where we have used the facts ( g∈G A g) 2 = |G A |( g∈G A g) and C E C E = C E⊕E .We observe that ρ 2 A is not proportional to ρ A because of the presence of the projectors P E .Therefore, ρ A is not a projector, which shows that the entanglement spectrum of our model is not flat, as it would be for the 3D toric code.Now we take the trace of ρ (C8) where we used that P E P E = P E .At this point, it is useful to consider the two boundaries of A separately.We can write, Observing that the final trace over ∂A F can further be factorized as Tr ∂A F = Tr ∂A L F Tr ∂A R F , and the left and right boundaries are disjoint and equivalent so we can simply square the result for the left boundary, we get, For each E this remaining trace is 0 unless g = e, or g = e∈E A e .In the latter case, we have P E e∈E A e = P E .Therefore, every term in the sum over E gets doubled, except when E = ∅ or E = ∂A L E .This gives, (C12) We now proceed by expressing the remaining sum in terms of the partition function of a 2D square lattice Ising model.To this end, we define an auxillary square lattice system with degrees of freedom σ i = ±1 on each vertex near the boundary.Each vertex i corresponds to an edge e i ∈ ∂A L E .Given E, we define a vector σ such that σ i = −1 if e i ∈ E, and σ i = 1 otherwise.With this notation, we can rewrite, P E = i,j P l(ei,ej ) with the trace given by, Tr P E = i,j This Ising model has a phase transition at inverse temperature β c = ln(1+ √ 2) 2 [67].Since ln √ 2 < β c , our partition function lies in the disordered phase.In this phase, the free energy is extensive, in the sense that F(β) := ln Z(β) = αN for large N .Importantly, there is no constant term in ln Z(β), as there would be in the ordered phase.The 2D Ising model has been solved exactly in the large-N limit by Onsager [67], and is given in Ref. [94] as, where k = 2 sinh(2β)/ cosh 2 (2β).Evaluating this integral numerically for β = ln √ 2, we find that α is equal to ln 2(ln 2 + Note that e ln 2(ln 2− 1 2 ) > 1, such that we can approximate ln e ln 2(ln 2− The final step is to evaluate |G A |. To do this, we need to count the number of indepedent generators of G A .This is done using the usual counting arguments for the 3D toric code [42].Suppose for simplicity that the region A contains L × L × L vertices, such that N = L 2 .We then have 3L 3 + L 2 edges contained in A, and 3L 3 + 2L 2 faces, giving |A| = 6L 3 + 3L 2 .Starting with the body center terms B c , we have L 3 − L 2 terms contained in A, all of which are independent.We have 3L 3 edge terms A e , but they are not all independent.Namely, the product of all edge terms around a given vertex is the identity, and so is the product of all edge terms on a non-contractible plane (such as the plane Σ z in Fig. 6 plus corrections which go to zero as L goes to infinity.

Figure 1 .
Figure 1.(a) Schematic diagram of a 2D SSPT order with line-like symmetries (dashed line) in the presence of a boundary (thick red line).Near a boundary, neighbouring line symmetries locally anti-commute (small circles), leading to an extensive degeneracy on the edge.(b) The ground state of the 3D SSET order has planar subsystem symmetries, and is a condensate of closed membranes of 2D SSPT orders.Left: a closed membrane of 2D SSPT order.Planar operators intersect the membrane along a loop, indicated by the dashed line, and reduce to the line symmetries of the membrane.Right: an open membrane carries a line-like excitation on its boundary, indicated by the thick red line, which transforms non-trivially under subsystem symmetries in the same way as the boundary of the 2D SSPT order.

Figure 3 .
Figure 3. (a) Top: a unit cell of the 3D SSPT state.The four triangles show a subset of the triangles appearing in Eq. (2).There are four more triangles like these for each face within this cell.The red lines represent the CZ gates contained in Uc (Eq.(4)).Bottom: the Hamiltonian term Ce, where e is the central vertical edge.(b) Decoration by 2D cluster states.The colored faces depict either domain walls, in the case of |SSP T , or membranes, in the case of |SSET .

Figure 4 .
Figure 4.The 3D cylinder with periodic boundary conditions in the x and y directions, and open boundaries in the z direction.The effective action of the bulk symmetries on the boundary is shown.The red lines denote CZ gates coming from the global symmetry action, while blue lines indicate the action of xz planar symmetries, with blue dots denoting X operators.The plane in the yz direction represents a domain wall that is inserted to create a global symmetry flux through the cylinder.

Figure 5 .
Figure 5. (a) Top: a unit cell of the 3D SSET state.The four triangles on the front face show a subset of the triangles appearing in UCCZ .There are four more triangles like these for each face within this cell.Bottom: The Hamiltonian term Ce, where e is the central vertical edge.(b) A loop-like excitation appearing on the boundary of the surface Σ, indicated by shaded faces.Excitations appear on edges marked by points.(c) A 2D cut of our 3D lattice with a loop excitation corresponding to the half-infinite membrane Σ ∞ , as viewed from above.There is an effective qubit degree of freedom for each edge on the boundary of Σ ∞ , indicated by dots.The horizontal dashed lines represent the intersection of two symmetry planes P with this 2D slice, and the action of XP on the edge degrees of freedom, according to Eq. 27, is shown beside each line.

Figure 6 .
Figure 6.(a) The geometry considered throughout this section.Opposite sides are identified, resulting in a 3D torus.The A and B subsystems are shown, as well as the noncontractible loops Λx,y and membrane Σz.(b) The dumbbell configuration for detecting the SPEE in 2D systems fromRef.[52].(c) The picture frame configuration proposed to detect the SPEE that arises due to planar subsystem symmetries in 3D systems.Here we see a view from the top.A cross section along the dotted line reveals a shape identical to that of (b).(d) The alternate picture frame configuration that does not detect SPEE due to line-like symmetries.

2 , Z sub1 2 , and Z sub2 2
Symmetry defects are defined to appear at the boundaries of open domain walls.The Z glob 2 symmetry defects are closed 1D loops, while the Z sub 1/2 2 symmetry defects are point-like.The symmetry defects carry important information about what happens after gauging the symmetry.After gauging, domain wall operators are proliferated, and the symmetry defects become deconfined topological excitations that we call gauge fluxes.As we saw in Section III, Z glob 2 gauge fluxes are mobile loop-like excitations.Gauge fluxes corresponding to generators of Z sub 1/2 2

Figure 7 .
Figure 7. Hamiltonian terms of the fracton model obtained from gauging subsystem symmetries of |SSP T .Top left: the minimal interaction term symmetric under all subsystem symmetries.The gauge qubit (blue dot) is placed at the nexus of the interaction.Bottom left: labelling of the 12 gauge qubits within a cube.Top middle: the flux term Ae, where e is the central vertical edge.Top right: the Hamiltonian term Ce, where the red line denotes CZ and connects a gauge qubit to a body center qubit (red).Dashed lines are used as visual cues that indicate which edge each gauge qubit is associated to.The interaction is only shown in one quadrant for clarity; it acts in an analogous (rotated) manner in all 4 cubes surrounding the edge e as indicated by the dashes.Bottom right: the flux terms within the cube.Z• is shorthand for Z (c,•) .The flux terms Ae and A (i) c satisfy a number of relations which determine the mobility of gauge fluxes, such as the product of all Ae and A (7) c in an xz dual-plane; Ae, A (4) c , A (5) c in a yx lattice-plane; and Ae, A

2 and Z glob 2 , 2 generator creates Z sub1 2 gauge
we end up with a model that is, in the absence of symmetry, equivalent to a stack of the 3D toric code and the fracton model from Section IV B 1. The ungauged Z sub2 2 symmetry couple the two models by permuting excitations between them.Specifically, a Z sub2 charges on a pair of planes wherever it intersects the loop-like Z glob 2 gauge fluxes, and it attaches a Z glob 2 gauge charge to Z sub1 2

2 and Z glob 2 2 gauge fluxes and by attaching Z glob 2 gauge charges to Z sub1 2 gauge
gives a model in the same topological phase as a stack of the 3D toric code and layers of 2D toric codes in all three directions.Similar to the previous model, the remaining Z sub1 2 symmetry generators permute excitations between the 2D and 3D toric codes by attaching Z sub2 2 gauge charges to the loop-like Z glob fluxes in adjacent planes .

Figure 8 .
Figure 8.(a) The Union-Jack lattice.Filled triangles represent the four CCZ unitaries associated to one face.The shaded circle contains the Hamiltonian term Cv, where X acts on v and red lines represent CZ unitaries.(b) The decorated domain wall structure of |SP T .Vertex qubits are initialized in the |+ state, then CZ unitaries are applied along thick red lines.
If we let m b/y and e b/y denote the anyons associated to violations of A b/y v and B b/y v , respectively, then the residual Z 2 symmetry permutes anyons by attaching to every m particle an e particle of opposite color, i.e. m b/y ↔ m b/y e y/b .

Figure 9 .
Figure 9.After tracing out qubits in the B subsystem, as, indicated by the slashes, the operator Ce (left) is transformed to C e Pe up to a constant factor of 2 4 = 16.

Where
|A| and |B| are the number of qubits in subsystems A and B. When we trace out the B subsystem, all elements in G which have non-trivial support in B have zero partial trace, since the Pauli operators are traceless.The only exception are the operators C e where e ∈ ∂A E .
P ∅ = Tr P ∂A L E = 2 |∂A L F | , we write,

2 .
the sum in such a way that it corresponds to the partition function Z(ln √ 2) of the 2D Ising model, where, Z(β) = σ e β i,j σiσj .(C17) For conceptual clarity, let us define the group G A which is obtained by adding the two generators e∈∂A L/R E C e , which we omitted earlier, into G A .Then we have |G A | = 4|G A |.By letting N = |∂A L E | be the number of spins in the Ising model, we get, S (2) A = − ln Tr A (ρ 2 A ) = |A| ln 2 − ln |G A | − 2 ln(2 −N Z(ln √ 2) − 1).(C18)

1 2 )N − 1 ≈ ln e ln 2 A
= |A| ln 2 − ln |G A | − 2 ln 2(ln 2 − 1 2 )N .(C21) (a), giving L 3 + 1 constraints, and hence 2L 3 − 1 independent terms.Finally, we have one C e term for each edge away from the boundary, all of which are independent, and the two non-local terms e∈∂A L/R E C e , giving 3L 3 − L 2 + 2 terms.Finally, we have the three non-local operators S m Σz and S e Λx,y .Altogether, this gives 6L 3 − 2L 2 + 4 terms, such that |G A | = 2 6L 3 −2L 2 +4 .The final result for the entropy is therefore, S (2) A (L) = ln 2(6 − 2 ln 2)L 2 − 4 ln 2 , (C22) (b).Let C ⊂ C be a set of body center qubits, and define the state |C such that each qubit in C is in the |1 state, while the rest are in |0 .Noting that CCZ|i |j |k = |i ⊗ CZ i |j |k , we can write the ground state in the following way, 2A in three steps, Tr A = Tr ∂A F Tr ∂A E Tr A−∂A .When taking the first trace, all non-trivial elements of G A with support outside of ∂A are traceless, except again for C e for e ∈ ∂A E .This gives, E after removing all CZ's.G ∂A contains all elements of G A which have act non-trivially only on ∂A, and is generated by the operators A e for e ∈ ∂A E .Now, the trace over ∂A E is 0 unless E = E , giving,Tr ∂A E Tr A−∂A (ρ 2 A ) = 2 −|∂A F | 2 |A| |G A |