Coupled Fredkin and Motzkin chains from quantum six- and nineteen-vertex
models

Zhao Zhang1$\star$ and Israel Klich2

1 SISSA and INFN, Sezione di Trieste, via Bonomea 265, I-34136, Trieste, Italy

2 Department of Physics, University of Virginia, Charlottesville, VA, USA

⋆<EMAIL_ADDRESS>
## Abstract

We generalize the area-law violating models of Fredkin and Motzkin spin chains
into two dimensions by building quantum six- and nineteen-vertex models with
correlated interactions. The Hamiltonian is frustration free, and its
projectors generate ergodic dynamics within the subspace of height
configuration that are non negative. The ground state is a volume- and color-
weighted superposition of classical bi-color vertex configurations with non-
negative heights in the bulk and zero height on the boundary. The entanglement
entropy between subsystems has a phase transition as the $q$-deformation
parameter is tuned, which is shown to be robust in the presence of an external
field acting on the color degree of freedom. The ground state undergoes a
quantum phase transition between area- and volume-law entanglement phases with
a critical point where entanglement entropy scales as a function $L\log L$ of
the linear system size $L$. Intermediate power law scalings between $L\log L$
and $L^{2}$ can be achieved with an inhomogeneous deformation parameter that
approaches 1 at different rates in the thermodynamic limit. For the $q>1$
phase, we construct a variational wave function that establishes an upper
bound on the spectral gap that scales as $q^{-L^{3}/8}$.

###### Contents

  1. 1 Introduction
  2. 2 Review of the colored Fredkin and Motzkin chains
  3. 3 Six-vertex construction of coupled Fredkin chains
  4. 4 Scaling of entanglement entropy
  5. 5 Scaling of spectral gap
  6. 6 Nineteen-vertex construction of coupled Motzkin chains
  7. 7 Conclusions
  8. A Ergodicity of the Hamiltonian and uniqueness of ground state

##  1 Introduction

Entanglement entropy (EE) and its scaling has been a central theme of quantum
many-body physics [1], not only because entanglement is a unique feature in
the quantum world by itself, but also for their crucial role in determining
the computational complexity of the numerical simulations of quantum many-body
systems [2], indication of topological order [3] and understanding of the
holographic principle and black hole entropy [4]. While EE of a generic
eigenstate in the Hilbert space is shown to scale with the systems size [5],
EE of the ground states of gapped local Hamiltonians are generally observed to
obey the so-called area-law, scaling with the size of the boundary. A
milestone in the study of area-law has been Hastings’ rigorous proof of the
result in one-dimensional systems [6]. Recently, a similar result in two-
dimension has been proven for frustration-free models [7, 8]. While area-law
has been ubiquitous in gapped systems, plenty of examples of area-law
violation has also been found in various gapless systems. (1+1)-dimensional
critical system described by a conformal field theory has EE of logarithmic
scaling [9]. EE of a system consisting of free fermions with a Fermi sea in
dimension $d$ scales as $L^{d-1}\log{L}$ [10]. On the other hand, violations
beyond logarithmic have only been known in one dimension so far.

Quantum vertex and height models are an invaluable tool for the description of
phases of strongly correlated systems [11, 12, 13, 14]. They often emerge as
an efficient description of quantum dimer models where strong local
constraints facilitate the existence of a well defined “height" degree of
freedom. In such models, a ground state may be well described in terms of a
height field and its fluctuations. When coupled to other local degrees of
freedom in such a way that the height field remains single valued, it is
possible to enrich the model to use the fluctuating height field in order to
further mediate correlations. One of the most spectacular examples of such a
behavior is exhibited in the colored Motzkin and Fredkin spin chains [15, 16,
17, 18, 19, 20, 21], where the height degree of freedom can assist in
generating an extensive entanglement entropy in the ground state. Such ground
states thus exhibit a maximal violation of entanglement “area law”. It is
important to note that the degrees of freedom associated with the height
field, due to continuity constraints, although playing a crucial role in
facilitating the entanglement contribution from the color degree of freedom to
be discussed in this work, cannot solely reproduce such area violation in
higher dimensions by itself [22]. In this paper we construct a bicolor six-
and nineteen-vertex models that admit exactly such behavior, in analogy with
the recent lozenge tiling based model we have presented [23]. Our models are
frustration free, with ground states being superpositions of surfaces with
colorings, when viewed along a horizontal or vertical direction, obeying the
coloring rules of arrays of colored Fredkin or Motzkin spin chains.

The paper is organized as follows. In Sec. 2, we quickly review the definition
of Fredkin and Motzkin chains in one dimension and their common entanglement
phase diagram. In Sec. 3, we first introduce the six-vertex construction of
coupled Fredkin chains, with the Hamiltonian and its ground state explicitly
written. In Sec. 4, the EE scaling of the ground state is extracted from a
field theory description of the random surfaces in the ground state
superposition, showing an entanglement phase transition of the $q$-deformation
parameter. In Sec. 5, an upper bound on the spectral gap of the highly
entangled phase is provided with a variational wave function. Sec. 6 sketches
a similar nineteen-vertex construction of coupled Motzkin chains with similar
EE scalings based on the previous sections. Finally, a summary and discussions
of future direction are given in Sec. 7.

##  2 Review of the colored Fredkin and Motzkin chains

Figure 1: Two configurations of a colored Fredkin chain of length $8$ that
differ by the two local configurations related by the projector
$\ket{F^{{\mathrm{r,r,b}}}_{2,6}}\bra{F^{\mathrm{r,r,b}}_{2,6}}$, which appear
in the ground state superposition with a relative weight of $q$.

The colored Fredkin spin chain [19, 20, 21] has a local Hilbert space of
spin-$\frac{1}{2}$ with two colors (red and blue) for up and down spins. It
has a unique ground state as a superposition of colored Dyck walks/paths,
which are random walks starting and ending at the origin and staying on one
side in between, when spin up and down are mapped to an up and down move
respectively, as depicted in Fig. 1. Its parent Hamiltonian consists of
projection operators designed to make the ground state orthogonal to the
projection onto the vectors

$\displaystyle\ket{F^{{c_{1},c_{2},c_{3}}}_{1,j}}=$ $\displaystyle
q^{-\frac{1}{2}}\ket{\uparrow_{j-1}^{c_{1}}\uparrow_{j}^{c_{2}}\downarrow_{j+1}^{c_{3}}}-q^{\frac{1}{2}}\ket{\uparrow_{j}^{c_{2}}\downarrow_{j+1}^{c_{3}}\uparrow_{j+2}^{c_{1}}},$
(1) $\displaystyle\ket{F^{{c_{1},c_{2},c_{3}}}_{2,j}}=$ $\displaystyle
q^{-\frac{1}{2}}\ket{\uparrow_{j-1}^{c_{1}}\downarrow_{j}^{c_{2}}\downarrow_{j+1}^{c_{3}}}-q^{\frac{1}{2}}\ket{\downarrow_{j-1}^{c_{3}}\uparrow_{j}^{c_{1}}\downarrow_{j+1}^{c_{2}}},$

in the spin or height sector to enforce weighted superposition of Dych paths
of different height between the $j$’th and $(j+1)$’th spin, and

$\displaystyle\ket{C_{j}}=$
$\displaystyle\ket{{\color[rgb]{1,0,0}\uparrow}_{j}{\color[rgb]{1,0,0}\downarrow}_{j+1}}-\ket{{\color[rgb]{0,0,1}\uparrow}_{j}{\color[rgb]{0,0,1}\downarrow}_{j+1}},$
(2)

in the color sector to enforce a balanced mixture of coloring of neighboring
up-down pairs. Together with the projectors that applies energy penalty on
color mismatching in the bulk and starting or ending the chain in the wrong
direction, the Hamiltonian of the colored Fredkin chain is

$\displaystyle H_{\mathrm{F}}=$
$\displaystyle\sum_{j=2}^{L-1}\sum_{a=1}^{2}\sum_{c_{1},c_{2},c_{3}=\mathrm{r,b}}\frac{1}{[2]_{q}}\ket{F^{{c_{1},c_{2},c_{3}}}_{a,j}}\bra{F^{{c_{1},c_{2},c_{3}}}_{a,j}}$
(3)
$\displaystyle+\sum_{j=1}^{L-1}\left(\frac{1}{2}\ket{C_{j}}\bra{C_{j}}+\ket{{\color[rgb]{1,0,0}\uparrow}_{j}{\color[rgb]{0,0,1}\downarrow}_{j+1}}\bra{{\color[rgb]{1,0,0}\uparrow}_{j}{\color[rgb]{0,0,1}\downarrow}_{j+1}}+\ket{{\color[rgb]{0,0,1}\uparrow}_{j}{\color[rgb]{1,0,0}\downarrow}_{j+1}}\bra{{\color[rgb]{0,0,1}\uparrow}_{j}{\color[rgb]{1,0,0}\downarrow}_{j+1}}\right)$
$\displaystyle+\sum_{c=\mathrm{r,b}}\left(\ket{\downarrow^{c}_{1}}\bra{\downarrow^{c}_{1}}+\ket{\uparrow^{c}_{L}}\bra{\uparrow^{c}_{L}}\right),$

where the q-deformed integer 2 is defined as $[2]_{q}\vcentcolon=q+q^{-1}$.
Its ground state can be written as

$|\mathrm{GS_{F}}\rangle=\frac{1}{\mathcal{N}_{\mathrm{F}}}\sum_{w\in\mathrm{colored}\
\mathrm{Dyck}\ \mathrm{walks}}q^{\frac{1}{2}\mathcal{A}(w)}|w\rangle,$ (4)

where $\mathcal{N}_{\mathrm{F}}$ is the normalization constant and
$\mathcal{A}(w)$ denotes the area underneath the Dyck path $w$.

The integer spin counterpart of colored Fredkin chain is called colored
Motzkin chain [16, 17], for its ground state is a superposition of Motzkin
paths/walks, which are Dyck paths diluted with spin-0’s or flat steps. The
bulk Hamiltonian projects onto the vectors orthogonal to the superposition in
the ground state

$\displaystyle\ket{M^{{c}}_{1,j}}=$ $\displaystyle
q^{-\frac{1}{2}}\ket{\uparrow_{j}^{c}0_{j+1}}-q^{\frac{1}{2}}\ket{0_{j}\uparrow_{j+1}^{c}},$
(5) $\displaystyle\ket{M^{{c}}_{2,j}}=$ $\displaystyle
q^{-\frac{1}{2}}\ket{0_{j}\downarrow_{j+1}^{c}}-q^{\frac{1}{2}}\ket{\downarrow_{j}^{c}0_{j+1}},$
$\displaystyle\ket{M^{{c}}_{3,j}}=$ $\displaystyle
q^{-\frac{1}{2}}\ket{\uparrow_{j}^{c}\downarrow_{j+1}^{c}}-q^{\frac{1}{2}}\ket{0_{j}0_{j+1}}.$

As color matching is automatically enforced, the colored Motzkin Hamiltonian
is just the sum over these bulk terms and the same boundary terms as the
Fredkin Hamiltonian.

$H_{\mathrm{M}}=\sum_{c=\mathrm{r,b}}\left(\sum_{j=2}^{L-1}\sum_{a=1}^{3}\frac{1}{[2]_{q}}\ket{M^{c}_{a,j}}\bra{M^{{c}}_{a,j}}+\ket{\downarrow^{c}_{1}}\bra{\downarrow^{c}_{1}}+\ket{\uparrow^{c}_{L}}\bra{\uparrow^{c}_{L}}\right).$
(6)

The ground state of the colored Motzkin chain is given as

$|\mathrm{GS_{M}}\rangle=\frac{1}{\mathcal{N}_{\mathrm{M}}}\sum_{w\in\mathrm{colored}\
\mathrm{Motzkin}\ \mathrm{walks}}q^{\mathcal{A}(w)}|w\rangle,$ (7)

where $\mathcal{N}_{\mathrm{M}}$ is the normalization constant and
$\mathcal{A}(w)$ denotes the area underneath the Motzkin path $w$.

When a cut in the middle separates the chain into two subsystems, both of the
ground states (4) and (7) can be Schmidt decomposed by the height $m$ of the
path in the middle at the cut, and the coloring of the $m$ spins in one of the
subsystems, which is to be matched with those in the other subsystem

$\ket{\mathrm{GS_{F(M)}}}=\sum_{m=0}^{L/2}\sqrt{p_{m}}\sum_{c_{1},...,c_{m}=\mathrm{r,b}}\ket{\mathcal{W}_{m}^{\vec{c}}}_{\mathrm{L}}\otimes\ket{\mathcal{W}_{m}^{\vec{c}}}_{\mathrm{R}},$
(8)

where $\ket{\mathcal{W}_{m}^{\vec{c}}}_{\mathrm{L,R}}$ are the normalized
superposition of all paths on the left and right subsystems that have height
$m$ in the middle point. Any such walk will have $m$ uncompensated up spins on
the left whose colors are exactly matched with $m$ down spins respectively on
the right subsystem. The vector $\vec{c}$ specifies the colors of the extra up
steps on the left. The Schmidt coefficient $p_{m}$ determines the EE of the
half chain

$S_{\frac{L}{2}}=-\sum_{m=0}^{L/2}2^{m}p_{m}\log p_{m}.$ (9)

Detailed analysis of the probability distribution of $m$ or its weighted
average among the Motzkin and Fredkin walks shows that the two model share the
same phase diagram 2, characterized by the scaling of half chain EE [16, 17,
19, 20, 21].

Figure 2: The common entanglement phase diagram of colored Fredkin and Motzkin
chains.

##  3 Six-vertex construction of coupled Fredkin chains

Figure 3: (a) A coloring of the minimal height configuration configuration of
a lattice of linear size 6. Three colors are used to manifest the matching
pattern. The shaded area with dashed line boundary marks the subsystem A. Then
numbers in the plaquettes indicate the height configuration in the dual
lattice. (In both (a) and (b), an additional green color is used to illustrate
the pairing of up-down spins at the same hight that match in color, even
though in the rest of the paper, the local Hilbert space of the model is
defined with only red and blue color.) (b) A coloring of the maximal height
configuration. The heights in the i’th square of plaquettes counting from
boundary inward alternates between $i$ and $i+1$.(c) The convention of
positive direction of the spins living on the array of horizontal chains
$S^{\mathrm{h}}$ and vertical chains $S^{\mathrm{v}}$ at vertex $(i,j)$, along
with the other 5 allowed vertex configurations.

The degrees of freedom of the model live on the edges or bonds between
vertices of a square lattice of linear size $L$. They can be decomposed into
two arrays of one-dimensional spin-$\frac{1}{2}$ chains, one horizontally and
one vertically aligned. Each of the edges in the array of horizontal (resp.
vertical) chains can have a spin $S^{\mathrm{h}}$ (resp. $S^{\mathrm{v}}$)
either up or down (resp. left or right) corresponding to $\pm\frac{1}{2}$.
These two sets of degrees of freedom are coupled to each other by the ice rule
[24, 25] in Fig. 3 (c), enforced by the bulk local Hamiltonian

$H_{0}=\sum_{x,y=2}^{L-1}(S_{x,y}^{\mathrm{h}}-S_{x,y+1}^{\mathrm{h}}-S_{x,y}^{\mathrm{v}}+S_{x+1,y}^{\mathrm{v}})^{2}.$
(10)

The global Hilbert space can be constrained to the subspace of six-vertex
configurations by making the coefficient of this term $V_{0}\gg 1$. The
boundary spins can be fixed by the Hamiltonian

$\begin{split}H_{\partial}=&\sum_{y=1}^{L}\left(S_{L,y}^{\mathrm{h}}-S_{1,y}^{\mathrm{h}}\right)+\sum_{x=2}^{L-1}(-1)^{x}\left(S_{x,1}^{\mathrm{h}}+S_{x,L-1}^{\mathrm{h}}\right)\\\
&+\sum_{x=1}^{L}\left(S_{x,L}^{\mathrm{v}}-S_{x,1}^{\mathrm{v}}\right)+\sum_{y=2}^{L-1}(-1)^{y}\left(S_{1,y}^{\mathrm{v}}+S_{L-1,y}^{\mathrm{v}}\right)+4L-4,\end{split}$
(11)

such that the boundary configurations in Fig. 3 has the minimal energy of 0,
and any other configurations will be penalized in proportion to the number of
local differences from them along the boundary.

The six-vertex condition allows a well-defined height function
$\phi_{x+\frac{1}{2},y+\frac{1}{2}}$ living on the dual lattice of plaquette
centers satisfying the rules according to the convention in Fig. 3 (c)

$\displaystyle\phi_{x+\frac{1}{2},y+\frac{1}{2}}-\phi_{x+\frac{1}{2},y-\frac{1}{2}}=2S_{x,y}^{\mathrm{v}},$
(12)
$\displaystyle\phi_{x+\frac{1}{2},y+\frac{1}{2}}-\phi_{x-\frac{1}{2},y+\frac{1}{2}}=2S_{x,y}^{\mathrm{h}},$
(13)

up to a global gauge transformation of shifting the heights by a constant. For
convenience we will fix the gauge so that the height $\phi_{0,0}=0$ at the
lower left corner of the lattice. The effect of boundary Hamiltonian amounts
to picking a (Dirichlet) boundary condition on the height for the ground state
wave-function, as shown in Fig. 3(a), Fig. 3(b) where the height function
alternates between $0$ and $1$ along the boundary.

To enrich the entanglement of the ground state, the local Hilbert space of
each spin is further enlarged to have either a red or blue color, along with a
local Hamiltonian $H_{C}$ between neighboring up-down spin pairs to match in
color

$\displaystyle H_{C}=$
$\displaystyle\sum_{x,y=1}^{L-1}\Bigg{(}\ket{{\color[rgb]{1,0,0}\boldsymbol{\uparrow}}_{x}{\color[rgb]{0,0,1}\boldsymbol{\downarrow}}_{x+1}}_{y}\bra{{\color[rgb]{1,0,0}\boldsymbol{\uparrow}}_{x}{\color[rgb]{0,0,1}\boldsymbol{\downarrow}}_{x+1}}+\ket{{\color[rgb]{0,0,1}\boldsymbol{\uparrow}}_{x}{\color[rgb]{1,0,0}\boldsymbol{\downarrow}}_{x+1}}_{y}\bra{{\color[rgb]{0,0,1}\boldsymbol{\uparrow}}_{x}{\color[rgb]{1,0,0}\boldsymbol{\downarrow}}_{x+1}}+\ket{\begin{subarray}{c}{\color[rgb]{0,0,1}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{1,0,0}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}_{x}\bra{\begin{subarray}{c}{\color[rgb]{0,0,1}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{1,0,0}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}+\ket{\begin{subarray}{c}{\color[rgb]{1,0,0}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{0,0,1}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}_{x}\bra{\begin{subarray}{c}{\color[rgb]{1,0,0}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{0,0,1}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}$ (14)
$\displaystyle+\frac{1}{[2]_{r}}\left(r^{-\frac{1}{2}}\ket{{\color[rgb]{1,0,0}\boldsymbol{\uparrow}}_{x}{\color[rgb]{1,0,0}\boldsymbol{\downarrow}}_{x+1}}_{y}-r^{\frac{1}{2}}\ket{{\color[rgb]{0,0,1}\boldsymbol{\uparrow}}_{x}{\color[rgb]{0,0,1}\boldsymbol{\downarrow}}_{x+1}}_{y}\right)\left(r^{-\frac{1}{2}}\bra{{\color[rgb]{1,0,0}\boldsymbol{\uparrow}}_{x}{\color[rgb]{1,0,0}\boldsymbol{\downarrow}}_{x+1}}_{y}-r^{\frac{1}{2}}\bra{{\color[rgb]{0,0,1}\boldsymbol{\uparrow}}_{x}{\color[rgb]{0,0,1}\boldsymbol{\downarrow}}_{x+1}}_{y}\right)$
$\displaystyle+\frac{1}{[2]_{r}}\left(r^{-\frac{1}{2}}\ket{\begin{subarray}{c}{\color[rgb]{1,0,0}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{1,0,0}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}_{x}-r^{\frac{1}{2}}\ket{\begin{subarray}{c}{\color[rgb]{0,0,1}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{0,0,1}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}_{x}\right)\left(r^{-\frac{1}{2}}\bra{\begin{subarray}{c}{\color[rgb]{1,0,0}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{1,0,0}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}_{x}-r^{\frac{1}{2}}\bra{\begin{subarray}{c}{\color[rgb]{0,0,1}\boldsymbol{\rightarrow}}\\\
{\color[rgb]{0,0,1}\boldsymbol{\leftarrow}}\end{subarray}_{y}^{y+1}}_{x}\right)\Bigg{)},$

where the up and down (resp. left and right) arrows are used to denote spin
$\rm\frac{1}{2}$ in the horizontal (resp. vertical) direction. The terms in
the first line energetically penalize adjacent up-down spin pairs mismatching
in color, so that spin configurations containing, say
$\ket{{\color[rgb]{1,0,0}\boldsymbol{\uparrow}}_{x}{\color[rgb]{0,0,1}\boldsymbol{\downarrow}}_{x+1}}_{y}^{\mathrm{h}}$
do not appear in the spin configuration of the zero energy ground state. Two
colored spin configurations that are not penalized by the mismatch penalty
terms in $H_{C}$ are examplified in Fig. 3 (a) and (b), where an additional
color green is employed to better illustrate the color matching between up-
down pairs of the same height. The terms in the next two lines enforce a
superposition of colorings of such adjacent up-down spin pairs when their
color is matched, tuned by the deformation parameter $0<r<1$. Indeed, whenever
a spin/color configuration
$\ket{{\color[rgb]{1,0,0}\boldsymbol{\uparrow}}_{x}{\color[rgb]{1,0,0}\boldsymbol{\downarrow}}_{x+1}}_{y}$
appears in the ground state, it must appear through the combination
$r^{\frac{1}{2}}\ket{{\color[rgb]{1,0,0}\boldsymbol{\uparrow}}_{x}{\color[rgb]{1,0,0}\boldsymbol{\downarrow}}_{x+1}}_{y}+r^{-\frac{1}{2}}\ket{{\color[rgb]{0,0,1}\boldsymbol{\uparrow}}_{x}{\color[rgb]{0,0,1}\boldsymbol{\downarrow}}_{x+1}}_{y}$
in order to a avoid an energy penalty from these projection operators. In this
way, these projection operators are necessary to provide color mixing and
ergodicity within the subspace of product states annihilated by the terms in
the first line.

The deformation parameter $r$ plays the role of an external color field, such
that when $r=1$, the ground state will have a uniform superposition of
different coloring, while when $r>1$, the configurations with more red colored
spins will be favored.

Since the color Hamiltonian only acts on up-down and left-right pairs, for it
to affect all the spins in the system, there must be a net surplus of up
(resp. left) spins in any sub-chain counting from left (resp. bottom). In
other words, the height function in the dual lattice must stay non-negative
and the spins form Dyck paths along the chains in both directions. This can be
enforced by the correlated swapping Hamiltonian

$H_{S}=\sum_{x,y=2}^{L-1}\sum_{f_{\mathrm{h}},f_{\mathrm{v}}=\pm}\sum_{c_{1},...,c_{6}=\mathrm{r,b}}\frac{1}{[2]_{q}}\ket{\pi_{x,y,f_{\mathrm{h}},f_{\mathrm{v}}}^{c_{1},...,c_{6}}}\bra{\pi_{x,y,f_{\mathrm{h}},f_{\mathrm{v}}}^{c_{1},...,c_{6}}}$
(15)

where

$\displaystyle\ket{\pi_{x,y,-,+}^{c_{1},...,c_{6}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 35.0pt 0.0pt
0.0pt,scale={0.6}]{ltfc.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
35.0pt 0.0pt 0.0pt,scale={0.6}]{ltec.pdf}},$ (16)
$\displaystyle\ket{\pi_{x,y,+,+}^{c_{1},...,c_{6}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 35.0pt 0.0pt
0.0pt,scale={0.6}]{rtfc.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
35.0pt 0.0pt 0.0pt,scale={0.6}]{rtec.pdf}},$
$\displaystyle\ket{\pi_{x,y,-,-}^{c_{1},...,c_{6}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 35.0pt 0.0pt
0.0pt,scale={0.6}]{lbfc.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
35.0pt 0.0pt 0.0pt,scale={0.6}]{lbec.pdf}},$
$\displaystyle\ket{\pi_{x,y,+,-}^{c_{1},...,c_{6}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 35.0pt 0.0pt
0.0pt,scale={0.6}]{rbfc.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
35.0pt 0.0pt 0.0pt,scale={0.6}]{rbec.pdf}}.$

The total Hamiltonian

$H_{cF}=H_{0}+H_{\partial}+H_{C}+H_{S}$ (17)

is a frustration-free sum of projection operators, meaning its zero energy
ground state is the simultaneous lowest energy eigenstate of each term. Since
each term in the Hamiltonian requires a superposition of locally different
height and coloring in a particular way, the ground state is therefore a
weighted superposition of bicolored six-vertex configurations with alternating
heights between 0 and 1 along the boundary, and non-negative heights in the
bulk

$\ket{\mathrm{GS}}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{\phi(\partial
P)=\frac{1}{2}}\prod_{x,y=1}^{L}\theta(\phi_{x,y})\sideset{}{{}^{\prime}}{\sum}_{C}r^{\frac{M(C)}{2}}q^{\frac{\mathcal{V}(P)}{2}}\ket{P^{C}},$
(18)

where for simplified notation, the height function of each spin is taken to be
the average between the heights of its two adjacent plaquettes, the first sum
is over all six-vertex configurations $P$ with boundary height $\frac{1}{2}$,
the second primed sum is over coloring patterns with spins in the same chain
on the same height matching. $\theta$ is the Heaviside step function,
indicating the sum is over configurations with non-negative height in the
bulk. The volume of a configuration is defined as
$\mathcal{V}(P)=\sum_{x,y=1}^{L-1}\phi_{x+\frac{1}{2},y+\frac{1}{2}}$, $M(C)$
is the “warmness” magnetization of coloring $C$, defined as the difference
between the number of pairs of red and blue spins, and $\mathcal{N}$ is the
normalization constant that depends only on $q$ and $r$. The uniqueness of the
ground state is guaranteed by the ergodicity of the Hamiltonian (15), which is
proven in the appendix.

##  4 Scaling of entanglement entropy

The model has an apparent $D_{4}$ lattice symmetry, so a cut across the middle
along either the horizontal or vertical direction gives the same bipartite
entanglement entropy between subsystems. Unlike a quasi-2D model of trivial
stacking an array of Fredkin or Motzkin chains, the ground state EE scaling
behavior of this coupled 2D model is the same for cuts in any direction.
Without loss of generality, we choose a vertical cut as shown in Fig. 3. Just
like the one-dimensional model described in Sec. 2, the entanglement comes
from the extra up spins in the left subsystem, or equivalently the surplus of
down spins in the right subsystem. Here, when performing the Schmidt
decomposition on the ground state (18), we need to keep track of the heights
of each chain in the array along the cut in the middle, denoted by the vector
$\vec{\phi}_{\frac{L+1}{2}}$, as well as the vectors $\vec{c}_{y}$’s denoting
the colors of the spins in each row to be matched between the two subsystem.
The number of components of each vector $\vec{c}_{y}$ is given by the value of
the component of height vector $\phi_{\frac{L+1}{2},y+\frac{1}{2}}$. Notice
that $\vec{\phi}_{\frac{L+1}{2}}$ itself is constrained to be a Dyck path in
the vertical direction, the bulk components of which are always larger than or
equal to the first and last components
$\phi_{\frac{L+1}{2},\frac{1}{2}}\equiv\phi_{\frac{L+1}{2},L+\frac{1}{2}}=0,1$.
Therefore the Schmidt decompostion can be written as

$\ket{\mathrm{GS}}=\\!\\!\\!\sum_{\vec{\phi}_{\frac{L+1}{2}}\in\mathrm{Dyck\
paths}}\\!\\!\\!\sqrt{\frac{\mathcal{M}^{2}_{\vec{\phi}}}{\mathcal{N}}}\\!\\!\sum_{\vec{c}_{1}\in\\{\mathrm{r,b}\\}^{\phi_{\frac{L+1}{2},3/2}}}\\!\\!\\!\cdots\\!\\!\\!\\!\sum_{\vec{c}_{L}\in\\{\mathrm{r,b}\\}^{\phi_{\frac{L+1}{2},L+1/2}}}\left(\prod_{y=1}^{L}r^{\frac{m(\vec{c}_{y})}{2}}\right)\ket{P^{\vec{c}_{1},...,\vec{c}_{L}}_{\vec{\phi}_{\frac{L+1}{2}}}}_{\mathrm{L}}\\!\\!\\!\otimes\ket{P^{\vec{c}_{1},...,\vec{c}_{L}}_{\vec{\phi}_{\frac{L+1}{2}}}}_{\mathrm{R}},$
(19)

where the “warmness” magnetization $m(\vec{c}_{y})$ of the unmatched colors of
the half system in the y’th chain is the difference between the number of red
and blue spin pairs among them, and

$\ket{P^{\vec{c}_{1},...,\vec{c}_{L}}_{\vec{\phi}_{\frac{L+1}{2}}}}_{\mathrm{L(R)}}=\frac{1}{\sqrt{\mathcal{M}_{\vec{\phi}_{\frac{L+1}{2}}}}}\sideset{}{{}^{\prime}}{\sum}r^{\frac{m_{\mathrm{L(R)}}}{2}}q^{\frac{\mathcal{V}_{\mathrm{L(R)}}(P^{C})}{2}}\ket{P^{C}}_{\mathrm{L(R)}}$
(20)

are normalized wave functions of the left (resp. right) subsystems, the primed
sum is a shorthand notation for summing over six-vertex configurations with
non-negative height in the bulk plaquettes and in particular of height
specified by $\vec{\phi}_{\frac{L+1}{2}}$ on the middle boundary.
$m_{\mathrm{L(R)}}$ is the redness magnetization of the pairs with color
matched within the subsystem, which takes value between 0 and
$(L^{2}-L-\frac{A(\vec{\phi}_{\frac{L+1}{2}})+3L}{2}+1)/2$, with
$A(\vec{\phi}_{\frac{L+1}{2}})=\sum_{y=1}^{L-1}\phi_{\frac{L+1}{2},y+1/2}$
being the cross-sectional area of the stepped surface outlined by the height
function, as depicted in Fig. 4.

Figure 4: Cross-sectional view of the stepped surface outlined by the height
function along the middle cut. Area $\mathcal{A}$ is defined as the sum of
heights at each step counting from the minimal heights in Fig. 3 (a).

The normalization constants are given by

$\mathcal{M}^{2}_{\vec{\phi}_{\frac{L+1}{2}}}=[2]_{r}^{L^{2}-\frac{5L}{2}-\frac{\mathcal{A}(\vec{\phi}_{\frac{L+1}{2}})}{2}+1}\sum_{\vec{\phi}_{\frac{L+1}{2}}(T|_{\mathrm{cut}})=\vec{\phi}_{\frac{L+1}{2}}}q^{\mathcal{V}(P)},$
(21)

and

$\displaystyle\mathcal{N}=$
$\displaystyle\sum_{\vec{\phi}}[2]_{r}^{\frac{A(\vec{\phi}_{\frac{L+1}{2}})+3L}{2}-1}\mathcal{M}^{2}_{\vec{\phi}_{\frac{L+1}{2}}}$
(22) $\displaystyle\equiv$ $\displaystyle[2]_{r}^{L^{2}-L}\sum_{\phi(\partial
P)=\frac{1}{2}}q^{\mathcal{V}(P)}.$ (23)

The Schmidt coefficients are given by the probability of height configuration
$\vec{\phi}$ with coloring $\\{\vec{c}_{1},...,\vec{c}_{L}\\}$ of the cross-
section between subsystems

$p(\vec{\phi}_{\frac{L+1}{2}},\\{\vec{c}_{1},...,\vec{c}_{L}\\})=\frac{\left(\prod_{y=1}^{L}r^{\frac{m(\vec{c}_{y})}{2}}\right)^{2}\mathcal{M}^{2}_{\vec{\phi}_{\frac{L+1}{2}}}}{\mathcal{N}}=p(\\{\vec{c}_{1},...,\vec{c}_{L}\\}\mid\vec{\phi}_{\frac{L+1}{2}})p(\vec{\phi}_{\frac{L+1}{2}}),$
(24)

with

$p(\\{\vec{c}_{1},...,\vec{c}_{L}\\}\mid\vec{\phi}_{\frac{L+1}{2}})=[2]_{r}^{-\frac{\mathcal{A}(\vec{\phi}_{\frac{L+1}{2}})}{2}-\frac{3L}{2}+1}\prod_{y=1}^{L}r^{m(\vec{c}_{y})},$
(25)

can be factorized as a product of probability of having a particular coloring
$\\{\vec{c}_{1},...,\vec{c}_{L}\\}$ of the unmatched pairs within the
subsystems, conditioned on having a Dyck path $\vec{\phi}_{\frac{L+1}{2}}$
along the cut, and the marginal probability
$p(\vec{\phi}_{\frac{L+1}{2}})\equiv\sum_{\vec{c}_{1}}\cdots\sum_{\vec{c}_{L}}p(\vec{\phi_{\frac{L+1}{2}}},\\{\vec{c}_{1},...,\vec{c}_{L}\\})$
of finding such a cross section among uncolored random height configurations.
The entanglement entropy decomposed into a piece given in terms of average
cross-sectional area of a random height configuration, and another subleading
contribution from the fluctuation of the random surface

$\displaystyle S_{\frac{L}{2}\times L}(q,r)=$
$\displaystyle-\sum_{\vec{\phi}_{\frac{L+1}{2}}}\sum_{\vec{c}_{1}}\cdots\sum_{\vec{c}_{L}}p(\vec{\phi}_{\frac{L+1}{2}},\\{\vec{c}_{1},...,\vec{c}_{L}\\})\log
p(\vec{\phi},\\{\vec{c}_{1},...,\vec{c}_{L}\\})$ (26) $\displaystyle=$
$\displaystyle\sum_{\vec{\phi}_{\frac{L+1}{2}}}p(\vec{\phi}_{\frac{L+1}{2}})S^{c}_{\frac{L}{2}\times
L}(\vec{\phi}_{\frac{L+1}{2}},r)+S_{\frac{L}{2}\times L}^{\phi}(q)$
$\displaystyle=$
$\displaystyle\frac{C_{r}}{2}\left(\langle\mathcal{A}\rangle+3L-2\right)+S^{\phi}_{L}(q),$

where $S^{\phi}_{\frac{L}{2}\times
L}(q)=-\sum_{\vec{\phi}_{\frac{L+1}{2}}}p(\vec{\phi}_{\frac{L+1}{2}})\log
p(\vec{\phi}_{\frac{L+1}{2}})$ and in third line we have used

$\displaystyle S^{c}_{\frac{L}{2}\times L}(\vec{\phi}_{\frac{L+1}{2}},r)=$
$\displaystyle-\sum_{\\{\vec{c}_{1},...,\vec{c}_{L}\\}}p(\\{\vec{c}_{1},...,\vec{c}_{L}\\}\mid\vec{\phi}_{\frac{L+1}{2}})\log
p(\\{\vec{c}_{1},...,\vec{c}_{L}\\}\mid\vec{\phi}_{\frac{L+1}{2}})$ (27)
$\displaystyle=$ $\displaystyle-(\log
r)[2]_{r}^{-\frac{\mathcal{A}(\vec{\phi}_{\frac{L+1}{2}})}{2}-\frac{3L}{2}+1}\sum_{\\{\vec{c}_{1},...,\vec{c}_{L}\\}}\left(\sum_{y=1}^{L}m(\vec{c}_{y})\right)\prod_{y=1}^{L}r^{m(\vec{c}_{y})}$
$\displaystyle+\log[2]_{r}\left(\frac{\mathcal{A}(\vec{\phi}_{\frac{L+1}{2}})}{2}+\frac{3L}{2}-1\right)$
$\displaystyle=$ $\displaystyle
C_{r}\left(\frac{\mathcal{A}(\vec{\phi}_{\frac{L+1}{2}})}{2}+\frac{3L}{2}-1\right).$

The sum in the second line can also be written as a sum over the total number
of pairs of red spin among $\\{\vec{c}_{1},...,\vec{c}_{L}\\}$,
$i\equiv\sum_{y=1}^{L}\frac{m(\vec{c}_{y})+\phi_{\frac{L+1}{2},y+\frac{1}{2}}}{2}$
which ranges from $0$ to
$N\equiv\frac{\mathcal{A}(\vec{\phi}_{\frac{L+1}{2}})}{2}+\frac{3L}{2}-1$:

$\displaystyle\sum_{\\{\vec{c}_{1},...,\vec{c}_{L}\\}}\left(\sum_{y=1}^{L}m(\vec{c}_{y})\right)\prod_{y=1}^{L}r^{m(\vec{c}_{y})}=$
$\displaystyle\sum_{i=0}^{N}\binom{N}{i}(2i-N)r^{2i-N}$ (28) $\displaystyle=$
$\displaystyle\frac{r-r^{-1}}{[2]_{r}}[2]_{r}^{N}N,$

which gives the coefficient

$C_{r}=-\frac{r-r^{-1}}{[2]_{r}}\log r+\log[2]_{r}.$ (29)

This kind of decomposition of entanglement entropy as a result of enlarging
the local Hilbert space has also been observed recently in the Bethe Ansatz
integrable excited states of a non-integrable one-dimensional multicomponent
spin chain [26], which emerges from certain phases of a quasi-2D spin ladder
[27].

For any finite $r$, $C_{r}$ is a finite constant independent of $L$, so the
problem is reduced to finding the scaling of the average area
$\langle\mathcal{A}\rangle$. That can be done in a field theoretic fashion, as
was previously used to study the dynamics of the one-dimensional Motzkin and
Fredkin chains [28, 29]. A continuous field of the height configuration can be
defined as a piece-wise linear function $\phi(x,y)$, which takes the value of
$\phi_{x+1/2,y+1/2}$ on the dual lattice. It is well known that the “entropy”
of random surface is captured by a surface tension $\sigma(\nabla\phi(x,y))$
as a function of the height gradient alone [30, 31, 32, 33]. Also taking into
account the “energy” contribution from volume weighting, we get the partition
function

$Z=\int\mathcal{D}\phi(x,y)e^{\iint dxdy(-\sigma(\nabla\phi(x,y))+(\log
q)\phi(x,y))},$ (30)

where $\mathcal{D}\phi(x,y)$ is a continuous version of

$\prod_{x,y}\int_{0}^{+\infty}d\phi_{x+1/2,y+1/2}\equiv\prod_{v}\int_{-\infty}^{+\infty}dh_{v}\theta(h_{v}),$
(31)

and where $\phi$ obeys a Lipschitz condition
$|\partial_{x^{\prime}}\phi,\partial_{y^{\prime}}\phi|\leq 1$, and $\theta$ is
the Heaviside step function.

The linear contribution in $\phi$ is dominant when $q>1$. To see this
explicitly, we substitute

$x=Lx^{\prime},\quad y=Ly^{\prime},$ (32)

which makes

$\nabla=L^{-1}\nabla^{\prime},\quad dx=Ldx^{\prime},\quad dy=Ldy^{\prime}.$
(33)

The free energy associated to a height configuration then becomes

$F[\phi]=L^{2}\iint_{0}^{1}dx^{\prime}dy^{\prime}\Big{(}\sigma\big{(}\frac{\nabla^{\prime}\phi(x^{\prime},y^{\prime})}{L}\big{)}-(\log
q)\phi(x^{\prime},y^{\prime})\Big{)},$ (34)

where $\phi(x^{\prime},y^{\prime})=\phi(x,y)$ now satisfy the Lipschitz
property of $|\partial_{x^{\prime}}\phi,\partial_{y^{\prime}}\phi|\leq L$
instead. The surface tension term counts the entropy of height configurations
associated with height variations in a small region with height differences
$\partial_{x^{\prime}}\phi,\partial_{y^{\prime}}\phi$ on the boundary of the
region, and is thus trivially bounded by the entropy density of ice. Thus, in
the thermodynamic limit, the surface tension term becomes irrelevant compared
to the linear term $(\log q)\phi(x^{\prime},y^{\prime})$ when $q\neq 1$.
Therefore $F[\phi]$ is minimized by the Lipschitz property for the $q>1$ case,
where minimization of $F$ is achieved with maximal gradient and maximal
volume; and by the positivity for the $q<1$ case, where $F$ is minimized
taking $\phi(x,y)=0$. Therefore, for $q$ larger and smaller than 1
respectively, we have $\langle\mathcal{A}\rangle=O(L^{2})$ and $O(L)$. (26)
then says $S_{L}(q,r)$ goes through a phase transition at $q=1$ from volume
scaling to area law.

At the critical point, the height field becomes a massless field conditioned
on staying positive. Given that the surface tension is a strictly convex even
function of the height variable [30, 31, 32, 34, 35, 36, 37], the average
height was rigorously shown to have the scaling $O(\log L)$, as a result of
being repelled by the hard-wall at zero height [38]. This gives the same EE
scaling of $O(L\log L)$ as in the recent quantum lozenge tiling model [23],
despite the height field of uniform weighted six-vertex model being
interacting and not described by a Gaussian free field. This entanglement
phase transition can be summarized in the phase diagram in Fig. 5.

Figure 5: The entanglement phase diagram with two distinct EE scaling
separated by the critical line at $q=1$.

The stark phase transition for any $\epsilon=q-1$ is a consequence of the
discontinuity of the partition function $Z$ when the thermodynamic/scaling
limit is taken. To obtain an intermediate scaling between $L\log L$ and
$L^{2}$, one can consider a varying $q=e^{\lambda L^{-\alpha}}$ that
approaches 1 as $L\to\infty$. Such scaling limits are of interest random
surface models, as they admit existence of non-trivial limit shapes [39]. For
$\alpha\in(1,2)$, simple scaling argument gives an EE scaling of
$L^{3-\alpha}$ with a $\lambda$ dependent coefficient. Whereas for $\alpha\geq
2$, it gives the $L\log L$ scaling, and for $\alpha\leq 1$, it gives the
$L^{2}$ scaling. Interestingly, one can think of this intermediate entropy
scaling as the scaling of entropy associated with a square neighbourhood of
size $L^{\prime}$ attached to the corner of a larger lattice where the
deformation parameter is inhomogenous, decaying as a function $e^{\lambda
d^{-\alpha}}$ of the distance $d$ from the corner of the lattice to the
center.

##  5 Scaling of spectral gap

Following the strategy in the proof of the gaplessness of the highly entangled
phase of the one-dimensional models [40, 21], we construct a variational wave
function that has both a small overlap with the ground state (18) and an
exponentially small expectation value of the Hamiltonian (17). Hence it
inevitably implies that the spectral gap of the $q>1$ phase is exponentially
small and hence gapless in the thermodynamic limit.

Figure 6: Two height configurations of a $14\times 14$ system that differ by
the height of one plaquette $(6+\frac{1}{2},7+\frac{1}{2})$. (a) belongs to
the set $\mathcal{E}$ but (b) doesn’t.

We start by defining a smaller set of six-vertex configurations than those
with non-negative height in the bulk, which will appear in the superposition
of the trial excited state. The criteria of this set $\mathcal{E}$ are that:
first, the lowest height in the bulk of a configuration must be either 0 or 1;
and second, the longest distance between the lowest height in the bulk (be it
0 or 1) and any of the four sides of the boundary is larger than
$\frac{L}{2}$. An example of a configuration in this set, incidentally also
one of the four such ones with lowest total volume, is shown in Fig. 6 (a).
Whereas the configuration with largest volume among those not belonging to
this set is given in Fig. 6 (b).

The trial excited state is defined by changing the color of the spin on the
endpoint inside the bulk along the said longest distance to the boundary from
the plaquette with lowest height

$\ket{\mathrm{ex}}=\frac{1}{\sqrt{\mathcal{N}^{\prime}}}\sum_{P\in\mathcal{E}}\sideset{}{{}^{\prime\prime}}{\sum}_{C}r^{\frac{M(C)}{2}}q^{\frac{\mathcal{V}(P)}{2}}\ket{P^{C}},$
(35)

where $\sqrt{N^{\prime}}$ is the new normalization constant, and compared to
(18), the double primed sum over colored refers the matching of all the other
pairs of spins at the same height in color except the one at the endpoint of
the longest nonzero height streak that is now forced to mismatch in color. Due
to this mismatching, the excited state must be orthogonal to the ground state

$\bra{\mathrm{ex}}\mathrm{GS}\rangle=0.$ (36)

Figure 7: The volume of the dent between the minimal volume configuration in
Fig. 6 and the maximal volume configuration scales as $\frac{L^{3}}{8}$ with
the system size.

Furthermore, since all the configurations in $\mathcal{E}$ are of larger
volume than the one with the longest streak containing the color mismatch, the
two mismatched spins never appear as neighbors in the superposition, we have

$H_{0}\ket{\mathrm{ex}}=0,\quad
H_{\partial}\ket{\mathrm{ex}}=0,\quad\mathrm{and}\quad
H_{C}\ket{\mathrm{ex}}=0.$ (37)

So the non-vanishing contribution to the energy expectation can only come from
$H_{S}$, precisely from the term acting on the plaquette, decreasing the
height on which would result in a configuration outside the set $\mathcal{E}$.
In the case of the configuration $P_{a}$ in Fig. 6 (a), the terms involved in
$H_{S}$ will be

$h_{6,7}=\sum_{c_{1},...,c_{6}=\mathrm{r,b}}\frac{1}{[2]_{q}}\left(\ket{\pi_{6,7,+,+}^{c_{1},...,c_{6}}}\bra{\pi_{6,7,+,+}^{c_{1},...,c_{6}}}+\ket{\pi_{6,7,+,-}^{c_{1},...,c_{6}}}\bra{\pi_{6,7,+,-}^{c_{1},...,c_{6}}}\right).$
(38)

Together they contribute $\bra{P_{a}}h_{6,7}\ket{P_{a}}=\frac{2}{1+q^{2}}$,
for each particular color configuration. The number of such height
configurations that can be brought out of set $\mathcal{E}$ can be very
roughly upper bounded by the total number of spin and color configurations
$(4\max\\{(1+r^{2}),(1+r^{-2})\\})^{L^{2}}$. However, all of them has a volume
$\frac{L^{3}}{8}$ smaller than the maximum configuration, as depicted in Fig.
7. Lower bounding the normalization $\mathcal{N}^{\prime}$ by the weight of
the largest volume configuration, we have the upper bound on spectral gap

$\bra{\mathrm{ex}}H_{cF}\ket{\mathrm{ex}}<\frac{2(4\max\\{(1+r^{2}),(1+r^{-2})\\})^{L^{2}}}{1+q^{2}}q^{-\frac{L^{3}}{8}},$
(39)

which is gapless in the thermodynamic limit for the $q>1$ phase.

##  6 Nineteen-vertex construction of coupled Motzkin chains

Figure 8: (a) A random coupled Motzkin lattice configuration with Aztec
diamond boundary. (b) The 19 vertices in the constrained Hilbert space
satisfying equal number of in and out arrows. Figure 9: One coloring
configuration for each of the 9 coupled Motzkin moves defined on plaquette
$(x,y)$ that generate the ergodic dynamics in the Hilbert space of random
surfaces of height function conditioned to be non-negative.

Building on the previous sections, we introduce a 2D generalization to the
Motzkin chain, where each spin takes value of either $\pm 1$ or $0$. This can
be mapped to solid edges with arrows and dashed edges without (which
correspond to spin $0$ in the one-dimensional chain) respectively, giving 19
vertex configurations in Fig. 8 (b) a full loop around each of which the net
height change is $0$, so that the height change is well-defined counting from
two different paths from one plaquette to another. Nineteen-vertex model is a
generalization to six-vertex model and is well-studied in the context of
classical statistical mechanics [41, 42, 43, 44, 45, 46, 47]. Note that
classical nineteen-vertex models are mapped to quantum spin-$1$ chains, by
transfer matrix method, which was studied in the context of integrability [48,
49]. However, in this section, we construct a different (2D) quantum
Hamiltonian that is frustration free, which enforce the ground state to be a
weighted superposition between pairs of locally different configurations in
Fig. 9. The boundary spins in the lattice shown in Fig. 8 (a) is enforced by
boundary Hamiltonians that penalizes $-1$ spins on the left and bottom side
and $+1$ on the right and top side. The bulk Hamiltonian can be defined as

$H_{cM}=\sum_{p\in\mathrm{bulk}}\sum_{h,v=1}^{3}\sum_{c_{1},c_{2}=\mathrm{r,b}}\frac{1}{[2]_{q}}\ket{\pi_{p,h,v}^{c_{1},c_{2}}}\bra{\pi_{p,h,v}^{c_{1},c_{2}}},$
(40)

where

$\displaystyle\ket{\pi_{p,1,1}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 27.0pt 0.0pt
0.0pt,scale={0.5}]{hm11f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
10.0pt 0.0pt 0.0pt,scale={0.5}]{hm11e.pdf}},$ (41)
$\displaystyle\ket{\pi_{p,1,2}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 22.0pt 0.0pt
0.0pt,scale={0.5}]{hm12f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
18.0pt 0.0pt 0.0pt,scale={0.5}]{hm12e.pdf}},$
$\displaystyle\ket{\pi_{p,1,3}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 22.0pt 0.0pt
0.0pt,scale={0.5}]{hm13f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
20.0pt 0.0pt 0.0pt,scale={0.5}]{hm13e.pdf}},$
$\displaystyle\ket{\pi_{p,2,1}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 27.0pt 0.0pt
0.0pt,scale={0.5}]{hm21f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
10.0pt 0.0pt 0.0pt,scale={0.5}]{hm21e.pdf}}$
$\displaystyle\ket{\pi_{p,2,2}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 22.0pt 0.0pt
0.0pt,scale={0.5}]{hm22f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
18.0pt 0.0pt 0.0pt,scale={0.5}]{hm22e.pdf}},$
$\displaystyle\ket{\pi_{p,2,3}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 12.0pt 0.0pt
0.0pt,scale={0.5}]{hm23f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
20.0pt 0.0pt 0.0pt,scale={0.5}]{hm23e.pdf}},$
$\displaystyle\ket{\pi_{p,3,1}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 25.0pt 0.0pt
0.0pt,scale={0.5}]{hm31f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
15.0pt 0.0pt 0.0pt,scale={0.5}]{hm31e.pdf}},$
$\displaystyle\ket{\pi_{p,3,2}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 22.0pt 0.0pt
0.0pt,scale={0.5}]{hm32f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
18.0pt 0.0pt 0.0pt,scale={0.5}]{hm32e.pdf}},$
$\displaystyle\ket{\pi_{p,3,3}^{c_{1},c_{2}}}$
$\displaystyle=q^{-\frac{1}{2}}\ket{\includegraphics[trim=0.0pt 22.0pt 0.0pt
0.0pt,scale={0.5}]{hm33f.pdf}}-q^{\frac{1}{2}}\ket{\includegraphics[trim=0.0pt
20.0pt 0.0pt 0.0pt,scale={0.5}]{hm33e.pdf}}.$

$H_{0}$ and $H_{\partial}$ are defined in exactly the same way as in the model
of coupled Fredkin chains enforcing local constraints on the Hilbert space in
Fig. 8 (b) and boundary configurations as in Fig. 8 (a), but the Hamiltonian
acting on the color sector is already encoded in $H_{cM}$.

##  7 Conclusions

In this paper we have shown how quantum height models may be enhanced to give
a range of exotic entropy scalings. Our models can be viewed as coupled
Fredkin and Motzkin chains. They provide another example of a local
Hamiltonian with volume scaling of EE. While the height degree of freedom can
be described via an appropriate field theory, the addition of the color
degrees of freedom within such a description is an interesting open question.
Moreover the field theory description only holds for the ground state, while
structure of excited states is a subject for additional work.

The equal time correlation functions of the ground state of our model are
given by the correlation functions of classical six-vertex model subject to
the constraint of positive height. Even in the absence of such a constraint,
the analytical result of its two-point correlation functions are only computed
for certain boundary conditions such as the domain wall boundary [50, 51, 52].
But it’s possible to compute them numerically using Markov Chain Monte Carlo
method [53]. Adding the non-negative height constraint would pose a challenge
to the application of worm or loop-building algorithms, as maintaining the
non-negativity would require checking a larger neighborhood as the loops get
longer in each update. Another interesting next step in that direction will be
the construction of a tensor network characterization for the state, as was
done for in the 1D case [54, 55]. Finally, our model in the absence of an
internal color degree of freedom is of interest as it promises anomalous slow
dynamics and fragmentation analogous to the classical and quantum Fredkin
chains in one dimension [56, 57, 58].

## Acknowledgements

ZZ thanks Filippo Colomo, Kari Eloranta, Christophe Garban, Hosho Kastura,
Yuan Miao, Henrik Røising and Benjamin Walter for fruitful discussions. ZZ
acknowledges the kind hospitality of the Galileo Galilei Institute for
Theoretical Physics during the workshops “Randomness, Integrability and
Universality” and “Machine Learning at GGI". We gratefully acknowledge support
from the Simons Center for Geometry and Physics, Stony Brook University at
which some of the research for this paper was performed.

### Funding information

The work of IK was supported in part by the NSF grant DMR-1918207.

##  Appendix A Ergodicity of the Hamiltonian and uniqueness of ground state

We now show that when the Hamiltonian $H_{S}$ acts on a properly colored
height configuration it generates another such configuration, and that
moreover by actions of $H_{S}$ we can get from any such configuration to any
other. Thus the set of non-negative weighted height configurations with
Dirichlet boundaries is closed under the operation, with the weighted
superposition of states a unique ground state. In complete analogy with the 1D
Motzkin and Fredkin chains, starting from a state which violates non-
positivity in the bulk, by applying the projectors we create a superposition
that will carry the negative region back to the end of the sample to get
penalized by the boundary terms. Just as in the Fredkin chain case, in a non-
negative height superposition involving a color violation, by reducing the
height of unmatched color pairs may be pushed closer until the violation can
be detected by local terms.

Figure 10: A snippet of plateaux contour of the maximal height $\phi_{M}$,
with the surrounding lower height plaquettes coded by sequentially lighter
shades. The heights of the plaquettes numbered \raisebox{-.9pt} {1}⃝ located
at the corners of the contour are ready to be lowered, while those numbered
\raisebox{-.9pt} {3}⃝ are not. The plaquette \raisebox{-.9pt} {2}⃝ is an
accidental mobile plaquette in this contour configuration despite not lying on
the corner.

Let us now check that we can get to the lowest height configuration from any
positive height configuration. Given any six-vertex configuration, there must
be a plaquette of maximal height $\phi_{M}$, which may not necessarily be
unique. Their nearest neighbor have height $\phi_{M}-1$, but the next-nearest
neighbors could either have height $\phi_{M}$ or $\phi_{M}-2$. In the former
case, we say the maximal height forms a plateau, while in the latter case, it
either lies on the boundary of a plateau, or is isolated. We note that a local
maximal height plaquette will have color matched pairs of edges, because of
the color rule Eq. (14), therefore it can be removable by one the four moves
in (16). Similarly, plaquettes that are on the the boundaries of plateaux are
removable if they are at the corner of boundaries (along a straight line of
boundary, both sides in the direction of the boundary are not in the right
configuration to allow one of the correlated swapping moves), since the next-
nearest neighbors are both of the same height. Thus, given any boundary of a
plateau, we can always reduce the volume of a surface by first removing the
height cubes on the (convex) corners of plateaux boundaries, after which new
corners will appear, so that the procedure keeps going. The only scenario such
a procedure terminates is when the boundary forms a straight line with the
plateau extending to the boundary of the lattice. In that case, both sides of
the straight line have the same constant height as the boundary, meaning we
have arrived at the lowest height configuration.

## References

  * [1] N. Laflorencie,  _Quantum entanglement in condensed matter systems_ ,  Physics Reports 646, 1 (2016).
  * [2] J. Eisert, M. Cramer and M. B. Plenio,  _Colloquium: Area laws for the entanglement entropy_ ,  Rev. Mod. Phys. 82, 277 (2010),  10.1103/RevModPhys.82.277.
  * [3] A. Kitaev and J. Preskill,  _Topological entanglement entropy_ ,  Phys. Rev. Lett. 96, 110404 (2006),  10.1103/PhysRevLett.96.110404.
  * [4] T. Nishioka,  _Entanglement entropy: Holography and renormalization group_ ,  Rev. Mod. Phys. 90, 035007 (2018),  10.1103/RevModPhys.90.035007.
  * [5] D. N. Page,  _Average entropy of a subsystem_ ,  Phys. Rev. Lett. 71, 1291 (1993),  10.1103/PhysRevLett.71.1291.
  * [6] M. B. Hastings,  _An area law for one-dimensional quantum systems_ ,  Journal of Statistical Mechanics: Theory and Experiment 2007(08), P08024 (2007).
  * [7] A. Anshu, I. Arad and D. Gosset,  _An area law for 2D frustration-free spin systems_ ,  arXiv e-prints arXiv:2103.02492 (2021),  2103.02492.
  * [8] J. Eisert, M. Cramer and M. B. Plenio,  _Colloquium: Area laws for the entanglement entropy_ ,  Rev. Mod. Phys. 82, 277 (2010),  10.1103/RevModPhys.82.277.
  * [9] P. Calabrese and J. Cardy,  _Entanglement entropy and conformal field theory_ ,  Journal of Physics A: Mathematical and Theoretical 42(50), 504005 (2009),  10.1088/1751-8113/42/50/504005.
  * [10] D. Gioev and I. Klich,  _Entanglement Entropy of Fermions in Any Dimension and the Widom Conjecture_ ,  Phys. Rev. Lett. 96, 100503 (2006),  10.1103/PhysRevLett.96.100503.
  * [11] E. Ardonne, P. Fendley and E. Fradkin,  _Topological order and conformal quantum critical points_ ,  Annals of Physics 310(2), 493 (2004),  https://doi.org/10.1016/j.aop.2004.01.004.
  * [12] S. Chakravarty,  _Theory of the d-density wave from a vertex model and its implications_ ,  Phys. Rev. B 66, 224505 (2002),  10.1103/PhysRevB.66.224505.
  * [13] O. F. Syljuåsen and S. Chakravarty,  _Resonating plaquette phase of a quantum six-vertex model_ ,  Phys. Rev. Lett. 96, 147004 (2006),  10.1103/PhysRevLett.96.147004.
  * [14] Z. Yan, Z. Y. Meng, D. A. Huse and A. Chan,  _Height-conserving quantum dimer models_ ,  Phys. Rev. B 106, L041115 (2022),  10.1103/PhysRevB.106.L041115.
  * [15] S. Bravyi, L. Caha, R. Movassagh, D. Nagaj and P. W. Shor,  _Criticality without frustration for quantum spin-1 chains_ ,  Phys. Rev. Lett. 109, 207202 (2012),  10.1103/PhysRevLett.109.207202.
  * [16] R. Movassagh and P. W. Shor,  _Supercritical entanglement in local systems: Counterexample to the area law for quantum matter_ ,  Proceedings of the National Academy of Sciences 113(47), 13278 (2016),  10.1073/pnas.1605716113,  https://www.pnas.org/content/113/47/13278.full.pdf.
  * [17] Z. Zhang, A. Ahmadain and I. Klich,  _Novel quantum phase transition from bounded to extensive entanglement_ ,  Proceedings of the National Academy of Sciences 114(20), 5142 (2017),  10.1073/pnas.1702029114,  https://www.pnas.org/content/114/20/5142.full.pdf.
  * [18] L. Dell’Anna, O. Salberger, L. Barbiero, A. Trombettoni and V. E. Korepin,  _Violation of cluster decomposition and absence of light cones in local integer and half-integer spin chains_ ,  Phys. Rev. B 94, 155140 (2016),  10.1103/PhysRevB.94.155140.
  * [19] O. Salberger and V. Korepin,  _Entangled spin chain_ ,  Reviews in Mathematical Physics 29(10), 1750031 (2017),  10.1142/S0129055X17500313.
  * [20] O. Salberger, T. Udagawa, Z. Zhang, H. Katsura, I. Klich and V. Korepin,  _Deformed fredkin spin chain with extensive entanglement_ ,  Journal of Statistical Mechanics: Theory and Experiment 2017(6), 063103 (2017),  10.1088/1742-5468/aa6b1f.
  * [21] Z. Zhang and I. Klich,  _Entropy, gap and a multi-parameter deformation of the fredkin spin chain_ ,  Journal of Physics A: Mathematical and Theoretical 50(42), 425201 (2017),  10.1088/1751-8121/aa866e.
  * [22] Z. Zhang and H. S. Røising,  _Hilbert space fragmentation in a frustration-free fully packed loop model_ ,  10.48550/ARXIV.2206.01758 (2022).
  * [23] Z. Zhang and I. Klich,  _Entanglement phase transition of colored quantum dimers on the honeycomb lattice_ ,  10.48550/ARXIV.2210.01098 (2022).
  * [24] E. H. Lieb,  _Exact solution of the problem of the entropy of two-dimensional ice_ ,  Phys. Rev. Lett. 18, 692 (1967),  10.1103/PhysRevLett.18.692.
  * [25] R. Baxter,  _Exactly Solved Models in Statistical Mechanics_ ,  Dover books on physics. Dover Publications,  ISBN 9780486462714 (2007).
  * [26] Z. Zhang and G. Mussardo,  _Hidden bethe states in a partially integrable model_ ,  Phys. Rev. B 106, 134420 (2022),  10.1103/PhysRevB.106.134420.
  * [27] G. Mussardo, A. Trombettoni and Z. Zhang,  _Prime suspects in a quantum ladder_ ,  Phys. Rev. Lett. 125, 240603 (2020),  10.1103/PhysRevLett.125.240603.
  * [28] X. Chen, E. Fradkin and W. Witczak-Krempa,  _Gapless quantum spin chains: multiple dynamics and conformal wavefunctions_ ,  Journal of Physics A: Mathematical and Theoretical 50(46), 464002 (2017),  10.1088/1751-8121/aa8dbc.
  * [29] X. Chen, E. Fradkin and W. Witczak-Krempa,  _Quantum spin chains with multiple dynamics_ ,  Phys. Rev. B 96, 180402 (2017),  10.1103/PhysRevB.96.180402.
  * [30] H. Cohn, R. Kenyon and J. Propp,  _A variational principle for domino tilings_ ,  Journal of the American Mathematical Society 14(2), 297 (2001).
  * [31] R. Kenyon and A. Okounkov,  _Limit shapes and the complex burgers equation_ ,  10.48550/ARXIV.MATH-PH/0507007 (2005).
  * [32] N. Destainville,  _Entropy and boundary conditions in random lozenge tilings_ ,  10.48550/ARXIV.COND-MAT/9804062 (1998).
  * [33] E. Granet, L. Budzynski, J. Dubail and J. L. Jacobsen,  _Inhomogeneous gaussian free field inside the interacting arctic curve_ ,  Journal of Statistical Mechanics: Theory and Experiment 2019(1), 013102 (2019),  10.1088/1742-5468/aaf71b.
  * [34] N. Reshetikhin and A. Sridhar,  _Integrability of limit shapes of the six vertex model_ ,  Communications in Mathematical Physics 356(2), 535 (2017),  10.1007/s00220-017-2983-x.
  * [35] N. Reshetikhin and A. Sridhar,  _Limit shapes of the stochastic six vertex model_ ,  Communications in Mathematical Physics 363(3), 741 (2018),  10.1007/s00220-018-3253-2.
  * [36] P. Lammers and M. Tassy,  _Macroscopic behavior of lipschitz random surfaces_ ,  10.48550/ARXIV.2004.15025 (2020).
  * [37] H. Duminil-Copin, K. K. Kozlowski, D. Krachun, I. Manolescu and T. Tikhonovskaia,  _On the six-vertex model’s free energy_ ,  Communications in Mathematical Physics (2022),  10.1007/s00220-022-04459-x.
  * [38] J.-D. Deuschel and G. Giacomin,  _Entropic repulsion for massless fields_ ,  Stochastic Processes and their Applications 89(2), 333 (2000),  https://doi.org/10.1016/S0304-4149(00)00030-2.
  * [39] A. Borodin, V. Gorin and E. M. Rains,  _q-distributions on boxed plane partitions_ ,  Selecta Mathematica 16(4), 731 (2010),  10.1007/s00029-010-0034-y.
  * [40] L. Levine and R. Movassagh,  _The gap of the area-weighted motzkin spin chain is exponentially small_ ,  Journal of Physics A: Mathematical and Theoretical 50(25), 255302 (2017).
  * [41] Y. M. M. Knops, B. Nienhuis, H. J. F. Knops and H. W. J. Blöte,  _19-vertex version of the fully frustrated xy model_ ,  Phys. Rev. B 50, 1061 (1994),  10.1103/PhysRevB.50.1061.
  * [42] Y. Honda and T. Horiguchi,  _Critical behavior of a 19-vertex model with full frustration_ ,  Phys. Rev. B 58, 11501 (1998),  10.1103/PhysRevB.58.11501.
  * [43] T. KOJIMA,  _The 19-vertex model at critical regime $|$q$|$=1_,  International Journal of Modern Physics A 16(09), 1559 (2001),  10.1142/S0217751X01003445.
  * [44] A. Lima-Santos,  _Bethe ansätze for 19-vertex models_ ,  Journal of Physics A: Mathematical and General 32(10), 1819 (1999),  10.1088/0305-4470/32/10/004.
  * [45] C. Hagendorf,  _The nineteen-vertex model and alternating sign matrices_ ,  Journal of Statistical Mechanics: Theory and Experiment 2015(1), P01017 (2015),  10.1088/1742-5468/2015/01/P01017.
  * [46] K. Eloranta,  _The bounded 19-vertex model_ ,  10.48550/ARXIV.1710.03609 (2017).
  * [47] A. Bossart and W. Galleas,  _Functional relations in nineteen-vertex models with domain-wall boundaries_ ,  Journal of Mathematical Physics 60(10), 103509 (2019).
  * [48] A. B. Zamolodchikov and V. A. Fateev,  _Model factorized s-matrix and an integrable spin-1 heisenberg chain_ ,  Sov. J. Nucl. Phys. 32:2 (1980).
  * [49] A. Klümper, S. I. Matveenko and J. Zittartz,  _Exact solution of new integrable nineteen-vertex models and quantum spin-1 chains_ ,  Zeitschrift für Physik B Condensed Matter 96(3), 401 (1995).
  * [50] A. G. Izergin, D. A. Coker and V. E. Korepin,  _Determinant formula for the six-vertex model_ ,  Journal of Physics A: Mathematical and General 25(16), 4315 (1992),  10.1088/0305-4470/25/16/010.
  * [51] F. Colomo and A. G. Pronko,  _On two-point boundary correlations in the six-vertex model with domain wall boundary conditions_ ,  Journal of Statistical Mechanics: Theory and Experiment 2005(05), P05010 (2005),  10.1088/1742-5468/2005/05/p05010.
  * [52] F. Colomo and A. G. Pronko,  _An approach for calculating correlation functions in the six-vertex model with domain wall boundary conditions_ ,  Theoretical and Mathematical Physics 171(2), 641 (2012),  10.1007/s11232-012-0061-2.
  * [53] P. Belov and N. Reshetikhin,  _The two-point correlation function in the six-vertex model_ ,  Journal of Physics A: Mathematical and Theoretical 55(15), 155001 (2022),  10.1088/1751-8121/ac578e.
  * [54] R. N. Alexander, A. Ahmadain, Z. Zhang and I. Klich,  _Exact rainbow tensor networks for the colorful motzkin and fredkin spin chains_ ,  Phys. Rev. B 100, 214430 (2019),  10.1103/PhysRevB.100.214430.
  * [55] R. N. Alexander, G. Evenbly and I. Klich,  _Exact holographic tensor networks for the Motzkin spin chain_ ,  Quantum 5, 546 (2021),  10.22331/q-2021-09-21-546.
  * [56] L. Causer, J. P. Garrahan and A. Lamacraft,  _Slow dynamics and large deviations in classical stochastic fredkin chains_ ,  Phys. Rev. E 106, 014128 (2022),  10.1103/PhysRevE.106.014128.
  * [57] C. M. Langlett and S. Xu,  _Hilbert space fragmentation and exact scars of generalized fredkin spin chains_ ,  Phys. Rev. B 103, L220304 (2021),  10.1103/PhysRevB.103.L220304.
  * [58] J. Richter and A. Pal,  _Anomalous hydrodynamics in a class of scarred frustration-free hamiltonians_ ,  Phys. Rev. Research 4, L012003 (2022),  10.1103/PhysRevResearch.4.L012003.