Strong planar subsystem symmetry-protected topological phases and their dual fracton orders

We classify subsystem symmetry-protected topological (SSPT) phases in $3+1$D protected by planar subsystem symmetries, which are dual to abelian fracton topological orders. We distinguish between weak SSPTs, which can be constructed by stacking $2+1$D SPTs, and strong SSPTs, which cannot. We identify signatures of strong phases, and show by explicit construction that such phases exist. A classification of strong phases is presented for an arbitrary finite abelian group. Finally, we show that fracton orders realizable via $p$-string condensation are dual to weak SSPTs, while strong SSPTs do not admit such a realization.

Introduction-Global symmetries, such as the Z 2 spinflip symmetry of the Ising model, act throughout the bulk of a system. Recently, there has been an emerging interest in symmetries that act on only part of a system. These include higher-form symmetries which act on deformable lower-dimensional manifolds of a system, 47 as well as subsystem symmetries, [5][6][7] which act on rigid lower-dimensional subsystems. It has also been realized that such subsystem symmetries may protect nontrivial symmetry-protected topological (SPT) phases 1-4 : gapped, disordered, short-range entangled phases which cannot be adiabatically connected to the trivial disordered phase in the presence of symmetry, but can be if the symmetry is not enforced. Examples of subsystem symmetries include those which act along straight lines 8,9 or planes, 10,16 or even fractal subsystems. [11][12][13][14][15]17 Such phases have been aptly named subsystem SPT (SSPT) phases, and the present work concerns their classification.
In 2 + 1D, such systems have gained interest due to the discovery that non-trivial SSPT phases may serve as a resource for universal measurement-based quantum computation (MBQC) [48][49][50][51]53 and also due to their unusual patterns of quantum entanglement. 9,[64][65][66][67] In attempting to classify such SSPT phases, one is faced with the issue that the total symmetry group is infinitely large in the thermodynamic limit, and consequently there are infinitely many distinct phases. However, many of these phases can be constructed by stacking (a process which we will define) 1+1D SPT phases along the subsystems. We call such phases weak SSPTs, whose nontriviality are a manifestation of lower dimensional physics. Ref. 9 defined an equivalence relation between phases wherein two phases that differ by stacking 1+1D SPTs belong to the same equivalence class. Phases not in the trivial equivalence class are, by definition, strong SSPTs. Ref. 9 found that there are a finite number of equivalence classes of phases, thus providing a sensible classification for the infinitely many SSPT phases. The present work is the natural extension of Ref. 9 to planar symmetries in 3+1D (henceforth, simply 3D).
The brief history of 3D planar SSPT phases begins with Ref. 10, which constructed a non-trivial 3D planar SSPT model. However, it was later discovered that its fracton dual belonged to the same foliated fracton phase as the X-cube model, 63 implying that it is weak. More recently, fracton phases were constructed in Ref. 62 which possess 'twisted' foliated fracton orders, raising the question as to the nature of their SSPT duals. We find that these phases, too, are weak. This prompts the question: do any strong planar SSPTs exist? We answer this in the affirmative.
We will first show how to construct weak 3D planar SSPT phases via a stacking process of 2D SPTs. We then ask whether there are SSPT phases which cannot be realized by this process. We identify mechanisms by which an SSPT may be strong, leading to a classification of such phases, and construct exactly solvable, zero-correlation length models realizing these phases. In the fracton dual picture, this construction corresponds to one in which 2D topological orders are stacked on to and strongly coupled to an existing fracton model in a certain way. 62 The duals of our strong SSPTs are novel fracton phases which cannot be attained via such a procedure, which also implies that they cannot be realized by a p-string condensation transition, 45,46 as we will show.
Planar subsystem symmetries-Throughout we will consider a system with degrees of freedom on each site of a cubic lattice. Each site r transforms under the finite abelian on-site symmetry group G under a unitary representation u r (g). An xy planar symmetry acting on plane z acts as S xy (z; g) = x,y u r=(x,y,z) (g) for g ∈ G. Similarly, we may define S yz (x; g) and S zx (y; g), which act on yz and zx planes respectively. Importantly, individual sites transform under the same on-site representation regardless of the orientation of the planar symmetrythere is therefore a redundancy: the product of all xy symmetries is identical to the product of all yz or all zx symmetries. We will refer to models which respect only one type of planar symmetry as 1-foliated, those with two as 2-foliated, and those with all three as 3-foliated. To construct explicit models, we choose the on-site degrees of freedom to be G-valued, |g r , which transform under the on-site symmetry as u r (g) |g r = |gg r .
Construction of weak SSPT phases-It is possible to construct non-trivial SSPT phases from known 2D global SPTs, as we will show in this section. Phases obtained in this way are 'weak', by definition, whose nontrivial properties are in some sense a manifestation of lowerdimensional physics. We emphasize here that we do not assume any translation invariance in our system. First, we briefly review the group cohomological classification of 2D bosonic SPTs with global symmetry group G. 3,20 For the purpose of being self-contained, we also include a more detailed review in Appendix A. The classification of such phases 3 is given by the third cohomology group H 3 [G, U (1)]. For simplicity, we may consider G = (Z N ) M , in which case an element of H 3 [G, U (1)] is specified by integers, p i I , p ij II (i < j), and p ijk III (i < j < k), all modulo N , called type-I, II, and III cocycles respectively. We will specify p i I and p ij II compactly in a single symmetric integer matrix M with M ii = 2p i I and M ij = M ji = p ij II . Upon gauging the global symmetries of a 2D SPT, one obtains a topologically ordered system with fractional quasiparticles carrying gauge charge or flux (or both). Nontrivial type-III cocycles give rise to non-abelian topological order -as we are only interested in SSPTs with abelian fracton duals, we will not consider them here. The elements of M characterize the self and mutual statistics of gauge flux excitations. 20 In particular, the type-I cocycles give rise to a self exchange statistic e πiMii/N 2 of the gauge flux m i , and type-II cocycles lead to a mutual braiding statistic of e 2πiMij /N 2 FIG. 1. (Left) Examples of our construction of 1-foliated or weak 2 or 3-foliated models, for G = ZN × ZN , in the graphical notation. 2D SPTs to be stacked, are shown in the blue boxes, and the arrow points to the resulting SSPT after stacking. The color of the edges connecting two vertices indicate its weight modulo N . (Right) Examples of M matrices that cannot be obtained by stacking 2D phases onto 2 or 3foliated models. The Type 1 phase is only strong for even N , and Type 2 strong phases can only be realized for 2-foliated symmetries.
between m i and m j .
It is always possible to view a 3D planar SSPT as a quasi-2D system in the xy plane with a subextensively large symmetry group G L by compactifying the z direction. We may then proceed to compute its classification in terms of H 3 [G L , U (1)], which is characterized by a subextensively large M matrix. We note that it is possible to define M matrices corresponding to yz or zx as well, but for reasons that will become clear we will always consider the xy symmetries only. It is useful to introduce a graphical notation for M, which is used in Fig. 1. The αth generator of G in a plane z is denoted by a vertex a i=(α,z) . Two vertices i and j are connected by an undirected edge with weight M ij , and a vertex i is connected to itself via a self-loop with weight M ii /2, where weights are defined modulo N .
Consider the 2D global symmetry group G 2D = G K for an integer K. For appropriate choice of the pure phase function f 2D , the wavefunction |ψ 2D = {gr} f 2D ({g r }) |{g r } 2D is a zero-correlation length ground state of a commuting-projector Hamiltonian with SPT order. All phases in the group cohomology classification can be realized in this way [24][25][26] (see Appendix A).
Suppose we start with the trivial disordered wavefunction |ψ 0 = {gr} |{g r } 2D . We can construct a nontrivial 1-foliated SSPT by identifying each factor of G in G 2D in the function f 2D ({g r }) with a planar G symmetry in the arbitrary collection of planes z 1 , . . . , z K (where z k are all within some finite range to ensure locality). The wavefunction |ψ 1-fol = U |ψ 0 with U = {gr} f 2D ({g r } rz∈{z k } ) |{g r } {g r }| is the ground state of a 1-foliated 3D SSPT, which is nontrivial only near the planes z k . We may then repeat this procedure arbitrarily many times, each time acting on the previous state with U for different choices of f 2D and {z k }. We will call this procedure "stacking" the 2D SPT |ψ 2D onto the planes {z k }.
More generally, we may define a stacking operation between two SSPTs in which the two systems, with on-site symmetry representations u (1) r (g) and u (2) r (g), are placed on top of each other to create a new SSPT with on-site representation u r (g) = u (1) r (g)⊗u (2) r (g). The group structure of the standard SPT classification is realized under such a stacking operation. Stacking a 2D SPT onto a 3D SSPT can be viewed as stacking two 3D SSPTs, in which the first is only nontrivial in the vicinity of a number of planes {z k }. We define any phase realizable by stacking 2D SPTs in this way to be weak. In the case of our 1foliated SSPT construction, each additional stacked 2D SPT simply adds to the corresponding elements of M, shown graphically in Fig 1. For 1-foliated symmetries, it is thus possible to realize any M by stacking 2D SPTs; hence all phases are weak.
On the other hand, for 2-or 3-foliated models, this procedure may not work because |ψ 1-fol is not guaranteed to be symmetric under the orthogonal planar symmetries (if it is, we can simply follow the same procedure). Instead, let us define variables d r = g r+z g −1 r , which transform under the xy planar symmetries but are invariant under all orthogonal symmetries. We may then define non-trivial SSPT wavefunctions as before, but in terms of d r instead using the unitary which is explicitly invariant under the orthogonal symmetries. However, in this case the M matrix of the 2D SPT does not map directly onto that of the SSPTinstead one should view the 2D SPT as living "in between" the planes of the SSPT, at {z k + 1/2}. To obtain the M matrix of the SSPT, one can compute the appropriate type-I and II cocycles of the 2D SPT in the basis of the xy planar symmetries. For details of this basis change, see Appendix F. This process is shown in Fig 1. The consequence is that, as opposed to 1-foliated symmetries, there are allowable M matrices that cannot be realized by stacking any number of 2D SPTs, due to certain constraints on M which we discuss in the next section.
Note that in this discussion we have implicitly ignored nontrivial SSPTs that have trivial M matrices. Such phases do in fact exist as discussed in Appendix C. However, we conjecture that all such phases are weak (they can be realized by stacking 2D linear SSPTs 9 ) and therefore irrelevant in the classification of strong phases.
General constraints and invariants-In the presence of orthogonal symmetries, there are general constraints that must be satisfied by M. Conceptually, these arise due to the aforementioned redundancy: the global symmetry S glob (g) = z S (xy) (z; g) = x S (yz) (x; g). Since yz symmetries do not contribute to M, the generator S glob (g) must therefore manifest trivially in M. This leads to two types of constraints on the elements of M: the global symmetry must have trivial type-I cocycle with itself and trivial type-II cocycle with any other symmetry.
In Appendix E, we prove that these constraints must hold generally by analyzing the possible action of the symmetry on the edges. 19 Let us label the αth generator of G on the zth plane by i = (α, z). Then, the two constraints are expressed as z M (α,z),(β,z ) = 0 mod N, ∀α, z, β and 1 2 z,z M (α,z),(α,z ) = 0 mod N, ∀α. These constraints have a natural interpretation in terms of our graphical representation. They define a restricted subgroup of H 3 [G L , U (1)] in which 2-or 3foliated SSPTs must reside. As we will show, there are now allowed phases which cannot be realized by stacking any number of 2D SPTs -these are precisely the strong phases we are searching for. This motivates us to define two types of strong invariants, F 1 and F 2 , which cannot be changed by stacking with 2D SPTs.
Strong SSPTs: Type 1-Consider G = Z 2N . Then M zz is an L × L matrix. We pick an arbitrary cut that divides the system into two halves z < z 0 and z ≥ z 0 . Then, is a Z 2 -valued global invariant. To see why, let us view M zz mod 2 as a Z 2 "flux" flowing from vertex z to z in the graphical representation. Then, Eq. 2 is a divergencefree constraint at each vertex. The invariant F 1 is simply the total Z 2 flux flowing through a cut at z 0 . It is therefore clear that F 1 does not depend on the choice of cut z 0 , nor can it be modified by stacking a 2D SPT which amounts to adding closed flux loops locally. Type 2-Consider G = Z N × Z N , so that M (α,z),(β,z ) is a 2L × 2L matrix. Again we pick a cut z 0 . Then, (1,z ) mod N (4) is a Z N -valued global invariant. To see how this arises, let us interpret M (1,z),(2,z ) as a Z N "flux" flowing from vertex (1, z) to (2, z ). Like before, Eq. 2 is a divergencefree constraint on this flux and F 2 measures the total flux flowing across a cut, which therefore does not depend on z 0 nor can it be modified by stacking with 2D SPTs.
In Appendix F, we prove three important statements. First, that the invariant F 1 or F 2 is the same regardless of whether we consider the M matrix obtained from xy symmetries or that obtained from yz (or zx) symmetries. Secondly, 3-foliated systems must have trivial F 2 = 0. Thirdly, F 1 and F 2 (which we also define for general G) completely characterize M up to stacking with 2D SPTs. Finally, in Appendix G, we provide an explicit construction of a 3-foliated model which realizes a nontrivial type 1 strong phase F 1 = 1, and a 2-foliated model which realizes arbitrary F 1 and F 2 , thereby demonstrating the existence of such strong phases. Examples of M matrices with non-trivial F 1 and F 2 are shown in Fig. 1 (right).
Let us define a 'strong' equivalence relation between SSPTs, under which two phases belong to the same equivalence class if they can be connected with one another by stacking of 2D SPTs 70 (along with, of course, symmetric local unitary transformations and addition/removal of disentangled degrees of freedom transforming as an on-site linear representation of G 23 ). For an arbitrary finite abelian group G, the set of equivalence classes is given by for 3-foliated and 2-foliated models respectively. The group structure is realized via the stacking operation between two SSPTs. We note that this equivalence relation can be naturally formulated in terms of planar-symmetric local unitary circuits, generalizing the definition of Ref. 9. Indeed the unitaries U used to construct weak SSPTs are examples of such circuits. Fracton duals-It is well known that, under a generalized gauge duality, [16][17][18] SSPT phases map onto models of fracton topological order. 10,62 The simplest and most well-studied fracton model is the X-cube model, 16 which is obtained by gauging the planar symmetries of the plaquette Ising paramagnet, and hosts fractional quasiparticle excitations with limited mobility including immobile fractons, lineons mobile along lines, and planons mobile within planes (which are either fracton dipoles or lineon dipoles). For our discussion, we will assume that the reader has a rudimentary understanding of the Xcube fracton model and its quasiparticle excitations (see Ref. 31 for a review).
Let us begin with 3-foliated SSPTs, which are dual to 'twisted' X-cube fracton topological orders with fractonic charge. 62 The gauge flux m (g,z) of an element g on the plane z is a planon: a composite excitation composed of a lineon anti-lineon pair on the planes z +1/2 and z −1/2, i.e. a lineon dipole. A single lineon can be regarded as a semi-infinite stack of lineon dipoles mobile in the x and y directions. For a more nuanced discussion of the mobility of such excitations, see Appendix D.
The constraints on the matrix M have a simple interpretation in this language: the infinite stack of lineon dipoles, which belongs to the vacuum superselection sector, 69 must have trivial braiding statistics with all other lineon dipoles, and a trivial exchange statistic with itself. The invariant F 1 also has a simple interpretation in this picture: the quantity e 2πiF1/N 2 corresponds to the braiding statistic 71 of a lineon and its anti-lineon on the same plane, modulo e 4πi/N 2 .
It is possible to construct fracton topological orders by strongly coupling intersecting stacks of topologically ordered 2D discrete gauge theories oriented along the xy, yz, and zx planes, inducing a type of transition called pstring condensation. 45,46 More generally, these stacks of 2D gauge theories can be replaced by arbitrary 1-foliated gauge theories. 62 The twisted X-cube models that emerge from this construction are dual to weak 3-foliated SSPTs constructed via the planar-symmetric local unitaries U in Eq. 1. We walk through this correspondence in more detail in Appendix H.
Conversely, strong 3-foliated SSPTs are dual to fracton models that cannot be realized through p-string condensation. This correspondence sheds light on the F 1 strong invariant -in p-string condensation, lineon crossing statistics are inherited from the self-braiding statistics of fluxes in the 1-foliated gauge theories, and are therefore the square of a flux exchange statistic, i.e. a multiple of e 4πi/N 2 for G = Z N with N even. In a strong phase, F 1 = 1 implies that this statistic is offset by e 2πi/N 2 . The fracton dual of the Type 1 strong G = Z 2 model, constructed in Appendix G, is an example of a novel such fracton order in which lineons satisfying a triple fusion rule have ±i mutual crossing statistic, and therefore cannot be realized via p-string condensation.
One can also consider the fracton duals of 2-foliated SSPTs, which are novel 'twisted' versions of the 2-foliated lineon-planon model introduced in Ref. 63. Furthermore, the X-cube model may be ungauged in two different ways, by regarding either the fracton sector or the lineon sector as gauge charge. The former procedure results in a paramagnet with G-valued degrees of freedom transforming under all 3 sets of planar symmetries as before, whereas the latter yields a model with two G-valued degrees of freedom per site, the first transforming under xy and yz planar symmetries, and the second under yz and zx planar symmetries. Both Type 1 and Type 2 strong SSPTs, as well as arbitrary weak SSPTs, may be constructed for the latter. Their fracton duals are novel variants of the X-cube model whose fracton dipoles exhibit non-trivial braiding and exchange statistics.
Conclusions-We have formulated a classification of strong 3D planar SSPTs. Each phase falls into one of a finite set of equivalence classes modulo stacking with 2D SPTs, which we have fully enumerated. For 1-foliated systems, all SSPT phases are weak. For 2-foliated systems, there are two mechanisms by which a phase may be strong, characterized by Type 1 and Type 2 invariants. For 3-foliated systems, only Type 1 strong phases exist. Under a generalized gauge duality, our classification has a natural interpretation in terms of p-string condensation, 45 and we have explicitly constructed strong SSPT models which are dual to fracton phases that cannot be realized via this mechanism.
There are various natural extensions of our work. A relevant and open question regards the structure of entanglement in strong SSPT phases. 56,[64][65][66][67] Another is the addition of non-trivial type-III cocycles, which are dual to non-abelian fracton topological orders. We leave the classification of such phases to future work. Finally, it would be interesting to study the foliation structure of the fracton duals.
Acknowledgments-T.D. thanks Fiona Burnell, Dominic Williamson, Abhinav Prem, and Shivaji Sondhi for many helpful discussions, especially in the early parts of this work. W.S. thanks Xie Chen and Sagar Vijay for helpful discussions. T.D. acknowledges support from the Charlotte Elizabeth Procter Fellowship at Princeton University. W.S. is supported by the National Science Foundation under award number DMR-1654340 and the Institute for Quantum Information and Matter at Caltech. J.W. was supported by NSF Grant PHY-1606531 and Institute for Advanced Study. This work is also supported by NSF Grant DMS-1607871 "Analysis, Geometry and Mathematical Physics" and Center for Mathematical Sciences and Applications at Harvard University.
In this Appendix, we review the group cohomological classification of SPTs in 2D, as well as some additional aspects which will prove useful for our arguments related to the SSPT. We will also review its interpretation as an anomalous action of the symmetries on the edges, as well as the connection to the braiding and exchange statistics of quasiparticle excitations in its gauge dual.

Group cohomological classification of 2D SPTs
In the presence of symmetry, the unique ground states of two gapped Hamiltonians belong to the the same phase if they can be transformed into each other via a symmetric local unitary (SLU) transformation. 3 That is, a finite depth local unitary circuit in which each gate commutes with the symmetry operation. A state describes a nontrivial 2D SPT phase if it cannot be connected to the trivial product state via an SLU, but can be trivialized if the symmetry restriction is removed. Two dimensional bosonic SPTs with on-site symmetry G, under this phase equivalence relation, are known 2 to be classified according to the third cohomology group H 3 [G, U (1)]. For the finite abelian group G = i Z Ni , this can be written out explicitly as where gcd denotes the greatest common denominator. The three factors are commonly referred to as type-I, type-II, and type-III cocycles. Type-III cocycles correspond to a gauge dual with non-abelian quasiparticle excitations; as our focus is on SSPTs with abelian fracton duals, we will be focusing only on type-I and II cocycles.
a. The Else-Nayak procedure Let us derive the group cohomological classification via a series of dimensional reduction procedures, introduced by Else and Nayak 12 , which will prove useful in our discussion of SSPTs. Although the original procedure observed a system with a physical edge, here we prefer to deal with a "virtual" edge, meaning: the full system has no edges, but we will consider applying the symmetry only to a finite region M of the system. At the edges of M , this symmetry will act non-trivially as if at a physical edge. The advantage of this approach is that it removes any ambiguity related to choice of how the model is defined at the physical edges (and will be useful in the case of SSPTs).
Let |ψ be the unique gapped ground state of our Hamiltonian H with on-site symmetry group G, and S(g) be the symmetry operation realizing the symmetry element g ∈ G. We have that [H, S(g)] = 0 and, without loss of generality, take the ground state to be uncharged under the symmetry S(g) |ψ = |ψ . Now, let S M (g) be the symmetry operation S(g), but restricted to a region M . S M (g) acting on the ground state will no longer leave it invariant, but will create some excitation along the boundary of this region, ∂M . Since |ψ is the unique ground state of a gapped Hamiltonian, this excitation may always be locally annihilated by some symmetric unitary transformation U ∂M (g) † , which only has support near ∂M . That is, It is straightforward to show that the matrices U ∂M (g) form a twisted representation of G, satisfying where B A ≡ BAB † denotes conjugation of A by B, and that they must commute with any global symmetry operation, [U ∂M (g), S(g )] = 0. We now perform a further restriction: from ∂M down to a segment C, U C (g). This is always possible. U C (g) need only satisfy Eq A3 up to some unitary operator V ∂C (g 1 , g 2 ) at the two endpoints of C, when acting on |ψ . The final restriction is from ∂C, which consists of two disjoint regions a and b, down to simply a: . V a (g) need only satisfy Eq A5 up to a U (1) phase factor, which can be cancelled out by the contribution from V b (g).
where ω : G 3 → U (1). This entire dimensional reduction process is shown in Figure 2.
One can further show that ω(g 1 , g 2 , g 3 ) satisfies the 3cocycle condition 12 and since V a (g 1 , g 2 ) is only defined up to a phase factor β(g 1 , g 2 ), we must identify is called a coboundary. The classification of functions satisfying Eq A7, modulo transformations Eq A8, is exactly the definition of the third cohomology group Suppose we have followed the Else-Nayak procedure on a system and obtained the cocycle function ω(g 1 , g 2 , g 3 ). How do we identify which class in Eq A1 it belongs to? One way to do so is to identify combinations of ω which are invariant under the transformation Eq A8, whose value can tell us about the class.
For simplicity, we focus first on G = (Z N ) M . Let us first write down an explicit form 9,10 for ω, where g i is an integer modulo N denoting the component of g in the ith Z N factor, g = (g 1 , g 2 , . . . , g M ), [·] means to take the modulo N , and p ij are integers mod N . It is straightforward to confirm that ω satisfies the 3-cocycle condition. As we will show, the different choices of p ij From Eq A1, p i I ≡ p ii specify the value of the type-I cocycles and p ij II ≡ p ij specify the type-II cocycles. Define and both of which one can readily verify are invariant under transformations of the type Eq A8. Given a choice of generators, G = a 1 , . . . , a M , an explicit calculation shows that and thus correctly identifying the value of the type-I and type-II cocycles. Thus, if we are given an unknown ω, we may simply compute Ω(a i ) and Ω II (a i , a j ) for all i and j to identify its class. We may define the symmetric matrix M ij = p ij II and M ii = 2p i I . Then, we have and for arbitrary elements g and h, where g = (g 1 , . . . , g M ).
c. Group cohomology models The group cohomology models are a powerful construction that allows us to explicitly write down models realizing SPT phases corresponding to an arbitrary cocycle 9,10 . Although these models have an elegant interpretation in terms of a path integral on arbitrary triangulations of space-time, we will simply be using them to define Hamiltonian models on regular lattices. We focus on the case of a square lattice.
Take G-valued degrees of freedom on each site r, |g r . The ground state of our model |ψ is an equal amplitude sum of all possible configurations where f ({g r }) is a U (1) phase for each configuration. The group cohomology model is defined by the choice The dimensional reduction procedure in the Else-Nayak procedure. We start with a truncated global symmetry operator, SM (g). This acts on the ground state as a unitary U ∂M (g) along the edge of M . We further restrict this unitary down to a line segment C, UC (g). Restricted to C, UC (g) behaves as a representation of G only up to unitaries V ∂C (g) at its endpoints. Finally, we restrict to a single endpoint Va(g), where associativity of the representation is only satisfied up to a phase ω(g1, g2, g3), defining our 3-cocycle.
where x, y are the two unit vectors, g * ∈ G is an arbitrary element which we can simply take to be the identity g * = 1, and we have defined a phase contribution f r for each plaquette. This arises from a triangulation of each square plaquette into two triangles, each of which contribute a phase; those interested in the details of the construction are directed to Ref 9. Performing the Else-Nayak procedure outlined in Appendix A 1 a on this ground state results in exactly the cocycle ω used to construct the state, up to a coboundary (Eq A8).
To obtain a gapped local Hamiltonian realizing this state as its ground state, we simply consider a set of local ergodic transitions {g r } → {g r } , multiplied by an appropriate phase factor, which by construction has |ψ as its unique ground state. We can simply choose {g r } to differ from {g r } by the action of a generator a i of G on a single site r. The Hamiltonian will then be a sum of mutually commuting terms consisting of a "flip" operator |a i g r g r | on each site, multiplied by an appropriate phase factor depending on the state {g r } near that site.

d. Gauge duality
The group cohomological classification of an SPT has an elegant interpretation in terms of braiding statistics of its gauge dual 11 . As our main interest is on the ungauged side of things, we will only very briefly outline the gauging process (as applied to the group cohomology models), and identify the relevant statistical processes.
Take as our ungauged SPT a group cohomology model on a square lattice. For each nearest neighbor pair (r, r ), we define gauge degrees of freedom g r,r = g r g −1 r . We then write the Hamiltonian (Eq A21) in terms of these degrees of freedom, which is always possible. In addition, we energetically enforce the constraint g r1r2 g r2r3 g r3r4 g r4r1 = 1 for the square plaquette with corners r 1...4 (labeled going clockwise or counterclockwise), by adding an appropriate projection term to the Hamiltonian.
The resulting model describes a topologically ordered phase, with characteristic properties such as a topological ground state degeneracy on a torus and quasiparticle excitations with anyonic braiding statistics. There are two main types of excitations: gauge charge, denoted by e g , and gauge flux, denoted by m g , for each g ∈ G. The former are created by gauged versions of operators of the form which creates a charge-anticharge pair, e g and e −1 g , at positions r 2 and r 1 . To create gauge flux excitations, instead consider the gauged version of the operator where S M (g) is a symmetry operator restricted to a region M and U ∂M (g) is the action on the boundary ∂M , as in the dimensional reduction procedure of Appendix A 1 a. The gauged version of S M (g) only flips g rr near at the boundary, and so the gauged L(g) operator has support only on ∂M . Now, if we further restrict L(g) → L C (g) to an open segment C, L C (g) creates two quasiparticle excitations at the two endpoints, which we identify as the gauge flux-antiflux pair m g and m −1 g . Note that there is an ambiguity in defining the gauged version of L(g), which may result in a different definition of the gauge flux excitation, m g ∼ m g e g . Thus, gauge fluxes are only well defined modulo attachment of charges.
The group cohomological classification of the ungauged SPT manifests in the self and mutual statistics of gauge fluxes in the gauged theory. Let a i be the generator of the ith factor of Z N in G, and e i and m i be its gauge charge and flux excitations. For two identical excitations, we can define an exchange phase via a process in which their two positions are exchanged. For two different excitations, we may instead define the full braiding phase, which is accumulated when one particle encircles another. with m j . Notice that the exchange and mutual braid of m i is only well defined modulo e 2πi N , since m i is only well defined modulo charge attachment. For a general gauge flux m g , its exchange phase is given by an N th root of Ω(g), which can be straightforwardly calculated from the M matrix (Eq A15).

Appendix B: The symmetry action in SSPTs
Let us briefly discuss how symmetries may act anomalously on the edges in an SSPT. Again, let us consider only virtual edges, as in our earlier 2D discussion.
We first discuss the case for 3-foliated phases. Let us take a cubic subregion M , and consider applying symmetry operations restricted to this subregion.
An xy planar symmetry restricted to this region, S However, there is more information contained in U h M (g) that is missed in this process. We know that U h M (g) must commute with all untruncated symmetries, such as S (xy) (z; g), when acting on the ground state. That is, U h M (g) has to be overall charge neutral under all symmetries. Consider the symmetry S (xy) (z; g ) The four intersection locations are spatially separated, thus one can sensibly define a charge on each of the four hinges, which do not have to be trivial. That is, S (xy) (z; g ) commuting with the first hinge may result in a phase e iφ1 , the second e iφ2 , and so forth, which is fine as long as e i(φ1+φ2+φ3+φ4) = 1. These charges pinned to the hinges cannot be removed by a symmetric local unitary transformation, and are therefore a sign of a non-trivial phase. Such charges arise due to the existence of 2D linear SSPTs: 2D phases with SSPT order protected by line-like subsystem symmetries 4 . We will later construct an explicit example of this.
Appendix C: Non-trivial SSPT phases with trivial H 3 Here, we highlight a mechanism by which an SSPT phase may be non-trivial, despite appearing trivial in our H 3 [G L , U (1)] picture along all planar directions (but still a weak phase overall). Let us begin with an example: the so-called semionic X-cube model (See Eq 14 and is a product of 8 Zs on the corners of a cube. We immediately notice that this model has a trivial H 3 [G L , U (1)], which one can confirm by simply noting that a planar symmetry actually acts trivially along the edges (except at a corner), and so cannot produce anything non-trivial under the Else-Nayak procedure. This is due to the presence of higher symmetry: this model actually is symmetric under line-like symmetries. For example, a product of Xs along the x direction, commutes with H sem . Acting on the ground state, the global symmetry operator truncated to a cube S M (g) acts by creating charges localized at its corners (as it must be for a 3D system with line-like subsystem symmetries). It is this pattern of charges which leads to the non-trivial lineon braiding phase.
Notice that there is no contradiction between the system having a trivial H 3 [G L , U (1)] classification and lineons having a −1 braiding statistic. This is due to the fact that the fundamental braiding process which H 3 [G L , U (1)] cares about is between lineon dipoles. Braiding two lineons in a plane z 0 − 1/2 is like braiding a stack of lineon dipoles on planes z < z 0 with another stack z ≥ z 0 . However, the braiding phase in a Z 2 theory is only defined modulo ±1, and so a braiding phase of −1 is the same as trivial from this perspective.
It is straightforward to show that the model described by H sem is weak. One can write the Hamiltonian as where U is a local unitary circuit consisting of CZ gates. The ground state is then where |ψ 0 is the trivial paramagnetic phase. It is possible to write U as where U z acts only between layers z and z + 1, and U z commutes with all planar symmetries (and each other) This is exactly the form of a planar-symmetric local unitary circuit (just a higher dimensional version of the linearly-symmetric local unitary circuit defined in Ref 4). Thus, U z |ψ 0 , where |ψ 0 is the trivial paramagnetic state X r |ψ 0 = |ψ 0 , describes a 2D phase on planes z, z + 1 (which is actually a 2D linear SSPT), and the ground state of H sem is simply the stack of these.
One of the consequences of the fact that H sem describes a weak phase is that there is no obstruction to constructing an SSPT phase which is described by H sem for z 0 and is completely trivial H triv for z 0. For example, we can define the Hamiltonian which is composed of commuting terms. If we look at the action of S M (g) where M is a large cube crossing z = 0, one finds that there is a single Z pinned to each hinge at z = 0, exactly as discussed in Appendix B. Indeed, by stacking 2D SSPTs, it is possible to realize phases with various choices of allowable charges pinned on each hinge. Finally, we note that we do not have a proof that all phases with a trivial H 3 [G L , U (1)] classification are weak. There may also exist other mechanisms by which a phase may be non-trivial. However, we are not aware of any counterexamples.

Appendix D: Mobility of single lineons
Here, we show that in the fracton dual of a 3-foliated phase, a "single lineon" need not actually be mobile along lines. In a slight abuse of nomenclature, we will still call this excitation the lineon, even though it may not be mobile along a line.
First, let us identify what is commonly referred to as the lineon. Again, consider the action of the global symmetry truncated to a cube, S M (g), which acts on the ground state as some unitary along the hinges h M , Then, the operator L M (g) ≡ S M (g)U † h M (g) acts trivially on the ground state. The gauged version of the operator L M (g) will define our lineon.
Let us review how the generalized gauging process 13,17,18 works for a 3-foliated model. For each xy plaquette, we define a plaquette variablẽ and similarly for yz and zx plaquettes. We may then write the Hamiltonian of any subsystem symmetric model in terms of these plaquette variables. In addition, we energetically enforce the constraints by adding terms to the Hamiltonian which projects onto this subspace. In terms of these plaquette variables, the symmetry S M (g) acts only on the hinges of the cube M . If the gauged version of U h M can be written such that it only acts along the hinges as well, then the gauged version of L M (g) = S M (g)U † h M also only acts along the hinges. In this case, if we truncate L M (g), we obtain an operator which creates a g lineon excitation at each of its truncated hinges. A single lineon is guaranteed to be mobile along a line, simply due to the fact that the operator L M (g) is line-like along the hinges. While this has been true in virtually all previously studied models, it is not true generally.
Indeed, consider H half from Appendix C. In that model, U h M (g) had a single charge Z r pinned at each place where h M crossed z = 0. However, there is no way to gauge U h M (g) in such a way as to keep the support of the operator only along the hinges. This means that one cannot construct a lineon which crosses the z = 0 plane alone.
From the perspective of the fracton order, we may consider a z-moving lineon at z > 0. Now, suppose we naively move this lineon down to z < 0, crossing the z = 0 plane. What one will find is that, upon crossing, there is a single fracton charge excitation stuck at the z = 0 plane. As it is a single fracton, which is immobile, one cannot simply move it along with the lineon (which would simply amount to a redefinition of the lineon for z < 0 vs z > 0). Thus, a z-moving lineon cannot cross the z = 0 plane without paying an energy penalty in the form of a fracton stuck at z = 0. Now, instead of having simply a single plane at z = 0 where charges are pinned, we can imagine constructing a model in which charges are pinned on every plane, or every other plane, for example. In this case, a single lineon moving would create fracton excitations as it moved along, which are unable to be annihilated or moved along with the lineon as a redefinition. A single lineon therefore cannot be moved along a line without creating additional excitations. However, a pair of lineon anti-lineon on adjacent planes (the gauge flux) is always guaranteed to be a planon.

Appendix E: Proof of constraints
Here, we prove the two constraints mentioned in the main text for 2 or 3-foliated models. Let us label by g z the group element g in the zth factor of G L , and g gl = z g z a global symmetry. Again, view the SSPT as a quasi-2D system with the large symmetry group G L . We will denote the representation of the symmetry g acting on the zth plane by simply S(g z ), rather than S (xy) (z; g) as in the main text. We take the system to also be symmetric under yz-planar symmetries.
Consider a square region M of this quasi-2D system (M would contain all sites x, y, z with x 1 < x < x 2 , y 1 < y < y 2 , and all z, for some choice of x 1,2 , y 1,2 ). The key fact is that the global symmetry truncated to M , S M (g gl ), acts trivially along the yz face of M (simply due to the fact that it acts identically to yz-planar symmetries). Thus, U ∂M (g gl ) acts trivially along the yz face. Now, we may perform the Else-Nayak procedure, further choosing a restriction to an open segment C which ends along the yz face. Going through the procedure with a trivial U ∂M (g gl ), we can always get ω(g gl , h gl , k gl ) = 1 for arbitrary g, h, k ∈ G (up to a coboundary).
This leads to our second constraint. Calculating the invariant Ω(g) = N n=1 ω(g, g n , g) for any global symmetry results in a trivial type-I cocycle with itself, Ω(g gl ) = 1. In terms of the M matrix, for each generator of G is exactly the global constraint from the main text. Next, consider the type-II cocycle between a global symmetry h gl and g z . This is given by the ratio We can calculate Ω(h gl g z ) using the Else-Nayak procedure, which we wish to show is simply equal to Ω(g z ).
First, note that if we have S M (g) defined on some larger M , which has a boundary action U ∂M (g), we may always use this to construct an edge action for S M (g) as which acts simply as S M (g)S † M (g) near ∂M , and has deferred all the non-triviality over to ∂M . Now, we may use this construction for U ∂M (g z ) in the Else-Nayak procedure, which, along the yz face, is equivalent to U ∂M (g z h gl ) (since U ∂M (g z ) = 1 is trivial along this edge). The procedure then continues, and since U ∂M (g z h gl ) is (by construction) invariant under conjugation by S M (h gl ), the process proceeds exactly the same regardless of whether we had chosen g z h gl or just g z . We can therefore always choose to have ω(g z h gl , (g z h gl ) n , g z h gl ) = ω(g z , (g z ) n , g z ) (E4) so that Ω(g z h gl ) = Ω(g z ), and therefore Ω II (g z , h gl ) = 1.
In terms of the M matrix, for h, g, being generators of G, is exactly the local constraint in the main text.
Appendix F: Various proofs for invariants F1 and F2 1. Independence of direction for F1 and F2, and triviality of F2 in 3-foliated model In this section, we prove the claims in the main text that 1) the invariants F 1 and F 2 must be the same regardless of which direction of planar symmetry we look at, and 2) that F 2 must be trivial in a 3-foliated model. We first introduce some ideas for a regular 2D SPT. First, let us make the simplifying assumption that U ∂M (g) is a purely diagonal operator. This is always possible to do in our class of models, where |ψ is an equal amplitude sum since if S M (g) sends {g r } → {g r }, then we may simply choose which one can verify satisfies U † ∂M (g)S M (g) |ψ = |ψ and will be only supported along ∂M as |ψ is symmetric. Note that although we have made this assumption, the spirit of our argument should remain the same even without it. In the Else-Nayak procedure, this means that U C (g) and V ∂C (g 1 , g 2 ) can also be chosen to be purely diagonal, and Eq A6 reads To measure Ω(g), consider the product which one can show using Eq F3 satisfies That is, the charge of Q a (g) under S M (g) is exactly the type-I invariant Ω(g). This procedure has the nice interpretation in the gauged language as measuring (half) the charge of N gauge fluxes m g . Next, consider a measurement of Ω II (g, h). One way to do so is by noting that we may use a region M 1 for S M (g) = S M1 (g), but instead consider a much larger region M 2 which fully contains M 1 for the symmetry S M2 (h), and also define S M (g n h n ) = S M1 (g n )S M2 (h n ) (formally, we would absorb some of the symmetry into U ∂M (h), like in Eq E3). Then, using the fact that U ∂M (g) commutes with all full symmetries S(h), and S M2 (h) ≈ S(h) when acting on U ∂M1 (g) since M 2 is much larger than M 1 , we have Next, one can always choose the truncation to a segment U ∂M1 (g) → U C1 (g) in a way that U C1 (g) also commutes with all full symmetries S(h), in which case as well. Using this choice, we have that Q a (gh) = Q a1 (g)Q a2 (h).
From this, one can readily compute the type-II cocycle And by symmetry, (note that these expressions are unambiguous since both numerator and denominator are diagonal). These have the nice interpretation on the gauged side of measuring the number of charges e h obtained as a fusion result of N gauge fluxes m g , or vice versa. Notice that while Q a (g) carries a charge under S M (g) and S(h), if we consider the contribution from the other endpoint of ∂C, Q b (g), then one must have that the phase factors cancel out from the two endpoints. This is simply due to the fact that the phase ω only appears when isolating V ∂C (g) to a single endpoint. Now, let us begin talking about the SSPT. Consider applying a symmetry S M (g) to a cubic region M , which acts non-trivially as U h M (g) along the hinges. Then, consider a symmetric truncation of U h M (g) → U C (g) which leads to V a (g 1 , g 2 ) in the Else-Nayak procedure, and consider Q a (x, z; g) on an upper hinge (see Fig 4), where we are now explicitly labeling the x and z coordinate of the hinge. Notice that if we had instead chosen to consider V a (g 1 , g 2 ) defined from the bottom hinge, we would end up with the conjugate Q * a (x, z; g) instead (as shown in Fig 4), which follows from the fact that the bottom hinge of S M (g) is related by a symmetry action to the top hinge of S M (g −1 ). Knowing Q a (x, z; g) is sufficient to calculate the invariants F 1 and F 2 .
Consider calculating F 1 , using H 3 [G L , U (1)] obtained from xy planar symmetries. Let us choose g to be the generator of G = Z 2N . Then, the invariant F 1 corresponds to for some arbitrary z 1 , with z 0 z 1 z 2 . Let us take the region M to be some region x < x 1 , such that the relevant edge is at x-coordinate x 1 . Then, applying Eq F5 For convenience, let us divide Q a (x 1 , z 1 ; g) into four quadrants, as shown in Fig 5, and denote its charge in each quadrant as Q ,Q ,Q , and Q . For example, Using this, we can express using Eqs F14 Alternatively, we could have used Eq F8 and Eq F9 to obtain and where Q ≡ Q Q , and similarly for others. Eq F16, F17, and F18 imply that the charge distribution in Q must satisfy Thus, there are two degrees of freedom for the charge distribution in Q, which we may call q 1 and q 2 , in which case Ω II (g < , g ≥ ) = e 2πi(q1+q2)/(2N ) . The invariant F 1 is then F 1 = q 1 + q 2 mod 2. Now, suppose we calculate the same quantity except using yz planar symmetries instead. We may perform the calculation using the same hinge Q a (x 1 , z 1 ; g), as shown in Fig 4. In this case, one finds that Ω (yz) where we have explicitly labeled everything with yz to avoid confusion (g (yz) < is the product of g (yz) x for x < x 1 , for example). In this case, one has F (yz) 1 = q 2 − q 1 mod 2. However, q 2 −q 1 = q 2 +q 1 mod 2, and so F (yz) 1 = F 1 is independent of whether we had chosen the xy or yz plane. In the 3-foliated case, we may use the same argument along a different hinge to show that F (zx) 1 is also given by the same quantity.
Next, consider the quantity F 2 . Take G = Z N × Z N , and choose g and h to be the two generators of G. Then, we wish to compute e 2πiF2/N = Ω II (g < , h ≥ )/Ω II (h < , g ≥ ) (F22) using the same set-up as before. Let us define to be the h charge in the quadrant of Q a (x 1 , z 1 , g), and similarly for the other quadrants. Then, using Eq F8 and Eq F9, (F26) such that is simply the total h charge of Q a (x 1 , z 1 , g). Clearly, if we were to perform this calculation for the yz plane using this same hinge (x 1 , z 1 ), we would find exactly the same result, e 2πiF (yz) 2 /N = Q h,g . Thus, F 2 is independent of whether we measure using the xy or yz planes.
Now, suppose our model is 3-foliated. We have shown that if we consider every endpoint (not just Q a ), the total charge must be zero under any untruncated symmetry (Eq F10). However, in a 3-foliated model, we may choose a symmetry operator which acts as a global symmetry near Q a , but does not act on the other endpoints at all (see Fig 5). This means that Q a must have trivial total charge under any symmetry. Thus, e 2πiF2/N = Q h,g = 1 must be trivial.
On the gauged side this has a natural interpretation: for the gauged 3-foliated model, N lineons at (x 1 , z 1 ) (which are mobile in the y direction) cannot carry any charge under S (zx) (y; h), otherwise they could not have been mobile in the y direction in the first place.

Completeness and basis change
Here, we first go through how to obtain M for a 3D SSPT after stacking by a 2D SPT, as described by the main text. Then, we prove that the invariants F 1 and F 2 are a complete classification of all matrices M modulo this stacking. We prove this by showing that all possible M may be brought into a canonical form M canon , determined solely by F 1 and F 2 , via stacking 2D SPTs. Recall that we wish to stack a 2D SPT with symmetry group G 2D = G K , and that we do so by identifying each factor of G in G 2D with a plane z k in the SSPT.
That is, let k = 1, . . . , K label the factors of G in G 2D , which we associate with the plane z k , and M 2D (α,k), (β,k ) the KM × KM matrix characterizing the 2D SPT. To ensure locality, all {z k } must reside within some finite O(1) interval. Then, define M 2D to be the matrix with the same dimensions as M, whose elements are obtained directly from M 2D , and all other elements with z / ∈ {z k } zero. If we were stacking on to a 1-foliated model described by M, we would simply modify M → M + M 2D . However, when stacking to a 2-or 3-foliated model, we instead define the 2D SPT in terms of d r degrees of freedom. Thus, one For example, suppose G = Z N and we have a single type-I cocycle on plane z 1 , where we show only the {z 1 , z 1 + 1} submatrix. Then, within this submatrix, and all other elements outside of this submatrix are 0. This therefore results in two type-I cocycles valued 1 (recall that the diagonal elements are M ii = 2p i I ) on planes z 1 and z 1 + 1, and a type-II cocycle valued −2 between the two planes. This is one of the examples shown in Fig 1 of the main text.

b. Completeness
Next, let us show completeness of the invariants F 1 and F 2 . In a general group G = M α=1 Z Nα , we may define F α 1 for each even N α , and F αβ 2 for each N αβ = 1. These are defined in reference to some plane z 0 , Our strategy is as follows: assume we are given a general M, with some longest range coupling d max , defined as the maximum |z − z | where M (α,z),(β,z ) is non-zero. We then show that by stacking M 2D we can reduce M down to one in which d max = 1, such that there is only couplings between planes z and z ± 1. Then, we finally reduce M down to some canonical M canon , which only depends on F α 1 and F αβ 2 . Thus, any two M with the same F α 1 and F αβ 2 can be related to one another by stacking various M 2D , and are therefore a complete set of invariants.
First, suppose we have some matrix M with some longest range coupling d max > 1 (which is always O(1) due to locality). This means there is some element M (α1,z1),(α2,z2) = 0 where |z 2 − z 1 | = d max . By symmetry of M, we may consider z 2 > z 1 without loss of generality.
Suppose α 1 = α 2 , then take where M 2D α1α2 is viewed as a matrix indexed by z, with fixed α 1 , α 2 , and we only show the relevant non-zero submatrix. We then have which has a single −1 as its ((α 1 , z 1 ), (α 2 , z 2 ))th element (along with its symmetric partner), and all other elements are of |z − z | < d max . Thus, we may take which now has M (α1,z1),(α2,z2) = 0. Note that although in writing the submatrix we have assumed d max = z 2 −z 1 > 2, this also works for d max = 2. If α 1 = α 2 , then we may instead use again only has a −1 as its ((α 1 , z 1 ), (α 2 , z 2 ))th element, and all other elements have range smaller than d max . We may repeat this on all non-zero elements of M with distance d max , after which we end up with some matrix with d max < d max . We can repeat this process until we have d max = 1, meaning M α1α2 is a tridiagonal matrix.
Let us now define a canonical form M canon , for a given set of F α 1 and F αβ 2 , by and for α < β. For β < α, we simply have M canon αβ = (M canon βα ) T . We have also simply set F α 1 = 0 for any odd N α , and F αβ 2 = 0 for any N αβ = 1. The strong examples shown in Fig 1 of the main text are both already in canonical form. We will now show that our tridiagonal M can always be brought into its canonical form.
First, for each α, examine the symmetric matrix M αα . Consider each 2 × 2 block coupling z 1 and z 1 + 1. We may stack with which we can add to M to modify the offdiagonal element to be 0 or 1 depending on its parity (if N α even) or 0 (if N α odd). We may do this for all the offdiagonal elements, bringing them all to F α 1 . The diagonal elements are automatically constrained by the local constraint (Eq E5) to be −2F α 1 . Next, we may do a similar thing to M α,β for each α < β. In this case, we stack which we can use to eliminate all the lower-diagonal elements M (α,z1+1),(β,z1) . Then, the upper-diagonal elements are M (α,z1),(β,z1+1) = F αβ 2 and the diagonal elements are all automatically fixed by the local constraint to be −F αβ 2 . We have therefore brought an arbitrary initial matrix M, via moves of the form W T M 2D W (stacking 2D SPTs), to a canonical form which only depends on F α 1 and F αβ 2 . From this, we conclude that F α 1 and F αβ 2 are a complete set of invariants for M.

Appendix G: Strong models
In this section, we introduce two strong models. The first is the 3-foliated Type 1 strong model with G = Z 2 , which we write down in the form of a Hamiltonian. The fracton dual is a novel fracton model which we explicitly write down. The second is the 2-foliated Type 1 and Type 2 strong model with G = Z N × Z N , which we write down the ground state wavefunction |ψ for. We may consider the 2-foliated model as part of a model with two sets of 2-foliated symmetries, in which case the fracton dual is again a novel model with unusual braiding statistics between fractons. Alternatively, we may examine the fracton dual of a single 2-foliated model by itself, which results in a 2-foliated fracton phase, with non-trivial braiding statistics between gauge fluxes. To obtain strong models for more general groups G, one may simply identify Z 2 or Z N ×Z N subgroups of G, and define the model in terms of those degrees of freedom. The Hamiltonian will be written as a sum of terms of the form where Z p are products of Z on the four corners of a plaquette p, and F r ({Z p }) is some function of these variables near the site r. The planar symmetries will act as products of Xs along xy, yz, or zx planes. As F r ({Z p }) only depends on the combinations Z p which commutes with all planar symmetries, this Hamiltonian is explicitly symmetry respecting. The function F r ({Z p }) consists of 6 Z p , 12 S p , and 12 CZ p1p2 operators on various plaquettes, and an overall factor of i. Fig 6 shows the model on the dual lattice, where plaquettes are represented by bonds, and the site r is mapped on to the red cube. Careful calculation will show that [B r , B r ] = 0 and B 2 r = 1. This Hamiltonian is is therefore simply a commuting projector Hamiltonian, and as every term is independent (only B r can act as X r ) and there are the same number of terms as sites, H has a unique gapped group state |ψ and describes a valid SSPT. We found it simplest to write a small computer script to confirm these commutation relations (and to compute the wireframe operator later), rather than doing so by hand.
The wireframe operator (Fig 7) obtained as a product of B r over a large cube, when ungauged, gives the action of the symmetry on the hinges of the cube. One may confirm using the Else-Nayak procedure that this model has M matrix and all other elements zero. This therefore realizes the Type 1 strong phase shown in Fig 1 for Taking a product of the cubic terms Bc (Fig 6) results in a wireframe operator with support along the hinges of the cube. This wireframe operator is shown here for a 7 × 7 × 7 cube, where the action of the Z, S, and CZ operators precede the X (which acts along the red cube).
puting the charges Q defined in Appendix F, one finds Q = Q = −1 and Q = Q = 1.
The fracton dual of this model is defined on the square lattice with qubit degrees of freedom on the bonds. The Hamiltonian is given by where c represents cubes, B c is the operator shown in Fig 6, v represents vertices, and A µν v is the product of Zs along the four bonds touching v in the µν plane (the usual cross term from the X-cube model). B c consists of Xs along the cube (the cube term from the X-cube model) but with an additional phase factor depending on the Z state around it in the form of S, Z, and CZ operators. Note that while the ungauged operator B r squares to 1, the gauged operator B c does not square to 1, it instead squares to a product of A v operators.
This model has the same fracton charge excitations as the usual X-cube model. However, the lineon excitations are modified. To find out what they are, consider the product of B c over a large cube, c B c , shown in Fig 7. This results in an operator with support only along the hinges of the cube. This operator, when truncated, is the operator which creates lineon excitations at its ends.
From this, the crossing (braiding) statistic of two lineon can be readily extracted. Reading off of Fig 7, a pair of x-moving lineons on line (y 1 , z 1 ) is constructed by the operator x,y1,z1 S (z) x,y1,z1 CZ (x↔y) x,y1,z1 where X (x) x,y,z is an X on the bond originating from the vertex at (x, y, z) going in the positive x direction, and similarly for S (z) x,y,z , and CZ (x↔y) x,y,z is a CZ between Z (x) x,y,z and Z (y) x,y,z . L x creates two lineons at x 0 and x 2 . Meanwhile, a pair of y-moving lineons is constructed by which creates two lineons at y 0 and y 2 . Note that depending on which hinge of the wireframe we obtain L x and L y from, there may be additional Z operators, which correspond to a choice of lineon or antilineon (and will affect the braiding phase by a ±1). It can be readily verified that when these two operators cross (i.e. y 0 < y 1 < y 2 and x 0 < x 1 < x 2 ), they only commute up to a factor of i, using the relations XSX = iZS and X 1 CZ 12 X 1 = Z 2 CZ 12 . Thus, the braiding phase of any two lineons in this model is ±i.

2-foliated strong model
In this section, we describe a 2-foliated model which realizes both Type 1 and/or Type 2 strong phases. The F 1 invariants are therefore simply q 11 and q 22 modulo 2, and the F 2 invariant is −q 12 . By the proof from our previous section, the invariants will also the same for the yz symmetries.
But before we can conclude that we have constructed a strong phase, we must show that this state is symmetric under yz symmetries. The purpose of the second term in f (SSP T ) r is to ensure that this is the case. Let us examine how f r ({g r }) transforms under a yz planar symmetry which sends {g r } → {g (yz) g r }, or, on the relevant degrees of freedom, (g r , g r+y , g r+z , g r+y+z ) → (gg r , gg r+y , gg r+z , gg r+y+z ) (g r+x , g r+x+y , g r+x+z , g r+x+y+z ) unchanged (G17) A calculation shows that which simplifies to f r ({g (yz) g r }) f r ({g r }) = P (g r+z , g r+x+y+z , g r+y+z , g) P (g r , g r+x+y , g r+y , g) P (g r , g r+x+z , g r+z , g) P (g r+y , g r+x+y+z , g r+y+z , g) P (g r+y , g r+x+y , g r+y , g) P (g r+z , g r+x+z , g r+z , g) (G19) where P (g 1 , g 2 , g 3 , g) = exp If one considers the contribution from neighboring cubes, one finds that the factors of P (. . . ) exactly cancel out between neighboring cubes. Repeating this calculation for a yz-planar symmetry which transforms the other four sites in Eq G17, one finds the same result. Thus, the wavefunction is indeed symmetric under yz planar symmetries and describes a strong SSPT phase for a 2-foliated model. If one wished, one could confirm that the matrix M (yz) obtained from yz planar is also strong with the same F 1 and F 2 invariants, by following the Else-Nayak procedure. Obtaining a gapped local Hamiltonian corresponding to this ground state is straightforward, and is done in the same way as for the standard group cohomology models, Eq A21.
the two qubit degrees of freedom become one effective Therefore, we have demonstrated that the Hamiltonian H 3-fol , obtained by strongly coupling three mutually perpendicular 1-foliated gauge theories, is dual to a weak SSPT Hamiltonian H un 3-fol . Moreover, any weak 3-foliated SSPT can be constructed in this way, since H un 3-fol describes a stacking of three 1-foliated SSPTs which can in principle by arbitrary. While our discussion has focused on the G = Z 2 case for simplicity, it can be straightforwardly generalized to arbitrary abelian group G.

Condensation transitions
As alluded to in the main text, and discussed in Ref. 8 , the procedure of stacking a 2D SPT onto a 3-foliated 3D SSPT is dual to the procedure of adding a 2D twisted gauge theory to a 3D twisted X-cube model, and condensing composite planon excitations composed of fracton dipole and 2D gauge charge pairs. This planon condensation process has the effect of confining the lineon dipoles and 2D gauge fluxes that braid non-trivially with these fracton dipoles and 2D gauge charges respectively, leaving deconfined only the composites of lineon dipoles and 2D gauge fluxes, which become the lineon dipoles of the condensed phase. The result is that the statistics of these lineon dipoles are now modified by the addition of the 2D gauge flux statistics.
Let us consider a simple example. Consider stacking an xy oriented 2D Z 2 SPT, with M matrix a single entrymatrix equal to 2, between layers z = 0 and z = 1 of a trivial 3D 3-foliated Z 2 planar SSPT, dual to a copy of the X-cube model. This procedure is dual to adding a 2D double semion layer, whose gauge flux has semionic exchange statistics, and condensing the planon composed of the 2D gauge charge plus the fracton dipole centered around z = 1/2. This condensation has the effect of confining both the 2D gauge flux and the lineon dipoles centered around z = 0 and z = 1, and leaving deconfined the composite of the z = 0 lineon dipole and the 2D gauge flux, and the composite of the z = 1 lineon dipole and the 2D gauge flux. Both of these composites therefore obtain semionic exchange statistics. This procedure corresponds to the addition of a single self-loop to a 3-foliated SSPT in our graphical notation (see Fig. 1).
This dual picture interpretation of the stacking construction of weak 3-foliated SSPTs sheds light on the correspondence with p-string condensation. The key point is that the p-string condensation procedure resulting in untwisted and twisted X-cube models, and the planon condensation procedure outlined above, commute with one another because they both involve condensation of pure gauge charge. Therefore, one can construct the dual phases of weak 3-foliated SSPTs by 1) starting with three decoupled stacks of 2D toric code layers, 2) adding 2D twisted gauge theory layers and identifying the added gauge symmetries with existing gauge symmetries by condensing pairs of gauge charges, and 3) driving a p-string condensation transition. Since all 1-foliated SSPTs are weak, and thus can be constructed by stacking 2D SPTs, step 2 allows for the creation of arbitrary 1-foliated gauge theories. Therefore, any fracton model which is obtained by performing p-string condensation on intersecting 1-foliated gauge theories, is dual to a weak 3-foliated SSPT.