Non-Abelian fracton order from gauging a mixture of subsystem and global symmetries

We demonstrate a general gauging procedure of a pure matter theory on a lattice with a mixture of subsystem and global symmetries. This mixed symmetry can be either a semidirect product of a subsystem symmetry and a global symmetry, or a non-trivial extension of them. We demonstrate this gauging procedure on a cubic lattice in three dimensions with four examples: $G=\mathbb{Z}_3^{\text{sub}} \rtimes \mathbb{Z}_2^{\text{glo}}$, $G=(\mathbb{Z}_2^{\text{sub}} \times \mathbb{Z}_2^{\text{sub}}) \rtimes \mathbb{Z}_2^{\text{glo}}$, $1\to \mathbb {Z}_2^\text {sub}\to G\to \mathbb {Z}_2^\text {glo}\to 1$, and $1\to \mathbb {Z}_2^\text {sub}\to G\to K_4^\text {glo}\to 1$. The former two cases and the last one produce the non-Abelian fracton orders. Our construction of the gauging procedure provides an identification of the electric charges of these fracton orders with irreducible representations of the symmetry. Furthermore, by constraining the local Hilbert space, the magnetic fluxes with different geometry (tube-like and plaquette-like) satisfy a subalgebra of the quantum double models (QDMs). This algebraic structure leads to an identification of the magnetic fluxes to the conjugacy classes of the symmetry.


I. INTRODUCTION
The discovery of topological quantum phases of matter revolutionizes our understanding of quantum many-body phases and leads to many theoretical developments and remarkable experimental results. Fractional quantum Hall states [1,2] are the most studied topological quantum phases that host fractional excitations with non-Abelian braiding statistics, robust gapless edge modes, and ground state degeneracies. The longrange entanglement and intrinsic strong correlations in these topological quantum phases bring a key challenge of studying these phases.
The concept of the emergent "gauge" degrees of freedom (DOF) in the topological quantum phases provides the breakthrough of tackling these systems. For example, the notion of the flux attachment in the composite electrons gives rise to the dynamical gauge fields in the fractional quantum Hall states. The fractional braiding statistics can be understood from the knot structures of these dynamical gauge fields described by the topological quantum field theories (TQFT) [3][4][5][6][7]. The topological properties of these phases at low-energy are governed by the TQFTs. The structure of TQFT allows us to extract their mathematical structure and enable us to apply the language of the modular tensor category [8,9] to describe these systems. On the other hand, various exactlysolvable models of the lattice gauge theories [10] also describe the fundamental excitations of the topological quantum phases and provide insights from studying these models. The famous models are the Kitaev's quantum double models (QDMs) [11], which are the exactly-solvable limit of the lattice gauge theories based on a discrete group G on a two-dimensional lattice. The Kitaev's QDMs (the toric code is one of them) have fundamental excitations like anyons whose fusion and braiding statistics are characterized by the modular tensor category.
A generalization of the topological quantum phases in three dimensions is believed to have more exotic phases beyond current theoretical developments. The "fracton order", is one of the novel lattice models discovered recently . These three-dimensional lattice models share common features with * ricktu256@gmail.com † pychang@phys.nthu.edu.tw the two-dimensional counterparts such as ground state degeneracies on the torus and fractional excitations. However, some of the excitations, which are referred to as fractons, cannot be moved by local operators, and are intrinsically immobile. This unique feature has no counterparts in two dimensions. It is desired to develop new concepts and mathematical tools for understanding the properties of the fracton order. Furthermore, the fracton order is related to tensor gauge theories [41] and is proposed to be applicable for quantum memory [42,43].
A natural extension of current existing theories of fractons is the non-Abelian generalization. These non-Abelian fractons have potential applications to topological quantum computation and quantum memory. A variety of the non-Abelian constructions, including gauging bilayer/permutation symmetries [44,45], cage-net models [20], the layer constructions [21,46,47] and the topological defect networks [48,49] are proposed. The former construction relies on gauging a global Z 2 /permutation symmetry of the "Abelian" fracton models and the second one requires flux-string condensation. These models are restricted in a certain geometry (e.g., the cage-net model and the layer construction) or a specific layer construction (e.g, gauging bilayer fractons). Due to the restrictions from these proposals, a generalization to other non-Abelian fracton orders is not straightforward. To overcome this difficulty, we propose a general gauging procedure of constructing the non-Abelian fracton from gauging a pure matter theory on a lattice with a mixture of a subsystem symmetry and a global symmetry. This mixed symmetry can be a semidirect product of a subsystem symmetry and a global symmetry, or a non-trivial extension of them. The gauge principle on the lattice has tremendous success in constructing many models with topological order. In particular, for a threedimensional paramagnetic state with Abelian subsystem symmetries, the gauged model in the exactly-solvable limit is the Abelian fracton lattice models [50]. Following our gauge principle, one can obtain the (non-)Abelian fracton from gauging the (non-)Abelian versions of the mixed symmetry. We demonstrate several examples of the gauged Hamiltonian in the exactly-solvable limit that host (non-)Abelian fractons and other excitations. To the best of our knowledge, gauging the mixed symmetry on a lattice has not been discussed in the literature [51].
Comparing with gauging the bilayer/permutation fracton orders, which the non-Abelian fractons are constructed from (1) gauging subsystem symmetry of the matter field, (2) taking the exactly-solvable limit, and (3) gauging the global symmetry, our gauging procedure only requires one-step gauging. This one-step gauging process is very general and allows us to construct various lattice models that host (non-)Abelian fracton orders. We demonstrate four examples by gauging The former two cases and the last one produce the non-Abelian fracton orders. The third case is the non-trivial Abelian fracton order which has a hybrid structure [52].
Besides the non-Abelian fracton orders can be systematically constructed from our one-step gauging process, a direct analogy to the Kitaev's QDMs can be obtained. The excitations in QDMs are characterized by the local operators around a vertex v and a neighboring plaquette p. The algebra of these local operators is generated by A g , g ∈ G and B h , h ∈ G, where A g is the gauge transformation at v, and B h is the projector onto the subspace that a product of group elements around a plaquette p equals h. They satisfy The pure electric charges are classified by the irreducible representations (irreps) of G in that A g acts on the Hilbert subspace with a specific charge according to that representations, while B h = δ h,e (e is the identity) on this subspace. The pure magnetic fluxes are classified by the conjugacy classes of G in that the Hilbert subspace with a specific flux is generated by the state with B h = δ h,g , where g runs over all elements of that conjugacy class. In our one-step gauging procedure, the local gauge transformations can generate a local symmetry G local , which allows us to identify the electric charges (including non-Abelian fractons) with the irreps of G local . In addition, the magnetic excitations together with the local gauge transformations form a subalgebra of the QDMs in the constrained local Hilbert space. This algebraic structure allows us to identify the magnetic excitations with the conjugacy classes of G local . With the one-step gauging, we construct explicitly the underlying local G local = S 3 , D 4 , Z 4 , and Q 8 symmetries, by gauging → 1, respectively. We further identify the charges with the irreps of these symmetries, and the magnetic fluxes with the conjugacy classes of G local in these systems. These excitations are summarized in Tables I, II, III, and IV. This paper is organized as follows: In section II, we present the details of the general gauging procedure. In section III, we demonstrate the construction of (non-)Abelian fracton models by our gauging procedure. We show that the electric and magnetic excitations satisfy the algebraic structure of the Kitaev's QDMs in the constrained Hilbert space. The algebraic structure of the electric charges and the magnetic fluxes allows us to classify the (non-)Abelian fracton orders. We construct two new non-Abelian fracton orders by gauging G = Z sub 3 ⋊ Z glo 2 and 1 → Z sub 2 → G → K glo 4 → 1. Our gauging procedure is very general. It also reproduces the non-Abelian fracton in the gauged bilayer X-Cube code [44,45] and also the nontrivial Z 4 fracton order which has a hybrid structure discussed in Ref. [52]. Moreover, it can be easily generalized to more exotic mixed symmetries. One interesting example will be a generalization of the Z sub N ⋊ Z glo 2 , which is related to the continuous field-theoretic approach in Refs. [53][54][55]. Finally, we summarize our results in section IV.

II. GENERAL GAUGING PROCEDURE
We start with an "ungauged" system, which can be described as consisting purely of the matter DOF on lattice sites. The ungauged systems is described by H = H o + H n , where H o is the onsite terms and H n = − c J c c + H.c. is the nononsite couplings with coupling strengths J c . The ungauged Hamiltonian has a symmetry G such that the non-onsite couplings c are the minimal symmetric couplings, e.g., the nearest neighbor interactions and the plaquette interactions. Here, G can be a global symmetry, a subsystem symmetry, or a mixture of them. For a subsystem symmetry, G is the group generated by the operation on each subsystem. For example, the subsystem Z 2 spin-flip symmetry group of a cubic lattice with spin-1 2 on each site is denoted by G = Z sub 2 . The group Z sub 2 is not actually isomorphic to the group Z 2 , but is a large group generated by the spin-flipping on every shifted coordinate planes of the cubic lattice.
The gauging principle is to promote the symmetry G to be local. The famous example is gauging the U (1) symmetry. The matter field ψ(x) has a global U (1) symmetry ψ(x) → ψ(x)e iφ such that the Hamiltonian H = dxψ † (x)(−i∂ x )ψ(x) is invariant. While promoting the global U (1) to be local, ψ(x) → ψ(x)e iφ(x) , one need to introduce a gauge field a(x) to ensure the Hamiltonian is invariant. The corresponding gauge transformation of the gauge field is a(x) → a(x)+∂ x φ(x). We also add the magnetic field and electric field terms in H and promote the derivatives to the covariant derivatives. With this example in mind, we describe the general gauging procedure as follows: 1. For each of the G-invariant minimal couplings c in H n , we define a gauge DOF τ c with Hilbert space dimension equals the number of distinct eigenvalues of c. We label the computational basis of the τ c with respect to the eigenvalues of the minimal couplings [56]. These gauge fields τ c are placed on the links or plaquettes for the corresponding minimal couplings c. We define supp c to be the set of vertices that c acts on. When all the matter sites are in the same state such that the matter state is the simultaneous eigenstate of all c, the corresponding state of all the gauge DOF constitute a natural flat connection, and is denoted by |0 .
by g as the original symmetry transformation, but leave other matter DOF on other sites unchanged. This action leads to changes of the eigenvalues of the non-onsite couplings c. To ensure the system is invariant under the gauge transformation, the gauge DOF must compensate the changes of the eigenvalues of c. The gauge transformations need to satisfy This means that G is promoted to a local symmetry G local of the gauged system. We demonstrate a systematical way for construction the gauge transformation in Appendix A [57].
3. Construct the gauge flux: The existence of flux means that the gauge field is "twisted", which corresponds to the violation of some "relations" that should be satisfied when the connection is flat. Those "relations" come from the dependence between the couplings c in the ungauged system. That is, some of c can be combined (product or sum) to the identity. Similar to the "minimal" coupling, we can find the "minimal" relations consists of local combination of the c on a plaquette or a tube that produce the identity. We denote these relations by R r ({c}) = 1, where r is the label of the local combination. We define supp r to be the union of supp c for all c in the combination. Then the corresponding flux term is where R r ({Z}) is the relation with c replaced by the computational-basis operators (generalized Z) of the corresponding τ c , and P r ′ is the projector onto the B r ′ = 1 subspace [58]. These projectors must be included when some smaller r ′ is contained in the geometry of r to ensure that the final Hamiltonian is gauge invariant [59].
4. Construct the gauged non-onsite terms H n (τ ) = − c J c c(τ ) + H.c.: The gauged coupling is defined to be a gauge-invariant (commutes with all A v,g ) operator acting on all the matter and gauge DOF within the geometry of c, When the gauge DOF of a state is trivial, c(τ ) reduces to the ungauged coupling c: For a general state |α , c(τ )|α can be computed by expanding |α with the computational basis of the gauged DOF. For each term, we use a series of A v,g to rotate the gauge DOF of each term to be |0 . We apply (4) to the terms that succeed, and drop the terms that fail.
5. Construct the "electric field" terms H e : The terms are the dynamical part of the gauge fields, which are gaugeinvariant operators.
The final gauged Hamiltonian is and the physical Hilbert space has the constraint A v,g |physical = |physical . The gauged Hamiltonian H g has a generalization of the confined-deconfined and the Higgsed-deconfined phase transition [60]. The exactly-solvable limit of the gauged Hamiltonian H g can be obtained by eliminating the non-onsite couplings and hopping terms of the fluxes. (taking both coupling and the hopping strengths to be vanishing).
If G is a global or subsystem Z 2 symmetry, then the above procedure reproduces the gauging process of Ref. [50]. Also, if G is a global symmetry with the regular representation on the matter on a 2D lattice and we eliminate the matter DOF by choosing the unitary gauge, then the above procedure gives the Kitaev G-QDM introduced in Ref. [11]. The idea of a semidirect product of global and subsystem symmetry in the continuous quantum field theories appeared in Refs. [53][54][55].
We argue that, if G is a mixture of planar subsystem and global symmetries, then the resulting model contains both fractons and mobile particles (a hybrid order). The argument is as follows: if G contains a planar subsystem part, there must be some minimal plaquette-couplings; if G also contains a global part not generated from the subsystem part, then there must also be minimal link-couplings. Now the strings and membranes of the electric charges can be constructed from the minimal couplings: If we pull the two ends of a linkcoupling apart, or the four corners of a plaquette-coupling apart, the resulting operators are non-minimal couplings invariant under G. The gauged version of these non-minimal couplings are the string (from link-couplings) and membrane (from plaquette-couplings) for the charges. Having both types of couplings means that we have a hybrid order containing both mobile particles and fractons. Moreover, if G is not a direct product of a subsystem and a global symmetry, then the fusion rule of the charges will not be trivial. It is because that the corresponding QDM will not be the direct product of the two types of charges corresponding to the global and subsystem ones. The correspondence with QDM will be discussed in the following examples.
In the following, we focus on the H g for G = Z sub In the exactly-solvable limit, these gauged lattice models describe (non-)Abelian fractons.

III. CLASSIFICATION OF EXCITATIONS WITH THE QUANTUM DOUBLE MODELS
Before discussing the examples, we highlight our main results. From the generic gauging principle, the (non-)Abelian fracton order can be obtained by gauging a mixed symmetry with a subsystem symmetry and a global symmetry in the exactly-solvable limit. Here we consider the case on a threedimensional cubic lattice with a mixture of a global and a two-dimensional subsystem symmetry. While promoting the mixed symmetry G to be local, the "gauge fields" are introduced on links and plaquettes. The local gauge transformations generate the local symmetry G local . The magnetic flux operators defined from the "gauge fields" have two kinds of geometry, one is plaquette-like and the other is tube-like. The operators are denoted as B p and B t , where p and t are the label of the plaquette and the tube, respectively. The gauged Hamiltonian in the exactly solvable-limit hosts point-like and string-like excitations. The point-like excitations are created at the corners of the membrane operators or the ends of the string operators. On the other hand, the string-like excitations are created at the boundaries of the membrane operators. The point-like excitations correspond to the electric charges are the excitations of the local gauge transformation A u v,g . Here the superscript u indicates the unitary gauge which eliminates the matter DOF. v is the location of the vertex, and g is the symmetry operation. These excitations are local with respect to the vertex v, and can be specified from the local operators which form a representation of G local on the Hilbert space. Hence we can identify electric charges of the G-fracton model with the irreducible representation of G local .
For the magnetic fluxes, the excitations of B p are loop-like and the excitations of B t are point-like. Due to such difference in geometry, a constraint of the local Hilbert space is needed for making possible comparison. We consider the geometry of the plaquette p 0 and the tube t 0 adjacent to a chosen vertex v 0 : where the constituent of that geometry is indicated with dots. The magnetic fluxes with the constraints of the local Hilbert space can be classified from the algebraic relation as the QDMs in Eq. (1). Define the projectors on this cube adjacent to v 0 : where P p denotes the projector onto the subspace without B p flux at the plaquette it is drawn on. Imposing the constraint P side = 1 means there should be no B p flux through the upper, lower, left, and right faces of the cube. This means that if there is B p flux through the front face (P p0 = 0), it must come out from the back face. Equivalently, if there is a loop-like excitation going through our region of discussion, the region must be at the side of the loop. On the other hand, imposing P corner = 1 means there should be no B p flux through the upper, lower, right, and back faces of the cube. This means that if there is B p flux through the front face (P p0 = 0), it must come out from the left face. Equivalently, if there is a loop-like excitation going through our region of discussion, the region must be at the corner of the loop. Here, we emphasize again that the excitation of B p0 can be either at the side of the string excitation which is specified by the protector P side = 1, or at the corner of the string excitation which is specified by the protector P corner = 1.
With certain constraint of the local Hilbert space, we find the flux operators can be mapped to a subalgebra of the QDM algebra in Eq. (1) with gauge group G local . From all the examples we demonstrated, the subalgebra is enough to distinguish all of the conjugacy classes of G local [61]. This allows us to identify the magnetic fluxes of the fracton model constructed by the mixed symmetry G with the conjugacy classes of G local .
It may happen that certain constraint of the local Hilbert space does the above, but some other constraint does not, and corresponds (at least for the fluxes) to the QDM of another gauge group instead. In the example we study (1 → Z sub 2 → G → Z glo 2 → 1), this situation leads to the fusion of two flexible loops into four lineons, a phenomenon first observed in Ref. [52]. In general, we expect that such constraintdependent QDMs in the Abelian cases [62] lead to the fusion of excitations of one geometry into that of another geometry. In particular, the parts of the original geometry satisfying each constraint (e.g. sides, corners) fuse as the corresponding QDM.
Here we demonstrate the classifications of charges and fluxes from the general gauging procedure and the algebraic structure of QDMs with four examples in detail.
The ungauged system is a three dimensional (3D) cubic lattice with a qutrit (labeled as 0) and a qubit (labeled as 2) on each site, corresponding to the quantum clock matter and the twist defect charge. The Hamiltonian is where the minimal couplings are Here X and Z are the (generalized) Pauli operators, [ψ] := |ψ ψ| is denoted as a projector, and the operators on the first and the second entries on the site act on the qutrit and the qubit, respectively. The other two directions are represented as These graphical notations of operators on the lattice define the minimal couplings and are included in the plaquettes in H n . All graphical notations in this article follow these orientations.
P and g (2) , where g (0) P acts on the matter by mapping |i to |i + 1 mod 3 on the qutrit on each site of a shifted coordinate plane P ; g (2) acts on the matter by flipping the qubit and mapping |i to | − i mod 3 on the qutrit on each site of the entire system. That is The procedure of gauging works as follows: 1. We put one qubit on each link (corresponds to c 2 ) and one qutrit on each plaquette (corresponds to c 0 ). The generalized Pauli operators for the qutrits are denoted as X 0 := X, X 1 := X † , Z 0 := Z, and Z 1 := Z † for notational convenience. I.e., the generalized Pauli operators have an additional subscript for the gauge qutrits on each plaquettes.
2. Based on the general gauging procedure we described in Sec. II.2, the gauge transformations are where v in A v,g indicates the site in the center. The gauge transformations for other g can be obtained by A v,g2 from the above two cases. One should be notified that i in X ( †) i is integer modulo two. Note that these gauge transformations generate a local Z 3 ⋊ Z 2 ∼ = S 3 symmetry for the gauged system: 3. There are two ways to combine the c's to produce the identity. The first type is the product of four c 2 on the links: The corresponding gauge field term is the plaquette operator of the link qubit.
The second type is the tube-like product and sum of the c i 's on the plaquettes (define c 1 = c † 0 for notational convenience.) where [a] c := 1 2 (c 2 + (−1) a ). Here we denote B t,1 = B † t,0 . and the corresponding field terms are where P p is the projector onto the subspace that B p = 1 and p runs over the six faces of the smallest cube containing the tube t. One should be notified that the projectors are essential. B t,0 +B t,1 commutes with all A v,g (so that the Hamiltonian is gauge invariant), but would not do so if we exclude the projectors P p .

The gauged couplings are
Note that in c 0 (τ ), the gauge DOF from c 2 along the boundary of the plaquette is also included. This is in contrast to gauging pure global or pure subsystem symmetry, in which c(τ ) is simply obtained by combining c with the corresponding τ c . This reflects the fact that when global and subsystem symmetry is mixed together, the branch cut created by the global part may split the coupling of the subsystem part.
The requirement b + c + d + e = 0 in the summation reflects the fact that if there is a Z glo 2 flux in the plaquette (B p = −1), then we cannot use A v,g to eliminate all gauge DOF along the plaquette boundary, giving the value 0 for the gauged coupling.

The electric field terms are
The decoratedX operator ensures the electric field is invariant under the gauge transformation A v,g . Similar to the toric code and the X-Cube code, in this example, the matter DOF can also be eliminate by choosing the unitary gauge. We regard the "physical state" satisfying the constraints as an equivalent class of computational basis state, and choose the representative such that all matter DOF are in the reference state |00 . In this case, the Hamiltonian becomes where A u v,g is the gauge part of A v,g by removing the operation on the center vertex of (13). The decoratedX operator becomes X in the unitary gauge. The Hilbert space contains only the gauge DOF without constraints. The above Hamiltonian H u describes the pure lattice gauge theory of Z sub Note that this model can also be constructed in the way similar to Refs. [44 and 45] by gauging the charge conjugation symmetry of the Z 3 X-Cube code [29]. Also note that the combination of global and subsystem symmetry of this model is similar to the field-theoretic approach in Ref. [53], where the Z sub 3 is the discrete version of U (1) x3 with constraints by the cubicity of the lattice.
In the exactly solvable limit (J 0 = J 2 = g 0 = g 2 = 0), the system is deeply in the deconfined phase and the fundamental excitations are fully mobile particles, fractons, and strings. These excitations can be categorized as electric charges or magnetic fluxes. The electric charges are the excitations from , and A u v,g (2) . The magnetic fluxes are the excitations from B p , B t,0 , and B t,1 . We discuss these excitations as follows.

Electric charge excitations
• [f 0 ]: The non-Abelian fracton, corresponding to the ex- where v is the upper left corner and i = 0, 1. Note that This means that the fracton created by M belongs to the same superselection sector, so is the same species [f 0 ]. Also note that this membrane operator only work when acting on a state with B p = +1 for all plaquette p on and near (at most one lattice spacing) the membrane.
To show that it is non-Abelian, we consider the following two states where L is large. The excitation patterns of the two states are both v,0 . On the other hand, they cannot be fused to the vacuum for |ψ 2 . This means that [f 0 ] has non-trivial fusion rules, and therefore is non-Abelian. If we continue this pattern to apply n M v,i 's to the ground state, it can be shown that there are asymptotically 2 n fusion channels, which implies the quantum dimension of [f 0 ] is 2.
• φ: The Z glo 2 charge, corresponding to the excitation of A u v,g (2) in (22). It is an Abelian quasiparticle that can move freely in 3D, created at the end point of the string operator We emphasize that these electric charge excitations are local with respect to the vertex v, and can be specified from the local operators which form a representation of S 3 on the Hilbert space. Hence we can identify electric charges of the G = Z sub 3 ⋊ Z glo 2 fracton model with the irreducible representations of G local ∼ = S 3 (left panel of Table I). This identification is the same as the QDM with S 3 symmetry [11].  (22). They are created at the endpoints at the string operator (the direction of the string is d) where v is the upper left vertex and i = 0, 1. Similar to so i = 0, 1 results in the same species. Also, this string operator only works when acting on a state with B p = +1 for all plaquette p on and near (at most one lattice spacing) this string.
To show that it is non-Abelian, we consider the following two states where L is large. The excitation patterns of the two states are both Now, for the state |ψ 1 , the two [e d ]'s on the left of (31) can be fused to the vacuum by applying an additional S L v,0 . On the other hand, they cannot be fused to the vacuum for |ψ 2 . This means that [e d ] has non-trivial fusion rules, and therefore is non-Abelian. If we continue this pattern to apply n S v,i 's to the ground state, it can be shown that there are asymptotically 2 n fusion channels, which implies the quantum dimension of [e d ] is 2.
• σ: The Z glo 2 flux, which is a flexible string-like excitation (referred to the σ string) corresponding to the excitation of B p in (22). It is created along the boundary of the membrane operator . . . . . . wherẽ The projectors on this decorated membrane operator ensures no additional excitations on the membrane. In the same spirit as the QDM, the magnetic fluxes would be the conjugacy classes of G local ∼ = S 3 . With the geometry defined in Eq. (7) and the constraints defined in Eq. (8), the flux operators can be mapped to a subalgebra of the QDM algebra for either P side = 1 or P corner = 1, where B g is the QDM flux operator in Eq. (1). We can also identify the star operators A u v0,g with the QDM star operator A g Eq. (1). The map of A's and B's together forms an injec-tive algebra homomorphism into the QDM algebra. In particular, they satisfy the relations (note that all the B's commute) In this way, the star and the flux operator of this fracton models is identified with a subalgebra of that of the corresponding QDM.
Note that this subalgebra is enough to distinguish all of the conjugacy classes of S 3 (more rigorously, if a state is the eigenstate of B g0 = 1, B g = 0, g = g 0 for some g 0 ∈ S 3 , then we can determine the conjugacy class of g 0 only by the eigenvalues of the image of B p0 and B t0,i .) This allows us to identify the fluxes of this fracton model with the conjugacy classes of S 3 . The corresponding fluxes are listed in the right panel of Table I.
The ungauged system is a 3D cubic lattice with three qubits (labeled as 0, 1, 2) on each site. 0 and 1 labels are corresponding to the Ising spin of the first layer and the second layer. 2 label is the twist defect charge in Refs. [44 and 45]. For convenience, we combine the first and the second qubits together and denote the operators by X 0 := XI, X 1 := IX, Z 0 := ZI, and Z 1 := IZ. With these definitions, the Hamiltonian and the gauging procedure we describe in this example will be similar as the first example, G = Z sub 3 ⋊ Z glo 2 . The Hamiltonian is where the minimal couplings are Sign representation Abelian fully mobile particle σ −1 0 {g (2) , g (2) g (0) , g (2) 2D representation Non-Abelian fracton {g (0) , (g (0) ) 2 } Non-Abelian lineon Here X and Z are the Pauli operators, [ψ] := |ψ ψ| is denoted as a projector, and the operators on the first and second entries on the site act on the Ising spin on the first and second layers, and the twist defect charge, respectively.
P and g (2) , where g (i) P acts on the matter by flipping (apply X) the ith qubit (i = 0, 1) on each site of a shifted coordinate plane P ; g (2) acts on the matter by flipping the 3rd qubit and swapping the first two qubits on each site of the entire system.
The procedure of gauging works as follows: 1. We put one qubit on each link (corresponds to c 2 ) and two qubits on each plaquette (corresponds to c 0 and c 1 , and labeled by 0 and 1). These qubits are the gauge fields.
2. Based on the general gauging procedure we described in Sec. II.2 [details are demonstrated in Appendix A 1], the gauge transformations are where v in A v,g indicates the site in the center. One should be notified that i in X i is integer modulo two.
On the other hand, the X and [a] on the link refer to the operators on the gauge field with respect to the twist defect charge. The gauge transformations for other g can be obtained by A v,g1g2 = A v,g1 A v,g2 from the above two cases. Note that they generate a local (Z 2 × Z 2 ) ⋊ Z 2 ∼ = D 4 symmetry for the gauged system: 3. There are two ways to combine the c's to produce the identity. The first type is the product of four c 2 on the links, which has the identical structure as the first example in Eq. (15). The corresponding gauge field term B p is the plaquette operator of the link qubit, which also has the identical structure as the first example in Eq. (16). The second type is the tube-like product and sum of the c i 's on the plaquettes where [a] c := 1 2 (c 2 + (−1) a ). The corresponding field term is where P p is the projector onto the subspace that B p = 1 and p runs over the six faces of the smallest cube containing the tube t. One should be notified that the projectors are essential. B t,0 +B t,1 commutes with all A v,g (so that the Hamiltonian is gauge invariant), but would not do so if we exclude the projectors P p . , Note that in c 0 (τ ) and c 1 (τ ), the gauge DOF from c 2 along the boundary of the plaquette is also included. This is in contrast to gauging pure global or pure subsystem symmetry, in which c(τ ) is simply obtained by combining c with the corresponding τ c . This reflects the fact that when global and subsystem symmetry is mixed together, the branch cut created by the global part may split the coupling of the subsystem part. In this case, we need to monitor the switching of the layer index created by the branch cut along the plaquette. The requirement b + c + d + e = 0 in the summation reflects the fact that if there is a Z glo 2 flux in the plaquette (B p = −1), we will assign the value 0 for the gauged coupling, which gives the correct gauge invariant c i (τ ).

The electric field terms are
Again, here [ab] := |ab ab| is the projector on the first and second qubits. The decoratedX operator ensures the electric field is invariant under the gauge transformation A v,g . Similar to the toric code and the X-Cube code, we can choose the unitary gauge to eliminate the matter DOF in this example. We regard the "physical state" satisfying the constraints as an equivalent class of computational basis state, and choose the representative such that all matter DOF are in the reference state |000 . In this case, the Hamiltonian becomes and the Hilbert space contains only the gauge DOF without constraints. The above Hamiltonian H u describes the pure lattice gauge theory of (Z sub 2 × Z sub 2 ) ⋊ Z glo 2 . In the exactly solvable limit (J 0 = J 1 = J 2 = g 0 = g 2 = 0), the system is deeply in the deconfined phase and the fundamental excitations are fully mobile particles, fractons, and strings. Similar to the previous case, the excitations can be categorized as electric charges (excitations of A u v,g (2) ) or magnetic fluxes (excitations of B t,0 , B t,1 , B p ). We discuss these excitations as follows.

Electric charge excitations
• [f 0 ]: The non-Abelian fracton, corresponding to the excitation of the first term in (46). It is created at the four corners of the membrane operator, which has the identical structure as the first example in Eq. (23). Note that where A u v,g (2) is the summand of the third term of (46) centered at v. This means that the fracton created by M belongs to the same superselection sector, so is the same species [f 0 ]. Also note that this membrane operator only work when acting on a state with B p = +1 for all plaquette p on and near (at most one lattice spacing) the membrane.
To show that it is non-Abelian, we consider the following two states where L is large. The excitation patterns of the two states are the same as the first example in Eq. v,0 . On the other hand, they cannot be fused to the vacuum for |ψ 2 . This means that [f 0 ] has non-trivial fusion rules, and therefore is non-Abelian. If we continue this pattern to apply n M v,i 's to the ground state, it can be shown that there are asymptotically 2 n fusion channels, which implies the quantum dimension of [f 0 ] is 2. • φ: The Z glo 2 charge, corresponding to the third term in Eq. (46). It is an Abelian quasiparticle that can move freely in 3D, created at the end point of the string operator, which has the identical structure as the first example in Eq. (27).
These electric charge excitations are local with respect to the vertex v, and can be specified from the local operators which form a representation of D 4 on the Hilbert space. Hence we can identify electric charges of the G = (Z sub 2 × Z sub 2 )⋊Z glo 2 fracton model with the irreducible representations of G local ∼ = D 4 (left panel of Table II). This identification is the same as the QDM with D 4 symmetry [11].

The magnetic flux excitations
• [e d ], where d = x, y, z: The non-Abelian lineons are constrained to move in the d direction, corresponding to the forth and fifth terms in Eq. (46). They are created at the endpoints at the string operator (the direction of the string is d). The string operator has the identical structure as the first example in Eq. (28). Similar to so i = 0, 1 results in the same species. Also, this string operator only works when acting on a state with B p = +1 for all plaquette p on and near (at most one lattice spacing) this string.
To show that it is non-Abelian, we follow the same discussion in the first example. We consider these two states where L is large. The excitation patterns of the two states have the identical structure as the first example in Eq. (31). Now, for the state |ψ 1 , the two [e d ]'s on the left of Eq. (31) can be fused to the vacuum by applying an additional S L v,0 . On the other hand, they cannot be fused to the vacuum for |ψ 2 . This means that [e d ] has non-trivial fusion rules, and therefore is non-Abelian. If we continue this pattern to apply n S v,i 's to the ground state, it can be shown that there are asymptotically 2 n fusion channels, which implies the quantum dimension of [e d ] is 2.
• σ: The Z glo 2 flux, which is a flexible string-like excitation (referred to the σ string) corresponding to the fourth term in (46). It is created along the boundary of the membrane operator {aij }, {bij } =0,1X (11) ⑧ ⑧ ⑧ ⑧X (12) wherẽ The projectors on this decorated membrane operator ensures no additional excitations on the membrane. In the same spirit as the QDM, the magnetic fluxes would be the conjugacy classes of G local ∼ = D 4 . With the geometry defined in Eq. (7) and the constraints defined in Eq. (8), the flux operators can be mapped to a subalgebra of the QDM algebra for either P side = 1 or P corner = 1, where B g is the QDM flux operator in Eq. (1). We can also identify the star operators A u v0,g with the QDM star operator A g Eq. (1). The map of A's and B's together forms an injective algebra homomorphism into the QDM algebra. In partic- (1, 1, −1) Abelian fully mobile particle σ −1 0, 0, 1  . (a, b, c) indicates the representation that g (0) → a, g (1) → b, g (2) → c. ular, they satisfy the relations (note that all the B's commute) (2) , g (0) g (2) , g (1) g (2) , g (0) g (1) g (2) , i = 0, 1.
In this way, the star and the flux operator of this fracton models is identified with a subalgebra of that of the corresponding QDM. Unlike the S 3 case, this subalgebra is not enough to distinguish all of the conjugacy classes of D 4 , since {g (2) , g (0) g (1) g (2) } and {g (1) g (2) , g (0) g (2) } both correspond to B p0 = −1, B t0,0 = B t0,1 = 0. To avoid this situation, we define another flux operator which does not exist in the original Hamiltonian This flux operator is interpreted as a tube operator wrapping around t twice as B t,0 B t,1 when B p0 = −1 (p 0 is the front plaquette). Then we have the additional correspondence and additional relations After including B t0,2 , this subalgebra is enough to distinguish all of the conjugacy classes of D 4 . To be more precise, if a state is the eigenstate of B g0 = 1, B g = 0, g = g 0 for some g 0 ∈ D 4 , then we can determine the conjugacy class of g 0 only by the eigenvalues of the image of B p0 and B t0,i . This allows us to identify the fluxes of this fracton model with the conjugacy classes of D 4 . The corresponding fluxes is in the right panel of Table II.
The ungauged system is a 3D cubic lattice with a d = 4 qudit on each site. The Hamiltonian is where the minimal couplings are where f (1) = f (i) = 1, f (−1) = f (−i) = −1. The construction of this function f , along with another choice of the coupling terms which are more similar to the previous cases, is discussed in Appendix B. Althoughc 0 is not a coupling in H n , it can be obtained from c 0 and c 1 : where the square root takes the eigenvalue −1 of the operator inside it to +i. G is generated by g (0) P and g (1) , where Note that g (0) P 's generate a normal subgroup Z sub 2 of G, but G is not a semidirect product of this Z sub 2 with G/Z sub 2 ∼ = Z glo 2 . Nevertheless, G can still be presented as the group extension In the notation of Ref. [52], this symmetry is called (Z 4 , Z 2 ).
The procedure of gauging works as follows: 1. We put one qubit on each link (corresponds to c 1 ) and one qubit on each plaquette (corresponds to c 0 ). Although we do not viewc 0 as a minimal coupling term in H o and therefore do not put the corresponding gauge qudit, we can still define the "dressed Z" operator associated withc 0 on the plaquettes, which acts on the gauge qubits on both the plaquettes and the links: In addition, we define the "dressed X" operator associated to the link They satisfy the commutation relations and all other combinations commute.
2. Based on the general gauging procedure we described in Sec. II.2, the gauge transformations are Note that these gauge transformations generate a local Z 4 symmetry for the gauged system: We will denote g (0) n simply by n for n = 0, 1, 2, 3, as the usual notation for Z 4 .
3. There are two ways to combine the c's to produce the identity. The first type is the product of four c 1 on the links, which has the identical structure as the first example in Eq. (15) [replacing c 2 by c 1 ]. The corresponding gauge field term B p is the plaquette operator of the link qubit, which has the identical structure as the first example in Eq. (16). The second type is the tube-like product ofc 0 The corresponding field term is Note that, unlike the non-Abelian cases, the projectors P p are not necessary. This is because that R t ({Z}) commutes with all A v,g .

The gauged couplings are
Again, in this case the projector P p is not necessary.

The electric field terms are
Similar to the toric code and the X-Cube code, we can choose the unitary gauge to eliminate the matter DOF in this example. We regard the "physical state" satisfying the constraints as an equivalent class of computational basis state, and choose the representative such that all matter DOF are in the reference state |0 . In this case, the Hamiltonian becomes and the Hilbert space contains only the gauge DOF without constraints. The above Hamiltonian H u describes the pure lattice gauge theory of G. Note that if we replaceX ( †) and Z ( †) by the usual Pauli operators X and Z on the corresponding links and plaquettes, the resulting model will be just the tensor product of an X-cube code and a 3D toric code. In the exactly solvable limit (J 0 = J 1 = g 0 = g 1 = 0), the system is deeply in the deconfined phase and the fundamental excitations are fully mobile particles, fractons, and strings. These excitations are electric charges which corre- , and the magnetic fluxes which are the excitations of B p , B t , B † t . We discuss these excitations as follows.
with its antiparticle e 3 created at the upper right and the lower left corner.
• e 2 : The combination of two e's, which is an Abelian fully mobile particle. It corresponds to the excitation of the first line of Eq. 72, with the eigenvalues of the four terms being 1, −1, 1 and −1, respectively. Besides the e 2 excitation can be created at the edge of the square of the membrane for e, it can also be created at the endpoints of the string operator, which has the identical structure as the first example Eq. (27). Indeed, the square of the membrane for e is just two copies of this string, one at the top and one at the bottom of the membrane.
The excitations of A u v,g (1) , A u v,g  Table III). This identification is the same as the QDM with Z 4 symmetry [11].

The magnetic flux excitations
• m: The flexible string-like excitation corresponding to the excitation B p = −1. It is created at the boundary of the membrane operator Note that this membrane also excites the B t operators. If we look at the R t ({Z}) (without the projectors), the excitation has eigenvalues ±i at the corners of the membrane, but still have 1 otherwise. The m 3 string can be created by the similar membrane operator with replac-ingX →X † .
• m 2 d , where d = x, y, z: The Abelian lineon constrained to move in the d direction, corresponding to the excitation having B t = −1. It is created at the endpoints at the string operator (the direction of the string is d) It is also created at the corners of the square of the membrane for m.
In the same spirit as the QDM, the magnetic fluxes would be the conjugacy classes of G local ∼ = Z 4 , that is, elements of Z 4 . With the geometry defined in Eq. (7) and the constraints defined in Eq. (8), the flux operators can be mapped to a subalgebra of the QDM algebra for only P corner = 1, where B g is the QDM flux operator in Eq. (1). We can also identify the star operators A u v0,g with the QDM star operator A g Eq. (1). The map of A's and B's together forms an injective algebra homomorphism into the QDM algebra. In particular, they satisfy the relations (note that all the B's commute) In this way, the star and the flux operator of this fracton models is identified with a subalgebra of that of the corresponding QDM. Unlike the S 3 case, this subalgebra is not enough to distinguish all of the conjugacy classes of Z 4 , since {1} and {3} both correspond to B p0 = −1, B t0 = 0. To avoid this situation, we note that in this Abelian case, the zero-flux projectors in the definition of B t0 is not necessary. Hence we can define the flux operator without the projectors. (this does not mean that we need to modified the Hamiltonian) Then we have the correspondence and additional relations After including B ′ t0 , the corresponding subalgebra is enough to distinguish all of the elements of Z 4 (more rigorously, if a state is the eigenstate of B g0 = 1, B g = 0, g = g 0 for some g 0 ∈ Z 4 , then we can determine g 0 only by the eigenvalues of the image of B p0 , B t0 and B ′ t0 . In this case, B ′ t0 alone is enough to distinguish them.) This allows us to identify the fluxes of this fracton model with the conjugacy classes of Z 4 in the situation that P corner = 1. The corresponding fluxes is in the right panel of Table III.
However, if we choose P side = 1 instead, then we have B ′ 2 t0 = 1 on the fracton side instead of B ′ 2 t0 = B p0 , which implies that the previous mapping into Z 4 QDM is no longer a homomorphism. Curiously, it is isomorphic to the Z 2 × Z 2 QDM algebra by the mapping Here the first Z 2 is associated to m 2 d , and the second Z 2 is associated to m. This explains why when we apply the square of the membrane for the string excitation m, we create four lineons at the corners but nothing at the edge of the membrane. The fusion rules of the fluxes behave like the Z 4 QDM (m×m = m 2 ) only at the corners, while at the edge of a string excitation, they behave like a Z 2 × Z 2 QDM (m × m = 0) instead. This phenomenon was first observed in Ref. [52], and we expect that it is general for gauged Abelian models with constraint-dependent QDM correspondence (see discussion before the first example).
The ungauged system is a 3D cubic lattice with three qubits on each site. We use the following correspondence of the qubits with the elements of Q 8 : For notational convenience, we combine the last two qubits together and define X 1 = XI, X 2 = IX, Z 1 = ZI, Z 2 = IZ, and II = I. The Hamiltonian is where a := a (1) a (2) with a (1) , a (2) = 0, 1, and Z a := Z a (1) 1 Z a (1) +a (2) 2 on the second entry acts on the second and the third qubits. The projector [a] := [a (1) a (2) ] = |a (1) a (2) a (1) a (2) | is defined on the second entry acting on the second and the third qubits. G is generated by g (−) where (· · · ) acts on the second and third qubits. These generators are the subsystem multiplication of −1 ∈ Q 8 and the global multiplication of i, j, k ∈ Q 8 , respectively. Note that P 's generate a normal subgroup Z sub 2 of G, but G is not a semidirect product of this Z sub 2 with G/Z sub 2 ∼ = K glo 4 . Nevertheless, G can still be presented as the group extension the links: The corresponding gauge field term is the plaquette operator of the link qubit.
The second type is the tube-like product The corresponding field term is six faces P p,1 P p,2 = a,b=00,01,10,11 4. The gauged coupling terms are c 0 (τ ) = a,b,c,d,e =00,01,10,11 b+c+d+e=00 5. the electric field terms are wherẽ [γb] The decoratedX (i) andX (j) operators ensure the electric field is invariant under the gauge transformation A v,g . Similar to the toric code and the X-Cube code, we can choose the unitary gauge to eliminate the matter DOF in this example. We regard the "physical state" satisfying the constraints as an equivalent class of computational basis state, and choose the representative such that all matter DOF are in the reference state |000 .
In this case, the Hamiltonian becomes where the dressed operators in the unitary gauge are [63] u . The Hilbert space contains only the gauge DOF without constraints. The above Hamiltonian H u describes the pure lattice gauge theory of G. In the exactly solvable limit (J 0 = J 1 = J 2 = g 0 = g 1 = g 2 ), the system is deeply in the deconfined phase and the fundamental excitations are fully mobile particles, fractons, and strings. These excitations are electric charges which are the excitations of and magnetic fluxes which are the excitations of B t , B p,1 , B p,2 . We discuss these excitations in the following.

Electric change excitations
• φ (q) , q = i, j, k: An Abelian quasiparticle that can move freely in 3D, corresponding to A v,g (q ′ ) = −1, q ′ = q. For q = i, j, k, it is created at the end point of the string operator respectively. • where v is the upper left corner Z = a=00,01,10,11 This means that the fracton created by M (m,n) v,b0 , b 0 = 00, 01, 10, 11 belongs to the same superselection sector, so is the same species [f 0 ]. Also note that this membrane operator only work when acting on a state with B p,i = +1 for all plaquette p on and near (at most one lattice spacing) the membrane. To show that it is non-Abelian, we consider the following two states where L is large. The excitation patterns of the two states are the same as the first example in Eq. (26). Now, for the state |ψ 1 , the four v,00 . On the other hand, they cannot be fused to the vacuum for |ψ 2 . This means that [f 0 ] has non-trivial fusion rules, and therefore is non-Abelian. If we continue this pattern to apply n M v,i 's to the ground state, it can be shown that there are asymptotically 2 n fusion channels, which implies the quantum dimension of [f 0 ] is 2.  Table IV). This identification is the same as the quantum double model (QDM) with Q 8 symmetry [11].
u are defined before and In the same spirit as the QDM, the magnetic fluxes would be the conjugacy classes of G local ∼ = Q 8 . With the geometry defined in Eq. (7) and the constraints defined in Eq. (8), the flux operators can be mapped to a subalgebra of the QDM algebra for either P side = 1 or P corner = 1, where B g is the QDM flux operator in Eq. (1). We can also identify the star operators A u v0,g with the QDM star operator A g Eq. (1). The map of A's and B's together forms an injective algebra homomorphism into the QDM algebra. In particular, they satisfy the relations (note that all the B's commute) In this way, the star and the flux operator of this fracton models is identified with a subalgebra of that of the corresponding QDM.
Note that this subalgebra is enough to distinguish all of the conjugacy classes of Q 8 (more rigorously, if a state is the eigenstate of B g0 = 1, B g = 0, g = g 0 for some g 0 ∈ Q 8 , then we can determine the conjugacy class of g 0 only by the eigenvalues of the image of B p0 and B t0,i .) This allows us to identify the fluxes of this fracton model with the conjugacy classes of Q 8 .The corresponding fluxes are listed in the right panel of Table IV.

IV. CONCLUSION
In this paper, we demonstrate a systematical gauging procedure on a lattice with pure matter fields. This general procedure works not only for a semidirect product of the global symmetry and the subsystem symmetry, but also for a nontrivial extension of them. We give four examples of gauging The former two cases and the last one produce the non-Abelian fracton orders. By using a one-step gauging, we give a transparent identification of electric charges with the irreducible representations of G local , which include the non-Abelian fracton orders. Furthermore, to compare the magnetic excitations with different geometry (tubes and plaquettes), we set a specific constraint on the local Hilbert space. We observe that the magnetic excitations satisfy the subalgebra of the QDMs, which allows us to identify the magnetic excitations as the conjugacy classes of G local . Our gauging procedure is very general, and can be easily extended to more exotic symmetries, and can produce different types of non-Abelian version of fractons.
Before we close the discussion, we would like to point out some future directions.
1. If one applies the symmetry of "(1D subsystem) ⋊ (2D subsystem) ⋊ global" or "fractal ⋊ global" form, we expect the gauged Hamiltonian in the exactly-solvable limit can produce new types of non-Abelian fractons.
2. The fully gauged Hamiltonian contains matter and gauge fields, which may contain partial confineddeconfined or Higgsed-deconfined transition for some of the excitations. One can also consider other types of matter fields such as the Majorana fermion.
3. The quotient superselection sector (QSS) [64] is a method to identify fracton species. We would like to understand the deep relation between a non-Abelian generalization of QSS and the corresponding electric charges with the irreps of G local .
4. Although we associate the magnetic fluxes with the conjugacy class of G local in the constrained local Hilbert space from the subalgebra of the QDMs, some magnetic fluxes (e.g., B t0,2 in the second example) cannot be obtained from our gauging procedure. A direct identification of the fluxes with structures of G and G local from the general gauging process is still desired.
5. The ungauging map [65] is the map from the zero-flux subspace of a pure lattice gauge system back to the matter Hilbert space with the symmetry. We expect that it can be generalized to non-Abelian symmetries.
Note added. While preparing the update of this paper, we recently learnt of a related work by Tantivasadakarn, Ji, and Vijay [66], which they construct a solvable lattice model labelled as (Q 8 , Z 2 ) in their terminology, which is our fourth example.
Here |m s(v) are the state of matter fields, and C v being a subset of the couplings associated with the site v (we will state the requirements for C v later). We now demonstrate a systematical way to find |τ ′ c after the gauge transformation, which defines A v,g acting on |τ c .
1. We construct a virtual state of matter fields |m = s =v |m s ⊗|m v associated with |τ c by the coupling c. This virtual state of matter fields |m is chosen to satisfy That means for a state of gauge field |τ c , we can find a corresponding virtual matter state |m v that satisfies Eq. (A2).
2. |τ ′ c associates with another virtual state |m ′ . Since |τ ′ c is the compensation of the changes of A v,g acting on the matter field, the virtual state |m ′ is obtained by applying g −1 to the state |m on v, |m ′ = s =v |m s ⊗ g −1 |m v .
3. The eigenvalue of c for the virtual state |m ′ c , c|m ′ c = τ ′ c |m ′ c then specifies the state |τ ′ c . This step determines the gauge transformation A v,g on |τ c . It is required that |m and |m ′ are eigenstates of c with the same eigenvalue, ∀c / ∈ C v (otherwise one should choose a larger C v ).
The above definition automatically implies A v,g1 A v,g2 = A v,g1g2 , and C v is chosen to be the minimal set such that A v,g is well-defined for all g and that [A v1,g1 , A v2,g2 ] = 0 when v 1 = v 2 . This means that G is promoted to a local symmetry of the gauged system. In the cases of gauging pure global or pure subsystem symmetries, we simply have C v = {c | v ∈ supp c}. However, when global and subsystem symmetries are mixed together, the choice of C v becomes nontrivial. See the example below.
1. Specific example for obtaining Av,g for (Z sub The reference state is |000 . At each vertex v [the center vertex in Eq. (A3)], the support C v of the gauge transformation is chosen to contain twelve c 0 and c 1 couplings on the plaquettes, and nine c 2 couplings on the links, (A3) Here we demonstrate the calculation of A v,g (0) P , where P is any of the three shifted coordinate planes containing v. Since g (0) P acts on v by the operator X 0 , A v,g (0) P also acts on v by X 0 . Also, by definition, A v,g (0) P acts on other sites trivially. Now we demonstrate how to calculate its action on the gauge DOF. Suppose we want to calculate A v,g (0) P acting on a computation basis state in which the part of the state around v is where v is the center vertex. Note that we only include the states of the τ c 's corresponds to couplings c ∈ C v in Eq. (A3). We will discuss this choice of C v below. Now the steps works as follows: 1. Consider the virtual situation that the matter is in the state (only the state relevant to C v is shown) This satisfies the requirements, since the state of v is the reference state |000 and one can check that the eigenvalues of the minimal couplings are in correspondence with |τ c . For example, |m is the (+1)-eigenstate of the c 0 coupling on the upper left blue plaquette and is the (−1)-eigenstate of the c 1 on that plaquette, which corresponds to the state |01 of the corresponding τ c0 and τ c1 shown in Eq. (A4) (Recall that we label the computational basis of τ c as |0 and |1 corresponding to the +1 and −1 eigenvalues of c). Of course, this is not the only virtual matter state that satisfies the requirement, but any of such states should give the same result. See the discussion on C v below.
2. The action of g on v is X 0 , which maps v on the virtual matter to |100 . Also, other sites of the virtual matter are unchanged. So the virtual matter state is changed to 3. Now we calculate the eigenvalues of the minimal couplings for |m ′ , and correspond them back to the state of τ c . The result is which is defined as the result of A v,g (0) P acting on the state in Eq. (A4). Note that for a minimal coupling c / ∈ C v [corresponding to a state not shown in Eq. (A4)], the eigenvalues for |m and |m ′ are the same. This is important when we discuss the choice of C v below.
The above steps give the result of A v,g (0) P acting on a specific computation basis state. To obtain the full expression of A v,g (0) P , one needs to consider its action on all possible states for gauge DOF (that is, change the state of Eq. (A4) to any possible combinations of 0 and 1. Then one should obtain the first line of Eq. (39). Now we discuss our choice of C v . Firstly, one can show that if we choose a different |m that also corresponds to (A4), then we still obtain the same Eq. (A7). Secondly, in the last step above, the eigenvalues of c / ∈ C v are not changed. This is true in general for this choice of C v , and therefore A v,g is well-defined. Also one can check [see Eq. (39)] that A v,g on different sites commute. Now we discuss why this choice is minimal. That is, if we delete some couplings from C v , then the above procedure will fail. Suppose we delete the c 2 in the upper left blue horizontal link in Eq. (A3). That is, we exclude the state on the horizontal blue link in Eq. (A4). Then Eq. (A5) still satisfies the requirements. In addition to that, if we modify the upper left corner from |101 to |100 in Eq. (A5), it also satisfies the requirements. However, if we use the latter choice of |m , than the upper left plaquette of the resulting Eq. (A7) will be |11 instead of |00 . That is, different choices of virtual matter lead to different result, so this C v does not work.
Next, suppose we delete the c 0 in the upper left blue horizontal link in Eq. (A3). Note that in the last step, the eigenvalue of this c 0 still changes in respond to the change of the virtual state of v. Then since this c 0 is not in C v , it violates the requirement that all eigenvalues of c / ∈ C v is not changed. Other deletion of elements in C v violates the requirements in similar ways. Note that some deletion still works when calculating A v,g (0) P but not when calculating A v,g (2) (we require that a single C v works for any g around a vertex.) Therefore, the chosen C v is minimal.
Note that this is not the only choice of C v , but other choice are essentially the same. If we require that A v,g is symmetric