Fragmentation-induced localization and boundary charges in dimensions two and above

We study higher dimensional models with symmetric correlated hoppings, which generalize a one-dimensional model introduced in the context of dipole-conserving dynamics. We prove rigorously that whenever the local configuration space takes its smallest non-trivial value, these models exhibit localized behavior due to fragmentation, in any dimension. For the same class of models, we then construct a hierarchy of conserved quantities that are power-law localized at the boundary of the system with increasing powers. Combining these with Mazur's bound, we prove that boundary correlations are infinitely long lived, even when the bulk is not localized. We use our results to construct quantum Hamiltonians that exhibit the analogues of strong zero modes in two and higher dimensions.

Fragmentation in particular, appears naturally in systems that conserve multipole moments [20,21].It manifests as an exponentially large number of dynamically disconnected sectors even after resolving the global conserved quantities.Such fragmentation (or, in the parlance of classical stochastic dynamics, reducibility [35]) comes in multiple flavors.One can, for example, distinguish weak and strong fragmentation, which give rise to distinct dynamical signatures [20,25].A particularly striking example of the latter was found in Refs.[20,21,36], where fragmentation was shown to result in localized dynamics that retains local memory of initial conditions.While such strong reducibility has long been known for kinetically constrained glassy models [35], the recently studied examples, originating from dipole-conserving systems, were restricted to one spatial dimension [25,36].
In this paper we introduce a set of models, which we name "discrete Laplacian models".Their defining feature is a correlated hopping term that distributes particles symmetrically among all neighboring sites, resembling a discrete second derivative and generalizing the 1D model of Refs.[20,21] to arbitrary lattices.We prove that the discrete Laplacian models all exhibit the same kind of localized behavior due to the strong fragmentation of their configuration spaces.This fragmentation originates from a combination of the correlated hopping of particles and a local constraint on the number of particles per site.Together these lead to a finite density of sites whose configuration remains frozen at all times, implying strong fragmentation.
Even away from the strongly fragmented limit, we find that these models exhibit localized behavior close to the boundaries of the system.We explain this by the presence of spatially modulated symmetries, a concept introduced in Ref. [37].We show the discrete Laplacian models posses conserved quantities that, related to solutions of the discrete Figure 1: Discrete Laplacian models.Family of models where the elementary process corresponds to a simultaneous hopping of particles from a site to all its neighbors.Illustrated here for different two-dimensional and three dimensional lattices.
Laplace equation, are localized near the boundaries.These include exponentially localized charges at the corners of the lattice, but also other quantities that are algebraically localized close to the boundary.In fact, we construct multiple families of such quantities, which decay with increasing powers of the distance away from the boundary.Just as in Ref. [37], we utilize Mazur's bound to show that these symmetries prevent the decay of correlations at the boundary, even when the bulk is not localized.
The idea of conserved quantities localized near the boundaries of a system bears resemblence to that of strong zero modes (SZM) in quantum spin chains [38].We will show that a modification of the discrete Laplacian models in the quantum setting can lead to phenomena analogous to SZM in two dimensions, with a system that exhibits degeneracies throughout its spectrum with open, but not with closed boundaries 1 .Nevertheless, our boundary modes differ in important ways from those of Ref. [38], as we discuss in detail below.
The remainder of the paper is organized as follows.In Sec. 2 we introduce a version of a discrete Laplacian model on a 2D square lattice and provide numerical results for the behavior of its bulk and boundary correlations, which will motivate our subsequent investigations.In Sec. 3, we first generalize the model to arbitrary lattices and graphs and then provide a proof of localization in the case when the number of local configurations is strongly restricted.In Sec. 4 we turn to the construction of conserved charges localized at the boundaries and discuss their implications for boundary correlations.Finally, in Sec. 5 we consider quantum versions of the discrete Laplacian models and we show how they can be modified to exhibit a phenomenology similar to strong zero modes.We conclude in Sec. 6.

Motivating example: Correlations in a 2D discrete Laplacian model
We begin by discussing numerical results for a particular 2D model.In subsequent sections, we will analytically explain the observed behavior, while also generalizing the model to arbitrary lattices.We consider models of classical stochastic dynamics, which allow for large-scale numerical simulations [16,17,37].The reason is that we mostly explore "classi-  (11).Inset: Spatial distribution along the xaxis for L = 250 and different times t.Unlike the weakly fragmented case, the correlation remains localized even for longer times.
cal" phenomena that relies entirely on symmetries that are common to both the classical and quantum cases.We discuss features specific to quantum Hamiltonians at the end of the paper, in Sec. 5. We consider a square lattice, where each site r hosts a classical spin s r taking values in s r ∈ {−S, −S + 1 . . ., +S} with (half)-integer S. We will equivalently refer to s r as the amount of "charge" on site r.The dynamics is generated by a set of local updates, or "gates", G + , acting on site r and its four neighbors.The effect of a gate is to change s r → s r − 4 and s r+δ → s r+δ + 1 where δ = ±e x , ±e y is any of the four elementary lattice vectors.In other words, it transfers one unit of charge from the central site r to all of its neighbors simultaneously (see Fig. 1a).This is well-defined for S ≥ 2, and it is applied only if it does not lead to a violation of the local constraint |s r | ≤ S. For future reference we summarize the effect of the gate G + by a set of integers as: such that s i,j → s i,j + n i,j under the effect of a gate applied at r = 0. We also consider the inverse gate, which is obtained by replacing each n i.j with −n i,j .
To turn these local updates into a stochastic evolution, in each time step we randomly tile the whole system with non-overlapping gates and apply the resulting layer (for additional details, see Ref. [37]).At each location, we choose randomly between the following three possibilities with equal probabilities: (i) apply the gate G + , (ii) apply its inverse, or (iii) no update is made.The defined model conserves several multipole moments of the charge, namely: the total charge, both components of the dipole moment, and the traceless part of the quadrupole moment.With open boundary condition, it also has a large number of additional conserved quantities, as we shall discuss below in Sec. 4.
We study the dynamics of the system by investigating the behavior of "infinite tem-perature" connected charge-charge correlations, which are defined as Here, we are uniformly sampling over all possible global initial configurations s(0) of the system, whose total number is given by D = (2S + 1) N for a system of N sites.s(t) is the time-evolved configuration corresponding to a particular trajectory of the system, and the brackets denote averaging over different circuit realizations.We use C(r, t) ≡ C r (r, t) to denote autocorrelations on the same site.The expected generic late-time behavior for a system with the above multipole symmetries is a subdiffusive decay of correlations of the form C(r, t) ∼ t −1/2 [14,18,41].This behavior is indeed observed in our model when S > 2 as Fig. 2a shows.However, Fig. 2a also shows that for S = 2, bulk correlations saturate to a finite value that remains finite even in the thermodynamic limit (see App.A for a careful finite-size scaling analysis).This is reminiscent of the behavior observed for S = 1 in the one-dimensional version of the same model [20], where it appeared as a result of Hilbert (or configuration in the classical model) space fragmentation.In Sec. 3 we argue that a similar explanation is valid here as well and it generalizes to a larger class of models, defined on arbitrary lattices and graphs.
While bulk correlations decay to zero for any S > 2 at infinite times, we find that correlations evaluated near the boundary of a system with open boundary conditions remain finite for any S, as shown in Fig. 2b.In Sec. 4, we will explain this by the presence of a set of additional conserved quantities localized at the edges of the system, similar to the behavior observed for certain 1D systems in Ref. [37].We construct a hierarchy of conserved quantities that are power-law localized towards the bulk and utilize Mazur's bound to show that these lead to the finite boundary correlations observed above.This is consistent with the fact that we do not only observe finite boundary auto-correlations, but that in fact, the spatial correlation along the boundary is also localized (see inset in Fig. 2b).

Localization from fragmentation
In this section, we first define a class of models, which we name "discrete Laplacian" models, that generalize both the 1D model studies in Ref. [20] and the 2D model discussed above in Sec. 2. We then go on to prove that all such models exhibit localized dynamics when the number of local configurations is sufficiently restricted.

General discrete Laplacian models
The model (1) introduced in the previous section can be thought of as a natural generalization of the 1D model studied in Refs.[20,36].In fact, as we now show, we can generalize this model to an arbitrary lattice in any spatial dimension, or even for an arbitrary graph.Consider a graph G(V, E) defined by vertices V and edges E. Let z v denote the degree (number of neighbors) of vertex v and we assume that the degrees are bounded: z v ≤ z 0 , ∀v ∈ V .We assign a spin variable s v = −S v , −S v + 1, . . ., +S v to each site, where we now allow the local spin to vary with the vertex v.In the proofs of Sec.3.2, we will find it easier to work with a shifted variable, m v = s v + S v , which therefore takes values 0, 1, . . ., M v ≡ 2S v .We will refer to this latter variable as the "particle number" on vertex v.
Let us generalize the model in Eq. ( 1) and define the following local gate acting on a vertex v and its neighbors N v , specified by the integers The effect of this gate is to remove z v particles from v and distribute them equally (one each) between its neighbors, leaving v completely empty when M v = z v .This is illustrated in Fig. 1 for a number of different graphs corresponding to 2D and 3D lattices.We also define the inverse gate, which takes one particle from each of the neighbors of v and collects them at v. As before, we only allow these gates to act if it does not lead to a violation of the constraint 0 ≤ m v ≤ M v for any vertex v.We refer to the set of models built from these elementary updates as discrete Laplacian models 2 .
While there are various ways in which these local updates rules can be turned into a stochastic evolution, the properties we discuss are largely independent of these details and follow directly from the constraints imposed by the structure of the gates (although we do assume detailed balance).We note that the discrete Laplacian updates are fairly natural to consider if one wants to impose a conservation of higher (in particular, dipole) moments while keeping interactions as local as possible.

Proof of localization
We now turn our attention to the case when the local constraints M v take on their minimal value, i.e., we set M v = z v throughout this section.For the 2D model in Fig. 1a, we observed that in this case autocorrelations exhibit a finite saturation value, a signature of localization.We now prove that this is a generic feature of the discrete Laplacian models with M v = z v .
As discussed above, there are two elementary processes that can occur: a site can 'fire', distributing a particle to each of its neighbors (i.e., G(v) = −z v ) or it can 'anti-fire', receiving a particle from all neighbors simultaneously (G(v) = +z v ).As we emphasized, these updates are allowed only if the do not lead to a violation of the constraints m v ∈ [0, z v ].To specify the dynamics, we consider stochastic evolutions where these two types of local updates are applied according to some random rules, at integer times, defining a Markov process with transition matrix T .We assume that this Markov process is reversible, T m,m ′ = T m ′ ,m .The space of particle configurations splits into various connected components; due to the reversibility of the dynamics, for an initial condition belonging to connected component C, there is a unique stationary distribution, which is uniform over all m ∈ C [43].
Our strategy to prove localization is as follows.We fix a vertex v and we identify connected components such that m v takes the same value for all m ∈ C; in particular the values m v = 0 or m v = z v (its minimum and maximum).Any initial configuration from such a C will lead to a finite contribution to the autocorrelation C v (t).What we then need to prove is that a finite fraction of all initial configurations belong to one of these connected components.
In particular, let D(C) = |C| denote the number of configurations belonging to C and D mv (C) denote the number of those where the vertex v is in the state m v (obviously, mv D mv (C) = D(C)).We can define the average value of m v within the sector C as We then have the following expression for the infinite timeaverage of the autocorrelator: where is the total number of configurations and m v = z v /2 is the overall ("infinite temperature") average particle number 3 (for a proof, see App.B).
The key point about this formula is that every term in the sum on the right hand side is non-negative; therefore calculating the contribution from any subset of the terms in the sum in Eq. ( 4) provides a strict lower bound for the time-averaged autocorrelations.
We will construct the appropriate set of connected components by identifying local spin configurations where certain spins remain frozen in the original state at all times, independently of configuration in the rest of the system.The same approach can be used to demonstrate strong fragmentation in kinetically constrained spin systems, such as certain variants of the Fredrickson-Andersen model [35].However, in our case, proving the existence of such frozen blocks is considerably more involved, and it will take up the rest of this section.We begin by proving a useful lemma: Lemma 1. Assume a vertex v fires twice at two different (discrete) times t ′ − 1 and t + 1 (t > t ′ ) such that it does not anti-fire (on net) between them.Then all neighbors of v neighbors must have fired at least once sometime in the time window [t ′ , t].
Proof.We want to keep track of how many times each site fires/anti-fires as the dynamics progresses.Let us define the net number of firings at a vertex F v (t, t ′ ) as the number of times v fires in a time interval [t ′ , t] minus the number of times it anti-fires in the same time interval.
In the situation above we have F v (t − 1, t ′ + 1) = +2 and F v (τ, τ ′ ) = 0 for any t ′ ≤ τ ′ < τ ≤ t.In particular, F v (t, t ′ ) = 0. Now we assume that v has some neighbor v ′ that does not fire in this time interval (more precisely, we only need the weaker condition that it does not fire on net, i.e., F v ′ (t, t ′ ) ≤ 0.) and derive a contradiction.
Note that between the two firings, v has to 're-charge', ∆m v (t, t ′ ) ≡ m v (t)−m v (t ′ +1) = z v .On the other hand, we can rewrite this as the total charge flowing into site v: That is, the charge needed for the second firing has to come from the firing of neighboring sites.If we assume that one of the neighbors does not fire, then another one has to fire at least twice: Now we can run this argument recursively.Let t 1 be the time when v 1 fires for the second time and t ′ 1 the last time it fired before that.Between these two times, v ′ needs to re-charge, which it can only do if its neighbors fire.However, since t ′ < t ′ 1 < t 1 < t we have from our previous assumption that F v (t 1 , t ′ 1 ) = 0. Therefore v 1 needs to have some different neighbor, v 2 that fires at least twice at some times t ′ 2 and t 2 between t ′ 1 and t 1 and so forth.We thus end up with an infinite regress: before our original site could fire for the second time, one of its neighbors has to fire twice, but then of its neighbors has to fire twice etc, making this process impossible (for a sketch of this on the square lattice, see Fig. 3a).The only resolution is to allow all neighbors of v to fire at least once between t ′ and t.We can then use this property to construct frozen blocks of sites that lead to fragmentation and localized dynamics.In particular, we have Corollary 1.Take some set of sites R with the property that for any v ∈ R there is some neighbor v ′ ∈ N v that is also in R. Let R denote the set of vertices in R along with all their neighbors, R ≡ v∈R N v .Consider a configuration where all sites in R have m v = 0. Then the sites in R will remain frozen forever, independently of what the initial configuration is outside of R.
Proof.An example of this situation for a square lattice is shown in Fig. 3b for the 2D square lattice, with dark blue sites belonging to R and light blue ones to the "padding region" P ≡ R \ R. We again use a proof by contradiction.Let us assume that at least one of the sites in R can change.Let us denote by v the first one to do so (site labeled '1' in Fig. 3b).The only way for v to change is if one of its neighbors in P fires.Let us denote this site by v ′ ∈ N v ∩ P (site labeled '0' in the figure).
Since v ′ starts in its lowest state, it must have been completely 'charged up' (i.e., m v ′ = z v ′ ) beforehand.Site v ′ could not have anti-fired up to this point, since it has a neighbor (site v) which has no particles.Therefore the charge must have come from the firing of its other neighbors.Since v is frozen, only the remaining z v ′ − 1 neighbors (sites 2, 3, 4 in the figure) have contributed to charge v up.Hence, one of them needed to fire at least twice.However, due to the previous theorem, this is in contradiction with the fact that v ′ did not yet fired, which finishes the proof.
We can now combine this corollary with Eq. ( 4) to arrive at Theorem 1.The time-averaged "infinite temperature" auto-correlations are finite in the thermodynamic limit.In particular Proof.Let us take R to be a set of two sites, including v and one of its neighbors and consider the set of configurations such that all sites v ′ ∈ R have , where we used that R contains at most 2z 0 sites.Let us now also include all other configurations that are connected to one of these by the dynamics, forming a set of connected components C. For all such C, it follows from the corollary that m v (C) = 0, since v is frozen.Therefore the total contribution from these to the RHS of Eq. ( 4) is (z v /2) 2 ′ C D(C) /D ≥ (z v /2) 2 /(z 0 + 1) 2z 0 .We hence conclude that the time-averaged autocorrelations must remain finite.In particular, for the square lattice model in Eq. 1 which has z v = 4 ∀v, we find C v ≥ 4/5 8 .Theorem 1 is the main result of this section.It shows that the dynamics is localized, in the sense that a finite amount of local memory of initial conditions is preserved forever.

Fragmentation of the configuration space
Let us now discuss what our results imply about the fragmented structure of the space of particle number configurations, i.e. the structure of the connected components C. For simplicity, we focus on a case where the system is defined on a d-dimensional lattice, so there is a notion of a linear size L, with the number of sites growing as ∼ L d .We also take periodic boundary conditions, so that there are no additional global symmetries beyond the multipole moments mentioned above.
Ref. [20] distinguished two types of fragmentation, labeled weak and strong.For the former, once the values of global conserved quantities (in our case, charge and its various multipole moments) are fixed, there is a dominant C that contains almost all configurations in the thermodynamic limit.In particular, this means that the size of the largest connected component has a size D/ poly(L).For strong fragmentation, on the other hand, the dimension of the largest C is an exponentially small fraction of D. The one-dimensional version of the discrete Laplacian model was found to exhibit strong fragmentation for S = 1 [20], and the localization of autocorrelations was derived explicitly from a complete classification of the fragmented components in Ref. [36].
Here, we proved localization more directly, without having to construct the full set of connected components.Nevertheless, our results do imply strong fragmentation.What we have proven above is that the system exhibits a finite frozen site density [44]: any given site has a non-vanishing probability of having its local state preserved by the dynamics at all times.It follows that the probability that a randomly chosen configuration contains a finite fraction of such frozen sites, approaches 1 in the limit of large system sizes 4 .The presence of a finite fraction of frozen sites means that the number of configurations connected to this initial state is exponentially small compared to D, hence implying strong fragmentation.
In fact, there are many more connected components beyond those which we used that we used in the proof of theorem 1, some of which we can also construct by using corollary 1.For example, we can take the region R to be a closed surface, meaning that it separates the graph into two regions -an "inside" and an "outside" -such that any path connecting them must go through R. Then the two sides of R will remain disconnected by the dynamics (an example of this behavior is shown in Fig. 4a).The region R in this case plays the same role as the "inert regions" discussed for the 1D model in Refs.[20,21].One could also put any frozen configuration in the inside of R to get some large frozen island, which would give additional contributions to the autocorrelations as well 5 .
Another quantity related to fragmentation is the total number of frozen states, i.e., spin configurations where every site is frozen.Numerically, we can estimate their number by randomly sampling configurations and checking if they are frozen: the results for the 2D model ( 1) with S = 2 are shown by the dots in Fig. 4b.We can also give an analytical estimate using Pauling's method, which yields that the number of frozen states, N F , scales as ln(N F ) ≈ ln(D) + (L − 2) 2 ln 1 − 2 9  5 5 .While this estimate is non-rigorous, it fits the numerical results very well (see black dashed line in Fig. 4b).Nevertheless, we can also derive a less tight but rigorous lower bound which reads ln(N F ) > 5  9 ln(2S + 1)L 2 for spin S.This is sufficient to show that N F grows exponentially with L 2 (see red line in Fig. 4b for S = 2) for any S, even away from the limit of S = 2 where our proof of strong fragmentation applies.For details of these calculations, see App.D).

Spatially modulated charges and boundary correlations
As observed in Fig. 2a, autocorrelations in the 2D model defined in Eq. ( 1) decay to zero in the bulk for any spin S > 2. However, as panel (b) of the same figure shows, this is not true for correlations near the boundary, when the system is defined with open boundary conditions.Here, we explain this fact in terms of additional conserved quantities that the model possesses in this case.

Recursion relation: Discrete Laplace equation
In general, we can look for spatially modulated conserved quantities of the form Q {αr} = r α r s r .For the model in Eq. ( 1), this ansatz leads to the following recurrence relation, whose solutions correspond to conserved quantities: With open boundaries, one can always solve the recurrence equation by specifying some set of boundary values α B i,j and then propagating them to the bulk.Such solutions give rise to a large number of additional conservation laws.As we show below, these are localized near the boundary and lead to the aforementioned long-lived correlations there, while their effect on bulk dynamics is negligible in the thermodynamic limit.
To prove that this is the case, we need to solve Eq. ( 5) for specified boundary values of α's for (i, j) ∈ B, i.e., α 0,j , α L+1,j , α i,0 , α i,L+1 , where we distinguish between interior sites i, j ∈ {1, . . ., L}, that belong to the bulk of the system D, and those at the boundary B.
To make the solution of the recurrence relation more apparent, we rewrite the equation in a slightly different form: Hence we see that the value at site (i, j) in the bulk, is given by the average value of the four neighboring sites.This equation is a discrete version of the Laplace equation ((∂ 2 x + ∂ 2 y )α(x, y) = 0) with Dirichlet boundary conditions for a square tiling.Its solutions are known as discrete harmonics (see e.g., Ref. [45]).There are two important properties of (discrete) harmonic functions that we will use: 1.A (discrete) harmonic function defined on L takes its maximum M and minimum m values at the boundary B. E.g., if α B i,j ∈ {0, 1} then 0 ≤ α i,j ≤ 1 for all points in the bulk.
2. Using this fact, it follows that the solution is unique: given two harmonic functions f, g on L such that f = g on B implies f i,j = g i,j for all (i, j) ∈ L.
One general way of solving a higher-dimensional linear recurrence equation, and in particular the Dirichlet's problem in Eq. ( 6), is using separation of variables α i,j = X i Y j .With this, Eq. ( 5) simplifies to solving the one-dimensional recurrences where we also indicated the corresponding continuum analogues.λ here is a constant whose possible values are restricted by the choice of boundary conditions.In particular, due to the linearity of the problem, we can decompose the Dirichlet's problem into four independent ones, with α vanishing along 3 out of the 4 boundaries in each case.This leads to a quantization condition on λ and the general solution will be a linear combination of such fundamental solutions.For now, we keep λ as an arbitrary parameter and consider what the form of the solutions admitted by Eqs. ( 7) is.See additional details in Appendix F and in Ref. [46].We can solve Eqs. ( 7) by finding the roots of the associated characteristic polynomials r 2 − (2 ± λ)r + 1 = 0 where the two signs correspond to the equations for X and Y respectively.For each equation, the two roots are inverses of each other and are either real (ξ x,y ) or pure complex phases (e ikx,y ), depending on the value of λ.Overall, we can identify three main types of solutions: (i) When λ = 0, we have ξ x = ξ y = 0; this case contains the multipole conserved quantities discussed earlier6 .(ii) When 0 < |λ| < 4, one of the solutions is real while the other one is complex, e.g., α i,j ∝ (ξ x ) i e ikyj .These correspond to conserved quantities that are exponentially localized in one dimension while being fully delocalized in the other.(iii) When |λ| > 4, both solutions are real α i,j = (ξ x ) i (ξ y ) j , and hence can be exponentially localized near one of the corners of the system, depending on the modulus of ξ x and ξ y .
Apart from giving insight into the family of possible conserved quantities, these fundamental solutions can be combined to obtain the unique solution corresponding to a choice of boundary values α B i,j .Instead of writing that general expression, we notice that close-form solution of equation ( 6) is already known in the context of two-dimensional (unbiased) random walks [47].This is given by where each sum corresponds to one of the four boundaries with the respective boundary values α B i,j and Hence, α i,j is given by the (discrete) convolution of T a,b (i, j) with α B i,j .T a,b (i, j) is the discrete analogue of the Poisson kernel, whose convolution with the boundary condition solves the continuum Dirichlet problem [48].We will use this fact to study the behavior of the solutions of the discrete recurrence relation Eq. ( 5) in the following section.
While Eqs. ( 8) and ( 9) provide a way to construct the exact solutions for any boundary condition, in practice we instead solve the recurrence relation numerically, applying an iterative method outlined in App.E. This has the advantage that it can be easily generalized to other cases where the exact kernel is not known, some of which are studied in Sec.4.3.

Boundary localized charges and Mazur's bound
A powerful tool to analytically prove finite boundary correlations is given by Mazur's bound [49] M i,j , which lower bounds the infinite time-average auto-correlations by the overlap with a set of conserved quantities Q a .More explicitly this is lim at any site (i, j).To define the bound, we introduce an inner product between observables as ⟨A, B⟩ ≡ ⟨AB⟩ where ⟨A⟩ = 1 |C| {s i,j }∈C A(s i,j ) is the "infinite-temperature" average.With this definition, the bound can be written as Here, K is a positive-definite matrix with elements K a,b = ⟨Q a , Q b ⟩.If one includes only a single conserved quantity Q {αr} = r α r s r , then the expression simplifies to where ⟨Q {αr} , Q {αr} ⟩ is simply given by i.e., proportional to the 2-norm of α.Hence, if one finds even a single conserved quantity that is strongly localized around some particular site at the boundary of the system and whose norm remains finite in the thermodynamic limit, this will give rise to finite boundary correlations.Alternatively, even if individual charges are not sufficiently well-localized, they can still be combined to give rise to a finite value on the RHS of Eq. (11).

Corner charges: Exponential localization
We first consider the situation near the corners of the square lattice.On a square lattice, the corner sites are completely decoupled under the dynamics, and hence provide a trivially conserved boundary charge.Let us ignore these, and consider instead the region close to the corner with coordinates r 0 = (0, 0).As we noted in Sec.4.1, the recursion relation has fundamental solutions of the form α ij ∝ (ξ x ) i (ξ y ) j .This suggests conserved quantities that are exponentially localized near one of the four corners of the system, similar to the exponentially localized quantities found for certain 1D systems in Ref. [37].Indeed, we can plug the above ansatz directly into the original recursion relation (5), which then turns into Solutions to this, restricted to the domain |ξ x |, |ξ y | < 1 exist (see Fig. 5a) and provide exponentially localized modes at the corner (0, 0) 7 .
Figure 6: Boundary charges with n > 1. Solutions of Eq. ( 6) with n = 2 (panel (a)) and n = 4 (panel (b)), which decay towards the bulk as α (2) 0,j ∼ j −3 and α (4) 0,j ∼ j −5 respectively.Red and blue correspond to positive and negative values of α We can plug these conserved quantities into Mazur's bound (12) at a site (x, y) close to the corner (i.e,.|x|, |y| not scaling with L) and including only this conserved quantity to find in the limit L → ∞.Therefore, as long as the site (x, y) is at finite distance from the corner, Mazur's bound will be finite, which in turn shows that the time-average boundary correlation saturates to a finite value.The optimal lower bound can be found by maximizing the expression on the right hand side of Eq. ( 15) over the set of solutions of the equation (14).The optimal bounds thus obtained for different coordinates (x, y) are shown in Fig. 5b.

Mid-boundary charges: Power-law localization
In Sec.4.1 we found that there exist fundamental solutions of the recursion equation of the form α i,j ∝ (ξ x ) i e ikyj , i.e. exponentially localized in one direction and plane-wave-like in the other.These are symmetries localized at the boundary.However, they are not normalizable and do not immediately yield a finite Mazur bound.Here we instead ask whether we can find solutions that are localized in all directions.To answer this, we return to Eq. ( 6), and recall that it is a lattice discretization of the continuum Laplace equation.
Motivated by this, we first consider the problem in the continuum, which will guide us in constructing new boundary charges and explaining their localization properties.
As we are interested in the long-distance behavior of boundary charges, we consider Laplace's equation on a semi-infinite plane, i.e., on the domain (x, y) Given a boundary condition f (x) ∈ L 1 (R), the solution is given by the convolution of f with the Poisson kernel P (x, y) = 1 π y x 2 +y 2 [48,50], i.e., α(x, y) = ∞ −∞ dz P (x − z, y)f (z), analogous to the discrete case.In particular, we want to study the decay of a solution localized around (0, 0).First, let us take f (x) = δ(x), the Dirac delta distribution.This leads to a diverging solution at the boundary.Nevertheless, our goal is to understand the long distance behavior where the solution is well-behaved.Alternatively, we could use a regularized version of the delta distribution instead.In any case, this boundary condition gives α(x, y) = P (x, y), which decays as α ∼ 1/r at large distances.Recall that r refers to the orthogonal distance to the boundary.Hence, while it decays towards the bulk, the decay is too slow for the norm dx dy |α(x, y)| 2 to converge.
Can we construct more localized solutions?Note that the above solution is reminiscent to the electrostatic potential generated by a point-like source at the boundary.This suggests a natural way of constructing solutions with a faster decay: replace the point-like source with a dipole or some higher multipole source.This can be achieved by using a boundary condition that is a higher derivative of the delta distribution: Using this boundary condition one finds, after integrating by parts, the solution α (n) (x, y) = (−1) n ∂ n z P (x − z, y)| z=0 .These decay asymptotically as α (n) (0, y) ∼ y −(n+1) , making them increasingly localized as we make n larger.In particular, for any n ≥ 1 they decay sufficiently quickly to make their norm convergent.
We now want to find an analogous set of solutions on the lattice.To do so, we first need to discretize the boundary condition f (x) = δ (n) (x).This can be simply done by replacing the derivatives by (central) finite differences of the Kronecker delta δ i,0 .In particular, the n-th derivative of a lattice function f i at site i is given by ∆ k , (here and below we focus on n even).Hence, for f i = δ i,0 , this requires fixing n + 1 non-zero boundary values given by for i ∈ {− n 2 , . . ., 0, . . ., n 2 }.We have normalized α (n),B such that α (n),B 0,0 = 1 for all n.Given this boundary condition, we can find the solution of Eq. ( 5) in the bulk.In Fig. (6) we plot the solution for n = 2 and n = 4.The two different colors emphasize the change of sign of the solution along the nodal lines where α (n) vanishes.E.g., α (2) is positive in the central lobe (red) and negative on the sides (blue).A higher n, gives n nodal lines with the corresponding changes of sign.The solutions α (n) give rise to a set of conserved quantities which we denote by Q (n) r 0 , where r 0 refers to the location around which they are localized (r 0 = (0, 0) in the discussion above).
The solutions α (n) inherit their properties from their continuum counterparts.For n = 0, where the boundary condition is simply a Kronecker delta, we find a 1/r decay in α (0) i,j , as shown in Fig. 7a.This means that the norm of Q (0) diverges logarithmically with the linear size of the system, as shown in Fig. 7b.Hence, these charges are not sufficiently strongly localized at the boundary for a single one of them to give a finite Mazur's bound.Instead, we can make use of the general expression for Mazur's bound (11) and include all of them simultaneously.We find that this is sufficient to obtain a lower bound that is tight to the result obtained in our numerical simulations for any S (see black dashed line in Fig. 1b).Additional details can be found in Appendix A and Ref. [46].
On the other hand, in agreement with the continuum case, the solutions for higher n decay faster, as r −(n+1) , as we confirm numerically in Fig. 8a.As a consequence, a single one of these charges is sufficient to give a finite Mazur's bound, using Eq. ( 12).However, while the lattice and continuum problems match at long distances, there are differences in their behavior close to the boundary that affect Mazur's bound.In particular, 2 ) i,j in the limit L → ∞.
while the charges become more strongly localized towards the bulk for higher n, they also become more spread out at the boundary, which leads to an increase in their 2-norm.As a consequence, the value of the bound (12) in fact decreases with n (while remaining finite for any finite n).Let us provide a rough estimate.For sufficiently large n, α (n) decays quickly away from the boundary.Motivated by this, we estimate the scaling of ∥α (n) ∥ 2 2 with n, by replacing it with the norm of the contribution from the boundary alone, . The latter can be evaluated from Eq. ( 17) and gives ∥α (n),B ∥ 2 2 ∼ √ n in the limit of large n.We evaluate Mazur's bound exactly for various n in Fig. 8b, and find that it agrees with this estimate.The same figure also shows Eq. ( 11) evaluated using all n = 0 charges taken together, which gives a stronger bound.In Appendix A we show the finite-size scaling of Mazur's bound with system size for several n.
All together, we managed to construct different families of conserved quantities {Q (n) r 0 }, which become more and more localized towards the bulk when considering higher and higher finite differences.Different Q (n) r 0 are localized around a boundary site with coordinates r 0 .For a given n, α (n),B vanishes everywhere except on n + 1 sites centered around r 0 .The set of charges where r 0 are at distance n 2 along the boundary are then trivially linearly independent.Hence, this proves that there exist at least O(L) of those.Nevertheless, different families parametrized by different n's are not linearly independent of each other.In particular, we can use the n = 0 charges to construct all other families.

Generalizations
So far we have focused on the particular square lattice model defined in Eq. ( 1).However, our construction can be extended to higher dimensions as well as to different lattices as long as the associated recurrence relation   (12) with n, when including a single higher-moment charge with modulation α (n) .Mazur's bound including all point-like charges at the boundaries is shown as a black triangle.For this data we used a linear system size of L = 300.example, while the models in Figs.1(a,b) are defined on different lattices, they both correspond to different discretizations of the same continuum Laplace equation.Hence the same long-distance decay for conserved quantities localized at the boundary of the system applies to them.
To extend the construction to higher dimensions, we now consider the continuum Laplace equation in d + 1 dimensions in the hyperplane H d+1 = {(x, y) ∈ R d+1 |y > 0}.The corresponding Poisson kernel is a generalization of the 2D one and reads P d (x, y) = c d y(y 2 + ∥x∥ 2 ) − d+1 2 , with some dimension-dependent constant c d [50].Once again the squared 2-norm of the (regularized) solutions with α B (x, 0) = δ(x) diverges logarithmically in the thermodynamic limit and hence these are not sufficient to prove finite boundary correlations.Instead we can again consider "multipole" boundary conditions.Fixing , the decay of the solution at large distances is given by α (n d ) (0, y) ∝ y − d i=1 n i −d .Therefore, we find that one can construct charges that are sufficiently localized at the d-dimensional boundaries of the system as to provide a finite Mazur's bound.
The previous discussion appears to suggest that as long as the continuum limit of the linear recurrence Eq. ( 18) is given by the Laplace equation, one should be able to show that the boundary correlations are finite.However, some subtleties need to be addressed.First of all, recall that we explicitly made use of the fact that solutions of Eq. ( 18) are discrete harmonic functions.This ensured that for any boundary condition (and any system size) there exist a unique solution.If this was not the case, even the existence of a solution is not ensured.One such example is the 2D model defined in Eq. ( 5) of Ref. [37].Our construction of boundary charged does not apply in that case; however, we have observed numerically that boundary correlations nevertheless fail to decay even in that model.
One set of models to which our previous discussion does directly apply are generated by the following local gates: with N = i n i .These correspond to the following recurrence relation, with p i = n i /N .These equations were solved exactly in Ref. [47].When n 1 = n 2 and n 3 = n 4 , the continuum limit of this recurrence reads where we defined p = 2p 1 .Solutions of this equation can be found by performing the change of variables x → x/ √ p, y → y/ √ 1 − p, and then are given by the evaluation of the solution of the isotropic problem ( 16) at α in (x, y) = α(x/ √ p, y/ √ 1 − p).Hence, our discussion about the localization properties of boundary charges also applies to these family of anisotropic models.

Quantum systems and strong zero modes
So far our discussion has focused on classical Markov chain dynamics.However, it is easy to map the systems we studied onto quantum Hamiltonians that share the same symmetries.Given a set of local gates G r = {n v } v∈Rr , acting on some region R r centered on a site r, one can construct a quantum spin-S Hamiltonian ĤG = r J r ĥr , on the same lattice as ĥr = v∈Rv Ŝ sgn(nv) where sgn(n v ) is the sign of n v , and J r is an arbitrary choice of real coefficients.By construction, when written in the computational basis, ĤG has the same block-diagonal structure as the original Markov generator built from the gates G r .In particular, it shares both its fragmentation properties and its spatially modulated symmetries (with r α r s r replaced by r α r Ŝz r ).These properties are also preserved by any additional diagonal terms V ({ Ŝz r }) that are function of the Ŝz r only.However, ĤG might also have additional symmetries, not present in the classical model.Of particular interest are Z 2 symmetries like R x = r e iπ Ŝx r and R y = j e iπ Ŝy j .Indeed, R x Ŝ± r R x = Ŝ∓ r and R x Ŝz r R x = − Ŝz r so the Hamiltonians constructed above are invariant under R x as long as V is built up from even products in the local Ŝz operators.The importance of these additional discrete symmetries stems from the fact that they anticommute with the spatially modulated symmetries and therefore lead to exact degeneracies in the many-body spectrum [38].
To have these degeneracies, any spatially modulated symmetry of the form r α r S z r is sufficient, along with one of the aforementioned Z 2 symmetries.For the discrete Laplacian models there are always such conserved quantities, independent of the choice of boundary conditions.However, we can modify the models in a way such that only the charges localized at the boundaries remain, which exist only for open boundaries.We do so by defining the local gates with associated recurrence relation reads Choosing N > 4 rules out the conservation of the global charge and hence of any of its higher-moments.In general, solutions of this recurrence equation are not discrete harmonics, but rather correspond to the eigenvalue problem △α = εα with ε ̸ = 0 in the continuum.In this case, one only finds solutions of the form α i,j = (ξ x ) i (ξ y ) j and α i,j = (ξ x ) i e ikyj , e ikxi (ξ y ) j ; both localized near the boundary.Importantly, this family of models does not show spatially modulated global conserved quantities for periodic boundary conditions.We thus expect them to feature similar phenomenology to that of the strong zero modes (SZM) introduced in Ref. [38], with degeneracies throughout the many-body spectrum for open, but not for closed boundaries.This construction can be extended to higher dimensions by for example considering a d-dimensional cubic lattice with n 0,0 = −N , where N > 2d, and the remaining contributions equals to +1 as in Eq. ( 23).These are higher-dimensional generalizations of the 1D models with exponentially localized symmetries introduced in Ref. [37].
Nevertheless, despite the similarities with SZM, some important differences remain.First, while the boundary modes of Ref. [38] are only approximately conserved for finite systems, ours are exact for any L.This somewhat changes the logic of the construction: Instead of using the zero mode to toggle between the two different symmetry sectors of the exact Z 2 symmetry, we can also classify energy eigenstates using the boundary symmetries.Arguably, the most important difference to standard SZM is that our boundary modes correspond to continuous, rather than discrete symmetries (i.e., there is no "normalization condition" of the form ( Q) m = 1 for any integer m).This condition appears to be "highly non-trivial and fundamental" [51] to ensure a non-zero radius of convergence in the perturbative construction of such modes (indeed, our boundary modes are presumably highly sensitive to any additional perturbations).While this condition appears to hold for previous constructions of SZM found in the literature, it is not clear whether it should be generally imposed [51].We therefore leave it to future work to determine whether the models introduced here can be meaningfully fit into the framework of SZM.

Conclusions and outlook
In this work we studied a family of models, which we named discrete Laplacian models, and which can be defined on an arbitrary lattice (or, more generally, bounded-degree graph).We proved two main results about these models.First, we proved that bulk autocorrelations saturate to a finite value when the on-site configuration spaces are chosen to take their minimal values.Our proof works by explicitly constructing spatial regions whose configurations are left unchanged by the dynamics, which can then be used to divide the rest of the system into disconnected regions.
Secondly, we constructed a hierarchy of linearly (in the linear system size) many conserved quantities, which are localized at the boundary of the system.These are new instances of spatially modulated symmetries whose modulations satisfy a discrete Laplace equation with Dirichlet boundary conditions.We showed that while these are only powerlaw (rather than exponentially) localized, the power of decay can be systemically increased by choosing appropriate boundary conditions.As a result, we are able to prove that boundary correlation are finite, by making use of Mazur's bound [49].
There are a number of questions left open by our work.One interesting aspect of our work is the construction of strongly fragmented models in higher dimensions, whereas previous examples in the recent literature were explicitly one-dimensional in nature, relying on the conservation of certain one-dimensional patterns due to hard-core constraints [25,36].While we proved the presence of strong fragmentation in discrete Laplacian models, constructing the conserved quantities to label all connected components, and understanding their algebraic structure in the spirit of Refs.[25,36], is an interesting challenge.This would shed more light on their behavior and would help in finding other strongly fragmented higher dimensional models.
While we proved strong fragmentation only in the case when we restricted the local spin (number of particles per site) to its smallest possible value, it is likely that even away from this limit, one would find a transition from a weakly to a strongly fragmented regime by tuning the density of particles per site, analogously to the one-dimensional case studied in Ref. [44] ( as well as other kinetically constrained models [35,52]).Understanding the higher dimensional versions of this transition is another interesting open problem.
The relative novelty of the new classes of modulated symmetries we uncovered also naturally leads to many open questions.The most pressing one is to understand their level of fine-tuning.In particular, what minimal mathematical structure is required to realize such symmetries, and to which extent it is sufficient for them to just be approximately conserved.Moreover, this work provides a stepping stone for a more ambitious goal: classifying all the possible types of spatially modulated symmetries that a physical system with local interactions can possess.While in 1D (with a finite number of modulated symmetries) it appears that only α j = r j with r an algebraic number are allowed, the infinitely many number of conserved quantities appearing in higher dimensions permit richer modulations.
On a different note, the constructions of quantum models in Sec. 5 with conserved quantities that are either localized in one direction and plane-wave like in the other, or localized at the corners is reminiscent of topological phases of matter where low-energy modes localized at the boundaries appear.The conserved quantities associated to corners in particular resemble the corner modes of higher order topological modes, which have been previously linked to systems with multipole symmetries [53].Whether these analogies can be pushed further is an interesting open problem.

A Finite-size Scaling Analysis
In this appendix we provide a more in-detail analysis of the finite-size scaling of the various numerical results we discussed in the main text.Fig. 9a shows that the finite saturation value of bulk auto-correlations for the model in Eq. ( 1) when S = 2, does not scale down with system size, being the data converged in system size for all shown times.This data has been obtained for periodic boundary conditions and averaged over N = 100 circuit realizations.Fig. 9b provides the finite-size scaling of boundary auto-correlations.Here, we required N = 400 simulations to decrease the late-time fluctuations.
In the following we show the scaling of Mazur's bound with system size.We start with its scaling when only including a single higher-moment charge Q (n) r 0 .This is shown in Fig. 10a for several n = 2, 4, 6, 8, 10.In addition, we also study the finite-size scaling when including O(L) point-like charges Fig. 10b.We first note, that while linearly independent, these conserved quantities are not orthogonal with respect to the infinite temperature inner product and thus, to compute M sx,y for spin correlations ⟨s x,y (t)s x,y (0)⟩ at the boundary, we need to use the general expression (11).Without loss of generality and due to the Z 4 symmetry of the lattice and local gates, we can restrict ourselves to focus solely on one boundary, e. g. on a site with coordinates (0, y) where 1 ≤ y ≤ L. As it turns out, we do not need to include all conserved quantities but only the ones of the form α B,a i,j = δ i,0 δ j,a for 1 ≤ a ≤ L, as all the remaining boundaries have exponentially small overlap with s x,y .Consequently, one finds  In the limit L → ∞, this can give a finite value depending on the scaling of the diagonal matrix elements (K −1 ) y,y with L. Deriving the scaling of a particular matrix element (K −1 ) y,y with system size, however, is an analytically difficult task.We provide a finitesize scaling in Fig. 10b  This allows us to limit our calculations to a more general property of the matrix K.The scaling with system size is shown in Fig. 10b.
B Derivation of Eq. (4) Consider Eq. ( 2) and split the sum over initial configurations into sums over connected components: m(0) = C m(0)∈C .Note that since we are considering connected correlations, it does not matter whether we are using the original variables s v or their shifted versions m v ; here we use the latter to make contact with the proofs of Sec.3.2.We now consider the contributions from each component C.
Let us write the average over random trajectories as ⟨m v (t)⟩ m(0) = m p m(0) (m, t)m v , where p m(0) (m, t) = T t m,m 0 is the probability distribution over possible spin configurations, conditioned on a given initial configuration m(0).Let us now focus on ⟨m v (t)⟩ m(0) for a particular initial configuration belonging to a specific connected component C. Restricted to this, the dynamics is by definition irreducible.Since we also assumed reversibility, this implies a unique stationary distribution, which is uniform over all m ∈ C [43].We therefore have lim Using this, the time-averaged expectation value is ⟨m v (t)⟩ m(0) = m v (C), as defined in the main text, which is independent of m(0) within the same C. Plugging this back into Eq.( 2) we get as promised.

C Strong fragmentation of the configuration space
We will argue that the constructions of Sec.3.2 imply a strong fragmentation of the configuration space in the sense that all connected components contain at most an exponentially small fraction of all configurations.The logic of our argument is as follows.Consider a system of N sites with M possible states on each.If a configuration has n frozen sites, i.e. sites whose state cannot evolve, than it is connected at most to M N −n other configurations; this is exponentially small if n/N is finite in the limit of large N .Therefore, any putative large component that is not exponentially small must contain only a vanishingly small fraction of frozen sites.If we can show that the total number of configurations where n does not grow proportionally to N is itself exponentially small compared to D then we have ruled out the possibility of a large (not exponentially small) connected component.This is what we will do below.We consider a d-dimensional lattice of linear size L and partition it into regions ("boxes") of linear size ℓ; the number of such boxes is N ℓ = (L/ℓ) d .Each site has M = z + 1 possible states, where z is the coordination number of the lattice.This means that each box has D ℓ = M ℓ d configurations, while the whole system has D = M L d = D N ℓ ℓ .Corollary 1 shows that we can construct finite-sized patterns of maximally polarized spins where some sites (those on the "inside", i.e., the region R in Fig. 3) are frozen independently of what the configuration is outside of the pattern.If we choose ℓ to be large enough (i.e., ℓ ≥ 4 for a hypercubic lattice), then each box can contain such a frozen pattern.This will happen with a finite probability p ℓ that depends on ℓ but not on the overall system size L. In particular, let F ℓ be the number of configurations of a box that contain a frozen pattern, then p ℓ = F ℓ /D ℓ .Let us denote by G ℓ = D ℓ − F ℓ the number of remaining configurations and q ℓ = G ℓ /D ℓ < 1 their probability.
We now ask the question: what is the probability that a randomly chosen initial spin configuration containts at most n frozen sites.Let us denote this quantity P ≤n .We want to show that P ≤n will be exponentially small in the large L limit, unless n is a finite fraction of L d (in particular, n ≳ p ℓ L d ).To do so, we note that we can upper bound it in the following way.Let us choose m boxes which contain a frozen pattern, while the remaining N ℓ − m boxes have no frozen patterns.Clearly, any such state has at least m frozen sites, so the only contributions to P ≤n come from m ≤ n.This gives the following Figure 11: Frozen lattice tiling.Tiling defined by the unit cell (grey box) such that the spin configuration on the entire lattice is frozen for all times.Blue sites are highest, red sites lowest charge.All sites are frozen for any choice of s r ∈ {−S, . . ., S} on the white sites.
upper bound: The expression on the right hand side is the cumulative distribution function of the binomial distribution.At large N ℓ , the binomial distribution is increasingly sharply peaked around its mean value, n = p ℓ N ℓ , and the cumulative probability at values of n much smaller than n (i.e., outside a window of size ∝ N 1/2 ℓ ) will be exponentially suppressed in N ℓ .In particular, the probability of configurations where the number of frozen sites is not a finite fraction of L d goes to zero exponentially in the thermodynamic limit.

D Counting of frozen states
To estimate the total number of frozen states of G + using Pauling's method, we consider the set of 5 sites on which a single gate acts and calculate the probability P F that the configuration cannot be changed by the effect of the gate.To do so, we first compute the probability that the gate can fire (or anti-fire) at the beginning and then take the complement.Given our five degrees of freedom s ∈ {−2, −1, 0, 1, 2} placed on a +-shape, we immediately see that in order for the gate to fire the central site has to be either ±2 and the outer sites are restricted to ∓{−1, 0, 1, 2}.With the total 5 5 possible configurations, we have a "firing" probability 2 • 1 1 • 4 4 /5 5 = 2 9 /5 5 .Taking the complement of the probability gives P F = 1 − 2 9  5 5 ≈ 84%, so more than 80% of all initial +-shaped configurations of five charges are frozen.Say we now have a square lattice with L 2 sites.There are in total (L − 2) 2 +-shaped configurations in the lattice all of which have the same probability P F to be frozen initially.Neglecting the effects of overlapping configurations, the probability that the whole lattice is initially frozen and thus frozen for all times is simply given by P . Therefore, the Pauling estimate for the total number of frozen states takes the While the Pauling estimate agrees well with the numerically calculated value of N F (see Fig. 4(b)), it is an uncontrolled approximation.We now derive a rigorous, albeit less tight, bound; it will also have the advantage that it applies to arbitrary spin S. Imagine again a square lattice of size L. We tile the lattice with some unit cell and fix the configuration of spins on some (proper) subset of sites in each unit cell in such a way that these prevent all sites from firing (or anti-firing), no matter what configuration we put on the remaining sites.If the fraction of fixed sites is ϕ, it directly follows that N F > (2S + 1) (1−ϕ)L 2 .There a several possible tilings that yield lattices frozen at all times; Fig. 11 shows one particular example, using a 3 × 3 unit cell with four fixed sites.The finite fraction is therefore ϕ = 4/9 and the lower bound becomes N F > (2S + 1) 5 9 L 2 .This is sufficient to prove that the number of frozen states scales exponentially with L 2 .

E Numerical solution of linear recurrence relations
Finding a close-form solution of linear recurrence relations in terms of (arbitrary) boundary conditions is in general difficult.Hence, given the wide variety of recurrence relations we study in this work, we instead follow a numerical approach based on Iterative Stencil Loops [54,55] we then compare with analytical results in the continuum limit.As in the case of G ⊞ , we can arrange the general recurrence relation v∈Nr n v α r+v = 0 (E.9) to yield an "averaging"-type recurrence relation: A route for solving the latter is recursively constructing a solution by obtaining the value at a site r as given by the RHS of Eq. (E.10).To achieve this, we (1) start with an empty lattice initialized with the imposed boundary conditions α B r for r ∈ B, (2) take the right sum as our "stencil" [54] and (3) iterating over the lattice updating every α r in the bulk with the respective weighted sum of its neighbors, in the sense that α are the values at the site before and after update, respectively.In particular for a lattice of size |L|, this results to each iteration step taking O(|L| 2 ) applications of the stencil.To terminate the algorithm, we compute the difference between two consecutive iteration steps L t and L t+1 and stop the algorithm whenever a certain threshold is reached r∈D |α r | < ε.In particular, we choose ε < 10 −10 .In the case of discrete harmonic functions, the algorithm is known to converge to the unique solution of the discrete Dirichlet's problem.

F Solution of Eq. (6) via separation of variables
Analogously to the continuum case, we can solve Eq. ( 6) via separation of variables, i.e., α i,j = X i Y j .Let us explicitly follow the main procedure which will teach us something about the structure of the solutions.Using this ansatz the recurrence relation becomes 4X i Y j = (X i+1 + X i−1 ) Y j + X i (Y j+1 + Y j−1 ) . (F.12) After dividing by X i Y j (assuming α i,j does not vanish in the bulk) we find 4 This implies that each term on the right hand side is a constant function of its argument and hence, they satisfy the one-dimensional recurrence relations along the horizontal and vertical lattice directions respectively, for certain values of λ ∈ R which are fixed by the boundary conditions.To proceed we notice that being Eq. ( 5) linear, the Dirichlet problem can be solved adding up the solutions of four different Dirichlet problems with vanishing boundaries conditions except at a given boundary.This means that a general solution can be written as i,j + α i,j + α i,j + α For example, solving problem (I) leads to Y n,j = sin(k n y j) and X n,i = sinh(κ n x i) with k n y = nπ/(N + 1) for n = 0, . . ., N + 1.This also restricts the values of 0 < λ < 4 to those satisfying cos(k n x ) = 1 − λ n /2, and in turn cosh(κ n x ) = 1 + λ n /2.Hence, we have found the fundamental solutions A n sinh(κ n x i) sin(k n y j), which by linear superposition lead to the general solution of problem (I) A n sinh(κ n x (N + 1)) sin(k n y j), (F.18) i.e., proportional to the Fourier coefficients of χ 2 (j).Solving the three remaining problems, which again involve products of sinusoidal and hyperbolic functions, one can find a general solution of Eq.( 6) for any choice of boundary functions.This implies that the fundamental solutions sinh(κ n x i) sin(k n y j), sinh(κ n x [i − (N + 1)]) sin(k n y j), sin(k n x i) sin(κ n y j) , sin(k n x i) sinh(κ n y [j − (N + 1)]) with n = 0 • • • , N + 1 form a basis for solutions of Eq. ( 6), showing that the number of independent symmetries scales with the linear system size.This approach not only gives us the fundamental solutions from where to obtain any other one, but also allows us to explicitly find particular spatial modulations that connect to the quasi-periodic and exponential modulations we for the first and second equations respectively.As the characteristic polynomial are palindromic, the two roots are inverse of each other.Whether these are real (η x,y ) or pure complex phases (e ikx,y ) is determined by the sign of the discriminants ∆ x = λ(λ+4) , ∆ y = λ(λ − 4): If ∆ x,y > 0, the corresponding solution can be written in terms of hyperbolic or exponential functions; while ∆ x,y < 0 correspond to sinusoidal or complex exponential modulations.
We can split the solutions into three main types: (i) λ = 0 contains all multipole conserved quantities we already identified except for the x 2 −y 2 quadratic moment.As this solution does not factorize in horizontal and vertical directions, but rather as α x,y = (x + y)(x−y), we can only recover it by superposing many fundamental solutions.Alternatively, these could be explicitly found when writing the recurrence relation in terms of center of mass s = i + j and relative r = i − j coordinates, i.e., α i,j = X s Y r .The second case (ii) corresponds to exponential (hyperbolic) solutions, which are localized near the corners of the 2D lattice when |λ| > 4.These can be directly found solving Eq.( 5) with the ansatz α i,j = (ξ x ) i (ξ y ) j .Finally, the third case (iii) corresponds to the product of sinusoidal and hyperbolic solutions along orthogonal directions, like e.g., α i,j = (ξ x ) i e −ikyj .Thus, while solutions with a finite momentum mode along one of the lattice directions exist (for |λ| < 4), these are exponentially damped along the orthogonal direction and do not contribute to spin correlations in the bulk.Moreover, we can also rule out solutions of the form α i,j ∼ e ikxi e ikyj which would have led to conserved finite modes in the bulk correlations.

Figure 2 :
Figure2: Spin-spin correlations of G + (a) Spin-spin autocorrelation C(0, t) for linear system size L = 300 and different spins S. Data is averaged over N = 1000 random initial configurations and circuit realizations.When S = 2, the autocorrelation saturates to a constant value independent of system size, hence compatible with a strong fragmentation of the configuration space.(b) Boundary autocorrelation for S = 3.The black dashed line shows Mazur's bound as computed in Eq.(11).Inset: Spatial distribution along the xaxis for L = 250 and different times t.Unlike the weakly fragmented case, the correlation remains localized even for longer times.

Figure 3 :
Figure3: Proof of localization.(a) Sketch of the argument leading to lemma 1.We want site 0 (dark red) to fire twice, while not allowing its neighbor to the right (gray) to fire in-between.This requires another neighbor (site 1) to fire twice, leading to an infinite regress shown by the red arrows.(b) Construction of a cluster of frozen sites.Blue sites are taken to have zero particles m v = 0.The light blue sites on the outside layer can change but the inner ones remain frozen forever, independently of the initial configuration on all other sites.

Figure 4 :
Figure 4: Fragmentation.(a): Construction of various sectors in configuration space using corollary 1, with inert walls splitting the lattice into disconnected regions.(b) Estimate for the number of completely frozen configurations, along with a Pauling estimate.Numerical simulations use system sizes with L ≤ 12.

Figure 5 :
Figure 5: Exponentially localized charges at the corners.(a) Solutions of Eq. (14).The black dashed line shows the unit circle.(b) Largest value of Mazur's bound in Eq. (15) over the set of solutions plotted in panel (a), at different coordinates (x, y).

Figure 7 :
Figure 7: Localization of n = 0 charges.(a) Decay of α (0, L 2 ) i,j towards the bulk along two different directions specified by the two planes in the 3D plot.In both cases, the solution decays as a 1/r with the distance r.(b) Divergence of the 2-norm of α (0, L 2 ) i,j v∈Nr n v α r+v = 0,(18) corresponds to a discrete Laplace equation, i.e., its solutions are discrete harmonic functions.Some examples of models satisfying this requirement are shown in Fig.1.For

Figure 8 :
Figure 8: Mazur's bound and higher-moment distributions.(a) Power-law decay of higher-moment charges with modulation α (n) 0,j ∼ j −(n+1) as predicted from the continuum Laplace equation.(b) Dependence of Mazur's bound(12) with n, when including a single higher-moment charge with modulation α (n) .Mazur's bound including all point-like charges at the boundaries is shown as a black triangle.For this data we used a linear system size of L = 300.

Figure 9 :
Figure 9: Finite-size scaling of correlations.(a) Finite-size scaling of bulk autocorrelations for S = 2 (a) and boundary auto-correlations (b) for model in Eq. (1).The data has been averaged over N = 100 and N = 400 circuit realizations respectively.

Figure 10 :
Figure 10: Finite-size scaling of Mazur's bound.(a) Finite-size scaling of Mazur's bound when including a single charge with modulation α (n) (panel a), and when including O(L) point-like charges.leading to (orange dots).Alternatively, we can instead look for a lower bound of the averaged boundary auto-correlations lim ⟨s r (t)s r (0)⟩ ≥ M B , (A.4) where Mazur's bound can be expressed in terms of the trace of the inverse matrix K −1

A
n sinh(κ n x i) sin(k n y j), (F.17)where the coefficients A n are fixed byχ 2 (j) = N +1n=0 encountered in the previous chapter.In general, we can solve Eqs.(F.14) finding the roots of the associated characteristic polynomials r 2 − (2 ± λ)r + 1 = 0 which take the values