Fractonic topological phases from coupled wires

In three dimensions, gapped phases can support"fractonic"quasiparticle excitations, which are either completely immobile or can only move within a low-dimensional submanifold, a peculiar topological phenomenon going beyond the conventional framework of topological quantum field theory. In this work we explore fractonic topological phases using three-dimensional coupled wire constructions, which have proven to be a successful tool to realize and characterize topological phases in two dimensions. We find that both gapped and gapless phases with fractonic excitations can emerge from the models. In the gapped case, we argue that fractonic excitations are mobile along the wire direction, but their mobility in the transverse plane is generally reduced. We show that the excitations in general have infinite-order fusion structure, distinct from previously known gapped fracton models. Like the 2D coupled wire constructions, many models exhibit gapless (or even chiral) surface states, which can be described by infinite-component Luttinger liquids. However, the universality class of the surface theory strongly depends on the surface orientation, thus revealing a new type of bulk-boundary correspondence unique to fracton phases.


I. INTRODUCTION
It is a common belief that gapped phases of quantum many-body systems can be described by topological quantum field theories (TQFT). There is strong evidence to support TQFT-based classifications in one and two spatial dimensions, but in three dimensions, a large family of gapped topological states have been discovered theoretically [1][2][3][4][5][6], whose properties do not easily fit within the TQFT framework. They all feature quasiparticle excitations with reduced mobility, and sometimes even completely immobile ones. When put on a three torus, these models exhibit a topological ground state degeneracy (GSD) which, asymptotically, grows exponentially with linear system size. These exotic quantum states are now known as fracton phases [7,8].
In searching for the unified structure underlying fracton order, it was realized that certain (type-I) fracton models can be built from coupling stacks of 2D topological phases [9,10]. In fact, a simple stack of layers of 2D gapped states already exhibits some of the characteristics of fracton topological order, e.g. size-dependent GSD, quasiparticles with reduced mobility (i.e. they can only move in planes). This observation has inspired more general constructions of fracton topological order [11][12][13][14][15][16][17], as well as shedding light on the precise meaning of fracton phases [18][19][20]. More recently, systematic constructions of fracton models from networks of 2D TQFTs (possibly embedded in a 3D TQFT) have been proposed [21][22][23], encompassing many existing models including type-II ex-amples. Ref. [21] further conjectured that all 3D gapped phases can be obtained this way.
In this work, we explore a different avenue to construct 3D fracton topological phases. The starting point of our construction is an array of 1D quantum wires, each described at low energy by a Luttinger liquid. Couplings between wires are then turned on to lift the extensive degeneracy. This approach is known as a coupled wire construction [24,25] and has been widely applied to construct explicit models for 2D topological phases (for a recent survey of these applications see Ref. [26]). The advantage of the coupled wire construction is that the correspondence between the bulk topological order and the edge theory can be made very explicit and chiral topological phases arise naturally. While being more general than other exactly solvable models (i.e. string-net models), the construction still remains analytically tractable. The method has previously been applied to 3D systems, however the examples have fallen under TQFT-type topological orders [27,28]. Related ideas of building fracton phases out of coupled spin chains have also been explored in Refs. [29,17].
In this work we systematically study 3D coupled wire construction. We find that the construction can easily give rise to both gapped and gapless phases. In both cases, there can exist gapped fractonic excitations. Notably, these excitations generally have fusion structures distinct from all previously known gapped fracton models: they are labeled by an integer-valued (internal) topological charge, much like electric charges in a gapless U(1) gauge theory, even though the system can be fully gapped. When it is gapped, we find that the excitations are mobile along the wires (although require a quasi-local string operator), but generally have reduced mobility in the transverse directions. We demonstrate in a class of examples that all excitations are immobile in the transverse directions, thus exhibiting a new type of lineon topological order.
Similar to the 2D case, the coupled wire construction allows direct access to surface states, at least when the surface is parallel to the wires. The surface states can be described by an infinite-component Luttinger liquid, defined by an infinite K matrix, which would ordinarily correspond to an (infinite-component) Abelian Chern-Simons theory in the bulk. However, the bulk-boundary relation is highly unusual in the gapped lineon phases: in most cases the surface K matrix strongly depends on the orientation of the surface, so in a sense different surfaces can host qualitatively different gapless states.
The paper is organized as follows: In Sec. II we provide a systematic overview of the coupled wire construction. We classify excitations in a general coupled wire model, emphasizing the important role of Gaussian fluctuations overlooked in previous studies. A polynomial formalism is introduced for these models which allows for the use of powerful algebraic methods. In Sec. III we consider a class of models built from wires hosting a single Luttinger liquid, which we term the chiral plaquette mod-els. There prove to be both gapped and gapless models within this class. Specific examples of each case are analyzed. In Sec. IV an example with each wire hosting a two-component Luttinger liquid is considered. These "CSS" models are shown to be gapped, with all quasiparticles being lineon. In Appendix A we give a general treatment of 2D coupled wire construction. Appendix B discusses the spectrum of Gaussian fluctuations for 3D coupled wire models. This analysis shows that models which naively appear gapped can actually posses gapless fluctuations. Appendix C gives an algorithm for computing the GSD of these models. In Appendix D we give an algorithm for finding the charge basis of the models. Finally, in Appendix E we prove that all of the excitations of the CSS models of Sec. IV are lineons.

II. COUPLED WIRE CONSTRUCTION
In this section we lay out the general theory of coupled wire construction. While our focus is on the 3D case, the formalism applies to 2D without much modification and in fact is more tractable there. We provide a detailed account of the general theory in 2D in Appendix A.
Consider quantum wires arranged in a square lattice, as illustrated in Fig. 1(a). Each wire is described by an M component Luttinger liquid, with a K matrix K w . If the wire is bosonic (fermionic), we take K w = σ x ⊗1 M ×M (K w = σ z ⊗ 1 M ×M ). The Lagrangian for one wire is Throughout this paper we choose the wires to extend along the x direction. V is the velocity matrix, which we take to be V = v1 for simplicity. Here Φ denotes the bosonic fields collectively Φ = (φ 1 , φ 2 , . . . , φ 2M ) T .
They satisfy canonical commutation relations Since K −1 w = K w , we do not distinguish the two throughout the paper. Denote the bosonic field of the wire at site r = (j, k) by Φ r (x). Bosonic fields on different wires commute. Importantly, all these fields are 2π periodic, which means that local operators in the theory are built out of derivatives of φ's, together with vertex operators of the form e il T Φr(x) with l being an arbitrary integer vector.
We add the following type of interactions to gap out the wires: Each of the Θ α r (x) is a linear combination of fields at nearby sites, with integer coefficients. We demand that the system is translation-invariant, including the continuous translation along x and the discrete translations along y and z. We write Θ α . The gapping terms should satisfy the null conditions [30][31][32][33]: This guarantees that the each Θ α r can be simultaneously frozen out in the large U limit.
Additionally, we demand that Θ α r are asymptotically linearly-independent, so there are sufficiently many of them to pin all bosonic fields. More precisely, consider the array of wires with periodic boundary conditions, the total number of independent Θ's should be M N w − c where N w is the number of wires, and c is bounded (c may vary with system size). This means that there should be no local linear relations for Θ's, so all possible linear relations involve infinitely many fields in the thermodynamic limit.
We further assume that the gapping terms are "locally primitive", which roughly says there are no local order parameter. More generally, we impose the condition that there exist no nontrivial local fields that commute with all Θ's, except the Θ's themselves and their linear combinations with integer coefficients. This is analogous to the topological order condition for stabilizer codes in qubit systems. Therefore, when U is sufficiently large, the gapping terms freeze all bosonic fields, at least at the classical level.

A. Energy spectrum and excitations
We now give a semi-classical description of the lowenergy excitations. In the limit of large U , the cosine potentials pin all Θ r at the minima and so the ground state manifold corresponds to the configurations with Θ r ∈ 2πZ. Excitations can be classified into the following two types.
The first type of excitations correspond to small oscillations of the fields Θ r around the minima. The spectrum of such oscillations can be found by a Gaussian approximation, i.e. expanding cos Θ r ≈ 1 − 1 2 Θ 2 r at Θ r = 0 and diagonalizing the resulting quadratic Hamiltonian (more details can be found in Appendix B). Physically, these excitations can be interpreted as density waves. We note that the spectrum of such oscillations may be gapless or gapped depending on the form of the gapping interactions. One may wonder whether the conditions on gapping terms listed in the previous section imply that the Gaussian fluctuations are gapped (which might be implicitly assumed in most coupled wire constructions in the literature, where such Gaussian fluctuations were not considered), but this is far from obvious and most likely incorrect.
The other type of excitations are known as "kinks" or "solitons", where some of the fields Θ r tunnel from one minima to another. More concretely, for the gapping term cos Θ l r (x), a n-kink where n ∈ Z is a configuration where Θ α r varies by 2πn over a short distance ξ. Such kink excitations are localized, and may be topologically nontrivial ( i.e. cannot be created locally). They represent gapped quasiparticle excitations, which are of most interest to us in this work.
To characterize the quasiparticle excitations, in particular their superselection sectors as well as the mobility, it is important to understand the structure of local excitations, i.e. how local operators act on the ground state.
As mentioned in the previous section, local operators come in two types: spatial derivatives of Φ and vertex operators. First, a local vertex operator e il T Φr(x0) acting on the ground state generally creates multiple kinks in the yz plane located at x = x 0 . From the commutation algebra, one can easily find that a (l T K w Λ r−r )-kink is created for the Θ r gapping term. It should be clear that such patterns of local creations of kinks determine the mobility of quasiparticles in the yz plane, i.e. the plane perpendicular to the wire direction.
The other type of local operators, namely derivatives of Φ, can be used to transport excitations along the wire direction. In particular, one can construct string operators of the general form Here w r can be real numbers. Without loss of generality, they can be restricted to [0, 1) as the integral part corresponds to a local vertex operator. For Eq. (6) to be a legitimate string operator, it should create allowed kink excitations (of strength 2πZ) at x 1 and x 2 . We further demand that w r should be "quasi-localized" in the yz plane, either strictly short-range, or exponentially localized.
For a coupled wire model in two dimensions, under quite general conditions we are able to completely classify superselection sectors (i.e. anyon types) of excitations and show that they are generally given by the determinant group of an integer K matrix determined from the gapping terms. We further show that lattice translations across wires can permute anyon types. Unfortunately for 3D models, we do not have results of similar generality and will have to work case-by-case.

B. Ground state degeneracy
In the limit of large U , the ground state of the coupled wire model on a torus is obtained by minimizing the gapping potentials cos Θ α r = 1, in other words Θ α r = 2πn α r where n α r ∈ Z. The 2π-periodicity of the underlying bosonic fields φ's induces equivalence relations between different configurations of Θ α r . For example, shifting φ i at site r by 2π leads to the following shifts of n's: Thus these two configurations are actually equivalent. Physically distinct ground states then correspond to equivalence classes of the integer configurations {n α r }. A more systematic algorithm to compute the ground state degeneracy is presented in Appendix C. The degenerate space is spanned by non-local string operators running along x and surface operators in the yz plane, referred to as logical operators, and therefore topologically protected as least when the system is fully gapped.

C. Surface states
With open boundary conditions, generally there exist boundary local fields unconstrained by the bulk gapping terms, which may form gapless surface states. This is most easily seen when the surface is parallel to the wire direction. For such a surface, "free" surface fields can be obtained as "incomplete" gapping terms, which by construction commute with the bulk ones. These incomplete gapping terms however do not commute with each other in general. Denote the corresponding fields on the boundary byΘ n (x), where n ∈ Z indexes the transverse direction on the surface. Their commutation relation generally should take the form Here K surf is an integer symmetric matrix. The algebra is formally equivalent to that of local bosonic fields in a multi-component Luttinger liquid with K surf as the K matrix. As a result, the surface degrees of freedom can be viewed as a 2D extension of a Luttinger liquid. Interestingly, in general K surf depends on the surface orientation, which is a very unusual form of bulk-boundary correspondence.

D. Polynomial representation
Here we introduce a polynomial representation for the coupled wire construction, inspired by Haah's polynomial formalism for Pauli stabilizer codes [34,35]. It gives a compact form of the gapping terms and allows applications of powerful algebraic methods.
We denote using y and z, the unit translations along the two directions perpendicular to wires, and byȳ = y −1 ,z = z −1 , the inverse translations along the same directions respectively. The Hamiltonian can be represented as a 2M × M matrix, where each column represents one Θ α . Rows of a given column give the Laurent polynomial for the corresponding bosonic field that shows up in Θ α . More explicitly, suppose Λ α rr ≡ Λ α r −r as required by translation invariance, the p-th row of the column is where 1 ≤ p ≤ 2M, 1 ≤ α ≤ M . Following Ref. [34], σ is called a stabilizer map (and each gapping term is a stabilizer). The null condition can then be summarized as where σ † ≡σ T . This is again reminiscent of Haah's formalism for stabilizer codes, except that now the polynomials have Z coefficients, instead of Z 2 . The other important difference is that K is symmetric, not symplectic.
We also define = σ † K w as the excitation map. The rows in the excitation map correspond to the different gapping terms. Acting with the excitation map on local operators reveals the excited gapping terms i.e. the kink excitations. In other words, it is a map between local operators and the kink excitations created by the action of those local operators, hence the name. The null condition in Eq. (10) implies Im σ ⊂ ker . We further require that Im σ = ker on an infinite system, which means that there are no gapless degrees of freedom left, which is the primitive condition discussed previously.
We study two families of stabilizer maps. For the first family, we consider M = 1 and K w = σ z . The stabilizer map is given by: σ = m 1 + n 2 y + n 1 z + m 2 yz m 2 + n 1 y + n 2 z + m 1 yz .
We refer to these stabilizer maps as chiral plaquette models.
The second general family of models is defined for even M and bosonic wires. We use M = 2 as example with K w = σ x , and denote the two bosonic fields on each wire as φ and θ. The stabilizer map is given by where f and g are finite-degree polynomials of y and z.
The Hamiltonian is similar to the "CSS" codes for Pauli stabilizer models, in the sense that one term only involves φ and the other only involves θ. The polynomial formalism, besides providing an economic representation of the Hamiltonian, allows for the use of powerful mathematical tools from the theory of polynomial rings. We use the representation to compute the basis of nontrivial superselection sectors or the "charge basis" of a given model. In other words, the charge basis is defined as the set of quotient equivalence classes of charges modulo trivial charge configuration. A charge configuration is trivial if and only if it can be created out of the ground state (e.g. two charges created by applying a string operator) using a local operator. For conventional topological phases, the charge basis is a finite set. In contrast, a fracton phase necessarily has an infinitely large charge basis. An efficient algorithm to compute the charge basis is described in Appendix D along with the calculation for some models.

III. CHIRAL PLAQUETTE MODELS
The chiral plaquette model is illustrated in Fig. 1(b). Specifically, the gapping term is given by Here we define m = (m 1 , m 2 ), n = (n 1 , n 2 ) and m = (m 2 , m 1 ), similarly for n . The components m 1,2 , n 1,2 are integers. Using the single wire K matrix K w = σ z , we define a dot product between vectors as x · y = x T K w y, for example m · m = 0. We also write In the following we assume gcd(m 1 , m 2 ) = gcd(n 1 , n 2 ) = 1. We further assume that m and n are linearly independent, the same with m and n . Equivalently, both m · n and m · n must be non-zero.
In Appendix B 3 we find the density wave spectrum where f k = m 1 + n 2 e iky + n 1 e ikz + m 2 e i(ky+kz) .
The spectrum is fully gapped if and only if Otherwise there are gapless points, near which the dispersion is linear.

A. Surface theory
First we study the surface states. Consider the (0, 1, 0) surface, shown in Fig. 2(a), which is parallel to the xz plane. The bulk occupies the y > 0 region. We follow the general procedure outlined in Sec. II C to find the K matrix. At y = 0, consider vertex operators supported on two adjacent wires e i(l T 1 Φz+l T 2 Φz+1) . If they commute with all bulk terms, l 1 and l 2 should satisfy l 2 · m = 0, l 1 · n = 0, l 1 · m + l 2 · n = 0.
We find that the only nonzero solution is l 2 = m , l 1 = n . Indeed, these terms can be viewed as "half" of a hypothetical plaquette term on the boundary. DenoteΘ z = n T Φ 0,z +m T Φ 0,z+1 . It is not difficult to show that any local vertex operators can be expressed as linear superposition of Φ z . Therefore they form a basis for local vertex operators. Their commutation relations define the K matrix on the surface, according to Eq. (8). The nonzero entries of the K matrix are One can similarly find the K matrix on the xy surface.
It is easy to see that all one needs to do is to replace n and n . Next we turn to the (0, 1, 1) surface, illustrated in Fig.  2(b). Due to the zigzag shape, there are two kinds of vertex operators surviving at low energy, which can again be viewed as "incomplete" plaquette terms, as illustrated in Fig. 2(b). The surface K matrix reads: which is generally distinct from the K matrices on the other surfaces. This example illustrates that in general surfaces of different orientations have very different K matrices. Since the surfaces are gapless, one may wonder whether they have any instabilities. The stability can be analyzed a la Haldane [30] by attempting to find a complete set of null vectors to gap out the edge, for example. We just point out that in some of our examples, the (0, 1, 0) and (0, 0, 1) surfaces are stable with respect to any local perturbations, and in fact fully chiral (i.e. all modes are moving in the same direction).

B. Mobility of fundamental kinks
Below we consider the mobility of the most fundamental excitation, a kink on a single plaquette. Following the discussions in Sec. II A, we first consider mobility in the yz plane, and then mobility along x, the wire direction.

Mobility in the yz plane.
Determining if a quasiparticle is mobile is equivalent to finding a string operator for the excitation. We now show that it is sufficient to consider string operators of a minimal width, i.e. one site. To see this, we consider a string operator of width 2 creating a single plaquette excitation at one end, illustrated by the following figure: r+ẑ r Suppose the operator at the corner site r+ẑ is e il·Φ r+ẑ . Because it should only excite the plaquette on the bottom left, but not on the top left, we must have l · n = 0, which means l = an for some a ∈ Z. We thus multiply the string operator by a stabilizer e −iaΘr to clean the operator at site r +ẑ. This "cleaning" can be repeatedly applied to the other sites on the same line, so now the width is reduced to just 1. It should be clear that a similar argument works when the width is greater than 2.
First we study a string operator directed along z direction, which can be written as z e il T z Φyz (20) where l z are integer vectors. If the string commutes with the Hamiltonian, the following conditions must be satisfied: We write l k = a k m + b k n, which is always possible when m · n = 0. The second equation gives b k+1 = a k . The first equation gives the following recursion relation which can be written as (when m · n = 0) A general solution can be expressed as powers of the roots of the characteristic polynomial: Here A and B can be fixed by initial conditions of the sequence {a k }. Notice that when |λ z | > 2, one of the two roots has absolute value greater than 1 and the other smaller than 1, so |a k | generally grows exponentially with k.
For concreteness, suppose that at the bottom of the string one finds a single kink excitation. This corresponds to the initial condition l 1 = m, which excites a (m · n)kink. Thus the whole string is fixed by . The coefficients a k should all be integers. This is possible only when λ z is an integer: 1. |λ z | = 0: the string repeats with period-4 m, n, −m, −n.
Only for λ z = 0, ±1, can the string operator actually move the excitation, at least for appropriate lengths. When |λ z | ≥ 2, the excitation created at the end of the string operator has a strength growing with the length of the string, thus costing more and more energy as the distance increases. Therefore we conclude that a (m · n)kink can move along z only when λ z = 0, ±1.
For the mobility along y, a similar calculation can be done and the ratio λ y = −m 2 +n 2 m ·n determines the mobility (basically m → m ).
For certain choices of m and n, e.g. when n = n , it is necessary to consider string operators along the z = ±y diagonal directions. We study an example of this type below in Sec. III C.

Mobility along x
We now consider the mobility along the wire direction. As discussed in Sec. II A, in general kinks can be moved along x with the following string operator: where w r are real numbers. Locality requires that these w r must have a (quasi-)localized profile, i.e. either completely short-range, or decaying exponentially away from the location of the excitation. Below we present a method to construct such a quasi-local string operator, whose profile in the yz plane is strictly short-range in one direction, but only quasi-localized in the other direction. In fact, using a cleaning argument similar to the one in Sec. III B 1, one can prove that no string operator of strictly finite support in the yz plane can move a single plaquette excitation in the wire direction, so a string operator, if it exists, must be quasi-localized. The construction is most easily explained when the system is compactified to a quasi-2D one. Such a process has been used in Ref. [36] to understand properties of fracton models. Without loss of generality, we choose the compactification direction to be z, i.e. a periodic boundary condition is imposed along z. Denote the length of the z direction by L z . Basically, a whole column of wires at a given y is viewed a "super"-wire, with L z bosonic degrees of freedom. Denote these fields on one "super"-wire collectively by Φ y , and we write the gapping Hamiltonian in the following way: It is straightforward to read off P z 's from Eq. (13). In particular, the K matrix K xz can be expressed as [K xz ] zz = P z · P z . We consider the following form of string operator: One can prove that it only creates excitations at the y-th "super"-wire. It is shown in Appendix A that to move an elementary excitation located at the (0, z 0 ) plaquette, we can set In order to have a quasi-localized string, u z must be sufficiently localized. For |λ z | > 2, u z decays exponentially where ξ ∼ cosh −1 |λz| 2 . Because of the exponential decay, the choice of of boundary condition is inessential.
If |λ z | < 2, the string operator is only algebraically localized (or even worse), and for certain system size the K matrix K xz can become degenerate. Therefore the construction does not yield a quasi-localized string operator for a fundamental, single-plaquette excitation (it is possible that for composite excitations a localized string operator still exists).
In principle, the construction can be applied to compactification along any direction, and a quasi-localized string operator can be defined as long as the corresponding K matrix is "gapped". Here, a gapped K matrix means that the eigenvalues of K are separated from 0 by a finite spectral gap in the infinite size limit. Otherwise the K matrix is said to be gapless. If for any compactification direction the K matrix is gapless, we believe that the bulk must be gapless (i.e. the condition Eq. (16) must be violated).

C. Example of gapped phases
We now consider two examples of the chiral plaquette model with fully gapped bulk. Both examples are not covered by the general discussions in the previous section due to the special choices of m and n, and in fact turn out to be planon models.

Planon phase I
Consider m = (p + 1, p), n = (1, 1). One can easily check that the bulk is fully gapped. The K matrix reads which is also gapped since λ z = 2p + 1 > 2. In this model, the 1-kink can move along y = z lines. The string operator is given by From the analysis in Sec. III B, they can also move along the wires, but not along y or z directions. Therefore all excitations are at least planons. We conjecture that there are no fully mobile excitations in this model. Consistent with the planar structure, since n 2 = n 2 = 0, the K matrices on all surfaces different from (0, 1, −1) are actually the same (up to an overall sign). On the other hand, the surface (0, 1, −1) can be made completely gapped: the K matrix becomes Therefore, boundary fieldsΘ 2j , j ∈ Z mutually commute and form a complete set of null vectors. The surface can be gapped by the following interactions: From these results, it is plausible to conjecture that the model realizes a pure planon phases in the (0, −1, 1) planes. Such a phase can be described by an infinite Chern-Simons theory [37] with a K matrix given in Eq. (30). To further check this conjecture, we compute the ground state degeneracy for a L y × L z grid of wires, with periodic boundary condition imposed. We find that Here l = gcd(L y , L z ) is the number of (effective) 2D planes with this boundary condition, which naturally explains the degeneracy except the special gcd(L y , L z ) = 1 case.

Planon phase II
Consider m = (p, 0), n = (0, 1), with p > 1. The bulk is fully gapped according to Eq. (16). The surface K matrices can be easily found: and The K matrix on (0, 1, 1) surface is block diagonalized. Each block is the following 2 × 2 matrix: It is not difficult to see that the kink charge on each plaquette is defined mod p 2 −1, namely charge (p 2 −1) can be created locally. An elementary 1-kink can move along y, with a period-2 string operator, and the mobility along x is also clear from K xz . We have computed the charge basis in Appendix D, and the result, The GSD is found to be Note that (p 2 − 1) Lz = det K xz (L z ). The reduction for odd L y can be understood as follows: under T y , n-kink becomes −pn, mod p 2 −1. T y -invariant kinks must satisfy n ≡ −pn mod p 2 − 1, or n is a multiple of p − 1, which form a Z p+1 subgroup. While the degeneracy and the mobility of bulk excitations seem to be compatible with the model describing decoupled xy planes, this picture is inconsistent with the fully chiral surface theory given by Eq. (36) in the xy plane (it would be gappable if the decoupled layer picture was correct).

D. Examples of gapless phases
In this section we study two examples of chiral plaquette model, whose Gaussian fluctuations are gapless. We will analyze the properties of the gapped sectors, and leave the effect of gapless modes for future work.

Gapless fracton phase
Consider m = (p, q), n = (−p, q). We assume p and q are coprime, |p| = |q| and both nonzero. The Gaussian spectrum is found to be gapless.
It is useful to first give the surface K matrices. For the xz surface, the nonzero elements are For the xy surface: (40) Similar to the K matrix in Eq. (32), the xy surface can be fully gapped.
One can easily see, from previous discussions that a 1-kink is immobile in the yz plane (a 2pq-kink can move along y). Since both surface K matrices have |λ| < 2, the construction in Sec. III B 2 fails to produce a quasilocalized string operator (only an algebraically-localized one). As mentioned in Sec. III B 2, since gcd(p, q) = 1 no string operator of bounded support in the yz plane can move this excitation. With this body of evidence we conjecture that the 1-kinks are fractons. As we show below, they can be created at corners of a rectangular sheet operator in the yz plane.
The GSD has interesting dependence on L y and L z . We will specifically consider the case when p + q is odd an L z is even, to illustrate the physics. From numerical computations we find the GSD is given by the following formula: . (41) An interesting feature of the GSD is that when L y is a multiple of 4, the additional factor | detK xz | appears in the GSD. HereK xz will be defined below, but it is related to K xz by an order one factor, and grows exponentially with L z . We will now explain the exponential factor (2pq) Ly as well as the dependence on L y mod 4. We note that the odd L z case has a similar L y dependence.
xz planons. The exponential dependence on L y can be understood in terms of planons in the xz plane. These excitations take the form of two 1-kinks 1 y,z 1 y+2,z . Here we label plaquettes by the coordinate of the down-left corner.
First we show that 1 y,z 1 y+2,z can move along z. We explicitly construct the string operator. Place the two kinks at (0, 0) and (2, 0), and the string is supported on y = 1 and y = 2: where It is straightforward to check that the string operator commutes with the gapping terms. The initial condition x 1 = (a, b) should satisfy qa − pb = 1, which is always solvable over Z if gcd(p, q) = 1. We notice that the string operator has period-2 along z. Furthermore one can show that (2pq) y,z (2pq) y+2,z is a local excitation. This should be evident from the following figure: Here (a, b) denotes a vertex operator e i(aφ1+bφ2) . Thus these planons satisfy Z 2pq fusion rules. This construction also immediately shows that 1 y,z 1 y+2,z can move along the x direction. From now on we refer to 1 y−1,z 1 y+1,z as the y-th planon. From the string operators one can easily determine the braiding statistics of these planons. We find that they are all bosons, and there is a e iπ pq mutual braiding phase between neighboring planons.
We note that with the planon string operators can be easily extended to a rectangular sheet operator that creates four 1-kinks on the corners, similar to what happens in the X-cube model.

Compactification.
To understand the L y dependence of the GSD, it is useful to consider compactification along z. Closed string operators wrapping around the z cycle for the xz planons become local order parameters, which must be fixed under compactification. We give their explicit expressions: We may view W 1 as moving the planon q, q around the z cycle, and W 2 as moving the −p, −p planon. Notice that these two string operators are not totally independent, due to the relation qW 1y + pW 2y ≡ 0 (mod local stabilizers). Diagonalizing these order parameters, the theory is partitioned into different sectors labeled by eigenvalues of W 1y and W 2y . They naturally form a Z 2q × Z 2p group. We now prove that with the additional condition it is in fact Z 2pq . In fact, consider the element (1, 1), whose order is obviously 2pq. Now for a general element (a, b) in Z 2q × Z 2p , we look for an integer n such that n ≡ a mod 2q, n ≡ b mod 2p. In other words, there must exist x 1 , x 2 ∈ Z such that n = a + 2qx 1 = b + 2px 2 , or a−b = 2(px 2 −qx 1 ). Since a+b is even and gcd(p, q) = 1, one can always find x 1 and x 2 to satisfy this equation. In this compactified system, we have to fix the local "order parameters", which are generated by W 1y and W 2y when summing over all z. In analyzing the quasi-2D system, it is convenient to just add another term − cos W 1y or − cos W 2y . For simplicity, consider p = 1, then we just need to add − cos W 1 . A linearly independent set of gapping terms are W 1y , Θ y,z , z = 1, 2, · · · , L z − 1. With this choice of gapping vectors, the local degeneracy is completely removed. We denote the new K matrix computed from this set of gapping vectors asK xz .
The GSD of the compactified system reads: which almost replicates the L y dependence of the GSD in Eq. (41), up to a factor of 2pq when L y is a multiple of 4 due to an additional relation among the compactified planon string operators. The L y dependence here can be traced to the ytranslation symmetry action on anyons in the compactified system. It turns out that | detK xz | is always a perfect square, so denote M = | detK xz |. The fusion group of Abelian anyons turns out to be Z 2 M . Denote the unit translation action on anyons by T y . We show quite generally that T 2 y = −1 (i.e. the charge conjugation). If we choose the basis to be (1, 0, . . . , 0) and (0, 0, . . . , 1), we have numerically found that the translation action along y is given by the following SL(2, Z) matrix: which satisfies T 2 y = −1, so T y is order-4. The only order-4 element for SL(2, Z) is in fact the S matrix, so T y is in the same conjugacy class. We assume that a basis transformation has been done to make T y = 0 −1 1 0 . For odd

Gapless planon phase
Consider the following gapless model motivated by the constructions in Ref. [25] and [39]. Let m = ( 1−q 2 , 1+q 2 ) and n = (q, −q) where q is odd. We choose a slightly different stabilizer map which corresponds to flipping the sign of the top-right and bottom-right corner terms in the plaquette of the chiral plaquette model: Using the the results of Appendix B 1 one can confirm that the Gaussian fluctuations above the ground state manifold are gapless for this model. The gapless points occur at momenta (0, ± 2π 3 , ∓ 2π 3 ). Denote l = gcd(L y , L z ). The GSD is given by To understand the size-dependence of the GSD, we need to know the mobility of certain elementary excitations. First of all, a 1-kink 1 yz is a lineon with mobility along the diagonal line z = y + α with α ∈ Z. The string operator j e i q n·Φy+j,z+j translates the 1-kink along this path. However, a q-kink can move along both y and z directions (in steps of 6). The composite, 1 y−1,z (−1) yz 1 y,z−1 is a planon with mobility in the planar surface with normal vector (0, 1, −1). The planon is moved in the discrete direction by combinations of the lineon string operator and is moved along the wire by the string operator e i q L 0 m·∂xΦ . For those familiar, we can interpret that the plaquette term in the formalism of Ref. [25]: Θ y,z =φ L y,z − φ R y+1,z+1 + qθ y,z+1 + qθ y+1,z where φ L/R r = φ r ± qθ r . This observation along with the mobility of the excitations gives the heuristic picture of a stack of ν = 1 q Laughlin states lying in the planes defined by normal vector (0, 1, −1) which are then coupled by the qθ y,z+1 +qθ y+1,z terms.
If PBC are imposed in the discrete directions, which have length L y and L z , one can see there are gcd(L y , L z ) distinct paths of slope z = y. This is shown in Fig  3. Since each plane hosts Laughlin-like planons, naively this suggests q gcd(Ly,Lz) superselection sectors of planons living on the 2D surfaces partitioning the three torus, therefore explaining the q l factor in GSD. Let α = 1, 2, . . . , gcd(L y , L z ) label these distinct diagonal planes.
With the picture of layers of Laughlin states in mind, we write down a naive basis of logical operators for the system. Fixing a diagonal plane labeled by α, we define where yz ∈ α for X α . The operator X α cycles the planon 1 y−1,z (−1) yz 1 y,z−1 around the wire direction while Z α moves the planon around a cycle in the discrete direction.
However one can check this pairing does not produce a diagonal commutation matrix. Following the procedure in Appendix C one may consider the commutation matrix of Computing the Smith normal form of this matrix allows one to construct a set of canonical logical operators such that a d a = |GSD|, where GSD is given in Eq. (48). This procedure is necessary because while the model may superficially resemble stacks of Laughlin states, the inter-planar couplings can induce relations amongst logical operators. As an example, consider the case when gcd(L y , L z ) = 6k. One may verify, using Eq. (49) the following relations between the operators Z α : This lack of linear independence only holds at system size gcd(L y , L z ) = 6k and is reflected in the value of the GSD which is q 6k−2 .

IV. CSS MODELS
We study an example of the CSS model, given by the following polynomials: f (y, z) = y + z + yz, g(y, z) = n + y + z.
For general CSS models, the Gaussian spectrum is found to be In this case, we find that the spectrum is fully gapped for any value of n. If n > 4, it is further shown in Appendix E that all excitations are lineons moving along the x direction (i.e. the direction of the wires). In other words, since the excitations can move only along x, it is a "type-II" model in the yz plane. This is somewhat similar to Yoshida's Sierpinski spin liquid [4], which is a Z 2 stabilizer model with only lineons. To prove this result we generalize the cleaning argument for Pauli stabilizer models [3,40] to the present case. The details of the proof can be found in Appendix E. Notice that unlike the proof in Ref. [3] showing that there are no string operators at all in Haah's cubic code, here the model actually has "string operators" in the yz plane. This string operator, if cut into a segment, however, does not create a charge and its inverse. In fact, if one fixes the charge at one end of the string operator, the magnitude of the charge on the other end grows exponentially with the separation between them, costing an exponentially large energy to create the configuration. Therefore, the charges are still immobile. In Appendix D we compute a charge basis of this model, and the result is given by a + by, where a, b ∈ Z. Physically it means that there are infinitely many types of excitations at the (arbitrarily chosen) origin (0, 0) or (1, 0), labeled by two integers, and any other excitation can be transformed to an excitation at the two sites by applications of local operators.

V. DISCUSSIONS
In this work we have uncovered new classes of 3D fracton models through coupled wire constructions. When gapped, they are found to be lineon models in general, exhibiting infinite fusion structures and some with highly unusual bulk-surface correspondence. All these features distinguish them from previously known fracton phases, which come in two varieties: either they arise in Pauli stabilizer models, or can be constructed from condensation transitions in stacks of two-dimensional topologically ordered phases [9-11, 13, 16]. The latter construction has been generalized to a "topological defect network" picture of fracton topological order [21][22][23]. Common to all these existing constructions is that quasiparticle excitations have finite-order fusion. In contrast, the excitations in our model naturally have a Z fusion structure, typically associated with gapless fractonic U(1) gauge theories [41,42]. Technically, the difference can be attributed to the use of continuous bosonic fields in our microscopic model, as compared to the other constructions typically starting from finite-dimensional local Hilbert space. However, we believe the same physics can be realized starting from e.g. spin chains whose low-energy theory are Luttinger liquid. A related class of planon phases was recently studied in Ref. [37], whose surface states are similar to the models in this work, but with clearer layered structures.
In addition, we conjecture that at least many of the coupled wire models studied in this paper have relatively simple mobility structure, namely every excitation is a lineon. Examples of pure lineon phase were constructed in Pauli stabilizer CSS models (e.g. Yoshida's Sierpinski spin liquid), and our coupled wire models provide a natural framework for translation-invariant lineon phases, since the wire direction is a continuum field theory to begin with, and nontrivial mobility structure in the transverse direction can be encoded in the interactions between wires. This is similar to the field-theoretical construction in Ref. [37], where excitations coupled to (2+1)d Chern-Simons gauge fields are naturally planons. It will be interesting to understand whether existing lineon models can be incorporated into the coupled wire construction, for example whether the Siepinski spin liquid can be realized with the "CSS"-type construction discussed in Sec. IV. More broadly speaking, together with results in Ref. [37], our work suggests that it may be fruitful to classify fracton phases based on mobility structure, e.g. those with only planons or lineons, and we hope to explore these questions in future works.
In this work we have focused on the fusion and mobility structures of excitations, but have not explored any other statistical processes, such as the generalization of exchange process for lineons [42,43]. A related issue is that the topologically protected ground state degeneracy is generally spanned by rigid string and membrane operators, where the string operator moves excitations along the line. On the other hand, it is not entirely clear what the physical interpretation of the membrane operator should be. The nontrivial commutation algebra between the string and membrane operators should be related to certain statistical phases.
Many of our models have gapless modes, and currently we do not have a clear physical understanding of the nature of these gapless excitations. One possibility is that they can be interpreted as photons of certain U(1) gauge fields [44][45][46][47][48][49]. It will be of great interest to identify an effective field theory for these phases, perhaps along the lines of Ref. [50]. It is also important to understand the interactions between the gapped quasiparticles and the gapless modes, e.g. whether the gapless modes mediate long-range interactions between gapped excitations. A related question is the stability of the gapless phase against perturbations.
Another possible direction for generalization is to consider more complicated, interacting conformal field theories, such as Wess-Zumino-Witten theories with higher levels, replacing the Luttinger liquid (essentially free bosons) in each wire [25]. This may lead to interesting non-Abelian lineon phases.

ACKNOWLEDGMENTS
MC is grateful to Xie Chen for inspiring discussions and collaboration on a related project, and Kevin Slagle for collaboration at the initial stage of this work. AD thanks Dominic J Williamson for related discussions. JS would like to acknowledge discussions with Thomas Iadecola and Dominic J Williamson on related work which proved helpful. JS thanks Chris Harshaw for patient discussions about some very old results in linear algebra. M. C. is supported by NSF CAREER (DMR-1846109) and the Alfred P. Sloan foundation, and thanks Aspen Center of Physics for hospitality and support under the NSF grant PHY-1607611, where the work was initiated.

Appendix A: Coupled wire construction in 2D
Here we consider the coupled wire construction in 2D and show that under fairly general assumptions the model is an Abelian topological phase. While this result is certainly expected and a well-known folklore, it has not been explicitly shown in literature.
The wires are labeled by a single index j. Each wire is a Luttinger liquid described by a K matrix K w . We define l 1 · l 2 = l T 1 K w l 2 . For the gapping term, without loss of generality we only include interactions between nearest-neighboring wires: We make the following assumptions about P, Q's: 1. They should satisfy the null conditions So that the gapping terms commute with each other.
2. All bosonic fields are gapped out when the system is closed. Since the number of gapping terms is the same as the fields, as long as the gapping terms are linearly independent the condition is satisfied.
3. Topological order condition: if a local vertex operator creates no excitations, it must be a linear superposition of the "stabilizers" (with integer coefficients). For example, on each site, if l ∈ Z 2N satisfies l · P α = l · Q α = 0 for all α, then l = 0. Therefore, viewed as vectors over R 2N , {P α , Q α } span a complete basis. Moreover, the subspace spanned by {P α } and that of {Q α } are orthogonal. Two useful corollaries follow: 1) K αβ = P α · P β is an invertible matrix. 2) if l · P α = 0 for all α, then l is a linear superposition of Q α 's (over Z), and vice versa. As a special but important case of the topological order condition, there should be no local degeneracy. In other words, there exist no integers m α such that α m α Θ α j,j+1 is a non-primitive vector. This leads to the following condition: let M denote the following N × 4N matrix then the Smith normal form of M must have all non-zero diagonals being ±1.
In fact one should allow superposition of Θ α j,j+1 's from a finite cluster of wires. Now we classify the superselection sectors of kink excitations, which give anyon types of the topological phase. They are defined as the equivalence classes of localized excitations, up to local ones. In the coupled wire model, first consider kinks of Θ α at the j, j + 1 bond. They can be labeled by a vector e = (e 1 , e 2 , . . . , e N ) ∈ Z N with a e α -kink in Θ α j,j+1 . Next we classify which kinks can be locally created. It is not difficult to show from the topological order condition that it is sufficient to consider a two-wire local operator e i(l T 1 Φj +l T 2 Φj+1) . In order for the operator to only create excitations on the bond j, j + 1, one must have From our non-degeneracy assumption, we see that Therefore, the equivalence class is given by Z N mod out vectors generated by row (or column) vectors of K. Formally this agrees with the superselection sectors of an Abelian Chern-Simons theory with the K matrix K.
We also need to understand how kinks on different bonds are related. Suppose there is a kink e (j−1) on j − 1, j bond. To locally transform it into a kink on j, j + 1 bond, apply a vertex operator e il T Φj at site j, where l must satisfy. (A5) Let Q denote the N × 2N matrix formed by Q α 's. Eq. (A5) is solvable for any e if and only if the Smith normal form of Q has only ±1 entries. It is not clear whether this follows from the conditions imposed on P and Q, but we do not know of any counterexamples. Assuming this is the case, then e il T Φ annihilates the kinks on j − 1, j bond and create new kinks on j, j + 1 bond, given by e = P · l. The superselection sector [e ] may be different from [e], but since there are only a finite number of them, after sufficiently many steps the kinks can be transported without changing its charge type. Now we consider moving excitations along the wire direction. Consider an excitation e on bond j, j +1, and an operator W j (x) = e i α wαP T α Φj (x) , where w α are rational numbers. W j commutes with the gapping terms at the j − 1, j bond and creates excitations at the j, j + 1 bond, in particular w α P α · P β for Θ β j,j+1 . Then if we choose w = K −1 e, W j defines a string operator to move e along the wire: We have essentially described the anyon string operators, and can compute their braiding statistics. However, the string operator that moves an anyon across wires does not have explicit form, so we do not have general expressions for the braiding statistics.
We also need to consider the spectrum of Gaussian fluctuations. While we do not have closed-form expressions for the general case, we expect that the Gaussian spectrum should be gapped when all the conditions on P and Q are satisfied and the K matrix is invertible.

Example with N = 1
We consider fermionic systems with N = 1. We take P = (p 1 , p 2 ) and Q = (p 2 , p 1 ), where gcd(p 1 , p 2 ) = 1. It is easy to check that all our conditions are satisfied. Local excitations are of the form ±(p 2 1 − p 2 2 ), so the group is just Z |p 2 1 −p 2 2 | . A period-1 string is given by l = (1, 1) T , which generates an excitation p 1 −p 2 . If p 1 −p 2 = ±1, to get a "unit" excitation one needs to consider l 1 = (x, y) T , l 2 = (y, x) T , where xp 1 − yp 2 = 1 (always solvable as gcd(p 1 , p 2 ) = 1). It implies that translation along y can act nontrivially on anyons: under translation T y , a kink of strength p 1 x−p 2 y becomes p 1 y−p 2 x. Notice that T 2 y = 1. As an example, if p 1 = 5, p 2 = 2, the anyons form a Z 8 group and T y takes a ∈ Z 8 to a 5 . Ref. [24] considered p 1 = m+1 2 , p 2 = m−1 2 to obtain K = (m). In this case, T y does not act. This kind of Laughlin states enriched nontrivially by lattice translation was also studied in Ref. [51]. Now consider the system has an edge at j = 0. It is easy to check that the only local vertex operator is e iQ T Φ0 . The edge theory is thus a chiral Luttinger liquid, with 1 × 1 K matrix: K = (p 2 2 − p 2 1 ).

Appendix B: Gaussian spectrum
In the U → ∞ limit of coupled wire models, the gapping terms pin the fields Θ = 0. Here we study small oscillations of Θ around the minima by expanding cos Θ ∼ −1 + Θ 2 2 and solve the resulting quadratic theory. Below we introduce a mode expansion for the various bosonic fields involved and review how one finds single-particle spectrum of a quadratic Hamiltonian of bosonic creation and annihilation operators. We work out the spectrum of the chiral plaquette model as an example.

Mode expansion of bosonic fields
Mean field theory gives the following translationally invariant effective Hamiltonian (B1) where the index q allows for more than one Luttinger liquid per wire. We use the following mode expansion where the index k = (k x , k) and N w is the number of wires. Canonical commutation relations are imposed on a and a † 's: [a k,q , a † l,q ] = δ kl δ qq , [a k,q , a l,q ] = [a † k,q , a † l,q ] = 0. Consider the Fourier representation of the kinetic part of H eff : The term Θ 2 will involve terms of the form φ r+∆ where ∆ is some vector in the yz plane. The mode expansion for these terms is as follows Using these expressions above one can construct a BdG type Hamiltonian for the corresponding quadratic bosonic theory.

Bogoluibov transformation for bosons
We will be studying theories which are quadratic in bosonic creation/annihilation operators. Here we describe how to find the spectrum for a general quadratic Hamiltonian of bosons: where the "first-quantized" Hamiltonian h is defined as: Here T is Hermitian and U is symmetric. a i 's satisfy the canonical commutation relations [a i , a † j ] = δ ij . We perform a canonical transformation to a new set of annihilation operators b in which the Hamiltonian is diagonalized: Here Λ is the diagonal matrix of single-particle energy eigenvalues. Note the requirement that b i satisfy the canonical commutation relations for bosons means that the Bogoluibov transformation W is symplectic: So Λ does not simply correspond to the eigenvalues of the "first-quantized" Hamiltonian matrix h. However, using the fact that JW † J = W −1 , we can rewrite the diagonalization equation as a more standard eigenvalue problem: So we can solve for the spectrum by diagonalizing the matrix

Spectrum of the chiral plaquette models
Here, as an example, we calculate the spectrum of the chiral plaquette models. These models were written in the chiral basis, where m · Φ = m 1 φ L + m 2 φ R . We work, because of the simplicity of the mode expansion, in the (φ, θ) basis with m · Φ = aφ + bθ, where a = m 1 + m 2 and b = m 1 − m 2 is clear. Similarly c = n 1 + n 2 , d = n 1 − n 2 .
Define the following functions of k f φ = a 2 + c 2 + a 2 cos(k y + k z ) + c 2 cos(k y − k z ) + 2ac(cos k z + cos k y ), Schematically, Θ 2 term involves terms of the form φφ, θθ and φθ + θφ. Using the results of Appendix B 1 one can check that (B12) The single particle hamiltonian h k then has the follow-ing form Diagonalizing the matrix gives the spectrum So to determine if the fluctuations are gapped one needs to check that min( Thus one just needs to find the zero locus of |f k |, given by the following equations: (m 1 + m 2 ) cos α + (n 1 + n 2 ) cos β = 0, (m 1 − m 2 ) sin α + (n 2 − n 1 ) sin β = 0.
Let us define Assume for now u = 0. We can easily find So for both expressions to be positive-definite, we must have which implies that either t ≤ u ≤ s or s ≤ u ≤ t. It is easy to see that the s = t case is included. Therefore, if (t − u)(s − u) > 0, there are no zeros for |f k | 2 , which implies that it must have a positive minimum. One can further check that this condition also covers the u = 0 case.

Appendix C: Ground state degeneracy on torus
When the model is fully gapped, an interesting quantity to consider is the ground state degeneracy (GSD) with periodic boundary conditions imposed. The GSD can be computed using a method introduced by Ganeshan and Levin [52]. In their approach, all fields are treated as real-valued, with the Hamiltonian still given by Eq. (4). Compactness is then imposed dynamically by adding Collecting all the pinned fields C = {Θ α r (x), Q a r }, we compute their commutation matrix Z. Notice that Θ α r commute with each other, so do the Q a r 's, thus the nonzero commutators only occur between Θ and Q, and the commutation matrix takes an off-diagonal form: We then find the Smith normal form of Z 1 : where A and B are unimodular integer matrices. Then define and one obtains We assume that diagonal elements of D are ordered such that the first I of them, d 1 , d 2 , · · · , d I are non-zero. Then the GSD is given by |d 1 d 2 · · · d I |. The matrix V in fact gives the logical operators that span the ground state space [52]. More precisely, the commutation matrix is "diagonalized" in the new basis This form of C suggests that the logical operators come in two conjugate groups, one being AQ (with additional 1/d i factors that we haven't included yet), physically string operators along x, the other being B T Θ, which can be generally interpreted as surface operators in the transverse directions.

Appendix D: Algorithm to find a charge basis
We first define the charge basis in terms of the excitation map discussed in the main text.
Definition D.1. Charge basis Any local operator is said to create a trivial charge configuration. In other words, any charge cluster that belongs to Im where is the excitation map, is trivial. We now denote the set of all excitations by E, also referred to as the excitation module. We use Theorem 1 of Ref. [35] which states that the equivalence class of excitations modulo trivial ones is a torsion element of the cokernel of the excitation map. In other words, any topologically nontrivial local charge is an element of T coker = T (E/Im ). Torsion submodule T (M ) of a module M is defined as T (M ) = {m ∈ M |∃r ∈ R\{0} such that rm = 0}.
In order to calculate the charge basis given by T coker = T (E/Im ), we first note that we consider the excitation map represented by a matrix with matrix elements belonging to a polynomial ring R[x, y, z] over the ring of integers Z i.e. each element is polynomials in variables y and z with coefficients of monomials in Z. We can always bring the excitation map to this form i.e. with non-negative exponents of translation variables since we can choose any translate of the stabilizer generators as our generating set to write down a polynomial representation of the excitation map. The same holds for the charge basis. Even though an arbitrary charge configuration is expressed as a Laurent polynomial, if it is finite, we can change our choice of origin to write it as a polynomial with non-negative exponents of translation variables i.e. over a polynomial ring. We will use this idea to compute the charge basis using the trivial charge polynomials expressed in the non-negative cone i.e. with non-negative exponents of translation variables. Any non-trivial charge configuration can be expressed using the elements of this charge basis up to a translation.
We now introduce some definitions and concepts needed in the calculation of the charge basis. These definitions are taken from Ref. [53]. .., s}, set monomials X J = LCM(X j |j ∈ J) where X j are monomials. We say that J is saturated with respect to X 1 , ..., X s provided that for all j ∈ {1, ..., s}, if X j divides X J , then j ∈ J. In other words, it is saturated if all the monomial from X 1 to X s divide the LCM of the smaller subset defined by J. For example, consider the set (X 1 = xy, X 2 = x 2 , X 3 = y, X 4 = x 4 ) and choose the subset (X 1 = xy, X 2 = x 2 ). The LCM of elements in the subset is x 2 y which is divisible by X 3 but X 3 / ∈ J and hence the subset (X 1 , X 2 ) is not saturated.
For each saturated subset J ⊆ {1, . . . , t}, we let I J denote the ideal of R generated by {lc(g i )|i ∈ J} where lc denotes the leading coefficient. C J be the complete set of coset representatives for R/I J . Assume that O ∈ C J and also for each power product X, let J X = {lm(g i )|lm(g i ) divides X} where lm(g i ) denotes the leading monomial of g i . Definition D.4. Totally reduced polynomial A polynomial r ∈ A is totally reduced provided that for every power product X, if cX is the corresponding term of r, then c ∈ C J X . For a given polynomial r ∈ A, a normal form for f provided that f ≡ r (modI) and r is totally reduced. Now we state the main theorem (Theorem 4.3.3. of Ref. [53]) that describes the result for the coset representatives of the quotient R/I J Theorem D.5. Let G be a Groebner basis for the nonzero ideal I of A. Assume that for each saturated subset J ⊆ {1, ..., t}, a complete set of coset representatives C J for the ideal I J is chosen. Then, every f ∈ A has a unique normal form. The normal form can be computed effectively provided linear equations are solvable in R and R has effective coset representatives.
The actual calculation is best understood through examples. We now show an example from different classes of models mentioned in the main text.

Charge basis for different models
• We first consider the CSS model.
where n is an integer. Since there is a duality between the φ and θ sectors, we can consider only one sector, let's say φ and calculate the charge basis in the φ sector. The excitation map (D1) implies that any excitation pattern that belongs to the im is a linear combination of the two polynomials as shown in the map i.e. it belongs to the ideal yz + y + z, n + y + z . The Groebner basis of the ideal with lexicographic ordering is given by {g 1 = y + z + n, g 2 = z 2 + nz + n}. The leading terms are then given by y and z 2 i.e. leading monomials y and z 2 with coefficients 1 and 1. Now we use the definition that for each power product X, We can also simply arrive at this result by writing down relations y = −n − z and yz = n from the relations y+z+n = 0 and y+z+yz = 0 in the ideal. Using these two relations, we get z 2 + nz + n = 0. Hence, an arbitrary polynomial in y and z can be expressed only in terms of monomials 1 and z since y and z 2 can be reduced to polynomials in 1 and z.
The choice of basis monomials is not unique. Notice that because our original ideal is symmetric in y and z, we can also use the relation y 2 + ny + n = 0 and express an arbitrary polynomial in terms of basis monomials 1 and y i.e. as {a + by|a, b ∈ Z}.
• We now consider a family of models described by = m 1 + n 2 y + n 1 z + m 2 yz m 2 + n 1 y + n 2 z + m 1 yz .
The Groebner basis of the ideal with lexicographic ordering is given by {g 1 = y + z + 2p + 1, g 2 = z 2 + (2p + 1)z + 1}. The leading terms are then given by y and z 2 i.e. leading monomials m 1 = y and m 2 = z 2 with coefficients 1 and 1. Then, we get the saturated subsets J 1 = ∅, J y µy = {m 1 }, J z = ∅, J z µz >1 = {m 2 }, J y µy z = {m 1 } and J y µy z µz >1 = {m 1 , m 2 } where µ y and µ z are non-zero integer exponents of y and z. Thus, the corresponding ideals I J are I J1 = 0, I Jz = 0, I J y µy = 1 , I J z µz ≥2 = 1 and I J y µy z µz = 1 . We get C J1 = C Jz = Z while all other coset representatives are 0. Thus, a complete set of coset representatives for Z[y, z]/I is the set {a + bz|a, b ∈ Z}. 3.
where p is an integer. The excitation map (D5) implies that the trivial charge configuration ideal is given by p + y, 1 + py . The Groebner basis of the ideal with lexicographic ordering is given by {g 1 = y+p, g 2 = p 2 −1}. The leading terms are then given by y and p 2 − 1 i.e. the leading monomials m 1 = y and m 2 = 1 with coefficients 1 and p 2 − 1. Then we get the saturated subsets J 1 = {m 2 }, J y µy = {m 1 , m 2 }, J z µz = {m 2 } and J y µy z µz = {m 1 , m 2 } where µ y and µ z are positive integer exponents of y and z. Thus, the corresponding ideals I J are I J1 = p 2 − 1 , I J y µy = 1 , I J z µz = p 2 − 1 and I J y µy z µz = 1 . We get C J1 = Z p 2 −1 and C J z µz = Z p 2 −1 while the other coset representatives are 0. Thus, a complete set of coset representatives for Z[y, z]/I is the set {a + bz µz |a, b ∈ Z p 2 −1 }.
• We now consider another family of models given by = m 1 − n 2 y + n 1 z − m 2 yz m 2 − n 1 y + n 2 z − m 1 yz .
Then, we get the saturated subsets J 1 = ∅, J y µy = {m 2 }, J z = ∅, J z µz >1 = {m 3 }, J y µy z = {m 1 , m 2 } and J y µy z µz >1 = {m 1 , m 2 , m 3 } where µ y and µ z are non-zero positive integer exponents of y and z. Thus, the corresponding ideals I J are I J1 = I Jz = 0, I J y µy = 3 , I J z µz ≥2 = 3 and I J y µy z µz = 1 . We get C J1 = C Jz = Z, C J y µy = Z 3 and C J z µz = Z 3 while all other coset representatives are 0. Thus, a complete set of coset representatives for Z[y, z]/I is the set {a + by µy + cz + dz µz>1 |a, c ∈ Z, b, d ∈ Z 3 }. Recall that the stabilizer map is given by Formally the stabilizers can be written as Here for brevity and in analogy with stabilizer codes we denote XI ≡ e iφ1 , IX ≡ e iφ2 , ZI ≡ e iθ1 , IZ ≡ e iθ2 , suppressing the x coordinate dependence. We consider cleaning of arbitrary pair creation operators to show that there is nontrivial logical string operator. Since the code is "CSS", cleaning the pair creation operators of one type would be enough. Thus, we consider Z pair creation operators.

Cleaning to a minimal box containing the excitation patches
We can clean an arbitrary pair creation operator that creates a pair of excitation patches to a minimal box that contains the two patches. This can be done by using the commutation constraints due to the corners shared with X stabilizers, i.e. where the independent vertices XI and XX of the X stabilizer operator hits the pair creation operator enclosing the excitation patches. There are two orthogonal edges with these type of independent vertices, XX and XI, in the X stabilizer. Such edges are called good edges for cleaning [3] and having two of them here implies that one can clean the Z pair creation operator down to a minimal box containing the excitation patches as shown in Figs. 4, 5 and 6 using commutation constraints with the corners of the kind XX and XI of the X stabilizer.  The horizontal and vertical strips can be reduced again to the lines. For the horizontal line, we can show the deformation result for recursion as follows. Suppose the operator on site i is (a i , b i ). Then they must satisfy It follows that a i+2 = nb i , and The characteristic polynomial is x 2 + nx + n = 0, with roots ω 1,2 = −n± √ n(n−4) 2 . So we can generally write It is easy to see that if n > 4, both ω 1,2 are real and |ω 1,2 | > 1, so a i or b i grows exponentially large with i.

FIG. 6. Cleaning of pair creation operators
Similarly, for the vertical line, we have the recursion relations It follows that a i = nb i+2 , and The characteristic polynomial is nx 2 + nx + 1 = 0, with roots λ 1,2 = − 1 2 ± 1 4 − 1 n . So we can generally write Then When n ≥ 4, both roots |λ 1,2 | < 1. As a result, a i and b i decays exponentially and the string can not extend to arbitrarily long length. We have shown that both horizontal and vertical string operators must create charges exponentially large in the length of the string at least on one end. Now we further prove explicitly that no string operators can create charges of opposite values, meaning nb i = −a 0 − b 0 and a i + b i = −a 0 . Consider the horizontal line. This implies n[u 1 (1 + ω 1 )ω i 1 + u 2 (1 + ω 2 )ω i 2 ] = u 1 (2 + ω 1 ) + u 2 (2 + ω 2 ), and They have no non-zero solution for any i. Similar holds for the vertical relations. We now consider string operators that could be formed from L shaped operators in Fig. 4(c) and Fig. 5(c). Using the cleaning done for the horizontal and vertical strips, we can reduce these operators to width 1 operators shown in Fig. 7. The lines do not join exactly at the corner in order to cancel out the excitations around it. The patches shown at the ends show the excitation strengths at those ends. We now show string operator formed from such joining L shaped width 1 operators cannot form a nontrivial logical operator. a. L shaped operators in Fig. 7(a) and (b) In Fig.7a, due to the commutation with stabilizer generators, the vertices O h0 and O v0 have constraints [O h0 , IX n ] = [O v0 , IX n ] = 0 which imply O h0 = (ZI) n h ≡ (n h , 0) and O v0 = (ZI) nv ≡ (n v , 0). In order to cancel the excitation shared by the two lines of the L shape as shown, we get n h + n v = 0. Along the horizontal line, we have the recursion constraints due to the commutation as follows which give the recursive equation This can be solved using the quadratic ny 2 + ny + 1 = 0 which has two roots λ 1 , λ 2 . Thus, we get for the other corner Using the constraint O h0 ≡ (a h 0 , b h 0 ) = (−n v , 0), we get u h 1 = −u h 2 = −nv n(λ2−λ1) . Similarly, for the vertical line, we have which gives the recursive equation with roots λ 1,2 of the same characteristic equation as before i.e. ny 2 + ny + 1 = 0. Hence, we again get Using O v0 ≡ (a v 0 , b v 0 ) = (n v , 0), we get u v 1 = −u v 2 = nv n(λ2−λ1) . Now, in order for the L shape to form a string operator, we need to cancel out the excitation created at the two ends. Hence, we require a v i = 0, a h i = 0 and b v i + b h i = 0. But We notice that both λ 1 and λ 2 are negative with |λ 2 | > |λ 1 | and thus a v i = 0 is not possible for n > 4. We can have the same L shape with different boundary operators at the corner for the horizontal and vertical segments, as shown in Fig.7b. In this case, in order to cancel the excitations in the plaquettes, we get O h0 = (nn v , 0) and O v0 = (n v , −n v ). The recursion relations are the same as the L shape in Fig.7a We again notice that both λ 1 and λ 2 are negative with |λ 2 | > |λ 1 |, and so it is impossible to have a v i = 0.
b. L shaped operators in Fig. 7(c) and (d) In Fig. 7c , XX] = 0 such that O h0 = (ZZ −1 ) n h and O v0 = (ZZ −1 ) nv . In order to cancel the common excitation, we require n h = nn v . Thus, O h = (ZZ −1 ) n h ≡ (nn v , −nn v ) and O v = (ZZ −1 ) n v ≡ (n v , −n v ). For the vertical line, we have the same constraints as (E18) and (E19), thus we get with roots λ 1,2 of the same characteristic equation as before i.e. nx 2 + nx + 1 = 0. Using O 2 ≡ (a v 0 , b v 0 ) = (n 2 , −n 2 ), we get u v 1 = n2λ2 λ1−λ2 and u v 2 = −n2λ1 λ1−λ2 . For the horizontal line, we get which leads to n(a h i + a h i+1 ) + a h i+2 = 0. This leads to where ω 1 , ω 2 are roots of the characteristic equation . We get The cancellation of excitations requires a v i + b v i = 0, a v i + nb h i = 0 and a h i + b h i = 0. From (E31), we see a h i + b h i = 0 is not possible for n > 4.
We can have the same L shape with different boundary operators at the corner for the horizontal and vertical segments, as shown in Fig.7d. In this case, in order to cancel the excitations in the plaquettes, we get the boundary which gives The condition a h i + b h i to cancel the excitation at the ends to form a string operator cannot be satisfied for n > 4.