The frustration-free fully packed loop model

We consider a quantum fully packed loop model on the square lattice with a frustration-free projector Hamiltonian and ring-exchange interactions acting on plaquettes. A boundary Hamiltonian is added to favor domain-wall boundary conditions and link ground state properties to the combinatorics and six-vertex model literature. We discuss how the boundary term fractures the Hilbert space into Krylov subspaces, and we prove that the Hamiltonian is ergodic within each subspace, leading to a series of energy-equidistant exact eigenstates in the lower end of the spectrum. Among them we systematically classify both finitely entangled eigenstates and product eigenstates. Using a recursion relation for enumerating half-plane configurations, we compute numerically the exact entanglement entropy of the ground state, confirming area law scaling. Finally, the spectrum is shown to be gapless in the thermodynamic limit with a trial state constructed by adding a twist to the ground state superposition.


I. INTRODUCTION
Ergodicity and its breaking lies at the foundation of modern statistical mechanics.It plays a key role in understanding of why the long-time average of an observable of a single system can be well-approximated by a statistical ensemble average.In quantum systems, any initial state being thermalized necessarily requires each eigenstate of the Hamiltonian to be thermalized, leading to the so-called eigenstate thermalization hypothesis (ETH) [1,2].Over the past few years, ETH violation has been realized outside the scope of the integrability and many-body localization paradigms [3,4], such as in quantum many-body scars (QMBSs) and Hilbert space fragmentation [4,5].
Two types of models have played important roles in understanding these novel mechanisms of weak ergodicitybreaking.The first type is kinetically constrained models, which can arise as low energy effective models through a Schrieffer-Wolf transformation [6,7].The dimensionality of constrained Hilbert spaces typically grows as an integer sequence, reflecting an underlying combinatorial structure.In one dimension, a prime example is the PXP model [8][9][10], which successfully explains Rydberg blockade experiments [11].The dimensionality of its Hilbert space grows as the Fibonacci numbers with the asymptotic scaling 1.618 N .In two dimensions, arguably one of the most studied models in classical statistical mechanics and combinatorics is the six-vertex model.With periodic boundary conditions, its Hilbert space dimension grows as 1.540 N 2 , following from Lieb's solution to the square ice problem at the ice point for which the weights of the six vertices are identical [12].Sophisticated results have been established when the model is subject to domain-wall boundary conditions (DWBCs), for which a bijection between configurations and alternating sign matrices (ASMs) [13] has been proven.The exact enumeration of ASMs is a celebrated result in combinatorics [13][14][15][16].Notably, progress has also been made in technologies and ideas for realizing classical and quantum spin ice models [17], with platforms ranging from arrays of ferromagnetic islands [18] to two-dimensional Rydberg atom arrays [19][20][21][22][23].
A second type of models studied in weak ergodicitybreaking are frustration-free (FF) Hamiltonians, which have a unique ground state being the superposition of configurations from a usually classical statistical mechanical or combinatorial ensemble.Such Hamiltonians are called frustration-free because their ground state is the simultaneous lowest energy eigenstate of all its local terms.Examples of frustration-free models in 1D include the Motzkin [24] and the Fredkin spin chain [25,26], for which ground state configurations reassemble combinatorial structures known as Motzkin and Dyck walks.Recently, it has been shown that by flipping the signs of some of the projectors, the FF eigenstate can be relocated to the middle of the spectrum, making it qualified for a quantum many-body scar.These systems, as well as the original models, also exhibit Hilbert space fragmentation [27][28][29], which refers to the emergence of exponentially many dynamically disconnected subspaces.A classification further distinguishes genuine from "local" fragmentation, as related to the scaling of the Krylov subspaces with system size [30,31].
Entanglement entropy (EE) plays an important role in the study of both types of aforementioned models of novel weak ergodicity-breaking.In the first type of models, the growth of EE is used to characterize the slow thermalization behavior of the so-called scarred initial state.In the second type of models in one dimension, the ground state provides an example of violation of area law.Area law here refers to a ground state EE scaling of S ∼ N d−1 where d is the spatial dimension.A milestone in the study of EE has been Hastings' proof of area law in gapped one-dimensional (1D) systems [32].Recently, a similar result in two dimensions (2D) has been proven for frustration-free models [33,34].The EE of free fermions generally violate area law, but only logarithmically, N d−1 ln N [35].In 1D there are multiple mechanisms that generate beyond logarithmic violation of area law, such as enlarging the degrees of freedom or the local Hilbert space [36], and strong inhomogeneities [37].A combination of these approaches also generalizes extensive entanglement growth to Hausdorff dimension one lattices embedded in higher dimensions [38].One route to generalize area law violation to higher dimensions is offered by the Motzkin and Fredkin spin chains, which are both translationally invariant and can violate the area law strongly, with up to volume-law scaling [39][40][41][42][43].A crucial ingredient of the Motzkin and Fredkin models is that they allow a height representation that can carry non-local information of the interactions.The first obstacle in generalizing this to 2D is to find a well-defined height function that does not give rise to any ambiguity when going around a closed loop in the lattice.This problem is intrinsically avoided in the context of dimer and fully packed loop (FPL) models [44][45][46][47].To further enforce the incremental height change between adjacent plaquettes to be ±1 on square lattice, we opt for FPLs.
In this manuscript, we combine the two mechanisms of weak ergodicity-breaking in a single model in two dimensions and explore its various features, including Hilbert space fragmentation, ground state entanglement entropy, and an upper-bound on the spectral gap.We find that the notion of a height function alone in two dimensions is not enough to generate beyond area-law EE, due to the strong constraint of gauge invariance.A further decoration of the model in the current manuscript, with enlarged local degrees of freedom and next-nearest neighbor interactions, was recently proposed in two models possessing up to volume EE scaling of ground state [48,49].In this manuscript we lay the foundation that make these two generalizations possible, including detailed discussions of ergodicity breaking and its dynamical consequences, why a height model alone is not enough to break area-law (contrary to the 1D case with Motzkin and Fredkin spin chains), and a proof of gaplessness in the thermodynamic limit.The current manuscript also complements the models of Refs.48 and 49, which are more focused on EE area-law violation and phase transitions, from considerations of the global structure of the Hilbert space and excited states.
The paper is organised as follows.In Sec.II, we introduce the Hilbert space and Hamiltonian of our model.In Sec.III we explore ergodicity breaking in detail (Sec.III A), in addition to constructing a series of exact eigenstates (Sec.III B).In Sec.IV we compute numerically the exact ground state bipartite entanglement entropy between half-systems, yielding area law scaling.In Sec.V we construct a trial state with vanishing en-ergy in the thermodynamic limit, demonstrating that the model is gapless.Finally, we provide a summary and propose future directions.In Appendix A we prove that FPL configurations with a fixed boundary are ergodic with respect to plaquette flipping in the bulk.In Appendix B we provide a sufficiency proof for exact unentangled exited states.In Appendix C we explain how half-system configurations can be counted recursively, which is used to calculate the entanglement entropy, and in Appendix D we devise a Monte Carlo algorithm to faithfully sample the classical configuration space.Finally, in Appendix E we provide an example of a 2D frustration-free FPL model that does not suffer from boundary constraints.

II. THE MODEL AND ITS DUAL LATTICE FORMULATION
We consider a square lattice G of N × N vertices, with binary degrees of freedom living on the horizontal and vertical bonds connecting neighbouring vertices.We constrain the Hilbert space to that of the fully packed loop (or 2-dimer) coverings of the lattice by the diagonal Hamiltonian 2 , where subscripts v and h refer to vertical and horizontal bonds respectively, and σ = +1 (σ = −1) for a covered (uncovered) bond.In the limit V 1 we can effectively work in the ground state manifold of H 0 .In the constrained Hilbert space FIG. 1. Bijections with domain-wall boundary conditions.Left: fully packed loop configuration, middle: six-vertex model configuration, right: alternating sign matrix.The mapping between the first two representations is explained in Fig. 2 (a), and the mapping between the latter two amounts to assigning a and b vertices to 0, and c vertices to ±1 [50].
we consider the low-energy effective Hamiltonian The Rokhsar-Kivelson projectors, P p , contain the diagonal potential term  [51].We will refer to plaquettes that contain two parallel covered bonds as flippable, and unflippable otherwise.Here bonds are either covered (black) or uncov-ered (light gray), with the bond-spin conversion rules in Fig. 2 (b).The sum above runs over the bulk plaquettes, which there are (N − 1) 2 of.The boundary terms x,1 ( impose a domain-wall boundary condition (labelled DWBC1) on the ground state, where every other bond is covered along the boundary (see Fig. 2 (c)).Above, subscripts denote the (horizontal or vertical) bond position (x, y), as counted from the lower left of the graph with row-major ordering.With these boundary conditions FPL coverings are also in bijection to alternating sign matrices [52], see Fig. 1.
The bulk Hamiltonian has the apparent Z 2 symmetry R of reversing the covering of all the bonds, satisfying which is broken by the boundary terms: The commutation relation of Eq. ( 3) can be understood by observing that all the projectors P p are invariant upon interchanging covered and uncovered bonds.The anticommutation relation of Eq. ( 4) becomes apparent by noticing that := 1 − = 1 2 (1 + − ), which makes H ∂ − 2N anticommute with R. The full Hamiltonian also has a hidden symmetry given by Wieland's gyration [53], which reverses the coverings around only the non-flippable plaquettes, while leaving the flippable ones unchanged.
The off-diagonal terms in Eq. ( 1) relate FPL configurations differing on a single plaquette with parallel bonds covered in different directions.This is conveniently expressed in the height representation on the dual lattice [44,46].On the square lattice the height field is integer-valued and changes in units of ±1 between neighbouring plaquettes.Using the bipartite rules summarized in Fig. 2 (b) the conversion between spins and height is compactly expressed as S = (S v , S h ) = ∇h, with ∂ x h(i, j) := h(i + 1, j) − h(i, j) and ∂ y h(i, j) := h(i, j + 1) − h(i, j).The FPL constraint then amounts to imposing ∇× S(i, j) = 0 around each vertex.It can easily be verified that these definitions make precisely the (flippable) plaquettes with two parallel covered bonds local extrema in the height, see Fig. 3. On the dual lattice, the model can be expressed as the kinetically constrained Hamiltonian

>(<) p
projects onto the state where the four neighboring plaquettes of p all have the same height above (resp.below) the height of p, and h +(−) p increases (resp.decreases) the height of plaquette p by 2. The boundary term is given by The height representation on the dual lattice reveals another symmetry of the Hamiltonian.The height difference along each row and column (or along any other path connecting two boundary plaquettes), which are equal to the sum of spins N j=1 S v (i, j) and N i=1 S h (i, j), respectively, is conserved.In fact, in the absence of kinetic (off-diagonal) terms on the boundary, which transform one boundary configuration into another [54], there are a linear number of local discrete symmetries: which are responsible for the ergodicity breaking to be discussed in the next section.This makes the model possess "local fragmentation" in the terminology of Ref. [31], to be contrasted with ergodicity breaking due to either discrete or continuous global symmetry such as the total magnetization in spin chains, or genuine fragmentation in models such as caused by pair flipping [30].In Appendix A we prove that the bulk Hamiltonian is ergodic within each Krylov subspace spanned by all the FPL configurations sharing the same boundary configuration.

III. ERGODICITY BREAKING AND EXACT EIGENSTATES
In this section we identify exact unentangled and finitely entangled eigenstates, by making simple observations of constraints on the boundary configurations.Among the exact eigenstates we construct the exact and unique ground state and its twin ceiling state.

A. Product eigenstates and finitely entangled eigenstates
If we follow the height counter-clockwise around the perimeter, the fact that it must come back to the same value after a full cycle implies that there must be at least one pair of height kink and anti-kink along the graph perimeter, as the height of neighbouring plaquettes must differ by ±1.Moreover, if anywhere on the lattice, there is a segment, such as those in the same color in Fig. 4, in which the height changes monotonically, and which has a right-angle turn, the height in the entire rectangle spanned by the two perpendicular sides must also change monotonically between diametrically opposing corners of the rectangle.We will henceforward refer to this as the convexity lemma, and use it repetitively.By "segment" we here refer to any path between two plaquettes with a single right-angle turn, i.e., any "L"-shaped path.The key point is that since flippable plaquettes correspond to local extrema in the height, a monotonic segment can not contain any flippable plaquettes.
An immediate consequence of the convexity lemma is that a monotonic segment on the perimeter cannot cover two corners.Thus, besides the special case of two kinks being located on diagonally opposing corners as described below, there must be at least two pairs of alternating kink and anti-kink, with the sum of lengths of every other segment being 2N , since the height along half of the perimeter must increase, and decrease along the remaining half, the net height change along the entire perimeter (of length 4N ) is zero.A four-segment case is depicted in Fig. 4.
As for the exceptional kink-antikink case, the height change must be monotonic in each row and column inside the bulk, so there are no flippable plaquettes in this case.We have thus found the first product eigenstate, as depicted in Fig. 5 (a) (for which the kink and the antikink are located at diagonally opposing corners).With two kink-antikink pairs, it is still easy to pick boundary configurations that only allow one bulk configuration.First, notice that the only way in which none of the four segments go around a corner is for each of them to cover exactly one side of the lattice, including one of two corner bonds in either end.This gives the ground state boundary configuration depicted in Fig. 2 (c).Otherwise, each segment must go around one corner and one corner only.
If the four rectangles formed by the four monotonic segments cover the entire square lattice, it also results in a product eigenstate.This gives the states shown in Fig. 5 (b), (c), and (d).Together, up to translations and rotations, they exhaust all product eigenstates with two pairs kink-antikink on the boundary.This can be seen by writing each segment length as l i = x i + y i for i = 1, . . ., 4, see Fig. 4. Monotonicity within each rectangular region and alternation of monotonicity between regions require the constraints min(x 1 +x 3 , y 1 +y 3 ) and min(x 2 + x 4 , y 2 + y 4 ) to be either N or N − 1 for the excited state to be a product state.The former case gives (b) and (c), up to translations, and the case latter gives (d), up to translations and π 2 rotations.By now it can be seen that in general, a rule in searching boundary configurations yielding product eigenstates is that each plaquette in the lattice is penetrated by at least one height gradient line.Here gradient line just means any directed path (among multiple choices) along FIG. 4. A boundary configuration with four monotonic segments (two kinks and two anti-kinks), demanding four rectangular monotonic regions near four corners, leaving a core with a single flippable plaquette.The right panel shows the corresponding height going around the lattice boundary counterclockwise.FIG. 5. Representatives of the two types of product eigenstates in isolated one-dimensional Krylov subspaces, with the opaque green color scale reflecting the height.(a), (e), and (f) have straight height gradient lines penetrating to opposite sides of the lattice, while (b), (c), and (d) have a saddle point of height gradients in the bulk.The energy of (a), (b), (c), (e), (f) is always 2N , while the energy of (d) is 2N − 2. Up to translations and rotations, the only remaining product states not depicted here have height functions being monotonic in one direction, and globally shifted parallel random walks in the orthogonal direction.
which the height changes monotonically and should not be confused with the direction of steepest descend.Gradient lines must start and end at boundaries, as loops are not allowed and endpoints imply a flippable plaquette.We now prove that the only possible full coverage of gradient lines, other than the special ones going around four corners as we discussed above, are parallel straight line coverings, corresponding to product states of the type Fig. 5

(e) and (f).
If there is at least one plaquette being traversed by a gradient line that starts and ends on opposite sides, we call this plaquette (i, j), otherwise, it becomes the one of special cases to be discussed in the next paragraph.Assume the plaquette (i, j) is traversed by a gradient line starting from side a and ending on side b.Since a and b are on opposite sides, without loss of generality, we can assume the coordinates of the two end points are (i a , 1) and (i b , N ), with i a ≤ i ≤ i b .We then immediately have that the entire region between i a and i b is traversed by straight gradient lines from the bottom to the top of the lattice, by the convexity lemma.But the lines in the region with first coordinate i ∈ [i a , i b ] will have no where to go except turning to the adjacent side, resulting in either a curl or violation of the convexity lemma.We thus conclude that if we have a product state in which one plaquette is traversed by a gradient line that starts and stops on opposite sides, then all plaquettes must be traversed by straight and parallel gradient lines starting and ending on opposite sides.
If however, all gradient lines turn to end on adjacent sides, an analogous argument will result in product states of the type in Fig. 5 (b), (c) or (d).Hence we have exhausted all the possible unentangled eigenstates.Alternatively, the above argument can be made elegantly in the six-vertex language.As such, a discrete version of electrostatics dictates that gradient lines, coming from strings of unflippable plaquettes, are curl-free and sourcefree.These are precisely what the FPL configurations shown in Fig. 5 correspond to when tracing their height profile.
Apart from product eigenstates, the lowest entangled frustration-free eigenstate corresponding to the boundary configuration in Fig. 4, which is a equal superposition of two FPLs differing only by the center plaquette, with EE of ln 2. One can further explicitly diagonalize a 2 × 2 plaquette system subject to DWBC, which is satisfied by the four-plaquette core surrounded by four monotonic rectangular wings.These states have area law entanglement entropy of the size of the core, which can be anywhere between 1 and N , if the cut is through it.

B. Exact eigenstates
As we prove in Appendix A, each consistent boundary configuration forms a Krylov subspace, in which a uniform superposition of all FPLs with that boundary configuration is an eigenstate.Among these would-be degenerate eigenstates, the boundary Hamiltonian picks the one with DWBC1 to be the unique and global ground state where the normalization constant is equal to the number of ASMs [55] due to the domain-wall boundary conditions realized by the ground state [13][14][15][16], cf.Fig. 1.
The Hamiltonian can be deformed locally with a large class of parameters.One can consider the graph Laplacian that gives the Hamiltonian.That graph turns out to have chordless cycles only of length four, originating from flipping two corner-sharing plaquettes in different orders.The projectors acting on each plaquette p can be deformed in a frustration-free manner independently by a spatially varying angle θ p not subject to any consistency relation of the type in Ref. [41]: Furthermore, it appears possible to deform the Hamiltonian in a neighborhood-resolved manner such that the ground state is a superposition of weighted six-vertex configurations.The detailed investigation of this quantum sixvertex model beyond the ice point, we intend to undertake in future work.
The necessary and sufficient condition for a boundary configuration allowing an integer energy eigenvalue is most succinctly expressed in the six-vertex language: Firstly, the total number of inward and outward pointing arrows must be the same; secondly, if the lattice is divided along any row or column, the difference between inward and outward pointing arrows on either side of the partition line can not exceed N , see Fig. 9 in Appendix B, in which the sufficiency condition is proven with backward induction.The necessity follows straightforwardly from the conservation of arrows at each vertex, so any partition of the lattice must leave the interior arrows along the cut capable of balancing the net arrow flow in the exterior.
An eigenstate with 2n flipped boundary spins relative to DWBC1 (in Fig. 2 (c)), and otherwise being annihilated by the bulk Hamiltonian in a construction similar to the ground state, will have energy E = 2n.As explained above and proved in Appendix B, a consistent boundary condition in the height representation requires color balancing, meaning that there is a balanced number of covered boundary bonds connected to vertices of either sublattice, i.e. n red and n blue vertices connected to covered boundary bonds.The excited state for a consistent boundary condition S n (the set of boundary spins) thus takes the form where the sum runs over FPLs with boundary condition S n .A simple estimate for the degeneracy of level 2n is d 2n ≤ 2N n 2 since there are 2N n ways in which both n blue and n red boundary spins can be flipped from DWBC1.This is an upper bound since it includes a handful of boundary configurations violating the second part of the sufficiency condition explained above, so that global conservation of arrows is respected, but where the arrow conservation in a subsystem including the partition line is violated [56].For n = 0 we recover the unique ground state with E = 0 and S 0 = DWBC1, and for n = 2N we get the unique excited state with E = 4N and S 2N = DWBC2.This E = 4N state is related to the ground state by the Z 2 operator R from Eq. ( 3) and (4).

C. Relation to quantum many-body scars
Here we comment on how our mechanism of weak ergodicity-breaking contrasts and complements the known paradigms of spectrum generating algebra [57] and the Shiraishi-Mori embedding formalism [58].Our Hamiltonian resembles the Shiraishi-Mori formalism, with a target space given by the space of all classical fully packed loops (FPLs).The diagonal boundary Hamiltonian lifts the a priori ground state degeneracy, making the ground state unique on the open square grid.Our exact frustration-free excited states are also equidistant in energy.Yet, the majority of our states are not related by ladder operaters, except for the ground state and the ceiling state, which are related by the gyration operator [53].Furthermore, the energy of our exact excited states are of order N while the spectrum ranges from 0 to order N 2 , so the eigenstates do not have finite energy density in the thermodynamic limit and hence do not qualify for QMBS states.The boundary puts stringent constraints on 2D models in constrained Hilbert spaces, to the point that even product state can become eigenstates, which is not observed in 1D, and could lead to potential applications in quantum technology.

IV. BIPARTITE ENTANGLEMENT ENTROPY
The bipartite entanglement entropy of the half-system defined by the central cut shown in Fig. 2 (c) is given by where m = (h(N/2, 1), h(N/2, 2), . . ., h(N/2, N )) is the height field across the cut.In what follows we will evaluate this for the ground state using an exact recursion relation.

A. Schmidt decomposition and exact recursion
The numbers p m are obtained by Schmidt decomposing the ground state: (12) where |P L, m (N/2) (resp.|P R, m (N/2) ) is the normalized sum of FPLs in the left (resp.right) half-system with heights m imposed along the right (resp.left) edge (and DWBC on the remaining three).The number of different height configurations along the (full system) central cut is given by the number paths with height changes of ±1 across neighbouring plaquettes, starting and ending at the same height of N/2: The coefficient p m above then takes the form where |P R, m (N/2)| is the number of FPLs in the right half-system with the left boundary m.

B. Scaling of the entanglement entropy
The combinatoric nature of the counting problem explained above, related to the rapid growth of A(N ) [55], and the exponential slowdown encountered when solving the recursion relation numerically still pose as a practical challenge.11), as obtained with the exact recursion relation (black squares) described in Appendix C and with the Monte Carlo method (red diamonds) described in Appendix D. The blue dashed line shows a fit to the last three data points.The inset shows the difference between consecutive points, ∆S(N ) = S(N )−S(N −2), which from the upper bound of Eq. ( 15) is expected to saturate to a constant for large N .
In Fig. 6 we show the outcome of calculating the entanglement entropy of the half-system from the exact recursion relation.A Monte Carlo algorithm, which does not rely on ergodicity and which faithfully samples FPL configurations, was devised and verifies the results with excellent accuracy.Details of the Monte Carlo method are provided in Appendix D. The entanglement entropy resulting from the exact recursion relation has area law scaling.Due to the shape of our lattice (having corners) one generally expects subleading corrections like ln N .On the cylinder and torus the general behaviour is a nonuniversal linear term and a universal constant related to the Lifshitz field theory relevant at the quantum critical point [46,61,62].
To shed further light on the scaling we can state a simple upper bound on Eq. ( 11) in the height field representation.The extreme case is obtained for a flat distribution with p m = 1/N m : This shows that the growth of the number of central height field configurations is only enough to make the half-system entanglement entropy upper bounded by N ln 2. The same argument applies to the exact excited states, while a generic excited state would not be subjected to the same Schmidt decomposition of Eq. ( 12) and therefore have EE upper bounded by Lieb's residue entropy of square ice instead [12].Recently, a 2D generalization of the Fredkin model was proposed on the hexagonal lattice with the dimer constraint [63].In Appendix E, we define a similar 2D generalization of the Fredkin spin chain to the square lattice with the FPL constraint.By the logic of the upper bound above this model presumably also faces a ground state bipartite entanglement entropy bounded by area law.The same conclusion is expected to hold for other height-conserving dimer models considered recently [64].
From a tensor network point of view, the ground states of both models can be view as given by projected entangled paired state (PEPS), with virtual legs corresponding to our physical degrees of freedom, and no degrees of freedom on the physical legs.In the tensor network language, the EE scaling is given by the number of bonds crossed by the boundary between subsystems.So to achieve an area law breaking entanglement entropy or a real 2D generalization to the Fredkin model, one may consider constructing a 3D holographic tensor network for which the physical legs only live on the bottom layer, while for the rest of the layers, the virtual leg from the previous layer plays the role of a physical leg [65].

V. UPPER BOUND ON THE SPECTRAL GAP
The classical six-vertex model at the homogeneous point, at which the Boltzmann weights of the a-, b-, and c-type vertices of Fig. 2 (a) are equal, is critical [66], with algebraically decaying bond-bond correlators.By construction, the ground state of our quantum Hamiltonian has the same equal-time correlation function as the classical model, which mean the spectral gap between ground state and first excited state must vanish in the thermodynamic limit [67][68][69][70].In this section, we construct a trial state (Eq.( 16)) to show that the excitation gap is upper bounded by a quantity that decays exponentially with system size.
The trail state needs to meet two requirements: it has to be both orthogonal to the ground state, and the Hamiltonian must have an expectation value that approaches zero in thermodynamic limit for it to provide an upper bound on the excitation gap.The first requirement can be satisfied by adding a "twist" to the ground state superposition [24], such that when taking the inner product with the ground state, the two parts with opposite signs cancel.The second requirement calls for a carefully chosen boundary between the two sets of configurations, such that the Hamiltonian annihilates the intra-set contributions, leaving only contributions from the interface of the sets in the energy expectation value.A convenient choice of the set boundary is halfway between the highest V max = 1  3 N (N + 1)(2N + 1) and lowest V min = 1  3 N (N + 1)(N + 2) volume configurations (as found by summing up the configurations in Eq. (D4)).The two sets contain the same number of configurations, and one has to traverse one of the many such configurations to go from configurations with a volume smaller than their average, , to one with a larger volume, see Fig. 7.We thus pick the trail state It immediately follows that GS|π = 0 and π|π = A(N ).Before evaluating its energy expectation value, we first count the numbers of configurations near the volume interface V 0 .A representative configuration with volume V 0 − 1 is shown in Fig. 7.The flippable plaquettes of this configuration are all located inside a diamond with half the size of the lattice.Outside this diamond the plaquettes are frozen.Inside, every other plaquette is flippable, giving a total number of M := N 2 4 flippable plaquettes.Among these there are M +1 2 plaquettes with height N 2 −1 and M −1 2 with height N 2 + 1 [71].The number of configurations with volume V 0 − 1 can thus be enumerated as since one can simultaneously flip any number of pairs of plaquettes of heights N 2 − 1 and N 2 + 1 from this reference state and remain in vicinity to the volume boundary.Each of these configurations can brought across the volume boundary by M +1 2 Hamiltonian terms flipping one of the plaquettes with height N 2 − 1.We have Using the asymptotic behaviour of A(N ) form Eq. ( 9) [55] we find proving that the Hamiltonian is gapless in the thermodynamic limit.A remark is in order here to point out the connection between our proof of gaplessness and the "arctic curves" of the six-vertex model [72].The vanishing asymptotics above strongly relies on the fact that near V 0 , the region containing flippable plaquettes only occupy half of the lattice and the corners outside the orange rhomboid in Fig. 7 (b) are frozen.

VI. CONCLUSIONS
We constructed a quantum fully packed loop model with a Rokhsar-Kivelson type Hamiltonian [51], in which configurations permit multiple equivalent formulations from the classical statistical mechanics and combinatorics literature.By making the model frustration-free the quantum ground state becomes an equal superposition of configurations from the classical space of configuration.We showed that the bulk configurations are heavily constrained by the boundary, to the point that certain boundary configurations imply product eigenstates.The bulk Hamiltonian is not ergodic in the entire Hilbert space, but only within each Krylov subspace, as dictated by the boundary configuration.Each ergodic subspace has its own lowest energy eigenstate, which are equidistant in energy across subspaces.Owing to enumerable half-system configurations by recursion, we performed an exact lattice calculation of the ground state bipartite entanglement entropy for systems of sizes up to 18 × 18 giving area law scaling.
Our methodology may turn useful in the study of other height models and in related studies of ergodicity breakdown induced by boundary terms.One possible generalization is to consider a Z n generalized model involving a boundary condition that alternates with period n instead of two.One may also expect the emergence of interesting phases and refined structures in the ergodicity breaking for fully packed loop models adopted to non-bipartite lattices.Lattices of interest include the triangular one [73][74][75], or more exotic ones such as the Kagome lattice [76], or even aperiodic tilings like the Penrose [77] and the Amman-Beenker tiling.
One could attempt to construct a microscopic Hamiltonian that makes the FPL constraint emerge, like what was done for the dimer model in Refs.78 and 79 using Klein terms of the Hamiltonian.There are also alternative ways to implement the DWBC, for example, by employing an antiferromagnetic interaction along the boundary.The outcome of this choice, other than making the ground state two-fold degenerate and the exact excited states reordered in energy, is that the entire Hamiltonian will have a Z 2 symmetry.
One can also introduce dynamic terms in the boundary Hamiltonian, so that the fragmentation is removed and the unique ground state becomes the superposition of FPLs with all boundary configurations.It would be of interest to explore the consequences of that on the EE scaling.In addition, it is also interesting to think of whether the Hamiltonian can be modified to reallocate the exact excited states to mid-spectrum QMBS states.
−1, or +1 such that the sum of each column and row is 1, and where the +1 and −1 elements alternate along rows and columns.The mapping to six-vertex model configurations is obtained by assigning c1 vertices +1, c2 vertices −1, and all the remaining four vertices 0 [13,16] In this appendix we prove by induction that each Krylov subspace specified by the boundary configuration (H ∂ ) is ergodic with respect to the bulk Hamiltonian (H bulk ).
For a 2 × 2 lattice, there are at most two choices of the bulk configuration, differing by the height on the central plaquette, for any fixed boundary configuration.These two are related by the flipping the central plaquette.
Suppose we have proven the ergodicity for (the interior of) a lattice of size k×k for any fixed boundary configuration, that is, if such a boundary allows multiple bulk configurations, they can all be transformed into each other by sequences of plaquette flips.Then for a lattice of size (k + 1) × (k + 1), for which the bottom left corner is a k × k lattice with established ergodicity, see Fig. 8, we just need to show that any configuration of the newly added strip above the top and to the right with identical bulk configurations, transform into each other by sequences of plaquette flips.
To show this, we only need to show that sequences of plaquette flips can transform any two strip configurations differing only by the height of any single plaquette into each other, because once that is established, configurations differing at multiple places can be connected step by step transitively.If the difference is located in the corner of the strip, then the two configuration are related by simply flipping the corner plaquette.Otherwise, if the heights h 0 ± 1 differ on a plaquette p 1 , see Fig. 8, on the right boundary, its neighbors above and below must both have height h 0 , or there are at least two plaquettes along the strip with different heights between the two strip configurations in question, contradicting our assumption.
FIG. 8. Proving ergodicity of the bulk Hamiltonian by induction: Two configurations differing only by the height of one plaquette along the strip added in step (k + 1) are either related by flipping that plaquette (p1), or if it is not flippable, and a sequence of flipping operations in the k × k bulk (by its proven ergodicity) will make it flippable.Here, the height configuration of one of the two configurations differing at p1 and p2 is shown, the other configuration has height h0 − 1 on p1 and h0 on p2.
All we need to do now is to show that there is another bulk configuration in the k × k lattice compatible with the strip configuration on its boundary, for which the adjacent plaquette p 2 to the left of p 1 has height h 0 , between the two plaquettes of the strip with height h 0 ± 1. Then these two configurations of the k + 1 × k + 1 lattice with identical heights in the k × k lattice are related by flipping plaquette p 1 .And this must be the case, because if not, p 2 must have h 0 + 2 (resp.h 0 − 2) if the height on p 1 is h 0 + 1 (resp.h 0 − 1), with neighbors above and below both of same height as p 1 , otherwise, the flippability of p 2 implies that there is another bulk configuration with height h 0 on p 2 making p 1 flippable.This argument goes on until we reach at the left boundary, with a monotonic gradient along this row, see Fig. 8.As we established in the exhaustion of product eigenstates, this means the entire lattice must be covered by horizontal straight gradient lines, making the ergodicity proof irrelevant for the boundary configuration under discussion.
In the periodic boundary case, ergodicity was proven in Ref. 84 for the "zero-flux" sector, which contains two configurations with every plaquette flippable, called "ideal states", and they are connected to each other by flipping half of the plaquettes.Since a global extremum implies a local extremum, each configuration has at least two flippable plaquettes which can be flipped to reduce the difference between maximum and minimum heights until one reaches one of the two ideal states.As elegant as this proof is, it does not work for fixed boundary conditions since the extrema can be located on the boundary.Modulo exceptions of configurations with no flippable plaquettes, which we provide the precise criteria for in Sec.III A, the ergodicity in other Krylov subspaces under periodic boundary conditions can be proved without too many modifications of the above proof.
In Ref. 85 it is stated that the "topological sectors" are classified by the height change along two perpendicular non-contractible loop directions in the case of periodic boundary conditions.There are 2N + 1 such sectors and the bulk Hamiltonian is ergodic within each of these.This can be seen by cutting the torus along two orthogonal non-contractible loops, due to periodicity, the boundary configuration are uniquely determined by the height along these two loops, which form an "L"-shaped strip.The kinetic terms in the Hamiltonian can change any configuration along the strip to any other with the same magnetization of arrows pointing in and out.The disconnected sectors can simply be enumerated with a representative where arrows pointing in come before all arrows pointing out, from each magnetization.The number of possible domain wall locations is 2N + 1 for 2N spins.Note that this is different from the the case of open boundary conditions considered above, for which the number of sectors grows exponentially with system size.We caution, however, that the meaning of "topological sectors" here is different from the more conventional usage, the latter in which there are 4 g disjoint sets of configurations for a genus g graph that can not be transformed into each other by sequences of plaquette flips.
Finally, ergodicity of plaquette flipping moves among FPLs or dimer coverings is trivial on the hexagonal lattice, but less understood on the (tripartite) triangular lattice [73,74].In Fig. 11 panel (a) we show the maximal height configuration of this model, and in panel (b) we show the plaquettes that are not flippable in this model, to ensure a positive height.For system sizes of N = 4, 6 we confirmed by explicit enumeration [87] that this model exhibits a larger ground state bipartite entanglement entropy than the model of Sec.II.In fact, this is the case for any straight horizontal or vertical system partition.And even though this Hamiltonian poses as a perhaps more direct 2D generalization of the 1D Fredkin chain than Eq. ( 1), it comes at the price of losing the combinatorial counting technology developed for FPLs with DWBC.And crucially, the simple bound of Eq. ( 15), or more precisely the growth of the number of height configurations, makes it immediately clear that also this model is bounded by area law bipartite entanglement entropy scaling in the ground state.

FIG. 2 .FIG. 3 .
FIG. 2. (a) Mapping between configurations of the six-vertex model (top row) and fully packed loop configurations (FPLs, lower two rows), with rules alternating on the even (cyan dots) and odd (red dots) sublattices [50].(b) Height (on plaquettes) and spin (on bonds) representation of FPLs, which are related as S = (Sv, S h ) = ∇h with discrete derivatives, see the main text.The FPL constraint amounts to imposing ∇ × S = 0 around each vertex.(c) Maximal height configuration for N = 6 with domain-wall boundary condition DWBC1 in the height representation.The dashed green line represents the cut between the left (L) and right (R) half-systems of which we calculate the ground state entanglement entropy in Sec.IV.(d) Boundary conditions labelled DWBC2 relevant for the unique and exact excited state with energy E = 4N described in Sec.III B.

FIG. 6 .
FIG.6.The bipartite entanglement entropy of the halfsystem from Eq. (11), as obtained with the exact recursion relation (black squares) described in Appendix C and with the Monte Carlo method (red diamonds) described in Appendix D. The blue dashed line shows a fit to the last three data points.The inset shows the difference between consecutive points, ∆S(N ) = S(N )−S(N −2), which from the upper bound of Eq. (15) is expected to saturate to a constant for large N .

FIG. 7 .
FIG. 7. (a) Fully packed loop configurations with domainwall boundary conditions partitioned into two sets, dictated by the volume of the configuration.Only configurations in the vicinity of the boundary between the two sets, as measured in a volume metric, contribute to the energy expectation of the trial state in Eq. (16).(b) One of the M M +1 2 configurations (with M = N 2 /4) with volume V0 − 1.The other configurations are obtained by simultaneously flipping an equal number of plaquettes with height 6 and 4.

FIG. 11 .
FIG. 11.(a) Maximal height configuration for the model defined in Eq. (E1).(b) The four plaquettes that are not flippable in the model of Eq. (E1) but that are flippable in the quantum dimer model.

Proof of ergodicity of plaquette flipping within boundary sectors
. [53] B. Wieland, Large Dihedral Symmetry of the Set of Alternating Sign Matrices, The Electronic Journal of Combinatorics 7, 10.37236/1515 (2000).