Area-law entangled eigenstates from nullspaces of local Hamiltonians

Eigenstate thermalization in quantum many-body systems implies that eigenstates at high energy are similar to random vectors. Identifying systems where at least some eigenstates are non-thermal is an outstanding question. In this work we show that interacting quantum models that have a nullspace -- a degenerate subspace of eigenstates at zero energy (zero modes), which corresponds to infinite temperature, provide a route to non-thermal eigenstates. We analytically show the existence of a zero mode which can be represented as a matrix product state for a certain class of local Hamiltonians. In the more general case we use a subspace disentangling algorithm to generate an orthogonal basis of zero modes characterized by increasing entanglement entropy. We show evidence for an area-law entanglement scaling of the least entangled zero mode in the broad parameter regime, leading to a conjecture that all local Hamiltonians with the nullspace feature zero modes with area-law entanglement scaling, and as such, break the strong thermalization hypothesis. Finally, we find zero-modes in constrained models and propose setup for observing their experimental signatures.

Introduction.-Eigenstate thermalization hypothesis (ETH) [1, 2] provides a specific mechanism for thermalization in isolated quantum many-body systems. ETH suggests that the eigenstates of the Hamiltonian at a given energy density are indistinguishable by local measurements and resemble random vectors. A particular consequence of ETH is that highly excited states of quantum system feature strong entanglement. Numerical studies demonstrated that ETH can describe the vast majority of quantum systems [3]. At the same time, possible mechanisms leading to violations of ETH are a subject of active research. Typically, ETH can be avoided due to the emergence of additional conserved quantities that may originate from special properties of Hamiltonian in integrable models [4] or from the presence of strong disorder in the many-body localized phase [5][6][7].
While in integrable and localized systems all eigenstates disobey ETH, recently the focus shifted to systems with weak ergodicity breaking, which have a small number of weakly entangled eigenstates coexisting with the bulk of "thermal" eigenstates that obey ETH. These weakly entangled and thus non-thermal eigenstates were dubbed quantum many-body scars (QMBS) and were reported in a number of different models [8][9][10][11][12][13][14][15][16], see also Ref.
[17] for a recent review. Interestingly, a large fraction of scarred systems features an exponentially large in system size nullspace which is protected by the symmetries of the model [18][19][20][21]. The relevance of such nullspaces to the weak ergodicity breaking was suggested by Ref. [22] which analytically constructed a particular eigenstate from the nullspace (zero mode) of socalled PXP model [19,23,24] as a matrix product state (MPS). Similar zero modes were also discovered in two-dimensional models [25,26] and models with larger blockades [27], while Ref. [28] proposed a systematic way of constructing parent Hamiltonian for MPS zero modes. The MPS form of these zero modes implies an area-law scaling of entanglement. Thus, such zero modes can be regarded as QMBS and provide an example of weak ergodicity breakdown. Moreover, in some cases they could be utilized as a "vacuum" for the construction of other QMBS states outside of the nullspace [22,29]. However, despite weakly entangled zero modes were established for certain models and in many-body localized systems [30], the general conditions for their existence remain unclear.
In this work we explore the structure of the exponentially degenerate nullspaces in a large class of spin chains. We analytically construct a MPS zero mode for a broad class of two-local Hamiltonians with symmetry-protected nullspace. For more generic Hamiltonians, we use a numerical algorithm to construct a basis in the nullspace that is ordered according to entanglement entropy [31]. We define a notion of a least-entangled zero mode, that is shown to obey area-law entanglement scaling in a family of generic Hamiltonians, even though the majority of states in this basis features volume-law entanglement scaling, signaling thermalization [20]. Thus, we conjecture that all local Hamiltonians with an exponentially degenerate nullspace feature a zero mode with area-law scaling of entanglement entropy, establishing a generic route to QMBS and weak ergodicity breaking. Finally, we find the MPS zero modes in a kinetically constrained Hamiltonian and propose an experimental scheme to observe their effects in Rydberg atom arrays [32].
Exponentially degenerate nullspace.-A generic mechanism leading to the exponentially degenerate nullspace of local Hamiltonians is provided by the existence of spatial inversion symmetry and the symmetry of the many-body spectrum around zero energy [19,20]. For concreteness, we consider one-dimensional, inversion-symmetric, spin-1/2 chains with Hamiltonian (1) parametrized by two constants, a and b. Operators X i , Y i , Z i correspond to σ x,y,z i Pauli matrices operating on the local Hilbert space of spin i spanned by ↑, ↓ states. We assume L to be even and use periodic boundary conditions identifying spins L + 1 and 1. The Hamiltonian (1) anticommutes with the operator Π =  Figure 1. (a) Nullspace basis from exact diagonalization (dark blue) has uniformly high entropy similar to the entropy of the eigenstates with energy near zero (green) and close to the Page value (dashed line) [33], SP = L/2 ln 2 − 1/2 . The LENB (light blue) reveals low-entangled zero modes. Inset: Density of states. Data is for ZXZ model with L = 18 and the parameters denoted by a star in (b). (b) Entanglement of |LE1 changes smoothly in the broad range of parameters of ZXZ model with L = 18. The dotted lines denote the regions where |LE1 is a product state. (c) Scaling of entanglement of the LEZM with system size L is consistent with area-law. Data is shown for four different points in parameter space, corresponding to symbols of the same shape in (b).
{H ZXZ , Π} = 0. Π ensures that for each eigenstate |E at energy E, a partner eigenstate Π |E = |−E exists at an energy −E, resulting to a symmetric energy spectrum around zero energy. We note that this property holds for any Hamiltonian which contains terms with odd number of X, Y operators and an arbitrary number of Z operators.
The existence of a degenerate nullspace in model (1) and its generalizations is guaranteed by the presence of spectral reflection and inversion symmetries. While this is basis-independent statement, it is easiest to understand in a computational basis, since product states with even (odd) number of ↓-spins correspond to eigenvalues Π = 1 (Π = −1) respectively. Organizing basis elements into blocks with Π = ±1, the relation {H ZXZ , Π} = 0 implies the block-off-diagonal structure of the Hamiltonian in such basis. Presence of inversion symmetry leads to mismatch in number of basis states between different blocks provided they are restricted to a particular inversion sector, which results in a non-empty kernel. In particular, inversion-symmetric product states necessarily have even number of ↓-spins, and there are 2 L/2 such states. An explicit calculation [34] shows that both inversion-even and -odd sectors feature at least 2 L/2−1 zero modes with even/odd number of ↓-spins. Summing up these contributions, we obtain a lower bound for the dimension of the nullspace, D 0 = dim(ker H ZXZ ) ≥ 2 L/2 [20,23,34].
Analytic MPS zero mode.-We find that all spin-1/2 Hamiltonians of the form H To construct the solution for this equation we use the two dimensional nullspace of the h i,i+1 operator, spanned by a singlet |↑↓ − |↓↑ and a state θ |↑↑ + |↓↓ , where θ = (2a − 1)/(2a + 1) for H ZXZ with b = 0 (see [34] for more generic Hamiltonians). The following choice of local matrices effectively combines these on-sites nullspaces allowing to satisfy the condition h i,i+1 |ψ = 0 for any i, thus giving a MPS zero mode. The existence of a MPS zero mode for H ZXZ with b = 0 and its two-spin generalizations that include Y matrices opens the question regarding the fate of zero modes in more general Hamiltonians.
Least entangled nullspace basis.-A systematic investigation of the nullspace is complicated due to its degeneracy and the absence of a natural basis. In order to overcome this limitation, we use entanglement to construct an unambiguous least entangled nullspace basis (LENB), in which vectors are ordered according to bipartite entanglement entropy, S = −Tr A ρ A ln ρ A . The reduced density matrix ρ A = Tr B |ψ ψ| is obtained by tracing the right half of the chain B. The LENB is constructed in an iterative procedure: first we calculate the least entangled zero mode (LEZM), i.e. a superposition of all zero modes |LE 1 = D0 n=1 c n |n , where n |c n | 2 = 1, that has the least possible amount of entanglement. This is achieved by employing algorithms [31,35] which minimize the entanglement of a vector in a subspace, here chosen to be the nullspace. The resulting state, |LE 1 may be viewed as an analogue of the ground state in the nullspace. Once obtained, the state |LE 1 is projected out of the nullspace and the entanglement minimization algorithm is applied again to the remaining states resulting in |LE 2 . The iteration of this process results in the LENB, {|LE n }, n = 1, . . . D 0 , in which the states are ordered according to their bipartite entanglement entropy.
We numerically construct the LENB for the H ZXZ model (1) at generic values of parameters. Figure 1(a) shows that the ZXZ model for L = 18 chain in the inversion-symmetric sector of zero total momentum has a Gaussian density of states with a small peak at E = 0 corresponding to nullspace. We explicitly check that the model is non-integrable and shows Wigner-Dyson level statistics [34]. The initial basis for the nullspace {|n }, obtained by an exact diagonalization algorithm, consists of states with approximately the same entropy, and almost coincides with the entropy of finite energy eigenstates, see Fig. 1(a). In contrast, the LENB construction results in a small number of weakly entangled states. In what follows we focus on the systematic analysis of LEZM and its entanglement scaling.
LEZM phase diagram.- Figure 1(b) shows the entanglement of the LEZM as a function of parameters of ZXZ model (1) for L = 18. The parameters space features two special lines b = −1 ± 2a for which the LEZM is a product state, |↓↓ . . . and |↑↑ . . . respectively. These lines include the point a = b = 1 which corresponds to the kinetically constrained PXP model [18]. When b = 0, the ZXZ model reduces to a sum of two-site operators for which we constructed an MPS zero mode with bond dimension χ = 2 in Eq. (2), thus implying an area-law entanglement bounded as S ≤ ln 2. These three lines in parameter space correspond to local minima in the entanglement of the LEZM as constructed by the numerical algorithm. The entropy changes smoothly around these minima which suggest the persistence of area-law entangled zero modes beyond the set of lines where analytical results are available.
We study the scaling of entropy with the system size for a wide range of model parameters in Fig. 1(c). For all simulated parameters the behavior of entanglement entropy is consistent with area-law scaling. In particular for parameters that are close to the special lines in the phase diagram entanglement does not change significantly with L. For other values of parameters, finite size effects are more pronounced, yet the finite size scaling is consistent with area-law and corrections decaying algebraically or exponentially with L [34]. Crucially, for all parameters, the state |LE 1 is locally similar between different system sizes as witnessed by the fidelity between local density matrices, and it features a large entanglement gap in entanglement spectrum [34].
The existence of area-law entangled LEZM in a broad parameter regime in ZXZ model raises the question if it is a simple consequence of the exponentially degenerate nullspace, or if the locality of the Hamiltonian is essential for its existence. This is addressed by comparing our local Hamiltonian to a random matrix Hamiltonian with similar symmetries. In this case, we observe that the least entangled zero mode follows a volume law scaling S ∝ L, see [34]. In addition, we show that the distribution of the entanglement spectrum of |LE 1 for a random matrix approaches Marcenko-Pastur distribution [36]. This result shows that for a random matrix the LENB construction does not lead to states which are qualitatively different from random vectors. The drastic difference in the behavior between random matrices and local Hamiltonians implies that the area-law LEZM is related to the locality of the Hamiltonian. Eigenstate thermalization in LENB.-We shift the focus from LEZM to characterizing the whole LENB from the perspective of thermalization. We compare the expectation value of the average magnetization, O = (1/L) i=1 Z i between the LENB states and non-zero energy eigenstates of H ZXZ . ETH suggests that expectation values of the operator O in eigenstates E|O|E are a smooth function of energy up to small fluctuations that are suppressed with system size [3]. Figure 2(a) shows that as energy E approaches zero, the expectation values of observables in eigenstates concentrate around zero, while entanglement entropy is rapidly increasing. In contrast, the states from the LENB defy the expectations from ETH despite having energy E = 0: a significant number of states from the LENB have both anomalously large expectation values of magnetization and small values of entanglement.
To understand the difference between zero modes we explore the finite size scaling of entanglement for all LENB states, see Fig. 2(b). We observe a two component behavior, where |LE n states with small n have a sub-volume law entanglement scaling, while the rest of the LENB states display a volume law scaling. This suggests that even though the nullspace as a whole will display a thermal behaviour in agreement with previous results [20], it generically hosts a small number of exceptional non-thermal states. LEZM in P P XP P model.-We illustrate an existence of area-law zero mode in a constrained spin-1/2 model where P i = (1 − Z i )/2 is the projector to the ↓-state and we restrict to the subspace where ↑ spins are separated by at least 2 sites. This Hamiltonian corresponds to the idealized description of Rydberg atom chains with range-2 blockade [37]: while in ↓ environment any given spin performs free Rabi oscillations, the presence of a nearest or next nearest neighbor ↑-spin arrests the dynamics. As the Hamiltonian has inversion symmetry and anticommutes with operator Π, it features a nullspace as described previously. The entanglement minimization in different momentum sectors reveals a number of lowentangled states in the LENB [34]. A particularly simple LEZM can be written analytically using two-site singlet Experimental signatures of LEZM.-The P P XP P model can be approximately implemented by the triangular ladder of Rydberg atoms, see Fig. 3 inset. Viewing the Rydberg atoms as spin-1/2 degrees of freedom, the system is governed by the Hamiltonian, where n i = (1 + Z i )/2 projects onto ↑ that corresponds to the excited state of a Rydberg atom. The long range interactions between excited Rydberg atoms decay with distance, allowing to find a nearest neighbor atom spacing r, such that V nn = V /r 6 ω. The zigzag geometry in Fig. 3 inset leads to the equal strength interaction V nn between atoms i and i + 1, i + 2. Therefore, for V nn ω the effective range-2 blockade condition emerges and the Hamiltonian (4) can be transformed into perturbed P P XP P Hamiltonian using Schrieffer-Wolff transformation [34].
The perturbations to the H P P XP P include the spin hopping terms and also longer range interaction terms δH = V 3 L i=1 n i n i+3 . Crucially, these terms do not anticommute with Π (since interaction terms in Eq. (4) contain only n, or equivalently, Z operators), thus lifting degeneracy of nullspace. While the preparation of |S is within the limits of current experiments in Rydberg arrays [32], it is an exact eigenstate only for Hamiltonian (3). Thus, while |S remains invariant under unitary dynamics generated by H P P XP P , we investigate its fate in time evolution under H Ry Hamiltonian.
To this end we use a Trotter-based algorithm to evolve |S with H Ry for experimentally realistic value V nn = 2.5ω, leading to a weak but yet considerable range-3 interaction terms, V 3 /ω = 2.5/( √ 3) 6 ≈ 0.1 that cause splitting of the nullspace [32]. Figure 3 shows the evolution of the fidelity F = | S|e −iHRyt |S) | 2 and local Rydberg excitation densities n(t) for a 24-atom system within the experimentally accessible timescale. We observe that although fidelity decreases from one, it remains of the order of 0.5 even at long times. Likewise, the local densities deviate from their original values, but remain far from their equilibrium value according to the microcanonical ensemble. Such dynamics, signals that despite the presence of significant perturbations that destroy nullspace, the initialization of the system in the LEZM of idealized P P XP P model results to a very slow thermalization.
Discussion.-We conjecture the existence of LEZM with area-law entanglement for generic local Hamiltonians with exponentially degenerate nullspace. This conjecture is supported by an analytic construction of a LEZM in the form of a MPS for a particular class of Hamiltonians, and by numerical constructions of LEZM states in the broad parameter regime. Moreover, we demonstrate the existence of LEZMs in kinetically constrained models, whose presence can be probed using Rydberg atom arrays.
These results suggest that slightly entangled zero modes are much more common than previously thought, suggesting that the presence of a nullspace in a local Hamiltonian may be sufficient for the existence of QMBS, thus inviting the systematic studies of nullspaces. It would be interesting to understand the general conditions for the existence of a LEZM that can be represented as MPS with finite bond dimension, and extend these results to higher dimensions using projected entangled pair states representation [38]. From a numerical perspective, the existence of zero-energy eigenstates with area-law entanglement invites the development of efficient numerical algorithms based on MPS that may be able to find such states for system sizes that are beyond the reach of exact diagonalization, or even directly in thermodynamic limit. Finally, area-law entangled LEZM may be used as "ground states" for creating anomalous eigenstates outside of the nullspace using local operators [22,29]. Understanding the conditions for a zero mode to provide a vacuum for stable quasiparticles could result to a novel mechanism of thermalization breakdown at finite energies.
Acknowledgments.-We acknowledge useful discussions with V. Gritsev and A. Garkun and suggestions on implementation of P P XP P model by D.

APPENDIX Appendix A: Lower bound on number of zero modes
In this section we calculate the lower bound in the number of zero modes for a one-dimensional Hamiltonian with spectral reflection and spatial inversion symmetries. We focus on even system sizes. The same calculation for odd system sizes results to a vanishing lower bound. We choose I to be the generator of inversion symmetry and Π = The lower bound in the number of zero modes is given by the mismatch in the number of columns and rows of H e and H o , To calculate this bound we evaluate the dimension of each sector. The basis states in the (e/o, ±) sectors are written as |O e± = |s ± + I |s ± , |O o± = |s ± − I |s ± , where |s ± is a computational state (and therefore an eigenstate of Π) with even/odd number of ↓-spins. To calculate the dimensions of each sector we need to take into account the set of inversion symmetric computational states |s sym = I |s sym . These states are structured as |s sym = |m 1 . . . m L/2 m L/2 . . . m 1 where m i =↑, ↓, e.g. for 4 sites, |↑↓↓↑ is a symmetric state. We observe that there are M = 2 L/2 such states and they always have an even number of ↓-spins. Therefore, the dimensions of the two sectors with even number of ↓-spins are In this section we analytically calculate exact zero mode states which can be represented as a matrix product states of bond dimension χ = 2 for a large class of two-local Hamiltonians including the ZXZ model presented in the main text for b = 0. We focus on spin-1/2, translation-invariant, so-called "two-local" Hamiltonians H = i h i,i+1 that can be written as a sum of operators acting on just two sites, for which the Hamiltonian density features spatial and spectral reflection symmetries, Due to the presence of spectral and spatial reflection symmetries the Hamiltonian density operator h i,i+1 has two zero modes. To understand the structure of these zero modes we write the Hamiltonian density in the eigenbasis of the generator of reflections I and Π i . Since [I, Π i ] = 0 and I 2 = Π 2 i = 1, the eigenbasis is labeled by two binary quantum numbers (i, π i ). We use the basis consisting of triplet and singlet states with quantum numbers, {(1, 1), |↑↑ }, Due to the algebraic constraints of Eq. (B1), the application of the Hamiltonian density to a vector with well-defined quantum numbers changes those numbers as (i, π i ) → (i, −π i ). This leads to the following general structure of the Hamiltonian density in the above defined basis, Solving for two (unormalized) zero modes in different I sectors we get |1 = − c2 c1 |↑↑ + |↓↓ and |2 = |↑↓ − |↓↑ . To calculate the zero mode of the full Hamiltonian we propose a translation invariant matrix product state (MPS) ansatz, where the local tensors A have two virtual indices a ∈ {1, χ} and a physical index s ∈ {1, d} with d = 2 being the local Hilbert space dimension. For simplicity we chose periodic boundaries but the results also hold for open boundary conditions. As it is mentioned in the main text, a sufficient condition for a state to be a zero mode of H is that the matrix hAA obtained from the action of local Hamiltonian on two sites in the MPS vanishes, Indeed the above condition corresponds to h i,i+1 |ψ = 0 and thus, leads to the whole MPS state being zero mode, H |ψ = 0 due to translation invariance. This implies that the tensor hAA (sisi+1) aiai+2 , where the physical indices are vectorized vanishes if its matrix elements AA a,b are superpositions of the zero modes of the Hamiltonian density (|1 , |2 ) or zero. The existence of a non-trivial solution depends on the structure of the zero mode subspace of the Hamiltonian density. For the Hamiltonian density of Eq. (B3), a solution has the form, This means that |ψ is an area-law entangled zero mode (S ≤ ln χ = ln 2). We leave the detailed analysis of this class of MPS, generalizations to higher spins, extensions to different Hamiltonian density nullspaces and more general constraints to future work.

Appendix C: Entanglement minimization algorithm
In this section we present the numerical algorithm used to generate the least entangled nullspace basis (LENB). The algorithm is an implementation of the entanglement minimization scheme developed in [31], which is based on the minimization of Renyi entropy,  1 Choose a state from the subspace |0 ∈ H 0 which consists of zero modes generated by the exact diagonalization algorithm.
2 Apply the singular value decomposition to the state, Eq. (C2), and replace S → S where S nn = λ 2α−1 n .
3 Project back the state to the subspace and normalize it. Repeat steps 2-3 until entanglement entropy converges to a fixed point.
To minimize the entanglement entropy (α → 1 + ) we start by minimizing the Renyi entropy for some α > 1 and slowly decrease α when the entanglement entropy saturates to a minimum. In particular, we find that α n → 1 + (α n−1 − 1)/2 with an initial α 0 = 2, where n ∈ {1, N } is the iteration index, provides an efficient formula to vary α. The saturation of entanglement entropy S 1 is determined from the standard deviation SD N loc (S 1 ) over the last N loc = O(10) iterations. An iteration dependent threshold n = 0.1/10 n , is used to identify whether the entropy saturated, SD N loc (S 1 ) ≤ n . The algorithm is considered to have converged when SD N loc (S 1 ) ≤ 10 −4 , which typically happens after N = O(100) iterations. For large Hilbert spaces ∼ 10 4 we fine tune the parameters to achieve optimal results and run the algorithm for many different initial states to be sure that it converges to the global minimum. We note that since we focus on the nullspace at specific momentum sectors, the entropy will be minimized automatically for all partitions which are translations of the A,B partition.
To calculate the LENB we add an additional step to the algorithm: Following the entropy minimization of a state, we project it out of the subspace. Running the algorithm using the new subspace will generate a least entangled state which is orthogonal to the previous one. Repeating this process will result to an orthonormal set of least engangled zero modes, i.e. the LENB. scaling of entanglement entropy as a function of the system size for the least entangled state in the nullspace of such random matrix ensemble. We observe that the entanglement entropy scales as S ∝ L for every state in the ensemble. In addition, the distribution of entropies becomes sharper as the system size increases, an indication of thermalization of the least entangled zero mode independently of the random matrix parameters.
To further examine the structure of entanglement we study the entanglement spectrum (i.e. the eigenvalues of the reduced density matrix) of a particular LEZM in Figure 4(c) and compare it to the Marchenko-Pastur (MP) distribution which is the distribution of the singular values of a random matrix. We observe that the entanglement spectrum flows towards a distribution which is close to MP distribution, further confirming the generic structure of entanglement of the least entangled zero mode of a random matrix Hamiltonian. These results imply that the presence of weakly entangled zero modes in our systems is not just na artefact of the size of the nullspace ∼ 2 L/2 and that the locality of the Hamiltonian is a critical ingredient for the existence of weakly entangled zero modes.
We observe that r ∼ r GOE for all generic parameters a, b of the ZXZ Hamiltonian. We also observe that the probability distribution function is sufficiently close to the GOE prediction, showing no signs of enhancement of P (r) for small r that typically stems from the absence of level repulsion and is characteristic of non-chaotic systems. In addition to the level-statistics we explore the bipartite entanglement entropy of the eigenstates of the ZXZ model, Figure 5(c). For large enough systems, all eigenstates follow an inverse parabola which is expected from thermalizing eigenstates and a Gaussian density of states. This illustrates the uniform thermalization of the model and the absence of scarred eigenstates [18] that would be visible as "entanglement outliers" and are sometimes present in chaotic quantum systems. For the rest of this section we focus on the least entangled zero modes shown in Fig. 1(c) of the main text and give additional numerical evidence of the area-law scaling of the least entangled zero modes. In Figure 6(a) we show that the finite size corrections to the area-law are well fitted by S(L) ∼ c 1 − c 2 /L 2 . We note however, that that exponential fit S(L) ∼ c 1 + c 2 e −c3L is also able to describe the saturation of entanglement.
In Figure 6(c), we compare the local structure of the least entangled zero mode for different system sizes. To achieve this we construct the reduced density matrices  Fig. 1(c) of the main text. The color gradient denotes the size of the system with L = 8 being the most transparent and L = 20 being the least transparent curves. The curves tend to collapse for large system sizes, indicating that the entanglement is generated only in the boundary between the subsystems. Moreover a jump that can be identified as "entanglement gap" separating large and small singular values is clearly visible in all plots.
chains of different lengths L. We calculate the fidelity of the density matrices for adjacent system sizes f ρ L−2 A , ρ L A , where f (σ, ρ) = tr[ √ ρσ √ ρ] 2 . We observe that the fidelity approaches one as a system size increases, which implies that the algorithm converges to the same local state for all system sizes.
To further examine the fixed points of these states, we shift our comparison from local subsystems to global subsystems. In Figure 7 we compare the different eigenspectra of the half-system reduced density matrix ρ A for different system sizes. We observe that for all parameters, the largest eigenvalues tend to flow towards some fixed point as we increase the system size. This is a strong indication that entanglement is only generated at the boundaries of the subsystems (area-law).
Finally, we have observed that for some parameter points, the entropy saturates for small system sizes and then starts increasing again, see Figure 8. A careful examination of the parameter regimes where this behavior is present, leads to a hypothesis that this increase of entropy is not a themodynamic feature as it leads to a second saturation point at higher entropy visible for several curves in Fig. 8.

Appendix F: Rydberg Hamiltonian and PPXPP model
In this section we derive the formal relationship between the Rydberg Hamiltonian and the PPXPP model defined in Eqs.(4)-(5) of the main text. The Rydberg blockade mechanism arises in the limit of strong nearest-neighbor interactions,V ω, such that the many-body Hilbert space is split into disconnected sectors distinguished by the total number of nearest-neighbor excitations. We employ Schrieffer-Wolff (SW) perturbation theory to address the connection between the Rydberg and kinetically constrained Hamiltonians (H R , H C ) given by Eq. (5) and Eq. (4) of the main text. The SW expansion of order l corresponds to approximations of O(ω l /V l−1 ). We perform the leading order expansion (l = 1) which implies an approximation H 1 = H R + O(ω 2 /V ). We find that for some parameters the entropy grows, flowing towards a larger value before saturating. In some cases the saturation may occur for system sizes beyond our numerical capabilities, thus suggesting that the growth of entropy does not rule out existence of area-law entangled zero modes.
To build the expansion we split the Rydberg Hamiltonian as, where H 0 is the unperturbed Hamiltonian, V is the perturbation and (i − 2, i − 1, i + 1, i + 2) are the four nearest neighbors of site i in the zig-zag lattice shown in Figure 3 of the main text. To derive the expansion, the perturbation is split using the generalized ladder operators, V = The entanglement spectrum of the black circled states for which the density matrix is full rank converges to a fixed point distribution for large system sizes.
illustrate the connection between the exact zero mode and the LEZMs we focus on the LEZM of |π/2| + sector as indicated by the red circles in Figure 9.a. This zero mode is a numerical approximation to the analytic zero mode, which is a zero mode of the k = |π/2| − sector if L = 8n and k = |π/2| + sector if L = 8n − 4, where n is a positive integer. For this state, the reduced density matrix for the half-system bipartition has two finite eigenvalues λ 2 1 = λ 2 2 = 0.5. In Figure 9(b) we show that the numerical algorithm converges with very high accuracy to the analytic zero mode |S which happens to be the least entangled state in that subspace. This result suggests that the entanglement minimization algorithm can also be used to find analytical zero modes as long as they are not highly entangled. In addition, the precise convergence indicates that the entanglement minimization algorithm converges to a global minimum for all available system sizes.
Zero modes with full rank reduced density matrices are also present in the system, see for example the black circles in Figure 9(b). Similarly to the ZXZ Hamiltonian we find that these zero-modes are area-law entangled and their entanglement spectrum flows towards some fixed point distribution, Figure 9(c).