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

TL;DR

This paper demonstrates that charge and dipole conservation in fracton systems cause Hilbert space fragmentation, leading to non-thermal behavior and ergodicity breaking in one-dimensional spin models, with implications for understanding thermalization failure.

Contribution

It reveals how charge and dipole conservation induce Hilbert space fragmentation and ergodicity breaking in dipole-conserving Hamiltonians, supported by concrete models and comparisons to random circuit systems.

Findings

Finite auto-correlation saturates at finite value, indicating non-thermal behavior.

Strong fragmentation into exponentially many invariant subspaces prevents thermalization.

Weak fragmentation persists with four-site terms, but allows thermalization for typical states.

Abstract

We show that the combination of charge and dipole conservation---characteristic of fracton systems---leads to an extensive fragmentation of the Hilbert space, which in turn can lead to a breakdown of thermalization. As a concrete example, we investigate the out-of-equilibrium dynamics of one-dimensional spin-1 models that conserve charge (total $S^z$) and its associated dipole moment. First, we consider a minimal model including only three-site terms and find that the infinite temperature auto-correlation saturates to a finite value---showcasing non-thermal behavior. The absence of thermalization is identified as a consequence of the strong fragmentation of the Hilbert space into exponentially many invariant subspaces in the local $S^z$ basis, arising from the interplay of dipole conservation and local interactions. Second, we extend the model by including four-site terms and find that this perturbation leads to a weak fragmentation: the system still has exponentially many invariant subspaces, but they are no longer sufficient to avoid thermalization for typical initial states. More generally, for any finite range of interactions, the system still exhibits non-thermal eigenstates appearing throughout the entire spectrum. We compare our results to charge and dipole moment conserving random unitary circuit models for which we reach identical conclusions.

Equations

$$H_3=$$

$$H_4=$$

$$Q\equiv \sum _n{S}_n^z\text{and}$$

$$\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^{T}dt⟨S_0^z(t)S_0^z(0)⟩$$

$$\left(\begin{array}{c}{F}_{k+1}^{+}\\ F_{k+1}^{-}\end{array}\right)=\left(\begin{array}{cc}2& 3\\ 3& 2\end{array}\right)\left(\begin{array}{c}{F}_k^{+}\\ F_k^{-}\end{array}\right)={\left(\begin{array}{cc}2& 3\\ 3& 2\end{array}\right)}^k\left(\begin{array}{c}1\\ 1\end{array}\right)=5^k\left(\begin{array}{c}1\\ 1\end{array}\right),$$

$$\displaystyle\Pi^{z}_{\text{odd}}=\exp\big{(}{i\pi\sum_{n\,\mathrm{odd}}{S_{n}^{z}}}\big{)}=\exp\big{(}{i\pi\sum_{n}{nS_{n}^{z}}}\big{)}=\exp\big{(}{i\pi P}\big{)},$$

$$S_{n_1}^{+}{S}_{n_2}^{-}{S}_{n_3}^{-}{S}_{n_4}^{+}\stackrel{{\Pi }^x,{\Pi }^y}{⟷}{S}_{n_1}^{-}{S}_{n_2}^{+}{S}_{n_3}^{+}{S}_{n_4}^{-}$$

$$\displaystyle\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}dt\,\left<S_{0}^{z}(t)S_{0}^{z}(0)\right>_{(0,0)}\geq\frac{1}{D_{(0,0)}}\sum_{\mathcal{H}_{i}\subset\mathcal{H}_{(q,p)}}{\frac{1}{D_{i}}\left[\mathrm{tr}\big{(}Z_{i}\big{)}\right]^{2}}.$$

$$H_{\text{PXP}}=\sum _{n=1}^L{P}_{n-1}{\sigma }_n^x{P}_{n+1},$$

$$|{\mathbb{F}}_2⟩=|\dots \boxed{-+-}-\stackrel{2(k-1)}{\boxed{-+-}}\stackrel{2k}{-}\stackrel{2(k+1)}{\boxed{-+-}}-\boxed{-+-}\dots ⟩,$$

$$h_{n}\equiv S_{n-1}^{+}\big{(}S_{n}^{-}\big{)}^{2}S_{n+1}^{+}+\text{H.c.},$$

$$|{\mathbb{F}}_2⟩⟶h_{2(k-1)}{h}_{2(k+1)}|{\mathbb{F}}_2⟩=|\dots \boxed{-+-}-0\stackrel{2(k-1)}{-}\stackrel{2k}{\boxed{0-0}}\stackrel{2(k+1)}{-}0-\boxed{-+-}\dots ⟩.$$

$$h_{2(k-1)}{h}_{2(k+1)}|{\mathbb{F}}_2⟩⟶h_{2k}{h}_{2(k-1)}{h}_{2(k+1)}|{\mathbb{F}}_2⟩=|\dots \boxed{-+-}-0\stackrel{2(k-1)}{-}\stackrel{2k}{\boxed{-+-}}\stackrel{2(k+1)}{-}0-\boxed{-+-}\dots ⟩.$$

$${H_3|}_{\text{PXP}}=\sum _{k=1}^{L/2}{h}_{2k}=4\sum _{k=1}^{L/2}|-{+}_{2k}-⟩⟨0{-}_{2k}0|+\text{H.c.}.$$

$$|0-0⟩\leftrightarrow |\downarrow \downarrow \downarrow ⟩,|-+-⟩\leftrightarrow |\downarrow \uparrow \downarrow ⟩,|0--⟩\leftrightarrow |\downarrow \downarrow \uparrow ⟩,$$

$$|--0⟩\leftrightarrow |\uparrow \downarrow \downarrow ⟩,|---⟩\leftrightarrow |\uparrow \downarrow \uparrow ⟩,$$

$${H_3|}_{\text{PXP}}=4\sum _{k=1}^{L/2}|\downarrow \uparrow \downarrow ⟩⟨\downarrow \downarrow \downarrow |+H.c.=4\sum _{k=1}^{L/2}|\downarrow ⟩{⟨\downarrow |}_{2(k-1)}\otimes |\uparrow ⟩{⟨\downarrow |}_{2k}\otimes |\downarrow ⟩{⟨\downarrow |}_{2(k+1)}+H.c.=4H_{\text{PXP}},$$

$$S_0^z(t)=\sum _{\vec{a}}{c}_{\vec{a}}(t){\lambda }^{\vec{a}}.$$

$${\rho }_R(n,t)\equiv \sum _{\text{RHS}(\vec{a})=n}|c_{\vec{a}}(t){|}^2.$$

$$H_4^{\prime }=$$

Full text

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 Gring et al. (2012); Hild et al. (2014); Brown et al. (2015); Kaufman et al. (2016); Tang et al. (2018); Brydges et al. (2018); D’Alessio et al. (2016); Gogolin and Eisert (2016); Meinert et al. (2017). Thermalization is believed to be characterized in terms of the Eigenstate Thermalization Hypothesis (ETH) Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016). 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 Rigol et al. (2008); Kim et al. (2014); D’Alessio et al. (2016).

Given its supposed generality, there has been much interest in systems that violate ETH. Two well-known instances are integrable systems Rigol et al. (2007); Kinoshita et al. (2006) and the many-body localized (MBL) phase Basko et al. (2006a); Nandkishore and Huse (2015); Altman and Vosk (2015); Schreiber et al. (2015), both of which avoid ETH due to the existence of extensively many conserved quantities Essler and Fagotti (2016); Huse et al. (2014); Serbyn et al. (2013). These conservation laws lead to non-ergodicity even at high energy densities. One important question concerns whether behavior similar to MBL can appear in systems without spatial disorder Schiulaz et al. (2015); Yao et al. (2016); Papić et al. (2015); Smith et al. (2017a, b); Brenes et al. (2018); van Nieuwenburg et al. (2018); Schulz et al. (2019).

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 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 Biroli et al. (2010); Mori (2016); Shiraishi and Mori (2017). It is important to stress that if only weak ETH is satisfied, then we can always find initial conditions which have narrow energy distributions but nevertheless fail to thermalize Gogolin and Eisert (2016). Indeed, generic systems are expected to exhibit the strong version of ETH Kim et al. (2014). Recently, however, several exceptions have been discovered Moudgalya et al. (2018a, b); Iadecola and Znidaric (2018); Iadecola et al. (2019); Ok et al. (2019). Systems with constrained dynamics van Horssen et al. (2015); Lan et al. (2018) are an especially promising avenue where non-thermal eigenstates, dubbed scar states, can occur Shiraishi and Mori (2017); Turner et al. (2018a, b); Choi et al. (2018); Lin and Motrunich (2018); Feldmeier et al. (2019). These are believed to be responsible for persistent oscillations observed in a recent Rydberg atom experiment Bernien et al. (2017).

Constrained dynamics occurs naturally in so-called fracton systems, which are characterized by the existence of excitations that exhibit restricted mobility Haah (2011); Vijay et al. (2015). On one hand, these systems have been studied in three-dimensional exactly-solvable lattice models with discrete symmetries, where fractons are created on the corners of a membrane or fractal operator Chamon (2005); Vijay et al. (2015, 2016); Hsieh and Halász (2017); Halász et al. (2017); Shirley et al. (2018). 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. Pretko (2017a, 2018, b); Williamson et al. (2018); Bulmash and Barkeshli (2018); Marcos et al. (2014); Gromov (2018) An analytical connection between the two approaches has been discussed in Refs. Bulmash and Barkeshli, 2018; Ma et al., 2018. The exotic behavior of fractons also gave rise to the study of their non-equilibrium physics Prem et al. (2017); Pretko (2017c), and it has been argued that fracton models with discrete symmetries show glassy dynamics Chamon (2005); Castelnovo and Chamon (2012); Prem et al. (2017).

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 systems, for which a numerical study is feasible. Apart from fracton systems, such charge and dipole conserving Hamiltonians also occur naturally in other contexts, for example in the quantum Hall effect Rezayi and Haldane (1994); Bergholtz and Karlhede (2008); Bergholtz et al. (2011); Wang et al. (2012); Nakamura et al. (2012) and in systems of charged particles exposed to a strong electric field Schulz et al. (2019); van Nieuwenburg et al. (2018). Interestingly, a recent work Pai et al. (2019) 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 discover that the source of non-ergodicity is an extensive fragmentation of the Hilbert space into exponentially many disconnected sectors in the local $z$-basis. In particular, based on the Hilbert space structure, we obtain a lower bound 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 circuits 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. App. D shows an explicit construction, relating the minimal Hamiltonian we consider to the PXP model Lesanovsky and Katsura (2012); Turner et al. (2018a) that appears in the context of quantum many-body scars Zla .

II Model and symmetries

We consider two spin-$1$ Hamiltonians on a chain of length $N$ of the form

[TABLE]

and

[TABLE]

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_{n_{0}}$:

[TABLE]

with respective eigenvalues $q$ and $p$ defining the symmetry sector $\mathcal{H}_{q,p}$ Gra . Since $[Q,P_{n_{0}}]=0$, the local $S^{z}$-basis, denoted by $|+\rangle,|0\rangle,|-\rangle$, is a common eigenbasis of $Q$ and $P_{n_{0}}$. The definition of the dipole symmetry $P_{n_{0}}$ also depends on the reference position $n_{0}$, except when $Q=0$. Unless specified otherwise, we choose open boundary conditions PBC and take $N=2m+1$ odd, labeling sites $n=-m,\dots,0,\dots,m$. We choose the reference site $n_{0}$ to be the center site, $n_{0}=0$, and denote $P\equiv P_{n_{0}=0}$. The operator $P$ does not commute with spatial translations and changes sign under inversion; thus, it is not an internal symmetry Gromov (2018). Dipole conservation is the relevant global symmetry appearing in the description of fracton phases of matter with U(1) symmetry group Pretko (2017a, 2018, b); Williamson et al. (2018); Bulmash and Barkeshli (2018); Marcos et al. (2014); Gromov (2018). Motivated by this, we use the following notations: we call the states $\left|\pm\right>$ on a given site a fracton with charge $q=\pm 1$, and a two-site configuration $\left|+-\right>$ ($\left|-+\right>)$ a dipole with zero charge and dipole moment $p=-1$ ($+1$). Notice that the dipole moment of a $(\pm)$-fracton on a site $n$ is $p=\pm n$. Thus, in order to conserve the total dipole moment, a fracton can only move by emitting dipoles Pretko (2017a); Pai et al. (2019).

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 Rezayi and Haldane (1994); Bergholtz and Karlhede (2008); Bergholtz et al. (2011); Wang et al. (2012); Nakamura et al. (2012). In addition, such dipole-conserving chains can arise naturally in the presence of strong electric fields, as we discuss in the outlook.

III Hamiltonian $H_{3}$

We start by investigating the three-site Hamiltonian $H_{3}$ in Eq. (1), as a minimal model that conserves both the total charge $Q$ and the 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.

III.1 Lack of thermalization

We first investigate the behavior of the auto-correlation function $C^{z}_{0}(t)\equiv\left<S_{0}^{z}(t)S_{0}^{z}(0)\right>$ at infinite temperature. Relying on quantum typicality Reimann (2007); Jochen Gemmer and Mahler (2009); Steinigeweg et al. (2015), we compute $C^{z}_{0}(t)$ for a random state on the full Hilbert space. For thermalizing and translational invariant spin-$1$ systems, $C^{z}_{0}(t)$ is expected to decay to $2/(3N)$ for a chain of length $N$, up to potential boundary contributions Czt . In Figure 1(a) we show $C^{z}_{0}(t)$, obtained via an iterative Krylov space based algorithm Saad (1992), for system sizes $N=13,15$. Instead of relaxing to the thermal expectation value, the auto-correlation saturates to a finite value $C^{z}_{0}(t)-2/(3N)\sim 0.2$ at long times. In App. B we confirm that this finite saturation value persists up to long times, $t\sim 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 $\left<S_{n}^{z}(t)S_{0}^{z}(0)\right>$, which exhibits a peak in the center site at all times. We conclude that the system exhibits non-ergodic behavior. 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

III.2 Hilbert space fragmentation

In this section, we demonstrate that the constrained dynamics of $H_{3}$ leads to a fragmentation of the many-body Hilbert space: most $(q,p)$ symmetry sectors split into many smaller invariant subspaces in the local $S^{z}$-basis, 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.

III.2.1 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 $\left|\bm{0}\right>\equiv\left|\cdots 0000\cdots\right>$, which is annihilated by all terms in $H_{3}$, due to $(S^{\pm})^{2}\left|0\right>=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., $\left|0\cdots 0++0\cdots 0---0\cdots\right>$. These are annihilated by all terms, since $S_{n}^{+}S_{n+1}^{-}\left|\pm\pm\right>=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 follow Pauling Pauling (1931); Onsager and Dupuis (1959) to estimate the 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). Treating the constraints on each triangle as independent, and using that there are $N-2$ triangles and 19 frozen states per triangle, we estimate their total number to be 3^{N}\times\big{(}19/27)^{N-2}\approx 2.02\times 2.11^{N}. 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. The numerical results together with an explicit computation for small system sizes, suggest that this estimate gives a lower bound of the actual number of frozen states. However, we do not have a proof of this general statement. Nevertheless, in Sec. IV.2 we provide an analytical lower bound for arbitrary finite-range dipole-conserving Hamiltonians.

III.2.2 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 ($\left|00\cdots 0\right>$) 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 left- and 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 $\left|++{0+0}++\right>$ and $\left|++{+-+}++\right>$. We can check that these two states can only evolve to each other under $H_{3}$, defining a small invariant subspace. More generally we can consider states of the form $\left|++0\cdots 0{+}0\cdots 0++\right>$: 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 Pretko (2017a); Pai et al. (2019). 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 annihilated by $H_{3}$; the wall therefore remains intact and the fracton bounces back harmlessly. Consequently, if the fractons on both sides of a $++$ block have positive charge, the chain is cut into two disconnected halves, 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 $\left|--0\cdots 0{-+}0\cdots 0++\right>$. 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\pm 1$ Dif . 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, each of which has its own conserved charge and dipole moment, 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.

III.2.3 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.2 we prove that exponentially many invariant subspaces exist for any extension of the model that only involves dipole-conserving interactions with finite range.

We close this section by noting that, apart from the overall fragmented structure of the Hilbert space, which is our main concern in this paper, there is also the possibility of interesting dynamics within certain connected components. For example, as we show in App. D, there are particular subspaces where the Hamiltonian $H_{3}$ maps exactly to the so-called PXP Hamiltonian Zla , studied in the context of quantum many-body scars Turner et al. (2018a, b). A similar mapping has been uncovered in a spin-1/2 version of this model in a recent preprint Moudgalya et al. (2019).

III.3 Saturation value of $C_{0}^{z}(t)$

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, let us define $\mathcal{P}_{i}$ as the projection onto the connected subspaces $\mathcal{H}_{i}$. These projectors form an orthogonal set of conserved quantities ($\mathcal{P}_{i}\mathcal{P}_{j}=\delta_{ij}\mathcal{P}_{i}$), such that one can use Mazur’s inequality Mazur (1969); Suzuki (1971); Caux and Mossel (2011) to lower bound the infinite time average of the charge auto-correlator as

[TABLE]

where $Z_{i}\equiv\mathcal{P}_{i}S_{0}^{z}\mathcal{P}_{i}=\mathcal{P}_{i}S_{0}^{z}$ is the projection of $S_{0}^{z}$ onto $\mathcal{H}_{i}$, and $D_{i}=\text{tr}\left(\mathcal{P}_{i}\right)$ is the dimension of the subspace.

The bound $C_{0}^{z}(\infty)$ is shown in Figure 1(a) for $N=15$ by the dashed horizontal line; we observe that it is close to being tight, indicating that the main cause of the lack of ergodicity is indeed the fragmentation of the Hilbert space. We computed the estimated value for $C_{0}^{z}(\infty)-2/(3N)$ for a variety of different system sizes, and found that the result appears to remain finite in the thermodynamic limit, even increasing slightly with $N$ for the system sizes available in our numerics (blue dots in Fig. 4).

Since the $\mathcal{H}_{i}$’s are invariant and disjoint subspaces, the weight of the operator $S_{0}^{z}$ within a given sector, \mathrm{tr}\big{(}Z_{i}^{2}\big{)}, remains constant under time evolution. Therefore, we introduce the operator weight W_{D}\equiv\sum_{D_{i}=D}\mathrm{tr}\big{(}Z_{i}^{2}\big{)}/\mathrm{tr}\big{[}\big{(}S_{0}^{z}\big{)}^{2}\big{]} as a function of the sector size $D$ for all invariant subspaces $\mathcal{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 $\sim 2.2^{N}$, the size of the largest sector in the entire Hilbert space scales as $\sim 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_{0}^{z}(t)$, even in the thermodynamic limit. We also confirm the same behavior when considering only the largest symmetry sector $(q,p)=(0,0)$ (see Appendices A and B); this emphasizes the relevance of the fragmentation within each $(q,p)$-sector.

IV Combined Hamiltonian $H_{3}+H_{4}$

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$.

IV.1 Thermalization for $H_{3}+H_{4}$

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_{0}^{z}(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, $\left<S_{n}^{z}(t)S_{0}^{z}(0)\right>$, becoming approximately homogeneous at long times, as shown in Fig. 1(c). The remaining small peak 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.2 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}$. 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.3. This effect is even more apparent in the distribution of the operator weight $W_{D}$ (defined in Sec. III.3) for the operator $S_{0}^{z}$, 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. (4): plugging in the invariant subspaces of $H_{3}+H_{4}$ we find that $C_{0}^{z}(\infty)$ approaches the thermal value, $2/(3N)$, exponentially in the thermodynamic limit, as shown in Fig. 4.

From these results we infer that including longer-range interactions makes the system sufficiently ergodic to thermalize. One possible reason for this qualitative difference is that the 4-site terms break the rule stated at the beginning of Sec. III.2.2, thus allowing for the desctruction of blocks of charges that would be inert under the dynamics of $H_{3}$. A different path to break the non-ergodicity 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 lower bound $C_{0}^{z}(\infty)$ for the charge autocorrelator using Eq. (4) 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 non-trivial model is $H_{4}$, which appears to be non-ergodic, while adding 5-site interactions restores ergodicity.

IV.2 Constructing 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 Deutsch (1991); Srednicki (1994); Rigol et al. (2008); Kim et al. (2014). 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 $\ell$, there exist at least $2\cdot 5^{N/\ell}$ frozen states. While for $\ell=3$ this lower bound is not as tight as the Pauling estimate discussed in Sec.III.2.1, it provides useful insight into generic 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 $\ell-1$ ($+$)-fractons on one side and $\ell-1$ ($-$)-fractons on the other. We now prove that this configuration is an eigenstate of any dipole-conserving term with range at most $\ell$, 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 center 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 [math] or $+$, the state is frozen.

One can 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 take two values (e.g $+$ or [math] between blocks of $+$ charge), 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 $\pm$ 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 $\ell$ sites, consisting of a wall of $\ell-1$ positive/negative charges, followed by a flippable spin. Let $F^{\pm}_{k}$ denote the number of different such configurations to the left of the $k$-th wall (but before the flippable spin), ending in a $(\pm)$-block. Then the considerations outlined above lead to the following recursion formula:

[TABLE]

where we have used that $F^{+}_{1}=F^{-}_{1}=1$. Since each step $k\to k+1$ corresponds to a shift by $\ell$ sites, we conclude that the number of frozen states we constructed scales as $2\cdot 5^{N/\ell}$. This is only a lower bound on the total number of frozen states, which can include other configurations not captured by this construction. In particular one could systematically improve the bound by allowing blocks to be separated by more then one site.

We compare the lower bound $\sim 5^{N/\ell}$ to the numerical results on the number of frozen states for Hamiltonians with interactions of range at most $\ell=3,4,5$ in Fig. 7, where we extract the asymptotic scaling. The comparison to numerical data in Fig. 7 shows that the scaling is relatively close to $\sim 10^{N/\ell}$ Thi : $10.73^{N/3}$ ($\ell=3$), $9.94^{N/4}$ ($\ell=4$) and $10.79^{N/5}$ ($\ell=5$). Thus, the lower bound is not tight but it proves the exponential scaling of frozen states.

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 Dim $(2S+3)^{N/\ell}$. Last, we note that in the limit $\ell\to\infty$ 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.

IV.3 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 dynamically 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’. Here we take the point of view of fixing only local symmetries, as non-local ones usually do not lead to distinct distributions for the diagonal matrix elements Santos and Rigol (2010); Sorg et al. (2014); Mondaini et al. (2016). In our case, this means fixing $Q$ and $P$, but not the additional symmetries that correspond to the invariant subspaces, since we expect these to be non-local ETH .

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.2.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 Biroli et al. (2010); Mori (2016); Shiraishi and Mori (2017) 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_{0}^{z})^{2}$, where [math] 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: $\langle{{(S_{0}^{z})^{2}}}\rangle$ takes similar values for 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 $\langle{{(S_{0}^{z})^{2}}}\rangle$ does not become narrower with increasing $N$, as shown in Fig. 8(a). This is in contradiction with 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 non-thermalizing 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, as 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., $\mathrm{max}[D_{(q,p)}^{i}]/D_{(q,p)}\propto\exp(-\alpha N)$ for some $\alpha>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: $\mathrm{max}[D_{(q,p)}^{i}]\approx D_{(q,p)}$, 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: although the majority of the eigenstates obey ETH, non-thermal eigenstates exist even in the bulk of the spectrum Shiraishi and Mori (2017); Moudgalya et al. (2018a, b); Turner et al. (2018a, b). However, while the number of these ‘scarred’ states is usually $\mathcal{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 $\mathrm{max}[D_{(q,p)}^{i}]/D_{(q,p)}$ either decays (exponentially) to zero or approaches unity. It is an interesting and open question, whether systems with intermediate behavior—with 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. Pai et al., 2019. These define a discrete time evolution, the building blocks of which are unitary gates acting on $\ell$ 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. Pai et al., 2019 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 $\ell=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 $\ell=4$, we find that the system thermalizes, also in complete agreement with our results on the Hamiltonian $H_{3}+H_{4}$.

The two circuit geometries, with gates of size $\ell=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 $\ell$ 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 $\ell=3$ ($\ell=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_{0}^{z}(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. E we also consider the spatial spreading of an initial $S_{n}^{z}$ operator and similarly find that for $\ell=4$ the operator is delocalized at long times. We therefore conclude that the localized behavior observed in Ref. Pai et al., 2019 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.3 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, especially for the $\ell=3$ case where we predict that the majority of eigenstates should be nonthermal.

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 auto-correlation 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, and while it is no longer sufficient to make the infinite-temperature dynamics non-ergodic, it nevertheless results in exponentially many non-thermal eigenstates. 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 $\sim 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 conserved operators look like is an interesting open problem.

Another interesting problem for future work is to investigate the equilibrium properties and dynamics of the system at finite or zero temperature, both within the whole Hilbert space and within individual subspaces, as well as to clarify the requirements for strong and weak fragmentation. It would be also worth investigating the relationship between the phenomenon of Hilbert space fragmentation and the concept of reducibility in classical constrained models Ritort and Sollich (2003). 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 further the connections to: many-body localization Basko et al. (2006b); Serbyn et al. (2013), quantum scars Moudgalya et al. (2018c); Turner et al. (2018a), gauged formulations of the considered systems Pretko (2018, 2017d) as well as Stark localization Schulz et al. (2019); van Nieuwenburg et al. (2018). The latter also provides a potential experimental platform for realizing dipole-conserving Hamiltonians of the type studied here in the limit of strong electric fields van Nieuwenburg et al. (2019), which perturbatively lead to the kind of terms we have considered above. Finally, it would be interesting to search for connections between the models discussed in this work and type I fracton models Vijay et al. (2016, 2015), 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 Kim and Haah (2016).

VII Acknowledgments

The authors thank Mari Carmen Bañuls, Claudio Castelnovo, Juan Garrahan, Johannes Hauschild, Daniel Hetterich, Vedika Khemani, Sanjay Moudgalya, Rahul M. Nandkishore, Shriya Pai, Michael Pretko, Nic Shannon, John Sous, Dávid Horváth and Giuseppe de Tomasi for discussions. P.S. specially thanks María del Carmen Giménez Aurusa, and acknowledges the Boulder school 2018, which inspired part of this work. We acknowledge support from “la Caixa” Foundation (ID 100010434) fellowship grant for post-graduate studies (P.S.), the Technical University of Munich - Institute for Advanced Study, funded by the German Excellence Initiative and the European Union FP7 under grant agreement 291763 (M.K.), the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2111-390814868 (F.P., M.K.), from the DFG grant No. KN 1254/1-1 (M.K.), the DFG TRR80, Project number 107745057 (F.P., M.K.), the DFG Research Unit FOR 1807 through grant no. PO 1370/2- 1 (F.P.), the DFG Collaborative Research Center SFB 1143 (R.V.), the Nanosystems Initiative Munich (NIM) by the German Excellence Initiative (F.P.), and the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program grant agreement No. 771537 (F.P.) and No. 851161 (M.K.). This research was conducted in part at the KITP, which is supported by NSF Grant No. NSF PHY-1748958 and the Heising-Simons Foundation.

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 $H_{3}$

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 $\Pi^{z}_{\text{odd}}$. However, one can easily prove that this quantity is fixed by the dipole moment as

[TABLE]

where in the second step we have used that for spin-1, a $2\pi$ rotation is equal to the identity. From this it is clear that the total parity \Pi^{z}=\exp\big{(}{i\pi\sum_{n}{S_{n}^{z}}}\big{)} is obtained as $\Pi^{z}=\Pi_{\text{odd}}^{z}\Pi^{z}_{\text{even}}$ and is related to the total charge as \Pi^{z}=\exp\big{(}i\pi Q\big{)}. In general, the terms in $H_{3}$ are also invariant under the parity transformations given by \Pi^{x}=\exp\big{(}{i\pi\sum_{n}{S_{n}^{x}}}\big{)} and \Pi^{y}=\exp\big{(}{i\pi\sum_{n}{S_{n}^{y}}}\big{)}, which map

[TABLE]

for all $n_{1},n_{2},n_{3},n_{4}$. Note that $\Pi^{x}$ and $\Pi^{y}$ do not commute with $Q$ or $P$.

Moreover, as stated in the main text, there exists an operator \mathcal{C}=\prod_{n}{e^{i\pi\big{(}S_{4n}^{z}+S_{4n+1}^{z}\big{)}}} which anti-commutes with $H_{3}$. Since $\mathcal{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, $\mathcal{C}$ does not anti-commute anymore with the resulting Hamiltonian. There is also at least one additional anti-commuting operator \tilde{\mathcal{C}}=\prod_{n}{e^{i\pi\big{(}S_{4n+2}^{z}+S_{4n+3}^{z}\big{)}}} but since $\mathcal{C}\tilde{\mathcal{C}}=\Pi^{z}$, they are not independent. Note that since $\mathcal{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 $\mathcal{C}_{4}=\prod_{n}e^{\pi S_{4n}^{z}}$, which does not anti-commute with $H_{3}$.

In Fig. 11(a) we show the density of states, $\rho(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\sum_{n}(S_{n}^{z})^{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}\equiv\sum_{D_{i}=D}\mathrm{tr}\big{(}Z_{i}^{2}\big{)}/\Big{(}\sum_{D}\sum_{D_{i}=D}\mathrm{tr}\big{(}Z_{i}^{2}\big{)}\Big{)} 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.

Appendix B Finite-size scaling of the autocorrelator

In this section we present in more detail the finite size scaling of the auto-correlation function and its lower bound.

First, we discuss the scaling of the auto-correlation function $\left<S_{0}^{z}(t)S_{0}^{z}\right>$ 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 $\left<S_{0}^{z}(t)S_{0}^{z}\right>_{(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, $\sum_{n}\left<S_{n}^{z}(t)S_{0}^{z}(0)\right>=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 lower bound $C_{0}^{z}(\infty)$ in Eq. (4) with system size $N$ restricted to the $(0,0)$ symmetry sector of $H_{3}$. In this case the lower bound takes the form

[TABLE]

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 $\sqrt{1/N}\cdot 2^{N}$; hence, the exponents are the same up to logarithmic corrections.

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 Page (1993) $S_{\text{Page}}=\frac{N-1}{2}\log{3}-\frac{1}{6}$.

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 than 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_{\text{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 Mapping to the PXP model

In this appendix we explain the relation between the dipole-conserving Hamiltonian $H_{3}$ introduced in Eq. (1) and the PXP model Zla , that appears in the context of quantum many-body scars Turner et al. (2018a, b); Choi et al. (2018); Shiraishi (2019). Such relation has been already obtained in Ref. Moudgalya et al. (2019) in the spin-$1/2$ version of the models we study (see e.g., Fig. 4) in the context of fractional quantum Hall. The PXP model describes a chain of interacting Rydberg atoms Bernien et al. (2017), which in the limit of strong nearest-neighbor interactions is effectively described by the spin-$1/2$ Hamiltonian

[TABLE]

where the projectors $P_{n}\equiv(1-\sigma_{n}^{z})/2$, ensure that not two adjacent Rydberg atoms become simultaneously excited into the $\left|\uparrow\right>$ state, a phenomenon known as Rydberg blockade. Restricted to the lowest energy subspace with no adjacent excited states, the dimension of the constrained Hilbert space can be shown to be Turner et al. (2018b) $d_{L}=F_{L+2}$ for open (OBC) —in the presence of additional boundary terms OBC — and $d_{L}=F_{L-1}+F_{L+1}$ for periodic boundary conditions (PBC), where $F_{n}$ is the $n$th Fibonacci number. Note that in particular, this subspace contains the Néel states $\left|\mathbb{Z}_{2}\right>=\left|\downarrow\uparrow\downarrow\uparrow\cdots\right>$ and $\left|\mathbb{Z}^{\prime}_{2}\right>=\left|\uparrow\downarrow\uparrow\downarrow\cdots\right>$, whose atypical real-time dynamics has been experimentally realized Bernien et al. (2017) and has been identified as a probe of the existence of quantum many-body scars Turner et al. (2018b).

In the following we show that the dynamics of certain connected subspaces $\mathcal{H}_{i}$ discussed in the main text, are governed by the PXP Hamiltonian in Eq. (9) and identify the analogs to the Néel states in the fractonic language. Let us consider states of the form

[TABLE]

with a $\left|+\right>$ state on every fourth site separated by three $\left|-\right>$’s. In the following we fix the length of the chain to be a multiple of four, such that we contain an integer number of unit cells. For OBC the dipole moment is given by $p(n_{+})=L/2(1-L/2+n_{+})$, where $n_{+}$ is the location of the first $\left|+\right>$ state starting from the left boundary. Thus, the location $n_{+}$ labels different symmetry sectors containing the same spin pattern. However, due to the periodicity of the configuration, there only exist four different dipole moments $p(n_{+})$ containing such configuration. When considering PBC, the dipole is defined modulo $L$.

Recalling that every local term $h_{n}$ in the Hamiltonian $H_{3}=\sum_{n}h_{n}$ takes the form

[TABLE]

it is clear that the only non-trivial local actions of $H_{3}$ on the state $\left|\mathbb{F}_{2}\right>$, are those contained within the blocks shown in Eq. (10). After applying two local terms $h_{n}$ centered around the location of $\left|+\right>$ states at sites $n=2(k\pm 1)$, $\left|\mathbb{F}_{2}\right>$ becomes

[TABLE]

Now the action on the intermediate site $2k$ becomes non-trivial

[TABLE]

One can then realize that only terms $h_{2k}$ centered around even sites generate non-trivial dynamics conditioned to the states on odd near sites, such that the only allowed local transition is $\left|-+-\right>\leftrightarrow\left|0-0\right>$. Then, the restriction of the Hamiltonian $H_{3}$ to subspaces containing configurations of the form given by Eq. (10) becomes

[TABLE]

Note also that there are never two $\left|+\right>$ states in adjacent even sites, i.e., the local configuration $\left|+\right>_{2k}\left|+\right>_{2(k+1)}$ is not generated under the evolution of $H_{3}$. This effectively implements the Rydberg blockade as imposed by the projectors in Eq. (9). With these observations in mind, we can construct a reversible map relating local spin-$1$ configurations centered around even sites $\{2k-1,2k,2k+1\}$, to spin-$1/2$ configurations on even sites $\{2(k-1),2k,2(k+1)\}$ in the PXP model via

[TABLE]

such that Eq. (14) becomes

[TABLE]

i.e., the restriction of $H_{3}$ into this family of connected subspaces becomes equivalent to a PXP model on a chain of length $L/2$ up to a factor of 4. Thus, there exist eight different symmetry sectors (the other four subspaces are obtained applying the $\Pi^{x}$ parity symmetry (see App. A) to the configuration $\left|\mathbb{F}_{2}\right>$), whose evolution is governed by the PXP Hamiltonian. This explicitly shows that quantum many-body scars appear in the dipole conserving Hamiltonian $H_{3}$, similarly to Ref. Moudgalya et al. (2019).

Appendix E Operator spreading of $S_{0}^{z}(t)$

Here we consider another measure of localization, that contains complementary information about the Heisenberg picture evolution of the charge density operator $S_{0}^{z}(t)$ compared to its auto-correlation function. In particular we look at how $S_{0}^{z}$ 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 $\lambda^{a}$ for $a=0,\ldots,8$, where $\lambda^{0}\equiv 1\!\!1$. A basis of operators on the entire chain is then given by products of such local basis elements of the form $\lambda^{\vec{a}}\equiv\bigotimes_{n=-N/2}^{N/2}\lambda_{n}^{a_{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 $\langle A,B\rangle\equiv\text{tr}(A^{\dagger}B)/3^{L}$ where $A$ and $B$ are two arbitrary operators.

Given such a basis, one can always expand the time evolved operator as

[TABLE]

The coefficients $c_{\vec{a}}(t)$ characterize how $S_{0}^{z}(t)$ spreads out in the space of all possible operators. In particular, focusing on spatial spreading, it is useful to classify the basis strings $\lambda^{\vec{a}}$ according to their right endpoints (assuming open boundary conditions), i.e., the rightmost site $n$ such that $\lambda^{a_{n}}\neq 1\!\!1$ but $\lambda^{a_{m>n}}=1\!\!1$. Denoting this site by $\text{RHS}(\vec{a})$ we can define the right endpoint density of $S_{0}^{z}$ at time $t$ as Roberts et al. (2015); von Keyserlingk et al. (2018); Nahum et al. (2018)

[TABLE]

At time $t=0$ this is a delta function at the initial position of the operator, $\rho_{R}(n,0)=\delta_{n0}$. During time evolution, as the support of $S_{0}^{z}(t)$ increases, $\rho_{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 von Keyserlingk et al. (2018), and as a power law when the operator is a conserved density Rakovszky et al. (2018); Khemani et al. (2018). 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 $\rho_{R}(t)$ in random circuits, first with $3$- and then with $4$-site gates. In order to evaluate $\rho_{R}(n,t)$ we represent $S_{0}^{z}(t)$ as a matrix product operator Verstraete et al. (2008) (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 Vidal (2003).

Our numerics only allow us to access small systems of size $N=6,8,10$. To compute the spreading of $\rho_{R}(n,t)$, we place an operator $S^{z}$ on the third site from the left end of the system and calculate $\rho_{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 $\rho_{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

[TABLE]

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}^{\prime}$ does not exhibit a presistent peak in $\rho_{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