Towards classification of Fracton phases: the multipole algebra

We present an effective field theory approach to the Fracton phases. The approach is based the notion of a multipole algebra. It is an extension of space(-time) symmetries of a charge-conserving matter that includes global symmetries responsible for the conservation of various components of the multipole moments of the charge density. We explain how to construct field theories invariant under the action of the algebra. These field theories generally break rotational invariance and exhibit anisotropic scaling. We further explain how to partially gauge the multipole algebra. Such gauging makes the symmetries responsible for the conservation of multipole moments local, while keeping rotation and translations symmetries global. It is shown that upon such gauging one finds the symmetric tensor gauge theories, as well as the generalized gauge theories discussed recently in the literature. The outcome of the gauging procedure depends on the choice of the multipole algebra. In particular, we show how to construct an effective theory for the $U(1)$ version of the Haah code based on the principles of symmetry and provide a two dimensional example with operators supported on a Sierpinski triangle. We show that upon condensation of charged excitations Fracton phases of both types as well as various SPTs emerge. Finally, the relation between the present approach and the formalism based on polynomials over finite fields is discussed.

We present an effective field theory approach to the Fracton phases. The approach is based the notion of a multipole algebra. It is an extension of space(-time) symmetries of a charge-conserving matter that includes global symmetries responsible for the conservation of various components of the multipole moments of the charge density. We explain how to construct field theories invariant under the action of the algebra. These field theories generally break rotational invariance and exhibit anisotropic scaling. We further explain how to partially gauge the multipole algebra. Such gauging makes the symmetries responsible for the conservation of multipole moments local, while keeping rotation and translations symmetries global. It is shown that upon such gauging one finds the symmetric tensor gauge theories, as well as the generalized gauge theories discussed recently in the literature. The outcome of the gauging procedure depends on the choice of the multipole algebra. In particular, we show how to construct an effective theory for the U (1) version of the Haah code based on the principles of symmetry and provide a two dimensional example with operators supported on a Sierpinski triangle. We show that upon condensation of charged excitations Fracton phases of both types as well as various SPTs emerge. Finally, the relation between the present approach and the formalism based on polynomials over finite fields is discussed.

I. INTRODUCTION
Fracton order is a class of gapped phases of matter that exhibits a system-size dependent groundstate degeneracy on a space of non-trivial topology. This degeneracy cannot be lifted by local perturbations. Another striking feature of the fracton order is the existence of local topologically non-trivial 1 excitations with "restricted mobility" [2][3][4][5][6][7] . The latter refers to the absence of stringlike operators that would allow the local excitations to move through the space without creating additional excitations. This is in sharp contrast with, for example, fractional quantum Hall states, where quasiholes can be freely transported (provided they were localized by an external potential) 8 . These models were originally introduced as an example of glassy behavior without disorder by Chamon and as a model for stable quantum memory by Haah. The term "fracton" has been previously used to refer to small scale thermal vibrations of fractal structures 9 . We hope that the present use of the term will not cause confusion.
According to the recent nomenclature 7 , the fracton phases come in two varieties: type-I and type-II. Type-I phases support completely immobile excitations -fractons -at corners of codimension 1 surface operators (such as membranes in d = 3) 2,5-7 . Various combinations of fractons can freely move around on lower dimensional sub-manifolds 6,7 . Type-II phases support only immobile fracton excitations, that exist at "corners" of fractal operators, i.e. the non-local operators, supported on a fractal.
Despite significant research effort it is presently not known what is the appropriate mathematical structure that encodes the exotic properties of the fracton order in a model-independent fashion. First substantial progress in model-independent description of fractons was made in a series of papers [10][11][12], where it was explained that restricted mobility of type-I models can be incorporated into an effective field theory by enforcing a certain set of Gauss law constraints. These constraints lead to the conservation of the dipole moment (or, generally various multipole moments) of the matter fields. It was also shown that lattice gauge theories with such Gauss law constraints have been previously studied in [13][14]. Degrees of freedom in these theories are described by a symmetric tensor gauge field and the Gauss law is enforced upon a symmetric tensor electric field. This type of effective theories does not describe a gapped phase "as is", since there are gapless excitations. A version of Higgs mechanism was developed to remove the gapless modes in [15][16].
In a somewhat surprising parallel development it was argued that certain type-I fracton models are related to a quantum theory of elasticity [17][18][19][20] . The relation can either be argued from a duality point of view 17 , where the gauge "symmetry" emerges from solving the momentum conservation equation, or starting with the observation that the gauge transformations in symmetric tensor gauge theories are identical to linearized diffeomorphisms, which is a symmetry in theories of elastic defects. Under this correspondence the immobile fractons map to disclinations, while the partially mobile fracton dipoles map onto dislocations (which satisfy the glide constraint). We note in passing that the duality in the context of elasticity has been previously studied in a series of papers by Kleinert 21,22 , where (Euclidean) symmetric tensor gauge theories (of vector charge type) were introduced, however the glide constraint was omitted. It was further noted in [18] (and later extended in Ref. [23]) that symmetric tensor gauge theories cannot remain gauge invariant in a arXiv:1812.05104v1 [cond-mat.str-el] 12 Dec 2018 general curved background, which is in striking difference compared to tranditional electrodynamics. This observation is in correspondance with a series of works [24][25][26][27], where the effective theory of the X-cube model (a particular representative of type-I fracton models) was derived from the microscopics, studied in the presence of disclinations, and generalized to be defined on arbitrary foliated manifolds. The physical properties of these models (such as groundstate degeneracy and restricted mobility) depend on the geometry of the underlying space, thus we find it to be more appropriate to refer to fracton systems as geometric order, as suggested in [25].
In yet another parallel development it was explained that the fracton phases can be further obtained by "gauging" subsystem symmetry [28][29][30][31] . This leads to a swarm of microscopic models. At the same time, the longdistance description of the subsystem symmetries is not well-understood. A close relative of such symmetries, known as sliding symmetry, appears in a theory of smectics (as well as other layered phases) that are studied in soft condensed matter physics 32,33 . On a lattice it is possible to define a subsystem symmetry that acts on a fractal set of lattice sitee. Gaugin such symmetry leads to the type-II fracton phases 28,29 . Parity breaking phases of fractons, with possible gapless boundary modes were studied in [34,35,18]. Further work on fracton phases and related topics can be found in [36][37][38][39][40][41][42][43][44][45][46][47], while a broad picture is given in the review [48]. Finally, we note that a similar formalism has been used to describe the geometric properties of fractional quantum Hall states. Namely, the foliated spacetime has been used in [49][50][51][52] to describe the transport of energy and momentum, while the symmetric tensor degrees of freedom have appeared in the context of nematic quantum Hall states and collective magnetoroton modes [53][54][55][56][57][58][59].

A. Summary of results
The main objective of the present manuscript is to introduce a language that allows to systematically construct a rich variety of effective field theories, which exhibit the phenomenology of both type-I and type-II fracton models. This will be done in such a way that both types appear as particular cases of the general framework.
Our construction rests upon an extension of a space(time) symmetry algebra, which we dub the multipole algebra. This algebra is a natural generalization of the symmetry algebras generated by the polynomial shift symmetries studied in [60][61][62][63]. These symmetries were originally introduced in the context of Galileon gravity 60 . In systems with a conserved U (1) charge these global symmetries lead to the conservation of the (various components of) multipole moments of the charge density. The aforementioned symmetries cannot be regarded as "internal" because they do not commute with spatial translations and rotations. Throughout the manuscript the symmetry generators will be denoted as P I a . These algebras are, in general, quite delicate objects as certain consistency conditions must be satisfied. These conditions arise from the intricate interplay between the spatial and multipole symmetries. If the multipole generators are picked "at random" the algebra will only close if all spatial symmetries are discarded, leading to trivial theories. After defining the multipole algebra we explain how to construct effective field theories, invariant under its action. We will restrict our attention to the matter described either by a real scalar "phase" field or by a charged complex scalar. These theories are introduced by first constructing all possible invariant derivative operators, consistent with the multipole algebra and then including all terms, allowed by the symmetries, in the effective action. Such theories usually break spatial rotations, have quite unusual scaling properties, and are generally quite exotic. Further, we discover that some of these theories in d = 3 exhibit an enhanced sliding symmetry, alluded to in the Introduction.
Gauging of the multipole algebra should lead to exotic theories of elasticity and/or gravity, because in order to consistently gauge these theories one must gauge the spatial rotations and translations. This explains why symmetric tensor gauge theories (which are a particular case of the present construction) are very sensitive to the background geometry. It is, however, possible to "partially" gauge these symmetry algebras, under the assumption that the space(-time) part of the curvature tensor is trivial (simply put, in flat space, without torsion). This partial gauging procedure leads to a very rich set of gauge theories, which includes all known symmetric tensor gauge theories (and various variations thereof: higher rank; scalar, vector or tensor charge; tracefull/traceless, etc.) as well as the "generalized gauge theories" introduced in [64]. These gauge theories naturally satisfy the exotic Gauss law constraints, which in the present formalism, is systematically derived by gauging the multipole algebra. These Gauss law constraints can be visualized, upon discretizing on a lattice, as prescribing the "allowed" charge configurations. These charge configurations specify which excitations are mobile, which are sub-dimensional and which are fractal. To be concrete, we derive the continuous model for the U (1) version of the Haah code [64] from the symmetry principles.
The gauge theories described above are gapless and do not correspond to the gapped fracton phases. To further advance our construction we show, building upon [15], that an extremely rich variety of phases emerges after the condensation of charge k objects, which leads to the reduction in symmetry from U (1) to Z k . It is particularly interesting that depending on value of k some of the "allowed" charge configurations become redundant, while in other cases a complex charge configuration turns into a hopping operator, which hops over several lattice spacings. Using this procedure we find a version of the Sierpinski triangle models in two dimensions and the Z 2 Haah code. We then explain how to translate the obtained results into the language of polynomials over fi-nite fields. We further explain how to define multipole moments over a finite field directly from the polynomials.
The paper is organized as follows. In Section II we introduce the generalized polynomial shift symmetries and use these symmetries to motivate the multipole algebra. Next we define the multipole algebra in abstract terms and illustrate the definition on a couple of examples. In Section III we explain how to construct the invariant field theories and explain how to partially gauge the multipole algebra, directly at the level of the matter theory. We investigate several examples of such gauge theories in two and three dimensions. In Section IV we discuss various extensions of the multipole algebra, most notably the charge condensation and crystalline symmetries. We also explain the relation between the present formalism and the approach based on the polynomials over finite fields. Finally, in Section V we present our conclusions and discuss open directions.

A. Polynomial shift symmetry
The conservation of the dipole moment and its relevance to the fracton order was emphasized in [11,45]. The symmetries that have to do with the conservation of the multipole moments have been extensively studied prior to these works [60][61][62][63] and are known as polynomial shift symmetries. To simplify the presentation we will first introduce these symmetries by specifying the action on the matter fields. To be concrete, consider a real scalar field ϕ. The action of the polynomial shift symmetry is defined according to where λ α is a symmetry parameter and P α (x) is an arbitrary polynomial; the sum over α is understood. We will further assume that there is a finite number of transformations that appear in (1). If the dynamics of ϕ is described by an action S[ϕ], which is invariant under (1), then we have the following set of conserved charges where q α (x) is the charge density that can be found via a direct application of Noether's theorem. If among P α (x) there is constant polynomial -corresponding to the global U (1) charge conservation -then we can write a more intuitive expression for the charges. Denoting the charge density as ρ(x) we find Which implies the conservation of various generalized multipole moments. If we further assume that the poly-nomials P α are homogeneous we can identify (3) with the (components of) proper multipole moments. To be specific, Notice, that since the coordinates x are dimensionalfull, then so are the parameters λ α . This ultimately leads to the generalized Mermin-Wagner theorem which states that if the power of the polynomials in (1) is no larger than n, then the symmetry cannot be spontaneously broken in d ≤ n + 1 spatial dimensions 63 . This is quite similar to the Mermin-Wagner theorem for the higher form symmetry 65 . The relation between these two types of symmetry warrants further exploration.
A particularly simple case of this structure is the symmetry under all polynomial shifts of the degree no greater than n. Such transformation takes form This leads to the conservation of all multipole moments of degree less or equal to n. The conserved charges are the arbitrary moments of the density Q ij..
. . The case of n = 1 corresponds to the conservation of the dipole moment and leads to the scalar charge theory 11 , while the restricted mobility of dipoles can be added by supplementing the n = 1 symmetry with δϕ = λ |x| 2 , which leads to the traceless scalar charge theory.

B. Multipole algebra
The transformations (1) commute with each other and form an unremarkable algebraic structure. However, these transformations do not commute with the spatial symmetries: translations and rotations. Instead, the polynomial shift symmetries extend the algebra of spatial symmetries to a bigger multipole algebra, m. Consequently, the symmetries responsible for the conservation of the multipole moments are not internal. This provides a general explanation to the observation made in Ref. [18] (and later extended in [23]) that the symmetric tensor gauge symmetry is "broken" on a general curved manifold.

Intuitive preamble
Before diving into the formal details of the multipole algebra we would like to illustrate the physical origins of the structure. Consider a transformation law of a quadrupole moment of the charge density under a translation by the vector r k where Q i is the total dipole moment and Q is the total charge. Usually, one would state that provided that all lower moments vanish, the quadrupole moment is invariant under translations. Imagine a situation when Q 1 = Q 2 = 0 enforced by the symmetry and all other components of the dipole moment are non-zero. Then the quadrupole moment is invariant under translations in the x 1 − x 2 plane and we can further constrain it by the symmetries in within this plane. Effective description of systems with such constraints will be translation invariant in x 1 − x 2 plane only. We will further encounter another interesting degenerate case, which will play the central role in the discussion of the fractal phases. Consider a particular component of the quadrupole moment Q µ = Q ij µ ij , where µ ij is a certain degenerate matrix with the property that its kernel equals to the orthogonal complement of the x 1 − x 2 plane. Then such component of the quadrupole moment is translationally invariant under all translations This equality holds for all translations since if we take translations outside of the x 1 − x 2 plane the matrices µ ij will project them down to the x 1 − x 2 plane, but in this plane the projection of the total dipole moment vanishes. Thus Q µ is translation invariant. The multipole algebra formalizes these simple ideas in the language of symmetry.

Formal definition of the multipole algebra
To get some insight into the algebraic structure we need to sort the polynomials by their degree. We introduce a set of generators of the polynomial symmetries, P Ia a , so that a corresponds to the degree of the polynomial, whereas I a runs through all polynomials of degree a. Then the general multipole algebra is defined via the following set of commutation relations where O [ij] denotes anti-symmetrization over i, j and f ij J b are the structure constants, that define the algebra. In the polynomial representation the content of (10) is quite simple: rotations and translations, when applied to the polynomials, should produce linear combinations of the polynomials within the multipole algebra, i.e. the set of polynomials P Ia a is closed under rotations and translations. The multipole algebra can be exponentiated to a multipole group, M. The multipole group is, in a way, a more fundamental construct since, unlike the algebra, it admits a crystalline analogue and survives the condensation of charge p objects (provided the charge is well-defined) as a discrete group.
If the polynomials, in the polynomial representation, depend on d variables, then the indices i, j run from 1 to k ≤ d. This happens because generally the vector space of polynomials corresponding to P Ia a is not compatible with translations or rotations; in such case the multipole algebra is trivial. However, if P Ia a are chosen carefully then very rich and intricate structures can emerge.

C. Maximal multipole algebra
Next we turn to consider a few explicit examples. The simplest case is the set of transformations (6) with n = 1. It gives rise to the following algebra (we omit the usual relations between translations and rotations (9)) where P 0 refers to the constant shifts. We emphasize that the internal index I was identified with the spatial index i. This is due to an "accidental" isomorphism between the set of linear shifts and set of translations. One distinguishing feature of this algebra is that is it consistent with all spatial translations and rotations. This happens because all polynomials of degree ≤ 1 were included, consequently such set is closed under any operation that does not increase the power of the polynomial. One simple extension of the algebra (11), consistent with all spatial symmetries, is the addition of one extra generator P 0 2 , corresponding to δϕ = λ |x| 2 . This generator leads to the commutation relations We note in passing that the generators {T i , R ij , P 0 , P i 1 , P 0 2 } form the Bargmann algebra 66 , where the spatial translations T i correspond to the Galilean boosts, linear shifts P i 1 correspond to the spatial translations, P 0 2 to the Hamiltonian and P 0 to the mass central charge. It would be interesting to explore this isomorphism to construct field theories invariant under this type of multipole algebra.
In the above set of examples it is amusing to note that the commutation relations between the R ij and P k 1 are the same as between R ij and T k . In the absence of the generator P 0 2 , which provides the asymmetry, we could swap T k and P k 1 and obtain the same algebra back. This observation has a nice elastic interpretation. If the theory of elasticity is viewed as a gauge theory of translations and rotations, T d SO(d), then we could use either real translations or shift symmetries to construct such a gauge theory. This ambiguity corresponds to the two distinct approaches to the relationship between fractons and elasticity 17,18 . If the generator P 2 0 is introduced and gauged, then one will end up with two different theories.
We will elaborate on this distinction in a forthcoming work.
These symmetric cases can be generalized to arbitrary multipole moments. In order to include the multipole moments up to n in d spatial dimensions we introduce a set of polynomial symmetries described by arbitrary symmetric tensors of degree a, P Ia a = P i1 ...ia a . Each such tensor has d+n−1 n independent components. Together with translations and rotations these symmetry generators form the following algebra Theories that are symmetric under this algebra conserve all the multipole moments up to order n. Such multipole algebras are completely characterized by the spatial dimension and a single integer n, thus we will refer to those as maximal multipole algebra of order n, m n max . It may be tempting to speculate 11 that the theory obtained by gauging the algebra (13) (perhaps, in the limit n → ∞), may be related to the Vassiliev theories of higher spin gravity. In view of the algebra (13)- (14), this speculation seems unlikely since the symmetric multipole generators all commute with each other, whereas in the Vassiliev higher spin theory they form a version of the W ∞ algebra (also referred to as shs(1) in earlier works 67 ).

D. Homogeneous multipole algebra
Next we consider a less symmetric, but still very palatable case. Consider a set of polynomial symmetries where all polynomials P Ia a are homogeneous, that is of the form (4). We will refer to such multipole algebras as homogeneous. We would like to demonstrate how certain rotation and translation generators drop form the algebra due to requirement that the set of polynomials P Ia a is closed under as many translations and rotations as possible. Consider, for example, where the five polynomials P Ia a (x) are given by where µ I1 i and µ I2 ij are fixed rank-1 and rank-2 tensors and i, j = 1, . . . , d. These tensors are fixed up to the freedom of replacing µ I1 i and µ I2 ij with the linear combinations thereof.
Depending on the circumstances, these polynomials can give rise to several multipole algebras. First, assume that µ I2 ij are all non-degenerate. Then only two out of d translations survive. These are translations in the directions t j (1) and t j (2) , which are explicitly determined from where α i , β i are some constants. Eqs. (18)-(19) Next we consider rotations. Clearly, the only rotation generator that has a chance to survive is the rotation in µ 1 i − µ 2 j plane, R 12 . Within this plane we can always arrange that via redefining the symmetry transformations in such a way that, within the µ 1 i −µ 2 j plane, µ 1 ij ∝ σ 1 and µ 2 ij ∝ σ 3 , where σ i are the Pauli matrices. Thus we obtain the multipole algebra where sum over repeated I a , J a is understood. Although, we have started in d spatial dimensions, only two translation generators and one rotation generator have survived.
There is another interesting possibility, which will arise in the study of the U (1) Haah code. When both µ I2 ij are degenerate in such a way that the kernels of µ I2 ij coincide with each other and with the orthogonal complement of µ I1 i . Then Eqs.(18)- (19) hold true for all translations since µ I2 ij are projectors to the µ 1 i − µ 2 i plane. In such algebra we get an additional set of trivial commutation relations where i runs over the orthogonal complement to µ I1 i : from 1 to d − 2 in the present example.

A. General constraints
We turn to the construction of the field theories invariant under the action of the multipole algebra. First, we fix a multipole algebra m and construct an irreducible representation. To start, we have to fix the transformation law under rotations. For simplicity we take a single real scalar field ϕ 68 . The transformation laws under the action of P Ia a are given by Eq.
(1). We will denote the highest power that appears in (1) as a max . The time derivative term will take the ordinary form,φφ 69 . To construct the kinetic term, we need an invariant derivative operator, i.e. a derivative operator, consistent with (1). To this end we consider a general differential operator where s ≤ a max . The coefficients q i1i2... will be chosen in such a way that Dϕ is invariant under the action of P Ia a . In the most general case these equations take form where P Ia a is the polynomial corresponding to the action of P Ia a . Eq.(26) must hold for all P Ia a . The solutions to these equations are the differential operators D α (the index α will label the solutions). The constraints (26) can be written more explicitly if we introduce the following parametrization for P Ia a (x) Then Eq. (26) turns into a set of linear equations on q i1...
which hold for all I a . It immediately follows that q = 0. Thus, the invariant derivatives must start with at least q i1 ∂ i1 . This implies an additional invariance of the effective theory under a global U (1) transformation δϕ = µ I0 . While the multipole algebra does not require to have a constant shift as a separate generator, it seems difficult to find a matter theory that represents such algebra polynomially without a symmetry enhancement.
The system (28)-(31) contains (2a)! a! equations for every polynomial, and only (2a)! a! unknowns. Thus, generally it is severely overdetermined and has no solutions. As we will see shortly, there are, indeed, "degenerate" cases when the system does admit solutions. This phenomenon is somewhat reminiscent to the existence of a solution to the (severely overdetermined) pentagon and hexagon equations, which are the consistency conditions for fusion tensor categories.
If there are no solutions for s ≤ a max then the system can always be solved by higher order differential operator that annihilates all polynomials of degree no grater than a max . Such operator takes form for any q i1...ia max+1 . The action constructed using these solutions will have an enhanced symmetry under all polynomial shifts and will represent the maximal multipole algebra m amax max . Presently it is not clear how to establish the existence of a solution to (28)-(31) without solving O(4 a a a ) equations.
A comment is in order. In the above construction, as well as in the remainder of the manuscript, we demand the complete invariance of the action under the symmetries. It is, however, possible to consider a weaker condition -the invariance of the action up to a total derivative. This condition allows for a wider array of the invariant Lagrangians. We leave exploration of this scenario for the future work.

B. Constraints in the homogeneous case
In the remainder of the manuscript we will focus on the homogeneous multipole algebras. The constraint equations take form . . . . (33)-(35) is more overdetermined than (28)-(31), since the contractions such as µ I ijk q ij = 0 must vanish separately. Nevertheless, as we will see shortly, these systems still admit solutions.

The system of equations
Finally, we note that the solutions of (33)-(35) are determined up to an overall scale (which is not the case for (28)-(31)). We will discuss how to fix the scale later.
To the lowest order in derivatives, the invariant effective action is where α runs over all solutions of (26).
To the lowest order in gradients the effective action is quadratic in the invariant derivatives. However, the latter are not necessarily of the same degree. The derivative expansion of such theories is quite unusual. In more traditional situations global symmetries restrict the terms that appear in the effective action to all orders in derivative expansion. In the present case the situation is drastically different. To the lowest order(s) in we have a few special derivative operators D α . The leading terms in the gradient expansion must be written using D α operators only. Such terms dominate up to the order 2a max . However, when we move further in the gradient expansion, say to the order 2a max + 1, we are allowed to use any operator of the form (32).
In the case of the maximal multipole algebra m n max it is possible to write a rotationally-invariant action where D i1i2...in = ∂ i1 . . . ∂ in . Rotationally invariant actions of the type (37) were studied in [60][61][62][63]. We will view the actions (36)-(37) within the framework of the effective field theory. In particular, this entails to including all possible terms, consistent with the postulated symmetries and organizing these terms according to some degree of "relevance".

C. Multipole gauge theory
In this Section we will explain how to gauge the multipole symmetries. With the invariant derivatives at hand it is now possible to introduce a local version of the multipole symmetry. In principle, we could do it directly from the commutation relations, however this entails to gauging the spatial symmetries as well. We will leave this program to a future work.
Present approach allows to "partially gauge" the symmetries: the multipole symmetries will become local, while rotations and translations will remain global. This amounts to setting the corresponding gauge fields to 0, so that Ricci curvature and torsion vanish. The gauging procedure is well known -we promote the global symmetry in (1) to a local one δϕ = ζ(x) and introduce a covariant derivative operator where a α (t, x) is the gauge field. It is importnant to note here that a α is the space-time scalar, in that it does not transform under rotations or translations, provided the coefficients q i1... transform as proper tensors. The time derivative is replaced with the ordinary ∇ 0 ϕ =φ + χ , with δχ = −∂ 0 ζ. With the gauge fields at hand we introduce a set of conjugate momenta (or, "electric fields") according to The curvature of the connection (or, "magnetic field") is guaranteed to be non-vanishing if we gauge the entire algebra m, provided that multipole constraints do not reduce the spatial symmetries to one or zero dimensions. In the latter case the curvature will vanish identically. In present approach we will construct the magnetic fields on the case-by-case basis as as gauge-invariant combinations of a α . The general invariant Lagrangian takes form where H[e, b] schematically denotes the Hamiltonian for the gauge fields. The action is supplemented with the Gauss law constraint (obtained by integrating out χ) where the conjugate derivative D † α is defined as The Gauss law of the type (41) was postulated in Ref. [64]. We find that (41) follows directly from the underlying structure of the multipole algebra. At the same time, not every Gauss law constraint is consistent with some multipole algebra. The invariant derivatives D α can be discretized on a lattice, leading to various charge configurations. To be concrete (see Ref. [64]), consider an operator e −iaα(x) , acting on a site with label x. This operator changes the value of the electric field e α (x) at the same lattice site, say raises it by 1. This, in turn, requires to introduce electric charges at all lattice sites x , that are connected to e α (x) by the Gauss law (41). We will make heavy use of this pictorial representation.
Such configurations of finite number of point charges are characterized by a set of multipole moments. These moments are determined by q i α , q ij α , . . .. Only the lowest of these moments is independent of the choice of coordinate origin, however the solution of (33)-(35) is independent of the choice of origin. Indeed, assume that the total charge is zero, then the dipole moment is welldefined. When the origin is shifted by r i the quadrupole moment transforms as δq ij α = r i q j α + r j q i α , which still satisfies (33)- (35). In the most general, non-degenerate case the solutions are invariant only under the translations in "allowed" directions, such as the ones specified by (18)- (19). This is not too surprising since other translations do not belong to the symmetry algebra.

D. Maximally symmetric gauge theory
Next we consider a few examples of the general formalism outlined above. First, we would like to make sure that the symmetric tensor gauge theories follow. This is indeed so, provided we gauge the maximal algebra m n max . The covariant derivative takes form In this case a rotationally invariant action is possible This type of (ungauged) actions has been studied in great detail in Refs. [61][62][63], in relation to "slow" Goldstone bosons.
Including, additionally, pure trace generators of one higher degree P In+1 n+1 necessitates the change in the covariant derivative where q i1i2...in α are symmetric traceless tensors and the gauge field a α can be related to the symmetric tensor gauge field via projecting out the tracefull parts a α = q i1i2...in+1 α a i1i2...in+1 . Both traceless and tracefull scalar charge theories exhibit a scaling symmetry with dynamical critical exponent z = n + 1 Thus, in d spatial dimensions, the field ϕ is dimensionless if n + 1 = d. As discussed in [63], this corresponds to a generalized version of the Mermin-Wagner theorem: the (maximal) multipole symmetry of degree n cannot be spontaneously broken in d ≤ n + 1 dimensions. In the case n = 1 we get the effective theory for the traceless scalar charge theories where q ij are symmetric traceless tensors. The charge configurations for these theories have been discussed previously 11 . We will use this case as the first benchmark of the present formalism. In two dimensions the Gauss law takes form since there are only two symmetric traceless tensors in d = 2, namely the Pauli matrices σ 1 ij and σ 3 ij . The corresponding charge configurations are illustrated in Fig.1.

E. Quadratic multipole algebras in two dimensions
In two spatial dimensions, while restricting ourselves to quadratic order, we can solve the problem of classification of multipole algebras. The constraint equations are We consider the solutions on the case by case basis. First case is when µ I2 ij are non-degenerate. Then q i = 0 and the µ I1 i are arbitrary, which implies that the dipole moment is conserved. We have to consider a few possibilities for q ij . To start, we parametrize the symmetric tensors using Pauli matrices and the identity matrix where τ ν = (σ 0 , σ 1 , σ 3 ), ν = 1, 2, 3, and σ 0 is the identity matrix. The corresponding polynomials take form The non-trivial constraints from (49) then take form where I 2 = 1, . . . , k. The number of solutions of (52) is 3 − k. If there are no quadratic symmetries, we find the usual symmetric tensor gauge theory, associated to the maximally symmetric multipole algebra of order 1. If k = 1 and γ 1 ν = δ 1,ν , we find the symmetric traceless gauge where R is the only rotation generator and T i are the translation generators. As mentioned previously, in these cases it is convenient to use δ iI1 to identify the spatial indices with the multipole index I 1 . With this identification, the corresponding gauge fields are proper tensors. In all other cases we find theories that break rotational symmetry since general solution for q ij α takes form where α = 1, 2. There are two invariant derivatives and two corresponding gauge fields a α . The multipole algebra takes form where f i,I1 are the structure constants that are determined by µ 1 ij . Next we consider k = 2. In this case there is a single invariant derivative and a single gauge field. After gauging the Gauss law eliminates all local gauge degrees of freedom. If the quadratic symmetries are both traceless, then the only allowed q ij ∝ δ ij and the corresponding gauge field a 1 is the pure trace a 1 = δ ij a ij . This is reminiscent of the linearized dilaton coupling.

Degenerate case
More interesting structure arises in the degenerate case. The multipole symmetry contains a single linear and single quadratic term, that take form The corresponding polynomials are The dipole moment q i has to be orthogonal to µ I1 and be a null-vector of µ I2 . To following vector satisfies these criteria Furthermore, there are two quadrupole matrices that satisfy (49), which are given explicitly by where 2 , 3 are the overall length scales, which will be determined when we discretize the theory on the lattice. The multipole algebra does not contain rotations and takes a simple form where T || and T ⊥ are translations in the direction parallel and perpendicular to q i . There are three invariant derivatives, and three corresponding gauge fields, a α . There is a nonlinear relation between the invariant derivatives, which takes form This relation allows to include only the terms linear in D 1 . The Gauss law takes form The corresponding charge configurations are illustrated in Fig.2. Notice, that due to the existense of the hopping operators for the (n, m) dipole, the latter is fully mobile. However, the charges themselves are immobile. To be complete we construct a set of "magnetic fields". For the purpose of this work, we will call a "magnetic field" any gauge-invariant combination of the gauge fields a α . Such construction turns out to be quite a non-trivial problem. To illustrate the sublety, we first proceed in a "naive" way. By inspection we find 4 "magnetic fields" where the first three are gauge invariant by the virtue of the commutativity of the partial derivatives, while the latter is invariant due to the non-linear relation (65). It turns out, that these magnetic fields satisfy several constraints. These constraints can be obtained by applying various invariant derivatives to b 4 . We find the following constraints These constraints leave a single independent component of the magnetic field. Finally, we could impose only a single linear polynomial symmetry and no quadratic symmetries. This will lead to a single choice of q i , but no restrictions on q ij . Such theories brake rotational symmetry.

F. Anisotropic scaling
Next, we would like to investigate the scaling properties of the two-dimensional theories. We will focus on the case when a single component of the dipole moment is conserved. The Lagrangian does not possess any peculiar scaling properties for generic values of m, n and the coupling constants. However, there are a few interesting cases to consider. To make the scaling properties apparent we introduce the variables q = q i x i /|q| and r = µ 1 i x i /|µ 1 |. In terms of these variables the Lagrangian takes form where λ α are linear combinations of λ α , with the coefficients determined by m, n. The lack of ∂ 2 r ϕ term is a manifestation of the degeneracy of the theory. If λ 2 = 0 then the theory exhibits the following scaling symmetry This symmetry implies scale invariance in the direction of the conservation of the dipole moment. If such scaling symmetry is enforced, it leads to a relation between the coupling constants λ 2 = − m 2 n 2 λ 3 . Next we relax the quadratic multipole symmetry. The multipole algebra, in this case, still closes and still does not contain the rotational symmetry since we have picked a preferred direction (namely, the conserved component of the dipole moment). For such theory the Lagrangian is no longer degenerate and is given by (73) plus an extra term λ 4 (∂ 2 r ϕ) 2 . Such theory acquires an interesting scaling symmetry if we set λ 2 = λ 3 = 0 where the scaling symmetry takes a highly anisotropic form If such symmetry is enforced, only two charge configurations are allowed. We leave the detailed investigation of such symmetries to future work. Clearly, anisotropic scaling is only possible if the systems in question break rotational symmetry, which is often the case in condensed matter physics. It is possible that such symmetries naturally emerge close to a quantum critical point that exhibit dimensional reduction 70,71 .

G. U (1) Haah code in three dimensions
Next we turn to the "U (1) Haah code" studied in [64]. We start by postulating the symmetries where The polynomials can also be represented by coefficient matrices as in (16)-(17) The multipole algebra takes exactly the form (21)- (24). The dipole and quadrupole vectors are found by solving the constraints (33)- (35). We find five solutions which are explicitly given bȳ where¯ α are the overall scales, among which¯ 0 is dimensionless, so we will set it to 1. Recall, that these scales cannot be determined from the constraints alone.
To lighten up the equations we will set all¯ α = 1, however the reader should be keenly aware that there is some freedom in overall scales. Thus we have five invariant derivativesD An invariant Lagrangian of the form (40) can be written already at this stage. It turns out, however, that such Lagrangian is not invariant w.r.t. rotations in the µ 1 i −µ 2 j plane. We can additionally enforce the invariance under rotations. Technically this is done by taking the linear combinations of the tensorsq ij α that are invariant under such rotations. There are two such linear combinations Thus we have three invariant derivatives that take form while the covariant derivatives are given by ∇ β ϕ = D β ϕ + a β . The most general Lagrangian, consistent with the (gauged) multipole algebra takes form where . . . stands for the terms higher in derivatives, g αβ is the matrix of coupling constants and H[e, b] ∝ β e 2 β + b 2 is the Hamiltonian for the gauge fields.
The Gauss law takes form The elementary charge configurations are illustrated in Fig.3. Finally, we note that there is a non-linear relation among the invariant derivatives, namely, which corresponds to These relations reduce the number of independent higher order terms and restrict the magnetic fields, but do not allow to eliminate any of the invariant derivatives.
Next we turn to the construction of the "magnetic fields". As before, the latter are defined to be the gauge invariant combinations of the gauge fields. There are 4 gauge invariant functions that we construct by inspection These magnetic fields are not all independent. The magnetic field b 4 is gauge invariant by the virtue of the nonlinear relation between the invariant derivatives (92). Furthermore, it can be explicitly checked that there are three relations between the magnetic fields where we have used (92). We again find that only one component of the magnetic field is independent. The above theory includes the "generalized gauge theory" for the U (1) Haah code of Ref. [64]. It appears that in Ref. [64] the authors only kept the following invariant derivatives We are not aware of an additional symmetry principle that would force us to discard D 2 . Since the Lagrangian (77) is effective, it must include all terms allowed by the general principles. Addition of an extra derivative, and, therefore an extra charge configuration does not contradict the conclusion about "fractal dynamics" observed in [64] as we will discuss in the next Section. The Lagrangian (90) has a hidden conformal sliding symmetry. To see it, we introduce new variables x = µ 1 i x i /|µ 1 | and y = µ 2 i x i /|µ 2 |. Then all invariant derivatives D α (and, consequently the Lagrangian) are also invariant under an infinite symmetry δϕ(z,z, x 3 ) = f (z) + g(z) , z = x + iy , where f (z) is holomorphic and g(z) is anti-holomorphic. This realization of conformal symmetry is an exotic example of a well-known "sliding" symmetry 32,33,72 , that appears in physics of smectics 73 and it can be understood as a continuous version of sub-systems symmetries. This symmetry is responsible for an infinite number of conserved charges noticed in [64]. Presently, it is not clear whether sliding symmetries are a generic feature of the models invariant under the multipole algebra, or it is an accident of the Haah code. Finally, the Lagrangian (90) exhibits an anisotropic scaling symmetry, which takes form We leave the investigation of the physical consequences of this symmetry to future work.

H. Coupling to charged matter
We will consider charged matter, represented by a complex scalar field. This is not the most general situation, since the matter fields will not transform under rotations, but it will serve a good illustrative purpose. The inspiration for the following construction is taken from [19,74].
In the previous Section we have explained how to construct the invariant derivates for an arbitrary charge conserving multipole algebra. Those derivatives were used to couple a phase field ϕ minimally to the multipole gauge fields a α . To introduce the charged matter we view the phase field as a phase of a charged scalar, according to We will concentrate on the homogeneous multipole algebras. In this case the invariant derivatives are also homogeneous and can be ordered by the degree as follows The covariant derivatives of the complex scalar are defined according to It immediately follows that using (102) in (104)-(105) leads to the invariant derivatives acting on ϕ that we discussed previously. The invariant Lagrangian then takes form where the terms are arranged in such a way that global U (1) invariance is preserved and g αβ is a matrix of coupling constants. It is an open problem to construct an analogue of such formalism for inhomogeneous multipole algebras as well as Lagrangians invariant up to a total derivative.

IV. EXTENSIONS
In this Section we discuss two extensions of the present formalism. One extension includes the point group symmetries of the lattice, whereas the other one includes charge condensation. We will also explain the relation between the present ideas and the formalism of polynomials over finite fields.

A. Crystalline multipole algebra
We have already observed that not all multipole algebras are consistent with continuous spatial rotations.
We have also noted that the symmetry parameters have the dimension of length and, ultimately, have to be determined by the lattice constant. Keeping the lattice physics in mind, we can relax the rotational symmetry from continuous to a point group symmetry. Presently there is no general theory of the multipole algebras combined with the crystalline symmetries. Instead, we consider an example -C 4 symmetry in two spatial dimensions. The polynomial symmetries compatible with C 4 are whereas if we were to require continuous rotational symmetry we would find a single quartic polynomial P 4 = P 1 4 + 2P 3 4 = (x 2 1 + x 2 2 ) 2 . Thus we find an interesting phenomenon: restricting continuous spatial symmetries to the crystalline ones (i.e. reducing the symmetry), allows for increasing the multipole symmetry. This will ultimately lead to more intricate constraints on the effective Lagrangian for the "spin-4" field.

B. Charge condensation
Lattice fracton models usually do not have a U (1) integer charge, rather they possess a Z p symmetry, which means that the charge lattice is reduced to Z/pZ, i.e. the excitations of charge p can disappear into vacuum and are equivalent to charge 0. This constraint is particularly effective if we have already introduced a lattice.
In light of this possibility, we will revisit the models from the previous Section. We start with a two dimensional model, characterized by the symmetry algebra (58). In this model condensation of charge-2 and of charge-3 objects leads to dramatically different macroscopic behavior. We will consider the case when n = m = 1. When charge-2 objects are condensed, we can modify Fig. 2, to the Z 2 -valued charges, as illustrated in Fig. 4. The (1, 1) dipole configuration, corresponding to D 1 , turns into an ordinary hopping operator, in the (1, 1) direction. The Z 2 charges can be easily separated as shown in Fig. 5, but only along the (1, 1) direction. Thus such charges are dimension-1 particles, capable of hopping in (1, 1) direction only.
In the Z 3 case, when charge 3n objects are equivalent to vacuum, the D 1 charge configuration is no longer a hopping operator since Q = −2 ∼ Q = 1. The charge configurations now take form illustrated in Fig. 6.
If we try to separate the charges created by e ia1 we find a fractal structure (see Fig. 7), which correspond to Z 3 version of Sierpinski triangle. This structure has appeared in [5]. We note, however that according to the general results of [3,75] such theories cannot be topologically ordered. Meaning, that they are not stable to perturbations that break the multipole symmetry. Rather, The charges can be hopped by applying an operator living on a fractal structure. Yellow stars correspond to the operator D1, while green stars correspond to the operator 2D1. It is clear that the ability to separate charges becomes sensitive to the system size since the linear size of the operators is 3 k . The next hopping operator will be of the size 27. Both charge 1 and charge 2 excitations can be hopped using similar operators. To hop charge 2 excitations one must replace all green stars by the yellow ones and vice versa. they should be viewed as SPTs.
Next we turn to the U (1) Haah code in three dimensions. In the previous section we have found that any charge configuration is generated by a combination of 3 basis charge configurations (see Fig. 8). In the original version of the Haah's code, which is based on the Z 2 charge, there are only 2 charge configurations. The two facts are reconciled by observing that after charge-2 condensation one of the charge configurations can be eliminated via applying the configuration corresponding to D 1 multiple times as shown on Fig. 8. In the Z 3 version, we find three independent charge configurations.

C. Polynomials over finite fields
The original works on the fractal phases 3,5,75 use the language of polynomials over the finite fields. In the present case, the "creation operators" originate in the structure of the invariant derivatives. The latter become differential operators with the coefficients in the same finite field upon the charge condensation. In this Section we describe the relation between the formalism of polynomials over finite fields and field theoretic approach.
To get some intuition about the possible relation we convert the graphical representation of charge configurations into polynomials as explained in [3,5,75]. In this construction one considers formal multivariate polynomials over a finite field, say Z p , for a prime p (or over Z in the U (1) case). The coefficients of the polynomials give the values of charges, while the powers of formal variables provide the coordinates. For example, qx n y m corresponds to a charge q located at position (n, m).
Polynomials H 2 and H 3 are divisible by H 1 over Z 3 . The coefficients in the polynomials sum up to 0 mod 3, which reflects the conservation of charge mod 3. The polynomials H i satisfy the same relation as the invariant derivatives (65) The use of these polynomials guarantees that all charge configurations will satisfy the conservation laws. The hopping of the dipole can be implemented in two different ways: additively and multiplicatively. The latter is accomplished via multiplication of the polynomial by either x or y. Indeed xH 1 and yH 1 correspond to the (1, 1) dipoles that hopped either in x or in y direction. The "additive" hopping polynomials can be constructed as −H 1 + xH 1 = 2H 1 + xH 1 = H 2 . Then, indeed, H 1 + H 2 = xH 1 . Note, that it is the additive hopping operators that correspond to the charge configurations Fig. 6. Since all the polynomials have a common factor, x + y + 1, it follows that the configuration corresponding to this common factor is mobile. Monomials x or y themselves are not allowed due to charge conservation. The pair creation, corresponding to x + 2y or y + 2x is also not allowed by the conservation of the dipole moment. The first allowed process is the triple creation, corresponding to x + y + 1. This triple combination is mobile since the operators which hop it are also allowed by the conservation laws. The fractal operators of the Fig. 7 are constructed by simply considering the powers H 3k 1 . It appears that there is enough information in the field theory to construct all of the polynomial and fractal structure after specifying the lattice and the condensation process. Alternatively, it should be possible to arrive to the same set of polynomials imposing the constraints (i.e. multipole moment "conservation") directly in the polynomial language.
The Haah code is specified by a similar data. In the Z 2 case the relevant polynomials are The redundancy of one of the charge configurations, in the Z 2 case discussed in the previous Section, corresponds to a relation H 2 1 = H 2 . One way to understand this nonlinear relation is to note that in (92) multiplication by 2 is the same as multiplication by 0 over Z 2 (which is no longer true over Z 3 ). This time, however, there are no mobile operators since multiplication by x, y or z takes us outside of the allowed set of polynomials. The fractal structures are generated by H 1 to various powers.

D. Multipole moments over finite fields
General multipole moments can be constructed using the formal derivatives over finite fields. We will illustrate the construction on the example of quadratic polynomials. Consider a polynomial with the coefficient either in Z p for a prime p or in Z.
The total charge of the corresponding charge configuration is given by the sum of the coefficients within the appropriate field The conservation of charge states that we consider the polynomials satisfy H({x i } = 1) = 0. The dipole moment can be evaluated as follows. Note that the power of a monomial indicates position of the charge. Thus to get a the value of the position we have to take a derivative. To this end we construct a vector of polynomials The dipole moment is determined by summing the coefficients in every component within the field. The sum over coefficients is formally evaluated by setting One has to be careful with the "multiplication". The symbol 2h kk really means 2h kk ≡ h kk + h kk mod p. This definition automatically allows to mod out by the equivalence relations between the values of the dipole moment. We now give an example of such relations in the case of Z 2 . Consider a charge configuration H = x 2 +xy +y 2 +1. The total charge is 0 over Z 2 (and would be 4 over Z).
The dipole polynomial takes form (117) The latter equivalence between the dipole moments is a consequence of the fact that charge-2 excitations can be pulled out of the vacuum and shift the total dipole moment by δ d = 2n 2m .
Arbitrary k-th multipole moment of the charge density can be defined in a similar fashion, if we restrict ourselves to configuration with all vanishing lower moments. For example, the quadrupole moment is constructed from the matrix of second derivatives. In general, we have for the k-th moment Q i1i2...i k To check that these relations are true one has to (i) construct the k-th moment according to the usual definition and (ii) use the constraints that all lower moments vanish. In the case when some of the lower moments are non-zero, the construction has to be appended. Finally, we would like to demonstrate that the multipole moments so defined are consistent with the multiplicative realization of the translations. For brevity we demonstrate this on the example of the second moment. Translation in the l-th direction by n lattice spacings is realized via multiplication by x n l . The change in the second moment is the given by δ l Q ij = ∂ i ∂ j (x n l H) − ∂ i ∂ j H = δ i l ∂ j H + δ j l ∂ i H + n(n − 1)δ i l δ j l H + (x n l − 1)∂ i ∂ j H , evaluating δ l Q ij at {x i } = 1 we find δ l Q ij = nδ i l d j + nδ j l d i + n(n − 1)δ i l δ j l Q = 0 , This variation vanishes provided all lower moments -total dipole and total charge in the present case -vanish.
In the language of polynomials the conservation laws are implemented as brute force constraints on the various moments of the charge density. It is not clear to us how to introduce the finite field version of the polynomial shift symmetries.

A. Conclusions
We have introduced the multipole algebra -an extension of space(-time) symmetries that enforce conservation of certain multipole moments of the charge density. This algebra contains both spatial symmetries and, in the simplest scalar representation, the polynomial shift symmetries. We have explained how to gauge the latter in the flat space and have shown that the corresponding gauge theory satisfies a set of Gauss law constraints. These constraints imply that the local excitations correspond to certain charge configurations with prescribed moments of the charge density. In such models one encounters a difficulty in trying to separate the U (1) charges away form each other. The recently studied symmetric tensor gauge theories of various kinds naturally fall into this structure and correspond to the maximally symmetric homogeneous multipole algebras. Crucially, the (gapless versions of) type-II fracton models also fall into the same category, and correspond to "less symmetric" multipole algebras.
We have discussed several concrete examples of the multipole algebras. The U (1) version of the Haah's code fits naturally into this structure. Upon charge condensation from U (1) to Z 2 we find exactly the charge configurations considered in the original Haah's work. It was found that for Z 3 charges there is an additional basis charge configuration that cannot be ruled out on the basis of the symmetry alone. We have also discussed a twodimensional example where the fractal structures naturally emerge upon the charge condensation from U (1) to Z 3 . Such 2D theories are gapped, but not topologically ordered. Thus, such theories should be viewed as SPTs of the multipole symmetry. Finally, an explanation relating the present construction to the formalism of polynomials over finite fields was provided.

B. Discussions
In this final part we discuss some open problems remaining after this work. First and foremost, we were able to formulate a general structure which, upon charge condensation may (or may not) lead to the fractal operators. It is important to find the necessary and sufficient conditions, in the field theory language, for the appear-ance of these operators. Currently, we can only establish their presence by inspection.
As our two-dimensional example illustrates it is possible to have fractal operators without topological order. Such theories are not stable to perturbations that break the multipole algebra. At the same time we were able to reproduce the Z 2 Haah code within the same framework. The latter, however, is topologically ordered and is stable against local perturbations, including the ones that break the multipole symmetry. It is not clear how to establish the existence of the topological order without going into details and comparing with the known commuting projector models.
The polynomial symmetries discussed above are clearly well-defined on an infinite plane. When the system is placed on a torus, i.e. is subject to the periodic boundary conditions, the polynomial symmetries become inconsistent with the boundary conditions. If the field ϕ is assumed to be compact, then we need to ensure that the exponents of these polynomials, e iP Ia a , are consistent with the boundary conditions. Restricting to such polynomials will lead to the reduction in the number of symmetries. It would be interesting to see if identifying such polynomials provides the information about degeneracy on a torus as well as an indicator that signals whether the "Higgsed" theory is topologically ordered or not.
On a more formal field theory side, it would be interesting to develop a general procedure that allows gauging of the entire multipole algebra. Such gauging should lead to very exotic theories of gravity and/or elasticity. Partial progress on this topic has been made in regards of gauging the Bargmann algebra, which we have encountered upon studying traceless scalar charge theory 66 . It will also be interesting to understand how the multipole algebra manifests itself in the theory of elasticity and its dual gauge theory along the ideas of [76]. We plan to address these and other questions in a forthcoming work.