Foliated fracton order from gauging subsystem symmetries

Based on several previous examples, we summarize explicitly the general procedure to gauge models with subsystem symmetries, which are symmetries with generators that have support within a sub-manifold of the system. The gauging process can be applied to any local quantum model on a lattice that is invariant under the subsystem symmetry. We focus primarily on simple 3D paramagnetic states with planar symmetries. For these systems, the gauged theory may exhibit foliated fracton order and we find that the species of symmetry charges in the paramagnet directly determine the resulting foliated fracton order. Moreover, we find that gauging linear subsystem symmetries in 2D or 3D models results in a self-duality similar to gauging global symmetries in 1D.


Introduction
Gauging is a powerful tool in the study of gapped quantum phases with global symmetry. When gauging the global symmetry of a system, gauge fields corresponding to the symmetry group are added to the system so that the global symmetry can be enhanced to a local symmetry. It is useful to consider such a procedure because different phases under global symmetry map into different phases of the gauge theory. Symmetric (e.g. paramagnetic) phases map into deconfined gauge theories while symmetry breaking phases map into a Higgsed gauge theory. Different symmetry protected topological (SPT)/symmetry enriched topological (SET) phases map into different deconfined gauge theories with different statistics among the gauge fluxes (see, e.g., Refs. [1,2]).
Recently, it has been realized that a similar gauging procedure can be applied to systems with subsystem symmetries as well [3][4][5][6][7][8][9]. Subsystem symmetries are symmetries with generators that act non-trivially only on a sub-manifold of the system. After gauging, the system is mapped to a model with 'fracton order' [3,. This relation has been demonstrated for various classical/quantum spin models, stabilizer codes, domain-frame condensate models, etc. In this paper, we summarize and make explicit the general gauging procedure. That is, we describe explicitly a systematic procedure for gauging models with subsystem symmetries which can be applied to any local quantum model with such symmetry. In particular, the gauge fields are added at the center of 'minimal' coupling terms which are not on-site symmetric and which generate all other non-on-site-symmetric coupling terms. A modified Hamiltonian can then be written with enhanced local symmetry and with dynamical terms for the gauge field, which defines the gauge theory. We focus on abelian symmetry groups only in this paper.
The next key question is: what is the relation between the ungauged order under subsystem symmetry and the gauged fracton order? To address this question, we study the mapping between ungauged and gauged phases (several of these examples have been studied in the previous literature [3][4][5][6][7][8][9]) and propose a way to interpret the correspondence. In 2D and 3D, gauging linear subsystem symmetries (which act on 1D lines) maps paramagnetic (trivially symmetric) phases and symmetry breaking phases into one another, while subsystem symmetry protected topological (SPT) phases [6] may map into themselves. This is similar to the case of global symmetries in 1D, where paramagnets are mapped into symmetry breaking phases, and SPT phases can map into SPTs. In 3D, gauging planar subsystem symmetries leads to foliated fracton order, as defined in Refs. [15,28]. In particular, symmetry charges that transform under planar symmetries in one, two or three directions map directly to planon, lineon and fracton charge excitations, which are restricted to move only in a plane, along a line, or which cannot move at all. The restricted motion of the charge excitations in the fracton model hence originates from the requirement to preserve subsystem symmetries in the ungauged model. By counting the species of symmetry charges in the ungauged model, we can make direct connection to the foliated fracton order after gauging. For example, it was shown in Ref. [3] that gauging the (paramagnet phase of) the plaquette Ising model and the tetrahedral Ising model results in the X-cube and the checkerboard model respectively. By counting symmetry charges, we can see that the checkerboard model should be equivalent to two copies of the X-cube model. We present the mapping between the two in Ref. [38] and in section 4, we explain how counting symmetry charges leads to the same conclusion. Given the analogous foliation (or layered) structure in 3D models with planar subsystem symmetry and 3D foliated fracton phases, there is a natural correspondence. As shown in Fig. 1, for 3D models with planar subsystem symmetry, to increase the system size by one lattice spacing in the direction of one set of planar subsystem symmetries, it is necessary to add degrees of freedom (DOFs) on an entire plane and increase the number of generators of subsystem symmetries by one. The added planar subsystem symmetry acts as a global symmetry on the added plane. On the other hand, as we discussed in Ref. [15,28], for 3D foliated fracton phases, to increase the system size by one lattice spacing along one of the foliation axes, it is necessary to add a layer containing a gapped 2D topological state as a resource. Thus, it is natural that subsystem symmetry symmetric states gauge into foliated fracton models since the added layer gauges into a deconfined 2D gauge theory with gapped topological order.
The paper is organized as follows: In section 2, we briefly review the procedure of gauging global symmetries using as an example the 2D paramagnetic state. Section 3 then discusses the generalized gauging procedure that can be applied to systems with subsystem symmetries in a systematic way. Multiple examples (including examples that have appeared in the previous literature) are discussed to show how the procedure works in different situations. Section 4 studies the correspondence between phases with subsystem symmetries and the phases of their gauged theories through multiple examples and the result is summarized in Table 1 in section 5.

Review: Gauging global symmetry
First, we give a brief review of the procedure for gauging global symmetries (for more careful discussions see, e.g., Refs. [1,39]). We consider the simplest example: the transverse field Ising model with global Z 2 symmetry, coupled to a Z 2 gauge field. The Hamiltonian takes the simple form of where the σ's are Pauli matrices on each lattice site (blue dots in Fig. 2) and vw denotes nearest neighbor pairs. The system has a global Z 2 symmetry of U = v σ x v . To couple the model to a Z 2 gauge field, we introduce gauge field degrees of freedom τ on each link of the lattice (green dots in Fig.2). τ x corresponds to (the exponential e iE of) the 'electric field' of the gauge field and τ z corresponds to (the exponential of) the 'vector potential' of the gauge field. The local symmetry, or the Gauss's law, is given by e where the product is over all edges e with v as one end point. Next, we couple H to the gauge fields such that the new Hamiltonian is invariant under the local symmetry transformations A v . The transverse field terms σ x i are already invariant under the local symmetries, so we do not need to modify them and simply include them in the new Hamiltonian. The Ising coupling terms σ z i σ z j need to be replaced with σ z i τ z ij σ z j in order to be gauge invariant (i.e. commute with the A v term). Besides that we add the vertex term A v = σ x v e v τ x e at every vertex v to enforce gauge symmetry (Gauss's law) and B p = e∈p τ z e , where the product is over all edges around a plaquette p, to enforce the zero flux constraint on every plaquette. The total Hamiltonian then reads When J z = 0, the Ising model H is in the symmetric paramagnetic phase. After gauging, it maps to the deconfined phase of the Z 2 gauge theory. This can be seen by noticing that when the energy of the v σ x v term is minimized, the gauged Hamiltonian reduces to which is exactly the toric code Hamiltonian representing the deconfined phase of the Z 2 gauge theory. The low energy excitations include a bosonic gauge flux, which corresponds to the violation of one e∈p τ z e term, and a bosonic gauge charge, which corresponds to the violation of one e v τ x e term. These two excitations can be created with string operators shown in Fig. 2b-c. They braid with each other with a phase factor of −1, which is the Aharonov-Bohm phase factor in the Z 2 case.
When J x = 0, the Ising model H is in the symmetry breaking ferromagnetic phase. After gauging, it maps to the Higgsed phase which lacks non-trivial topological order. This can be seen by noticing that when J x = 0, H g has a unique ground state and no fractional excitations.
This gauging procedure can be applied to any local quantum Hamiltonian on any lattice satisfying a global symmetry G by introducing gauge fields on the links of the lattice, enforcing gauge symmetry (Gauss's law), modifying interaction terms to be gauge invariant, and finally including a flux term for the gauge field. By doing so, we obtain a gauge theory of group G. The properties of the gauge theory can be determined from the ungauged model in the following ways: 1. If the symmetry is spontaneously broken in the ungauged model, then the gauge theory is Higgsed with trivial topological order. 4. The braiding statistics between a gauge charge and a gauge flux is independent of the original order; it is given by the Aharonov-Bohm phase factor, which is determined by the symmetry group. For example, in a Z N gauge theory, the phase factor between an elementary charge and an elementary flux is e i2π/N .

5.
In 1D, gauge theories are not topologically ordered. Symmetry breaking and trivial SPT phases map into each other upon gauging. Non-trivial SPT phases can map to themselves upon gauging. (We briefly review the gauging of 1D phases in appendix B.) 3 Gauging subsystem symmetry: general procedure How do we gauge models with subsystem symmetries? The simplest example of a system with subsystem symmetry is an Ising paramagnet on a cubic lattice (corresponding to the plaquette Ising model in Ref. [3]). Consider a cubic lattice with spin 1/2 degrees of freedom at each lattice site (blue dots in Fig. 3). The Hamiltonian is simply given by H = v σ x v . This Hamiltonian is invariant under planar subsystem symmetries where P XY n labels the XY plane with Z direction coordinate n and similarly for P Y Z n and P ZX n . Throughout this paper, we use X, Y , Z to label spatial directions and x, y, z to label spin directions. The red vertex is involved in the twelve minimal coupling terms highlighted by red squares. The gauge symmetry term is a product of a σ x at the red sphere and twelve τ x on the green spheres. (c) The product of four minimal coupling terms around the four blue plaquettes is the identity. The corresponding flux term is a product of four τ z on the green spheres.
This model (with additional plaquette terms) was originally considered in Ref. [3]; however, we are not including the Ising coupling term here for simplicity of discussion. To gauge it, Ref. [3] proposed to add a gauge degree of freedom τ at each face-center of the cubic lattice (green dots in Fig. 3). The gauge symmetry is then given by which is the product of a symmetry charge σ x v at a site v and the (twelve) electric gauge fields τ x f on the neighboring faces f . The gauge flux terms, which are minimal pure vector potential terms that satisfy the gauge symmetry, now involve the product of four τ z 's as shown in Fig. 3. The gauged Hamiltonian takes the form Since the symmetry charges are fixed by the transverse field σ x (in the ground state), the zero temperature phase of the gauged Hamiltonian becomes equivalent to that of the X-cube model [3]. However, for generic systems with subsystem symmetry the degrees of freedom may be located at different places in the lattice and may transform under the subsystem symmetry in different ways. For example, in Ref. [3], an example was discussed where the ungauged model contains DOFs at the vertices and at the face centers of a cubic lattice, where the subsystem symmetry acts on planes with integer and half integer coordinates (in units of the cubic lattice constant). Ref. [8] discussed an example where the DOFs lives both at the vertices and body centers; the ones at vertices transform under subsystem symmetry in one direction only. For a generic configuration of lattice structure and DOFs, where should the gauge fields be added and how should the gauge symmetry of the gauged model be defined?

General procedure
We will now outline a gauging procedure that is consistent with the gauging procedure for global symmetry [1,39] and various previous works for gauging subsystem symmetries. The input to the procedure is a lattice of degrees of freedom (in a Hilbert space), a set of symmetry operators, and a model H = h that is symmetric under the symmetry. We will focus on abelian groups only in this paper.
Suppose that the on-site symmetry charge at each site is measured by σ x v (in general the charge does not have to be a Z 2 charge, although we use the σ notation without loss of generality). The procedure is as follows: 1. Find the minimum coupling terms c that a) are not on-site symmetric; b) are a tensor product of operators carrying elementary symmetry charges at each site; and which, c) together with on-site symmetric terms, can be composed into any coupling term satisfying the symmetry. (Note that these minimum coupling terms are not necessarily included in the Hamiltonian; they are used only to locate the gauge degree of freedom in the next step.) 2. Assign a gauge degree of freedom τ c at the center of each minimum coupling term.
(τ x c can be thought of as the exponential e iE of the electric field E, while τ z c is the exponential of the vector potential. τ can be a general gauge field, not just a Z 2 one.)

The gauge symmetry is given by
where the product is over all minimum coupling terms c that contain v.
4. All symmetric coupling terms h can then be made into gauge symmetric terms h g by multiplying each minimal coupling factor in h by a τ z c .
5. The minimum coupling terms will usually not be independent of each other. Or sometimes, gauge fields are added for non-minimum coupling terms as well. In such cases, we then find independent minimum sets C of coupling terms c ∈ C whose product is either the identify or a product of on-site symmetric terms σ x . 1 Correspondingly, the product B C = c∈C τ z c becomes the flux term of the gauge field if it is a local term.
In this way, we can gauge a model H = h with global or subsystem symmetry into a gauge theory Note that a large part this procedure, such as determining the gauge degrees of freedom and the gauge symmetry, is completely independent of the original Hamiltonian and depends on the action of the symmetry operators on the ungauged Hilbert space. The only step that depends on the original Hamiltonian is step 4 where the original Hamiltonian is made gauge symmetric.
Let us consider some examples to see how this works.

Example: global symmetry
For global symmetry, the minimum symmetric coupling term is a nearest neighbor two-body term of the form O i O j where O i carries charge e and O j carries charge −e. Other symmetric coupling terms, including non-nearest-neighbor two-body terms and multi-body terms, can all be constructed as composites of the nearest-neighbor two-body terms and on-site symmetric terms. Therefore, the gauge DOFs are assigned to each link of the lattice. The gauge symmetry term involves one lattice site and all the emanating links. The set of two-body terms around the same plaquette combine into on-site symmetric terms; therefore we have one flux term per plaquette. This is exactly the gauging procedure we reviewed in Sec. 2.
Changing the lattice structure corresponds to choosing a different set of minimum coupling terms, which does not affect the nature of the gauge theory obtained.

Example: 3D planar symmetry on a cubic lattice
For the subsystem symmetry example discussed above (DOFs at vertices of cubic lattice, transforming under planar symmetry in three directions), the minimum symmetric coupling term is the four-body plaquette term v∈p σ z v , as shown in Fig. 3a. All other symmetric coupling terms can be obtained as composites of such plaquette terms and on-site symmetric terms. Therefore, as suggested in Ref. [3], we can add one gauge field per plaquette. Each vertex is involved in 12 minimum coupling terms; therefore the gauge symmetry term is a product of one σ x and twelve τ x (Fig. 3b). Four minimum coupling terms around the same cube combine into identity as shown in Fig. 3c; therefore we have the corresponding flux terms. This is exactly the gauging procedure we reviewed at the beginning of this section [Sec. 3].

Example: 3D planar symmetry on a FCC lattice
Consider the situation corresponding to the tetrahedral Ising model discussed in Ref. [3], as shown in Fig. 4. Besides the DOF σ v at vertices of the cubic lattice, there are DOF σ f at the faces of the cubic lattice. Subsystem symmetry acts on each XY , Y Z and ZX direction plane either with integer or half integer coordinates.: The minimum coupling terms, as shown in Fig. 4a, are the tetrahedral terms involving one σ z v and three σ z f 's. All other symmetric coupling terms, including four-body terms of σ z v 's and four-body terms of σ z f , can all be constructed from this minimum coupling term. Therefore, as discuss in Ref. [3], one gauge DOF τ is to be assigned to each tetrahedron. The gauge symmetry terms are the product of one σ x together with the eight τ x 's in the surrounding tetrahedrons (Fig. 4b). The product of the same eight tetrahedron minimum coupling terms also happens to be the identity; therefore, we have the product of eight τ z 's as the flux term (Fig. 4b). If the σ DOF are all polarized by −σ x , the gauged model becomes exactly the same as the checkerboard model. The red vertex is involved in the eight minimal coupling tetrahedron terms centered at the green spheres. The gauge symmetry term is thus a product of a σ x at the red sphere and eight τ x on the green spheres. The product of the eight minimal coupling tetrahedron terms is the identity. The corresponding flux term is a product of eight τ z on the green spheres.

Example: 3D planar symmetry on a BCC lattice
Now consider the situation described in Ref. [8], where there is one DOF σ 0 at each cube center and three DOFs σ a , σ b , σ c at each vertex, as shown in Fig. 5. σ 0 transforms under subsystem planar symmetries in all three directions while σ a , σ b , and σ c transform only under symmetries in Y Z, ZX, and XY planes, respectively. An XY -plane subsystem symmetry generator is a product of all σ x 0 in a particular XY plane (P XY m+1/2 ) with Z coordinate m + 1/2 and all σ x c in the two neighboring XY planes (P XY m and P XY m+1 )with Z coordinate m and m+1: U Y Z and U ZX are defined in similar ways. The minimum coupling terms are the triangular terms shown in Fig. 5a. All other symmetric coupling terms can be composed from these minimum coupling terms. Therefore, to gauge the model, we need to assign one gauge DOF τ per triangle. The gauge symmetry terms are then the product of one σ x 0 with 24 τ x 's around it (Fig. 5b), and the product of one σ x a (or σ x b , σ x c ) with four τ x 's around it (Fig. 5c). The product of four triangular coupling terms is the identity, therefore we have the product of the corresponding four τ z 's as the flux term (Fig. 5d). This is the minimum gauging scheme for such a distribution of symmetry charges.
We could add gauge fields corresponding to non-minimum coupling terms as well. This is what was done in Ref. [8], where a gauge field is added for each four-body plaquette coupling term of the σ 0 's. Since this four-body term can be obtained by composing two triangular terms, this results in one more type of gauge flux term corresponding to the product of the τ z associated with these three coupling terms (one plaquette and two triangular terms).

Correspondence before and after gauging
Using the general gauging procedure, in this section we are going to explore the correspondence between models with subsystem symmetry (before gauging) and the gauged model with (potential) foliated fracton order. We refer to such a correspondence as the 'gauging correspondence'. While the following discussion is mostly based on specific examples, we expect several features of the gauging correspondence to apply generically, as specified below. In Appendix A, we will also show that the gauging procedure can be applied to global dipole conservation symmetries to produce a symmetric tensor gauge theory.

Planar symmetry and foliated fracton order
First, let's discuss models with subsystem planar symmetries. We are going to study models of increasing complexity -paramagnets with subsystem planar symmetries in one direction, two directions, three directions and four directions respectively -as well as models where the symmetries are spontaneously broken. We expect the following features to be generically true in the gauging correspondence: 1. when the planar symmetries are spontaneously broken, the gauged model does not have nontrivial order; 2. when the planar symmetries are not spontaneously broken, the gauged model has foliated fracton order 3. symmetry charges transforming under planar symmetries in one direction, two directions, and three or more directions turn into planon excitations, lineon excitations, and fracton excitations respectively upon gauging. The first feature is analogous to the Higgs mechanism in usual gauge theories. For the second one, we gave an intuitive understanding in the introduction section. Let us briefly discuss the third one before moving onto examples.
In Ref. [53], we proposed to characterize fractional excitations in foliated fracton phases using the notion of quotient superselection sectors (QSS). In particular, two fractional excitations are considered as equivalent (i.e. they belong to the same QSS class) if they differ only by local excitation and planons -a fractional excitation that moves in a 2D plane. Among the foliated fracton phases that we have studied, there are two types of QSS: 1) fracton sectors where the fractional excitation is fully immobile as an individual quasiparticle, and 2) lineon sectors where the excitation can only move along a straight line. In terms of the gauging correspondence, it is easy to see how the fracton/lineon QSS can emerge after gauging subsystem symmetries. Before gauging, if a symmetry charge transforms under planar subsystem symmetries in three directions, then to preserve subsystem symmetry, this charge cannot move freely in any direction. It is pinned at the intersection point of the three planes, as shown in Fig. 6, and such fracton symmetry charges have to be created four at a time. Upon gauging, they become the fracton gauge charges. If a symmetry charge transforms under planar symmetries in two directions, then this charge can move but only along the intersection line of the two planes. Such lineon symmetry charges become the lineon gauge charge upon gauging. Finally, if a symmetry charge transforms under planar symmetries in one direction only, then this charge can move along the plane. Such planon symmetry charges become the planon gauge charge upon gauging. Composites of fracton charges can become lineon or planon charges. For example, composing two Z 2 fracton charges in the same plane and displaced by a diagonal direction results in a lineon charge because the composite carries nontrivial symmetry charge in the two orthogonal planes only. By analyzing how the symmetry charges and their composites transform under subsystem symmetry, we can see how the gauging correspondence emerges. Let us see how this works through the following examples.

3D paramagnet with planar symmetry in one direction
We start with a simple and almost trivial case where the subsystem symmetry acts only in XY planes. Consider again the cubic lattice with DOF at vertices and the paramagnetic model The subsystem symmetry is given by Upon gauging, this model should naturally map to a stack of 2D (untwisted) deconfined gauge theories in the XY plane. The symmetry charges become the planon gauge charges in each 2D layer. The gauged theory is a trivial foliated fracton phase. Of course, this result does not depend sensitively on the lattice structure or details of the Hamiltonian, as long as the planar symmetries are preserved.

3D paramagnet with planar symmetry in two directions
A less trivial example is the 3D paramagnet H = − v σ x v with two sets of planar symmetries Each symmetry charge transforms under planar symmetries in two directions and hence becomes a lineon gauge charge upon gauging. The combination of two symmetry charges separated in the X or Y directions transform under planar symmetry in one direction only and hence is a planon. The combination of two symmetry charges separated in the Z direction does not transform under subsystem symmetry at all and hence is a not a fractional excitation. Therefore, in the gauged theory, we expect only one lineon QSS in the charge sector. This can be seen explicitly by applying the gauging procedure described in section 3. The two minimum coupling terms are 1) four σ z 's around a plaquette in the same XY plane (Fig. 7a), and 2) two σ z 's along the Z axis (Fig. 7b). Correspondingly, gauge fields are placed in each XY plane plaquette and on each link in the Z direction. The gauge symmetry term involves the product of one σ x v , four τ x XY 's and two τ x Z 's, as shown in Fig. 7c. The product of two plaquette coupling terms and four link coupling terms is identity, giving rise to the flux term as shown in Fig. 7d. The gauge charge, which corresponds to the violation of the gauge symmetry term, is a lineon that moves in the Z direction. It turns out that the flux excitation is also a lineon that moves in the Z direction. This is the anisotropic model introduced in Ref. [53].  The red vertex term in the center is included in four plaquette minimal coupling terms (red plaquettes) and two Z-axis terms (red edges). Therefore, the gauge symmetry term is a product of a σ x at the center (red sphere) and six τ x at the green spheres. (d) The flux term is a product of six τ z at the green spheres.

3D paramagnet with planar symmetry in three directions
Now let us consider the case where the planar subsystem symmetries lie along three directions. We have discussed the gauging procedure of three different cases (with different distributions of symmetry charges) in section 3. Now we will examine how the symmetry charge becomes a gauge charge through the gauging process and how the corresponding foliated fracton order emerges after gauging.

A. Cubic lattice
In the case discussed in section 3.3, where symmetry charges live at the vertices of a 3D cubic lattice and transform under planar symmetries in all three directions, each symmetry charge is a fracton and cannot move (since the charge is conserved on every plane). If two symmetry charges separated in the X, Y or Z direction are combined, then the composite transforms under planar symmetry in one direction only and hence is a planon. Therefore, upon gauging, the gauge charge sector of the gauge theory should contain only one quotient superselection sector -a fracton QSS. This is indeed the case for the corresponding gauge theory of X-cube model. As discussed in Ref. [53], the X-cube model contains three elementary QSSs: one fracton QSS and two lineon QSS. The one fracton QSS is the gauge charge sector of the gauge theory while the two lineon QSSs are the gauge flux sector of the gauge theory.

B. Cubic lattice: dual model
In fact, the X-cube model can be obtained through gauging a different model. Consider a 3D cubic lattice with two DOFs σ r and σ b (red and blue) at each lattice site. The red σ r transform under planar symmetry in XY and Y Z directions; the blue σ b transform under planar symmetry in Y Z and ZX directions; and their composite at each lattice site transforms under planar symmetry in ZX and XY directions. That is, the symmetries act as The minimum coupling terms are two-body terms σ z v,r σ z v+ŷ,r in the Y direction, two-body terms σ z v,b σ z v+ẑ,b in the Z direction, and four-body terms σ z v,r σ z v,b σ z v+x,r σ z v+x,r in the X direction, as shown in Fig. 8a. Therefore, according to the general procedure, a gauge field is added to each link of the cubic lattice. The gauge symmetry term is the product of σ x v,r (σ x v,b ) with four τ x on neighboring links in the XY plane (ZX plane), as shown in Fig. 8b-c. The combination of twelve minimum coupling terms around a cube is identity, therefore the flux term is the product of twelve τ z around a cube as shown in Fig. 8d.
If the σ spins are all polarized by Hamiltonian H = − v σ x v,r + σ x v,b , then the gauged model is exactly the X-cube model, but as the electromagnetic dual of the previous case. The symmetry charges transform under two planar symmetries, and therefore gauge into two independent lineon gauge charges (that move in the Y and Z directions). Their combination is a lineon charge that transforms under the XY and XZ planar symmetries and therefore moves only in the X direction. If two red charges separated in the X, Y , or Z directions are combined, then they form either a planon or a local excitation, and similarly for the blue charges. Therefore, the gauge charge sector contains two independent lineon QSSs. The gauge flux in this case makes up the fracton QSS.

C. FCC lattice
In the second case discussed in section 3.4, symmetry charges live both at vertices and face centers and transform under planar symmetry in all three directions. Again, each symmetry charge (both the vertex and face-center charges) is a fracton and cannot move. The combination of two vertex charges separated in X, Y , or Z directions transforms under planar symmetry in one direction only, and hence is a planon. Therefore, the vertex charge alone makes one fracton QSS after gauging. The combination of a vertex charge and a facecenter charge separated by half of a face diagonal transforms under two planar symmetries and are hence lineons. Similarly, the combination of two face-center charges separated by half of a face diagonal (out of the plane of the face) are also lineons. Taking into account neutral excitations -excitations carrying no symmetry charges -involving one vertex charge and three face-center charges, we can see that there are all together two independent lineon sectors. Therefore, upon gauging, the charge sector should contain one independent fracton QSS and two independent lineon QSSs. This corresponds exactly to the combination of the original and dual cubic lattice examples discussed above. Therefore, the gauged theory -the checkerboard model [3] -should be equivalent to two copies of X-cube model combined in a electromagnetic dual way. This is exactly what we show in Ref. [38].

D. BCC lattice
Now we come to the case discussed in section 3.5, where symmetry charges at cube center transform in three directions while symmetry charges at vertices transform in one direction only. The vertex charges are planon charges so they can be omitted when counting QSSs. The cube center charge is a fracton. Two fracton charges separated in the X, Y or Z direction combine into a planon. Therefore, upon gauging, the gauge charge sector contains only one fracton QSS. If the ungauged Hamiltonian is in the trivial phase (given for example by , then the gauged model would belong to the same foliated fracton phase as the X-cube model. In Ref. [8], a twisted version of the ungauged Hamiltonian is discussed. Upon gauging, the charge sector remains the same, while the flux sector may have different statistics compared to the X-cube model. Ref. [8] discussed the difference in statistics in terms of the self rotation of lineons. In Ref. [53], we show that this difference can be removed if 2D layers of twisted gauge theories are added to the 3D fracton model. Therefore, the gauged model has the same foliated fracton order as the X-cube model. Correspondingly, the difference between the twisted and non-twisted versions of the ungauged Hamiltonian can be removed by adding 2D layers of twisted SPTs. Therefore, the twisted ungauged model is equivalent to a 'weak SSPT', i.e. a stack of 2D SPTs, as defined in Ref. [6].

3D paramagnet with planar symmetry in 4 directions
It is also possible to construct a paramagnet in which every DOF transforms under a planar subsystem symmetry in 4 different directions. The model is constructed as follows: first, a lattice is constructed out of a fourfold foliation structure. To be precise, given four stacks of parallel planes such that no four planes intersect at a single point, a natural cellulation structure is defined in which each elementary 3-cell is a polyhedron bounded by these planes. Then, a σ DOF is placed in each 3-cell. The planar subsystem symmetries act on all 3-cells between neighboring parallel planes. The minimal symmetric coupling terms are the four-body terms v∈p σ z v with a σ z operator on each of the four 3-cells adjacent to a given edge (which is along the intersection between two planes). In the dual cellulation (or lattice), this edge is dual to a quadrilateral plaquette p, and the 3-cells are dual to vertices v. Upon gauging, the subsystem symmetric paramagnet defined on this type of lattice yields a generalized X-cube model as discussed in Ref. [16]. For example, using this type of construction, one can obtain the stacked kagome lattice X-cube model.

3D symmetry breaking state with planar symmetry
In all previous examples, for the ungauged model, we considered the simplest symmetric Hamiltonian of the form H = − v σ x v where the ground state is symmetric under all subsystem symmetries. For global symmetry, it is known that when the matter field undergoes spontaneous symmetry breaking, the gauge field is Higgsed and the gauge theory become nontopological. For subsystem symmetry, a similar Higgs mechanism applies, as first discussed in Ref. [3]. Let us repeat the exercise and see how Higgsing occurs in the cubic lattice example of section 3.3.
The minimum Ising coupling term that can be added to the system is the plaquette term involving four σ z 's (Fig. 3a). To make the term gauge invariant, we attach a τ z term in the middle of the plaquette. The total gauged Hamiltonian hence takes the form The B c terms are actually redundant for determining ground state because they can be composed out of the plaquette terms. Therefore, the Hamiltonian can be simplified into This is a cluster state [54] Hamiltonian where the σ and τ DOFs are connected through face diagonals. It has a unique ground state, and hence no topological or fracton order.

Linear symmetry and duality
Now let's consider subsystem linear symmetries in 2D and 3D models. We find that the gauging correspondence works in a very similar way to that of linear symmetries in 1D. It is well known (and we review it in Appendix B) that upon gauging the linear (global) symmetry in 1D, the gauged model also has an emergent global linear symmetry at low energy which comes from the zero flux constraint around the 1D ring. The gauging procedure leads to a duality between trivial symmetric paramagnets and symmetry breaking phases and a (self)-duality among nontrivial symmetry protected topological phases. From the examples discussed in this section, we find a similar correspondence in 2D and 3D with subsystem linear symmetries: 1. the model after gauging has linear subsystem symmetries at low energy which comes from the zero flux constraint around nontrivial loops; 2. symmetry breaking phases are mapped to trivial paramagnets; 3. trivial paramagnets are mapped to symmetry breaking phases; 4. non-trivial subsystem symmetry protected topological phases are mapped to non-trivial subsystem symmetry protected topological phases. We expect these features to apply generically to all models with linear subsystem symmetries.

2D paramagnet/symmetry breaking state with linear symmetry
It is possible for 2D systems to have linear subsystem symmetries. As we will see, gauging 2D systems with linear subsystem symmetries bears great similarity to gauging global symmetries in 1D. In particular, in both cases, trivial paramagnet and symmetry breaking phases are dual to each other through gauging. Consider a 2D square lattice with a σ DOF at each vertex. The subsystem symmetries acts along each row L X m and each column L Y n of the square lattice: which is a 2D cluster state model with unique ground state that is symmetric under the subsystem symmetries B X and B Y . Moreover, this state can be mapped to a symmetric product state through a symmetric local unitary transformation, indicating that it is equivalent to a trivial paramagnet. The symmetric local unitary is given by where is the controlled-X operation from a vertex spin to its neighboring gauge field and the Hadamard operator H = 1 1 1 −1 maps between σ x and σ z . When J = 0, corresponding to the trivial paramagnet phase before gauging, the gauged model is which can be reduced to if the −B x σ x v terms are all satisfied. This corresponds to the symmetry breaking phase of the gauge field under subsystem symmetries B X and B Y .

2D linear symmetry protected topological model
We now discuss an example of a 2D model with linear SSPT order, which is self-dual under gauging the subsystem symmetries. The system contains a σ DOF at each vertex of two interlocking square lattices labelled α and β. The linear symmetries act on all spins in a given row or column of either the α or β lattice. Explicity, the symmetry generators are As discussed in Ref. [6], the 2D cluster state model is a strong SSPT, which exhibits a protected edge degeneracy that grows exponentially with the length of the boundary. The Hamiltonian (also shown in Fig. 9) is where i(a), j(a), k(a), and l(a) refer to the four β lattice vertices neighboring vertex a, and vice versa for i(b), j(b), k(b), and l(b).
The minimal coupling terms satisfying the subsystem symmetry are the four-body terms around each elementary plaquette of either the α or β lattice. Thus, to gauge the model, gauge fields τ v are placed at every vertex v of both the α and β lattices (on top of each matter DOF), as shown in Fig. 9. The gauge symmetries then take the form Figure 9: The 2D cluster state model. The two stabilizer terms in Eq. (23) are circled in green above. The black and gray lattices are the α and β lattices. After gauging, gauge fields τ are placed on both the red and blue vertices.
in the previous example, there are no local gauge-symmetric flux operators; the only allowed flux terms act along an entire row or column: These operators correspond to symmetry generators of the gauge theory. Upon gauging the Hamiltonian takes the form This gauged model is actually a linear SSPT and is dual to the original SSPT. To see this, note that the matter DOFs can be decoupled from the gauge DOFs via the symmetric local unitary operator where as before, C σ V τ is the controlled-X gate from the vertex spin σ to an adjacent gauge field τ . Then which is a 2D cluster state model residing on the gauge DOFs. Here the relation H ∼ = H indicates that H and H have coinciding ground spaces and thus represent the same gapped phase.

3D models with linear subsystem symmetry
It is also possible for 3D systems to have linear subsystem symmetries. For example, suppose a system has a σ DOF at every vertex of a cubic lattice and symmetries which act along lines of spins along the X, Y , or Z direction. In this case, the minimal coupling terms that commute with the symmetries are eight-body terms v∈c σ z v involving the 8 qubits at the corners of a cube c. Therefore, to gauge such models, gauge fields are placed at the centers of each cube.
The correspondence before and after gauging of linear subystem symmetries in 3D bears similarities to the case of linear symmetries in 2D and global symmetries in 1D. For instance, the cubic Ising Hamiltonian is self-dual under gauging: the weak-coupling paramagnetic phase maps into the strongcoupling subsystem symmetry breaking phase and vice versa. Furthermore, the linear SSPT given by the the 3D cluster state Hamiltonian [6] is self-dual under gauging, in analogy with the 2D cluster state linear SSPT and the 1D cluster state global SPT.

Discussion
The gauging correspondence revealed in the previous examples is summarized in the table below. Fracton charges are acted upon by planar symmetry in three directions, whereas lineon charges are acted upon by planar symmetry in two directions. The fracton and lineon charges in the table are counted up to the attachment of planon charges, which are acted upon by planar symmetry in one direction only.

Before Gauging After Gauging Planar
One fracton charge X-cube with lineon flux symmetry Lineon charges in X, Y , Z directions X-cube with fracton flux in 3D One lineon charge in Z direction Anisotropic model with lineon flux Symmetry breaking Topologically / fractonically trivial state Linear Trivial paramagnet Symmetry breaking symmetry Symmetry breaking Trivial paramagnet in 2D/3D Non-trivial SSPT Non-trivial SSPT Table 1: Correspondence between phases with subsystem symmetries and gauge theory phases. The X-cube and anisotropic model listed refer to the corresponding foliated fracton phase, not to the specific model. Therefore, by counting the types of symmetry charges before gauging, we can determine the gauge charge and correspondingly gauge flux quotient superselection sectors in the gauge theory. A highly interesting and open question is whether there are non-trivial SPT phases with planar subsystem symmetry in 3D. The model discussed in Ref. [8] we now know to be equivalent to a weak SSPT. Hence upon gauging, it gives the same foliated fracton order as the X-cube model [53]. For a truly non-trivial SSPT, upon gauging, we expect the gauge charge and gauge flux to correspond to the same quotient superselection sectors while the gauge flux has non-trivial statistics compared to the X-cube model.

4)
The minimal coupling term can be made gauge symmetric by coupling it to a gauge field: ∂ a ∂ b φ → ∂ a ∂ b φ − A ab . 5) We now need to find linear combinations of the minimal coupling terms ∂ a ∂ b φ that result in zero. Equivalently, we want to find linear combinations of derivatives of A ab that are invariant under the replacement A ab → A ab + ∂ a ∂ b λ, which is often referred to as a gauge transformation. Thus, we want to find the smallest possible basis of gauge invariant operators, which is given by the magnetic tensor B i j = iab ∂ a A bj [40]. Therefore, gauging the matter Hamiltonian [Eq. (29)] results in the following gauged Hamiltonian (E ab ) 2 is added at the end since the above model is a gapless gauge theory. Traditionally, the (π − ∂ a ∂ b E ab ) 2 is not explicitly written, but is instead imposed as a gauge constraint or is considered irrelevant (under RG) at long length scales.

B Gauging global symmetry in 1D systems
In this section, we review the process of gauging 1D symmetric, symmetry breaking and SPT phases and see how symmetric and symmetry breaking phases map into each other upon gauging while SPT phases can map into themselves. Consider the 1D transverse field Ising model with Hamiltonian and global symmetry U = i σ x i . To gauge the model, we put gauge fields τ on every link. The gauge symmetry term is A i = τ x i−1,i σ x i τ x i,i+1 . The only flux term that satisfies all the gauge symmetries is a global term B = i τ z i,i+1 . Therefore, the flux term effectively becomes a Z 2 global symmetry of the gauged model.
Coupling H to the gauge field, we obtain the gauged Hamiltonian When J = 0, in the ground state, all the σ spins are polarized in the X direction and the gauge fields couple effectively through τ x i−1,i τ x i,i+1 . With respect to the effective global symmetry of B = i τ z i,i+1 , the gauge field ground state spontaneously breaks the symmetry. On the other hand, if B x = 0, the Hamiltonian becomes a 1D cluster state [54] model with unique ground state which is symmetric under the global B = i τ z i,i+1 symmetry. Now let us discuss an SPT example. Consider the 1D cluster state model This model has a global Z 2 × Z 2 symmetry generated by and the model has symmetry protected topological order under this symmetry [55].
To gauge the Z 2 × Z 2 symmetry, we put gauge fields τ between neighboring gauge charges. That is, we place one gauge DOF per site. The ones on the even sites are gauge fields of g 2 . The ones on the odd sites are gauge fields of g 1 . The Gauss law terms are The flux terms, which are pure gauge terms that satisfy the Gauss law, are They become the global Z 2 × Z 2 symmetry of the gauged model. To make the original Hamiltonian terms gauge invariant, we modify them to be Now the total Hamiltonian is All the terms commute, are independent, and are symmetric under the global symmetry. Therefore, on a closed ring, the ground state is unique. On an open chain, the terms no longer commute with the symmetry and need to be removed, leaving a two fold degeneracy at the edge as the symmetry protected edge state.