Ergodicity-breaking arising from Hilbert space fragmentation in dipole-conserving Hamiltonians

Pablo Sala, 2, ∗ Tibor Rakovszky, 2 Ruben Verresen, 3 Michael Knap, 2, 4 and Frank Pollmann 2 Department of Physics, Technical University of Munich, 85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 München, Germany Max-Planck-Institute for the Physics of Complex Systems, 01187 Dresden, Germany Institute for Advanced Study, Technical University of Munich, 85748 Garching, Germany (Dated: April 10, 2019)


I. INTRODUCTION
Recent years have seen a great deal of effort-both theoretical and experimental-to understand quantum thermalization: the question of how closed quantum systems, evolving under unitary dynamics, reach a state of thermal equilibrium [1][2][3][4][5][6][7][8]. Thermalization is believed to be characterized in terms of the Eigenstate Thermalization Hypothesis (ETH) [7,[9][10][11]. According to this, each eigenstate of a thermalizing Hamiltonian essentially behaves like a thermal ensemble as far as expectation values of local observables are concerned. While no proof of ETH exists, there are many cases where it has been shown numerically that indeed all eigenstates satisfy this hypothesis [7,11,12].
Another key question is about the possibility of systems that exhibit interesting intermediate behavior, neither localized, nor fully ergodic. In particular we can distinguish between strong and weak ETH: the former says that all eigenstates in the bulk of the spectrum become * pablo.sala@tum.de FIG. 1. Thermalization and its absence in the autocorrelation function. Panel (a) shows the auto-correlation function C z 0 (t) ≡ S z 0 (t)S z 0 (0) in the full Hilbert space at infinite temperature for N = 13 (transparent curves) and N = 15 (opaque curves) spins. For Hamiltonian H3 in Eq. (1), C z 0 (t) saturates to a finite value at long times, closely matching our prediction in Eq. (5) (dashed line). The auto-correlation function of the combined Hamiltonian H3 + H4 decays to zero at long times. Panels (b) and (c) show the spatially resolved correlator S z n (t)S z 0 (0) for H3 and H3 + H4 respectively.
thermal in the thermodynamic limit, while the latter allows for the presence of outlying non-thermal states, as long as their ratio is vanishingly small at any given energy [30][31][32]. It is important to stress that if only weak ETH is satisfied, then we can always find initial condi-tions which have narrow energy distributions but nevertheless fail to thermalize [8]. Indeed, generic systems are expected to exhibit the strong version of ETH [12].
Recently, however, several exceptions have been discovered [33,34]. Systems with constrained dynamics [35,36] are an especially promising avenue where non-thermal eigenstates, dubbed scar states, can occur [32,[37][38][39][40]. These are believed to be responsible for persistent oscillations observed in a recent Rydberg atom experiment [41]. Constrained dynamics occurs naturally in so-called fracton systems, which are characterized by the existence of excitations that exhibit restricted mobility [42,43]. On one hand, these systems have been studied in threedimensional exactly-solvable lattice models with discrete symmetries, where fractons are created on the corners of a membrane or fractal operator [43][44][45][46][47][48]. On the other hand, different approaches for fractons with U(1) symmetry have shown that their mobility constraints are related to the conservation of the dipole moment which localizes isolated charged excitations. [49][50][51][52][53][54][55] An analytical connection between the two approaches has been discussed in Refs. 56 and 57. The exotic behavior of fractons also gave rise to the study of their non-equilibrium physics [58,59] and it has been argued that fracton models with discrete symmetries [43,45] show glassy dynamics [58].
In this paper, we study the consequences of dipole conservation associated with a global U(1) charge (i.e., the conservation of total spin S z ) in one-dimensional (1D) spin models for which a numerical study is feasible. Interestingly, a recent work [60] has argued that random local unitary dynamics with such symmetries fails to thermalize. We find the same non-ergodic behavior in a minimal Hamiltonian that contains only three-site interactions. We demonstrate that the non-ergodicity is due to an extensive fragmentation of the Hilbert space into exponentially many disconnected sectors in the local zbasis. In particular, based on the Hilbert space structure, we provide a semi-analytical prediction for the long-time auto-correlation, which remains finite in the thermodynamic limit. This is a novel type of non-ergodic behavior, arising in a translation invariant system, but nevertheless sharing certain features of MBL, which we denote by strong fragmentation of the Hilbert space.
However, we find that this strongly non-ergodic behavior disappears once we add longer-range interactions, such as a four-site term. In this case, the dipole constraint is no longer sufficient to violate ergodicity, and the infinite temperature autocorrelator decays to zero. Nevertheless, the model still violates the strong version of ETH and exhibits exponentially many non-thermal eigenstates, disconnected from the bulk of the spectrum, and co-existing with thermal eigenstates at the same energies. We term this behavior, which is reminiscent to quantum many-body scars, weak fragmentation and give an analytical lower bound on the number of product eigenstates for arbitrary finite range of dipole-conserving interactions. We compare our results to random unitary circuit dynamics, and find the same behavior: while cir-cuits constructed from three-site gates fail to thermalize, adding four-site gates is sufficient to delocalize the system and lead to thermalization for typical initial states. We numerically verify that the invariant subspaces for Hamiltonian and random circuit dynamics coincide exactly.
The remainder of the paper is organized as follows. In Sec. II we introduce the Hamiltonians we study, and describe their relevant symmetries. In Sec. III we investigate the minimal model containing only three-site interactions and show that it fails to thermalize. We prove that the Hilbert space fragments into exponentially many invariant subspaces, some of which we construct analytically, and connect these to the finite saturation value of the auto-correlation function. In Sec. IV we extend the model by adding four-site interactions and argue that while these are sufficient to make the majority of eigenstates thermal-leading to ergodic behavior for typical initial states-the system still violates strong ETH. In Sec. V we compare our results to random unitary circuit dynamics and find similar behavior. We conclude in Sec. VI with a summary and outlook. The appendices provide further comparisons of our numerical results on auto-correlations for different system sizes, as well as other dynamical quantities, such as entanglement growth and operator spreading.

II. MODEL AND SYMMETRIES
We consider two spin-1 Hamiltonians on a chain of length N of the form and acting on three and four consecutive sites, respectively. Apart from being translation and inversion symmetric, both Hamiltonians share two additional global symmetries: they conserve a U(1) charge Q and its associated dipole moment P n0 : with respective eigenvalues q and p defining the symmetry sector H q,p . Since [Q, P n0 ] = 0, the local S z -basis, denoted by |+ , |0 , |− , is a common eigenbasis of Q and P n0 . The definition of the dipole symmetry P n0 also depends on the reference position n 0 , except when Q = 0. Unless specified otherwise, we choose open boundary conditions [61] and take N = 2m + 1 odd, labeling sites n = −m, . . . , 0, . . . , m. We choose the reference site n 0 to be the center site, n 0 = 0, and denote P ≡ P n0=0 . The operator P does not commute with spatial translations and changes sign under inversion; thus, it is not an internal symmetry [55]. Dipole conservation is the relevant global symmetry appearing in the description of fracton phases of matter with U(1) symmetry group [49][50][51][52][53][54][55]. Motivated by this, we use the following notations: we call the states |± on a given site a fracton with charge q = ±1, and a two-site configuration |+− (|−+ ) a dipole with zero charge and dipole moment p = −1 (+1). Notice that the dipole moment of a (±)-fracton on a site n is p = ±n. Thus, in order to conserve the total dipole moment, a fracton can only move by emitting dipoles [49].
There also exists an operator that anti-commutes with H 3 , but commutes with Q and P (see App. A for details). Consequently, the spectrum of H 3 is symmetric around zero in every (q, p)-sector. The same is also true when H 4 is considered separately, but not for the combined Hamiltonian H 3 + H 4 . These anti-commuting symmetries can also be broken by adding terms diagonal in the S z basis, which would not change any of the physics observed in the following.
We note in passing that similar charge and dipole conserving Hamiltonians can be written for any spin representation, in any spatial dimension, as well as for fermionic systems. For the latter, the dipole symmetry becomes the center of mass of the particle number operator and the corresponding Hamiltonian consists of a symmetric redistribution of charges with respect to the center sites. A similar fermionic Hamiltonian appears in the study of fractional quantum Hall on a torus in the Tao-Thouless limit [62][63][64][65][66].

III. HAMILTONIAN H3
We start by investigating the three-site Hamiltonian H 3 in Eq. (1) as a minimal model conserving both total charge Q and dipole moment P . We detail its unusual non-ergodic dynamics and identify it as a consequence of extensive fragmentation of the Hilbert space into invariant subspaces. In Sec. IV we will add longer-range terms to this minimal model and describe their effect on the dynamics.

A. Lack of thermalization
We first investigate the behavior of the autocorrelation function C z 0 (t) ≡ S z 0 (t)S z 0 (0) at infinite temperature in the full Hilbert space. Relying on quantum typicality [67][68][69], we compute C z 0 (t) for random initial states. For thermalizing and translational invariant spin-1 systems, C z 0 (t) is expected to decay to 2/(3N ) for finite chains of length N up to potential boundary contributions [70]. In Figure 1(a) we show C z 0 (t), obtained via an iterative Krylov space based algorithm [71], for a system size N = 15. Instead of relaxing to the thermal expectation value, the auto-correlation saturates to a finite value C z 0 (t)−2/(3N ) ∼ 0.2 at long times. In App. B we confirm that this finite saturation value persists up to large times, t ∼ 10 10 , with no sign of decay. Moreover, the long-time values appear to be largely independent of N , indicating truly localized behavior that persists even in the thermodynamic limit. Figure 1(b) shows the spatially resolved correlation function S z n (t)S z 0 (0) , which exhibits a peak in the center site at all times. We conclude that the system exhibits non-ergodic behavior, similar to the random circuit results of Ref. 60. This is also supported by calculating the growth of entanglement starting from a random product state, which saturates to a sub-thermal von Neumann entropy density as we show in App. C

B. Hilbert space fragmentation
In this section, we demonstrate that the constrained dynamics of H 3 leads to a fragmentation of the manybody Hilbert space: most (q, p) symmetry sectors split into many smaller invariant subspaces in the local S zbasis, such that the total number of such subspaces grows exponentially with system size. These disconnected sectors come in a variety of different sizes; they include 'frozen' states (product eigenstates of H 3 ) and finite dimensional subspaces, where the chain splits into spatially disconnected regions.

Frozen states
We begin by constructing a family of exponentially many exact eigenstates of the Hamiltonian, which are all product states in the local z-basis. We will refer to these as frozen states. The simplest example is the vacuum state |0 ≡ |· · · 0000 · · · , which is annihilated by all terms in H 3 , because (S ± ) 2 |0 = 0. We can easily construct other frozen states by adding blocks of at least two contiguous charges of equal sign on top of the vacuum, e.g., |0 · · · 0 + +0 · · · 0 − − − 0 · · · . These are annihilated by all terms, since S + n S − n+1 |±± = 0. We conclude that any configuration where charges always occur in blocks of at least two consecutive sites are zero energy (mid-spectrum) eigenstates of H 3 . It is clear from the construction that their number is exponentially large in system size.
We can use a Pauling estimate [72] to derive a simple lower bound on total number of frozen states. To do this, we map the spin chain of length N to a triangular ladder with spins placed on the vertices as shown in the inset of Fig. 2(a). Since there are N − 2 triangles and 19 frozen states per triangle, we estimate their total number to be 3 N × 19/27) N −2 ≈ 2.02 × 2.11 N . This is a lower bound of the actual scaling, since constraints on neighboring triangles are not linearly independent as assumed in the Pauling estimate. In Fig. 2(a), we numerically verify that the estimate is quite close to the actual number of frozen states, as obtained by exact diagonalization.

Larger dimensional sectors
Above we saw that blocks of two or more consecutive charges of equal sign are annihilated by the local terms in H 3 that act on them. Let us now consider the empty region (|00 · · · 0 ) between two such frozen blocks and fill it with a random configuration of charges. These charges can now move around and potentially destroy the blocks on the two sides. However, we argue that there are initial configurations where this cannot happen: when the sign of the rightmost charge within the region matches the charge of the frozen block to its right, then this block remains inert at all times. The same holds for the frozen block on the left when its charge is of the same sign as the leftmost charge within the region. When the charges match on both sides, then both blocks are stable and the charges in the middle bounce back-and-forth between them, disconnected from the rest of the chain. This appears as a direct consequence of the general rule: For a region surrounded by empty sites, the signs of the leftand rightmost charges are invariant under the dynamics generated by H 3 .
The simplest example where we can observe this behavior, is as a 2-level system shown in Fig. 2(b), defined by the states |+ + 0 + 0 + + and |+ + + − + + + . We can check that these two states can only evolve to each other under H 3 , defining a small invariant sub-space. More generally we can consider states of the form |+ + 0 · · · 0+0 · · · 0 + + : an isolated fracton surrounded by two 'walls' of positive charge. Acting on this state with H 3 , maps the configuration 00+00 in the middle to 0+ − +0, showing that the (+)-fracton can move by emitting a dipole +− (or −+) in the opposite direction [49]. The fracton can then move forward by emitting further dipoles, until it reaches one of the walls. However, when it eventually hits the wall, we end up with the configuration + + +, which is annitilated by H 3 ; the wall therefore remains intact and the fracton bounces back harmlessly, as shown in Fig. 2(c). To destroy the wall, we would need to flip the charge of the isolated fracton to get a (−)-fracton: the resulting − + + configuration can then peel off a freely-moving −+ dipole, eventually melting the walls that surround it as shown in Fig. 2

(d).
A similar situation occurs for the initial configuration |− − 0 · · · 0−+0 · · · 0 + + . In this case the walls on the two sides have opposite charges and a single dipole is placed between them. For a single dipole surrounded by empty sites, the Hamiltonian H 3 acts as a free hopping term, moving the dipole from site n to n ± 1 [73]. Eventually it reaches one of the surrounding walls, but since the charges in the dipole are aligned with those of the walls, it always bounces back, effectively defining a single particle hopping problem on a finite region. If, on the other hand, the initial dipole in the middle was of the form +− it could again peel off charges from the two walls, eventually melting them.
The previously stated general rule, together with the fact that blocks with a given charge are frozen, allow us to construct more general spatially disconnected regions in the chain: take an arbitrary configuration in some finite interval and surround it with walls that have the same charge as the one closest to them on the inside. One can then cover the entire chain with such regions, giving rise to many invariant subspaces within each global (q, p) symmetry sector. The resulting eigenstates clearly break translation invariance and have small amounts of entanglement, limited by the size of the largest connected spatial region.
These constructions highlight the intertwined relation between dipole conservation, spatial translations, and locality. After the dipole quantum number is fixed, translation/inversion symmetry is generically broken, which allows us to derive conservation laws within different spatially disconnected regions. Our construction also shows that in order to determine which invariant subspace a given initial configuration belongs to, one has to consider it on the entire chain: even if a certain region looks initially frozen, it can eventually be melted by additional charges coming from the outside. This indicates that it might not be possible to systematically label all invariant subspaces in terms of quantum numbers of local conserved quantities.

Distribution of dimensions of invariant subspaces
Above we explicitly constructed invariant subspaces of H 3 of various dimensions within given (q, p) symmetry sectors. The distribution of these invariant subspaces can be studied by numerically identifying the connected components of the Hamiltonian written in the S z basis. The resulting distribution is plotted in Fig. 3(a), showing exponentially many sectors with a broad distribution. We point out that since the sectors are obtained in the local z-basis, they remain invariant under any perturbation that is diagonal in this basis. However, such diagonal perturbations would have the effect of changing the energy of the different frozen states, moving them away from zero energy, and distributing them throughout the entire spectrum.
Based on the constructions in the previous section, we infer that the existence of these invariant subspaces is a consequence of the interplay between the conservation of dipole moment (which fails to commute with translation and inversion) and the locality of interactions. In particular, in Sec. IV B we prove that exponentially many invariant subspaces exist for any extension of the model that only involves dipole-conserving interactions with finite range.
Equipped with the knowledge of the fragmented Hilbert space structure, we are now able to explain the long-time value of the auto-correlation function observed in Fig. 1(a). To this end, we can rewrite the autocorrelation evaluated in the full Hilbert space as where the sum is performed over all invariant subspaces H i within all (q, p) sectors, and Z i (t) is the projection FIG. 4. Saturation value of the autocorrelator. Finitesize study of the semi-analytical prediction (5) for C z 0 (t → ∞) as a function of system size. We have substracted the ther- S and a chain of length N . For the minimal spin-1 model H3 (blue dots) the prediction is slightly increasing with system size. On the other hand, it decays to zero exponentially for the combined Hamiltonian H3 + H4 (blue squares). For comparison, we also show results for other local spin S: the larger the on-site Hilbert space dimension, the easier it is for the system to thermalize [74].
of the operator S z 0 (t) onto H i . Approximating the longtime evolution within each subspace H i of dimension D i by independent Haar random unitaries, we obtain This prediction is shown in Figure 1(a) by the black dashed horizontal line. Our estimation works quantitatively well, giving a slightly lower value than the one observed numerically. We compute this predicted value for C z 0 (∞) − 2/3N (shown in Fig. 4) for a variety of different system sizes. This value appears to remain finite in the thermodynamic limit, even increasing slightly with N for the system sizes available in our numerics (blue dots).
Since the H i 's are invariant and disjoint subspaces, the weight of the operator S z 0 within a given sector, tr Z 2 i , remains constant under time evolution. Therefore, we introduce the operator weight W D ≡ 2 as a function of the sector size D for all invariant subspaces H i . This defines a probability distribution, shown in Fig. 3(b). We find a wide distribution with significant weight on small sectors. While the number of frozen states scales as ∼ 2.2 N , the size of the largest sector in the entire Hilbert space scales as ∼ 1.9 N , both much smaller than the total dimension 3 N . This suggests that sectors of all sizes have significant contributions to the evolution of S z 0 (t), even in the thermodynamic limit. We also confirm the same behavior when considering only the largest symmetry sector (q, p) = (0, 0); this emphasizes the relevance of the fragmentation within each (q, p)-sector.
IV. COMBINED HAMILTONIAN H3 + H4 So far we have only considered the 'minimal model', defined by the Hamiltonian H 3 in Eq. (1). We will now investigate to which extent the features found above are robust against local perturbations that preserve the symmetries Q and P .
A. Thermalization for H3 + H4 In the following, we add the four-site terms defined in Eq. (2) and consider the combined Hamiltonian H 3 +H 4 . We find that, while this Hamiltonian shares certain features with H 3 -in particular, it has exponentially many invariant subspaces-it nevertheless thermalizes at infinite temperature. Indeed, the auto-correlation function C z 0 (t) for the Hamiltonian H 3 +H 4 decays to zero at long times, in contrast to the dynamics governed by H 3 alone; see Fig. 1 for a comparison. This is accompanied by the spatially resolved correlation function, S z n (t)S z 0 (0) , becoming approximately homogeneous at long times, as shown in Fig. 1(c). The remaining small peak in the autocorrelator is due to finite size effects, as we show in App. B. Moreover, as we discuss in App. C, for a random product state evolving under H 3 + H 4 , the entanglement entropy approaches its thermal value at long times, providing an additional indication that the system thermalizes.
This qualitative difference suggests that the Hilbert space structure uncovered in Sec. III B should also be modified by adding H 4 to the Hamiltonian. Figure 5(a) compares the distribution of sector sizes D for H 3 + H 4 (blue stars) with the minimal Hamiltonian H 3 (red dots). While exponentially many invariant subspaces still exist, their total number is drastically reduced, as many previously disconnected sectors are coupled to each other by the perturbation H 4 [75]. Thus the number of sectors of small dimension D decreases and there are new larger blocks appearing; in fact, the largest global symmetry sector, q = p = 0, becomes almost (but not exactly) fully connected, as we discuss in Sec. IV C. This effect is even more apparent in the distribution of the operator weight W D (defined in Sec. III C) for the operator S z 0 , which we show in Fig. 5(b). Most of the weight is now concentrated around the largest sector, similarly to the case of a single global U(1) symmetry. Thus, even though invariant subspaces within symmetry sectors still exist, they do not appear to be sufficiently relevant to make the system non-ergodic. This is also reflected in the long-time value of the auto-correlation function as predicted in Eq. (5): plugging in the invariant subspaces of H 3 + H 4 we find that C z 0 (∞) approaches the thermal value, 2/(3N ), exponentially in the thermodynamic limit, as shown in Fig. 4.
From these results we infer that allowing for longerrange interactions, makes the system sufficiently ergodic to thermalize. A different path to break the nonergodicity of H 3 , would be to increase the local Hilbert space dimension, making the dynamics less constrained. Consequently, we expect that for larger spin, even a three-site Hamiltonian of the form (1) would lead to thermalization. Indeed, computing the long-time prediction of the autocorrelator C z 0 (∞) using Eq. (5) for H 3 acting on a spin-2 chain, we find that it decays to zero in the thermodynamic limit, as shown in Fig. 4. Similarly, if we consider spin-1/2 chains, the shortest range nontrivial model is H 4 , which appears to be non-ergodic, while adding 5-site interactions restores ergodicity.

B. Contructing frozen states
While the combined Hamiltonian H 3 + H 4 appears thermalizing at infinite temperature, it nonetheless violates the strong version of the Eigenstate Thermalization Hypothesis [9][10][11][12]. In particular, certain frozen states continue to exist not only for H 3 + H 4 , but even in the presence of longer-range local interactions. In fact, as we now prove, for a spin-1 chain that conserves charge and dipole, and involves only local terms with range at most , there exist at least 2 · 5 N/ frozen states. While this lower bound is not as tight as the Pauling estimation discussed in Sec.III B 1, it provides useful insights into longer-range Hamiltonians and can be generalized to any spin representation. We begin our construction by considering the configuration shown in Fig. 6(a), with a center site surrounded by a block of − 1 (+)-fractons on one side and − 1 (−)-fractons on the other. We now prove that this configuration is an eigenstate of any dipole-conserving term with range at most , where without loss of generality we can measure the dipole moment relative to the center site. It is sufficient to consider off-diagonal terms (in the z basis), consisting of spin raising and lowering operators. Due to the way we constructed the state, the only such terms that do not annihilate it are those that have only S − on one side and S + on the other. However, any such term would lead to a change in the dipole moment and is thus prohibited. Terms only acting on the center site do not change P but they are also excluded due to charge-conservation. We conclude that this configuration is frozen, independently of the state of the center spin, as promised. Next, we consider a similar configuration, but one where the central spin is surrounded by blocks of the same, rather than opposite, charges, as shown in Fig. 6(b). Let these blocks be made out of (+)-fractons. Then the only off-diagonal operators that can act on them are powers of S − , decreasing the total charge Q. One has to compensate for these charges by adding additional charges on the center site. Therefore the only allowed terms that could change this configuration are of the form S − −n (S + 0 ) 2 S − n , and only when the central spin is occupied by a (−)-fracton. When it is either 0 or +, the state is frozen.
We can then combine these two types of 'frozen patches' we constructed above to cover the entire 1D chain, resulting in a globally frozen state. These states are made up by blocks of + or − charges, with a single site between any two consecutive blocks, as shown in Fig. 6(c). As we showed above, these sites host flippable spins: the ones separating blocks of equal (±)-charge can be ± or 0, while those that separate blocks of opposite charge can be in any of the 3 possible spin states. This construction then results in exponentially many frozen states, coming from both the possible arrangements of ± blocks and from flipping the spins between blocks within a given arrangement.
We can count the total number of frozen states resulting from this construction iteratively, starting from the left edge of the system (assuming open boundaries). We cut the systems into blocks of sites, consisting of a wall of − 1 positive/negative charges, followed by a flippable spin. Let F ± k denote the number of different such configurations to the left of the k-th wall (but before the flippable spin), ending in a (±)-block. Then the considerations outlined above lead to the following recursion formula:

FIG. 7.
Scaling of the number of frozen states for Hamiltonians with at most range terms. We consider Hamiltonians with all possible combinations of charge and dipole conserving quartic in spin operators of range at most for = 3, 4, 5. The number of frozen states is exponential, with an exponent that decreases with , but is larger then the analytical lower bound 2 · 5 N/l . that the number of frozen states we constructed scales as 2·5 N/ . This is only a lower bound on the total number of frozen states, which can include other configurations not captured by this construction. We compare this bound to the numerical results of the number of frozen states for Hamiltonians with interactions of range = 3, 4, 5 in Fig. 7, where we extract the asymptotic scaling. The comparison to the numerical data shown in Fig. 7  We conclude this section with some comments about the construction we presented. First, while above we did not distinguish between different overall (q, p)-sectors, one could similarly construct frozen states with a given q and p. For example one can apply the construction on only the left half of the chain and for each state repeat the same configuration on the right half to obtain a state with p = 0. Second, the bound can be easily extended to chains with local spin S > 1. For example one can consider blocks that have maximal positive/negative charge. Repeating the same arguments then gives a scaling [76] (2S + 3) N/ . Last, we note that in the limit → ∞ the lower bound tends to one, consistent with the fact that for all possible charge and dipole conserving infinite-range interactions every (q, p) sector becomes completely connected.

C. Strong vs. weak fragmentation
As the previous section shows, the combination of dipole conservation and strictly local interactions is sufficient to lead to an emergence of exponentially many dy- namically disconnected sectors in the many-body Hilbert space, even after fixing q and p. While we only showed this rigorously for the case of one-dimensional sectors, we find numerically that others with larger dimension also exist (see Fig 5). While both H 3 and H 3 + H 4 share this feature, their behavior with respect to thermalization appears to be quite different, as we already observed in Fig 1. This motivates us to distinguish two cases, dubbed weak and strong fragmentation, which violate strong and weak ETH, respectively.
Let us first make precise what we mean by violation of ETH. In defining ETH we consider expectation values of few-body observables for all eigentates of the Hamiltonian within a fixed global (q, p) symmetry sector (we do not consider off-diagonal matrix elements here). By strong ETH we then mean the statement that the expectation values are the same for all eigenstates at the same energy density in the thermodynamic limit. Weak ETH, on the other hand, means that this statement only holds up to a small number of outlying states, where 'small' means here 'measure zero in the thermodynamic limit'.
Our construction in the previous section then proves that any dipole-conserving, strictly local Hamiltonian has weak fragmentation in the above sense, i.e., non-thermal eigenstates are present in the middle of the spectrum. Apart from the aforementioned frozen states, these also include other low entanglement eigenstates, stemming from small invariant subspaces, analogous to the ones discussed in Sec. III B 2. Generically, however, their ratio compared to thermal ones is vanishingly small within any energy shell in the thermodynamic limit; this is the case of H 3 + H 4 as we argue below. Thus the weak version of ETH [30][31][32] is still obeyed, and the system thermalizes for typical initial states, provided they have narrow energy distributions. On the other hand, we argue that the Hamiltonian H 3 , discussed in Sec. III, has strong fragmentation in the sense that at least a finite fraction of the eigenstates is non-thermal, leading to the manifestly non-thermalizing behavior we observed.
The difference is illustrated in Figs. 8(a,b), which shows the expectation value of a simple observable ((S z 0 ) 2 where 0 is the central spin) for all energy eigenstates within the q = p = 0 symmetry sector, for H 3 and H 3 + H 4 . For the combined Hamiltonian H 3 + H 4 , the majority of eigenstates, which all belong to the same invariant subspace, behave as predicted by ETH: (S z 0 ) 2 takes similar values for all states within a narrow energy shell, with the width of its distribution decreasing with system size. Nevertheless, we also observe outlying eigenstates, stemming from small invariant subspaces, that do not approach this line, violating strong ETH. The minimal Hamiltonian, H 3 , on the other hand, violates even the weak version of ETH: the distribution of (S z 0 ) 2 does not become narrower with increasing N , as shown in Fig. 8(a). This is in contradiction with the ETH, which predicts a vanishing width in the thermodynamic limit. Similar behavior occurs in the half-chain entanglement entropy of the eigenstates, shown in Fig. 8(c): the nonthermalizing nature of H 3 is reflected by the fact that the entropies of its eigenstates do not fall on a line when plotted as a function of the energy, instead being distributed over values much smaller than what is predicted at infinite temperature, realized by a random state in the (0, 0) sector.
The above discussion suggest that the difference between strong and weak fragmentation can be diagnosed by considering the sizes of the connected subspaces, in comparison with the size of the global (q, p) symmetry sector they belong to. In the strongly fragmented case of H 3 studied above, for a typical (q, p) symmetry sector, the dimension of the largest connected subspace is exponentially smaller than the dimension of the full symmetry sector, i.e., max[D i (q,p) ]/D (q,p) ∝ exp(−αN ) for some α > 0. In Fig. 9 we verify that this is indeed the case for the largest symmetry sector (0, 0) of H 3 . We propose that this decay indicates strong fragmentation, naturally leading to the absence of thermalization for physical observables such as the auto-correlation function considered above.
In the weakly fragmented case, the symmetry sectors can still split into many subspaces. However, the largest of these spans almost the entire (q, p)-sector: , with the ratio approaching 1 in the thermodynamic limit. Figure 9 shows that this is the case for H 3 + H 4 . Consequently, the vast majority of eigenstates within any energy shell in a given (q, p) symmetry sector belong to the same large invariant subspace, and look thermal as a consequence. Thus, while weakly fragmented systems violate strong ETH-due to outlying non-thermal eigenstates-they nevertheless thermalize for typical (but not all) initial states. This weak fragmentation is reminiscent to what has been observed in other models in the context of many-body quantum scars: while the majority of the eigenstates obey ETH, non-thermal eigenstates exist even in the bulk of the spectrum [32-34, 37, 38]. However, while the number of these 'scarred' states is usually O(N ), in our case we find exponentially many such states. Note that the non-thermal eigenstates belonging to low-dimensional invariant sectors in our system have finite overlap with simple product states which can potentially be prepared in experimental settings. This implies that a lack of thermalization up to infinite times could be observed, even in the weakly fragmented case, for appropriately chosen initial states.
In the above we observed that the ratio max[D i (q,p) ]/D (q,p) either decays (exponentially) to zero or approaches unity. It is an interesting and open question, whether systems with intermediate behaviorwith either slower than exponential decay or convergence to a finite fraction-can exist, and whether they exhibit strong or weak fragmentation.

V. COMPARISON TO RANDOM UNITARY CIRCUITS
We now argue that our findings are not specific to the Hamiltonians we considered so far, and generalize to arbitrary systems with the same global symmetries and a fixed range of interactions. In particular, we compare with random unitary circuits of the form originally introduced in Ref. 60. These define a discrete time evolution, the building blocks of which are unitary gates acting on sites, each of which is required to be block diagonal in Q and P , but is otherwise chosen randomly (i.e. every block is independently Haar random). In Ref. 60 it was argued that such circuits always lead to localized behavior. Here we argue that this is in fact only the case for gates with = 3, where the circuit exhibits exactly the same Hilbert space structure as the Hamiltonian H 3 above, and is therefore indeed similarly localized. When introducing larger gates of size = 4 we find that the system becomes delocalized, also in complete agreement with our results on the Hamiltonian H 3 + H 4 .
The two circuit geometries, with gates of size = 3 and 4, are shown in Fig. 10. In both cases we compute the connectivity of the Hilbert space. Instead of the Hamiltonian, we consider the unitary operator defined by the first layers of the circuit. This is a matrix with random entries, but its connected components are independent of the particular realizations. We find numerically that the connected components for = 3 ( = 4) coincide exactly with those of the Hamiltonians H 3 (H 3 + H 4 ), shown previously in Fig. 5. This follows from the fact that the allowed local transitions are the same in the Hamiltonian and the random unitary circuit. The fact that the invariant subspaces coincide supports the idea that the additional invariant subspaces are a consequence of dipole conservation and locality alone, and do not depend on any additional structure that might be present in the Hamiltonian case. Based on our previous analysis, we therefore expect that the three-site circuit does not thermalize, but the four-site circuit does. This is confirmed by calculating the autocorrelator C z 0 (t), which (after subtracting its thermal value) goes to a constant in the former case, while it decays to zero in the latter, as shown in Fig. 10. In App. D we also consider the spatial spreading of an initial S z n operator and similarly find that for = 4 the operator is delocalized at long times. We therefore conclude that the localized behavior observed in Ref. 60 is particular to the case of the circuit with three-site gates, contrary to what is suggested there.
The fact that the Hilbert space fragmentation coincides exactly between the random circuit and Hamiltonian cases also means that the conclusions we drew regarding non-thermal eigenstates in Sec. IV C also generalize to time-periodic (Floquet) models built out of similar local gates. In particular this implies the presence of exponentially many frozen eigenstates for such models, which appears to contradict the results of Ref. 77, where a polynomial number or scar states were claimed, even for the = 3 case where we predict that the majority of eigenstates should be non-thermal.

VI. SUMMARY AND OUTLOOK
In this work, we studied the out-of-equilibrium dynamics of spin chains conserving a charge and its associated dipole moment. For the minimal spin-1 Hamiltonian which is restricted to only three-site interactions, we found non-ergodic behavior in the charge-charge autocorrelation function. We explained this finding in terms of a strong fragmentation of the Hilbert space into exponentially many disconnected sectors which all contribute significantly to the dynamics even at infinite temperature. We found that a weaker form of fragmentation survives for more general, longer-range Hamiltonians, which is no longer sufficient to make the infinite-temperature dynamics non-ergodic though. Furthermore, we showed numerically that the fragmentation of the Hilbert space exactly matches that of random circuit dynamics with the same range of interactions, giving rise to similar dynamical behavior.
The observed fragmentation lies in-between the known cases of systems with a few global symmetries and that of integrable or many-body localized systems. The former have at most polynomially many symmetry sectors -most of which are exponentially large-while the latter have ∼ N independent local conserved quantities. In our case, however, sectors of all sizes co-exist and in the case of strong fragmentation they all are relevant for the dynamics, even at infinite temperature. Understanding how these different sector can be consistently labeled and what the corresponding operators look like is an interesting open problem.
It is also an interesting problem for future work, to investigate the equilibrium properties and dynamics of the system at finite or zero temperature, as well as to clarify the requirements for strong and weak fragmentation. Moreover, the extension of the current analysis to higher spatial dimensions seems promising, where disconnected sectors also appear but different roads to thermalization can be present. It is also worth exploring the connections to: Stark localization [28,29], many-body localization [21], quantum scars [37] or gauged formulations of the considered systems [50,78]. Finally, it would be interesting to search for connections between the models discussed in this work and type I fracton models [43,45], as well as for the possible connection between the discussed fragmentation of the Hilbert space and the emergence of superselection rules in the cubic code model [79]. Note added. While we were finalizing our manuscript, we have learned about related work by V. Khemani and R. Nandkishore which will appear in the same arXiv posting.
Appendix A: Symmetries of the 'minimal' Hamiltonian H3 Here we discuss some additional symmetries possessed by the 'minimal model' represented by the Hamiltonian H 3 in Eq. (1). One of these is the sublattice parity symmetry Π z odd . However, one can easily prove that this quantity is fixed by the dipole moment as where in the second step we have used that for spin-1, a 2π rotation is equal to the identity. From this it is clear that the total parity Π z = exp iπ n S z n is obtained as Π z = Π z odd Π z even and is related to the total charge as Π z = exp iπQ . In general, the terms in H 3 are also invariant under the parity transformations given by Π x = exp iπ n S x n and Π y = exp iπ n S y n , which map for all n 1 , n 2 , n 3 , n 4 . Note that Π x and Π y do not commute with Q or P .
Moreover, as stated in the main text, there exists an operator C = n e iπ S z 4n +S z 4n+1 which anti-commutes with H 3 . Since C commutes with Q and P , the spectrum of H 3 is symmetric around zero in every (q, p)-sector. However, when additional terms diagonal in the S z basis are considered, which by construction do not change the fragmentation of the Hilbert space as we discussed in the main text, C does not anti-commute anymore with the resulting Hamiltonian.
There is also at least one additional anti-commuting operatorC = n e iπ S z 4n+2 +S z 4n+3 but since CC = Π z , they are not independent. Note that since C commutes rather than anti-commutes with H 4 , the spectrum of H 3 + H 4 is no longer symmetric as can be seen e.g., in Fig. 8. Nevertheless, there also exists a separate anti-commuting operator for H 4 (2) taking the form C 4 = n e πS z 4n , which does not anti-commute with H 3 . In Fig. 11(a) we show the density of states, ρ(E), of H 3 for a chain of length N = 13, which has a divergent delta peak at zero energy. This peak contains all the frozen states described in the main text, among other zero energy eigenstates that arise as a consequence of the aforementioned anti-commuting symmetry. One could remove the peak at zero energy by adding e.g. a finite mass term of the form m n (S z n ) 2 to H 3 , which breaks the anti-commuting symmetry. This term also has the effect of shifting the energy of the frozen states to finite values and distributing them throughout the spectrum.
In Fig. 11(b) we show the size of the symmetry sectors with different global quantum numbers q and p. Note that this distribution is independent of the specific Hamiltonian under study. Each curve corresponds to a fixed value of the charge quantum number q. The dimension D (q,p) decreases with increasing absolute value of the charge. The distributions for +q and −q coincide due to time reversal invariance, the way we have chosen the reference site n 0 , and labeling the sites in the chain. A different labeling of sites, would simply shift the mean value of both distributions symmetrically with respect p = 0. We also observe that the distribution attains a maximum at the (0, 0)-sector, as claimed in the main text. In addition, we obtain symmetric distributions because P changes sign under inversion, while Q is invariant.
Finally, we show the sector size and the operator weight distributions with W D ≡ for invariant subspaces within the largest (q, p)-sector, q = p = 0. Fig. 12(a) shows qualitatively the same sector size distribution as in Fig. 2(c) in the full Hilbert space. Fig. 12(b) also reflects the main properties of the operator weight distribution, featuring a wide distribution with significant weight on small sectors.  In this section we present in more detail the finite size scaling of the auto-correlation function and the semi-analytical prediction.
First, we discuss the scaling of the auto-correlation function S z 0 (t)S z 0 at infinite temperature in the full Hilbert space in Fig. 13 for both the minimal model H 3 in Eq. (1) and the combined Hamiltonian H 3 + H 4 . On the one hand, the minimal model realizes a finite saturation value at long times which slightly grows with system size as can be seen in Fig. 13(a). On the other hand, when the combined Hamiltonian H 3 + H 4 is considered, the auto-correlation decays to zero with system size. This agrees with the discussion in the main text, where it was argued that for longer range Hamiltonians, the system thermalizes and the correlation decays to zero at long times in the thermodynamic limit.
Moreover, as we discussed in the main text, not only the auto-correlation function in the full Hilbert space shows a non-thermal (thermal) behavior for H 3 (H 3 + H 4 ). We can also realize this behavior within a specific restricted symmetry sector. In Figs. 14(a) we show the behavior of the auto-correlation function S z 0 (t)S z 0 (0,0) in the largest (q, p)-symmetry sector, q = p = 0, and size N = 15 showing the same qualitative behavior: a finite saturation value at long times for H 3 (panel (a)) and thermalization for the combined Hamiltonian H 3 + H 4 (panel (b)). Note that since charge is conserved and we evaluate the correlation within the q = 0 sector, n S z n (t)S z 0 (0) = 0 at all times and thus the surface under the peak must add up to zero.
Moreover, in Fig. 15, we show the persistence of the non-thermalizing behavior for H 3 at longer times t = 10 10 for smaller system size N = 13 and within the (0, 0)-sector. The space resolved correlation function is also shown in the inset showing the absence of thermalization even at long time scales.
In Fig. 16(a) we show the scaling of the prediction lim t→∞ S z 0 (t)S z 0 (0) (0,0) in Eq. (5) with system size N restricted  to the (0, 0) symmetry sector of H 3 . In this case the prediction takes the form We observe that the value increases with N and realize an even-odd dependence on N decreasing with system size. In addition, Figs. 16(b-c) show, respectively, how the number of frozen states and the size of the largest invariant subspace within the (0, 0) symmetry sector grow with system size. Since the largest sector does not scale with the size of the entire Hilbert space, the lower dimensional sectors become thermodynamically important. Compare for example with a spin 1/2 chain with charge conservation only. The dimension of the full Hilbert space is 2 N and the largest (zero charge) subspace scales as 1/N · 2 N ; hence, the exponents are the same up to logarithmic corrections. FIG. 17. Half chain entanglement entropy (EE) growth for an initial random product state for different system sizes. The dashed line signals the Page value [80]. Panel (a) shows the behavior of the EE for the minimal model H3. The entanglement reaches a size dependent saturation value below the Page value. This can be understood from the exponential fragmentation of the Hilbert space. However, when the combined Hamiltonian H3 + H4 is considered in panel (b), the EE almost reaches the Page value. We associate the offset between them to the existence of still exponentially many invariant subspaces. (c) Scaling of the time-averaged saturation value for EE reached at long times. While for H3 the slope is different from that of the Page value, signalling a sub-thermal entropy density in the steady state, the off-set for the combined Hamiltonian appears to be constant.

Appendix C: Entanglement growth from random product states
In this appendix we complement our results on auto-correlations with a different measure of thermalization: entanglement growth from a (random) product state. By choosing the initial state Haar randomly on each site and averaging, we ensure that the dynamics explores all (q, p) symmetry sectors. For an ergodic system, the long-time state is then expected to resemble a global random state in the entire Hilbert space. In particular, the entanglement between two halves of a bi-partition is expected to be given by the Page formula, which in our case (maximal bi-partition of a spin-1 chain with odd lengths) reads [80] We evaluate the time evolution starting from the aforementioned random product states exactly, for both the minimal Hamiltonian H 3 and the combined Hamiltonian H 3 + H 4 . In the former case, shown in Fig. 17(a), we find that while the entanglement quickly saturates to a volume law, the associated entropy density is smaller then the expected Page value, indicating a non-thermal state. This is consistent with our results on auto-correlation functions in Fig. 1, as well the entanglement of eigenstates in Fig. 8, all consistent with non-ergodic behavior. The entanglement growth for H 3 + H 4 is shown in 17(b) where we observe that the entanglement saturates to a value close to S Page . There is still a constant offset, which we associate to the influence of the remaining non-thermal eigenstates, but this does not affect the entropy density in the thermodynamic limit.
Appendix D: Operator spreading of S z 0 (t) Here we consider another measure of localization, that contains complementary information about the Heisenberg picture evolution of the charge density operator S z 0 (t) compared to its auto-correlation function. In particular we look at how S z 0 spreads out in the space of all possible operators, becoming a complicated superposition of many operators, and how its spatial support grows in time.
In order to do this, we first need to introduce a local basis in the space of operators acting on a single site of the spin chain. For the spin-1 models we consider, such a basis constists of 9 linearly independent operators that span the entire space of on-site operators. A possible choice is given by the 8 Gell-Mann matrices, together with the identity 1 1. Let us denote these as λ a for a = 0, . . . , 8, where λ 0 ≡ 1 1. A basis of operators on the entire chain is then given by products of such local basis elements of the form λ a ≡ N/2 n=−N/2 λ an n , labeled by a list of N indices {a n }. These operator string form an orthonormal basis in the Hilbert space of operators with respect to the Frobenius inner product A, B ≡ tr(A † B)/3 L where A and B are two arbitrary operators.
Given such a basis, one can always expand the time evolved operator as The coefficients c a (t) characterize how S z 0 (t) spreads out in the space of all possible operators. In particular, focusing on spatial spreading, it is useful to classify the basis strings λ a according to their right endpoints (assuming open boundary conditions), i.e., the rightmost site n such that λ an = 1 1 but λ am>n = 1 1. Denoting this site by RHS( a) we can define the right endpoint density of S z 0 at time t as [81][82][83] ρ R (n, t) ≡
At time t = 0 this is a delta function at the initial position of the operator, ρ R (n, 0) = δ n0 . During time evolution, as the support of S z 0 (t) increases, ρ R (n, t) moves to the right, ballistically for generic clean systems. At the same time, its value near the origin decays to zero, exponentially when symmetries are not present [82], and as a power law when the operator is a conserved density [84,85]. A possible alternative measure of localized behavior is therfore to look at the spreading of the right endpoint density and look for a finite weight remaining near the origin at infinite times, even in the thermodynamic limit.
We first consider the evolution of ρ R (t) in random circuits, first with 3-and then with 4-site gates. In order to evaluate ρ R (n, t) we represent S z 0 (t) as a matrix product operator [86] (MPO) and apply the random gates to that to evolve it in time. In order to simplify the calculations, we consider slightly modified circuit geometries, which allow us to use the well known time-evolving block decimation (TEBD) algorithm, after blocking pairs of sites together [87].
Our numerics only allow us to access small systems of size N = 6, 8, 10. To compute the spreading of ρ R (n, t), we place an operator S z on the third site from the left end of the system and calculate ρ R (n, t) at different positions and times. For a circuit made out of 3-site gates, we find a persistent peak near the original position, whose size decays only slightly with system size (Fig. 18(a)-(c)). For the circuit with gate-size l = 4 on the other hand, we observe a much smaller peak, which keeps decreasing until finite size effects kick in, similar to the behavior observed for the autocorrelator in the main text, and consistent with the perdiction that in the thermodynamic limit the peak would eventually disappear (Fig. 18(d)-(f)). We also observe a larger peak at the rightmost site, where most of the operator weight accumulates at long times.
The same difference in behavior between 3-site and 4-site interactions is also present in the Hamiltonian case. For H 3 we find that the peak in ρ R (0, t) is almost independent of system size, in agreement with the non-ergodic behavior observed in the autocorrelator in the main text. This is shown in the left panel of Fig. 19. This behavior changes, however, once we add 4-site terms to the Hamiltonian. In particular we consider the perturbation This is the same as in Eq. (2), except that only terms with even n are present. This is done in order to simplify numerical calculations (making the Hamiltonian nearest neighbor after blocking pairs of neighboring sites together). We expect that if H 3 + H 4 does not exhibit a presistent peak in ρ R , then neither should H 3 + H 4 , therefore it is enough to show its absence in the former case. This is indeed what we find as shown in the right panel of Fig. 19