Abelian combinatorial gauge symmetry

Combinatorial gauge symmetry is a principle that allows us to construct lattice gauge theories with two key and distinguishing properties: a) only one- and two-body interactions are needed; and b) the symmetry is exact rather than emergent in an effective or perturbative limit. The ground state exhibits topological order for a range of parameters. This paper is a generalization of the construction to any finite Abelian group. In addition to the general mathematical construction, we present a physical implementation in superconducting wire arrays, which offers a route to the experimental realization of lattice gauge theories with static Hamiltonians.


Introduction
Physically realizing topological ordered [1] gapped quantum spin liquids in the laboratory remains a problem at the frontier of our knowledge in both quantum materials and quantum information sciences.Control over such states would constitute a major step forward in building stable quantum memories, material simulation, and ultimately scalable quantum computation.Over the past decades, major theoretical advances have been made by formulating models within the framework of lattice gauge theories [2], for example the well-known 2 toric code [3].Recently, 2 quantum spin liquids have been demonstrated by preparing a state like the quantum ground state of the static system and its elementary excitations as non-equilibrium states in quantum simulators [4,5].However, as these states are prepared dynamically, they do not exist in the absence of a an external drive or modulation, and thus they do not possess the topological stability of the ground states of the static Hamiltonians.
A set of proposals utilizing the "combinatorial gauge symmetry" (CGS) construction plots a potential alternative path to overcoming these obstacles [22,23].CGS is a framework for constructing Hamiltonians out of realistic two-body interactions where the gauge symmetry is exact and non-perturbative for the full range of parameters.As such there is hope that a suitable range of parameters can be found where the topological phase is sufficiently gapped to be observable and stable.Earlier work has centered on special cases of small Abelian groups ( 2 and 3 ).By generalizing the framework to all Abelian groups we further enlarge the scope of possible implementations and parameters where one can look for such states.The mathematical construction developed in this paper to handle the generalization is interesting in its own right.
From a mathematical point of view, central to the construction is a classification of twobody interactions that can be written as matrices W i j that couple two quantum variables indexed by i and j (e.g., spins or Josephson phases).The allowed sets of interactions W must possess several properties in order to satisfy the requirements of a lattice gauge theory exactly.This paper is a systematic exposition of the general requirements.
As a high level preview, the gauge symmetry arises from automorphisms of the form W = L ⊤ W R where the matrices L and R must be monomials in order to preserve the commutation relations of the underlying quantum variables.Further, because of the geometry of the lattice configuration, R will be required to be diagonal, but L will not.Effectively, the problem of classifying the allowed W matrices will reduce to a study of how group actions relate to permutations of invariant sets.
The outline of the paper is as follows.First we introduce in Sec. 2 a simple example of a model with a combinatorial gauge symmetry corresponding to a 2 gauge theory on the triangular lattice to illustrate key elements of our construction.We summarize the main results of this paper in Sec. 3. In Sec. 4 we develop in detail the general construction for any Abelian gauge group, as well as modifications to the models in specific cases, and then provide illustrations for these in Sec. 5. Finally, we discuss experimental implementations in Sec. 6.

A motivating example
In this section we consider a 2 gauge theory on a triangular lattice with only two-body interactions.Before that, it should be made clear what kind of gauge symmetry such a theory possesses.In an ordinary 2 gauge theory, e.g., the classical toric code model, the degrees of freedom are located on the links of the lattice and the Hamiltonian is where the plaquette operator B p is the product of σ z operators acting on the links bounding the plaquette p.In the full toric code model, this Hamiltonian has gauge symmetries generated by star operators, that is, products of σ x operators along a path through the dual lattice that forms a domain wall around a subset of vertices on the lattice.At the level of individual degrees of freedom, such symmetries act by the transformation σ z → σ x σ z σ x = −σ z .Thus, whenever the support of a plaquette operator intersects with that of a gauge transformation, the degrees of freedom involved will be affected by the gauge transformation, but the plaquette operator will be invariant overall.(See Fig. 1(a).)Since in the Abelian case such a symmetry operation is always realized by a phase factor, we may represent its effect on a plaquette term by a diagonal matrix multiplying the vector of operators contained in this term.For example, if we collect the z-components of the three spins around a star onto a vector σ z ⊤ ≡ (σ z 1 , σ z 2 , σ z 3 ) ⊤ , then the action of the symmetry R1 that flips the first two spins is where R is the diagonal matrix Under the action of R the product , because the phase factors introduced by the two spin flips multiply to one.
Having understood the action of the gauge symmetry on individual degrees of freedom and on gauge invariant operators, we now construct a two-body Hamiltonian that possesses the same gauge symmetry.To this end, in addition to the gauge spinoperators σ z as before, we introduce matter spin operators µ z associated with each plaquette term in the pure gauge theory.More precisely, in this example we introduce four matter spins µ z a , a = 1, 2, 3, 4, which we also collect into a vector We would like the Hamiltonian to (i) have only one-and two-body interactions (here of the Z Z type) and (ii) possess a 2 gauge symmetry when acting on the gauge spins.Condition (i) is consistent with an interaction term of the form −J 4 a=1 3 i=1 W ai µ z a σ z i .Using the vector notation, this interaction can be written compactly as (µ z ) ⊤ W σ z .As for condition (ii), when the gauge symmetry acts on the σ i 's via σ z → Rσ z , the interaction is transformed into (µ z ) ⊤ W Rσ z , which is equivalent to multiplying the coupling matrix W on the right by R. If there exists a matrix L that we can multiply on the left of W so that L ⊤ W R = W , then L is the action of the gauge symmetry on µ z , so that the transformation (σ z , µ z ) → (Rσ z , Lµ z ) leaves the system invariant.The L and R matrices arise from gauge transformations, which ultimately come from conjugating spins operators by some symmetry operator, so the transformation represented by L should preserve the commutation relations between the spin operators.Therefore, we restrict our considerations only to L's that permute the spins or flip µ z 's.In other words, the matrix L should be a monomial matrix with entries ±1.
In the triangular lattice example, we can construct an interaction that satisfies all of the above conditions with four matter spins.The coefficient matrix is From the lattice geometry, we see that gauge transformations acting on the gauge spins around a triangular plaquette are generated by the following two diagonal matrices, The superscripts denote the sites where the spins are flipped.The corresponding L matrices that leave W invariant under W → L ⊤ W R are labeled by the same superscripts as the corresponding R matrices.(See Fig. 1(c) for an illustration of the effect of R 12 and L 12 on the operators σ z i and µ z a .)Next we take a closer look at the properties of these matrices.
By design W in (4) has the symmetry L ⊤ W R = W .The rows of W are the triples of ±1 that multiply to 1.These are precisely the triples that can be obtained from (+, +, +) by flipping the signs of two entries at a time.More formally, the rows of W form the orbit of (+, +, +) under the action of the group of gauge symmetries.Orbits under a group action are examples of invariant subsets, which means that if we flip the sign of any two columns of W via W → W R, i.e. act on all elements of the invariant subset by a group element, the effect is the same as permuting the rows of W , i.e. the elements of the invariant subset.Therefore, for each R, there will always be a permutation matrix L that "undoes" its effect by permuting the rows back into their original positions, so that W R = (L ⊤ ) −1 W is satisfied.For example, with the pair of transformations R 12 and L 12 , we have and Flipping the signs of all entries in the first two columns of W has the same effect as exchanging the first row of W with the last and the second with the third.This interaction term needs to be supplemented with additional terms, so the Hamiltonian of the entire system takes the form where the subscripts p on the operators µ x,z p and σ x,z p denote the plaquette that the operators are located in, the indices a ∈ p label the matter spins in plaquette p, and the indices i label the gauge spins on the set of links L of the lattice.If a large transverse field −Γ p a∈p µ x a is applied on the matter spins, it will effectively integrate out the matter degrees of freedom, so that the interaction term −J p µ z p ⊤ W σ z p will perturbatively generate a pure gauge theory.The lowest order term in perturbation theory will be of the form p i∈p σ z i .In addition, a star term can also be generated perturbatively through the application of the transverse field on the gauge spins −Γ ′ p i∈p σ x i , thus fully reproducing the pure gauge theory we started with.We stress that the gauge symmetry of the system described by Eq. ( 9) is preserved throughout, not merely emergent in a perturbative limit in which we justify integrating out the matter fields, meaning that a topological phase may exist within a large area of the parameter space of couplings J, h, Γ and Γ ′ .
For this system to be equivalent to a 2 gauge theory whose ground state is a spin-liquid state, it is not sufficient to only have the gauge symmetry.We would like ground state in absence of the σ x term to be an eigenstate of the plaquette operators B p = i∈p σ z i with eigenvalue 1, because this enforces a zero-flux condition that corresponds to the confined phase in the pure gauge theory. 2 However, in our motivating example on the triangular lattice, the B p = 1 state is always degenerate with the B p = −1 state, because a global spin flip of both matter and gauge spins commutes with the Hamiltonian but anticommutes with B p .In order to split this degeneracy, we apply a field −h p a∈p µ z a to the matter spins, with the sign of h chosen so that the B p = 1 sector is preferred.Without this term, the ground state manifold is spanned by the flux-0 state {σ, µ} = {(+, +, +), (+, −, −, −)} and its gauge symmetry partners, along with the flux-1 state {(−, −, −), (−, +, +, +)} and its symmetry partners.Since these two sets of states have distinct matter spin magnetization (indeed, the CGS preserves the magnetization of matter spins because it acts as a permutation on the matter spins), we are able to distinguish the flux states by applying a uniform field on the matter spins.
In the above example we outlined how to construct a combinatorial gauge theory whose ground state has topological order.In the following sections we will generalize these ideas.

Summary of results
The construction presented in the previous section can be repeated in full generality for all finite Abelian gauge groups by considering cyclic gauge groups.This is because any finite Abelian group G can be decomposed into a direct product of cyclic groups i k i .Focusing on the on the cyclic group k , we work with operators σ Z and σ X acting on a k clock variable, which satisfy the relation Schematically, the one-and two-body Hamiltonian for such a combinatorial gauge theory takes the following form: where we have extended the vector notation of the example of the previous section to write the Z Z term in the first line of Eq. (11), by collecting all q gauge spins and m matter spins on each plaquette into vectors (11) is the gauge-invariant interaction term.A gauge transformation is represented by a pair of matrices L p and R p acting on the vectors of operators associated with each plaquette, i.e. σ Z,X p → L p σ Z,X p and µ Z,X p → R p µ Z,X p .The diagonal R p matrices are determined by the behavior of the pure gauge theory under gauge transformations, and the L p matrices are monomial and constructed according to the requirement L ⊤ p W R p = W .This last requirement ensures the gauge invariance of this interaction term.It is possible to construct the coupling matrix W for any gauge group.Moreover, by exploiting internal symmetries of matter degrees of freedom, the number of matter variables, or equivalently, the size of the matrix W can be reduced.A CGS Hamiltonian on a lattice with coordination number q will have at least q parameters that we will be able to tune to adjust the spectrum, so that the correct flux sector becomes the ground state.The transverse field −Γ p a∈p µ X a induces dynamics in the matter degrees of freedom.When this term is strong (Γ ≫ J), the matter degrees of freedom are integrated out, and the effective Hamiltonian to the lowest order in perturbation theory is exactly the plaquette term in the pure gauge theory, since that will be the lowest order term in perturbation theory that respects the gauge symmetry.Next if we apply a weak transverse field −Γ ′ i∈L σ X i , the generators of gauge transformations, that is, the star terms in the pure gauge theory, will be produced.Note that as long as the transverse field on the matter spins is uniform in each plaquette, the gauge symmetry is preserved.Thus the Hamiltonian (11) will always possess an exact gauge symmetry.This means that even though the effective Hamiltonian only becomes exactly equal to the pure gauge theory in the Γ → ∞ limit, we should expect the existence of a topological phase that can be connected adiabatically to the confined phase of the pure gauge theory even at finite values of Γ .Indeed, in the case of 2 gauge theory on a square lattice, it has been shown that such a quantum spin liquid phase exists in this limit.[24] This provides evidence that under the Γ ′ ≪ J ≪ Γ energy hierarchy, the topological phase is robust.Since the gap is finite, as long as symmetry breaking perturbations are small compared to this gap, the splitting between the topologically degenerate ground states, or indeed within any flux sector, is exponentially small in the system size, so the topological phase should be preserved.The presence of the term −h p a∈p µ Z a in Eq. ( 11) is required to break accidental degeneracies when certain gauge-violating transformations leave the interaction term invariant.This is another measure to ensure the energetic separation of different flux sectors.

Combinatorial gauge symmetry
Suppose we want to construct a gauge-invariant two-body interaction whose behavior under gauge transformations mirrors that of a pure gauge theory where the first group of terms represents charges, i.e., generators of gauge transformations defined on a star labeled by s, and the second group enforces the constraint that the ground state is flux-free on each plaquette labeled by p.To describe this behavior algebraically, we need to consider three algebraic constructs.First, the gauge fields can be represented by variables σ that take value in the gauge group G and this group naturally acts on these variables. 3When an element of the gauge group g ∈ G acts on a gauge variable σ, we write it schematically as g • σ.
On the entire lattice, the total Hilbert space is a tensor product of all the local gauge variables i σ i , and is acted on by the direct product of the gauge groups i G i , where the factor G i acts only on the variable σ i .Not all elements of this direct product are valid gauge transformations, but only the ones that are generated by the charge operators A s .We call this subgroup G the group of gauge transformations.An operator is gauge invariant if it is invariant under the action of this group.Nevertheless, in order to test if a local operator O is gauge invariant, we need not test it against all possible gauge transformations, but only the ones that act nontrivially on the support S of O.This means that we are in fact interested in the quotient group G under the relation ∼ that identifies two gauge transformations T = i g i and Now we assume that in the pure gauge theory, the flux operator B p is a product of q operators σ 1 , . . ., σ q .However, since our aim is to construct a two-body interaction term involving a set of auxiliary, or "matter", degrees of freedom µ 1 , . . ., µ m , the resulting operator should not be a product of all relevant operators.The most natural form is a quadratic direct coupling between the gauge and matter variables, i,a µ a W ai σ i .Thus we write the gauge and matter degrees of freedom as vectors of operators σ and µ, so that the action of a gauge transformation T ∈ G p on σ can be represented by a diagonal matrix R, where and its action on the matter variables by a matrix L, i.e. µ → Lµ.Now we define that the term i,a µ a W ai σ i = µ ⊤ W σ has a combinatorial gauge symmetry that reproduces the gauge symmetry of the flux operator B p of the pure gauge theory H pure = −λ B p B p , if for every gauge transformation T in the local group of gauge transformations G p of the flux operator B p = i∈p σ i that acts on σ p via the matrix R, we can find a matrix L acting on a vector of matter variables µ, such that (Lµ) ⊤ W (Rσ) = µ ⊤ W σ, or L ⊤ W R = W .The symmetry L of the matter variables should be a physical symmetry of the corresponding physical degrees of freedom.
The overall invariance of µ ⊤ W σ is equivalent to the invariance of W under the multiplication by R on the right and a corresponding L on the left.Since an R matrix multiplies W on the right, it can be thought of as the same transformation applied to each row of W . Thus if the rows of W form an invariant set under all R-transformations then multiplying W on the right by a matrix R has the same effect as permuting the rows, so we can "undo" this transformation by multiplying W on the left by a permutation matrix L ⊤ .This gives us the invariance and consequently the invariance of µ ⊤ W σ under the transformations µ → Lµ and σ → Rσ.
The number of matter degrees of freedom µ is then determined by the number of elements in this invariant set.Since invariant sets are unions of orbits, such a coupling matrix W can be constructed by choosing one or more rows of initial entries, applying all R matrix representations of transformations in G p to obtain the orbits, and combining the results to form W .We can lay out the construction more concretely using k as the gauge group.σ could be represented by a clock variable . The action of the gauge group is generated by the X operator, associated to a k × k permutation matrix that is, take any generator x of the group G and represent it by the X operator, then for any group element g, which can be written as x c for some integer c, its action σ → g • σ is represented by Z → (X † ) c Z X c .In particular, since X Z X † = ζ k Z, such an action on Z reduces to multiplication by the phase ζ c k = e 2πic/k .Thus, the action of G p that we have represented schematically as a matrix multiplication in fact is a matrix multiplication by a diagonal R matrix whose entries are k-th roots of unity.For the flux term B p = i∈p X i , gauge invariance requires that the entries of R multiply to 1.If B p is supported on q sites, the local group of gauge transformations G p will be a quotient of q−1 k .In many cases, G p will be precisely q−1 k , and when k is prime, G p will always be a direct sum of fewer than q copies of k . 4This includes the case of k toric codes, but also models with much more complicated geometry, for example the Haah's code, which will be discussed in the following section.
Supposing G p is of the form q−1 k , we can construct the coupling matrix W by choosing the entries in the first row and generating the other rows under the action of all R ∈ G p .Suppose we take the first row of W to be (1, 1, . . ., 1), then the entries of W will also be k-th roots of unity, and the rows of W consist precisely of all q-tuples of k-th roots of unity that multiply to 1.There are N = k q−1 such combinations, so the CGS model will require k q−1 matter degrees of 4 Each charge operator will generate a subgroup of the local group of gauge transformation isomorphic to the gauge group.The relation among them is of the form R , where R i is the local gauge transformation generated by the i-th charge operator.When k is prime, this relation can be used to express R m in terms of the other generators, but if not, we may only be able to eliminate a power of R m , resulting in a local group of gauge transformation whose direct summands may include non-trivial quotients of k .freedom.Explicitly, where a n, j are integers such that q j=1 a n, j ≡ 0 (mod k), or equivalently q j=1 ζ a n, j k = 1, for all n.This matrix will satisfy the CGS invariance condition (15).In this case a general gauge transformation acts on the gauge degrees of freedom via a multiplication on the right of W k,q by the following diagonal matrix: where c j are integers satisfying the condition q j=1 ζ c j k = 1 or q j=1 c j ≡ 0 (mod k), so that the product of each row of W k,q is preserved.By construction, this right multiplication on W by R has the effect of permuting its rows, so that an N × N permutation matrix L exists to rearrange the rows back in their initial order, thereby satisfying the condition L ⊤ W R = W .
For a specific example, consider the k = 5, q = 3 case, that is, take 5 as gauge group and a triangular lattice.Starting with (1, 1, 1) as the initial row of W , the remaining rows can be obtained by the action of the local group of gauge transformations 2  5 , generated by ), acting on W by right multiplication, producing a 25-element orbit.(See Appendix B for the full form of this matrix.)In the following discussion we will denote such a combinatorially symmetric coupling matrix by W k,q [r 1 , . . ., r n ], where r 1 , . . ., r n are representatives of distinct orbits forming the set of its rows, so the 5 example above is W 5,3 [(1, 1, 1)].(When context makes it clear, we may also omit some or all of the arguments of W .)

Reducing the number of matter spins
In general, we do not impose anything other than permutation symmetry for the matter variables, so in principle the µ's in the Hamiltonian can be any local operator.However, if the matter variables are realized as spin- 1  2 or n degrees of freedom, they could possess an internal symmetry, which generalizes L's to monomial matrices.In some cases, this allows us to reduce the number of matter variables.To be precise, this reduction is possible whenever G p contains a subgroup consisting of elements of the form R = g I, which we may call "uniform" gauge transformations, forming a subgroup G U p of G p .Such uniform gauge transformations multiplies all entries by the same group element g from the right and permutes rows of the coupling matrix differing by an overall phase factor from the left.We can factor out all such uniform transformations and let them act as an internal symmetry of the matter variables, or equivalently quotienting G p by G U p .This reduces the size of G p and therefore the number of matter variables by a factor of G U p .Starting from a fully constructed CGS model with uniform gauge symmetries, we can partition rows of W into sub-orbits consisting of rows related by a uniform gauge transformation and choose only one representative from each to constructed a reduced W red .For this coupling matrix, when the action of R maps a row of W red into rows that are not present in the original, it will always be compensated for by a non-identity entry in the L matrix.We will see examples of this below.

Accidental symmetries
A subtle issue can arise related to uniform gauge transformations.When the matter variables possess internal symmetries, and when such symmetries coincide with elements of the gauge group, they may combine into spurious symmetries of the Hamiltonian that are not members of the local group of gauge transformations.In the language of the coupling matrix W and L, R matrices, this situation is described by a pair of transformations (L, R) = (g −1 I, g I), when g I itself is not a valid gauge transformation.For example, this could occur when the matter variables are spin- 1  2 variables, the gauge group contains an element that is represented by a phase shift by −1, and the flux operator B p contains an odd number of factors.The motivating example discussed in Section 2 falls into this case.Such a spurious symmetry will collapse flux sectors that should remain distinct, so we need to break it with an additional field −h p a µ z p,a on the matter spins to remove the internal symmetry.(See Appendix C for a demonstration that applying such a field is sufficient to lift the spurious symmetry.) Note that uniform spurious symmetries arise in a situation complementary to the situation discussed above where uniform gauge transformations allow for a reduction of the coupling matrix, so breaking the internal symmetry in the former case will not come at a cost of increasing the number of matter variables.

Product gauge groups
The framework for CGS is applicable to general Abelian gauge groups, but it is sufficient to consider the more restricted case of cyclic groups, because any finite Abelian group can be decomposed into a direct product of cyclic groups.A CGS Hamiltonian whose gauge group is a direct product can be built up from CGS Hamiltonians with the factor groups as the gauge group.When the gauge group G is a direct product G 1 ×G 2 , the group of gauge transformations and the local groups of gauge transformations are also direct products.This means that if p has combinatorial gauge symmetry with respect to a G 1 gauge theory defined on a lattice and so does H 2 = −J p (µ 2 p ) ⊤ W 2 σ 2 p with respect to a G 2 gauge theory on the same lattice, we can take the tensor product of the two systems, and write where id 1,2 are identity operators acting on the two tensor sectors.This Hamiltonian will then be a G 1 × G 2 combinatorial gauge theory.
Note that even in cases where the gauge group is already cyclic, the tensor product construct could reduce the number of matter spins.For example, a 6 gauge theory on the triangular lattice can be constructed with 36 matter spins, and the number can be reduced to 12 if the matter degrees of freedom have a 3 internal symmetry.If the gauge group 6 is treated as 2 × 3 so that the gauge theory is constructed through the tensor product construction, only 13 matter spins are needed, with 4 matter spins for the 2 sector (see Section 2) and 9 for the 3 sector (see Section 5.2).The coupling matrix for the 3 sector can be further reduced (also see Section 5.2), so that in total only 7 matter degrees of freedom is required, of which 3 need to be 3 variables.

Ground state flux and tuning the spectrum of the CGS coupling term
We note that there is no fundamental constraint on the initial row(s) of the coupling matrix W , so each of these entries can be regarded as a tunable parameter, which gives us a lot of freedom to adjust the spectrum of the CGS coupling term.This is fortuitous, because while the combinatorial gauge symmetry construction guarantees that any model we generate has an exact gauge symmetry, this does not mean that the ground state of our two-body model also corresponds to the ground state of the gauge theory we are constructing, as we saw in the case of a 2 gauge theory on a triangular lattice.To find the energy of a particular flux sector, we can fix a set of value of the σ Z operators that satisfy the flux constraint, then (supposing that the matter variables are spin-1 2 's) choose the signs of the µ z 's, so that the energy is minimized, Due to the combinatorial gauge symmetry, this expression only depends on the flux i∈p σ Z i and the initial coupling values that generate W .Given any specific combinatorial gauge symmetry, these initial couplings can be tuned to find a coupling matrix that has the zero flux sector as the lowest energy one.This argument also shows that when the matter degrees of freedom have internal symmetries, transforming them under such symmetries does not change the spectrum of the CGS Hamiltonian.In other words, if the matter degrees of freedom are spin-1 2 's, for example, flipping the sign of entire rows of the coupling matrix W will not affect the spectrum.
Take the coupling matrix (17) as an example, if a more general set of coupling constants (b 1 , b 2 , . . ., b q ) are taken as the starting point of the construction, the resulting coupling matrix will become W k,q with column j rescaled by b j , and the ground state energy will now become In the W 5,3 [(1, 1, 1)] example mentioned earlier, the ground state contains both flux-0 and flux-1 and sectors.If instead r = (−1, 1, 1) or r = (0.5, 1, 1) is the first row, which corresponds to rescaling the first column of the coupling matrix (40) shown in the appendix by −1 or 0.5, respectively, the ground state will be a flux-0 state.

Examples of combinatorial gauge symmetric Hamiltonians
In examples to follow, we will examine a few CGS Hamiltonians.Through these examples we shall further illustrate both the flexibility of these models and the necessary conditions for realizing a topological state.

Uniform gauge transformations and reduction of number of matter spins
Consider a 2 gauge theory on a square lattice, e.g. the toric code.Using the initial row (−1, 1, 1, 1), we will get the 8 × 4 interaction matrix On the other hand, the 4 × 4 Hadamard matrix with the same initial row was used to construct a 2 gauge theory [22].Both matrices possess combinatorial gauge symmetries corresponding to a 2 gauge theory on a square lattice, but (22) has twice as many rows as (23).We will show that the two constructions are equivalent, and W 4,2 can be reduced to W ′ procedurally.(Thus we may refer to the latter as W which not only permutes the rows of W 4,2 red but also flips the signs of two of them.However, if we relax the restriction that the L matrices be permutations, but also allow for −1 in addition to 1 in its non-zero entries, i.e., by letting L be monomial matrices, we can recover the combinatorial gauge symmetry.Physically, this corresponds to not only permuting, but also flipping the matter spins.This operation commutes with a transverse field, so the 2 gauge symmetry is preserved in the full Hamiltonian.What makes this manipulation possible is that while W 4,2 red does not include every element of the orbit under the action of the gauge transformations, it does include all such entries up to overall signs, so that the action of an R matrix permutes the rows of W 4,2 red up to signs, which can be compensated for by the signs of the entries in the L matrix.In this particular case, the corresponding L matrix is The full and reduced versions of the coupling matrices, W 4,2 and W 4,2 red , are mapped to each other as follows.To reconstruct W 4,2 , append to each row r a ∈ W 4,2 red a row −r a .Conversely, observe that W 4,2 consists of four pairs of rows that differ by an overall sign, so by choosing one representative from each pair, we form W 4,2 red .This procedure also produces a correspondence between the full and reduced versions of L matrices.A reduced L matrix is augmented into the full version by replacing each entry by a 2 × 2 block.Zeros and ones become zero and identity matrices, while −1 becomes a 2 × 2 permutation matrix that swaps the two rows that differ by an overall sign.For example, from L 12  red in (25), we get Compare the transformation of the full W matrix under the action of R 12 , This is a permutation of the rows 1 and 3, 2 and 4, 5 and 8, and 6 and 7, which precisely represented by L 12 .To reduce from a full L matrix, we first sort the rows of the coupling matrix so that rows related by a uniform gauge transformation are adjacent to each other, so that the full L matrix will take a block form.Using a permutation representation of 2 , each block in the full L matrix is mapped to +1 or −1.

Tuning the spectrum of a CGS interaction
Next we shall analyze the case that has been studied in [23], namely a 3 gauge theory on a triangular lattice (see Appendix A for the change in underlying lattice) in the context of the general method of building a CGS model as well as techniques for reducing the number of matter degrees of freedom and adjusting the spectrum.
A plaquette in a triangular lattice has coordination number 3 and we are working with the gauge group 3 , so starting from the initial row (1, 1, 1), we construct a W matrix where ω = e 2πi/3 .Note that the rows can be partitioned into three sets, and within each set the rows differ by an overall phase factor.Following the reduction procedure illustrated by the previous example, we pick one representative from each set to construct a reduced coupling matrix, and take the matter degrees of freedom to be 3 clock variables.However the ground state of this CGS Hamiltonian does not lie in the 3 flux 0 sector, instead it is in the flux 1 and flux 2 sector, which can be verified from equation (20).To fix this issue, we can adjust the initial row of the W matrix, for example to (−1, 1, 1).An alternative solution is presented in the appendix of [23], that is, a W matrix that combines two orbits under the gauge transformations, In Sec. 4 we noted that as long as the rows of W forms an invariant set under the action of the gauge transformations, it possesses a combinatorial gauge symmetry.The simplest invariant sets are orbits of a single row, but they can also combine multiple distinct orbits.This is not obviously advantageous, because the size of such a W matrix is not minimal.However, the above W for the 3 CGS turns out to have the correct flux state as the ground state, while the simple W matrix formed out of the orbit of (1, 1, 1) does not.Thus we have another method of adjusting the spectrum of CGS Hamiltonians.

Gauge theories with more complex structure
To illustrate the power of the CGS construction, we now construct a model with the symmetry of Haah's code.Recall that Haah's code is defined on three-dimensional cubic lattice, where there are two spins on each site.The Hamiltonian is the sum of two types of generators of a stabilizer code, each supported on each cubic cell, c, of the lattice, Each operator A c and B c is a product of σ z 's and σ x 's, respectively, acting on a subset of the eight spins on the cubic cell, as illustrated in Fig. 3. Reformulated as a gauge theory, we may treat the A's as charge operators and the B's a generators of the gauge symmetry.Thus to construct a CGS model with the same type of symmetry, we need to compute the group of local gauge transformations on an A-type operator generated by conjugating it by B-type operators.This is computed in Appendix D. Despite the complex arrangement of sites, the local group of gauge symmetries is very simple, 7  2 .This means that we could emulate an A-type operator by a CGS model of the (2, 8)-type.Figure 4: wire arrays that implement combinatorial gauge symmetries for (a) a 2 gauge theory on a square lattice [25] and (b) a 3 gauge theory on a honeycomb lattice [23].(a) The array is a grid of "waffles" in which horizontal wires extending across neighboring waffles (in black) represent the gauge degrees of freedom, and the vertical wires (in gray) are the matter degrees of freedom.Within each waffle, the wires are coupled by Josephson junctions with no phase shift that corresponding to +1 in the coupling matrix (black cross on white square) and junctions with a π phase shift that correspond to a −1 in the coupling matrix (white cross on black square).(b) In the 3 case, the grid of waffles is constructed by the same principle, but within each waffle, the wires are coupled via Josephson junctions, and the phase shift required by the CGS Hamiltonian is realized by a uniform magnetic field that generates a flux of Φ = (n + 1 3 )Φ 0 within each elementary plaquette of the waffle, where Φ 0 = h/2e is the flux quantum.

Implementations
Proposals have been made for realizing a 2 CGS model on a square lattice with a superconducting circuit [25], and a classical 2 spin liquid with the same combinatorial gauge symmetry was built and observed [26].These implementations suggest generalizations.Following the approach in [25], superconducting wires can be arranged into an array of "waffles", or grids of intersecting wires, as illustrated in Fig. 4(a).Generally, such wire arrays are described by the Hamiltonian (see also [23]) : : : with the Josephson coupling and a quadratic capacitance term H C that depends on the conjugate charges to the phase variables θ i , φ a .When the coupling coefficients W ai are those of a CGS model, the Josephson possesses a combinatorial gauge symmetry, where the action of gauge groups corresponds to shifts in the location of the minima of Josephson energies.(One can ensure that H C does not break the symmetry by requiring permutation invariance of the capacitance matrix.)In the regime where H J is dominant, the minima of the Josephson energy correspond to the minima of the CGS model.
To get the phase factors W ai in the Josephson Hamiltonian, two strategies are available.The junctions can be built as Josephson junctions or π-Josephson junctions [27][28][29][30] when the W matrix consists only of +1 and −1 entries, or by a more general approach given in [23], shown in Fig. 4(b), where the W matrix can be realized by threading a uniform flux of 1  3 flux quantum in every elementary plaquette of the waffle.The latter strategy can in fact be applied to other CGS Hamiltonians.
Generically, the Hamiltonian of a waffle of superconducting wires is the sum of the Josephson energies at the intersections of the wires, whose minima will acquire phase shifts in the presence of a magnetic field.Let θ i be the superconducting phase of the wire representing the i-th gauge variable, and let φ a be the phase of the a-th matter wire, both measured at the edges of the waffle (the left edge for φ a and the top edge for θ i ).The Josephson junction energy at the intersection of the θ i and φ a wires is equal to where is the phase accumulated by integrating the vector potential A, which originates from magnetic fluxes, along the path from the top of gauge wire i to the left of matter wire a (See Fig. 5(a)), and Φ 0 is the flux quantum.A gauge transformation A → A + ∇ϕ shifts the value of A ai , up to a factor of 2π/Φ 0 , by Since these two terms on the right-hand side of (33) only depend on the position where θ i and φ a are measured, they can be absorbed into θ i and φ a , respectively.Thus fixing a gauge amounts to choosing the reference point for the phases in the matter and gauge wires.We choose the gauge such that in which case A ai is tied to the total magnetic flux Φ ai within the area enclosed by the left and top edges of the waffle and gauge wire i and matter wire a: These shifts in the Josephson energy minima can be used to implement the coupling matrix W of a CGS model.To obtain the most general W matrix, all junctions are regular Josephson junctions with the exception of those involving the top wire of the waffle.These latter junctions need to be tunable, so that the Josephson phase is that of W 1i .This can be achieved using the asymmetric DC SQUID design in [25].Then an appropriate magnetic field should be applied, so that at the junction between the θ i and φ a wires a phase shift Φ ai is induced, satisfying This relation gives us a one-to-one correspondence between W matrices and flux configurations in waffles, because excluding the first column and first row, an n × m coupling matrix contains (n−1)×(m−1) entries, and there are precisely (n−1)×(m−1) elementary plaquettes in a waffle with m gauge wires and n matter wires.
For example, the 2 triangular CGS coupling matrix (4) corresponds to the following flux configuration.
, where Φ 0 = h/2e is the flux quantum.Note that a coupling matrix with its rows reordered or its rows transformed according to an internal symmetry of the matter degrees of freedom will have the same combinatorial gauge symmetry and spectrum, but in general it will correspond to a different flux configuration in the waffle realization.For example, if we swap the second and third rows of ( 4) and flip the sign of the second and last row of the resulting matrix, the coupling matrix will correspond to the following flux arrangement (as illustrated in 5(b)), Similarly, while the waffle corresponding to ( 30) is threaded by 1 3 Φ 0 and 2 3 Φ 0 fluxes, switching the fourth and fifth rows allows us to implement it as (see also Fig 5(c)). .
Given the total Josephson energy of a waffle that corresponds to a CGS model with W as its coupling matrix, the minimum is achieved at by a phase configuration (θ 1 , . . ., θ q ) on the gauge wires that corresponds to the ground state of the CGS Hamiltonian.(Compare to (20).)And as argued in [23], the phases (φ 1 , . . ., φ N ) on the matter wires are tethered to the θ i 's via in the same way that the signs of the matter spins in the ground state are fixed by the configuration of the gauge spins.This shows that the ground state degeneracy of the waffle system corresponds exactly to that of the CGS Hamiltonian, so when the Josephson energy is dominant, the superconducting wire array will exhibit a topological phase described by the CGS Hamiltonian.

Summary and outlook
We have presented a general method for constructing a class of one-and two-body-interaction Hamiltonians with an exact gauge symmetry for any given Abelian group.The work systematically encompasses earlier sporadic examples of combinatorial gauge symmetry.We illustrated the principles and the construction for a variety of groups and lattice geometries.
The core physical concept was to build a lattice with matter and gauge degrees of freedom, which can be spins or superconducting wires.The core mathematical concept was to classify the allowed interactions between the matter and gauge degrees of freedom.This problem reduced to ensuring consistency between the interaction matrix and action of the gauge group on the underlying quantum variables.
The Hamiltonians with the exact gauge symmetry should contain gapped topological phases for a wide range of parameters.We are hopeful that the construction here is sufficiently general that one or more of these examples can lead to realizations of topologically ordered quantum spin liquids in engineered or programmable quantum devices.At the very least, the mathematical construction here may offer a direction for further theoretical explorations.gauge theory) and by the NSF Grant DMR-1906325 (work on interfaced topological states in superconducting wires).

A Lattice gauge theory conventions
In this section we make a comment on the lattice gauge theory convention adopted in this paper and contrast it with the conventions in previous work on combinatorial gauge theory [22,23].Briefly, in previous examples of 2 and 3 combinatorial gauge theories, the flux variables of the gauge theory are located on vertices, so the flux operators are supported on stars, of a square lattice and a honeycomb lattice, respectively.Under the conventions of this paper, these constructions would be for a square plaquette-based 2 flux operator, and a triangular plaquette-based 3 flux operator.
In 2D Abelian lattice gauge theories, the flux and charge variables are dual to each other.Under a transformation that maps each face of graph to a vertex of its dual graph and vice versa, the flux variables can be interpreted as the charge variable of a theory defined on the dual graph.This means that, in terms of physical properties, it makes no difference if the flux variables are supported on the plaquettes of a graph or on the vertices.The more popular convention is to put flux variables on plaquettes, which is also the convention adopted in this paper.See Fig. 6 for an illustration of the 2 gauge theory in Section 2. Note that the position of the degrees of freedom and the Hamiltonian do not depend on the underlying lattice chosen.
We also note that the convention of situating flux variables on vertices of the lattice is the source of the terminology "matter" degrees of freedom.These matter variables are directly coupled to the gauge degrees of freedom, and are analogously supported on the vertices of the lattice.While strictly speaking individual "matter" variables don't carry representations of the gauge group, their Hilbert spaces may be rearranged into a direct sum of subspaces that do.Regardless of the true nature of these matter variables, in this paper we are interested in the confined phase of the gauge theory, which is most straightforwardly achieved when the matter variables are integrated out, so we need not dwell on their role in the lattice gauge theory.

B Further examples of CGS coupling matrices
Here we show two further examples of CGS coupling matrices.
For a 5 gauge theory on a triangular lattice, the following coupling matrix is constructed from the initial row (1, 1, 1), Here ζ = e 2πi/5 is one of the fifth roots of unity.The group of gauge transformations are generated by where perm{σ} denotes the permutation matrix corresponding to the permutation σ.
For a 2 gauge theory on a triangular lattice, using the initial row (1, 1, 1, 1, 1, 1) we can construct a 2 6−1 × 6 = 32 × 6 CGS coupling.However, since gcd(2, 6) = 2, we can reduce the number of rows of this matrix by a factor of 2, to the following 16 × 6 matrix, The generators of the gauge transformation can be taken to be  This example corresponds to the color code [31] on a honeycomb lattice.

C Breaking spurious degeneracies for 2 CGS theories
For a 2 CGS theory on a lattice with an odd coordination number q, the basic construction outlined in Sec. 4 gives us a 2 q−1 × q W matrix.The rows of W consist of q-element vectors of ±1 that contain an even number of −1's.To find the ground state of this CGS Hamiltonian, we need only minimize the energy under the constraint σ = 1, since all other flux 0 states are gauge symmetry partners of this state, and the flux 1 states can be obtained by a global spin flip.In this case, the Hamiltonian reduces to − a µ a i W ai .This is minimized when each µ a takes the sign that makes µ a i W ai negative, which is in turn determined by the sign of i W ai .Among these, there are q 2 j rows with 2 j negative signs.When 2 j < q/2, the sum of all entries in the row is positive, otherwise it is negative.Thus in the ground state, the total µ-magnetization is ⌊(q−1)/4⌋ j=0 q 2 j − (q−1)/2 j=⌈(q−1)/4⌉ q 2 j .
Note that the first sum is over the even binomial coefficients up to (q − 1)/2, and the latter sum is over the remaining even binomial coefficients.Since q 2 j = q q−2 j and q is odd, the latter sum is equivalent to a sum over the odd binomial coefficients up to but not including the (q − 1)/2-th one.This means that we can rewrite the total µ-magnetization as (q−1)/2 j=0 (−1) j q j .
This shows that the µ-magnetization in the CGS Hamiltonian for k = 2 and q odd is always non-zero, thus guaranteeing that a uniform field on the matter spins will serve to split the ground state degeneracy between flux 0 and flux 1 states.Since B is a product of Pauli x-operators, the matrices R are diagonal matrices with ±1 along the diagonal.Thus for A to be under this transformation, there must be an even number of −1 along the diagonal.In other words, the local group of gauge symmetries must be a subgroup of the group consisting of operations flipping an even number of the spins in (Z 1 , . . ., Z 8 ), which we shall denote by G.A generating set of G can be taken to be {P i,i+1 } i=1,..., 7 , where P i,i+1 can be represented by the diagonal matrix This also shows that G ∼ = 7 2 .Now by showing that every element in this generating set can be expressed in terms of symmetry operations generated by type-B operators, we can prove that the local group of symmetries G A is exactly equal to G. The generators of G and their decompositions in terms of B-operators are tabulated in Table 1.Note that there are 14 type-B operators that generate non-trivial symmetries of A, so this computation also shows that they are not all independent of each other.

Figure 1 :
Figure 1: Illustration of the 2 gauge symmetry on a triangular lattice.(a) A gauge transformation flips certain spins, but plaquette terms remain invariant overall, since there is always an even number of spin flips.(b) A two-body operator satisfying the combinatorial gauge symmetry corresponding to the gauge symmetry depicted in (a) requires four matter degrees of freedom.(c) Under the action of the gauge symmetry in (a), the combinatorial gauge symmetric operator in (b) stays invariant by flipping the gauge spins just as in the pure gauge theory, while also permuting the matter spins.

Figure 2 :
Figure 2: Two gauge transformations T and T ′ that represent the same element of the local group of gauge transformation (magenta) for a plaquette operator (highlighted in black).
that is, their actions coincide on the support of the operator O.We define this group as the local group of gauge transformations G[O] for the operator O. (See Fig. 2 for an illustration of the relation between the group of gauge transformations and the local group of gauge transformations in the 2 gauge theory of Section 2.) For compactness, when O is the flux operator B p , we shall write G p instead of G[B p ].
incomplete orbit under gauge transformations might suggest that W 4,2 red cannot have a combinatorial gauge symmetry.Indeed, multiplying W 4,2 red on the right by R 12 = diag{(−, −, +, +)} we get W

Figure 3 :
Figure 3: Illustration of the two types of operators in the Hamiltonian of the Haah's code.

Figure 5 :
Figure 5: Flux-based superconducting wire array realizations.The general configuration of gauge and matter wires in a waffle and the phase shift through a Josephson junction induced by a magnetic field is depicted in (a).Correspondences between CGS coupling matrices and flux configurations are shown in (b) for the 2 CGS model on a triangular lattice and a waffle with a uniform flux Φ 2 = 1 2 Φ 0 except for one plaquette, and (c) for the 3 CGS theory with the coupling matrix W 3,3 red [(1, 1, ω), (1, 1, ω)] and a waffle with a uniform flux Φ 3 = 1 3 Φ 0 except for one plaquette and one junction (represented by a dotted box) that has a 2π 3 phase shift.

1 Figure 6 :
Figure 6: A gauge theory on a triangular lattice where the flux operators are supported on plaquettes (blue lattice) can be reinterpreted as a theory on the honeycomb lattice with flux operators supported on the vertices (yellow lattice).

Figure 7 :
Figure 7: Notations used for writing down the generators of the local group of gauge symmetries.The sites that are acted on non-trivially by type-A operators are numbered from 1 through 8 as in the diagram.The type-B operators that generate nontrivial transformations of one particular A-operator are labeled by their coordinates relative to it, in the coordinate system illustrated on the right.

D
Group of local gauge symmetries of the A-type operators in Haah's codeTo describe the local group of gauge symmetries G A of a type-A operator in the Haah's code Hamiltonian, it is convenient to label sites that it acts on as in Fig.7.The gauge symmetry generated by a type-B operator can then be written as the a multiplication by a matrix R,

Table 1 :
Generators of the group of local gauge symmetries.