Lattice Clifford fractons and their Chern-Simons-like theory

We use Dirac matrix representations of the Clifford algebra to build fracton models on the lattice and their effective Chern-Simons-like theory. As an example we build lattice fractons in odd $D$ spatial dimensions and their $(D+1)$ effective theory. The model possesses an anti-symmetric $K$ matrix resembling that of hierarchical quantum Hall states. The gauge charges are conserved in sub-dimensional manifolds which ensures the fractonic behavior. The construction extends to any lattice fracton model built from commuting projectors and with tensor products of spin-$1/2$ degrees of freedom at the sites.


I. INTRODUCTION
A major goal of condensed matter physics is to understand and to classify all possible phases of matter; another one is to uncover phases outside contemporary paradigms. While these two goals are evidently contradictory, together they move the field forward. An example of a new class of systems whose complete understanding is still in progress is that of what is now commonly referred to as fractons in general, or more precisely, systems with fracton excitations. These systems have peculiar properties, including ground state degeneracies that depend both on topology and geometry of lattice discretizations, and excitations with restricted mobility that, in turn, make the dynamical relaxation to the ground states slow [1][2][3][4][5].
Recent reviews of fractons can be found in [6] and [7]. Thus far, they are classified into two types: in Type I phases a single fracton excitation cannot move alone, but a pair can bind into mobile dipoles; in Type II phases, all excitations are immobile [3,5,8,9].
It is this inherent immobility of isolated excitations that lead to slow dynamical behavior [1, 10,11]. The same restricted mobility and slow dynamical relaxation of excitations might be useful for building quantum memories [12][13][14]. In addition, fractons possess connections to elasticity theory [15,16] and gravity [17].
Fracton phases were originally constructed in lattice models; while their peculiar properties might appear unnatural for continuum descriptions, the construction of effective field theories that capture their low-energy properties is possible, as shown by Slagle and Kim [18] in the X-cube model [5,19]. The construction of effective field theories enables much further progress [18,[20][21][22][23][24][25]. Some features of fracton excitations are captured by the field theories in simple ways. Restricted mobility, for example, is encoded in additional charge conservation laws along sub-dimensional manifolds, such as planes for 3-dimensional (3D) 1 models, besides the conservation of total charge in the whole volume. The conservation of charges in planes implies that a dipole in the perpendicular direction is conserved. Hence charge conservation in sub-manifolds is equivalent to the conservation of vector charges (dipoles), a feature of higher-rank gauge theories [26][27][28][29][30][31][32][33][34][35][36][37][38][39], which, in general, are gapless. Nevertheless, gapped fracton models can be obtained from higher-rank gauge theories via the Higgs mechanism [40,41]. Gapped 3D fractons can also be obtained by either stacking [34,[42][43][44][45][46][47] 1 When referring to the dimensionality of the spacetime in this work we will use the notation (D + 1), with D the number of spatial components and the +1 refers to the time direction.
You et al. present a different route to a fracton field theory that is not cast as a higherrank gauge model. They present a Chern-Simons-like action with vector gauge fields that contains the sub-manifold conservation laws, hence also conserving dipoles. Their theory is gapped, and it can be discretized to a lattice to arrive at the Chamon model of Ref. [1].
The connection to Chern-Simons-like theories is appealing in that one would hope they can be generalized to describe classes of gapped fractons, much like Chern-Simons theories can describe classes of quantum Hall states [51].
In this work, we construct families of Chern-Simons-like theories of gapped fractons.
These theories have multiple gauge charges, and are described by an anti-symmetric K matrix and associated charge vectors. We arrive at these theories starting from microscopic lattice models, where we place a number n of spin-1/2 degrees of freedom (or qubits) at the sites. Such starting point is rather generic, and encompasses models such as the Chamon and Haah codes. Instead of tensor products of Pauli operators, we use the Dirac representation of the Clifford algebra to describe the site degrees of freedom. We show that the Dirac representation with 2 n -dimensional matrices is a natural mathematical framework to build the lattice models, and makes the connection to the field theory, a bosonization of sorts, rather simple. In the lattice theory, the fracton nature of the models are simple consequences of the lattice connectivity and the Clifford algebra, for example the immobility of single defects. In the continuum theory, these properties translate into charge conservation laws in sub-manifolds.
For the sake of giving a concrete but yet general example of the construction of these Clifford fractons, we build fracton models in any odd D = 2n + 1 spatial dimensions. This example allows one to track more easily the use of the 2 n ×2 n anti-commuting Dirac matrices γ I with I = 1, 2, . . . , 2n + 1, where γ 1 γ 2 . . . γ 2n+1 = i n . The model of Ref.
[1] corresponds to the simplest case, with D = 3 and 2 × 2 representations of the Dirac matrices. The 2n + 1 Dirac matrices form a maximal set of anti-commuting operators, and no operator (any product of Dirac matrices) other than the identity commutes with less than two of the γ's; it is this algebraic property that impedes the propagation of single fracton excitations.
We encode the anti-commutation relations of the Dirac matrices in a 2n × 2n antisymmetric matrix K for constructing a model in the continuum. This "bosonization"-type scheme is a generalization of that in Ref. [18]. The generic bosonic formulation in terms of K matrices and charge vectors allows us to take a continuum limit, and arrive at a (D+1)-dimensional Chern-Simons-like action L = 2n a, b=1 where the differential D = 0, again a condition ensured by relations between the microscopic lattice charge vectors T (I, α) and the K matrix. The conservation of the n currents not only on the full volume, but also on sub-manifolds, follows from a linear dependence of the charge vectors, D I=1 T (I, α) = 0, which in turn affects the differential operators D (α) a . These fracton models have a ground state degeneracy for systems of linear size L (hypervolume L D ). In the case of the "integer" fractons we constructed on the lattice, the K-matrix Pfaffian equals 1, and the degeneracy is 2 (D−1) 2 D−3 L .
The paper is organized as follows. In Sec. II, we construct, as example, a microscopic Clifford fracton model in D = 2n + 1 spatial dimensions. In Sec. III, we construct the corresponding effective theory in the continuum. In Sec. IV we examine several properties of the effective field theories. We close in Sec.V, with a brief summary and final remarks.
Details of several computations as well as additional relevant discussions are presented in the appendices.

II. MICROSCOPIC MODEL IN ARBITRARY ODD DIMENSIONS
The 3D model of Ref.
[1] uses the simplest representation of the Clifford algebra, where the Dirac γ matrices are 2 × 2: γ I = σ I , I = 1, 2, 3, with the σ I the Pauli matrices. The corresponding Hilbert space of the local degrees of freedom is 2-dimensional.
The local Hilbert space is 4-dimensional in this case. (This representation is obtained from a tensor product of two sets of Pauli matrices.) In D = 2n + 1 dimensions, we work with 2 n × 2 n representations of the Clifford algebra, i.e. the Dirac matrices γ I , I = 1, . . . , D, (We build these matrices explicitly in appendix A.) The construction of the fractons in odd-dimensional D = 2n + 1 space proceeds as follows. We start with an face-centered hypercubic lattice, that can be thought as the even sublattice Λ e of a hypercubic lattice with orthogonal basis vectorsâ I , I = 1, . . . , D.
We place the degrees of freedom on this even sublattice, as well as operators Γ (I,α) with α = 1, . . . , n acting on these degrees of freedom. The operators Γ (I,α) are built as products of the γ-matrices (in turn built from tensor products of Pauli matrices, see appendix A).
We take Γ (I,1) ≡ γ I , which we call principal configuration. The need for the additional Γ (I,α) with α = 2, . . . , n comes because the local Hilbert is 2 n dimensional, and consequently n operators are necessary to gap the theory.
We define O (α) operators centered on the odd sublattice Λ o , Notice that O (α) x 2 = 1 1 follows because the Γ (I,α) are products of Dirac matrices. Using these operators we construct the Hamiltonian where all coupling constants g α are chosen to be positive. We can further choose the oper- In this case, the Hamiltonian is a sum of commuting projectors and there are as many commuting projectors (up to constraints that we shall see in a moment give a topological degeneracy) as the number of degrees of freedom in the problem. The connection between the choice of operators Γ (I,α) and the commutations between the O (α) x stems from the geometry imprinted via the definition Eq. (7), and is depicted in Fig. 1. The O (α) 's trivially commute both when they are defined at the same site x ∈ Λ o or when they do not share any sites (in Λ e ); there just remains two cases to be checked: when they share one and two sites. The neighboring O (α) 's, defined at sites x and x + 2â I of Λ o , share the Λ e site at x +â I , and they commute if Neighboring O (α) 's, defined at sites x and x +â I +â J of Λ o , share the two Λ e sites at x +â I and x +â J . The operators on those sites either commute or anti-commute, which can be cast as with η (αβ) IJ = 0 or 1. The desired commutation relations Eq. (9) are guaranteed if In particular, the commutation (10) implies η (α,β) II = 0.
All these conditions can be satisfied using Dirac matrix representations of the Clifford algebra. The simplest example is the D = 3 contained in Ref.
[1], where one uses the 2 × 2 representation Γ (I,α) I = 1 2 3 In D = 5 we use the 4-dimensional representation of the Dirac matrices and take the following Γ (I,α) operators: that satisfy Eqs. (10) and (12) and hence yield the set of commuting projectors O (1) x and O (2) x , x ∈ Λ o . As we shall see in the next section, the field theory formulation provides a systematic way to construct operators of the sets α ≥ 2 in arbitrary odd dimensions and satisfying all the required commutation rules.
The ground state of these models (for any D = 2n + 1) correspond to all O (α) x having eigenvalue +1. Excitations or defects correspond to those operators having instead eigenvalue −1. That these models are fractonic requires that there is not a single local operator whose net effect is to generate a defect pair or equivalently move a single isolated defect. In x = 1 1 , k = 1, . . . , 2 D−1 , α = 1, . . . , (D − 1)/2 (15) where k accounts for all the 2 D−1 sub-lattices Λ o,k . Eq.(15) already give us a hint of the degeneracy of the model, which will be at least 2 [(D−1) 2 D−2 ] , but it happens that this number is a lower bound to the degeneracy and in fact, the degeneracy can also depend on the system size, as was shown in [2]. For a hypercube of volume L D , the degeneracy dependence on the linear size L is 2 (D−1) 2 D−3 L (see appendix C for details).

III. LOW-ENERGY EFFECTIVE FIELD THEORY
In this section we shall derive an effective field theory capturing the low-energy physical properties of the lattice models given by (8). To connect the operators that act on the microscopic degrees of freedom with the suitable operators that possess a well-defined continuum limit, we define a map parametrized by the vectors T (I,α) a of Eq. (5): where the repeated matrix indices a, b = 1, . . . , 2n are summed. We need only D − 1 = 2n independent fields A a to construct all the required operators. In addition, we see that this parametrization introduces a symmetry with Q an arbitrary matrix.
For the case of the principal configuration, where Γ (I,1) ≡ γ I and, accordingly, T Properties of the fields A a and of the matrix K can be obtained from the analysis of the principal configuration. Indeed, we start by computing where we have used the BCH formula and assumed that the commutator appearing in this expression is a c-number. Since γ I x and γ J x ′ must commute if x = x ′ and anti-commute if x = x ′ and I = J, we impose which will be interpreted as an equal-time commutation relation in a field theory formulation with a canonical pair A a ( x) and Π a ( x) = 1 π (K T ) ab A b . We shall return to this point later. Using the commutation relation (20) in (19) leads to the further following conditions to match the anti-commutation relations among the γ I x : The condition (21) is particular to the principal configuration, and reflects that the building blocks of the theory are anti-commuting objects (Dirac matrices); it is not needed for the other T -vectors with α ≥ 2, since the associated operators (products of Dirac matrices) may either commute or anti-commute. Notice that the conditions in (21) imply that the fields are compact, since the shifts do not change γ I x because exp 2iπ The most general solution (independent of t) for the condition in the second line of (21) corresponds to the case where K is an anti-symmetric matrix. In writing the commutation relation (20) we assumed that the inverse K −1 exists, which requires det K = 0, a condition only possible to satisfy for even-dimensional anti-symmetric K matrices. Recall that K is a 2n × 2n = D − 1 × D − 1 matrix, so the construction works for odd dimensions D.
A useful relation involving the t-vectors emerges when we consider the product of all matrices γ in the same site. Suppressing the matrix indices for simplicity, we have where we have used (21) in the last step. In order for the right-hand side to be proportional to the identity, we require which we refer to as the neutrality condition. Notice that this is satisfied with the choice in (6). Moreover, we shall require the neutrality condition for all sets of operators Γ (I,α) , which corresponds to We now proceed by analyzing the field theory counterpart of the lattice operators (7).
Using the representation (16), it follows that Given the commutation relation (20), we have to determine the conditions on T (J,α) and K so as to produce commuting operators O (α) x , i.e., so that the field theory representation on the right-hand side reproduces the commutations in Eq. (9). As discussed in the previous section, there are two situations where a nontrivial commutation rule may arise: when the operators share one or two sites. When they do not share any site, the commutation is trivially satisfied. To take into account these two situations, we just have to consider two operators O where Notice that C IJ . The symmetry in the IJ indices follows directly from the way C (αβ) IJ is defined, whereas the anti-symmetry in the αβ indices follows from the anti-symmetry of the matrix K. In particular, the condition (28) is automatically satisfied if α = β, which is consistent with the fact that O (α) x operators of the same kind commute with each other. A systematic procedure for constructing T -vectors satisfying the condition (28) is presented in appendix B. Next we consider the continuum limit of the relation (27). The expansion of the field A b reads As the unit vectorsâ J have the components a I J = δ I J , we get We see that the neutrality condition (26) ensures that the first term in the exponential vanishes, so that the operator O The Hamiltonian in (8) becomes with We see that the ground state corresponds to the case where all the cosines in (33) are We can enforce this in a corresponding field theory description of the ground state through a Lagrange multiplier, as we will discuss in a moment.
Before going to the field theory it is convenient to express the operator M (α) ( x) in a way that solves the constraint of the neutrality condition (26). Thus, we single out one of the directions, say the last one J = D, and write where the derivative operator D J is defined as It is also convenient to define another differential operator as In terms of D which makes evident its invariance under gauge transformations with ζ (α) = ζ (α) ( x, t) being a set of arbitrary functions of spacetime coordinates. In fact, notice that Therefore, the condition above, needed for gauge invariance, is precisely the condition for commutation of the cosine operators (28).
With all these elements in place, we can write down a field theory which describes the ground state of the microscopic fracton model, The first term is responsible for the commutation relation (20) 2 , whereas the second one enforces the ground state constraints, with A  (20), since for each pair of coordinates we always have two contributions because of the antisymmetry of the matrix K, for example, . This pair of terms must be brought into a single term through integration by parts before computing the canonical momentum.
transform as Thus, we end up with a bona fide gauge theory, which resembles the Chern-Simons description of topologically ordered systems. The gauge-invariant "electric" and "magnetic" fields can be defined as where ǫ a 1 a 2 ···a D−1 is the Levi-Civita tensor of rank D − 1.

A. Level Quantization
Now we will explore some properties of the effective field theory (40). Firstly, it is interesting to understand whether there is a notion of quantization of the "level" of the theory, which in the present case is given by the matrix K. To address this question we consider the principal configuration T (I,1) = t (I) . In this case, the t-vectors must satisfy the conditions in (21). Then, we use the symmetry transformations in (17) to make a specific choice for the t-vectors. For example, if we pick up the canonical form (6), we obtain the following level quantization condition: i.e., all the off-diagonal elements must be odd integers, and consequently nonvanishing. Of course, different representations of the t-vectors yield different quantization of the elements of the matrix K, but in all the cases we end up with some notion of quantization due to the conditions in (21).
From the field theory alone, the quantization of the level can be understood as follows.
Consider a manifold M = S 1 × M D , with S 1 representing the time direction with period [0, τ ) and M D a spatially closed manifold. Due to the compact nature of the fields A a , we have a quantized flux Consider large gauge transformations that wind around the time direction. The A (α) 0 field transform as and the corresponding variation of the action under these transformations is where Z ab an integer valued quantity obtained from the summation above. From this, it is straightforward to note that in order for the quantum theory to be invariant under the large gauge transformations (45), the K-matrix elements have to be integer valued. Further details of this calculation are found in appendix D

B. Three dimensional case
It is instructive to compare the effective field theory that we have obtained for the particular case of three spatial dimensions, D = 3, with the result of Ref. [20]. For D = 3 there is only one configuration, the principal configuration α = 1. The matrix K in this case is The action, in terms of electric and magnetic fields, reduces to With the canonical choice (6), the derivative operators entering the electric and magnetic , whereas the coefficient k must be an odd integer, in accordance with (43). In order to compare with [20], we just need to rename the fields and the derivative operators according to [20]), which leave both the action and the commutation relation in (48) unchanged. In this form, we can immediately compare with the results of [20] with the following identification between the parameters k = s/2. In that work, the original Chamon model (with full cubic symmetry) is recovered for s = 2 (in [20], the level quantization is s ∈ Z), which in our normalization corresponds to k = 1. This choice describes the 2-state system at each site, as expected. Also, this choice of k is allowed by the level quantization (43) associated with the canonical choice for the t-vectors. For a general discussion of how the K-matrix elements are determined by the microscopic theory, we refer the reader to appendix E.

C. Conservation Laws
A gauge-invariant coupling to matter can be introduced in the action (40) through the , provided that the current satisfies the continuity equation By integrating over the whole space and assuming periodic boundary conditions along all directions, it follows that charge is conserved in the whole system, In addition, given the form of the derivative operators D (α) a , we also have more restrictive conservation laws. These extra conservation laws require that charge is also conserved on a set of sub-manifolds of the system. It is due to these extra conservation laws that the fracton behavior of the excitations emerges.
To find the sub-manifolds where charge is conserved, we use the definition of the derivative operators D (α) a in (36) to write the continuity equation as where we have defined J it is convenient to introduced the directionsx σ I ID ≡â I + σ IâD , with σ I = ±1. In this notation, (51) can be written as This form of the continuity equation induces 2 D−1 extra conservation laws, explicitly, that charge must be conserved in each of the (D − 1)-dimensional sub-manifolds labelled by over any of these sub-dimensional manifolds, we obtain the following conserved charges These conservation laws, in turn, imply that the dipole moment in the direction perpendicular to those manifolds is conserved. Naturally, such conservation laws impose several restrictions on the mobility of the particles. We build in detail the form of the excitations in appendix F.

D. Ground State Degeneracy
Here we discuss the computation of the ground state degeneracy using the effective field theory. Naturally, in the continuum limit the degeneracy is infinite so that we shall adopt some kind of discretization (regularization) of the theory. The form of the conservation laws in the sub-dimensional manifolds provides a very natural way to discretize the theory in a layered structure. The basic idea is to consider the system as a stack of layers corresponding to the sub-dimensional manifolds where charge is conserved.
Let us start with the case D = 3. The action (40) reduces to where we keep explicitly only the part relevant for the computation of the degeneracy. In this case, charge is conserved in 2 2 sub-spaces labeled by σ 1 , σ 2 = ±, with the corresponding The strategy is to write the action (54) in terms of the coordinates x σ 1 13 , x σ 2 23 , x ⊥ , where x ⊥ is the coordinate perpendicular to the plane defined by the directions x σ 1 13 and x σ 2 23 . Upon this change of variables, where J is the Jacobian of the transformation. As this transformation is linear, J is just a constant and can be absorbed in dx ⊥ . The transformation from the coordinates x 1 , x 2 , x 3 to x σ 1 13 , x σ 2 23 , x ⊥ will change the limits of integration. However, as the ground state degeneracy in each plane with periodic boundary conditions (forming a torus T 2 ) does not depend on the area of the plane (torus), so we can ignore the area of integration in our computation as long as we assume periodic boundary conditions along the plane x σ 1 13 -x σ 2 23 . The next step it to discretize the coordinate x ⊥ . We consider that the perpendicular direction is composed by a stack of N layers, where 2a is the separation between the planes, twice the lattice spacing of the microscopic model. This discretization ties the number of layers to the linear size: N = L/2a. (Equivalently N = L/2 given we set a = 1). The gauge fields A a need to be rescaled properly The action (54) becomes Thus we end up with N copies of (2+1)-dimensional theories.
where l i is the size of each cycle of the 2-torus, and the arguments of the exponentials are properly dimensionless. These objects are gauge-invariant. In fact, under a gauge transformation, the fields transform as Let us analyse, say, the first holonomy in (60). Under a gauge transformation, it changes by a factor exp i The above condition is satisfied with the general periodic boundary condition Infinitesimal gauge transformations correspond to n i 1 = 0, whereas n i 1 = 0 are associated with large gauge transformations. Similarly, for the second holonomy in (60), we obtain A large gauge transformation satisfying all these conditions can be constructed explicitly, This implies an equivalence for the gauge fields Now we consider the ground state configuration, which corresponds to solutions depending only on the time, Plugging this equation into the action (59) we obtain The holonomies become From the action (67) it follows the commutation rule leading to the commutation relation between the holonomies which implies a 2k-fold degeneracy for each plane i. The degeneracy of the layered system is then Finally, taking into account that we have 4 sub-dimensional manifolds where charge is conserved, the total degeneracy is Using that N = L/2, we recover the degeneracy of the lattice model in D = 3.
In this case, charge is conserved in the following 2 4 sub-dimensional manifolds with the corresponding measures, We proceed similarly to the previous case, i.e., we write the action in terms of the coordinates of a sub-manifold where charge is conserved plus a perpendicular direction x ⊥ , which is then discretized. With this, the action (74) becomes where the fields A i a were rescaled as in (58). Now, the key point is that we can rotate the matrix K according to (17) to bring it to the block-diagonal form where k 1 and k 2 are real and positive. In this basis, the fields A i a decouple pairwise, Thus, we can construct the following pairs of holonomies the 4-dimensional torus T 4 in T 4 = T 2 × T 2 . Therefore, by proceeding in the same way as in the case D = 3, we see that these holonomies lead to a (2k 1 × 2k 2 )-fold degeneracy in each layer. For N layers, we get Finally, considering the 2 4 sub-dimensional manifolds, it follows that the total ground state degeneracy is which is expressed in a basis-independent way in terms of the Pfaffian of the original matrix K.
The generalization to the odd D-dimensional case is immediate. We decompose the space in a (D − 1)-dimensional sub-manifold corresponding to one of the 2 D−1 sub-spaces where charge is conserved, and a perpendicular dimension which is then discretized. Next, we make the transformation (17) to bring the matrix K to the block-diagonal form where all k's are real and positive. In this basis, the fields A a decouple pairwise, which is equivalent to decomposing the (D − 1)-dimensional torus as The corresponding degeneracy is Taking into account the N layers, we have Finally, considering all the 2 D−1 sub-dimensional manifolds, we obtain the total ground state degeneracy GSD = 2 For the case of Clifford fractons, where k 1 = k 2 = · · · = k D−1 2 = 1 or, equivalently, Pf(K) = 1, the ground state degeneracy reduces to where we have again used that N = L/2. This is precisely the result shown in the end of Sec. II obtained directly from the lattice model.

V. FINAL REMARKS
In this work we constructed fracton models on the lattice and identified their continuum description in terms of Chern-Simons-like theories. The construction is generic in that it applies to any system whose microscopic Hamiltonian is a sum of commuting projectors built from tensor products of spin-1/2 operators. Instead of working directly with tensor products of Pauli operators that represent the local variables, we utilize the Dirac representation of Clifford algebras. This representation makes a connection between the lattice model and the field theory simple. Our formalism can, in principle, be used to analyze other lattice models, such as those that exhibit subsystem symmetry protected topological (SSPT) phases [45,52] or type II fracton phases [3]. Applying this formalism to these problems is a natural direction for future work.
In the field theory, the algebraic structure of the Dirac matrices is encoded in an antisymmetric matrix K. The details about an specific lattice model enter via this matrix K (whose dimension depends on the size of the representation), the charge vectors T (that specify the operators that are placed on the sites), as well as the lattice vector positions of the sites themselves. Given these data, one can follow the prescription here presented and derive an effective field theory for any type of Clifford-like fracton, such as the 3D Chamon (with a 2×2 Dirac representation) or the 3D Haah (with a 4×4 Dirac representation) codes.
As a concrete example, we built fracton theories in odd D spatial dimensional spaces. We discussed the properties of the resulting Chern-Simons-like theory, such as their currents, which are conserved in sub-manifolds, and the topological degeneracy of the ground states, which formally depends on the Pfaffian of the matrix K and, as usual in fracton systems, on the linear size of the system.
Properties such as the mutual statistics of the quasiparticles where not explored in the present work. The restricted mobility of fractons makes it unnatural to speak of standard braiding. However, the authors in [53] were able to develop a theory of fusion and statistical processes that incorporates the mobility restrictions common in fracton models. An interesting question for future exploration is how our formalism could incorporate their notion of statistics.
For readers familiar with the K-matrices and charge vectors T appearing in the description of Abelian fractional quantum Hall states [51], as well as their quantum wire constructions [54][55][56][57], it is tempting to expect that the description here presented -for "integer" fractons given our K and T 's -could possibly lend itself to the analysis of fractional fractons.
This is an intriguing possibility that merits further investigation, but keeping the following points in mind.
The approach of this paper resembles quantum wire constructions of topological phases, but instead of wires we deploy (0 + 1)-dimensional degrees of freedom, i.e., ours is a "quantum dot" construction. Like in the wire constructions, we identify families of commuting operators that can be simultaneously pinned and gap the system. In the wire systems, fractionalization already takes place in the (1 + 1)-dimensional building blocks, and it is carried over to higher dimensions by coupling the wires, notably using only integer charge transfer operators. However, there is no fractionalization in the quantum dots of the construction of this paper. Of course, one may generalize the construction presented here to start with wires instead of dots, in which case fractionalization may appear more easily.
Added note : It has been brought to our attention that the word "fracton" has been used in physics in other contexts before. An early use was in [58] in reference to fractional charges in quantum chromodynamics. In our construction, we adopt the modern meaning of the word as stressed in the main text.

Appendix A: Euclidean Dirac matrix representations of Clifford algebras
We construct fracton models in odd D = 2n + 1 dimensions using representations of the Clifford algebra. Specifically, we use the Euclidean Dirac matrices. Below we construct these representations and show properties that these matrices satisfy. These properties are used, for example, to argue that there is no operator that can move defects in the corresponding fracton models.
Let us work with matrices defined as the tensor products of Pauli matrices: with µ i = 0, 1, 2, 3 and σ 0 ≡ 1 1. We shall obtain a set of 2n + 1 mutually anticommuting matrices for any n. We construct this set inductively.
For n = 1, the set contains the matrices γ These 2n + 1 matrices multiply to the identity up to a prefactor: where the ± simply depends on the order that the matrices are multiplied (the choice of order of the indices I ∈ S (n) ). This relation can also be proved by induction. Notice that it holds for n = 1. If it holds for n − 1, then it follows that This property means that the last, or (2n + 1)th, γ-matrix can be obtained from the product of all the other 2n matrices. It also follows that any matrix that is a tensor product of Pauli matrices can be written as products of these 2n γ-matrix. (Notice that there are 4 n possible tensor products of Pauli matrices, and 2 2n = 4 n choices of whether a γ-matrix enters or not the product of γ's.) The construction above yields a set of 2n + 1 matrices γ The set of indices I ∈ S (n) can be interchanged to I = 1, . . . , 2n + 1, which is the notation we use in the main text for the Euclidean Dirac matrices.

Properties of the Euclidean Dirac matrices
Let us now show three useful properties of the 2n + 1 matrices γ (n) I with I ∈ S (n) .
1. The identity is the only tensor product of Pauli matrices that commutes with all the Dirac matrices. In other words, only the matrix γ To show this property, suppose that there is a matrix γ (n) J that commutes with all the 2n + 1 matrices. This J must be of the form J = j0 for γ i3 , for i ∈ S (n−1) . But because S (n−1) is maximal, there is no new j / ∈ S (n−1) to add.
We thus conclude that the set of 2n + 1 matrices γ (n) I with I ∈ S (n) is maximal.
3. There is no matrix γ  I with I ∈ S (1) , because no one Pauli matrix commutes with two Pauli matrices. Now suppose that the statement is true up to n − 1; let us analyze the consequences for when we consider n.
Let us break the problem in four cases: In this case, the commutation with γ (n) 0...0 1 and γ (n) 0...0 2 comes for free. Therefore we reduce the problem to finding γ (n−1) j that commutes with 2(n − 1) matrices γ (n−1) i with i ∈ S (n−1) . Since there is no solution for this problem (the statement is true for the case with n − 1), then there is no solution for the case with n either.
This is the simplest case; γ i3 with i ∈ S (n−1) . This is equivalent to finding 2n − 1 matrices that anticommute with γ (n−1) j among the γ (n−1) i with i ∈ S (n−1) . This is impossible since the set S (n−1) is maximal (see above).

Appendix B: Fracton models build from the Clifford algebra representations
We can construct a fracton model in D dimensions if D is odd. In this case we take n = (D − 1)/2 and we can use the matrices γ (n) I , I ∈ S (n) in the construction. Here we shall label these 2n + 1 matrices simply γ I , I = 1, . . . , 2n + 1, as we did in the main text.
The construction can be made in the D-dimensional hypercube, with orthogonal basis vectorsâ I , I = 1, . . . , D, as presented in the main text. We place the degrees of freedom on the even sublattice Λ e . The dimension of the local Hilbert space at each site is 2 n , or equivalently, that associated with the n spins or gradings of Pauli operators used to construct the γ-matrix representations. At these even sublattice sites we place operators Γ (I,α) , α = 1, . . . , n built as products of the γ-matrices (in turn built from tensor products of Pauli matrices).
The first set of operators, with α = 1, is the set of Dirac matrices γ I , I = 1, . . . , 2n + 1, or explicitly The other sets are needed to gap the model.
operators centered at sites x on the odd sublattice Λ o , and using these the Hamiltonian We can choose the operators Γ (I,α) such that As stated in the main text, in this case i ) the Hamiltonian is a sum of commuting projectors, and ii ) there are as many commuting projectors as the number of degrees of freedom in the problem (up to constraints tied to the topological degeneracy).
Let us first focus on the operators O (1) for simplicity. These are defined as The operators γ I x satisfy the following commutation relations: (We remark that these models are bosonic, and not fermionic; the Dirac matrices represent the tensor product of local Pauli matrices, that in turn represent spin degrees of freedom on the lattice.) Given these commutation relations, it follows that all distinct O (1) x and O (1) x ′ that share common sites commute: 1) they either share a single site along the line that connects them, in which case the same operator (same I), or 2) they share two sites with different components I and J entering in each of O (1) x and O (1) x ′ , and hence there is a factor of −1 from the anti-commutation relation of each common site, and hence in total a factor (−1) 2 , leading to the commutation of the two operators. x , and as demonstrated above in Sec. A 1, there is no operator that anti-commutes with one and only one of the γ I . Therefore, it is not possible to construct a local operator whose sole effect is to create a pair of defects, or move a single defect. Defects are only created in at least quadruplets in any dimension D = 2n + 1, much as in the D = 3 model in Ref. [1]. This property that defects cannot be created in pairs, but only in at least quadruplets, underscores the fracton nature of these odd D models.
Let us now discuss the other operators O The commutation relations between the γ-matrices can be encoded in an integer-valued anti-symmetric K-matrix via where repeated index summation over the a and b are used. The correct commutation relations follow from requiring that It follows that the commutation relations or equivalently, using Eq. (B8), Then, condition Eq. (B9) is equivalent to where While we just need C (αβ) IJ to vanish mod 2, we can simply demand that it vanishes, and still solve the problem as we show below.
Let us now construct vectors T (I, α) a that satisfy C It follows, for I, J = 1, . . . , 2n, that or equivalently, that Hence the condition that the commutation relations C (αβ) vanish require that the sets of (α, β)-indexed matrices (L (α) K L (β) ⊤ ) be anti-symmetric (in the indices I and J). Let then where the A (αβ) are anti-symmetric matrices for any of the α, β pairs. For given choices of matrices A (αβ) , we can solve sequentially for i.e., start with β = 1 and L (1) = 1 1, obtain L (2) for some arbitrary choice of A (21) , then for some choice A (31) obtain L (3) , and so on. In other words, we can determine the L (α) from using β = 1 and L (1) = 1 1 in Eq. (B22): for anti-symmetric choices of A (α 1) . Notice that the L (α) cannot be equal, otherwise two sets of T (I,α) 's would be identical. The number of solutions (number of α ′ s) depend on the dimension D − 1 of the matrices, for example the matrix K. Notice that if K is 2 × 2, any anti-symmetric matrix is proportional to iσ 2 , and therefore it follows from Eq. (B23) that one cannot get a non-trivial solution other than L (1) ∝ 1 1.
The following anti-symmetric matrices A (α 1) commute with both K 2n and K −1 2n ; using this property and the anti-symmetry of both the A (α 1) and the K 2n , one can show that the A (αβ) obtained through Eq. (B25) are antisymmetric, as required.
With these A (α 1) , one can proceed to find the L (α) matrices and then the vectors T (I,α) , and finally the operators Γ (I,α) with the desired commutation relations.

Appendix C: Degeneracy
In the main text we argued that the topological degeneracy of the model is at least The degeneracy can be greater, and can depend on the system size. Here we follow Bravyi, Leemhuis, and Terhal's calculation in their appendix A of [2].
It follows from D I=1 T (I,α) a = 0 that the Γ (I,α) multiply to the identity (up to a factor of magnitude 1, that also depends on the order of multiplication). Using this property, we arrive at the equivalent of their parity checks: (C2) (The t above conforms to their notation; they are not related to our t-vectors.) Notice that, for each α, the D equations are linearly dependent, and that the sum over all of their left hand side is identically zero. Suming any pair of these equations yield The solutions of these equations for the case when L 1 = L 2 = · · · = L D in a similar way as in 3-dimensions: first use two lines (with 2L/2 = L sites on a sublattice) and generate the solution on a plane, then two planes and generate the solution in the 3rd dimension, and after that proceed accordingly, use 2 3-dimensional hyperplanes to generate the solutions in 4-dimensions, and so on. The number of logical quibts generated in this way is C D = L/2 × 2 × 2 × · · · × 2, with D − 1 2's, i.e., C D = 2 D−2 L. When we take into account all the α = 1, . . . , (D − 1)/2, we have (D − 1) 2 D−3 L logical qubits. Therefore, the ground state degeneracy is This implies The corresponding variation of the action is By using the flux quantization (D2), we obtain In order for the quantum theory to be invariant under large gauge transformations, we must By integrating over M D , we obtain By choosing again the configuration in (D10), which is an integer.

Appendix E: K-matrix and microscopic theory
The matrix K is determined in three steps: i) we determine its dimensionality; ii) we find constraints on the possible values of the elements; iii) we fix the elements in order to match the physical properties of the lattice model. The steps i) and ii) follow directly from the algebra of the operators in (16).
Step iii) is more subtle and follow from the algebra of the ground state holonomy gauge-invariant operators of the effective field theory.
Step (i): Dimensionality of K-matrix. In the class of models we have constructed in the manuscript, the dimensionality of the "spin" operator acting at each site is tied to the spatial dimensionality. Indeed, In this way, the algebra of the spin operators in D spatial dimensions can be written in terms of the Clifford algebra of Dirac matrices of dimensionality 2 n × 2 n . In this case, we have 2n Dirac matrices. Therefore, according to the representation of (16) of the paper, we need 2n distinct fields A to reproduce properly the algebra of the Dirac matrices. This fixes the dimensionality of the matrix K to be 2n × 2n.
Step (ii): Elements of the K-matrix. The possible values of the elements of the K-matrix are determined from the algebra of operators at each site. Starting in D = 3, with the K-matrix given by we have the two operators γ 1 = e ikA 2 and γ 2 = e −ikA 1 .
They will anticommute if k is odd. For the D = 2n + 1 dimensional case, the reasoning is similar. We consider a basis of fields so that the matrix K is block-diagonal (eq.(83) of the manuscript) Then, the representation (16) of the paper will provide the properly commutation rules between the operators only if all k's are odd.
Step ( In other words, the ground state will possess a Z 2k symmetry that is not present in the lattice model unless k = 1. The same reasoning goes in higher dimensions, where we have more pairs of operators satisfying the algebra (71) of the manuscript . Thus, the K-matrix corresponding to the lattice model is so that its block-diagonal form has k 1 = k 2 = . . . = k n = 1. In this sense, the K-matrix is determined by the microscopic system since it carries information about the symmetries of the lattice model.
The equation (17) of the manuscript suggest a class of equivalence for the K-matrices in our description, in a similar fashion as usual Chern-Simons theories that can lead to the same description with two different K-matrices. To make this point clear, consider a redefinition of the basis field in the effective action according to where W is a matrix with integer entries. This leads to a theory with new parameters Thus, two effective theories with parameters (K, T ) and (K,T ) related through (E6), with the matrix W possessing integer entries and det W = 1, describes the same fracton system.
Indeed, this implies and also leaves unchanged the quantization condition (for the principal configuration) given in (21) of the manuscript: which is still an even integer if I = J and an odd integer if I = J.
which corresponds to the creation of a fracton localized at x = −(a 1 + b 1 )/2, y = 0 and z = (b 1 − a 1 )/2. Notice, however, that this configuration corresponds to a process where charge is not conserved (Q = 0 → 1). Indeed, Consequently, it is not a full-fledged solution of the continuity equation.
We could try to avoid the violation of charge above by inserting a charge of opposite sign in a distinct point, which corresponds This is compatible with charge conservation in the whole system, Q = dx dy dz J 0 , but we still have to inspect the conservation in the sub-manifolds. Let us consider, for example, the following charge We need to be careful in computing the integrals, since the directions x + 13 and x + 23 are not orthogonal, whereas the directions x + 13 and x − 13 are orthogonal. This means that we can carry out the integration over x + 13 keeping x − 13 fixed. Thus, we proceed by integrating over x + 13 , letting x − 13 untouched: The computation of the remaining integral is a little trick because the directions x + 13 and x + 23 are not orthogonal, but actually we do not need to compute it to extract useful information.
Indeed, this expression shows that in order that the charge Q (++) to be conserved we need to require b 1 = d 1 . Similarly, by considering the charge we see that a 1 = c 1 in order that this charge to be conserved. The charges Q (+−) and Q (−−) do not provide additional conditions. Taking into account that a 1 = c 1 and b 1 = d 1 in (F4), we see that the density of charges J 0 trivially vanishes. In conclusion, the process of creation of a dipole is not compatible with the several conservation laws and, consequently, it is not allowed.
Let us try to find a different type of configuration, which is compatible with the whole set of conservation laws. Consider the density, and the corresponding flux which are compatible with the continuity equation (F1). It follows immediately that the charge is conserved in the whole three-dimensional manifold. Next, let us examine the conservation laws in the sub-manifolds. We start with the following charges, We have two possibilities ensuring charge conservation: Similarly, the remaining charges are which leads also to two possibilities i) a 1 = c 1 and e 1 = g 1 ii) a 1 = e 1 and c 1 = g 1 .
From these possibilities, it is clear that if we select choice i) of (F11) and i) of (F13), or ii) of (F11) and ii) of (F13), the density in (F8) will trivially vanish. However, we obtain a non-vanishing density if we choose crosswise i)/ii) of (F11) and ii)/i) of (F13). Let us choose, say, i) from (F11) and ii) from (F13). In this case, the density becomes which corresponds to the creation of four charges at the following positions: Let d(q i , q j ) be the distance between two charges. The above expressions ensure that d(q 1 , q 2 ) = d(q 3 , q 4 ) and d(q 1 , q 3 ) = d(q 2 , q 4 ), which physically means that the sum of all dipole moments of the configuration vanishes (see figure 2). This guarantees conservation of dipole moment, which is a consequence of the conservation of charges in sub-manifolds (planes). In fact, charge conservation along a plane implies that the dipole moment perpendicular to the plane is conserved.
We see that the location of the four charges in (F15) are specified by the set of arbitrary points a 1 , b 1 , c 1 , f 1 . By varying the values of these points we change both the size of the dipoles and their positions, in a way that preserves the structure depicted in figure 2, i.e., the charges are always localized at the corners of a parallelogram. Physically, this means that the dipoles can move freely in the system, but cannot be created or annihilated (remembering our previous discussion, the creation of a single dipole is not compatible with all the conservation laws).
Before closing, it is instructive to consider a simple symmetric choice, a 1 = b 1 = a and c 1 = f 1 = −a. In this case, the density reduces to while the flux can be written as where to write in terms of a single term we have used the property θ(x) + θ(−x) = 1. The density J 0 describes the creation of a set of four charges located at the points x = ±a, y = 0, z = 0 ⇒ positive charges x = 0, y = 0, z = ±a ⇒ negative charges .
This configuration is depicted in figure 3.

Case D = 5
In this case we have two conservation laws given by (51), and With this, the conservation law (F20) becomes We can construct a current J 1 that creates excitations simultaneously in planes x 1 − x 2 and x 3 − x 4 by using the four-charge configurations of the case D = 3 (F9) with the positions of the charges subject to (F15), since these configurations also live in planes. In this way, we write the generalization for the five-dimensional case as where Θ(x + 12 , x − 12 ) is defined as the θ-dependent part of (F9): Plugging J 1 in (F23) and (F25) gives the densities and J (2) 0 = θ(t) δ(x 1 ) δ(x 2 ) δ(x 5 ) ∆(x + 34 , x − 34 ) + θ(t) δ(x 5 ) Θ(x + 12 , where There are some important points to notice in the densities J 1 ) describe the same physical situation, since the redefined currents also satisfy ∂ 0J (1) and dx σ 1 15 dx σ 2 25 dx σ 3 35 dx σ 4