Quantum East model: localization, non-thermal eigenstates and slow dynamics

We study in detail the properties of the quantum East model, an interacting quantum spin chain inspired by simple kinetically constrained models of classical glasses. Through a combination of analytics, exact diagonalization and tensor network methods we show the existence of a fast-to-slow transition throughout the spectrum that follows from a localization transition in the ground state. On the slow side, we explicitly construct a large (exponential in size) number of non-thermal states which become exact finite-energy-density eigenstates in the large size limit, and -- through a"super-spin"generalization -- a further large class of area-law states guaranteed to display very slow relaxation. Under slow conditions many eigenstates have large overlap with product states and can be approximated well by matrix product states at arbitrary energy densities. We discuss implications of our results for slow thermalization and non-ergodicity more generally in systems with constraints.


I. INTRODUCTION
The dynamics and thermalization of interacting quantum systems are extremely challenging problems that attract considerable attention due to their fundamental and practical relevance to many areas of physical sciences, including condensed matter, quantum information, statistical mechanics and beyond (see [1,2] for reviews). Despite many advances in the last couple of decades [1,2], even for models as simple as one-dimensional spin chains with local interactions it has not been possible to reach a fully satisfactory and general understanding of these problems, neither through the use of analytical tools nor through numerical methods. The main obstacles [3][4][5][6][7] relate to the growth of quantum correlations, the spreading of information, and the highly entangled nature of the excited eigenstates that dominate the dynamical evolution.
Here we address several of these questions by studying in detail the properties of the quantum East model, introduced in Ref. [42] as a candidate disorder-free system displaying breakdown of ergodicity at long times. This model is inspired by the classical stochastic East model [61], a prototypical kinetically constrained model (KCM) of classical glasses (For reviews on classical KCMs and their application to the glass transition problem see [62,63]). The numerical simulations of Ref. [42] suggested a possible transition in the quantum East model from a thermalizing phase where relaxation is fast, to a phase of slow relaxation where dynamics retains memory of initial conditions for long times indicating the possible absence of ergodicity. However, as it is often the case with numerics for the small systems accessible to exact diagonalization, it is difficult to make convincing extrapolations from the results of [42] for the asymptotic behavior for large system sizes in the quantum East model.
In this paper we provide what we believe is robust evidence of non-thermal behavior in the quantum East model that goes beyond that of other quantum constrained systems such as the PXP model [54] or quantum dimers [45,46]. For technical convenience we consider the case of open boundary conditions. We employ a combination of analytical arguments, exact diagonalization and tensor network methods to show that the model displays a fast-to-slow phase transition throughout its spectrum, by which we mean a change from a dynamical phase where thermalization is fast to a phase where dynamics is slow and even non-ergodic depending on initial conditions. The transition occurs when changing the parameter that controls the balance between kinetic and potential energies in the Hamiltonian across a "Rokhsar-Kivelson" (RK) point [64,65]. In particular, we demonstrate that the slow dynamical phase is characterized by the following: (i) the ground state is ex-ponentially localized at one of the boundaries; (ii) there is number -scaling exponentially with system size -of eigenstates at finite-energy density that are non-thermal, which we show how to construct analytically exploiting the localization of the ground state and the kinetically constrained nature of the Hamiltonian; (iii) of these, at least a number which is linear in size has area-law entanglement, while for the rest their bipartite entanglement is spatially heterogeneous; (iv) these non-thermal eigenstates have large overlap with product states and can be approximated well by matrix product states (MPS) at arbitrary energy densities; (v) it is possible to construct an even larger number of area-law states, each with an extensive number of localized super-spins, that are guaranteed to display very slow growth of entanglement and long term memory of initial conditions. The remarkable range of non-thermal behavior that we uncover in the quantum East model underlines the potential richness of non-equilibrium behavior of quantum KCMs with appropriately tailored constraints.
The paper is organized as follows. In Sec. II we introduce the quantum East model and describe its basic properties. Section III considers the localization transition in the ground state of the model as one crosses the RK point. In Sec. IV we propose an analytic ansatz in the localized/slow phase which allows to construct an approximation to the ground state of a larger system starting from the exact ground state of a smaller system, and which becomes exact in the large size limit. As a generalization of this procedure, we show how to analytically construct approximate non-thermal eigenstates of finite energy-density using as ingredients eigenstates of smaller systems. These become exact eigenstates in the large size limit, and have area-law entanglement. Furthermore, we confirm the existence of low entanglement exited eigenstates under slow conditions by showing they can be efficiently approximated by MPS for larger system sizes.
In Sec. V we construct a large class of area-law states with small energy variance in terms of localized "superspins". While these are not strict eigenstates, unitary dynamics starting from these states is very slow, and we provide bounds to the decay of time-correlation functions and the growth of entropy with time. In Sec. VI, we analyze in detail the statistical properties of the spectrum of small systems accessible to exact diagonalization, showing that the fast/slow transition is manifested in a change of eigenstate characteristics -including their entanglement, localization and closeness to product states -throughout the spectrum. Finally, we conclude in Sec. VII with a discussion about the implications of our results to quantum constrained dynamics more broadly and to the general problem of quantum non-ergodicity in the absence of disorder. Figure 1. Fast vs. slow dynamics in the quantum East model. Relaxation to equilibrium of (time-averaged) density autocorrelator (2) starting from all possible product initial states for N = 10. For s < 0 (left) equilibration is fast and memory of the initial conditions is rapidly lost. For s > 0 (right) relaxation is slow and memory of initial conditions is preserved throughout the simulation window.

II. QUANTUM EAST MODEL
The quantum East model was originally introduced in [42] in order to consider slow quantum dynamics akin to (quasi-)MBL in the absence of disorder with kinetic constraints as the mechanism behind slow evolution. The model is defined in terms of spin-1/2 degrees of freedom on a one dimensional lattice with Hamiltonian [42], where the operator n i = (1 −σ z i )/2 is a projector onto the state |1 in the local z-basis, and σ α i is the Pauli-α operator at site i. When s = 0, the operator in Eq. (1) is (up to a sign) the same as the continuous-time Markov generator of the classical East model, a stochastic KCM much studied in the context of the classical glass transition [61,[66][67][68][69]. For s = 0, it corresponds to the "tilted" generator studied in the context of the dynamical large deviations of the stochastic model, see e.g. [63,70,71]. When considered as a quantum Hamiltonian, s = 0 is a so-called RK point [64,65], where the ground state corresponds to the equilibrium probability of the stochastic model.
When interpreted as a stochastic generator, the operator Eq. (1) corresponds to the "infinite temperature" classical East model. Note that this terminology does not refer to the temperature of the quantum system, but to the characteristics of the equilibrium probability, i.e., the ground state of Eq. (1) at the stochastic point s = 0, which is given by the equal superposition of all configurations. Such as maximum entropy probability is that of an "infinite temperature" equilibrium state. For the "finite temperature" version of the East model, i.e. one where the equilibrium state is not the equal weight combination of all configurations, see e.g. [71].
The factor n i at the front of each term in the Hamiltonian (1) is the kinetic constraint. It represents an operator valued rate, which in the case of H above makes the action of the local Hamiltonian at site i + 1 non-trivial only when n i projects into the state |1 . In the KCM jargon, when this constraint is satisfied, the site i is said to "facilitate" dynamics of its i + 1 neighbour (i.e., the one to the East, thus the name of the model) [62,63]. In contrast to Ref. [42], here we will study the properties of the Hamiltonian (1) with open boundary conditions. We do this for technical convenience, as we do not expect the physics we uncover below to be very different for the case with periodic boundaries.
The key numerical observation in Ref. [42] was the change in the dynamical behavior when the parameter s is changed from one side of the RK point, that is from s < 0, to the other side, that is s > 0. In Fig. 1 we reproduce this observation for the case of open boundaries: we show the relaxation to equilibrium of the normalized two-time density autocorrelator 2O(t) − 1, defined as the time average where n i (t) is the occupation operator in the Heisenberg picture under unitary evolution generated by the Hamiltonian Eq. (1), and Z ≡ i n i (0) is a normalization factor for the initial occupation. The figure shows results for initial states which are product states in the occupation basis (i.e., local z-basis) at different initial fillings (note that magnetization is not conserved in this model). Notice that, for finite systems, the energy is determined not only by the initial polarization Z, but also by the occupation of the last site. This is the reason why we observe two different thermal values for the same polarization. This effect vanishes in the thermodynamic limit. We observe two fundamentally different behaviors of the autocorrelator depending on the sign of s. For s < 0 dynamics is fast and most of the information about the initial state is quickly erased, as expected from thermalization and compliance with ETH [2]. In contrast, for s > 0 dynamics is slow and for a large class of initial product states, memory of the initial conditions is retained at arbitrarily long times. This is indicative of a transition in the quantum dynamics of the system.
Motivated by these results, in the following we will analyze the structure of the eigenstates of the Hamiltonian in order to collect information about the dynamical properties of the model both for finite system sizes and in the thermodynamic limit.

A. Symmetries of the quantum East model
Since the Hamiltonian is identically zero on the empty string |0 . . . 0 , for open boundary conditions the Hilbert space splits in blocks that are not connected by the dynamics. Each block is determined by the position of the first occupied site, i.e. the k-th block corresponds to the subspace spanned by all classical configurations that start with a string of k − 1 zeroes followed by a 1.
In the following, we will mostly focus on the dynamics of a single block, with N (dynamical) sites to the right of the first occupied one. The position of the latter naturally introduces an edge, and the effective Hamiltonian on the N dynamical sites to its right reads Since [H N , σ x N ] = 0, the Hamiltonian in Eq. (3) can be further divided in the sum of two commuting terms In the rest of the paper we will study and discuss the properties of the Hamiltonians in Eq. (1), (3) and (4).
B. The special case s = 0 At the RK point, s = 0, the Hamiltonian (1) has an additional symmetry. It can be written as a sum of projectors H = i n i ⊗Π − i+1 which, in addition to the empty string, annihilates also a string of |+ states. Thus the Hilbert space splits further in blocks determined by the lengths m and n of, respectively, the leading empty string and the trailing string of |+ . Hence the eigenstates have the form |0 ⊗m |1 |ψ N B |− |+ ⊗n , where N B is the length of the dynamical part of the (m, n)-block, and |ψ N B is an eigenstate of the corresponding effective Hamiltonian,

III. GROUND STATE LOCALIZATION PHASE TRANSITION
We now show that the ground state of the quantum East model (4) is localized when s > 0. Namely, in the ground state of a block of (dynamical) size N , the probability of finding an occupied site decreases exponentially with the distance to the left edge introduced by the first occupied site. The localization length ξ can be extracted already at small sizes, accessible by exact diagonalization, by analyzing the expectation value in the ground state of the local operator n i as a function of the position i. This is shown in Fig. 2 for the ground state of H N + , with N = 15. For s < 0 we observe an almost homogeneous occupation, independent of the system size and the value of s. For s > 0, in contrast, the occupation decays fast with the distance to the edge, with faster decay as we increase s. We find that the results can be fitted assuming an exponential decay, and the localization length ξ from the fit captures the phase transition at s = 0. Indeed, we find that the value of ξ diverges as s = 0 is approached from the positive side, according to ξ ∝ s −ν with ν = 0.533 ± 0.006 (see inset of Fig. 2). These results hold for the ground state of H N + in Eq. 4. We observe the same qualitative behavior for the ground state of H N − . Indeed, both Hamiltonians differ only in the last site, H N + − H N − = −e −s n N , with the difference decreasing fast for s > 0.
The form in Eq. (6) provides a reasonable fit of the numerical data for the occupation, but a more detailed look at our numerical results suggests in fact a fasterthan-exponential decay, as shown in Fig. 3. Indeed, in appendix A we show that for large s perturbation theory provides an approximate decay of the form As can be seen in figure 3, this is in good agreement with the numerical data.
In Fig. 4 we provide a cartoon picture of the ground state which for s > 0 is localized near the edge. The spatial structure of the GS revealed by these studies can be understood in light of the adiabatic theorem. Away from the phase transition, which happens at s = 0 [72], the system is gapped, and we can apply the adiabatic theorem to connect the ground state to the non-interacting one at s → ∞. The latter corresponds to the product state with only the first site occupied, |10 . . . 0 . Within the gapped region, the evolution with the adiabatically changing Hamiltonian will dress the initial site with an exponential tail like the one shown in our numerical results and depicted in Fig. 4 (see also appendix A).

IV. EIGENSTATES FOR LARGE SYSTEM SIZES
Given the localization properties of the ground state discussed above, and the peculiar form of the Hamiltonian, in this section we provide an ansatz for the ground state and some excited states of finite energy density at arbitrarily large system sizes.
A. The ground state for large system sizes Consider the normalized state where |GS N −1 + is the ground state of the Hamiltonian . We want to show that, in the localized phase, |Ψ 0 (L; N ) is close to the ground state of H L in Eq. (3), supported on L sites. In appendix B, we demonstrate that the only contribution to the energy variance comes from the boundary term between |GS N −1 + and the string of empty sites. By using H N −1 , it can be easily seen that neither the mean value, nor the variance of the energy evaluated in |Ψ 0 (L; N ) depend on L, and they take the simple form where we have defined δ = GS N −1 (9) show that both the mean energy and the variance of the state |Ψ 0 (L; N ) (supported on L sites) can be estimated from the knowledge of |GS N −1 + (supported on N − 1 < L sites). For small values of δ, namely when the last spin of |GS N −1 + is close to |0 , the state |Ψ 0 (L; N ) is close to an eigenstate of H L for any L. As can be seen in Fig. 2, this is precisely the case when s > 0. Eq. (9) also shows that the quantity δ fully quantifies the energy variance of the extended state. Accordingly, as long as the variance is smaller that the gap (which is sizable already for small positive values of s and for all system sizes), we expect that the state |Ψ 0 (L; N ) approximates the ground state of the Hamiltonian, independently of L.
Notice that the form of Eqs. (8), (9) are also valid (with the E 0 + in (8) replaced by the appropriate energy) if the factor |GS N −1 + in Eq. (7) is replaced by any other eigenstate of H N −1 + . In Sec. VI we will use δ as a figure of merit for quantifying the number of eigenstates that admit an extension as the one in Eq. (7), with small variance. We will show that, for positive values of s, the property above is shared by several eigenstates of the model and not only by the ground state.

B. Excited states for large system sizes
As we have shown above, by combining the ground state of small systems and strings of empty sites, it is possible to approximate ground states for large system sizes. The construction utilizes two particular ingredients: the localization properties of the ground state, and the fact that the Hamiltonian annihilates a string of empty sites. In this section we will construct an ansatz for excited states based on similar ideas. Suppose |φ M is an ex- (such that L = N + M + 1) exhibits similar properties as the one defined in Eq. (7). More precisely, as in the previous case, the only contribution to the energy variance comes from the boundary term between the ground states and the empty site and is given by Eq. (9). The corresponding expectation value of the energy is E = E 0 + + E M + e −s δ/2. Notice that the states in Eq. (10) can be arbitrarily close to an eigenstate of H L in Eq. (3) as long as δ is small enough. Since the typical energy gap between two neighboring eigenstates in the middle of the spectrum for a generic Hamiltonian supported on L sites scales as 2 −L , in order to provide accurate approximations, δ needs to decrease at least as fast. As illustrated by Fig. 3, δ decays super-exponentially, δ ∼ exp(−N log N ), which implies that N log N > ∼ L will be enough to satisfy that condition. For very large system sizes (L → ∞) this can be achieved if the ground state occupies a fraction of the sites N/L approaching zero. Therefore, the fraction M/L of sites that can be occupied by an excited state approaches one as we increase the system size. As M becomes larger, the states |φ M can reach higher energies leading to any finite energy density for the states |Ψ (L; N ) .
It is worth stressing that the approximate eigenstates |Ψ (L; N ) are non-thermal and, as long as M = O(L), they are exponentially many in system size L. More precisely, for any given N , there are 2 L−(N +1) states of that form: a fraction 2 −(N +1) of the total number of states in the Hilbert space.
Exploiting the maximally excited state The construction we just described provides an explicit way of addressing excited states at large system sizes by using eigenstates from smaller sizes. In general, nevertheless, states of the form (10) do not need to fulfill an area law of entanglement, even if the leftmost N sites are always in a product state with respect to the rightmost M +1 sites of the system, because a highly excited eigenstate |φ M may have volume law entanglement. Thus, the description of |Ψ (L; N ) may require exponential resources. However, there is at least one interesting exception to this situation, when the excited state |φ M corresponds to the maximally excited state of the Hamiltonian H M in Eq. (3), or equivalently, the ground state of −H M which also admits a MPS approximation.
If we choose |φ M in Eq. (10) to be the maximally excited state |φ M max , we obtain an area-law state , as long as M = O(L), the resulting |Ψ max (L; N ) has finite energy density. Moreover, its energy variance is ∆H L < δ, so that in the localized phase it can be made arbitrarily small by increasing N , and the construction can provide approximate eigenstates.
From the exact diagonalisation results above we know that even for small system sizes δ quickly reaches machine precision at least exponentially fast in N . This means that even for modest N its value becomes negligible in the construction above. The construction above immediately suggests an efficient numerical algorithm to construct quasi-exact highly excited eigenstates for system sizes much larger than the ones allowed by exact diagonalization, since we can use variational MPS methods to find the ground states of H N and −H M for chains of several hundred sites with extremely good precision [71]. Fig. 5 illustrates the construction for a chain of size L = 30. In particular, we show the energy variance and occupation distribution of MPS approximations to excited states, found numerically as described in Sec. VI D. For s > 0 and small energy densities, for which the MPS provide almost exact eigenstates, we observe that their spatial profile indeed agrees with that of the analytical construction presented in this section. Moreover, for s > 0 the construction yields energy variances close to machine precision over practically the whole range of energies, where the direct MPS search is far from reaching an exact eigenstate.

V. THE SUPER-SPIN PICTURE
Here we exploit the results from previous sections to engineer a large class of states with small variance. The basic idea is concatenating several blocks of N + 2 sites, each of them in one of two mutually orthogonal states, We identify the subspace spanned by these two vectors with the Hilbert space of a super-spin. key dynamical property that we exploit is their energy variance, which can be easily computed using the same procedure as in section IV A. Since |0 blocks do not contribute to the variance, the only contributions come from blocks in |1 , and correspond to the value computed in Eq. (9) where the index k runs over the positions of the occupied super-spins and M is the Hamming weight of |s . It is important to stress that M is potentially unbounded in the thermodynamic limit, in which case the variance becomes unavoidably large. The energy can also be easily computed, Equations (12) and (13) show that if M is chosen appropriately we can construct states with high energy and exponentially small variance in N . Notice however that, as we want states with small energy variance, we need to introduce limitations on the values of M. From Eq. (12), note that the variance of the super spins cannot exceed the value δL/(N + 2) < ∼ 2 −N log(N ) L/(N + 2) < 2 −N L/(N + 2), since for any given super spin we can accommodate at most M = L/(N +2) occupied sites. Clearly, we have the freedom of choosing N at will. However, the choice will affect the variance of the super spins, and the dimension of the Hilbert space spanned by them. It is illustrative to mention a few interesting cases. (i) If N ∼ log L β with β > 1 then, the state with M occupied super spins can have a large variance, exponentially larger than 2 −L : H s < L 1−β / log L β . The dimension of the corresponding Hilbert space is of the order 2 L/ log L β . (ii) An opposite scenario is when N ∼ L/ log L. In this case, the variance is small H s < 2 −L/ log L log L, but the dimension of the Hilbert space is linear in L. (iii) An interesting intermediate example consists in N ∼ L α , with 0 < α < 1. Here the variance is H s < L 1−α 2 −L α and the dimension of the Hilbert space scales as 2 L 1−α , which is sub-exponential in system size. In the following section we will show how the variance of an initial state can be use to quantify its slowness. The superspin picture provides a flexible platform where one can choose the appropriate trade-off between the dimension of the Hilbert space and the dynamical activity of the super-spin vectors that span it.

A. Dynamical properties of the super-spin states
The memory of the initial state during time evolution admits a general bound based on the initial energy variance. For an initial state |ψ 0 , we define the overlap [73] where |ψ t is the state at time t and ρ t = |ψ t ψ t | the corresponding density matrix. Using the Cauchy-Schwarz inequality, where · F denotes the Frobenius norm. In the second line we used the fact that for the commutator with the Hamiltonian this norm does not depend on time. Exploiting Eq. (15) we can compute the memory of the initial state as which leads to the bound where we used a(0) = 1, da dt | t=0 = 0 and a(t) ≤ a(0). Eq. (17) is a general bound on the memory of a time evolved state based on the energy variance of the corresponding initial state.
The bound in Eq. (17) can be used to bound the growth of the entanglement entropy of an arbitrary subsystem. According to the Fannes inequality, for any pair of density matrices M 1 , M 2 , of dimensions D × D [74], where where in the second inequality we used and Eq. (17). By plugging Eq. (19) in Eq. (18), we can bound the growth of the entanglement entropy as The bounds on the memory of the initial state in Eq. (17) and the growth of the entropy in Eq. (21), can be straightforwardly applied to the super-spins |s , defined in Sec. V. In the particular case when ρ 0 = |s s| (supported on L sites) we can bound the memory of the initial conditions by using Eq.(17) and ∆H s < Mδ. Namely, Accordingly, if we take σ 0 to be the corresponding reduced density matrix for a region of N L sites, D = 2 N . The bound in Eq. (21) then reads In the previous sections we showed that in the localized region δ decreases exponentially with N . As a consequence, if M is sufficiently small, Eq. (22) and Eq. (23) provide tight bounds on the dynamics of |s . Specifically, in order to erase half of the memory of the initial state, i.e. | s(t * )|s(0) | 2 ≤ 1/2, the dynamics needs at least exponentially long times in N , t * > ∼ (Mδ/2) −1/2 . At the same time, for entangling a sub-region of size N , i.e. 2t √ MδN ∼ 1, the time evolution necessitates exponential times of the order t * ∼ (2N √ Mδ) −1 . We conclude that the dynamics of the states |s , in order to entangle a sub-region, requires exponentially long time in the subsystem size.
The states |s can then be seen as an orthonormal set of quasi-conserved area-law vectors, and any superposition of them will result in a state whose dynamics at short times is governed by dephasing only. The superspin picture thus provides an effective description of a subset of the Hilbert space in the thermodynamic limit which evolves slowly in time, is weakly entangled, and efficiently simulable.
The results in Ref. [42] can be reinterpreted from a super-spin point of view. It was numerically argued there, for the case of periodic boundary conditions, that for certain product states the dynamics exhibits a slow growth of the entanglement entropy, exponential in system size. The slowness of the state was quantified by the number of empty sites following an occupied one. Since the previous statements about the energy variance of |s do not depend on the boundary conditions and, as argued in section III, the block ground state for s > 0 is very close to the product state |1000 . . . , the bound in Eq. (23) gives a rigorous interpretation of the previous numerical observations.
Extensions The super-spin construction described above can be made more general in several ways. On the one hand, a larger set of states can be constructed by allowing not only the ground state, but also (sufficiently localized) excited states as building blocks |φ M . In Sec. VI we show that such excited states do actually exist. On the other hand, by combining the super-spin picture with the excited state construction in Sec. IV B we can also construct states with finite energy density. Namely, we can construct states |S = |s ⊗ |φ M max with energy H S = M(E 0 + +e −s δ/2)+E M max and energy variance ∆H S < Mδ. By increasing M the energy density can be increased, but at the cost of reducing the dimension of the super-spin subspace to 2 (L−M )/(N +2) .

VI. NON-THERMAL EXCITED STATES IN SMALL SYSTEM SIZES
In the following we explore the properties of the whole Hamiltonian spectrum using exact diagonalization for small system sizes. The results indicate a substantial change in the properties of eigenstates across the spectrum in the region s > 0. In particular, in this region many localized eigenstates can be found, beyond the ground state, which can be used in the constructions of the previous sections.

A. Entropy of the exact eigenstates
Since eigenstates of the Hamiltonian incorporate the whole information about the dynamics of the system, their entanglement entropy is often used as an indicator of the associated dynamical behavior. In Fig. 6, we show the entanglement entropy in the middle of the chain of N = 14 spins from exact diagonalization, for two values of s. For negative s, the entanglement entropy of eigenstates exhibits behavior compatible with a thermalizing system -apart from the extremes of the spectrum, most of the eigenstates have large entanglement, almost saturating the upper bound given by system size. In contrast, for positive values of s a considerable number of excited eigenstates have low entanglement. This is an indication of non-thermal eigenstates and reminiscent of the quantum scars found in the PXP model [54], but here we observe this behavior for a much larger number of states.
In order to collect detailed information about the distribution of the entanglement along the spin chain, we compute, for each eigenstate, the entanglement entropy with respect to all possible cuts of the chain, S i = S (ρ i ), where ρ i is the reduced density matrix obtained when tracing out all but the leftmost i spins. In Fig. 7, we plot the entanglement entropy S i and single site occupation n i as a function of the position of the cut (respectively the site) i for all eigenstates in the case s = 0.5 and N = 15. The figure suggests a peculiar heterogeneous entanglement structure of a significant number of eigenstates, for which both quantities decay exponentially as the cut moves to the right. In other words, for many eigenstates, the spins far from the left edge are almost in a product state with the rest of the system, and the corresponding sites are almost empty. These results are quali-   tatively similar to the ones discussed in Sec. III where we analyzed the localization properties of the ground state.

B. Small-δ eigenstates
We diagonalize the Hamiltonian H N + in Eq. (4) for different system sizes N and values of s. Given the set of eigenvectors, we consider the probability distribution of the last site occupation δ = n N , the parameter which, as discussed in Sec. IV A, quantifies the variance of the extended eigenstates |φ N ⊗|0 . . . 0 . Fig. 8 shows the histogram of the corresponding probability density function PDF(δ). Notice that, for positive values of s, many eigenstates exhibit surprisingly small values of δ. Namely, there are several eigenstates |φ N such that the energy variance of the state |φ N ⊗ |0 . . . 0 can be bounded by extremely small values.
In order to quantify the number of eigenstates with small δ, in Fig. 9 and 10 we consider the cumulative distribution function CDF(δ). In particular, in Fig. 9 we observe an abrupt change from negative values of s, where most of the eigenstates have large values of δ to positive s, where more and more eigenstates have very small values. For the sizes accessible by exact diagonalization, the fraction of eigenstates with small δ does not seem to depend on the size of the system. In Fig. 10 we show the energy-resolved CDF. In particular, we divide the spectrum in four intervals of equal energy width, which we number in order of increasing energy. The figure shows that most of the small-δ eigenstates are concentrated in the lower part of the spectrum, in agreement with the results in Fig. 7. As s increases, we observe that the number of eigenstates with small δ values grows for all energy regions, as we indeed expect from the discussion in the previous sections and the smaller localization length.

C. Geometric entanglement
The geometric entanglement of a state, defined as its minimum distance to a product state, also provides interesting insights about the properties of the eigenstates. Given a pure state, the geometric entanglement can be found by maximizing its overlap with a product state. Although it is possible to solve this optimization problem with exact or approximate numerical algorithms, in our case this is unpractical, since we need to repeat the calculation for each eigenstate. Instead, we apply a simpler one-sweep truncation strategy to construct an approximation to the closest product state. Namely, for each eigenstate, we sequentially perform a singular value decomposition with respect to each cut of the chain and keep only the largest singular value for each of them. The resulting product state, once normalized, provides a lower bound to the maximum overlap.
In Fig. 11 we plot the probability density function (over all energy eigenstates) of this estimate for the maximum overlap.   Figure 10. Energy-resolved cumulative distribution function of δ. We split the eigenstates in four equal-size energy intervals, ordered by growing energy. Most of the eigenstates with small δ are concentrated at low energy. For the largest s = 2, the count for eigenstates with small δ increases at all energies.

CDF (δ)
dicating that many eigenstates have a large overlap with product states.
An alternative view of this feature is demonstrated in Fig. 12  uct states. They are mostly concentrated at small energy densities, but Fig. 12 shows that large overlaps can actually be found at arbitrary energy densities.

D. Numerical approximation of non-thermal excited states for large sizes
The discussion in Sec. IV B indicates the existence of highly excited states with small entanglement, that can be written as MPS. We can thus try to find them with numerical methods. There are several possible variations of DMRG to try and target excited states [76]. The simplest one attempts to find the MPS that minimizes the expectation value of the operator W = (H − λ) 2 , where λ is the target energy of the desired states. Since H is a matrix product operator, also W has that form and the minimization can be run efficiently with standard MPS algorithms. We use this tool to probe the whole energy spectrum for eigenstates that can be approximated by MPS.
In the numerical study we fix the system size to L = 30 and the bond dimension to D = 50. For several values of s, we then collect 300 data points, uniformly distributed in energy (excluding the edges of the spectrum). In Fig. 5 (left) we show the energy variance as a function of the energy density. We observe that, for s > 0 and low energy, the algorithm produces MPS with variance close to machine precision. Fig. 5 (right) shows the profile of the expectation value of the single-site occupation number n i as a function of the site i. At low energy densities, and positive s, the optimization finds states with a structure that resembles |Ψ max , with exponentially decreasing occupation from the left edge to the right until a certain site, where the occupation abruptly increases to stay close to one until the right edge. In Fig. 5 we marked with a black circle the states with energy variance smaller than 10 −8 . We find that all the states with small energy variance have an exponential tail which starts from the left edge. According to our analytical construction in Sec. IV, the position of the jump should correspond to the energy of the state. For large energies such construction becomes harder to capture, and the optimization is forced to search for a trade-off between accurate target energy or small energy variance. Notice that our optimization is not tailored to search for this specific construction, as each run starts from a random initial MPS.

VII. CONCLUSIONS AND OUTLOOK
We provided an in-depth analysis of the dynamical properties of the quantum East model. The model is known [70,71] to have a first-order quantum phase transition at the critical point s = 0. Here we showed that, in correspondence to the phase transition point, the ground state undergoes a localization transition from completely delocalized to super-exponentially localized (s > 0). We showed how this ground state transition leads to a singularity throughout the spectrum from a fast dynamical phase at s < 0, where ergodicity is established quickly under unitary evolution, to a slow dynamical phase at s > 0 where thermalization is impeded. We provided rigorous results about the dynamical consequences of this transition focusing on the behavior in the slow nonthermalizing side.
By combining analytical arguments, exact diagonalisation and tensor network methods, we made three key findings. First, the localized nature of the ground state for s > 0 allows for a systematic construction of the ground state for arbitrary system sizes. This construct is very simple, that of a tensor product of the ground state of a small system with a completely empty state. Since the second factor is annihilated due to the constraints, all cost is concentrated at the juncture, which the localization in the first factor makes vanishingly small in the large size limit.
Secondly, this procedure can be extended to obtain exact large-size eigenstates of finite energy density. By replacing the right factor by an excited state, one can systematically construct a large number of non-thermal excited states. If the right factor is that of the eigenstate of maximal energy, the ensuing large-size eigenstate has area law entanglement. This means that there are (at least) polynomially (in system size) many area law eigenstates of finite energy density. Exact diagonalisation of small systems suggests furthermore that many other eigenstates have small entanglement, so repeating the construction with these as the rightmost factor should allow us to obtain an even larger number of area law eigenstates.
Thirdly, by generalizing the tensor product construction to many junctions we can define an even larger class of non-thermal states in terms of what we call superspins. A state composed of super-spins is the tensor product of several ground states for a finite system of a fixed size, possibly separated by empty blocks of the same size, and thus corresponds to a dressed occupied spin localized at each occupied juncture. From the arguments above, if the number of super-spins scales sub-extensively with system size, and the distance between junctions is large enough, such states become area-law eigenstates in the large size limit. States with extensive number of superspins in contrast, while may still have small energy variance, are not guaranteed to be eigenstates. Nevertheless, it is easy to prove that under unitary evolution they retain memory of initial conditions and take exponentially long times to entangle a small sub-region. These are the states that underpin the slow dynamics of the model.
We supported our claims with extensive numerical results for small systems obtained with exact diagonalization, as well as for large systems using tensor networks. In particular, we considered several quantities of interest, such as the entanglement entropy of the eigenstates, the distributions of their last site occupation and of their maximal overlap with a product state. The statistical analysis shows that many eigenstates exhibit atypical behavior signaling the presence of non-thermal dynamical properties. These properties confirm for small sizes the (apparent) singular change throughout the whole spectrum as one changes the parameter s from negative to Figure 13. Single-site occupation of the quantum North-or-East model ground state with s = 0.01 on a 10 × 10 lattice. The low-entanglement structure allow to perform accurate DMRG calculations [77]. Even for small values of s, the ground state exhibits strong localization.
positive. Although all our analytical constructions rely on the localization of the ground state, we provided numerical evidences that many other eigenstates have similar localization properties. Our results suggest the existence of an emergent non-ergodicity of the quantum East model in its slow dynamical phase, stemming from the localization properties of many of its eigenstates, which is already evident at small system sizes. This provides a robust analytical framework for the investigation of the breakdown of ergodicity in this system.
Our findings suggest broader implications for constrained systems and opens interesting possible research directions: (i) It is natural to compare the large number of nonthermal finite-energy-density eigenstates we found here for the quantum East model to quantum scars [54] in the PXP model. The PXP [78,79] is a constrained system known to thermalize [80] for typical states, but possessing a linear in size number of low entanglement excited eigenstates across the spectrum -the scar states -responsible for non-thermalization if starting from particular initial conditions [54]. The PXP [78,79] is a constrained system that thermalizes [80] for typical states, but contains a number of low-entanglement exited eigenstates. These scar states -linear in system size -are responsible for the lack of thermalization when starting from particular initial conditions [54]. For the quantum East model we find an even larger number of non-thermal eigenstates.
(ii) There is one key difference with the PXP model. The PXP constraint is similar to the classical "2-spin facilitated" Fredrickson-Andersen (FA) model [62], i.e., spin-flips can only occur when two nearest neighbors are in the down state, in contrast to the "1-spin facilitation" (albeit directional) of the East model. The PXP constraint is thus stronger than that of the East model, and as a consequence, in one dimension the state space breaks into exponentially many dynamically disconnected blocks. In contrast, as we explained above, in the East model (in periodic boundary conditions) all states are connected except the fully empty one. Thus a weaker constraint in the East model gives rise to stronger dynamical features, associated with the fact that regions devoid of occupied spins are locally frozen -yet still dynamically connected -a feature we exploited in the construction of non-thermal states.
(iii) The slow dynamics of the East model is a firstorder phenomenon -cf. the transition in the ground state -is a consequence of having spatial coexistence of two very different kinds of dynamics. That is, since a region with no occupied spins is locally stable, it can only be relaxed starting from the interface with a more active region. While here we studied explicitly the spectral properties of the East model with open boundaries, we expect to find similar slow characteristics in the case of periodic boundaries, and for other one-dimensional models with similar constraints, such as the quantum FA model [43] (a 1-spin facilitated model but with a symmetric constraint).
(iv) Finding new and broader classes of area-law eigenstates might give insightful information about the interplay between localization in disorder-free models and ergodicity breaking. Along the lines of our constructions, one possible direction includes the characterization of excited states for smaller sizes in order to promote them to fundamental building blocks for eigenstates at larger sizes in other systems.
(v) Perhaps the most interesting research direction is the extension of the ideas here to higher dimensions. From classical KCMs [63] we know that qualitative features are not very dependent on dimensionality. It would be natural to study two quantum generalizations of KCMs in particular, the quantum North-or-East model (see e.g. Ref. [7]), which extends the East model constraint to two dimensions, and the quantum 2-spin facilitated FA model (see Ref. [62] for the classical version), both of which have the essential feature of local "inactive" regions annihilated by the constraint, which is an essential building block in the constructions above. Furthermore, a DMRG calculation for the quantum Northor-East model [81] on a 10 × 10 lattice shows that the ground state for s = 0.01 is strongly localized, too (see Fig. 13). This suggests that these two-dimensional KCMs should also display slow dynamics and a prominence of non-thermal states like the ones uncovered here for the quantum East model, therefore being candidates for nonergodic and non-disordered quantum systems beyond one dimension.  Figure 14. Transition graph of states in perturbation theory. The states contained in the first seven orders of perturbation theory are illustrated as a graph, where each nodes groups states with the same k (the number of occupied sites) and (the position of the rightmost occupied site). The number of states for fixed (k, ) is indicated on each node. The thickness of each edge represents the number of transitions between states contained in the nodes via an application of RV . We notice that the graph has the structure of a Pascal triangle, where its shallow diagonals are highlighted with a dashed line. Notice that with an even (odd) number of hops along the edges from the topmost node |1000 . . . we end up on an even (odd) shallow diagonal, colored in yellow (red).

Appendix A: Perturbation theory
We split Eq. (3) into a bare Hamiltonian and an interaction term where We treat ε = e −s as a small parameter, so we can look at the first order corrections in ε of the ground state deep in the localized phase, i.e. s → ∞.
The bare Hamiltonian H 0 is diagonal in the computational basis and has as unique ground state |ϕ (0) = |000 . . . with a unit gap above it. The excited states are labelled by |k α , where k = 1, . . . , N counts the number of occupied sites and α labels the degeneracy: H |k α = k |k α .
Using the notation in Ref. [82] we can express the corrections to the ground state as a Taylor expansion in the parameter ε H |ϕ = E 0 |ϕ , |ϕ = |ϕ (0) + ε |ϕ (0) + ε 2 |ϕ (0) . . . , By defining the resolvent matrix the n-th order terms can be computed as From a more physical point of view, the perturbation should be seen as as the introduction of a local impurity that has an effect that is localized at the boundary, rather than a global perturbation giving an extensive energy contribution. This can be seen by looking explicitly the first orders in Eq. (A5): The corresponding energy is straightforward to compute We now comment on the convergence radius of the perturbation expansion. Firstly, we notice that in a Banach space, if the series n=0 |ε n | ϕ (n) converges, then the series for |ϕ converges too. We shall now prove that the norm ϕ (n) 2 grows as the n-th power of the golden ratio ϕ. In order to understand the structure of each perturbation theory order, it is instructive to construct a graph of transition between the states, as shown in Fig. 14. The states are grouped in cluster labelled by the number k of occupied sites and the length of the non-trivial string. Each cluster contains k−1 states, and applying RV to a state |k, generates up to k + 1 states in at most 4 distinct ways: |k ± 1, or |k ± 1, ± 1 . Then the n-th perturbation term can contain at most the first n+1 2 even or odd "shallow" diagonals, depending on whether n is even or odd. In other words, the number of states generated by applying RV to the n-th perturbation order is equal to the number of states in the order (n − 1) plus the number of states in the n-th shallow diagonal, which is the corresponding Fibonacci number F n , since n−1 2 j=0 n−j−1 j = F n . Hence the total number of states # (n) in the n-th order is # (n) = n/2 j=0 F 2j+1 = F n+1 n odd n/2 j=0 F 2j = F n+1 − 1 n even (A8) which asymptotically tends to ϕ n+1 . The amplitude of each state will be the sum of all the paths of length (n − 1) leading to it from |1000 . . . , weighted by k 1/k of all the visited states |k α along the path. The weighted sum cannot exceed (1/2) (n−1)/2 , since this is the path connecting |ϕ (1) and |ϕ (2) . Furthermore, the number of distinct paths cannot exceed 4 n−1 , since in principle we can move in the 4 directions for (n−1) moves, as shown in Fig 14. While these estimates are not particularly tight, they are sufficient to bound the norm of |ϕ (n) as ϕ (n) 2 < 2 2(n−1) 2 (n−1)/2 2 ϕ n+1 ∼ (2 3 ϕ) n (A9) Thus the perturbation series in Eq. (A3) is at least convergent in the radius or equivalently s > ∼ 1.28. While this result is used to prove a positive radius of convergence, in practice we can monitor the rate of convergence of the norm of each order to have confidence in the method. Indeed, numerically, we observe a rapidly converging norm in the regime s > ∼ 0.9.
The perturbation theory picture not only gives us a method to compute the first order corrections to the ground state, but also an understanding of its structure. Indeed, looking at the first terms in the expansion, we notice that the terms remain localized close to the left boundary. As we apply powers of V , higher order terms in Eq. (A5) become progressively delocalized, but these contributions get damped at least by a factor ε n .
Let us define the domain-wall state |Θ m = |1 ⊗m |000 . . . . From Eq. (A6), we notice that Θ n |ϕ (n) = 1/n! and that Θ m>n |ϕ (n) = 0. At each order n there is a contribution |Θ n coming from a |Θ n−1 in the previous order. One can think of this term as the fastest possible "excitation" created by applying V repeatedly. Since all other terms in |ϕ (m) have a smaller support, clearly n m |ϕ (m) = 1/m!, and more generally, at the first order n r |ϕ = ε r r! |Θ r + . . . .
We can then conclude that We have dropped the normalization factor, since ϕ|ϕ = 1 + O ε 2 leads to sub-leading corrections. Checking with finite-order expansions, Eq A12 gives a qualitatively accurate prediction of the decay of the occupation number, as shown in Fig. 15. Using Stirling's formula, the asymptotic behavior of Eq. (A12) at large distances is n r ∼ exp[−2r log r + 2r(1 − s)]. Therefore, in the perturbative regime, the ground state exhibits superexponential localization around the first site.

Appendix B: Computation of the variance
Consider the Hamiltonian defined in Eq. (1) on a one dimensional lattice of N sites, and the eigenstates, We consider the normalized state in the |+ sector, similar considerations hold in the |− sector. The state |Ψ + admits a simple computation of its energy expectation, and its variance. In order to compute the variance, we need two ingredients. The expectation value of the energy Ψ + |H L |Ψ + = Ψ + H N −1 The expectation value of the square of the Hamiltonian is where we can explicitly calculate We notice that, in the limit of weak perturbation, the two methods give similar results. Close to the boundary of the convergence radius, higher order term contributions become relevant, giving corrections only at larger distances.