Multipole groups and fracton phenomena on arbitrary crystalline lattices

Multipole symmetries are of interest in multiple contexts, from the study of fracton phases, to nonergodic quantum dynamics, to the exploration of new hydrodynamic universality classes. However, prior explorations have focused on continuum systems or hypercubic lattices. In this work, we systematically explore multipole symmetries on arbitrary crystal lattices. We explain how, given a crystal structure (specified by a space group and the occupied Wyckoff positions), one may systematically construct all consistent multipole groups. We focus on two-dimensional crystal structures for simplicity, although our methods are general and extend straightforwardly to three dimensions. We classify the possible multipole groups on all two-dimensional Bravais lattices, and on the kagome and breathing kagome crystal structures to illustrate the procedure on general crystal lattices. Using Wyckoff positions, we provide an in-principle classification of all possible multipole groups in any space group. We explain how, given a valid multipole group, one may construct an effective Hamiltonian and a low-energy field theory. We then explore the physical consequences, beginning by generalizing certain results originally obtained on hypercubic lattices to arbitrary crystal structures. Next, we identify two seemingly novel phenomena, including an emergent, robust subsystem symmetry on the triangular lattice, and an exact multipolar symmetry on the breathing kagome lattice that does not include conservation of charge (monopole), but instead conserves a vector charge. This makes clear that there is new physics to be found by exploring the consequences of multipolar symmetries on arbitrary lattices, and this work provides the map for the exploration thereof, as well as guiding the search for emergent multipolar symmetries and the attendant exotic phenomena in real materials based on nonhypercubic lattices.

Prior explorations of multipole symmetries have largely been limited to either systems in the continuum or on hypercubic lattices.It is worth noting that formulating the problem in the continuum introduces certain pathologies -for instance, the dipoles and multipoles are not quantized and there appear a continuous infinity of particle types and superselection sectors.How to properly define the problem in the continuum is a program of ongoing research [45].In practice, most works implicitly assume an underlying lattice.When lattice symmetries are treated seriously, the problem is almost always formulated on a square or cubic lattice, with rare exceptions (e.g., [46]).However, given the wide range of crystal structures, it is natural to wonder what happens if we formulate the problem on an arbitrary crystal structure that is not a hypercubic lattice.Multipole symmetries are not purely internal, and mix with spatial transformations [6], so the extension to arbitrary lattices is decidedly nontrivial.As a simple example to illustrate this nontriviality, let us work in one dimension and imagine imposing dipole conservation, but not monopole conservation.This automatically implies that the theory must lack translation invariance -the dipole moment is defined with respect to an origin, and monopole charge may be freely added or removed at that origin without changing the dipole moment.Conversely, if we wish to retain translation symmetry, and conserve dipole moment, then we must also conserve monopole moment (charge).What happens if we move up from one dimension to higher dimensional lattices?How does the set of multipole symmetries that can be consistently imposed depend on the choice of lattice?Are there qualitatively new phenomena that arise once we move away from continuum systems or hypercubic lattices?These questions remain largely open, and could provide routes to the design of new fracton phases on non-hypercubic lattices, to the realization of new kinds of nonergodic dynamics, and to the identification of new hydrodynamic universality classes.They would also guide our search for such phenomena in real materials -since while multipolar symmetries may emerge in the low-energy description of a real material, which multipolar symmetries emerge would depend on the crystal structure of the material, which may not be hypercubic.
In this work, we undertake a systematic exploration of multipole symmetries on arbitrary lattices.We explain how, given a space group symmetry and a set of occupied Wyckoff positions (which together determine the crystal structure), one may systematically construct all possible consistent multipole symmetry groups to any desired order.While we work in two dimensions for simplicity, the methods we develop are general and should extend to three (or higher) dimensions mutatis mutandis.For all two-dimensional Bravais lattices, we exhaustively classify all possible consistent multipole groups at order  = 1 (dipole),  = 2 (quadrupole),  = 3 (octupole) and  = 4.We explain how the classification may be extended beyond Bravais lattices to deal with bases and nonsymmorphic symmetries, thereby allowing us to access arbitrary wallpaper groups.We also explain that the space group itself is not sufficient to fully specify the problem -one needs the full crystal structure, which also involves knowledge of the occupied Wyckoff positions.We illustrate the general framework by an explicit computation of the consistent multipole groups up to order  = 2 on the kagome and breathing kagome crystal structures.This gives an in-principle classification of all possible multipole groups in any space group.
While the first part of this manuscript is essentially mathematical in nature, classifying the possible consistent multipole groups on various crystal structures, the second part of this manuscript discusses the physical consequences.We begin by discussing how knowledge of the multipole group may be used to write down effective low-energy field theories, and how these may be discretized to yield effective Hamiltonians on the lattice in question.We then discuss how this framework may be used to generalize certain results originally obtained on hypercubic lattices.For instance, we present a general understanding of a minimal set of symmetries that must be imposed to yield localization on arbitrary crystal lattices, generalizing the results of Ref. [30,31].Then we discuss two seemingly novel phenomena that arise when we move beyond hypercubic lattices (a) an emergent robust subsystem symmetry arising on the triangular lattice, and (b) an unusual situation arising on breathing kagome, where one can obtain multipolar conservation laws without conservation of monopole charge, but while retaining translation symmetry (and also an emergent vector conserved charge).These two examples are not exhaustive, but illustrate that new physics can arise when one generalizes away from hypercubic lattices.The systematic exploration of new physics arising from multipole symmetries on arbitrary crystal structures therefore promises to be a fertile territory for exploration, and this manuscript provides the map.
This manuscript is structured as follows.We start by introducing the methodology for deriving multipole groups that are compatible with the space group of the lattice in Sec.II.For readers not interested in the technical details, we begin this section by presenting an intuitive overview of the general procedure.We then apply the formalism to the five Bravais lattices in two dimensions in Sec.II B, and to generic wallpaper groups in Sec.II C, where a number of additional complexities arise.The second half of the manuscript is concerned with exploring the consequences of the multipole groups found in Sec.II.First, we describe how to construct local Hamiltonians that are invariant under space group operations and conserve multipole moments belonging to a particular multipole group in Sec.III.Then, in Sec.IV, we extend some classic results derived on hypercubic lattices to arbitrary crystal structures, and also present two examples of interesting phenomena that can arise when one goes beyond simple hypercubic lattices.We close with a discussion of our results in Sec.V. Finally, some clarification on notation.Throughout the manuscript, we use the physics convention for the dihedral group: D  is the group of symmetries of a regular -gon.Similarly, our notation for the irreducible representations (irreps) of D  is set out in Appendix A.

II. FORMALISM FOR DERIVING LATTICE MULTIPOLE GROUPS
The goal of this section is to build a procedure that algorithmically determines the constraints that space group symmetries place on lattice multipole groups.More precisely, suppose the multipole group contains a particular polynomial shift symmetry f.Then we wish to determine: what additional polynomial shift symmetries must appear in the multipole group in order to preserve the space group symmetry?
We will first describe the general procedure.Then, we will implement it for all five Bravais lattices in 2D, explaining several details and subtleties that arise in the different examples.Finally, we will explain through examples several important concepts that appear when classifying multipole groups on a crystalline lattice with a basis.Using Wyckoff positions, we will give an in-principle classification of all possible multipole groups in any space group.

A. General procedure
The input to our procedure is a space group , its action on a given basis of the crystalline lattice L, and a nonnegative integer .The output is a list of all multipole groups M compatible with the space group  (in the sense the quotient of M by its subgroup of pure polynomial shift symmetries is ) that contain polynomial shift symmetries of degree at most .
We place a scalar field   (r) at each basis location in the crystalline lattice.Here  labels a basis element of L, and, depending on the physics, r may either be the position of a lattice site or a basis site; we will discuss the distinction in Sec.II C 3. In what follows, we assume for simplicity that, given an element s ∈ , the field transforms as where s() is some permutation of the lattice basis elements and s(r) is some action on the spatial coordinates.More generally, the field  need not transform as a scalar under space group operations.This more general case is easily incorporated into the formalism; we give an example in which the field transforms as vector and discuss some consequences in Appendix B. Under an infinitesimal polynomial shift symmetry f, we have where the   (r) are a collection of polynomials that we assume to have degree at most , and  is the symmetry parameter.A symmetry under such a polynomial shift implies conservation of the corresponding multipole moment.We demand that the set of polynomial shift symmetries be closed under the action of the space group -a defining property of the multipole group [6].Commuting s through the polynomial shift symmetry, we see that if f is in M, then M must also contain Given f, then, we must find its orbit under the space group.While this problem is well-posed, the space group has infinite order, although it is finitely generated (e.g., there are an infinite number of symmetry translations, but they are all generated by a finite set of basis vectors), and the classification problem becomes awkward1.However, in general, space group elements act on spatial coordinates as where  s is a matrix acting on the components of r and t s is some constant vector.Therefore, if we replace s(r) with a modified transformation s that acts as that is, we completely drop the translation, then applying this transformation to homogeneous polynomial of degree  produces another homogeneous polynomial of degree .We can therefore sort the set of homogeneous polynomials   (r) into irreps of the finite group generated by the transformations which is a straightforward task that can be done algorithmically.
There are several ways to accomplish this using standard grouptheoretic methods; see Appendix A for a numerical method using character theory and an analytical method using Clebsch-Gordon coefficients for discrete groups.
We call the group of modified transformations the "extended point group"; the reason for the terminology is that this group evidently contains the point group but, for a non-symmorphic symmetry group, is in general larger than the point group.The extended point group is in fact isomorphic to the quotient /, where  ⊂  is the set of translations, but the extended point group is implemented slightly differently because it discards 1 The basic technical problem is the following.Consider the vector space spanned by monomials, equipped with the natural inner product where different monomials are orthogonal.Then translations do not act as orthogonal operators on that vector space.We would thus be forced to consider nonorthogonal representations of the space group, which is a more challenging problem than the representation theory considered in this paper.both integer and fractional translations.This does not yet solve the problem, because under a true space group operation (including translations), a representation of degree  will generically produce terms of all lower degrees.However, the terms of degree  will stay within the representation we found earlier because translations only produce terms of degree less than .Hence, under a true space group operation, the representations of degree  "mix" with representations of degree less than , but do not mix with other representations of degree .In fact, we show in Appendix C that for each representation of degree , we may take an infinitesimal translation by t s , which produces homogeneous polynomials of degree  − 1 only.By iterating this procedure on all representations of degree less than , we find the full set of polynomials that appear under these translation operations.
To summarize, we sort homogeneous polynomials into irreps of the extended point group (6), and then see how they "mix" under the translation piece of each space group element.Under (true) space group transformations, an element fof a fixed irrep produces other elements of its irrep and linear combinations of its "descendant" irreps.If f ∈ M, then the entire irrep to which f belongs, and all of the descendants of that irrep, must appear in M as well.
This procedure can be generalized to the case where f is an arbitrary linear combination of elements from different irreps of the extended point group.We will see in the context of various examples that this generalization can be subtle.
Before proceeding to examples, we comment that, for convenience, we have assumed that the multipole group contains polynomial shift symmetry for scalar fields, but the approach would work equally well for, say, spins, e.g., rotating Ŝ(r) →   Ŝ (r)  (r) Ŝ(r) − Ŝ (r)  (r) (7) where  (r) is some polynomial.We will elaborate further on this context in Sec.III A.

1
x y x 2 2xy y 2 x 3 3x 2 y 3xy 2 y 3 x 4 4x 3 y 6x 2 y 2 4xy 3 y 4 FIG. 2. Mixing of monomials under lattice translations for the monoclinic and orthorhombic Bravais lattice types.Every arrow means that some (infinitesimal) translation acting on the monomial of higher degree generates the monomial(s) of lower degree.That is, if a particular monomial is included in the multipole group, compatibility with space group operations requires the multipole group to contain all lower-order monomials that can be reached by following an arrow from the original monomial.Linear combinations of degenerate irreps must be dealt with separately, as described in the main text.

B. Bravais lattices
When working with Bravais lattices, there is no basis index  in, e.g., Eq. ( 1).As a result, the procedure outlined above becomes straightforward.Extended point group operations (6) are exactly point group operations, and so we just need to sort homogeneous polynomials into irreps of the point group .The "mixing" of irreps then comes from translations by lattice vectors.We will find the possible multipole groups compatible with all five Bravais lattices in two dimensions (depicted in Fig. 1), starting with the lattice that possesses the fewest symmetries.As we increase the size of the symmetry group, various additional complexities will appear; while the formalism will remain unchanged, we will highlight some subtleties in its implementation.

Monoclinic
We begin with the monoclinic Bravais lattice [Fig.1(a)], whose point group is C 2 , with generator (, ) ↔ (−, −).In crystallographic notation, which will be used throughout the paper, the space group is 2.By inspection, each monomial     (with ,  ≥ 0) of even degree forms a trivial irrep of C 2 and each monomial of odd degree forms the nontrivial one-dimensional irrep of C 2 .
After sorting polynomials into irreps of the point group, we must determine what constraints translations put on the multipole group.Suppose that we include     in the multipole group for a given pair , .As shown in Appendix C, we must also include shift symmetries obtained by infinitesimal translations of     along the lattice directions.Since there are two linearly independent translations and polynomial shift symmetries can have any real (not just integer) coefficients, it suffices to choose the translations to be along  and  (denoted Iterating this procedure, we find that including  −1   in the multipole group requires  −2   (if  > 1) and  −1  −1 (if  > 0) to appear in the multipole group as well.Each irrep of degree  therefore has a set of "descendants" of degree  − 1, and the descendants have their own descendants, and so on.We can organize this information into a graphical tree-like structure, presented in Fig. 2. The meaning is that if we include any one monomial in M, then all of its descendants (i.e., any monomial that can be reached by following a series of arrows starting at the original monomial) must also appear in M.
We now address an important subtlety.Each irrep of the point group is highly degenerate, in the sense that any linear combination of even-(or odd-) degree monomials also forms an irrep of C 2 .This is unimportant if we include in M a linear combination of polynomials transforming under different irreps; in that case, we have to include both polynomials or neither in the multipole group.For example, including  4 + 3 2 as a generator in M also forces us to include  4 −3 2 since one term is odd under rotations and one term is even.The polynomials { 4 + 3 2 ,  4 − 3 2 } span the same set of polynomials as { 4 , 3 2 }, so we can just include both monomials.To ensure the multipole group is closed under translations, we should then include the descendants of both monomials.
However, if we include a linear combination of two identical representations, we may only need to include corresponding linear combinations of the descendants.For example,  3  +  3 only forces us to include the span of {3 2  +  3 ,  3 + 3 2 } (and lower-degree descendants), not the span of { 3 ,  3 , 3 2 , 3 2 }.
Mixing of irreps of the point group D 6 by translation symmetry for the triangular lattice.Every arrow means that some (infinitesimal) translation acting on the polynomial(s) of higher degree generates the polynomial(s) of lower degree.As for the square lattice, all multipole groups that contain a polynomial of degree  ≥ 2 must also include both components of dipole, and charge.The representations of the polynomials under the point group are given in Table II.
The tree in Fig. 2 is therefore still useful for finding the possible multipole groups, but one must be careful when including linear combinations of polynomials transforming under the same irrep of the (extended) point group.

Orthorhombic
Next, we consider Bravais lattices of orthorhombic type, for which the point group is D 2 [Figs.1(b), (c)].For simplicity we choose coordinates so that the point group is generated by reflections about the and -axes.Again, each monomial forms an irrep of the point group.Choosing the generating mirror of D 2 to send  ↔ −, the monomial     forms the 1D irrep  (−1) + , (−1)  where the notation is set in Appendix A; to briefly summarize, the first index is the eigenvalue under two-fold rotations [i.e., (, ) ↔ (−, −)] and the second index is the eigenvalue under the generating mirror.
For the rectangular Bravais lattice (space group ), the monomial     transforms by  −1   (provided  > 0) under   and    −1 (provided  > 0) under   .For the rhombic Bravais lattice (space group ), the monomial     transforms by  −1   under   and a linear combination of  −1   and    −1 (again, assuming ,  > 0) under translations parallel to the other lattice direction.In both of these Bravais lattices, including     in the multipole group requires both  −1   and    −1 to appear in the multipole group as well.Hence, the same sets of polynomial shift symmetries are allowed in the multipole group as in the monoclinic Bravais lattice.We note, however, that there are now four one-dimensional irreps for the orthorhombic lattices.At a given order, there are hence fewer irrep degeneracies that need to be taken into account when considering mixing of irreps under translations than for the monoclinic Bravais lattice.

Square
The square lattice's point group is D 4 [Fig.1(d)], and its space group is 4.In contrast to the lower-symmetry Bravais lattices, sorting the square lattice polynomials into irreps from first principles requires the use of systematic methods like those in Appendix A. See Table I for the representations under the point group.Mixing of these polynomials under translations is depicted in Fig. 3.
Here, we find a further subtlety in dealing with degenerate representations.The four independent cubic polynomials sort into two copies of the two-dimensional representation  1 of D 4 .Under translation, the quartic polynomial  −− irrep  3  +  3  mixes with the cubic  2 } span all cubic polynomial shift symmetries, so it seems natural to make a tree structure like Fig. 2 using these two  1 irreps as a basis for the cubic polynomial shift symmetries.
The conclusion here is that when there are degenerate irreps of a given degree, there is not generally a "canonical" choice of basis for those irreps for which a tree structure like Fig. 2 is unambiguous.For a given set of polynomials that we wish to include in M, compatibility with the space group symmetry requires us to find the "tree" of descendants and include them in M, which is a well-posed problem that can be solved algorithmically.However, that descendant information cannot be neatly encoded for all choices of polynomials into a tree of the sort shown in Fig. 2, because different "arrows" in the tree may require different basis choices.One can check that, in the present case, it is not possible to build such a tree, which is why Fig. 3 contains two equivalent bases shown in the shaded circles.Note, for example, that if both  4 +  4 and  3  +  3  are in the multipole group, then the former forces us to include { 3 ,  3 } and the latter forces us to include { 3 + 3 2 ,  3 + 3 2 }, which do span all cubic polynomials.However, including only one of these quartic polynomials in the multipole group only requires us to include one  1 irrep of cubic polynomials, not all cubic polynomials.

Triangular
The final Bravais lattice to consider is the triangular lattice.The point group is D 6 [Fig.1(e)] and the space group is 6.There is no additional subtlety compared to the square lattice.See Table II for the results, along with Fig. 4 for the translation mixing tree.a 1 a 2 FIG. 5. Kagome lattice with space group 6., , and  sublattices are labeled in red, blue, and green, respectively, but the sites are identical.The generating mirror, which is chosen to send  ↦ → −, is denoted by the double gray line.The lattice translation vectors a 1 and a 2 are denoted by the black arrows.

C. Beyond Bravais: Wallpaper groups
In the presence of a basis for the lattice, several things change compared to a pure Bravais lattice.First, the space group can be any of the 17 wallpaper groups.Second, space group symmetries can in general permute basis sites.Finally, polynomial shift symmetries do not in general need to act in the same way on each basis site.The interaction of all of these features leads to significant changes in the possible multipole groups.We will now give a few physically interesting or physically motivated examples which illustrate some key points about multipole groups.In particular, we give examples to illustrate the following facts: 1.The set of allowed multipole groups does not depend only on the space group, and instead depends on the details of the action of the space group on the lattice basis.Using the concept of Wyckoff positions, this action can be classified; we will give a procedure to generate this classification in principle, but will not perform the explicit classification for all 17 wallpaper groups.
2. In the presence of a basis, some multipolar symmetry groups may look very unnatural (particularly those generated by inhomogeneous polynomials), but are actually well-motivated in certain physical contexts.
3. The extended point group may in general be larger than the point group, specifically when the wallpaper group is non-symmorphic.

Multipole groups are not a function of space group alone
The kagome and triangular lattices both have space group 62.One might reasonably ask if the allowed multipole groups on the kagome and triangular lattices are "the same." The naïve answer is immediately "no," simply because there are more polynomial shift symmetries for the kagome lattice than the triangular lattice; on the kagome lattice, one can choose an independent polynomial shift symmetry on each of the three sublattices, whereas there is only one choice of polynomial on the triangular lattice.However, the same argument holds for three layers of the triangular lattice, and it is obvious that one can simply choose an allowed triangular lattice multipole group on each of the layers independently (with the space group operations acting simultaneously on all layers).We claim that the kagome lattice does not have this property; there are allowed multipole groups that are fundamentally distinct from all triangular lattice multipole groups.To show this, we classify multipole groups on the kagome lattice.
The wallpaper group 6 is generated by four operations: a six-fold rotation  (we are slightly overloading the letter , which also labels the  site in the lattice basis; the meaning should be clear from context), a reflection , and two translations  1 and  2 .With the sublattices labeled as in Fig. 5, space group operations act as where the six-fold-rotation and reflection matrices,   and   , respectively, and the two translation vectors, are given by with  the lattice constant.The point group is D 6 , the symmetries of a regular hexagon.The polynomial shift symmetry irreps of the point group are given in Table III, and the translation dependences are shown in Fig. 6.
Indeed, every triangular lattice multipole group has a corresponding kagome lattice multipole group, which acts identically on all three sublattices.However, even constant polynomials show new features.There are constant polynomials that form a 2D representation of the point group, which cannot happen on the triangular lattice or several copies of the triangular lattice.There are therefore only four possible multipole groups with only constant polynomial shift symmetries on the kagome lattice; there are independent choices of whether to include each of the two irreps, but those are the only possibilities.Compare this to three copies of the triangular lattice, where any constant shift  = (, , )  is a polynomial shift symmetry consistent with the space group.Physically, this is very interesting; including {  ,   } but not  0 means that conserving a vector charge but not the total scalar charge is consistent with the point group For the breathing kagome lattice (for which the point group is D 3 ), the irrep labels are modified according to Eq. ( 16), but the polynomials and their descendants remain unmodified.

Polynomial(s) Irrep Descendants
Polynomial(s) Irrep Descendants Mixing of irreps by translational symmetry for the kagome (and breathing kagome) lattice(s).Every arrow means that some (infinitesimal) translation acting on the polynomial(s) of higher degree generates the polynomial(s) of lower degree.Unlike previous trees, there exist two disjoint hierarchical structures.Physically, this means that there exist valid multipole groups that do not include total charge.The shorthand is used for the vectors that span the three-dimensional space that corresponds to the sublattice index.The representations of the polynomials under the point group are given in Table III.symmetry.We discuss the physical implications of this fact further in Sec.IV D. Similarly, the "dependency tree" in Fig. 6 does not decouple into multiple copies of the triangular lattice tree.
We conclude that, while the kagome lattice admits multipole groups identical to those of the triangular Bravais lattice, it also allows multipole groups that are fundamentally different from any multipole group on the Bravais lattice.The allowed multipole groups are therefore not simply a function of the space group.

A classification procedure
We saw above that the way that space group symmetries permute different basis sites in the lattice can substantially change the structure of the multipole group, and so the space group alone does not determine the allowed multipole groups.For a given space group, any set of basis sites can be decomposed into sets of independent, decoupled symmetry orbits; each orbit type is called a Wyckoff position, and the number of atoms in the orbit is called the multiplicity of the Wyckoff position.Wyckoff positions have been classified exhaustively for all space groups in both 2D and 3D [47][48][49][50][51][52].Since space group operations never mix lattice sites with different Wyckoff positions, one can independently choose allowed polynomial shift symmetries for each Wyckoff position.
Given a crystalline lattice, then, we can classify all possible multipole groups as follows.First, identify the space group of the crystalline lattice.Then, decompose the set of atoms in the unit cell into distinct Wyckoff positions.Next, for each Wyckoff position, compute the allowed multipole groups as we did above.A multipole group is generated by the space group operations and various polynomial shift symmetries that transform as direct sums of irreps of the extended point group, one summand per Wyckoff position.If a polynomial shift symmetry appears in M, so must its "descendants" under translations, which will be direct sums of the descendants of each individual irrep summand.
Given a Wyckoff position, the problem of sorting polynomial shift symmetries into irreps of the extended point group can be broken down further.First, one can determine the irrep of the extended point group formed by the permutation action of the space group on the basis sites.This is purely geometric information determined completely by the Wyckoff position.Second, one can take pure polynomials, without any reference to the Wyckoff position, and sort them into irreps of the extended point group.Finally, the irreps obtained in the prior two steps are combined using the Clebsch-Gordon coefficients of the extended point group; see Appendix A for details.
As an example of the above procedure, we can consider the classification problem for the wallpaper group .Assuming, without loss of generality, that the point group D 1 Z 2 is generated by reflections  ↔ −, there are three Wyckoff positions, shown in Fig. 7a.Wyckoff position  is an atom placed at (0, ) for any .Wyckoff position  is an atom at (1/2, ) for any .Wyckoff position  consists of a pair of atoms at (±, ) for any  (here,  is interpreted modulo the lattice constant in the  direction).
For atoms on either the  and  Wyckoff positions, polynomial shift symmetries do not permute basis sites within the unit cell.The multipole group classification for these atoms can be read off by inspection.The monomial     forms a trivial (resp.nontrivial) irrep of Z 2 if  is even (resp.odd), and translations mix irreps in the same way as Fig. 2.
For atoms on a  Wyckoff position, the point group operation  interchanges the two basis elements, that is, It is straightforward to see that irreps of the point group are given by transformations with  0 = (1, 1)  and  1 = (1, −1)  .The irrep is trivial if  +  is even and nontrivial if  +  is odd.Note that the permutation action on the basis sites is a direct sum of a trivial irrep  0 and a nontrivial irrep  1 , so the polynomial shift symmetries are just the (tensor) product of one of these   and the irrep formed by polynomials     in the spatial coordinates.The corresponding Clebsch-Gordon coefficients are trivial because the irreps involved are all one-dimensional.
For a generic arrangement of atoms with space group , then, we choose a set of polynomial shift symmetries that we wish to include.Each shift symmetry should be decomposed into a linear combination of direct sums of irreps of the extended point group, with one (possibly zero) summand in the direct sum for each Wyckoff position in the lattice.The complete symmetry orbit of the original polynomial shift symmetries then follows from the symmetry orbit and translation mixing of each direct summand.
For example, suppose we have atoms at Wyckoff position  (call the corresponding field  () ) and Wyckoff position  (call the corresponding fields  ()  , ).We could choose to include the shift symmetry that transforms the  Wyckoff position by  2  (trivial irrep of Z 2 ) and transforms the  Wyckoff position by  1 (nontrivial irrep of Z 2 ), where  is some constant with dimensions of length 2 : Under a mirror (as shown in Fig. 7a) since the first entry forms a trivial irrep of Z 2 and the last two entries form a nontrivial irrep of Z 2 .Hence, the multipole group must also contain ( 2 , −, )  .We know the translation descendants of each irrep ( 2  on the  Wyckoff position descends to , , , 1 and  1 on the  Wyckoff position descends to  1 ), but we must translate both irreps simultaneously to get the descendants of the combination.Namely, ( 2 , ±, ∓)  are descendants, but ( 2 , 0, 0)  is not a descendant.
We comment that the difference between the triangular, multi-layer triangular, honeycomb, and kagome lattices are exactly which Wyckoff positions are occupied; all of them have space group 6.According to the nomenclature in the Bilbao crystallographic server [48][49][50], the triangular lattice consists of an atom at the  Wyckoff position, and the multi-layer triangular lattice consists of multiple atoms at the  Wyckoff position.The honeycomb lattice has atoms in the  Wyckoff position (which has multiplicity two), and the kagome lattice has atoms in the  Wyckoff position (multiplicity three).Space group 6 also has two multiplicity-six Wyckoff positions, occupying one of which would generate the ruby lattice, and a multiplicity-12 Wyckoff position; the classification procedure for these Wyckoff positions using the techniques developed herein is straightforward but tedious.The , , and  sublattices are labeled in red, blue, and green, respectively, but the sites are identical.The generating mirror, which is chosen to send  ↦ → −, is denoted by the double gray line.The lattice translation vectors a 1 and a 2 are denoted by the black arrows.The C 3 rotation center to which each basis site is associated is depicted by a solid black circle.The zoomed circular inset illustrates the three vectors   that connect Bravais sites to basis sites.

Scope of the classification and physical interpretation of symmetries
There is an important subtlety in the notation when we write Eq. ( 2).In the absence of basis sites, there is no  index, and in the continuum, it is physically sensible to restrict our attention to continuous  (r).On the lattice, however, there is no a priori physical reason to restrict to continuous   (r).Rather than writing something like Eq. ( 2), we could instead consider shift symmetries where the field in question is shifted by an arbitrary number chosen independently at each basis site.This is the most general possibility, but attempting to classify such functions in any useful way is beyond the scope of this paper.
One natural way to restrict our attention to polynomial shift symmetries is to imagine defining a polynomial shift symmetry  (r) on the entire plane (without reference to basis sites), and then defining a shift symmetry on the lattice by restricting the domain of  (r) to r that belong to the crystalline lattice.This is a completely reasonable possibility, and it is included in our formalism.However, one could also define one polynomial   (r) per basis site, and then defining a shift symmetry on the lattice by restricting the domain of   (r) to those r that are an th basis element.The latter is the approach we take in this paper because it is more general; for a generic lattice with two basis elements  and , one could not obtain, e.g.,   (r) =  and   (r) =  by restricting the same finite-order polynomial to both the  and  sublattices.This subtlety in notation can sometimes obscure the physical meaning of certain multipolar symmetries.As an example which will become relevant in Sec.IV D, we consider the classification of multipole groups on the breathing kagome lattice, shown in Fig. 8.The breathing kagome lattice has space group 31, with point group D 3 .The allowed multipole groups for the breathing kagome lattice are almost identical to those for the kagome lattice.One can check that the table for breathing kagome can be obtained from Table III by simply replacing the irrep labels: with identical translation mixing.The only difference in the allowed multipole groups is that some irreps that are nondegenerate on the kagome lattice become degenerate on the breathing kagome lattice, which allows us to take linear combinations of these newly degenerate irreps without including both individual irreps.For example, ( 2 +  2 ) 0 transforms as  ++ on both the kagome and breathing kagome lattice, while   +   transforms as  −+ on kagome and  ++ on breathing kagome.On the kagome lattice, then, including ( 2 +  2 ) 0 + (  +   ) forces both ( 2 +  2 ) 0 and (  +   ) to appear in the multipole group, along with their descendants {, } 0 ,  0 , and {  ,   }.However, on the breathing kagome lattice, ( 2 +  2 ) 0 + (  +   ) transforms trivially under the point group, so only its descendants {2 0 +   , 2 0 +   } and  0 need to appear in the multipole group.Now consider the following polynomial shift symmetry: where   is the displacement of the th basis site from an associated Bravais lattice site, as shown in Fig. 8, and the   transform as vectors under symmetry operations.This shift symmetry is not, a priori, obviously physically meaningful.However, by inspection we can see that this amounts to shifting each field by R • x, where R is the closest Bravais lattice site.Physically, the presence of this symmetry means that the  component of the total dipole moment is conserved, where each field's charge is weighted by the closest Bravais lattice site's position instead of the physical position of the charge itself.One situation where this is physically meaningful is as follows.Imagine placing a spin-3/2  -like electron at the center of each small triangle and a localized  electron at each vertex of the breathing kagome lattice.With a strong Ising-like interaction between the  and  electrons of the sort on each triangle, one could imagine that the  electron's spin is equal to the sum of the surrounding  electron spins.When allowing  electrons to weakly interact between triangles, the symmetry Eq. ( 17) means that (a component of) the dipole moment associated to the  electrons is conserved.

Nonsymmorphic symmetries
When the space group is nonsymmorphic, namely when the space group is not a direct product of a point group and translations, the extended point group becomes distinct from the point group.
As an example, consider the space group , generated by translations and a glide consisting of reflection about the -axis (without loss of generality) combined with a half-translation.An example lattice and its symmetries is given in Fig. 7b, where the two types of lattice site are assumed to transform into each other under reflections about the -axis only.
Using the same formalism and applying a glide, we obtain where  = diag(1, −1).Now there is a nontrivial matrix acting on the coordinates rather than simply a translation.The same formalism as earlier applies, but now instead of finding irreps of the point group, we should find irreps of the extended point group Z 2 , which is generated by the transformation This extended point group is isomorphic as a group to / where  is the space group and  ⊂  is the subgroup consisting of pure translations.It is not hard to check that the irreps of the extended point group are       , where   = (1, (−1)  )  with  = 0, 1.The irrep is trivial if  +  is even and nontrivial if  +  is odd.After incorporating translations, it is straightforward to check that including       in the multipole group forces us to include   ′   ′   for any 0 ≤  ′ < , 0 ≤  ′ ≤ .These irreps form two decoupled copies of the tree in Fig. 2, one for each   .The formalism we have developed is therefore capable of constructing the allowed multipole groups to any desired order on any lattice.

III. CONSTRUCTING MODELS
Given a crystalline multipole group M, one would generally wish to construct models, both in the continuum and on the lattice, with M symmetry in order to study the effect of these exotic symmetries on interesting properties like their dynamics.In fact, as we will describe shortly, there is reason to expect that some multipolar symmetries, particularly when the multipole group is sub-maximal, can lead to highly interesting consequences like a sub-extensive number of emergent conserved quantities at long wavelengths.Such models are also of interest to guide the search for emergent multipolar symmetries in real materials.In this section we explain how one may construct effective Hamiltonians and field theories given a multipolar group M.

A. Spin models
We can also view the polynomials that comprise the multipole groups M that we have constructed as corresponding to conserved multipole moments of a local charge density.The construction of the multipole group ensures that we are able to write down Hamiltonians that are invariant under the space group of the lattice and conserve these multipole moments.Specifically, given spin- degrees of freedom living on the sites of a lattice L, we would like to construct tranlation invariant Hamiltonians that conserve the multipole moments of the local "charge density" Ŝ  for all polynomials  (r) belonging to the multipole group M. As we discussed in Sec.II C 3, in the presence of a basis, the weighting function  (r) corresponds to an in-principle independent multipole moment on each basis site.We now describe a systematic way to construct local Hamiltonians that conserve the moments (21) (with further details provided in Appendix D).The Hamiltonians are composed of local "gates" ĥx and their Hermitian conjugates.Gates are labeled by the index  (there may be more than one type of gate compatible with the imposed conservation laws) and are centered on x (which may or may not coincide with one of the basis sites).Hamiltonians built from such gates take the general form where the   are coupling constants.The gates ĥx are local, acting on spins belonging to a 'cluster' of spins  ⊂ L by incrementing or decrementing the  component of the spins by integers   () in the vicinity of position x: We will work exclusively with clusters of strictly finite support.This generic gate structure, or a specific variant thereof, has previously been utilized in Refs.[14,30,36,42,43] to construct Hamiltonians that conserve various moments of Ŝ  .When evaluating the time evolution of the charges (21) through their Heisenberg equation of motion, , Ĥ], the coefficients   () that define the gate effect a discrete derivative   acting on the function  .Importantly, the charge Q[  ] will be a conserved quantity if the derivative   annihilates  for all x, with this property being preserved under arbitrary perturbations (or operator insertions in ĥx , see Appendix D) that are diagonal in the Ŝ  basis.Hence, we can systematically construct Hamiltonians of the form (22) that conserve a finite list of multipole moments  (r) that form a multipole group M by identifying the possible discrete derivatives, specified by {  (  )}, that annihilate all  (r), i.e.,    = 0. Producing a Hamiltonian that is invariant under the space group from a given discrete derivative requires an extra step: One must find the orbit of the discrete derivative under the action of the space group and include all such operators as "gates" in Eq. ( 22) with equal weight   (all such operators generated in this manner from a valid discrete derivative are guaranteed to annihilate all polynomials belonging to the multipole group).In this way, space group operations will then just permute the gates, leaving Ĥ unchanged.

Constructing discrete derivatives
For a given multipole group, there are two handles that we can use to constrain the search for possible discrete derivatives: (i) the local Hilbert space dimension through the spin, , and (ii) the size (and shape) of the cluster .The size of the spin directly constrains the maximum absolute integer change in Ŝ  to be ≤ 2, while reducing the size of the cluster  reduces the dimension of the parameter space for the search; the total number of gates within a region of size | | is (2 + 1) | | .
For convenience, we arrange the integers that define the gate into a 'vector' |  ), where 0 ≤ |  | ≤ 2.In this language, a judicious choice of origin3 allows us to write the action of the discrete derivative as while the group status of M ensures that (   ) (x) will vanish for all x.Equation ( 24) therefore states that any gate that is compatible with all conservation laws, specified by the integer vector |  ), must be orthogonal to all vectors |  ) ∈ R | | corresponding to polynomials that belong to the multipole group.To touch base with continuum derivative notation natural for field theoretic treatments, we can find the continuum derivative to which   coarse grains by writing where, in the second line, we performed a Taylor expansion of the infinitely differentiable function (x) and replaced the sum over sites in the cluster with an inner product according to Eq. (24).Hence, to find the overlap of   with continuum derivatives of the form       , one simply has to evaluate the overlap between the vector |  ) that defines the discrete derivative and the vectorized monomial     (weighted by the appropriate combinatorial factor).Note that the reverse process -finding discrete derivatives from continuum derivatives -will not in general produce valid discrete derivatives that annihilate all  .This failure can occur when the multipole group is sub-maximal and a derivative of order , strictly less than the maximum polynomial degree , annihilates all polynomials belonging to M. Requiring that the discrete derivative coarse grains to the order- derivative does not constrain its overlap with derivatives of order ℓ > .This spurious overlap with derivatives of order ℓ can then prevent the order-ℓ polynomials from being annihilated by the discrete derivative.Indeed, this subtlety will be responsible for the interesting physical consequences that we discuss in Sec.IV.One should therefore always solve Eq. ( 24) in order to find the correct discrete derivative operators.

IV. PHYSICAL CONSEQUENCES
We have explained how the possible multipolar groups consistent with lattice symmetries may be constructed on arbitrary lattices.We now illustrate the power of the above formalism by exploring some physical consequences of multipolar symmetries.To this end, we pick one striking phenomenon that has previously been discussed, and generalize it to arbitrary crystal lattices, and also identify two remarkable phenomena that have not previously been discussed, as far as we are aware, but which we encounter when we explore multipolar groups on triangular and breathing kagome lattices respectively.A common theme in many of these discussions is the emergence of a (generally sub-extensive) set of conserved quantities at sufficiently long wavelengths for certain sub-maximal multipole groups.We work in the continuum to motivate the existence of these conserved quantities, and then identify examples with analogous behavior on lattices that host degrees of freedom with strictly finite local Hilbert space dimensions, subject to locality requirements.

A. General considerations -emergent symmetries
Given a scalar field  (for simplicity) subject to some polynomial shift symmetries, the field  can only appear in a symmetry-respecting Lagrangian as D  , where D  is some derivative operator that obeys for all polynomials  (x) in the multipole group.The notation D  is reserved for continuum derivatives, in order to distinguish them from their discrete counterparts   .Consider a multipole group whose highest-degree polynomials are degree ; then any D  of order ( + 1) or higher will solve the above equation.However, if the multipole group is sub-maximal, there may exist some D  of order  ≤  that annihilate all polynomials belonging to the group.As we will see, in some cases it may happen that the lowest-order D  annihilate additional polynomials (x) that are not in the multipole group, although higher-order solutions will generally not annihilate the additional polynomials.In this case, the additional polynomials lead to additional quantities where (x) is the density operator of the scalar field, which are (emergently) conserved at long wavelengths.One of the simplest examples of this phenomenon is a theory on the square lattice (see Sec. II B 3) that exhibits the following polynomial shift symmetries (see Ref. [43]): This is a valid multipole group; the polynomials  2 −  2 and 2 form distinct one-dimensional irreps of the point group D 4 ( −+ and  −− , respectively), and translation invariance requires that charge and both components of dipole are also conserved (see the hierarchy in Fig. 3).While this list of functions is annihilated by any third-order generalized derivative, all functions are also annihilated by the two-dimensional Laplacian ∇ 2 =  2  +  2  , which is second order in derivatives and the only such operator.However, we may immediately observe that any harmonic function (x), i.e., ∇ 2  = 0, will also be annihilated by the lowest-order generalized derivative ∇ 2 .Each such harmonic function  defines a quasi-conserved quantity via Eq.( 27).This infinite family of conservation laws is broken by higherorder, dangerously irrelevant derivative corrections that cease to annihilate the harmonic functions.

B. Localization in discrete Laplacian models
It was recently shown in Ref. [30] that if the dynamics on a generic lattice is given entirely by certain "discrete Laplacian" operators4, then the system exhibits strong fragmentation leading to localization [13,14].In our language, such gates are of the form   2 (0) =  and   2 () = −1 for all vectors  that connect a site to its  nearest neighbors (and the gate with the signs of   2 flipped).We will see examples of these gates momentarily.Using the formalism developed thus far, we are able to write down minimal valid multipole groups for which discrete Laplacians are the smallest allowed gates.Therefore, if the gate sizes are restricted to be the size of the discrete Laplacian and smaller, these symmetries are sufficient to enforce strong shattering/fragmentation in two dimensions.We illustrate this by explicitly working out the minimal necessary multipole group on the square and triangular lattices, although our formalism could obviously be extended to obtain the minimal sufficient multipole group on any lattice, in arbitrary dimensions.

Square lattice
For the square lattice, the minimal set of conserved multipole moments that is required to give rise to localization and is compatible with space group symmetry is: which produces the discrete Laplacian gate as the unique solution with smallest range.In (30), the color of a site denotes the sign of (  ), and the integer determines |(  )|, the combination of which determines the gate through Eq. ( 23).Since  and  2 −  2 transform as one-dimensional irreps of D 4 , we are able to remove just one of them whilst maintaining a valid multipole group.However, removing  allows gates that correspond to discretizations of  2  and  2  separately.On the other hand, removing  2 −  2 permits lattice discretizations of     , which can be discretized on the sites that surround a single plaquette, permitting operators with smaller range.The linear-order polynomials are not required to eliminate the undesired gates, but are required by translation symmetry.

Triangular lattice
On the triangular lattice, the set of multipole moments ( 29) is insufficient to fully constrain the gate of smallest range to be uniquely determined; instead, there are two D 3 -symmetric solutions with (0) = /2.This can be remedied by including an additional polynomial in the multipole group: which leads to the desired discrete Laplacian gate as the unique solution.The polynomial  3 − 3 2  can also be included in the multipole group (31) without changing the fundamental solution (32).If  3 − 3 2  is instead the only third-order polynomial in the multipole group, then the allowed gates are no longer unique; the same two solutions are permitted as for the set of multipole moments in Eq. (29).

C. Subsystem symmetries in 2D from 𝑶 (1) symmetries
We now discuss the first of the apparently new phenomena that we discover by applying our formalism, in this case to particular multipole groups on the triangular lattice.In particular, we find that, for certain choices of the multipole group M, all gates up to a certain size -parametrized by discrete derivatives -will additionally conserve (nonorthogonal) subsystem symmetries along lines of the triangular lattice.In this way, if we restrict the support of the possible gates, a finite list of conserved multipole moments on the triangular lattice will give rise to a much stronger emergent subsystem symmetry, involving an  () number of emergent conserved quantities for a lattice of linear dimension .
We begin by considering the fourth-order polynomial  (x) = ( 2 +  2 ) 2 , which transforms according to the trivial representation,  ++ , of the point group D 6 (see Table II).The 'descendant' polynomials of lower order  < 4 that must additionally be included for M to be closed can be found in Fig. 4. We further suppose that the third-order polynomial  (x) =  3 − 3 2 , which again transforms according to a one-dimensional irrep ( −+ ), is included in M. Note that the addition of this extra polynomial does not require any new descendants beyond those already included.Hence, our multipole group can be summarized as Despite the highest-order polynomial being of degree four, we are able to find a unique third-order derivative that annihilates all polynomials: D 1 =  3  − 3   2  .If we had not included  (x) =  3 − 3 2  in M, we would additionally be able to introduce a second third-order derivative, D 2 =  3  − 3 2    .We note in passing that the operator D 1 also naturally appears in the context of hydrodynamics in the presence of the point group D 3 , since it may be written compactly as D 1 =           , with     the third-order D 3 -invariant tensor [53,54].
At the level of continuum derivatives, the operator D 1 on its own annihilates a much larger family of functions than those belonging to M. This can be illustrated by acting on functions  k•x , we observe that D 1 acts multiplicatively as That is, any linear combination of oscillatory functions that satisfy   = 0 or   = ± √ 3  (i.e., those that do not vary parallel to the triangular lattice directions, e  or e  ∓ √ 3e  , respectively) will be annihilated by D 1 .This leads to the same conservation laws as a subsystem symmetry: since D 1 will annihilate () and analogous functions for the other two lattice directions, charge is conserved along every line that is parallel to each of the three lattice directions.

Lattice Hamiltonian
On the lattice, we can systematically construct discrete derivatives according to the procedure outlined in Sec.III A 1. We will restrict our attention to regions of the triangular lattice constructed by finding all sites contained within a circle of radius ℓ.The origin will be arbitrary, but circles centered on sites will generate D 6 -symmetric clusters while those centered on triangular plaquettes will generate D 3 -symmetric clusters, etc. Further, we will examine the most highly constrained case corresponding to minimal on-site Hilbert space dimension: spin-1/2 degrees of freedom.Subject to these constraints, the solution with the smallest range is unique and consists of spin raising and lowering operators acting around a hexagon (ℓ = 1 + , with  a positive infinitesimal): where the red sites (say) correspond to spin raising operators, Ŝ+ , and the blue sites to Ŝ− .We omit the integer labels utilized in the gates ( 30) and ( 32) since we are working with spin-1/2 degrees of freedom.Hence, if (x, 1), . . ., (x, 6) label the spins around a hexagon surrounding a lattice site x in a counterclockwise direction [the colored sites in ( 34)], we have the "ring-exchange" gate centered on sites, an interaction that is commonly found in the context of frustrated magnetism [55,56].The discrete derivative operator in Eq. ( 34), like the continuum derivative to which it coarse grains (D 1 =  3  − 3   2  ), conserves charge along lines.From (34), we observe that  ∈  1 (  ) = 0 for all lines  that are parallel to the three triangular lattice directions.Therefore, for every closed line connecting sites of the triangular lattice that is parallel to one of the lattice directions, there exists a corresponding conserved charge.For a lattice of  ×  primitive unit cells (see Fig. 1), there are 3 such conserved charges, although not all of these are independent.Hence, the behavior of the discrete derivatives mirrors that of the continuum derivatives: in the discrete (continuous) case, the operator that annihilates all polynomial moments with smallest range (with fewest derivatives) conserves charge along lines parallel to lattice directions.
We may now enlarge the radius of the circle to look for discrete derivatives that annihilate the moments in (33) with larger range.Clearly, if the enlarged region includes multiple hexagons of the form (34), we can form valid discrete derivatives by taking linear combinations of the motif in (34) (as long as the resulting coefficients all satisfy the constraint |  (r)| ≤ 1).For instance, the operator can be formed by taking a linear combination of solutions (34) on the central three sites forming an 'up' triangle, and will exhibit precisely the same conservation laws as (34).As a result, to find derivatives that do not preserve the constraint, we should restrict our search to operators that do not belong to the span of derivative operators of the form (34).The first such operator appears at the radius ℓ satisfying 2ℓ = √ 19 + : , which ceases to conserve charge along lattice directions that are not parallel to .Note that derivative operators belonging to the orbit of ( 36) under symmetry operations will also be valid solutions.Further note that we present the derivative operator with the smallest support by adding or subtracting operators of the form (34).The operator in (36) coarse grains to the fourthorder continuum derivative ∝   ( 3  − 3 2    ).While the exact conservation law is broken at the lattice scale by gates such as (36), in the long-wavelength limit, the subsystem-symmetrybreaking gates are less relevant and lead to an infrared (IR) description in which subsystem symmetry is broken only by higher-order, dangerously irrelevant corrections.
While we originally discovered this emergent subsystem symmetry on the triangular lattice, it turns out that an analogous scenario can also obtain on the square lattice.If we work on the square lattice, with spin-1/2 degrees of freedom, and conserve  2 ±  2 or { 2+1 ,  2+1 }, for even and odd orders, respectively, then these polynomials (and their descendants) are all annihilated by     , which has subsystem symmetry along lines.The gate of minimal size that is compatible with the chosen multipolar symmetries is then just ring exchange around a square plaquette, the consequences of which were discussed in detail in Ref. [24].Similar to the triangular lattice, there exist gates analogous to (36) that break the microscopic subsystem symmetry, but in principle we could always add additional multipolar conservation laws to forbid the smaller subsystem-symmetry-breaking gates, thereby ensuring that the smallest gate that breaks subsystem symmetry is larger than any desired size.

Haar-random circuits
To test our prediction that subsystem symmetry emerges at sufficiently long length and time scales, we perform simulations of Haar-random circuits that preserve the relevant conserved quantities (33).In particular, we work with Haar-random gates that are large enough for subsystem symmetry to be broken at the microscopic level, i.e., such that gates analogous to (36) are included.For two-point correlation functions, we perform the Haar average exactly, which gives rise to an effective stochastic automaton dynamics that can be simulated efficiently [scaling as poly()] for large systems (similar mappings exist for other quantities evaluated in Haar-random circuits, see, e.g., Refs.[57][58][59][60]).The mapping to automaton dynamics is outlined in detail in Appendix E. To summarize briefly, the Haar-random circuit is mapped to automaton dynamics that permits all symmetryallowed transitions (with equal probability) within a local region defined by the gate applied at each step.
As shown in Ref. [35], the Ward identity for charge conservation on the triangular lattice in the presence of subsystem symmetry along lattice directions is where the scalar current  is related to the charge density  via the constitutive relation  = − 1  2  3  at leading order, and the derivatives  1 ,  2 , and  3 are directed along the three triangular lattice directions.Equation (37) gives rise to the highly anisotropic decay rate Γ(k) =  6 cos 2 (3)/16, for density modulations of wave vector k, with (, ) the polar coordinates of k.That is, the decay rate Γ(k) has flat directions along  = /6 + /3 (arising directly from subsystem symmetry) and scales with the sixth power of .The Ward identity (37) straightforwardly determines the correlation function of Ŝ where the sum is over wavevectors k compatible with the periodic boundary conditions, and the overline denotes an average over circuit geometries and of the gates over the circular unitary ensemble (CUE), i.e., a Haar average.The function ( 38) is contrasted with the output of the Haar-random circuit simulations in Fig. 9, which exhibit excellent agreement.There is just one free parameter: the phenomenological subdiffusion constant .In the thermodynamic limit, we find from (38) that   (r; ) decays slowly with distance as |r| −1/2 for r parallel to one of the three triangular lattice directions, leading to distinctive sharp features that are reproduced by the numerical simulations.The coincidence of the theoretical prediction (38) and the random quantum circuit result verifies that the late-time behavior of correlation functions under generic dynamics that conserves the moments in ( 33) is governed by an equation of motion that exhibits subsystem symmetry.That is, even though there is generically no subsystem symmetry at the microscopic level, it nevertheless emerges at late times and long wavelengths if we conserve the  (1) list of multipole moments in Eq. ( 33).

Discussion
It is known that sub-maximal multipole groups can exhibit additional emergent conservation laws leading to unexpectedly slow dynamics controlled by dangerously irrelevant perturbations.The U(1) generalization of Haah's code is a central example of this [43].Our formalism can be used to construct many more models exhibiting such exotic physics on arbitrary lattices.We have illustrated this by a particularly striking construction, where imposition of a finite number of multipole moments on the triangular lattice leads to a robust emergent subsystem symmetry broken only by gates acting on twenty four or more sites (a similar construction can retrospectively be performed on the square lattice).If we limit the range of the gates such that the subsystem symmetry becomes exact (albeit accidental), then we get the exotic lattice consequences discussed in Ref. [24], including, for example, finite thickness 'shields' that disconnect the system.If we allow for terms in the Hamiltonian with arbitrary range, then the subsystem symmetry is eventually broken, but the subsystem symmetry breaking perturbations are irrelevant such that long wavelength hydrodynamic behavior is still consistent with subsystem symmetry (unless we specifically examine relaxation of subsystem symmetry charges, in which case it will be controlled by dangerously irrelevant perturbations).
The constructions we present herein are interesting for multiple reasons.Firstly, of course, there is the exotic quantum dynamics and hydrodynamics discussed above.Next, Haah's code is a particularly interesting example of a fracton phase [61][62][63][64][65], which continues to challenge several emerging paradigms (see, e.g., Refs.[66][67][68]).The U(1) generalization of Haah's code appears to inherit its exotic properties from a sub-maximal multipole group.We have provided a route to the construction of a multitude of sub-maximal multipole groups on arbitrary lattices.Going back from U(1) to Z 2 might then provide a whole family of models analogous to Haah's code, opening a new chapter for the field.Finally, subsystem symmetry itself is of interest [45,[69][70][71][72][73][74], but is extremely unlikely to arise exactly in a microscopic Hamiltonian.We have provided a general construction for how subsystem symmetry may be emergently obtained from a finite number of conservation laws, which may also be useful for uncovering subsystem symmetries in real materials.

D. Vector conserved quantities
In Sec.II C 3 we identified multipole groups for the breathing kagome lattice (Fig. 8) that did not require conservation of the  component of total magnetization  Ŝ  (i.e., did not require conservation of monopole charge).Here, we show that these multipolar conservation laws instead produce an emergent two-component vector conserved charge density, and that the smallest lattice derivatives additionally conserve moments of this vector density related to holomorphic complex functions.
Consider the following list of multipole moments, which, according to Table III, generate a valid multipole group on the breathing kagome lattice: where   ≡  −   , with   shorthand for the sublatticedependent vector shift that appears in Eq. ( 17) (analogously for   ).While the list of multipole moments in Eq. ( 39) looks somewhat unnatural, we are able to rewrite all position dependence of the polynomials in terms of the nearest Bravais lattice sites as where R corresponds to the position of the C 3 rotation center associated with the three sites (see Fig. 8), which are labeled by the index .That is, when evaluating multipole moments, sites are weighted according to the position of the rotation center to which they are associated.This point is discussed in further detail and motivated physically in Sec.II C 3. Unlike the examples considered thus far, this multipole group does not require conservation of total charge, which would correspond to conserving the moment specified by   =  0 = (1, 1, 1)  .In what follows, it will be convenient to introduce the following set of unit basis vectors {e () } for the triangular lattice Since the three vectors are related to one another by C 3 rotations, they satisfy  e () = 0. Note that this choice permits us to rewrite the conservation laws in terms of a two-component quantity   (R) via   (R) =  ()    (R).In terms of   (R), the conservation laws are simplified to Recall that the functions   (R) define conserved charges via Rewriting these conserved charges in terms of   (R), we have where we have introduced the operators ρ (R) ≡  ()  Ŝ R, , which live on the Bravais lattice site R.The conservation laws in (39) and ( 40) therefore imply that, while the total scalar charge Q Ŝ R, is not conserved, it is nevertheless possible to define an operator ρ (R) living on Bravais lattice sites with four conserved multipole moments: This follows from substituting (42) into (43).For convenience, we defined the two-dimensional cross product as the  component of the conventional three-dimensional cross product.Note that ρ (R) indeed transforms as a vector under point group operations.Specifically, C 3 rotations, under which ρ ↦ →  ℓ ρℓ , since three-fold rotations permute the spins associated to a given Bravais site ( ℓ is the rotation matrix for  = 2/3), and the generating mirror  ↦ → −, under which ρ ↦ → − ρ and ρ ↦ → ρ .
As shown in Sec.III A, we can construct lattice Hamiltonians that satisfy these conservation laws by finding (discrete) derivatives that annihilate the list of conserved functions in Eq. (39) [or, equivalently, Eq. ( 42)].While any second-order derivative will annihilate the polynomials in Eq. ( 42), there also exists a particular solution composed of first-order derivatives.Namely, the derivative operators (  , −  ) and (  ,   ) will also annihilate all moments   (R).However, these lowestorder derivatives additionally conserve a much larger family of functions beyond the finite list in Eq. (42).Any two-component function   (r) that satisfies will be annihilated by the first-order derivative operators that we have identified.These equations are simply the Cauchy-Riemann equations whose solutions define the holomorphic functions.Hence, any holomorphic function will be annihilated by the pair of first-order derivatives (  , −  ) and (  ,   ).In particular, the conserved first-order moments R R • ρ(R) and R R × ρ(R) appearing in Eq. ( 42) correspond to the holomorphic functions () =  and () = , respectively, where the  () component is recovered from the real (imaginary) part of the complex-valued function ().More generally, the first-order derivatives will annihilate all functions spanned by () =   and () =   for  ∈ N 0 .Such holomorphic conserved charges also arise in the context of hydrodynamics in the presence of a triangular point group when the current tensor transforms in the vector representation of D 3 [53].

Lattice Hamiltonian
To put the theory on a lattice comprised of spin-1/2 degrees of freedom, we are tasked with finding discrete derivatives defined by a set of integer coefficients   () satisfying |  ()| ≤ 1 that annihilate the list of functions in Eq. (39).Explicitly, for a cluster  composed of Bravais sites R and basis sites , we require that (   ) (x) =  ∈    (  )    (x+  ) = 0, where   is the displacement between the Bravais site associated with the site  and the center of the cluster x (see Appendix D for further details).As before, we will consider clusters defined by finding all sites contained within a circle of radius ℓ.Considering first a single 'down' triangle, we find the solution which adds or subtracts charge from the three sites associated to a given Bravais lattice site.The derivatives can also be expressed in terms of their action on discrete   (R) via (   ) (x) =  ∈,  ()    (  )  (x +   ), where the index  labels the Bravais lattice sites (we assume that the cluster  encompasses all  associated with each included Bravais site).We may therefore define the discrete derivative ( D  ) (x) =  ∈ ñ (  )  (x +   ), i.e., we define a set of coefficients ñ () ≡   ()     () which are associated to Bravais lattice sites (the coefficients ñ may no longer be integer valued).Note that the gate  0 in (46) leaves ρ unchanged since the basis vectors e () sum to zero (equivalently, D0 is the trivial derivative operator with ñ = 0).Next, we enlarge our search to include regions of that lattice that include up to three 'down' triangles.We find nontrivial solutions centered on both 'up' triangles and on hexagonal plaquettes The gates related to (47) by C 3 rotations are also valid solutions, which allows us to construct a Hamiltonian that is invariant under the space group.Up to these symmetry-equivalent solutions [and addition or subtraction of the solution in (46)], the gates in (47) are unique, although considering clusters of increasing size will eventually yield solutions that are linearly independent from (46) and (47).Upon coarse graining, the derivative operators D1 and D2 become precisely the operator (  , −  ) identified previously, i.e., D1 ∼ where the sites depicted now correspond to the triangular Bravais lattice formed by the centers of 'down' triangles.Their C 3 -rotated variants coarse grain to a linear combination of (  , −  ) and (  ,   ).Hence, as shown in Appendix D, the equations that define which functions are conserved by the Hamiltonian are precisely discrete versions of the Cauchy-Riemann equations ( 45), whose long-wavelength solutions will give rise to (approximate) holomorphic conserved charges.In a similar manner to Sec.IV C, even allowing for terms in the Hamiltonian with arbitrary range, the derivative operators that break holomorphic charge conservation are irrelevant.As a result, long-wavelength relaxation is determined by an equation of motion that preserves the holomorphic charges discussed herein.

Discussion
We have identified an exotic and (as far as we know) hitherto unknown possibility -a conserved multipole group that does not conserve monopole charge -and have identified a concrete realization on the breathing kagome lattice.We have pointed out that this multipole group conserves an emergent vector charge, plus all holomorphic functions of this vector charge.This constitutes a new and exotic class of problem, for which exploration of the thermodynamics and dynamics promises to be a particularly interesting challenge for future work.It also provides a 'proof of principle' of qualitatively new physics that is only accessible on non-hypercubic lattices, and thus can only be accessed through our formalism, or something analogous.

V. CONCLUSIONS
We have presented an algorithmic procedure by means of which one may construct all consistent conserved multipole groups to any desired order on arbitrary crystal lattices, and may further construct the minimal continuum field theory and lattice Hamiltonian consistent with said conservation laws.The procedure for constructing consistent multipole groups is as follows: given a space group, compute the extended point group and sort polynomials in the spatial coordinates into irreps.Then, for each Wyckoff position, decompose the permutation action of the extended point group on the basis sites into irreps.Using Clebsch-Gordon coefficients, combine the above irreps into irreps of the full extended point group.Finally, compute the translation mixing.This method is dimension-and latticeindependent, and can thus be applied to arbitrary crystal lattices.As such, it provides a complete in principle classification of the multipolar problem on arbitrary lattices.The procedure is labor intensive, so we have only carried it out for a representative set of lattices in two dimensions (all two dimensional Bravais lattices, plus kagome and breathing kagome).However, extension to, e.g., all 230 space groups in three spatial dimensions would be straightforward, albeit tedious.An explicit classification for all known crystal structures would be a worthwhile challenge for future work.
We have explored some interesting physical consequences using our construction.As a warmup, we have identified the minimal set of symmetries required to get localization from strong fragmentation on the square and triangular lattices respectively (using a method that could be generalized to any lattice).We have also identified two phenomena that do not appear to have been appreciated before.One involves an emergent subsystem symmetry arising from imposition of a finite number of multipolar conservation laws.The second new phenomenon has no known analog on hypercubic lattices whatsoever -it turns out that on the breathing kagome lattice, one can define a consistent multipole group that does not include monopole.Thus, one can write down consistent translation invariant Hamiltonians that conserve certain multipole moments of charge, but do not conserve total charge itself!Nevertheless, these models do conserve an emergent 'vector' charge, as well as all holomorphic functions thereof.This constitutes a striking and novel scenario, deserving of more detailed exploration in future work.
Our formalism opens up new directions for multiple lines of research.At the most basic level, it provides a guide to the construction of new fracton models on non-hypercubic lattices.5Particularly interesting, it allows us a way to identify consistent sub-maximal multipole groups.Given that the sub-maximal group appears to be the 'secret ingredient' that endows Haah's [U(1)] code with its uniquely exotic properties, it offers a route to the construction of a whole family of models analogous to Haah's code, on lattices other than cubic.Haah's code is the least well understood fracton model, challenging many paradigms [66][67][68], and construction of a family of models analogous to Haah's code might open up new directions for research into fracton phases.Given that Haah's code has interesting properties as a quantum memory [75,76], such models may also be useful for quantum information.On another front, subsystem symmetries are of considerable current interest [45,[69][70][71], but given that subsystem symmetries involve infinitely many conservation laws, it is not clear how such symmetries could arise in nature.We have provided constructions through which imposition of a finite number of multipolar conservation laws can lead to the emergence of subsystem symmetry, which could help guide the search for realizations of subsystem symmetries in nature.(It should be noted, however, that multipolar symmetries other than dipolar are already challenging to realize, as discussed in Ref. [13]).On a third front, our explorations have led us to discover new phenomena that have no known analog on (hyper)cubic lattices -for instance the possibility of having systems that do not conserve charge, but do conserve certain multipolar moments of charge.This opens a new direction for the study of fracton phases and associated phenomena on general lattices, and demonstrates that there is qualitatively new physics to be found.And finally, our work can guide the search for fracton physics and associated phenomena in real materials -many of which are not based on square or cubic lattices.Note that crisp experimental diagnostics for fracton phases have been identified in Refs.[77][78][79].
Finally, we have thus far limited ourselves to systems where all symmetries of the Hamiltonian are preserved.It would be very interesting to examine spontaneous symmetry breaking of crystalline multipole groups.Either the space group part or the polynomial shift part of the symmetry could be broken, and the interplay between discrete spatial symmetries and continuous polynomial shift symmetries could lead to interesting effects.One could even consider exotic possibilities like breaking a spatial symmetry and a polynomial shift symmetry but preserving the product.An exploration of such exotic symmetry breaking phenomena is an(other) interesting challenge for future work.

A.2. Sorting polynomials into irreps
There are two algorithmic ways to sort polynomial shift symmetries into irreps of the extended point group.The first method is a simple one to perform on a computer.Fix a Wyckoff position of multiplicity .By construction, each extended point group operation maps a collection of  homogeneous polynomials of degree  to another collection of homogeneous polynomials of degree , so we may restrict our attention to the case where all basis sites transform via a homogeneous polynomial of degree .The set of all such transformations is spanned by the transformations   (r) =  , 0    − for all  0 = 1, 2, . . .,  and  = 0, 1, . . ., .Viewing these transformations as a basis | 0 , ⟩ for the set of polynomial shift symmetries under consideration, it is easy to compute directly how extended point group operations act on these transformations.For example, for  = 2 on the honeycomb lattice, rotations take This gives us a matrix representation of each extended point group operation.Given such a matrix representation  () for  ∈ , where  is the group in question, and given the characters  ℓ () where ℓ labels a representation of , one can construct the projector  ℓ onto the representation ℓ via the formula The eigenvectors of this projector with eigenvalue 1 are a basis for the polynomials which transform under the given representation ℓ.One can then convert this into any convenient basis.
The second method is analytical and uses the Clebsch-Gordon coefficients of the extended point group.First, one can decompose the polynomials {, } into irreps of the pure coordinate transformations r → sr.All higher-order polyno-mials are formed as (tensor) products of some of these irreps with themselves; using Clebsch-Gordon coefficients, one can decompose those tensor products into irreps of these pure coordinate transformations.For example, under the extended point group D 4 , {, } forms the two-dimensional representation  1 .Hence quadratic polynomials all appear in the tensor product  1 ⊗  1 =  ++ ⊕  +− ⊕  −+ ⊕  −− .Using the Clebsch-Gordon coefficients, we find that the  ++ representation is  2 +  2 , the  −− representation is 2, the  −+ representation is  2 −  2 , and the  +− representation is 0 (and thus not present).This procedure can then be iterated; cubic polynomials are formed from (tensor) products of the linear polynomials and quadratic polynomials, and so on.A similar, related approach is to restrict polynomial representations of O(2) to representations of D  via "branching rules" [80,81].This is the complete classification procedure if the lattice is a Bravais lattice.In the presence of a basis, we notice that the permutation action of the extended point group on the fields, and in particular on each Wyckoff position, also produces a (not necessarily irreducible) representation of the extended point group.This representation is also the representation formed by constant polynomial shift symmetries (since the coordinate transformation does nothing to constant shift symmetries) and can be decomposed straightforwardly into irreps, potentially using the projector method as above.The overall action of the extended point group on the fields is the tensor product of this permutation action with the action of pure coordinate transformations.Hence the decomposition into irreps consists of choosing a permutation irrep (i.e., constant polynomial shift symmetry), tensoring with a coordinate transformation irrep (i.e., a polynomial), and using Clebsch-Gordon coefficients to decompose the tensor product into irreps.For example, on the honeycomb lattice the permutation action is represented as the  ++ ⊕  −− representation of D 6 , where the  ++ represents sublattice-even shift symmetries and  −− represents sublatticeodd symmetries.Hence, each irrep of the extended point group is just the tensor product of the irreps formed by pure polynomials (which are identical to those of the triangular lattice) with one of these permutation irreps.Since the permutation irreps are all 1D, the Clebsch-Gordon coefficients are very simple.

A.3. Clebsch-Gordon coefficients for D 𝑴
For convenience, we list the Clebsch-Gordon coefficients for D  in the reflection eigenbasis.One derivation is given in Ref. [82] in the rotation eigenbasis.We make a notation change for this appendix; instead of referring to the one-dimensional irreps of D  as  , , we will refer to them as  , for  = 0, /2, where  = 0 corresponds to  = +1 and  = /2 corresponds to  = −1.This notation makes several of the Clebsch-Gordon coefficients significantly simpler.
We first give the rules for which representations appear in the tensor product: text, namely where we formally take ℓ infinitesimal, generate a collection of polynomials of degree  − 1, and then repeat the process for the newly generated polynomials.
Without loss of generality, let a  be oriented along the -axis; for a general direction, all partial derivatives may be replaced by directional derivatives.The power series expansion of a polynomial  of finite degree yields The highest term  is the maximum power of  appearing in  , which is, of course, bounded from above by the degree of  .Observe that    /  always contains a term involving  −    (for some fixed power ) but never any terms involving a larger power of .Therefore, the  + 1 polynomials    /  are linearly independent as elements of the vector space of polynomials over R. We claim that this set of  + 1 polynomials spans the same subspace as the set of  (r + ℓa  ) for all integer ℓ.Equation (C.2) immediately shows that the span of the latter set is contained in the span of the former.To see the other way around, consider the  + 1 polynomials given by ℓ = 1, 2, . . .,  + 1. Writing these polynomials in the basis given by    /  , we obtain a set of vectors This is a Vandermonde matrix with nonzero determinant since no two of its rows are equal.Therefore, it is invertible; its inverse is therefore a basis transformation that expresses the    /  as a linear combination of the  (r + ℓa  ) with our chosen values of ℓ, so the span of the derivatives is contained in the span of the translates.In particular,   / is in the span of the translates.By symmetry,   / is as well.We summarize the above argument as follows: given a lattice translation a  , the directional derivative (a  • ∇)    for all  span the same space of polynomials as the translates   (r + ℓa  ) for all ℓ.It remains to show that given another lattice translation a ′  , the mixed derivatives (a  • ∇)  (a ′  • ∇)    spans the same space of polynomials as   (r + ℓa  + a ′  ).We can simply repeat the above argument, but applied to (a ′  • ∇)    .Then we can conclude that (a  • ∇)  (a ′  • ∇)    spans the same polynomials as (a ′  • ∇)    (r + ℓa  ).Applying the same argument again, we see that the latter spans the same polynomials as   (r + ℓa  + a ′  ), as desired.) || ( ∈ Z) for spin- degrees of freedom, we find that That is, the commutator between Q[  ] and the Hamiltonian in (D.1) composed of local gates effects a discrete derivative on  .Hence, any function  satisfying    = 0 defines a corresponding conserved charge   Q[  ] = 0.In particular, if the coefficients   () satisfy  ∈   (  ) = 0, then   will annihilate the unit function  (r) = 1 and total charge Q [1] will be conserved by dynamics generated by (D.1).Note that the arguments can also be applied in reverse: Given a Hamiltonian of the form (D.1), one can identify a family of conserved charges Q[  ] by solving the equations    = 0.

D.2. Introducing a basis
Now consider introducing a -spin basis at each site of the Bravais lattice labeled by an index  = 1, . . ., .The action of a discrete derivative is still (   ) (x) =  ∈   (  )  (x +   ), where   is now the displacement of the basis site  from the center of the cluster .However, we can alternatively index the derivative coefficients and the function  according to Bravais sites and a basis index.In this case the discrete derivative may be written (  () )    (x +   () ) , (D.5) where x +   () corresponds to the position of a Bravais site and  () labels the Bravais lattice site associated to site .That is, the vectors   are displacements of the Bravais sites from the center of the cluster (now defined as  ∈   , which is identical to  −1  ∈   () if  contains all basis sites associated with each Bravais site).Note that this is merely a reparametrization of the discrete function.Specifically, we are not assuming that the function depends only on the position of the Bravais site.For example, in the labeling scheme employed in (D.5), the function  (r) = r 2 would become   (R) = (R +   ) 2 , where   is the vector that connects Bravais site R to basis site r  .In this language, the putatively conserved quantities Q[  ] now depend on a function   (R) that itself depends on both the Bravais lattice site R and the index :

D.3. Additional symmetries
Note that the Hamiltonians in Eqs.(D.1) and (D.7) possess additional symmetries, as pointed out in, e.g., Refs.[14,30].For instance, the parity operator Π ≡     Ŝ  commutes with the Hamiltonians (D.1) and (D.7) since Π has the effect of interchanging Ŝ+  with Ŝ−  , i.e., Π Ŝ±  Π = Ŝ∓  , therefore interchanging the gates ĥx ↔ ĥ † x .We may, however, add any perturbation that is diagonal in Ŝ  basis to the Hamiltonian while maintaining the conserved operators (D.3) and (D.6), since this will not affect the commutators in Eqs.(D.4a) and (D.8).Because the discrete symmetry Π anticommutes with Ŝ  , it only remains a conserved quantity if the perturbation consists of an even number of Ŝ  operators.The anticommutation of Π and Ŝ  also implies that Π and the multipole moments Q[  ] anticommute.Hence, if Π does commute with the Hamiltonian, the sectors with quantum numbers { [  ]} and {− [  ]}, for all  that are conserved by Ĥ, will have precisely the same spectrum.

D.4. More general Hamiltonians
As noted in the previous subsection, the Hamiltonians (D.1) and (D.7) will conserve the same family of charges Q[  ] if they are subjected to arbitrary perturbations that are diagonal in the Ŝ  basis.In fact, the statement is more general: even the gates ĥx themselves can be 'decorated' by operator insertions that commute with Ŝ  .To illustrate this, consider a one-dimensional lattice (with no basis) that hosts a theory conserving only total charge  Ŝ  .The smallest gate of the form (D.2) that can be written down is simply ĥ,1 = Ŝ−  on the central site (physically, a hop of two units to the right cannot be exactly decomposed into two sequential hops, since charge on the central site may interfere with the sequential hopping process).Since Ŝ+

Ŝ−
commutes with Ŝ  , the operators ĥ−1,1 ĥ,1 and ĥ,2 effect exactly the same discrete derivative (∼   ) and therefore possess the same conserved quantities.As a result, the gates that we identify using the methods described in the main text are 'canonical' in the sense that they conserve all relevant  (x) and all such diagonal operator insertions are absent; more complex gates can then of course be constructed by introducing such diagonal operators.standard "brickwork" geometry of gates, the location of each gate is chosen at random from a uniform distribution over the lattice (and one unit of time is defined as an extensive number of such random gate applications).For a given cluster acting on a region ℓ, the unitary gate has the following structure  .The quantity that is averaged over in Eq. (E.3) is interpreted as the probability that the system transitions from state s to s ′ after evolving the system for a time : Note that, since each gate is drawn from an independent distribution, the CUE average has decomposed into two separate averages.Since the gate (E.1) acts as the identity outside of ℓ, we immediately observe that the matrix elements enforce that the states s 1 and s 2 must coincide in l.We now perform the average on the top line exactly using the one-fold Haar channel: for each block, Φ[ ] ≡ Û Â Û † =  −1 Tr[ Â]1 [85] (with  the dimension of the random matrix Û) where the subscripts ℓ and l denote the region of the lattice on which the inner products are evaluated, and   is the number of states in the symmetry block .Since the projectors can also be chosen to be diagonal in the Ŝ  basis, each term under the summation vanishes if |s 1 ⟩ ℓ or |s 2 ⟩ ℓ does not belonging to the block .If both belong to the block , then the two states must coincide.We can therefore eliminate s 2 from (E.The evolution from s → s ′ can therefore be decomposed into a sequence of such transfer matrices, each of which incorporates the effect of a gate application.The transfer matrix (E.8) is the state version of the operator transfer matrix derived in Ref. [60]; it checks which block the input state belongs to and sends it to a mixture of all other local states belonging to the same block with uniform probability determined by the size of the block: 1/  .The correlation function (E.2) can then be evaluated efficiently by performing a stochastic automaton time evolution, where Ŝ  eigenstates |s⟩ are sent to other eigenstates |s ′ ⟩ according to the transition probabilities determined by the transfer matrix (E.8).To make sure that the gate (36) (which microscopically breaks the subsystem symmetry) is included in the circuit, we work with clusters of sites of radius ℓ satisfying 2ℓ = √ 19 + , centered on bonds of the lattice.

FIG. 7 .
FIG.7.Left: An example of a lattice belonging to the space group .Axes of reflection are denoted by the double gray lines.The labels , , and  correspond to the three types of Wyckoff position, with multiplicities one, one, and two, respectively.Right: A lattice with space group .The dashed gray lines represent axes of glide reflection since the two types of lattice site are (only) mapped onto one another under the mirror  ↦ → − (supposing they were spatially coincident).The solid gray lines represent lattice translations, which are parallel to the and -axes.

a 1 a 2 FIG. 8 .
FIG.8.Breathing kagome lattice, which belongs to the space group 31.The , , and  sublattices are labeled in red, blue, and green, respectively, but the sites are identical.The generating mirror, which is chosen to send  ↦ → −, is denoted by the double gray line.The lattice translation vectors a 1 and a 2 are denoted by the black arrows.The C 3 rotation center to which each basis site is associated is depicted by a solid black circle.The zoomed circular inset illustrates the three vectors   that connect Bravais sites to basis sites.

FIG. 9 .
FIG. 9. Comparison between the autocorrelation function   (r; ) = ⟨ Ŝ (r; ) Ŝ (0; 0)⟩ obtained for a Haar-random circuit that conserves the multipole moments (33) (left) and the corresponding hydrodynamic prediction (38) (right).The correlation function is evaluated using an effective classical automaton evolution that allows us to reach large systems and times.The spatial profile is illustrated at a fixed time,  = 5 × 10 3 , for a system of linear size  = 512.

Ûℓ = 1
l ⊗    , (E.1) where l is the region of the lattice complementary to the gate region ℓ, and   is an   ×   random unitary matrix.The gate (E.1) is decomposed into blocks labeled by  according to their symmetry quantum numbers.That is, all states |,   ⟩ (where   = 1, . . .,   ) are eigenstates of the multipole charges Q[  ] for  belonging to the multipole group, Q[  ] |,   ⟩ =  [  ] |,   ⟩, and have the same eigenvalues  [  ] for all  .Since all Q[  ] are diagonal in the Ŝ  basis, we may take the decomposition (E.1) to be in the basis defined by product states of the form ⊗  |  ⟩, with |  ⟩ ∈ {|0⟩ , |1⟩} the eigenstates of Ŝ  on site , i.e., Ŝ |⟩ = (−1)  |⟩.We now show how Haar-averaged two-point correlation functions map onto stochastic automaton dynamics for operators that are diagonal in the Ŝ  basis.The derivation is similar to that of Ref.[60], except that we work with states rather than vectorized operators, which makes the correspondence with automaton dynamics more crisp.Consider an infinite temperature two-point correlation function of the form   () = Tr ρ Ŵ † () Ô Ŵ () Ô  , (E.2)where the time evolution operator Ŵ () is given by a product of microscopic random unitaries (E.1), and ρ = 1/ is the infinite temperature density matrix ( being the total dimension of the many-body Hilbert space).The overline denotes 'Haar averaging', i.e., averaging each block belonging to the gate (E.1) over the unitary group U (  ) with respect to the Haar measure.The corresponding ensemble of matrices is the circular unitary ensemble (CUE).If the operators are diagonal in the Ŝ  basis, it is convenient to evaluate the trace in this basis:    () = 1  ∑︁ s,s ′   (s ′ )  (s)⟨s ′ | Ŵ () |s⟩ ⟨s| Ŵ † () |s ′ ⟩ , (E.3) where s and s ′ represent eigenstates of Ŝ s ′ s () ≡ ⟨s ′ | Ŵ () |s⟩ 2 .(E.4)Suppose that each microscopic gate application is associated with a time  (i.e.,  2  = 1 defines one unit of time).If a gate was applied on the region ℓ to evolve the system from time  −  to , then, inserting two resolutions of identity, s ′ s () = ∑︁ s 1 ,s 2 ⟨s ′ | Ûℓ |s 1 ⟩ ⟨s 2 | Û † ℓ |s ′ ⟩× ⟨s 1 | Ŵ ( − ) |s⟩ ⟨s| Ŵ † ( − ) |s 2 ⟩ .(E.5)
and   , respectively), even if the discrete lattice translations are not orthogonal.Hence, including     requires us to include  −1   (if  > 0) and    −1 (if  > 0) in the multipole group.

TABLE I .
Polynomials up to degree  = 4 and their representations under the point group D 4 of the square lattice.The irrep labels  , and   are defined explicitly in Appendix A. The descendants of each polynomial under translation mixing are also listed.

TABLE II .
Polynomials up to degree  = 4 and their representations under the point group D 6 of the triangular lattice.We choose the generating mirror  to send  ↦ → −.The irrep labels  , and   are defined explicitly in Appendix A. The descendants of each polynomial under translation mixing are also listed.

TABLE III .
Polynomials up to degree  = 2 and their representations under the point group D 6 of the kagome lattice.We choose the generating mirror  to send  ↦ → −, and use the notation Consider a given discrete derivative   on a Bravais lattice defined by its action on discrete functions (   ) (x) =  ∈   (  )  (x +   ).The real coefficients   (  ) are asso-ciated with site , which is displaced by   from the center of the cluster of sites  (which we take to be  ∈   ).Note that x, the displacement of the center of the cluster, might not be centered on sites of the Bravais lattice (e.g., x may correspond to the plaquettes or the bonds of the lattice).The positions x+  always correspond to Bravais lattice sites, however.Suppose that the coefficients   () are all integer valued (perhaps by removing a common factor).If the integer-valued coefficients satisfy |  ()| ≤ 2, we can then define a Hamiltonian using the discrete derivative   acting on spin- degrees of freedom:The sign of the coefficients, sgn(  ()) ∈ {+, 0, −} (where sgn 0 = 0), determine whether the site  in the cluster is associated to a spin raising or lowering operator.Turning the problem around, given spin- degrees of freedom, the restriction | ()| ≤ 2 can impose strong constraints on the discrete derivatives that form valid Hamiltonians.In the most highly constrained case -spin-1/2 degrees of freedompermitted derivative operators must satisfy   ∈ {−1, 0, 1} only.The utility of the construction in (D.1) is that the Hamiltonian conserves the moments of "charge" defined by functions  that are annihilated by the discrete derivative   .Explicitly, consider the putatively conserved operator Q[  ], parameterized by the function  (r) Q[  ] = ∑︁ r  (r) Ŝ is a choice that is determined by the physics of the problem, but both cases are handled by the notation in (D.6).Suppose that the Hamiltonian now comprises multiple gates labeled by the integer    .Repeating the calculation that led to (D.4a), we find that the time evolution of the charge Q[  ] is determined by  Q ∝ ∑︁    (  () )    (x +   ()) ĥx − ĥ †  ] is a conserved quantity under dynamics generated by (D.7) if all discrete derivatives annihilate the function   (R).