# Variational manifolds for ground states and scarred dynamics of blockade-
constrained spin models on two and three dimensional lattices

Joey Li  Institute for Theoretical Physics, University of Innsbruck, 6020
Innsbruck, Austria  Institute for Quantum Optics and Quantum Information of
the Austrian Academy of Sciences, 6020 Innsbruck, Austria  Giuliano Giudici
Institute for Theoretical Physics, University of Innsbruck, 6020 Innsbruck,
Austria  Institute for Quantum Optics and Quantum Information of the Austrian
Academy of Sciences, 6020 Innsbruck, Austria  Hannes Pichler  Institute for
Theoretical Physics, University of Innsbruck, 6020 Innsbruck, Austria
Institute for Quantum Optics and Quantum Information of the Austrian Academy
of Sciences, 6020 Innsbruck, Austria

###### Abstract

We introduce a variational manifold of simple tensor network states for the
study of a family of constrained models that describe spin-1/2 systems as
realized by Rydberg atom arrays. Our manifold permits analytical calculation
via perturbative expansion of one- and two-point functions in arbitrary
spatial dimensions and allows for efficient computation of the matrix elements
required for variational energy minimization and variational time evolution in
up to three dimensions. We apply this framework to the PXP model on the
hypercubic lattice in 1D, 2D, and 3D, and show that, in each case, it exhibits
quantum phase transitions breaking the sub-lattice symmetry in equilibrium,
and hosts quantum many body scars out of equilibrium. We demonstrate that our
variational ansatz qualitatively captures all these phenomena and predicts key
quantities with an accuracy that increases with the dimensionality of the
lattice, and conclude that our method can be interpreted as a generalization
of mean-field theory to constrained spin models.

##  I Introduction

Arrays of neutral atoms trapped in optical tweezers [1, 2, 3, 4] have recently
emerged as a promising platform for the quantum simulation of many-body spin
models. These Rydberg atom quantum simulators have been used for the
experimental study of a variety of quantum phenomena in equilibrium, including
zero temperature quantum phase transitions between trivially disordered states
and states that break a variety of spatial symmetries [5, 6, 7, 8], as well as
spin liquids with topological order [9, 10, 11, 12, 13, 14] and symmetry
protected topological phases [15, 16]. They can also be used to study
thermalization dynamics [17, 18] and were instrumental in the experimental
discovery [19, 20] of the out-of-equilibrium phenomenon now known as quantum
many-body scars (QMBS) [21], which have sparked interest as an example of non-
thermalizing behaviour in quantum many-body systems, violating the eigenstate
thermalization hypothesis [22, 23]. Remarkably, often the physics of these
systems simply emerges from the Rydberg blockade mechanism [24, 25, 26] that
forbids neighbouring Rydberg excitations, introducing a constraint between
otherwise freely rotating spins. It is this constraint that renders the
problems in general non-trivial and thus underlies the difficulty of
describing these systems in a many-body setting.

Figure 1: Schematic depiction of the setups considered in this work,
consisting of PXP models in 1D, 2D and 3D hyper-cubic lattices. The blockade-
constraint in these systems does not allow two neighboring spins to be
simultaneously in the Rydberg state. After imposing unit cell translation
invariance, the system is parameterized by two variational parameters
$\theta_{A}$ and $\theta_{B}$, corresponding to the two sub-lattices.

Nevertheless, we expect that it should be possible to capture some of the
above phenomena by an appropriately constructed mean-field theory for these
constrained systems. This includes in particular phase transitions breaking
the sublattice symmetry, as well as the dynamics associated with QMBS. On a
technical level such an attempt faces the challenge of formulating a proper
variational manifold of states that satisfies the blockade constraint, and at
the same time allows for an efficient calculation of physical observables. In
this work we address this challenge and introduce a manifold of constraint-
satisfying states that can be elegantly represented as minimally entangled
tensor networks. Our ansatz can be seen as a generalization of minimal matrix-
product wavefunctions that have been employed for 1D chains [27, 28] and the
related tree networks in Ref. [29].

Going beyond these previous results, we show that our variational ansatz
allows for efficient calculation of expectation values, in the relevant
parameter regimes, in 2D square and 3D cubic lattices directly in the
thermodynamic limit. In particular, we recover the quantum phase transition
between a disordered and a symmetry broken phase in the ground states of
generalized PXP models, as well as the periodic dynamics of special initial
states associated to quantum scars in PXP models, via a time-dependent version
of the variational principle [30, 31]. In 1D, the phase diagram has been
studied analytically [32], and in 1D and 2D, these phenomena have been
observed experimentally [19, 6, 20] and studied numerically with density
matrix renormalization group techniques [33]. Here our method provides a
complementary viewpoint and adds the benefit of an analytical description. In
addition, it also treats the so far unexplored case of 3D lattices, and allows
for quantitative predictions, such as the location of the quantum phase
transition point in equilibrium, or the period of revival in scarred dynamics.

The paper is structured as follows. We first introduce the class of PXP models
analyzed in this work in Sec. II and review the basic phenomena that this
model is expected to display. In Sec. III we introduce our variational
manifold and we present a method that allow us to perform tensor network
contractions in two and three dimensions in the form of a power series. In
Sec. IV we discuss the ground state phase diagram of our model in 1D, 2D and
3D using our variational ansatz, and compare its predictions with exact
diagonalization results. In Sec. V we discuss time evolution obtained via the
TDVP and contrast its predictions with the exact dynamics of the PXP model.

##  II Model

###  II.1 Constrained Hilbert space

In this work we are interested in constrained spin models that are often
referred to as PXP models. Specifically, we consider $N$ spin-1/2 particles,
each with basis states $\ket{\downarrow}$ and $\ket{\uparrow}$, that are
defined on the sites of a $D$-dimensional lattice, as depicted in Fig. 1. We
consider states of this spin systems that satisfy a constraint that prevents
two spins on neighbouring lattice sites to be simultaneously in the
$\ket{\uparrow}$ state. In physical realizations this constraint arises
dynamically from strong state-dependent nearest neighbour interactions.
Specifically, this constraint captures the essence of the Rydberg blockade
mechanism in Rydberg atom arrays. The spin configurations that satisfy this
constraint form a subspace of the $2^{N}$ dimensional Hilbert space. For large
system sizes, the dimension of this constrained Hilbert space grows as
$x^{N}$, where the value of $x<2$ depends on the lattice. For a 1D chain,
$x=(1+\sqrt{5})/2\simeq 1.618$ 111The Hilbert space dimension of $N$ spins
$d_{H}(N)$ satisfies the equation $d_{H}(N+1)=d_{H}(N)+d_{H}(N-1)$, which
implies $d_{H}\sim\phi^{N}$ for large $N$, where $\phi$ is the golden ratio..
In higher dimensions the larger coordination number increases the effect of
the constraint, leading to a decrease of $x$ with higher $D$.

In the following we denote the projector onto the constraint-satisfying
subspace by $\mathcal{P}$. It can be composed from local projectors
$\mathcal{P}_{ij}$ that act on pairs of neighboring spins $i$ and $j$, i.e.,
$\mathcal{P}=\prod_{\langle ij\rangle}\mathcal{P}_{ij}$. Here
$\mathcal{P}_{ij}=\ket{\downarrow\downarrow}\bra{\downarrow\downarrow}+\ket{\downarrow\uparrow}\bra{\downarrow\uparrow}+\ket{\uparrow\downarrow}\bra{\uparrow\downarrow}$
simply annihilates components that do not satisfy the constraint on the link
between $i$ and $j$.

###  II.2 Hamiltonian

Let us now write down the Hamiltonian of a spin model that lives in this
constrained Hilbert space; the simplest example is the PXP model

$\displaystyle
H_{\text{PXP}}=\Omega\sum_{i}\mathcal{P}\sigma^{x}_{i}\mathcal{P},$ (1)

where
$\sigma^{x}=\ket{\uparrow}\bra{\downarrow}+\ket{\downarrow}\bra{\uparrow}$.
The Hamiltonian (1) can be interpreted as the projection of a non-interacting
model of independently rotating spins into the constraint-satisfying subspace.
It induces a conditional dynamics, where a spin only precesses if all its
neighbours are in the $\ket{\downarrow}$ state. Crucially this projected
Hamiltonian is interacting, and in general non-integrable [35, 36]. Note that
the geometry of the lattice enters in (1) implicitly via the projection
operation, which manifestly depends on the lattice structure.

In this work we study the physics of what we shall call generalized PXP
models, the family of Hamiltonians of the PXP model (1) with the addition of a
transverse-field term and a next-nearest-neighbour interaction:

$\displaystyle H=\sum_{i}\mathcal{P}(\Omega\sigma^{x}_{i}-\Delta
n_{i})\mathcal{P}+V\sum_{\langle\\!\langle
ij\rangle\\!\rangle}\mathcal{P}n_{i}n_{j}\mathcal{P},$ (2)

where $n=\ket{\uparrow}\bra{\uparrow}$.

This Hamiltonian may be interpreted as a toy model for an array of coherently
driven Rydberg atoms. In these systems the state $\ket{\downarrow}$ represents
an internal electronic ground state of an atom and the state $\ket{\uparrow}$
represents a highly excited Rydberg state. When driven by a laser that couples
these two states, the system is described by a Hamiltonian of the form

$\displaystyle H_{\text{Rydberg}}=\sum_{i}(\Omega\sigma^{x}_{i}-\Delta
n_{i})+\sum_{ij}V_{ij}n_{i}n_{j}.$ (3)

Here $\Omega$ and $\Delta$ denote the Rabi frequency and the laser detuning,
and $V_{ij}=C/|i-j|^{6}$ is the van der Waals interaction between two atoms
$i$ and $j$ that are both in the Rydberg state. The strength of the
interaction depends on the geometric distance between the atoms, with a length
scale set by the blockade radius $R_{b}=(C/\Omega)^{1/6}$. In the so-called
blockade approximation one considers the limit $V_{ij}\rightarrow\infty$ if
two atoms are less than one blockade radius apart, and $V_{ij}\rightarrow 0$
otherwise. In particular, when $R_{b}$ equals one lattice spacing $a$, the
model reduces to Eq. (2) with $V=0$, which becomes the PXP model for
$\Delta=0$. If one sets $V_{ij}=V$ for next-nearest-neighbour sites, i.e. for
$|i-j|=2a$ in 1D or $|i-j|=\sqrt{2}a$ in $D>1$, then one recovers the
Hamiltonian in Eq. (2). While it is natural to have $V>0$, it is also
theoretically and experimentally [37] interesting to consider the case $V<0$,
i.e. with attractive next nearest neighbour interactions. In what follows, we
will set $\Omega=1$.

###  II.3 Physics of the PXP model on bipartite lattices

Hamiltonians of the form (2) host a remarkably rich variety of physical
phenomena. For instance, depending on the choice of the lattice, which selects
the form of $\mathcal{P}$, the ground state of (2) can host a plethora of
symmetry broken ordered phases [32, 33, 38] as well as topologically ordered
phases [39]. In this work we focus on bipartite lattices, in particular on
hypercubic lattices, and denote the two sublattices by $A$ and $B$
respectively. For such bipartite lattices we now briefly review the expected
qualitative features of the system in and out of equilibrium.

####  II.3.1 Equilibrium

The ground state phase diagram of the Hamiltonian in Eq. (2), at fixed finite
values of $V$, depends only on $\Delta$. For $\Delta\to-\infty$, the ground
state is unique and given by $\ket{\downarrow}^{\otimes N}$ independent of the
lattice geometry. For $\Delta\to+\infty$ and any finite $V$, the ground state
maximizes the number of spins in the $\ket{\uparrow}$ state that is consistent
with the constraints. For a bipartite lattice this gives a twofold degenerate
ground state, which spontaneously breaks the sublattice symmetry: the two
ground states are given by configurations where all spins on sublattice $A$
are in the $\ket{\uparrow}$ state and all spins on sublattice $B$ are in the
$\ket{\downarrow}$ state, or viceversa. We denote these states by
$\ket{\mathbb{Z}_{2}}=\ket{\uparrow}^{\otimes A}\ket{\downarrow}^{\otimes B}$
and $\ket{\mathbb{Z}_{2}^{\prime}}=\ket{\downarrow}^{\otimes
A}\ket{\uparrow}^{\otimes B}$, respectively. In the limit
$\Delta,V\to+\infty$, with constant $\Delta/V=c$, the ground state can be
different. In 1D for $c>3$ a period-3 density-wave state with $1/3$ density of
up spins is favored. In 2D and 3D the states
$\ket{\mathbb{Z}_{2}},\ket{\mathbb{Z}_{2}^{\prime}}$ cease to be the classical
ground states for $c>4$ and $c>6$, respectively, in favor of states with
density $1/4$, leading to a four-fold degenerate striated phase which arises
perturbatively due to quantum fluctuations [33]. The latter can stabilize
these non-$\mathbb{Z}_{2}$-ordered phases even at finite $V$. In fact, in 1D,
the interplay between $\mathbb{Z}_{2}$ and $\mathbb{Z}_{3}$ order gives rise
to several interesting phenomena such as floating phases and chiral
transitions [40, 41, 42, 43]. In $D>1$, the precise nature of the phase
transition along all phase boundaries of (2) is less clear [33]. Here, we
focus on the regime where $V$, when positive, is sufficiently small, such that
as $\Delta$ is tuned from $-\infty$ to $+\infty$ the system undergoes a direct
quantum phase transition from the gapped, disordered phase to the gapped,
$\mathbb{Z}_{2}$ ordered phase. When $D=1$ and for $V\geq 0$, the phase
transition is second order. For sufficiently strong attractive next-nearest-
neighbour interaction $V<0$, the phase transition becomes discontinuous and
first order. A tricritical point for $V<0$ separates these two regimes. The
most accurate numerical study of the entire transition line was performed in
Ref. [40]. The location of the tricritical point is known analytically, since
it lies on the integrable line $\Delta=V-1/V$ for $V=-(\sqrt{5}+1)/2)^{5/2}$
[32]. A sketch of the $(\Delta,V)$ phase diagram is shown in Fig. 2. We expect
a similar phase diagram to hold in $D>1$ and, as we show below, our
variational method predicts that this is indeed the case; in particular, we
will use it to locate the phase boundary and the tricritical point in
$D=1,2,3$ and we will compare the outcome with exact diagonalization results
on 2D and 3D lattices with up to $48$ and $64$ sites, respectively. We note
that, whereas in 1D the transition is always surrounded by a disordered regime
for any finite $\Delta,V$, in 2D and for large enough $\Delta$ the system
undergoes a direct second order transition from the $\mathbb{Z}_{2}$ phase to
the striated phase, with $V$ as driving parameter [33, 44]. Since our ansatz
only captures $\mathbb{Z}_{2}$ order, we will restrict our discussion to the
regime where the effect of the $1/4$-density ordered phase is negligible.

Figure 2: Sketch of the phase diagram for the generalized PXP model with next-
nearest-neighbour interaction Eq. (2). A disordered phase and a
$\mathbb{Z}_{2}$ ordered phase are separated by second-order (dotted line) and
first order (solid line) transitions, which meet at the tricritical point (red
dot). Quantitative versions of this for 1D, 2D, and 3D are shown in Fig. 6.

####  II.3.2 Quantum many-body scars

When $\Delta=V=0$, the Hamiltonian in Eq. (2) reduces to the PXP model which
is known to host quantum many-body scars. These are special eigenstates that
do not satisfy the eigenstate thermalization hypothesis [22, 23] and have been
proposed in Ref. [45] as the explanation of the anomalous dynamics observed in
a Rydberg atom simulator initialized in one of the two $\ket{\mathbb{Z}_{2}}$
states [19]. The anomaly consists in the fact that despite the Hamiltonian (1)
being non-integrable, and the initial $\ket{\mathbb{Z}_{2}}$ state having the
energy density of an infinite temperature ensemble, the system does not
rapidly thermalize, but rather exhibits (approximate) periodic dynamics with
alternating revivals of the two $\ket{\mathbb{Z}_{2}}$ states. It is now
believed that this phenomenon arises in the PXP model in any dimension
whenever the lattice is bipartite [46], although numerical studies in $D>2$
are absent in the literature. The name quantum scars comes from the analogy,
put forward in Ref. [45], with single-particle physics, where the quantum
system possesses special eigenstates that are “scars” of unstable periodic
orbits of their, otherwise chaotic, classical counterpart [47]. In this
context, the one-dimensional version of the variational ansatz which we will
introduce in Sec. III was employed in Ref. [27] to map the constrained
dynamics of the PXP model into classical equations of motion where the
$\ket{\mathbb{Z}_{2}}$ states indeed lie on unstable periodic orbits. Below we
generalize these findings to $D>1$ and demonstrate the predictive power of our
variational ansatz in two and three dimensions. A characteristic quantity that
we can compute in our framework is the revival time, which we find to increase
with the coordination number of the lattice. By comparing the result with the
exact many-body quantum dynamics on small systems, we will show that our
method predicts the period with an accuracy increasing with the lattice
dimension.

Figure 3: Schematic illustration of the periodic revivals in quantum scarred
systems, and how they are captured by the variational method. The 2D surface
represents the variational manifold, and the loop on it represents the exact
periodic orbits in the variational dynamics connecting the two
$\ket{\mathbb{Z}_{2}}$ states. The dotted line represents the true time
evolved state $e^{-iH_{\text{PXP}}t}\ket{\mathbb{Z}_{2}}$.

##  III Variational ansatz

We now introduce a variational manifold of states that is designed to capture
some of the physics of PXP-type models on bipartite lattices. The manifold of
states $|\psi(\vec{\theta},\vec{\phi})\rangle$ can be represented as tensor
network states with bond dimension $\chi=2$ on the network given by the
lattice of interest. The key feature of our variational manifold is that the
states satisfy the blockade constraint on the corresponding lattice by
construction, i.e., they do not contain configurations where two neighboring
spins are both in the $\ket{\uparrow}$ state. The variational manifold is
parameterized by two angles, $\theta_{i}\in[-\pi,\pi)$ and
$\phi_{i}\in[0,2\pi)$, for each lattice site $i$. We emphasize that for the
applications presented in this paper, we only work with two-site unit-cell
translationally invariant systems (i.e. systems invariant under a translation
by an even number of sites in any direction) with parameter space
$(\theta_{A},\theta_{B},\phi_{A},\phi_{B})$, although in the discussion here
we will refer to the general case $(\vec{\theta},\vec{\phi})$ when it presents
no additional difficulty.

The variational state $|\psi(\vec{\theta},\vec{\phi})\rangle$ is a projected
entangled pair state (PEPS)[48] with a PEPS tensor $M(\theta_{i},\phi_{i})$
defined at each lattice site $i$. In 1D, the state is a matrix product state
(MPS), and the MPS tensor is

$\displaystyle\quad\vbox{\hbox{\includegraphics[width=59.75095pt]{figs/MPS01.pdf}}}=M(\theta,\phi)^{(j)}_{ab}$
$\displaystyle=$
$\displaystyle\begin{pmatrix}\cos(\frac{\theta}{2})\ket{\downarrow}_{j}&-ie^{i\phi}\sin(\frac{\theta}{2})\ket{\uparrow}_{j}\\\
\ket{\downarrow}_{j}&0\end{pmatrix}_{ab}$ (4)

where $\ket{\downarrow}=(1\;\;0)^{T}$ and $\ket{\uparrow}=(0\;\;1)^{T}$, and
where $j$ is the physical index and $a$ and $b$ are the virtual indices. In
the graphical notation above we assigned a direction to the virtual legs; the
usefulness of this should be readily apparent when we consider higher
dimensional generalizations.

We generalize this ansatz as follows. On a 2D square lattice, the PEPS tensors
take the form

$\displaystyle\vbox{\hbox{\includegraphics[width=54.06006pt]{figs/Mtensor.pdf}}}\;=\;M(\theta,\phi)^{(j)}_{(ab)(cd)}$
$\displaystyle=$
$\displaystyle\,\begin{pmatrix}\cos(\frac{\theta}{2})\ket{\downarrow}_{j}&\;0\;\;&\;\;0\;&-ie^{i\phi}\sin(\frac{\theta}{2})\ket{\uparrow}_{j}\\\
\ket{\downarrow}_{j}&\;0\;\;&\;\;0\;&0\\\
\ket{\downarrow}_{j}&\;0\;\;&\;\;0\;&0\\\
\ket{\downarrow}_{j}&\;0\;\;&\;\;0\;&0\end{pmatrix}_{(ab)(cd)},$ (5)

where the matrix is written in the basis $(00,01,10,11)$. Fig. 4 depicts the
2D PEPS for a system with a two-site unit cell translation symmetry imposed.

Figure 4: PEPS for the 2D square lattice

The above construction suggests a generalization that can be applied to
arbitrary dimensional lattices: in general each tensor will have some “in”
virtual indices, and some “out” virtual indices. For a hypercubic lattice in
$D$ dimensions, each PEPS tensor will have $2D+1$ total legs, with $D$
incoming virtual indices and $D$ outgoing virtual indices. Diagrammatically,
we write

$\displaystyle\vbox{\hbox{\includegraphics[width=79.6678pt]{figs/PEPSnD.pdf}}}=M(\theta,\phi)^{(i)}_{(\text{in})(\text{out})},$
(6)

where a thick line collectively represents all incoming or outgoing legs.
Introducing a bra-ket notation on the virtual indices, the tensor $M$ reads

$\displaystyle M=\sum_{i=\downarrow,\uparrow}\sum_{\text{in,
out}}M^{(i)}_{(\text{in})(\text{out})}\ket{i}\ket{\text{in}}\bra{\text{out}}.$
(7)

Let $I\in\\{0,1\\}^{D}$ be a binary string of length $D$, let $\sum_{I}$ be
the sum over all $2^{D}$ such strings, and denote by $\ket{I}$ the
corresponding product state in the virtual Hilbert space. We further use the
short-hand notation $\vec{0}=(0,0,0,\dots)$ and $\vec{1}=(1,1,1,\dots)$, and
define $\ket{\alpha}\equiv\sum_{I}\ket{I}$ and
$\ket{\beta}\equiv\sum_{I\neq\vec{0}}\ket{I}=\ket{\alpha}-|\vec{0}\rangle$.
Using this notation, the tensor $M$ can be written in general dimensions as

$\displaystyle
M(\theta,\phi)=\ket{\downarrow}\big{(}\cos(\theta/2)|\vec{0}\rangle\langle\vec{0}|+\ket{\beta}\langle\vec{0}|\big{)}$
$\displaystyle-
ie^{i\phi}\sin(\theta/2)\ket{\uparrow}|\vec{0}\rangle\langle\vec{1}|.$ (8)

It is straightforward to see that this reproduces Eqs. (III) and (III) for
$D=1,2$, respectively.

The variational ansatz $\ket{\psi}$ has several interesting properties. As
already noted, it satisfies the blockade constraint. Remarkably, it is always
normalized, i.e.
$\langle\psi(\vec{\theta},\vec{\phi})|\psi(\vec{\theta},\vec{\phi})\rangle=1$,
in the thermodynamic limit; this follows from a simple but rather technical
calculation which we present in Appendix A. An important consequence of the
normalization is that it allows one to calculate also expectation values and
time evolution of states using the TDVP. For bipartite lattices, the
variational manifold contains the important state
$\ket{\mathbb{Z}_{2}}=\ket{\downarrow}^{\otimes A}\ket{\uparrow}^{\otimes B}$;
this can be checked by setting $\theta_{i}=\pm\pi,\theta_{j}=0$ for all sites
$i\in A,j\in B$. This is in fact a special case of a more general property:
product states with all spins on one sublattice in the state
$\ket{\downarrow}$, and with all spins on the other sublattice being any state
on the Bloch sphere, are contained in the variational manifold. In $D=1$ the
variational manifold is known to be equivalent to the manifold of all product
states projected into the constraint-satisfying Hilbert space [27]. This is
not the case in $D>1$. In fact, the projected product state has a different
tensor network representation that is not gauge-equivalent to Eq. (III) and
does not yield a normalized state in the thermodynamic limit. Since this
property provides remarkable simplifications in the following variational
calculations, we will stick to the ansatz in Eq. (III) and provide a more in-
depth discussion about its relation to the projected product state manifold in
Appendix D.

####  III.0.1 $(\theta_{A},\theta_{B})$ space for systems on bipartite
lattices

For the applications that follow, we only consider two-site unit-cell
translational invariant states. The $\phi$ variables are a priori relevant,
but we will show that for the variational dynamics, the trajectories contain
only states with $\phi_{A},\phi_{B}=0$, and that the variational ground states
have $\phi_{A},\phi_{B}=\pm\pi/2$. Thus, the value of $\phi_{A}$ and
$\phi_{B}$ can be considered fixed, and we have a two dimensional parameter
space $\theta_{A},\theta_{B}\in[-\pi,\pi)$, which contains the states relevant
for variational energy minimization, and in which the variational time
evolution takes place. This space is shown in Fig. 5, where we label the
$\ket{\mathbb{Z}_{2}}$ states and the product states along the $\theta_{A}$
and $\theta_{B}$ axes, as well the non-perturbative and perturbative regimes
for the series expansion calculation of expectation values (discussed below).

Figure 5: Depiction of the $(\theta_{A},\theta_{B})$ space that our
variational state $\ket{\psi(\theta_{A},\theta_{B})}$ is in. Points along the
axes are all product states, where we have defined
$\ket{\theta\phi}=\cos(\theta/2)\ket{\downarrow}-ie^{i\phi}\sin(\theta/2)\ket{\uparrow}$,
and $\ket{\mathbb{Z}_{2}}=\ket{\uparrow}^{\otimes A}\ket{\downarrow}^{\otimes
B}$ and $\ket{\mathbb{Z}_{2}^{\prime}}=\ket{\downarrow}^{\otimes
A}\ket{\uparrow}^{\otimes B}$. (For the applications here $\phi$ is irrelevant
and can be fixed to $0$ or $\pm\pi/2$.) The shaded regions correspond to the
non-perturbative regime (which is slightly larger in 3D than in 2D). In 1D and
2D we can access the entire space via explicit contraction, while in 3D we are
limited to the unshaded region. However, we find that the perturbative regime
captures all of the physics discussed in this work.

###  III.1 Calculation of expectation values

The structure of the variational state
$\ket{\psi}=|\psi(\vec{\theta},\vec{\phi})\rangle$ allows us to efficiently
calculate quantities of the form $\langle\psi|\hat{\mathcal{O}}|\psi\rangle$,
for product operators $\hat{\mathcal{O}}$ in up to three dimensions. This
includes in particular local observables but also $n$-point functions,
although in this work we only consider one- and two-point functions. In this
section we give an outline of the methods we use; a full discussion is in
Appendix A.

Figure 6: $(\Delta,V)$ phase diagrams for the generalized PXP model (2) with
the magnetization
$\langle\sigma^{z}_{A}\rangle_{GS}-\langle\sigma^{z}_{B}\rangle_{GS}$ as the
order parameter. The red dots denote the tricritical points (as found from the
variational method) while the lines show, for comparison, the phase boundaries
found via exact diagonalization on various system sizes. Finite-size critical
points are extracted from the peaks in the ground state fidelity
susceptibility (see Appendix G for details). In 1D, we also show the
analytically known tricritical point.

In general, for a bond-dimension $\chi$ tensor network state $\ket{\Psi}$,
evaluating quantities like $\langle\Psi|\hat{\mathcal{O}}|\Psi\rangle$
requires the contraction of a tensor network of bond dimension $\chi^{2}$.
However, even though the tensor networks states introduced above have bond
dimension $\chi=2$, it turns out that the evaluation of quantities of the form
$\langle\psi|\hat{\mathcal{O}}|\psi\rangle$ reduces to the contraction of a
tensor network, whose bond dimension is not 4 but instead only 2. For $D=1$
the tensor network can be contracted analytically for infinte systems and for
$D=2$, one can calculate $\langle\psi|\hat{\mathcal{O}}|\psi\rangle$ via
explicit tensor network contraction on infinite cylinders with finite
circumference. The reduced bond dimension allows one to perform the
calculation for up to relatively large circumference. We performed our
numerical calculations with cylinders of circumference $L=10$ sites, but we
verified that compared to larger cylinders the results are unchanged to
several significant figures.

For $D\geq 2$, we introduce a method to express the tensor network contraction
as a perturbative expansion. For the unit cell translation invariant case that
we consider in the applications below, the variational ansatz becomes
$\ket{\psi(\theta_{A},\theta_{B},\phi_{A},\phi_{B})}$. We show that it is
possible to express $\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$ as

$\displaystyle\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$
$\displaystyle=h(\theta_{A},\theta_{B},\phi_{A},\phi_{B})\Sigma(f,\theta_{A},\theta_{B}),$
$\displaystyle\Sigma(f,\theta_{A},\theta_{B})=$
$\displaystyle\sum_{n,m=0}^{\infty}\\!(-1)^{n+m}f_{n,m}\sin^{2n}\\!\\!\Big{(}\frac{\theta_{A}}{2}\Big{)}\sin^{2m}\\!\\!\Big{(}\frac{\theta_{B}}{2}\Big{)},$
(9)

where the prefactor $h(\theta_{A},\theta_{B},\phi_{A},\phi_{B})$ is a function
that depends on the operator $\hat{\mathcal{O}}$, and where $f$ is a matrix of
“counting factors”, that also depends on $\hat{\mathcal{O}}$, but is
independent of the variational parameters. The central idea underlying this
method is based on a decomposition of each tensor in the network into a sum of
two tensors. This seems counterproductive at first, since it expresses the
quantity $\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$ not as a contraction of a
single tensor network, but instead as sum over contractions of $2^{N}$ tensor
networks. However, with a clever choice of the decomposition, one can show
that only few of those contractions give a nonzero contribution. Moreover,
those contractions that give a finite contribution, can be easily evaluated.
These terms can be counted and collected into a series expansion with the
$\sin^{2}(\theta/2)$ as the expansion parameter. For more details on this
method we refer the reader to Appendix A.3.

We note that Eq. (III.1) should be understood as an asymptotic series. This is
because the counting factors $f_{n,m}$ scale super exponentially (that is,
there do not exist any $a,b\in\mathbb{R}$ such that $f_{n,m}<a^{n}b^{m}$ for
all $n,m\in\mathbb{N}$), causing the series to eventually diverge at high
enough order. However, we find that in regions of the
$(\theta_{A},\theta_{B})$ plane where the effective expansion parameters
$\sin^{2}(\theta_{A}/2)$ and $\sin^{2}(\theta_{B}/2)$ are small enough, i.e. a
star shaped region around the $x$ and $y$ axes, the terms of the series
decrease quickly, and accurate results are obtained at finite order. The
perturbative regimes in Fig. 5 were defined as the region where the expansion
is accurate to an error of $10^{-3}$. We verified that this perturbative
method is valid for the two applications we work with (energy minimization and
time evolution): for the former, we verify that our variational ground states
lie inside the perturbative region, and for the latter, we considered time-
evolved trajectories from different orders of the expansion, and verified that
they converge.

##  IV Ground state phase diagram

As a first application of the variational manifold, we calculate the energy
$E=\bra{\psi}H\ket{\psi}$ and study the properties of the variational ground
state. We consider the generalized PXP model including the detuning term and a
next-nearest-neighbour (NNN) term Eq. (2).

We use our variational manifold as a minimal ansatz beyond mean-field theory
to study the ground state of this Hamiltonian and to reproduce the phase
diagram. With translation invariance we have a four dimensional parameter
space with variational states
$\ket{\psi(\theta_{A},\theta_{B},\phi_{A},\phi_{B})}$. However, we can
restrict the manifold to real states by setting $\phi_{A}=\phi_{B}=\pi/2$, as
the Hamiltonian is real. We consider the energy function
$E(\theta_{A},\theta_{B})$ on the variational space
$\theta_{A},\theta_{B}\in[-\pi,\pi)$ and find the variational ground state by
minimizing $E(\theta_{A},\theta_{B})$. Fig. 6 shows the resulting phase
diagrams.

Let us first qualitatively describe the landscape of the optimal variational
parameters $\theta_{A}^{GS},\theta_{B}^{GS}$ as the Hamiltonian is tuned.
Suppose we fix $V$ and allow $\Delta$ to vary: when $\Delta=-\infty$, the
ground state is the $|\vec{0}\rangle$ state, and when $\Delta$ remains large
and negative, we find that $\theta_{A}^{GS}=\theta_{B}^{GS}$ at the minimum
energy state of the variational manifold. As $\Delta$ reaches a critical value
$\Delta_{c}$, the variational ground state changes from the symmetric phase to
a symmetry-broken phase where $\theta_{A}\neq\theta_{B}$. When
$\Delta>\Delta_{c}$, the variational ground state becomes twofold degenerate,
with nonzero $\theta_{A}^{GS}-\theta_{B}^{GS}$. For $V=0$, we find that the
transition occurs at $\Delta_{c}^{1D}\simeq 0.77$ in 1D, at
$\Delta_{c}^{2D}\simeq-0.45$ in 2D, and at $\Delta_{c}^{3D}\simeq-1.0$ in 3D.
In 1D it is known from exact diagonalization studies that $\Delta_{c}\simeq
1.31$ [49]. In 2D, exact diagonalization on small systems sizes combined with
finite-size scaling indicate that $\Delta_{c}\simeq-0.2$ (see Appendix G). In
3D, the critical point for 64 sites is at $\Delta_{c}\simeq-0.61$. Thus we see
that in all cases the variational method approximates, but underestimates the
transition point.

Figure 7: Plots of the minimum energy (per site) along lines of constant
$\theta_{A}-\theta_{B}$ for $D=1,2,3$, showing the second order transition
(top row) and first order transition (bottom row).

When considering different values of the NNN interaction $V$, we find that in
all three cases ($D=1,2,3$) the variational method predicts a change from a
second order transition to a first order transition as $V$ is decreased. For
$V>V_{c}$ the variational ground state $(\theta_{A}^{GS},\theta_{B}^{GS})$
evolves smoothly with $\Delta$, while for $V<V_{c}$ it changes abruptly at the
transition point. We find that $V_{c}\simeq-1.0,-0.6,$ and $-0.25$ in 1D, 2D,
and 3D, respectively. In 1D, this significantly overestimates the tricritical
point compared with the known analytical value of $V_{c}\simeq-3.33$ [32],
although it is still notable that the effect is captured by the variational
method. For $\Delta$, we find that
$\Delta_{\text{tricritical}}\simeq-0.51,-1.1,$ and $-1.5$ in 1D, 2D, and 3D,
respectively.

We present the data here in two series of plots. In Fig. 6 we plot the
$(\Delta,V)$ phase diagrams, using the order parameter
$\langle\sigma^{z}_{A}\rangle_{GS}-\langle\sigma^{z}_{B}\rangle_{GS}$. In Fig.
7 we present the one-dimensional plots of minimum energies along lines of
constant $\theta_{A}-\theta_{B}$, i.e.
$E_{\text{minimum}}(\zeta=\theta_{A}-\theta_{B})$ where
$E_{\text{minimum}}(\zeta)=\min_{\theta_{A}}[E(\theta_{A},\theta_{A}-\zeta)]$
, for various values of $\Delta$ at both $V=0$ and $V<V_{c}$. In these plots
one can clearly observe the continuous and discontinuous change of the ground
state in the second-order and first-order transition regimes, respectively.
Moreover, we extract the critical exponent of the order parameter by fitting
the numerical data close to the transition point at $V=0$, where the
transition is continuous, and find that it agrees with the Ising mean-field
critical exponent $\beta=1/2$ (see Fig. 8(c)).

Figure 8: (a) Correlation lengths $\xi$ as function of $\Delta-\Delta_{c}$ of
the variational ground state, for $V=0$. $\Delta_{c}=0.77,-0.45,-1.0$ for 1D,
2D, and 3D, respectively. (b) Illustration of how $\xi$ is computed: plots of
the correlation function at different values of $\Delta$, in 1D. (c) Plot of
$|\theta_{A}^{GS}-\theta_{B}^{GS}|$ as a function of $\Delta-\Delta_{c}$,
showing scaling with the critical exponent $\beta=1/2$.

###  IV.1 Correlation length from two-point functions

The energy optimization presented above only requires the calculation of one-
point functions. However, as described in Appendix A.2, the perturbative
expansion can be employed to compute also two-point functions, and can be in
principle extended to the calculation of any $n$-point function. Here, we
focus on the density-density connected two-point functions $f(|i-j|)=\langle
n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle n_{j}\rangle$, which we compute
on the optimal variational ground state for various values of $\Delta$ and
$V=0$. We plot a few examples of $f(|i-j|)$ for 1D in Fig. 8 (b). By fitting
the exponential decay via $f(|i-j|)=A\exp(-|i-j|/\xi)$ we can extract the
correlation length $\xi$, which we plot as a function of $\Delta$ (for the
case $V=0$) in Fig. 8 (a). As expected from a mean-field ansatz, the
correlation length is always finite. Note however, that our ansatz is a tensor
network with bond dimension 2; such states can be gapless in $D>1$ [50], but
we find that our variational manifold only includes gapped states.

##  V Variational dynamics

We now use our variational manifold combined with the time-dependent
variational principle (TDVP) [30, 31] to study the time evolution of the PXP
model Eq. (1) with $\Omega=1$, and show that it captures the periodic revivals
characteristic of quantum scarred systems. The 1D results were already
obtained in Ref. [27] and are presented here for comparison with the $D>1$
case.

###  V.1 TDVP equations of motion

In general, the variational manifold is parameterized by the variables
$\vec{\theta}$ and $\vec{\phi}$. However, since we are mainly concerned with
time evolution from the initial states
$\ket{\mathbb{Z}_{2}}=\ket{\uparrow}^{\otimes A}\ket{\downarrow}^{\otimes B}$
and $\ket{\mathbb{Z}_{2}^{\prime}}=\ket{\downarrow}^{\otimes
A}\ket{\uparrow}^{\otimes B}$, it suffices to consider unit cell
translationally invariant states. Also, since the TDVP is energy conserving
[30, 31], it can be shown that $d\vec{\phi}/dt=0$ when starting from the
$\ket{\mathbb{Z}_{2}}$ states, and we can set all $\vec{\phi}=0$ (see Ref.
[27] and Appendix C.3.2). Thus our variational manifold is two dimensional,
parameterized by the position vector $\vec{\theta}=(\theta_{A},\theta_{B})$.
We calculate the dynamics within this variational space by minimizing the
leakage rate, or “quantum leakage”, given by
$\Gamma^{2}=\langle\delta|\delta\rangle$, where $\ket{\delta}$ is the
component of the exact time evolution that is orthogonal to the variational
space.

The application of the TDVP to this variational manifold (described in
Appendix C), results in equations of motion
$\dot{\vec{\theta}}=(\dot{\theta}_{A},\dot{\theta}_{B})$ that, due to
symmetry, may be expressed in the form
$\dot{\theta}_{A}=f(\theta_{A},\theta_{B})$ and
$\dot{\theta}_{B}=f(\theta_{B},\theta_{A})$, where $f(\theta_{A},\theta_{B})$
can be calculated by exact tensor network contraction in $D\leq 2$ and via
perturbative expansion in any $D$.

Along the axes of the $(\theta_{A},\theta_{B})$ space, the equations of motion
take on the following, particularly simple form (derived in Appendix F)

$\displaystyle\dot{\vec{\theta}}(\theta_{A},0)$
$\displaystyle=2(\hat{\theta}_{A}+\cos^{D}\\!(\theta_{A}/2)\hat{\theta}_{B})$
(10a) $\displaystyle\dot{\vec{\theta}}(0,\theta_{B})$
$\displaystyle=2(\cos^{D}\\!(\theta_{B}/2)\hat{\theta}_{A}+\hat{\theta}_{B}),$
(10b)

where $\hat{\theta}_{A}$ and $\hat{\theta}_{B}$ are the unit vectors. In 1D
there also exists a closed-form expression for the equations of motion in
general; it was shown [27] that these tensor network calculations can be
performed analytically to show that
$f(\theta_{A},\theta_{B})=2(\cos(\theta_{B}/2)+\sin(\theta_{A}/2)\cos(\theta_{A}/2)\tan(\theta_{B}/2))$.

Figure 9: Flow diagrams and paths for 1D (left) and 2D (right). The heatmap is
the leakage rate $\gamma=\sqrt{\innerproduct{\delta}{\delta}/N}$; note that
the leakage rate is exactly zero along the $x$ and $y$ axes.

###  V.2 Periodic trajectories

In the TDVP framework, the interacting spins problem is approximated by a two
dimensional nonlinear dynamical system with equations of motion
$d\vec{\theta}/dt=\dot{\vec{\theta}}(\theta_{A},\theta_{B})$. These equations
can be integrated numerically to obtain trajectories starting from any initial
condition, and the vector field $\dot{\vec{\theta}}(\theta_{A},\theta_{B})$
can be straightforwardly analyzed to infer the stability of these
trajectories. In Figs. 9-10 we plot
$\dot{\vec{\theta}}(\theta_{A},\theta_{B})$ and the leakage rate in the
$(\theta_{A},\theta_{B})$ plane of the exact equations of motion obtained
analytically in Ref. [27] in 1D, on finite cylinders via exact contraction of
the tensor network in 2D (Fig. 9), and the equation of motions obtained from
the perturbative expansion in 2D and 3D (Fig. 10). We have marked the
$\ket{\mathbb{Z}_{2}}$ and $\ket{\mathbb{Z}_{2}^{\prime}}$ states; the all
spin down state $\ket{\downarrow}^{\otimes N}$ is at the origin. Note that all
points on the $x$ and $y$ axes correspond to product states. For example,
points on the $x$ axis are states of the form $\bigotimes_{i\in
A}(\cos(\theta_{A}/2)\ket{\downarrow}_{i}-i\sin(\theta_{A}/2)\ket{\uparrow}_{i})\otimes(\bigotimes_{j\in
B}\ket{\downarrow}_{j})$.

In each case (1D, 2D, and 3D), we observe a periodic path connecting the
$\ket{\mathbb{Z}_{2}}$ and $\ket{\mathbb{Z}_{2}^{\prime}}$ states; these non-
ergodic trajectories correspond to the revival behaviour observed in the
quantum spin system and it is easy to see from the vector field plot of
$\dot{\vec{\theta}}(\theta_{A},\theta_{B})$ that they are unstable orbits
alike the one-dimensional case. The path periods, i.e. the time it takes to go
$\ket{\mathbb{Z}_{2}}\rightarrow\ket{\mathbb{Z}_{2}^{\prime}}\rightarrow\ket{\mathbb{Z}_{2}}$,
are $T_{\mathrm{TDVP}}=4.820,5.168,5.345$ for 1D, 2D, 3D, respectively. These
have to be compared to the revival times
$T_{\mathrm{exact}}=4.786,5.154,5.340$ that we obtained from the exact
calculation of the many-body dynamics on periodic systems of up to $36$, $48$,
and $64$ sites in 1D, 2D, and 3D, respectively. Despite these numbers being
extracted from small systems, they exhibit no noticeable finite-size effects
(see Appendix G for details) and they demonstrate that the TDVP combined with
our variational ansatz provide periods that improve with increasing $D$
($T_{\mathrm{TDVP}}-T_{\mathrm{Exact}}=0.035,0.014,0.005$ for $D=1,2,3$).

Figure 10: Flow diagrams and paths for 2D (left) and 3D (right), from the
perturbative method. Note that the corners have been removed: in these regions
the perturbative expansion behavior is uncontrolled at the order
$(\mathcal{N}=12$ (2D) and $(\mathcal{N}=9$ (3D) at which the calculation is
performed. The heatmap is the leakage rate
$\gamma=\sqrt{\innerproduct{\delta}{\delta}/N}$.

We also plot the leakage rate $\gamma=\sqrt{\innerproduct{\delta}{\delta}/N},$
normalized to be an intensive quantity. In regions where the leakage rate is
small, the variational dynamics closely approximates the true evolution of the
quantum system. We observe that the periodic paths lie within regions of low
leakage, thus giving validity to the results. We calculated the integrated
leakage rates $\oint\gamma dt\simeq 0.17$ in 1D, $\oint\gamma dt\simeq 0.16$
in 2D, and $\oint\gamma dt\simeq 0.13$ in 3D. The leakage rate is exactly zero
along the lines $\theta_{A}=0$ and $\theta_{B}=0$; we show this analytically
(in any dimension) in Appendix F.

It is nontrivial that the variational dynamics remains accurate along the
entire path. In general, a starting point in an area of low leakage may travel
to an area of high leakage, given our variational manifold of low-entanglement
states. A simple example is the path starting at the origin (i.e. with initial
state $\ket{\downarrow}^{\otimes N}$) and ending at $(\pi,\pi)$, which has
leakage $\oint\gamma dt\simeq 1.28$ and $0.46$ in 1D and 2D, respectively.

Observing the trend from one to three dimensions, we conjecture that as one
goes to higher dimensional hypercubic lattices, the path period will
monotonically increase, asymptotically approaching $2\pi$ in the infinite
dimensional limit, with the trajectory becoming more and more square shaped,
i.e. more closely following the $x$ and $y$ axes. We can give an intuitive
argument for this. Consider starting in the state
$\ket{\mathbb{Z}_{2}}=\ket{\uparrow}^{\otimes A}\ket{\downarrow}^{\otimes B}$.
Initially the $A$ spins will rotate freely while the $B$ spins are frozen. As
the dimension increases, so does the number of neighbours of each site, thus
increasing the effect of the Rydberg blockade. In the infinite dimensional
limit, the $B$ spins will remain nearly frozen until the $A$ spins are nearly
pointing down, so we expect the time of half a periodic path to approach
$\pi$.

This is corroborated by the variational equations of motion Eq. (10): in the
limit $D\rightarrow\infty$, the velocity is $2\hat{\theta}_{A}$ along the $x$
axis and $2\hat{\theta}_{B}$ along the $y$ axis. An alternative picture of the
higher-dimensional limit can be given via the related idea of PXP models on
complete bipartite graphs [51].

##  VI Conclusion

We have introduced a manifold of low-entanglement states that accurately
captures much of the essential physics of the Rydberg arrays in the blockade
regime on hypercubic lattices. Our approach provides a general way to study
such constrained systems by a variational mean-field method, which is
otherwise challenging due to the blockade, that precludes direct application
of product state ansatzes.

We studied the equilibrium properties of the generalized PXP model with the
variational ansatz in cubic lattices up to three dimensions, and found that it
predicts phase boundaries and transitions in good agreement with exact
diagonalization and previously known results.

We applied the TDVP to calculate the approximate time evolution of the PXP
model by considering only the component that lies within the variational
manifold. We not only obtained important evidence of the non-ergodic behavior
that originate from the $\mathbb{Z}_{2}$ initial state in $D>1$ PXP models,
but also demonstrated that TDVP on this manifold of states yields quantitative
predictions that improve with increasing lattice dimensionality, as validated
by leakage rate and accuracy of the revivals period $T$. Moreover, our results
allowed us to infer the infinite-dimensional limit of the PXP non-ergodic
dynamics, that consists of a periodic orbit exactly supported on the cartesian
axes in the parameter space of our variational manifold, with period $T=2\pi$.

The are several potential directions for future work. At equilibrium, our
ansatz can be generalized to larger blockade radii, to encode lower density
symmetry broken phases of Rydberg arrays. Moreover, thanks to its tensor
network representation, it can be employed to approximate low-energy
excitations of Rydberg atom Hamiltonians [52]. Out of equilibrium, the TDVP
can expended to calculate the variational dynamics for constrained systems
with a time-dependent Hamiltonian, which is especially relevant in
experimental platforms [20].

###### Acknowledgements.

We acknowledge financial support from the ERC Starting grant QARA (grant no.
101041435), the European Unions Horizon 2020 research and innovation program
under Grant Agreement No. 101079862 (PASQuanS2), and by the EU-QUANTERA
project TNiSQ (N-6001). G.G. acknowledges support from the European Union’s
Horizon Europe program under the Marie Sklodowska Curie Action TOPORYD (Grant
No. 101106005).

##  Appendix A Methods

The simple structure of the variational ansatz allows one to perform tensor
network calculations relatively easily. Here we provide the details of how
these calculations are done.

Specifically, we discuss methods to calculate
$\langle\psi|\hat{\mathcal{O}}|\psi\rangle$, where
$\ket{\psi}=|\psi(\vec{\theta},\vec{\phi})\rangle$ is the variational state
and $\hat{\mathcal{O}}$ is a local operator. We will show that it is possible
to calculate $\langle\psi|\hat{\mathcal{O}}|\psi\rangle$ by contracting a bond
dimension 2 tensor network. This allows one to perform explicit tensor network
contractions for up to relatively large quasi-2D cylinders. We also introduce,
for $D\geq 2$, a method to express the tensor network contraction as a
perturbative expansion. As a corollary, we will show that the variational
state is always normalized, i.e. $\innerproduct{\psi}{\psi}=1$.

###  A.1 Properties of the variational ansatz

For clarity, let us temporarily work in two dimensions. We define the tensor
$T$ as a contraction of the PEPS tensor $M$ with its adjoint:

$\displaystyle
T=\;\vbox{\hbox{\includegraphics[width=54.06006pt]{figs/MM2D.pdf}}}\;=\;\;\vbox{\hbox{\includegraphics[width=51.21504pt]{figs/T2D.pdf}}}\;,$
(11)

where parallel in and out legs have been combined, so the legs of $T$ have
dimension $4$. The quantity $\innerproduct{\psi}{\psi}$ is thus a contraction
of an infinite network of $T$s.

It turns out that the network of $T$s is equal to a network of bond-
dimension-$2$ tensors $t$. That is, for the purposes of calculating
$\innerproduct{\psi}{\psi}$ or $\langle\psi|\hat{\mathcal{O}}|\psi\rangle$, we
can perform a “bond dimension reduction” operation

$\displaystyle\vbox{\hbox{\includegraphics[width=51.21504pt]{figs/T2D.pdf}}}\quad\longrightarrow\quad\vbox{\hbox{\includegraphics[width=45.52458pt]{figs/tt2D.pdf}}}\,,$
(12)

where $t_{abcd}=T_{(aa)(bb)(cc)(dd)}$ for $a,b,c,d\in\\{0,1\\}$, and all other
elements of $T$ are discarded. This works because for all the nonzero elements
of $T$, the outgoing legs are $(00)$ or $(11)$. Thus when considering the
entire network contraction, we can ignore the indices $(01)$ and $(10)$.

Let us now consider the general ($D$ dimensional) case. Using the thick-line
notation introduced in Eq. (6), we write the tensor $T$ as

$\displaystyle
T=\vbox{\hbox{\includegraphics[width=54.06006pt]{figs/T_tensor.pdf}}}\;.$ (13)

Using the matrix elements of $M$ from Eq. (III), we can find the matrix
elements of $T=\innerproduct{M}{M}$, which is given by

$\displaystyle
T=\cos^{2}(\theta/2)|\vec{0}\rangle\langle\vec{0}|\otimes|\vec{0}\rangle\langle\vec{0}|$
$\displaystyle+\sin^{2}(\theta/2)|\vec{0}\rangle\langle\vec{1}|\otimes|\vec{0}\rangle\langle\vec{1}|$
$\displaystyle+\cos(\theta/2)\big{(}|\beta\rangle\langle\vec{0}|\otimes|\vec{0}\rangle\langle\vec{0}|$
$\displaystyle+|\vec{0}\rangle\langle\vec{0}|\otimes|\beta\rangle\langle\vec{0}|\big{)}$
$\displaystyle+|\beta\rangle\langle\vec{0}|\otimes|\beta\rangle\langle\vec{0}|,$
(14)

where the $\ket{\alpha},\ket{\beta}$ notation is the same as originally
introduced for Eq. (III). The exact meaning of Eq. (A.1) is as follows: $T$ is
a tensor with $D$ incoming legs and $D$ outgoing legs, and each leg has
dimension 4. Let $a,b,c,d$ each be a string of ones and zeros of length $D$.
When the term
$x|\vec{a}\rangle\langle\vec{b}|\otimes|\vec{c}\rangle\langle\vec{d}|$ shows
up in the right hand side of Eq. (A.1), it means that the element of $T$ with
incoming indices set to $(a_{1}c_{1}),(a_{2}c_{2}),\dots,(a_{D}c_{D})$ and
outgoing indices set to $(b_{1}d_{1}),(b_{2}d_{2}),\dots,(b_{D}d_{D})$ has the
value $x$.

Like in the 2D case, we can perform a bond dimension reduction, and then
$\innerproduct{\psi}{\psi}$ is equal to the contraction of a network of bond-
dimension-2 tensors $t(\theta)$, which we express diagrammatically as

$\displaystyle
t_{(\text{in})(\text{out})}=\vbox{\hbox{\includegraphics[width=79.6678pt]{figs/t_reduced.pdf}}}$
(15)

and are given by

$\displaystyle
t=\cos^{2}(\theta/2)|\vec{0}\rangle\langle\vec{0}|+\sin^{2}(\theta/2)|\vec{0}\rangle\langle\vec{1}|+\ket{\beta}\langle\vec{0}|$
(16)

We can write $t$ as $t(\theta)=p-\sin^{2}(\theta/2)q$, where the constant
tensors $p$ and $q$ are

$\displaystyle p=\ket{\alpha}\langle\vec{0}|$ (17a) $\displaystyle
q=|\vec{0}\rangle\langle\vec{0}|-|\vec{0}\rangle\langle\vec{1}|$ (17b)

The decomposition $t=p-\sin^{2}(\theta/2)q$ turns out to be very useful in
performing calculations. Observe that in the network of $t$s, following a
directed line never results in a closed loop. (In the 2D square lattice, for
example, this follows trivially if we define all the arrows to point
rightwards or downwards.)

A simple but important result (which follows directly from Eq. (17)) is the
fact when $q$ is contracted with a $p$ on each of its outgoing legs, the
result is a tensor with all elements equal to zero. That is,

$\displaystyle\vbox{\hbox{\includegraphics[width=93.89418pt]{figs/qnContraction.pdf}}}=0.$
(18)

Thus, each $q$ must be attached via an outgoing leg to at least one other $q$;
otherwise the contraction of the network will be zero.

####  A.1.1 Normalization

An immediate corollary is that the variational state is normalized in the
thermodynamic limit, i.e. that $\innerproduct{\psi}{\psi}=1$. There is one
additional condition that we must impose: we must assume that at most a finite
number of the $\theta_{i}$s in $\vec{\theta}$ are integer multiples of $\pi$.
That is, we assume that $\sin(\theta_{i}/2)\neq\pm 1$ for all but a finite
number of $\theta_{i}$s.

In the expansion of $\innerproduct{\psi}{\psi}$ using
$t=p-\sin^{2}(\theta/2)q$, there are no terms with a finite number of $q$s due
to Eq. (18): there are either no $q$s, or at least an infinite number of them,
since each $q$ must attach to at least one other $q$, ad infinitum. But the
tensor networks with an infinite number of $q$s give zero contribution,
because each $q$ comes with a factor of $\sin^{2}(\theta/2)$. Thus
$\innerproduct{\psi}{\psi}=(\text{network of all }p\text{s})=1$.

####  A.1.2 Calculation by explicit contraction

In $D\leq 2$, the most straightforward way to calculate
$\langle\psi|\hat{\mathcal{O}}|\psi\rangle$ is to explicitly contract the
tensor network. In 1D this can be done for arbitrary system sizes, and in 2D,
we can define the system on a cylinder. Due to the reduced effective bond
dimension, we only need to work with a $2^{L}\times 2^{L}$ transfer matrix
(rather than $4^{L}\times 4^{L}$), where $L$ is the circumference of the
cylinder; in other words, the computational cost is equivalent to performing
exact diagonalization on a 1D spin chain of length $L$.

Due to the limited ability of the ansatz to encode long range correlations,
the results do not depend strongly on the size of the cylinder. We performed
numerical calculations on $L=10$ cylinders, but we verified that increasing
the circumference to $L=12$ or $L=14$ leaves all results unchanged well beyond
the relevant accuracy of a few significant figures.

###  A.2 Perturbative expansion

Using the results of the previous section, we show how one can express
$\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$ as an infinite series.

In the interest of presenting the ideas here compactly, we write things down
as algebraic equations rather than diagrams, but at the expense of using
somewhat schematic notation. Here, $\Tr(\cdots)$ means a contraction over the
infinite network, and when tensors inside a $\Tr(\cdots)$ are “multiplied”,
they are being placed on lattice sites and then contracted. In the one
dimensional case, the tensors are matrices, $\Tr(\cdots)$ is indeed a trace,
and the equations here may be taken literally.

Recall that $\innerproduct{\psi}{\psi}$ is the contraction of a network of $t$
tensors. Let $\Theta$ be the tensor that results when $\hat{\mathcal{O}}$ is
sandwiched with the PEPS tensors: then $\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$
is the contraction of a network with $\Theta$ on a local region and $t$ on all
other sites. That is,

$\displaystyle\Theta=\Big{(}\prod_{i\in
I}M(\theta_{i})^{*}\Big{)}\hat{\mathcal{O}}\Big{(}\prod_{i\in
I}M(\theta_{i})\Big{)},$ (19)
$\displaystyle\bra{\psi}\hat{\mathcal{O}}\ket{\psi}=\Tr\Big(\Theta\prod_{i\notin
I}t(\theta_{i})\Big{missing}),$ (20)

where $I$ is the region (i.e. set of sites) that $\Theta$ covers. Note that
due to the mechanics of the bond dimension reduction operation, $I$ may be
slightly larger than the region that $\hat{\mathcal{O}}$ itself acts on.

Using $t(\theta_{i})=p-\sin^{2}(\theta_{i}/2)q$, we can expand:

$\displaystyle\bra{\psi}\hat{\mathcal{O}}\ket{\psi}=\Tr\Big(\Theta_{I}\prod\limits_{i\notin
I}(p_{i}-\sin^{2}(\theta_{i}/2)q_{i})\Big{missing})$ (21)
$\displaystyle=\sum\limits_{J}\Tr[\prod\limits_{j\in
J}(-\sin^{2}(\theta_{j}/2)q_{j})\Big{(}\prod\limits_{i\notin J\cup
I}p_{i}\Big{)}\Theta_{I}\Big{]}$ (22)
$\displaystyle=\sum\limits_{J}\Bigg{(}(-1)^{|J|}\Big{(}\prod\limits_{j\in
J}\sin^{2}(\theta_{j}/2)\Big{)}\Tr[\Theta_{I}\prod\limits_{\begin{subarray}{c}j\in
J\\\ i\notin J\cup I\end{subarray}}q_{j}p_{i}\Big{]}\Bigg{)}.$ (23)

We have expressed $\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$ as a sum over
regions $J$. Note that $p$ and $q$ are constant tensors: the subscript denotes
which site they are on. Likewise we wrote $\Theta_{I}$ to emphasize that it
lies on the region $I$. $|J|$ is the number of sites in $J$.

The key point now is that $\Tr(\Theta\prod q_{i}p_{j})$ will be zero, due to
Eq. (18), unless the region $J$ is “anchored” by $\Theta$. When
$\Tr(\Theta\prod q_{i}p_{j})$ is nonzero, it can usually be computed quickly
by reading off the tensor elements of $\Theta$.

The perturbative expansion is thus a sum over allowed regions $J$. This is the
general idea; the numerical implementation is described below.

###  A.3 Numerical implementation

Now let us discuss the practical implementation. Variants of the calculation
we discuss here are used for both the ground state and the TDVP calculations.

Let us now restrict ourselves to the unit cell translation invariant case,
where the variational ansatz becomes
$\ket{\psi(\theta_{A},\theta_{B},\phi_{A},\phi_{B})}$. Using Eq. (23), we can
express $\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$ as

$\displaystyle\bra{\psi}\hat{\mathcal{O}}\ket{\psi}=h(\theta_{A},\theta_{B},\phi_{A},\phi_{B})\Sigma(f,\theta_{A},\theta_{B}),$
(24)

where the prefactor $h(\theta_{A},\theta_{B},\phi_{A},\phi_{B})$ is some
function that comes from the aforementioned $\Tr(\Theta\prod q_{i}p_{j})$, and
is usually not hard to evaluate, and where $\Sigma(f,\theta_{A},\theta_{B})$
is the perturbative expansion:

$\displaystyle\Sigma(f,\theta_{A},\theta_{B})=\sum_{n,m=0}^{\infty}\\!(-1)^{n+m}f_{n,m}\sin^{2n}\\!\\!\Big{(}\frac{\theta_{A}}{2}\Big{)}\sin^{2m}\\!\\!\Big{(}\frac{\theta_{B}}{2}\Big{)},$
(25)

where $f$ is a matrix of “counting factors” that depends on
$\hat{\mathcal{O}}$ and is obtained by counting the ways that the $q_{i}$s may
be arranged in $\Tr(\Theta\prod q_{i}p_{j})$. To see how exactly Eq.(25)
arises, it is instructive to work out the calculation for a specific example.

####  A.3.1 Example: $\bra{\psi}n_{i}\ket{\psi}$

Consider the case where $\hat{\mathcal{O}}$ is the number operator $n_{i}$,
for some site $i\in A$. Note that since the state $\ket{\psi}$ lies inside the
constrained subspace, this is equivalent to calculating
$\bra{\psi}\mathcal{P}n_{i}\mathcal{P}\ket{\psi}$. In this case, the tensor
$\Theta$ occupies only a single site, and the prefactor $h$ is simply
$h=h(\theta_{A})=\sin^{2}(\theta_{A}/2)$.

Let us assume that we have a 2D square lattice with all arrows pointing
rightwards or downwards, though the general idea is the same for other
lattices. (In particular, Eq. (25), as written, applies for any lattice, as
long as a suitable version of $f_{n,m}$ is used.)

We want to calculate the contraction of a tensor network with $\Theta$ on one
site, and $t$ on all other sites. We use the fact that
$t=p-\sin^{2}(\theta/2)q$, and then expand in powers of $\sin^{2}(\theta/2)$.
That is, we want to express the result of the contraction of the tensor
network as sum = ($\Theta$ with $p$ on all other sites) + (configurations with
one $q$) + (configurations with two $q$s) + ….

The allowed configurations are constrained by where the $q$s can be placed:
each $q$ must be above or to the left of another $q$, or $\Theta$. Also, note
that each $q$ comes with a factor of $-\sin^{2}(\theta/2)$. We can write the
expansion as in Eq. (25), where, in each term in the sum, $n$ is the number of
$q$s on sublattice $B$ and $m$ is the number of $q$s on sublattice $A$, and
where $f_{n,m}$ is a coefficient that can be obtained by counting.

Specifically, $f_{n,m}$ is defined as: the number of possible configurations
on a bipartite lattice (with the origin on sublattice $A$) with: $n$
$q$-tensors on $B$, $m$ $q$-tensors on $A$, a single $\Theta$ tensor at the
origin, and a $p$-tensor on all the other sites; with the constraint that each
$q$ must have another $q$ (or the $\Theta$) immediately “downstream” of it
(i.e. to the right of and/or below).

The $f$ matrix elements are obtained via brute force counting using a
computer. For a 2D square lattice, the first few values of $f_{n,m}$ are
(where the top left element of the matrix is $f_{0,0}$):

$\displaystyle f^{(2D)}=\begin{pmatrix}1&0&0&0&0&0&0&\dots\;\\\
2&4&2&0&0&0&0\\\ 1&11&25&21&6&0&0\\\ 0&10&72&174&192&100&20\\\
0&3&87&510&1281&1680&\\\ 0&0&48&732&3780&&\\\ 0&0&10&560&&&&\\\
&&\vdots&&&&&\ddots\end{pmatrix}.$ (26)

The perturbative expansion only is valid in the region where $\theta_{A}$ or
$\theta_{B}$ or both are small enough that
$f_{n,m}(\sin^{2}(\theta_{A}/2))^{n}(\sin^{2}(\theta_{B}/2))^{m}$ decreases
with larger values of $n$ and $m$.

####  A.3.2 Two point functions

The method described here also allows for the calculation of two point
functions, which in practice corresponds to calculating
$\bra{\psi}\hat{\mathcal{O}}\ket{\psi}$ where $\hat{\mathcal{O}}$ is a non-
local operator consisting of two parts. The simplest example would be
$\bra{\psi}n_{i}n_{j}\ket{\psi}$, which may be calculated by a generalization
of the procedure described above for $\bra{\psi}n_{i}\ket{\psi}$.

###  A.4 Asymptotic nature of the expansion

Here we continue the discussion on the perturbative and non-perturbative
regimes from Sec. III.1. Since the exact form of the series expansion depends
on the operator being calculated, it is difficult to discuss this in a
completely general way. We use the simplest case where the operator is
$\hat{\mathcal{O}}=n_{i}$; that is, we use $f_{n,m}=f^{(2D)}_{n,m}$ from Eq.
(26). However, the results are valid for any local operator
$\hat{\mathcal{O}}$ as the scaling of $f_{n,m}$ with $n$ and $m$ would be
similar.

We now analyze the quantitative behaviour of the asymptotic series, Eq. (25).
We define the order of a truncated expansion to be
$\mathcal{N}\equiv\max(n+m)$. Let $s_{n,m}$ refer to a term in the sum, and
let $S_{\mathcal{N}}\equiv\sum_{n+m\leq\mathcal{N}}s_{n,m}$ be the partial sum
up to order $\mathcal{N}$. To probe the accuracy of the expansion at a certain
order, we consider the quantity $|S_{\mathcal{N}}-S_{\mathcal{N}-1}|$. In Fig.
11 we show $|S_{\mathcal{N}}-S_{\mathcal{N}-1}|$ for the highest orders of the
expansion as used in the calculations in this paper, $\mathcal{N}=12$ in 2D
and $\mathcal{N}=9$ in 3D. The darker region of this plot therefore
corresponds to the perturbative regime. In this work we define it as the
region where $|S_{\mathcal{N}}-S_{\mathcal{N}-1}|<10^{-3}$ at order
$\mathcal{N}=12$ (2D) and $9$ (3D).

The boundaries of the perturbative regimes in Fig. 5 in the main text are the
points $|S_{\mathcal{N}}-S_{\mathcal{N}-1}|=10^{-3}$ from the same data as in
Fig. 11.

Figure 11: Plot of $|S_{\mathcal{N}}-S_{\mathcal{N}-1}|$ in log scale, across
the parameter space, for $\mathcal{N}=12$ in 2D (left) and $\mathcal{N}=9$ in
3D (right). Note that the colormap scale has been truncated in the sense that
numbers $<10^{-5}$ have been neglected.

####  A.4.1 Proof of non-convergence

Here we prove that the series expansion Eq. (25) divertges. In Section III.1
it was claimed that the counting factors $f_{n,m}$ scale super exponentially,
i.e. that there do not exist any $a,b\in\mathbb{R}$ such that
$f_{n,m}<a^{n}b^{m}$ for all $n,m\in\mathbb{N}$. This would imply that the sum
(25) diverges everywhere. Here, we prove that for any $N\in\mathbb{N}$ there
exists $n,m>N$ such that $f_{n,m}>\sqrt{(n+m-1)!/(n+m)}$, which is a strong
enough condition to imply super-exponential scaling. We present the proof for
2D, but the generalization to 3D is straightforward.

Instead of $f_{n,m}$, it will be easier to work with $f_{k}$, where $k=n+m$.
That is, define $f_{k}\equiv\sum_{n+m=k}f_{n,m}$, so we are counting all
allowed configurations with $k$ sites occupied. Let us state the definition of
$f_{k}$ with the aid of a picture. $f_{k}$ is the number of ways of arranging
$k$ dots on the vertices of the grid shown in Fig. 12, where the first dot is
placed at the bottom point, and then each additional dot is placed with the
condition that at least one of the two sites below it must be occupied. Now
consider how $f_{k+1}$ is related to $f_{k}$. The number of ways to place the
$(k+1)$th dot depends on the configuration of the $k$ dots already placed: in
particular, it is equal to the “perimeter” of the region occupied by the $k$
dots. This is clearly smallest when the $k$ dots are arranged compactly, such
as in the $k=9$ example shown in Fig. 12, and in this case the perimeter is
$2\sqrt{k}$. Thus $f_{k+1}>\sqrt{k}f_{k}$, and by induction,
$f_{k}>\sqrt{k!}$.

Now let us recast the result in terms of $n$ and $m$. Here let $\bar{f}_{n,m}$
be the largest $f_{n,m}$ such that $n+m=k$. Clearly $k\bar{f}_{n,m}>f_{k}$.
Then $f_{k}>\sqrt{k!}$ implies that $(n+m)\bar{f}_{n,m}>\sqrt{(n+m)!}$, which
is the claimed result.

Figure 12: A visualization of the counting problem: one counts the number of
ways of arranging $k$ red dots on this grid, subject only to the constraint
that for each dot, at least one of the two sites immediately below it must be
occupied (for sites on the edge, the condition is that the one site below it
must be occupied). Shown here is one allowed configuration for $k=9$.

##  Appendix B Energy calculation

Here we elaborate on the specifics of the tensor network calculations for the
variational energy $\bra{\psi}H\ket{\psi}$. Let us consider the parts of the
Hamiltonian (2), $H=H_{\text{PXP}}+H_{\Delta}+H_{\text{NNN}}$ one by one.

From these calculations we obtain the energy $E(\theta_{A},\theta_{B})$,
plotted in Fig. 13 and Fig. 14, and find the variational ground state by
finding the $\theta_{A},\theta_{B}$ that minimizes $E(\theta_{A},\theta_{B})$.

###  B.1 PXP term

Let us first consider $E_{\text{PXP}}=\bra{\psi}H_{\text{PXP}}\ket{\psi}$.
Using $\mathcal{P}\ket{\psi}=\ket{\psi}$, we have

$\displaystyle E_{\text{PXP}}=\sum_{i}\bra{\psi}\sigma^{x}_{i}\ket{\psi}$
$\displaystyle=\frac{N}{2}\big{(}\bra{\psi}\sigma^{x}_{a}\ket{\psi}+\bra{\psi}\sigma^{x}_{b}\ket{\psi}\big{)},$
(27)

where here (and in the sections below) we introduce the convention that $a$
and $b$ each refer to any single site on the (respectively) $A$ and $B$
sublattices.

To calculate $\bra{\psi}\sigma^{x}_{i}\ket{\psi}$, we must consider the
following tensor:

$\displaystyle\qquad\qquad\qquad
S(\theta,\phi)\equiv\vbox{\hbox{\includegraphics[width=56.9055pt]{figs/S_tensor.pdf}}}$
$\displaystyle=ie^{i\phi}\sin(\theta/2)\cos(\theta/2)(|\vec{0}\rangle\langle\vec{0}|\otimes|\vec{1}\rangle\langle\vec{0}|-|\vec{0}\rangle\langle\vec{0}|\otimes|\vec{0}\rangle\langle\vec{1}|)$
$\displaystyle=ie^{i\phi}\sin(\theta/2)\cos(\theta/2)(|\vec{0}\rangle\langle\vec{10}|-|\vec{0}\rangle\langle\vec{01}|).$
(28)

Note that in diagrams like the above, each in (out) leg on the top is paired
with an in (out) leg on the bottom, such that for a $D$ dimensional lattice we
have $D$ in-pairs and $D$ out-pairs. Since this tensor $S$ has outgoing
indices of the form $\vec{01}$, we cannot perform the bond dimension reduction
described above. The solution is to consider the object as a single unit, a
tensor occupying $D+1$ sites, on which we can perform the bond dimension
reduction.

The quantity $\bra{\psi}\sigma^{x}_{a}\ket{\psi}$ is the contraction of the
tensor network with $S(\theta_{A},\phi_{A})$ on a single site, and with
$T(\theta_{A})$ or $T(\theta_{B})$ on all other sites: the calculation may be
performed via either explicit contraction ($D\leq 2$) or the perturbative
method ($D\geq 2$). Upon performing this calculation, we find that the
$\theta,\phi$ dependence decouples in a way that we may write
$\bra{\psi}\sigma^{x}_{a}\ket{\psi}=\sin\phi_{A}\,\mathbb{F}(\theta_{A},\theta_{B}),$
and thus for the energy we have

$\displaystyle E_{\text{PXP}}(\theta_{A},\theta_{B},\phi_{A},\phi_{B})=$ (30)
$\displaystyle\frac{N}{2}\big{(}\sin\phi_{A}\,\mathbb{F}(\theta_{A},\theta_{B})+\sin\phi_{B}\,\mathbb{F}(\theta_{B},\theta_{A})\big{)},$
(31)

It is easy to check that due to the symmetry of
$\mathbb{F}(\theta_{A},\theta_{B})$, one can without loss of generality set
$\phi_{A}=\phi_{B}=\pi/2$.

###  B.2 Detuning and next-nearest-neighbour terms

For the detuning and NNN terms, we have, similar to the above:

$\displaystyle E_{\Delta}$
$\displaystyle=-\Delta\sum_{i}\bra{\psi}n_{i}\ket{\psi}$
$\displaystyle=-\frac{\Delta
N}{2}\big{(}\bra{\psi}n_{a}\ket{\psi}+\bra{\psi}n_{b}\ket{\psi}\big{)}$ (32)
$\displaystyle E_{\text{NNN}}$ $\displaystyle=V\sum_{\langle\\!\langle
ij\rangle\\!\rangle}\bra{\psi}n_{i}n_{j}\ket{\psi}$ (33)

We thus need to consider the tensor:

$\displaystyle
J(\theta)\equiv\vbox{\hbox{\includegraphics[width=56.9055pt]{figs/J_tensor.pdf}}}=\sin^{2}(\theta/2)|\vec{0}\rangle\langle\vec{1}|.$
(34)

It follows that the quantity $\bra{\psi}n_{i}\ket{\psi}$ is the contraction of
a tensor network with $J(\theta_{i})$ on a single site and $T$ on all other
sites, and that $\bra{\psi}n_{i}n_{j}\ket{\psi}$ would be the contraction of a
tensor network with $J(\theta_{i})$ on site $i$ and $J(\theta_{j})$ on site
$j$ and $T$ on all other sites.

Figure 13: Plots of the energy $E(\theta_{A},\theta_{B})$. Top row: 1D, at
$\Delta=0$ (left) and $\Delta=1$ (right). Bottom row: 2D, at $\Delta=-1$
(left) and $\Delta=0$ (right). Thus the plots on the left are for
$\Delta<\Delta_{c}$ and on the right for $\Delta>\Delta_{c}$. In each plot,
the red dot(s) indicate the minimum of $E(\theta_{A},\theta_{B})$. The line
$\theta_{A}=\theta_{B}$ is included as a visual aid. The red dot is at the
origin at $\Delta=-\infty$. As $\Delta$ is increased, it travels down the
line, until it splits into two at $\Delta=\Delta_{c}$ when the variational
ground state becomes degenerate.

Figure 14: plots of the energy $E(\theta_{A},\theta_{B})$ for 3D. at
$\Delta=-1$ (left) and $\Delta=-0.6$ (right). Note that we zoomed in on a
smaller window near the origin. The red dot is the minimum.

##  Appendix C TDVP calculations

Here we present the technical details of the variational dynamics
calculations. We first give a brief but self-contained review of the TDVP.

###  C.1 TDVP

The TDVP (see Fig. 3 for a visualization in the context of quantum scarred
dynamics) is a way of calculating the approximate time evolution of a quantum
state by projecting the true unitary time evolution onto a variational
manifold, which in our case is the family of tensor network states introduced
in Sec. III. As explained in the main text, we may impose unit cell
translation invariance, i.e. have all spins on each sublattice be in the same
state, and due to energy conservation we may set $\vec{\phi}=0$, so that
ultimately our variational manifold is parameterized by only two variables,
$\theta_{A}$ and $\theta_{B}$.

Let the variational state be $\ket{\psi(\theta_{A},\theta_{B})}$. Define the
vector $\ket{\delta}$ as the difference between the time evolution within the
manifold, and the true evolution within the full Hilbert space. The latter is
given by the Schrodinger equation, $i\frac{d\ket{\psi}}{dt}=H\ket{\psi}$. In
the following we use Greek letters as placeholders for the sub-lattices $A$
and $B$, and use the notation
$\ket{\partial_{\mu}\psi}\equiv\frac{\partial\ket{\psi}}{\partial\theta_{\mu}}$.
The leakage rate vector is given by

$\displaystyle\ket{\delta}=\sum_{\mu}\ket{\partial_{\mu}\psi}\dot{\theta}_{\mu}+iH\ket{\psi}$
(35)

The leakage rate is defined as $\Gamma^{2}=\langle\delta|\delta\rangle$. To
derive the TDVP equations of motion, we minimize the leakage rate with respect
to $\dot{\theta}_{A}$ and $\dot{\theta}_{B}$. That is, we have
$\partial\Gamma^{2}/\partial\dot{\theta}_{A}=\partial\Gamma^{2}/\partial\dot{\theta}_{B}=0$,
which gives us the equations of motion

$\displaystyle\dot{\theta}_{A}=-i\frac{\langle\partial_{A}\psi|H|\psi\rangle}{\langle\partial_{A}\psi|\partial_{A}\psi\rangle}\,,\quad\dot{\theta}_{B}=-i\frac{\langle\partial_{B}\psi|H|\psi\rangle}{\langle\partial_{B}\psi|\partial_{B}\psi\rangle}.$
(36)

In deriving these equations, we used the fact that
$\langle\partial_{A}\psi|\partial_{B}\psi\rangle=0$, i.e. the Gram matrix
$G_{\mu\nu}\equiv\langle\partial_{\mu}\psi|\partial_{\nu}\psi\rangle$ is
diagonal. This is nontrivial, but we will derive this below.

After performing the calculations for the numerators and denominators of Eq.
(36), it turns out that we can write the equations of motion as
$\dot{\theta}_{A}=f(\theta_{A},\theta_{B})$ and
$\dot{\theta}_{B}=f(\theta_{B},\theta_{A})$, where

$\displaystyle
f(\theta_{A},\theta_{B})=\frac{\mathbb{K}(\theta_{A},\theta_{B})+D\mathbb{S}(\theta_{A},\theta_{B})}{\mathbb{G}(\theta_{A},\theta_{B})},$
(37)

where $D$ is the dimension of space, and $\mathbb{K}$, $\mathbb{S}$, and
$\mathbb{G}$ are each the contraction of a single (infinite) tensor network,
the technical details of which are described below.

###  C.2 Gram matrix calculation

Here we describe the calculation of the Gram matrix, i.e. quantities of the
form $\langle\partial_{\mu}\psi|\partial_{\nu}\psi\rangle\equiv G_{\mu\nu}$.
We have, from the product rule of calculus, that
$\ket{\partial_{A}\psi}=\sum_{a\in A}\ket{\partial_{a}\psi}$, so that
$\langle\partial_{\mu}\psi|\partial_{\nu}\psi\rangle=\sum_{i\in\mu,j\in\nu}\langle\partial_{i}\psi|\partial_{j}\psi\rangle$,
where $\ket{\partial_{i}\psi}$ means that the derivative is only applied on
the site $i$. In other words, $\ket{\partial_{i}\psi}$ is a tensor network
state with the tensor $\partial M(\theta_{i})$ on site $i$, and the usual
$M(\theta_{j})$ tensor on all other sites $j\neq i$. The tensor $\partial M$
is defined as

$\displaystyle\partial
M\equiv\frac{dM}{d\theta}=\vbox{\hbox{\includegraphics[width=48.36958pt]{figs/dM.pdf}}}\qquad\quad$
$\displaystyle=-\frac{i}{2}\sin(\theta/2)\ket{\downarrow}|\vec{0}\rangle\langle\vec{0}|-\frac{1}{2}\cos(\theta/2)\ket{\uparrow}|\vec{0}\rangle\langle\vec{1}|$
(38)

Define the tensor

$\displaystyle\partial
T=\vbox{\hbox{\includegraphics[width=56.9055pt]{figs/dT.pdf}}}\;.$ (39)

It turns out that, after performing the bond-dimension reduction, that
$\partial T\propto q$. This allows us to conclude that
$\langle\partial_{i}\psi|\partial_{j}\psi\rangle=0$ for all pairs of sites
$i\neq j$. This leads to the conclusion that the Gram matrix is diagonal, and
that $G_{AA}=(N/2)\langle\partial_{a}\psi|\partial_{a}\psi\rangle$, and the
analogous equation for the $B$ sublattice. Here $a$ is any single site in $A$,
and $N$ is the total number of sites ($N$ is infinite, but it always cancels
out and never shows up in any final results).

Let us define

$\displaystyle
g=\vbox{\hbox{\includegraphics[width=56.9055pt]{figs/gtensor.pdf}}}\longrightarrow\;\;\vbox{\hbox{\includegraphics[width=51.21504pt]{figs/greduced.pdf}}}\;,$
(40)

performing a bond dimension reduction as before. The elements of $g$ are given
by
$g(\theta)=(1/4)\sin^{2}(\theta/2)|\vec{0}\rangle\langle\vec{0}|+(1/4)\cos^{2}(\theta/2)|\vec{0}\rangle\langle\vec{1}|$.

Thus $\innerproduct{\partial_{a}\psi}{\partial_{a}\psi}$ is the contraction of
a tensor network with $g(\theta_{A})$ on a single site and $T$ on all other
sites. We define the aforementioned $\mathbb{G}$ as
$\mathbb{G}(\theta_{A},\theta_{B})=\innerproduct{\partial_{a}\psi}{\partial_{a}\psi}=G_{AA}/(N/2)$.
By symmetry, the contraction of a single $g(\theta_{B})$ tensor on a $B$
lattice site is
$\mathbb{G}(\theta_{B},\theta_{A})=\langle\partial_{b}\psi|\partial_{b}\psi\rangle=G_{BB}/(N/2)$.
Similar to the calculation of the energy, these tensor network calculations
can either be performed explicitly or perturbatively, using the methods
introduced above.

###  C.3 Calculation of $\langle\partial_{\mu}\psi|H|\psi\rangle$

Here we discuss the calculation of $\langle\partial_{\mu}\psi|H|\psi\rangle$,
the numerator of the equations of motion, Eq. (36). Since $\ket{\psi}$ is
inside the Rydberg-blockaded subspace we have

$\displaystyle\langle\partial_{\mu}\psi|H|\psi\rangle=\Omega\sum_{j\in\mu}\sum_{i}\langle\partial_{j}\psi|\sigma_{i}^{x}|\psi\rangle.$
(41)

We will see that $\langle\partial_{j}\psi|\sigma_{i}^{x}|\psi\rangle$ is
nonzero only if either $i=j$, or $j$ is immediately “downstream” of $i$, i.e.
in 1D, they are nearest neighbours with $i$ to the left of $j$ (and in 2D,
either to the left or above, and so on).

####  C.3.1 Calculation of $\mathbb{K}$

Let us first consider the $i=j$ contribution, which we call $\mathbb{K}$.
Define the tensor $K$ as

$\displaystyle
K=\vbox{\hbox{\includegraphics[width=62.59596pt]{figs/Ktensor.pdf}}},$ (42)

where we find that the elements are given by
$K(\theta)=(i/2)\cos^{2}(\theta/2)|\vec{0}\rangle\langle\vec{10}|+(i/2)\sin^{2}(\theta/2)|\vec{0}\rangle\langle\vec{01}|$.
The quantity
$\mathbb{K}(\theta_{A},\theta_{B})=\langle\partial_{a}\psi|\sigma_{a}^{x}|\psi\rangle$
is the contraction of a network with $K(\theta_{A})$ on one site $a\in A$ and
$T$ on all other sites. $A\leftrightarrow B$ symmetry gives the analogous
quantity
$\mathbb{K}(\theta_{B},\theta_{A})=\langle\partial_{b}\psi|\sigma_{b}^{x}|\psi\rangle$.

Similar to the calculation for $S$ for the energy in Appendix B.1, we cannot
perform the bond dimension reduction on $K$ alone, so we consider the tensor

$\displaystyle\vbox{\hbox{\includegraphics[width=93.89418pt]{figs/KTcontraction.pdf}}}.$
(43)

We then contract this tensor with $T$ on all other sites using the same
methods as for the other tensor network contractions above, to obtain
$\mathbb{K}(\theta_{A},\theta_{B})$.

####  C.3.2 Calculation of $\mathbb{S}$

Next we consider $\langle\partial_{j}\psi|\sigma_{i}^{x}|\psi\rangle$ when
$i\neq j$. To start, consider the tensor $S$ as defined in Eq. (B.1), but we
now have $\phi_{A}=\phi_{B}=0$:

$\displaystyle\qquad\qquad\qquad
S(\theta)\equiv\vbox{\hbox{\includegraphics[width=56.9055pt]{figs/S_tensor.pdf}}}$
$\displaystyle=i\sin(\theta/2)\cos(\theta/2)(|\vec{0}\rangle\langle\vec{10}|-|\vec{0}\rangle\langle\vec{01}|).$
(44)

In this case, we find that

$\displaystyle\vbox{\hbox{\includegraphics[width=93.89418pt]{figs/STcontraction.pdf}}}=0.$
(45)

Note that this implies that $\bra{\psi}\sigma^{x}_{i}\ket{\psi}=0$ which
implies that $E=\langle H\rangle=0$ when $\ket{\psi}$ is a variational state
with $\vec{\phi}=0$. (Eq. (45) is true only if all the external legs are $0$
or $1$, i.e. if we perform bond dimension reduction. This is sufficient for
our purposes here, but the subtlety will matter when calculating the leakage
rate.)

Thus the only nonzero $i\neq j$ contributions are the cases where $i$ is next
to $j$. By symmetry, we only need to consider the diagram

$\displaystyle\vbox{\hbox{\includegraphics[width=99.58464pt]{figs/SdTcontraction.pdf}}},$
(46)

multiplied by a factor of $D$. We contract this object with the $T$ tensor on
all other sites, to get $\langle\partial_{j}\psi|\sigma_{i}^{x}|\psi\rangle$.
We define
$\mathbb{S}(\theta_{A},\theta_{B})=\langle\partial_{a}\psi|\sigma_{b}^{x}|\psi\rangle$
and
$\mathbb{S}(\theta_{B},\theta_{A})=\langle\partial_{b}\psi|\sigma_{a}^{x}|\psi\rangle$.

##  Appendix D Projected product states

Consider the set of projected (i.e. Rydberg blockaded) product states:

$\displaystyle|\psi(\vec{\vartheta},\vec{\varphi})\rangle=\mathcal{P}\bigotimes_{i}\ket{\vartheta_{i},\varphi_{i}},$
(47)

where the state of each site is a spin coherent state
$\ket{\vartheta_{i},\varphi_{i}}=\cos(\vartheta_{i}/2)\ket{\downarrow}-ie^{i\varphi_{i}}\sin(\vartheta_{i}/2)\ket{\uparrow}$.
These states resemble the variational ansatz, but are not normalized.

As noted earlier, $\mathcal{P}$ is the application of the two-site operator
$\mathcal{P}_{ij}=\ket{\downarrow\downarrow}\bra{\downarrow\downarrow}+\ket{\downarrow\uparrow}\bra{\downarrow\uparrow}+\ket{\uparrow\downarrow}\bra{\uparrow\downarrow}=\text{diag}(1,1,1,0)$
to all nearest-neighbour pairs of sites. The two-site projector can be
decomposed in a “matrix product operator” manner: diagrammatically, we have

$\displaystyle\vbox{\hbox{\includegraphics[width=19.91684pt]{figs/P.pdf}}}\;\;=\;\;\vbox{\hbox{\includegraphics[width=59.75095pt]{figs/PLPR.pdf}}}\;,$
(48)

where the internal line has bond dimension 2. Specifically, $\mathcal{P}_{R}$
and $\mathcal{P}_{L}$ are given by: $\mathcal{P}_{R}^{0}=\begin{pmatrix}1&0\\\
0&1\end{pmatrix}$,
$\mathcal{P}_{R}^{1}=\mathcal{P}_{L}^{0}=\begin{pmatrix}1&0\\\
0&0\end{pmatrix}$, and $\mathcal{P}_{L}^{1}=\begin{pmatrix}0&0\\\
0&1\end{pmatrix}$, where the superscript refers to the horizontal leg.

We can use this to write the state
$|\psi(\vec{\vartheta},\vec{\varphi})\rangle$ as a tensor network. In one
dimension, the MPS tensor is

$\displaystyle\vbox{\hbox{\includegraphics[width=56.9055pt]{figs/MPS1.pdf}}}\;\;=\;\;\vbox{\hbox{\includegraphics[width=88.2037pt]{figs/MPS.pdf}}}.$
(49)

This can be easily generalized to PEPS in higher dimensions. We can read off
all the nonzero elements of the PEPS tensors. Using the same notation as we
used for the variational ansatz (A.1), the PEPS tensor at each site is given
by

$\displaystyle\mathcal{M}(\vartheta,\varphi)=\cos(\vartheta/2)\ket{\downarrow}\ket{\alpha}\langle\vec{0}|-ie^{i\varphi}\sin(\vartheta/2)\ket{\uparrow}|\vec{0}\rangle\langle\vec{1}|.$
(50)

Comparing this expression to Eq. (III) one can see that the only difference
lies in the coefficient of the component $\ket{\beta}\langle\vec{0}|$. In 1D,
a gauge transformation exists that maps the PEPS state resulting from the
tensors contraction into each other [27]. We checked numerically that this is
not the case in $D>1$. We did so by taking the exact contraction on periodic
2D systems of the PEPS in Eq. (50) and Eq. (III). For a fixed value of the
variational parameters in the former, we optimized the overlap between the two
states over the variational parameters of the latter. The failure in finding
unity overlap within numerical precision signals the non-equivalence of the
two states. This can also be inferred from the fact that the projected product
state 2D PEPS is known to exhibit a second-order phase transition along the
line $\theta_{A}=\theta_{B}$ [53], whereas this phase transition is absent in
our ansatz (III), as we verified numerically from infinite-cylinder
calculations. We stress that the main advantage of our ansatz is the fact that
is normalized in the thermodynamic limit. This property ensures a dramatic
simplification of the perturbative expansion which allowed us to compute with
high accuracy the overlaps and expectation values required for the ground
state and TDVP analysis of the PXP model in two and three dimensions.

##  Appendix E Leakage rate

Here we discuss the calculation of the leakage rate for the variational
dynamics, the rate at which the state evolving under the full unitary
evolution leaves the variational manifold at a given instant. Specifically, we
will show how the evaluation of the leakage rate reduces to quantities of the
form $\bra{\psi}\hat{\mathcal{O}}_{\text{local}}\ket{\psi}$, which can then be
calculated via a tensor network contraction similar to the calculations
described above.

The leakage rate is the magnitude of the leakage rate vector $|\delta\rangle$:

$\displaystyle\langle\delta|\delta\rangle=\sum_{\mu\nu}\dot{\theta}_{\mu}\dot{\theta}_{\nu}\langle\partial_{\mu}\psi|\partial_{\nu}\psi\rangle+i\sum_{\mu}\dot{\theta}_{\mu}\langle\partial_{\mu}\psi|H|\psi\rangle$
(51)
$\displaystyle-i\sum_{\mu}\dot{\theta}_{\mu}\langle\psi|H|\partial_{\mu}\psi\rangle+\langle\psi|H^{2}|\psi\rangle$
(52)
$\displaystyle=\sum_{\mu}\Big{[}\dot{\theta}^{2}_{\mu}\langle\partial_{\mu}\psi|\partial_{\mu}\psi\rangle+2i\dot{\theta}_{\mu}\langle\partial_{\mu}\psi|H|\psi\rangle\Big{]}$
(53) $\displaystyle+\langle\psi|H^{2}|\psi\rangle\equiv\Gamma^{2}$ (54)

where to get to the second line, we used the fact that (for the first term)
the Gram matrix is diagonal, and (for the second term) the fact that
$\langle\partial_{\mu}\psi|H|\psi\rangle$ is pure imaginary, so that
$\langle\partial_{\mu}\psi|H|\psi\rangle=-\langle\psi|H|\partial_{\mu}\psi\rangle$.

We can plug our derived expressions for $\dot{\theta}_{A}$ and
$\dot{\theta}_{B}$ into the leakage rate expression to get

$\displaystyle\Gamma^{2}=\langle\psi|H^{2}|\psi\rangle-\sum_{\mu}\dot{\theta}^{2}_{\mu}\langle\partial_{\mu}\psi|\partial_{\mu}\psi\rangle$
(55)

The new work that we have to do is to calculate
$\langle\psi|H^{2}|\psi\rangle$.

First, note that there are two equivalent ways of writing down the PXP
Hamiltonian, which we call the “global” and “local” versions:

$\displaystyle\text{global:}\quad
H=\sum_{i}\mathcal{P}\sigma^{x}_{i}\mathcal{P}$ (56)
$\displaystyle\text{local:}\quad H=\sum_{i}\bigg{[}\Big{(}\bigotimes_{\langle
ji\rangle}P_{j}\Big{)}\sigma^{x}_{i}\bigg{]}$ (57)

Here, $P=\ket{\downarrow}\bra{\downarrow}$, and $\mathcal{P}$ is the whole-
system Rydberg-blockade projection, i.e. $\mathcal{P}=\prod_{\langle
ij\rangle}\mathcal{P}_{ij}$, where $\mathcal{P}_{ij}=\text{diag}(1,1,1,0).$ In
1D, the local version can also be written as
$H=\sum_{i}P_{i-1}\sigma^{x}_{i}P_{i+1}.$ Working with the global version, let
us expand $\langle\psi|H^{2}|\psi\rangle$:

$\displaystyle\langle\psi|H^{2}|\psi\rangle=\sum_{ij}\langle\psi|(\mathcal{P}\sigma^{x}_{i}\mathcal{P})(\mathcal{P}\sigma^{x}_{j}\mathcal{P})|\psi\rangle$
(58)

At this point, it is convenient to use the local $H$ for the $i=j$ terms, and
the global $H$ for the $i\neq j$ terms. For the $i=j$ terms, we have

$\displaystyle\sum_{i}\bra{\psi}\Big{[}\Big{(}\bigotimes_{\langle
ji\rangle}P_{j}\Big{)}\sigma^{x}_{i}\Big{]}^{2}\ket{\psi}=\sum_{i}\bra{\psi}\bigotimes_{\langle
ji\rangle}P_{j}\ket{\psi}$ (59)
$\displaystyle=\frac{N}{2}\bra{\psi}\bigotimes_{\langle
ia\rangle}P_{i}\ket{\psi}+\frac{N}{2}\bra{\psi}\bigotimes_{\langle
ib\rangle}P_{i}\ket{\psi}.$ (60)

For the $i\neq j$ terms, making use of the fact that
$\mathcal{P}\ket{\psi}=\ket{\psi}$, we have

$\displaystyle\sum_{i\neq
j}\langle\psi|(\mathcal{P}\sigma^{x}_{i}\mathcal{P})(\mathcal{P}\sigma^{x}_{j}\mathcal{P})|\psi\rangle=\sum_{i\neq
j}\bra{\psi}\sigma^{x}_{i}\mathcal{P}\sigma^{x}_{j}\ket{\psi}$ (61)

If $i$ and $j$ are not next to each other, all of the $\mathcal{P}$ operators
in $\mathcal{P}$ can be commuted out (this is easiest to see by drawing
$\mathcal{P}$ as a brickwork-style “quantum circuit”, and noting that all the
bricks are diagonal and therefore commute with each other, and that the bricks
are annihilated upon hitting $\ket{\psi}$). So, for non-nearest-neighbour
sites,

$\displaystyle\bra{\psi}\sigma^{x}_{i}\mathcal{P}\sigma^{x}_{j}\ket{\psi}=\langle\psi|\sigma_{i}^{x}\sigma_{j}^{x}|\psi\rangle=0.$
(62)

where in the last step we used the fact that we know from the explicit tensor
network calculation that
$\langle\psi|\sigma_{i}^{x}\sigma_{j}^{x}|\psi\rangle$ is nonzero only for
adjacent sites. Specifically, by nonadjacent, we mean neither directly
adjacent nor diagonally adjacent. The quantity
$\langle\psi|\sigma_{a}^{x}\sigma_{a^{\prime}}^{x}|\psi\rangle$ is nonzero if
$a$ and $a^{\prime}$ are diagonally adjacent $A$ sites that are downstream of
the same $B$ site.

For quantities involving the projector, we only have to deal with the cases
where $i$ and $j$ are next to each other. In those cases, all of the
$\mathcal{P}$ operators can be commuted out, except for the one acting on $i$
and $j$, so we have
$\bra{\psi}\sigma^{x}_{i}\mathcal{P}\sigma^{x}_{j}\ket{\psi}=\bra{\psi}\sigma^{x}_{i}\mathcal{P}_{ij}\sigma^{x}_{j}\ket{\psi}$
Using symmetry, we can write the sum as

$\displaystyle\sum_{\langle
ij\rangle}\bra{\psi}\sigma^{x}_{i}\mathcal{P}\sigma^{x}_{j}\ket{\psi}=2D\sum_{i}\bra{\psi}\sigma^{x}_{i}\mathcal{P}\sigma^{x}_{i+\hat{x}}\ket{\psi}+\frac{N}{2}{D\choose
2}\langle\psi|\sigma_{a}^{x}\sigma_{a^{\prime}}^{x}|\psi\rangle+(a\leftrightarrow
b)$ (63)
$\displaystyle=\frac{N}{2}\bigg{[}2D\langle\psi|\sigma_{a}^{x}\mathcal{P}_{a,a+\hat{x}}\sigma_{a+\hat{x}}^{x}|\psi\rangle+{D\choose
2}\langle\psi|\sigma_{a}^{x}\sigma_{a^{\prime}}^{x}|\psi\rangle+(a\leftrightarrow
b)\bigg{]}$ (64)

In the first line, we gained a factor of $2D$: the factor of 2 comes from
having considered each nearest-neighbour pair only once, and the factor of $D$
comes from the symmetry between the different directions (recall the fact for
each PEPS tensor, all of the “in” legs are interchangeable, likewise for the
“out” legs). In the second line, we used translation invariance. Using
$\langle\partial_{A}\psi|\partial_{A}\psi\rangle=\frac{N}{2}\langle\partial_{a}\psi|\partial_{a}\psi\rangle$
(and the analogous equation for $B$) for the quantity
$\sum_{\mu}\dot{\theta}^{2}_{\mu}\langle\partial_{\mu}\psi|\partial_{\mu}\psi\rangle$
, we now have the full leakage rate expression:

$\displaystyle\begin{aligned}
\frac{\Gamma^{2}}{N}\equiv\gamma^{2}=\frac{1}{2}\bra{\psi}\bigotimes_{\langle
ia\rangle}P_{i}\ket{\psi}+D\Big{(}\langle\psi|\sigma_{a}^{x}\mathcal{P}_{a,a+\hat{x}}\sigma_{a+\hat{x}}^{x}|\psi\rangle\Big{)}-\frac{1}{2}\Big{(}\dot{\theta}_{A}^{2}\langle\partial_{a}\psi|\partial_{a}\psi\rangle\Big{)}+\frac{1}{2}{D\choose
2}\langle\psi|\sigma_{a}^{x}\sigma_{a^{\prime}}^{x}|\psi\rangle+(a\leftrightarrow
b),\end{aligned}$ (65)

where $a$ and $a^{\prime}$ have a specific relationship as mentioned above,
and where the $(a\leftrightarrow b)$ applies to all of the preceding terms. At
this point, each of the terms in the above equation can be calculated via a
tensor network contraction.

##  Appendix F Explicit TDVP calculations for product states

Here we analytically perform the TDVP calculations for states along the
$\theta_{A}$ and $\theta_{B}$ axes in Fig. 9 and Fig. 10. We derive the
equations of motion Eq. (10) and demonstrate that the leakage rate is zero,
i.e. that $\gamma(\theta_{A},0)=\gamma(0,\theta_{B})=0$.

Let us consider the case $\theta_{B}=0$; by symmetry, the case $\theta_{A}=0$
will be the same but with $\theta_{A}\leftrightarrow\theta_{B}$. When doing
the TDVP calculations, we work with the variational state
$\ket{\psi(\theta_{A},\theta_{B})}$, and the states
$\ket{\partial_{a}\psi(\theta_{A},\theta_{B})}$ and
$\ket{\partial_{b}\psi(\theta_{A},\theta_{B})}$. Recall that the latter states
are defined as applying a derivative $\partial\\!/\\!\partial\theta_{A}$ or
$\partial\\!/\\!\partial\theta_{B}$ to $\ket{\psi(\theta_{A},\theta_{B})}$ at
any one particular site $a\in A$ or $b\in B$. In general, these are tensor
network states; however, when $\theta_{B}=0$, these become product states,
allowing one to perform the TDVP and leakage rate calculations by hand.

The state $\ket{\psi(\theta_{A},0)}$ is the following product state:

$\displaystyle\ket{\psi(\theta_{A},0)}=\bigotimes_{i\in A,j\in
B}\ket{\theta_{A}}_{i}\ket{\downarrow}_{j}$ (66)

where
$\ket{\theta_{A}}\equiv\cos(\theta_{A}/2)\ket{\downarrow}-i\sin(\theta_{A}/2)\ket{\uparrow}$.
The state $\ket{\partial_{a}\psi(\theta_{A},0)}$ is obtained by taking the
derivative at one site; it is given by

$\displaystyle\ket{\partial_{a}\psi(\theta_{A},0)}=\ket{\theta_{A}^{\prime}}_{a}\otimes\Big{(}\bigotimes_{\begin{subarray}{c}i\in
A\backslash\\{a\\}\\\ j\in
B\end{subarray}}\ket{\theta_{A}}_{i}\ket{\downarrow}_{j}\Big{)}$ (67)

where $\ket{\theta_{A}^{\prime}}\equiv
d\\!\ket{\theta_{A}}\\!\\!/\\!d\theta_{A}=-(1/2)\sin(\theta_{A}/2)\ket{\downarrow}-(i/2)\cos(\theta_{A}/2)\ket{\uparrow}$.

The state $\ket{\partial_{b}\psi(\theta_{A},0)}$ is more tricky to write down,
as one must start with $\ket{\partial_{b}\psi(\theta_{A},\theta_{B})}$ and
then derive the product state that results when $\theta_{B}=0$. A careful
inspection yields the result:

$\displaystyle\ket{\partial_{b}\psi(\theta_{A},0)}=\Big{(}\bigotimes_{a\in
N_{L}(b)}\cos(\theta_{A}/2)\ket{\downarrow}_{a}\Big{)}\otimes\big{(}\frac{-i}{2}\ket{\uparrow}_{b}\big{)}\otimes\Big{(}\bigotimes_{a\in
N_{R}(b)}\ket{\downarrow}_{a}\Big{)}\otimes\Big{(}\bigotimes_{\begin{subarray}{c}i\in
A\backslash N(b)\\\ j\in
B\backslash\\{b\\}\end{subarray}}\ket{\theta_{A}}_{i}\ket{\downarrow}_{j}\Big{)}$
(68) $\displaystyle=\frac{-i}{2}\cos^{D}(\theta_{A}/2)\Big{(}\bigotimes_{a\in
N(b)}\ket{\downarrow}_{a}\Big{)}\otimes\ket{\uparrow}_{b}\otimes\Big{(}\bigotimes_{\begin{subarray}{c}i\in
A\backslash N(b)\\\ j\in
B\backslash\\{b\\}\end{subarray}}\ket{\theta_{A}}_{i}\ket{\downarrow}_{j}\Big{)}$
(69)

where $N_{L}(i)$ is the “left neighbourhood” of $i$, i.e. the set of sites
immediately upstream of $i$, and $N_{R}(i)$ is the “right neighbourhood” ,
i.e. the set of sites immediately downstream of $i$, and $N(i)=N_{L}(i)\cup
N_{R}(i)$ is the full neighbourhood. Note that $|N_{L}(i)|=|N_{R}(i)|=D$.

Note that $\bra{\theta_{A}}\sigma^{x}\ket{\theta_{A}}=0$,
$\innerproduct{\theta_{A}^{\prime}}{\theta_{A}}=0$,
$\bra{\theta_{A}^{\prime}}\sigma^{x}\ket{\theta_{A}}=i/2$, and
$\bra{\downarrow}\sigma^{x}\ket{\theta_{A}}=-i\sin(\theta_{A}/2)$. This
implies that $\bra{\partial_{a}\psi}\sigma^{x}_{i}\ket{\psi}=0$ for $i\neq a$
and $\bra{\partial_{a}\psi}\sigma^{x}_{a}\ket{\psi}=i/2$. Thus we find that

$\displaystyle\innerproduct{\partial_{a}\psi(\theta_{A},0)}{\partial_{a}\psi(\theta_{A},0)}=1/4$
(70)
$\displaystyle\dot{\theta}_{A}(\theta_{A},0)=-i\frac{\langle\partial_{A}\psi|H|\psi\rangle}{\langle\partial_{A}\psi|\partial_{A}\psi\rangle}=-i\frac{\langle\partial_{a}\psi|H|\psi\rangle}{\innerproduct{\partial_{a}\psi}{\partial_{a}\psi}}=2$
(71)

Now let us consider $\dot{\theta}_{B}(\theta_{A},0)$. Since
$\ket{\partial_{b}\psi}$ has the site $b$ in the state $\ket{\uparrow}$,
$\bra{\partial_{b}\psi}\sigma^{x}_{i}\ket{\psi}=0$ unless $i=b.$ Note that
$\innerproduct{0}{\theta_{A}}=\cos(\theta_{A}/2)$; thus when taking the scalar
product in $\bra{\partial_{b}\psi}\sigma^{x}_{b}\ket{\psi}$ we gain a factor
of $\cos^{2D}(\theta_{A}/2)$. We find that

$\displaystyle\innerproduct{\partial_{b}\psi(\theta_{A},0)}{\partial_{b}\psi(\theta_{A},0)}=\frac{1}{4}\cos^{2D}(\theta_{A}/2)$
(72)
$\displaystyle\bra{\partial_{b}\psi}\sigma^{x}_{b}\ket{\psi}=\frac{i}{2}\cos^{3D}(\theta_{A}/2),$
(73)

which yields

$\displaystyle\dot{\theta}_{B}(\theta_{A},0)$
$\displaystyle=-i\frac{\langle\partial_{B}\psi|H|\psi\rangle}{\innerproduct{\partial_{B}\psi}{\partial_{B}\psi}}$
(74)
$\displaystyle=-i\frac{\langle\partial_{b}\psi|H|\psi\rangle}{\innerproduct{\partial_{b}\psi}{\partial_{b}\psi}}=2\cos^{D}(\theta_{A}/2).$
(75)

Now let us calculate the leakage rate using Eq. (65). For the case that we
consider here where $\ket{\psi}$ is a product state, only the
$\bra{\psi}\bigotimes P_{i}\ket{\psi}$ terms and the
$\dot{\theta}^{2}_{\mu}\innerproduct{\partial_{\mu}\psi}{\partial_{\mu}\psi}$
terms contribute to the leakage rate; the second and fourth terms in Eq. (65)
are zero. Since $\bra{\theta_{A}}P\ket{\theta_{A}}=\cos^{2}(\theta_{A}/2)$ and
$\bra{\downarrow}P\ket{\downarrow}=1$, we have

$\displaystyle\bra{\psi(\theta_{A},0)}\bigotimes_{i\in
N(a)}P_{i}\ket{\psi(\theta_{A},0)}$ $\displaystyle=1$ (76)
$\displaystyle\bra{\psi(\theta_{A},0)}\bigotimes_{i\in
N(b)}P_{i}\ket{\psi(\theta_{A},0)}$
$\displaystyle=\big{(}\cos^{2}(\theta_{A}/2)\big{)}^{2D}$ (77)

Thus for the leakage rate at the state $\ket{\psi(\theta_{A},0)}$, we have
$2\gamma^{2}=1+\cos^{4D}(\theta_{A}/2)-1-\cos^{4D}(\theta_{A}/2)=0$.

##  Appendix G Details on exact diagonalization

Figure 15: Ground state fidelity susceptibility for $V=0$ as function of
$\Delta$ in 2D for $N=24,36,48$ sites and in 3D for 64 sites. The peak
position signal the finite-size phase transition point from the disordered to
the $\mathbb{Z}_{2}$-ordered phase. The inset shows the extrapolation of the
thermodynamic limit critical value of $\Delta$ resulting in
$\Delta_{c}\simeq-0.19$.  Figure 16: Revival fidelity employed to extract the
first period of the approximately-periodic many-body dynamics generated by the
$|\mathbb{Z}_{2}\rangle$ initial state. The upper panel shows a comparison
between the symmetric and non-symmetric dynamics starting from the cat state
$(|\mathbb{Z}_{2}\rangle+|\mathbb{Z}^{\prime}_{2}\rangle)/\sqrt{2}$ and from
the symmetry-breaking $|\mathbb{Z}_{2}\rangle$ state, in 1D for $N=18$. The
former case produces a revival at a time $t\simeq T/2$ due to the fact that
the periodic orbit from the $|\mathbb{Z}_{2}\rangle$ state back to itself goes
through the $|\mathbb{Z}^{\prime}_{2}\rangle$ state. The second revival in the
symmetric time evolution exactly corresponds to the end of the orbit. The
lower panel shows the symmetric dynamics for the largest systems diagonalized
in all dimensions. Vertical lines correspond to the revival times.

In this appendix we provide details on the exact diagonalization techniques
employed for the ground state and dynamics calculations presented in the main
text. For diagonalizing system sizes of up to 48 and 64 sites in 2D and 3D
periodic lattices, respectively, we numerically computed the Hamiltonian
directly in the sector which is invariant under the symmetry group of the
lattice. In particular, we combined translations along and reflection with
respect to the Cartesian axes. The dimensions of the neutral symmetry sectors
for the largest system sizes considered in this work are 1682382 (48 sites in
2D) and 13292545 (64 sites in 3D), which can be easily tackled with sparse
matrix techniques to obtain the lowest-energy eigenstates and to apply the
exponential of the Hamiltonian for the time evolution of the
$|\mathbb{Z}_{2}\rangle$ state.
To extract the phase boundaries depicted in Fig. 6 we computed the ground
state fidelity susceptibility

$\mathcal{F}=\frac{1-|\langle\mathrm{GS}(\Delta)|\mathrm{GS}(\Delta+d\Delta)\rangle|}{d\Delta^{2}},$
(78)

which exhibits a peak at the finite-size phase transition point. We show in
Fig. 15 the finite size scaling of $\mathcal{F}$ for $V=0$ in 2D, from which
we extrapolate the thermodynamic limit critical value $\Delta_{c}\simeq-0.19$,
and for $V=0$ in 3D for the $4\times 4\times 4$ cube of 64 sites.
For the computation of the revival periods reported in Sec. V, we need to
compute the fidelity $|\langle\mathbb{Z}_{2}|e^{-iHt}|\mathbb{Z}_{2}\rangle|$.
Although the $|\mathbb{Z}_{2}\rangle$ state breaks the lattice symmetry of the
Hamiltonian, the revival time can be extracted from the time evolution of the
cat state
$|\psi_{0}\rangle=(|\mathbb{Z}_{2}\rangle+|\mathbb{Z}^{\prime}_{2}\rangle)/\sqrt{2}$,
whose dynamics is symmetric and exhibits twice the revivals occurring in the
non-symmetric case (see Fig. 16). We also note, as can be seen by comparing
the upper and bottom panels of Fig. 16, that finite-size effects on revival
times are negligible, since the position of fidelity maximum is unchanged from
$L=18$ (upper panel) to $L=36$ (lower panel, blue line) in 1D. We verified
that the same occurs in 2D for $N\geq 16$ ($4\times 4$ square). For the 3D
revival time a cube with 4 atoms in each direction ($N=64$) is required,
otherwise the system is quasi-two-dimensional, and the resulting revival time
is the same as in 2D.

## References

  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, A rydberg quantum simulator, Nature Physics 6, 382 (2010).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Morgado and Whitlock [2021] M. Morgado and S. Whitlock, Quantum simulation and computing with rydberg-interacting qubits, AVS Quantum Science 3, 10.1116/5.0036562 (2021).
  * Wu  _et al._ [2021] X. Wu, X. Liang, Y. Tian, F. Yang, C. Chen, Y.-C. Liu, M. K. Tey, and L. You, A concise review of rydberg atom based quantum computation and quantum simulation, Chinese Physics B 30, 020305 (2021).
  * Keesling  _et al._ [2019] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum kibble–zurek mechanism and critical dynamics on a programmable rydberg simulator, Nature 568, 207 (2019).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227 (2021).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms, Nature 595, 233 (2021).
  * Chen  _et al._ [2023] C. Chen, G. Bornet, M. Bintz, G. Emperauger, L. Leclerc, V. S. Liu, P. Scholl, D. Barredo, J. Hauschild, S. Chatterjee, M. Schuler, A. M. Läuchli, M. P. Zaletel, T. Lahaye, N. Y. Yao, and A. Browaeys, Continuous symmetry breaking in a two-dimensional rydberg array, Nature 616, 691 (2023).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Verresen  _et al._ [2021a] R. Verresen, M. D. Lukin, and A. Vishwanath, Prediction of toric code topological order from rydberg blockade, Phys. Rev. X 11, 031005 (2021a).
  * Giudici  _et al._ [2022] G. Giudici, M. D. Lukin, and H. Pichler, Dynamical preparation of quantum spin liquids in rydberg atom arrays, Phys. Rev. Lett. 129, 090401 (2022).
  * Sahay  _et al._ [2023] R. Sahay, A. Vishwanath, and R. Verresen, Quantum spin puddles and lakes: Nisq-era spin liquids from non-equilibrium dynamics (2023), arXiv:2211.01381 [cond-mat.str-el] .
  * Slagle  _et al._ [2022] K. Slagle, Y. Liu, D. Aasen, H. Pichler, R. S. K. Mong, X. Chen, M. Endres, and J. Alicea, Quantum spin liquids bootstrapped from ising criticality in rydberg arrays, Phys. Rev. B 106, 115122 (2022).
  * Giudice  _et al._ [2022] G. Giudice, F. M. Surace, H. Pichler, and G. Giudici, Trimer states with ${\mathbb{z}}_{3}$ topological order in rydberg atom arrays, Phys. Rev. B 106, 195155 (2022).
  * de Léséleuc  _et al._ [2019] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected topological phase of interacting bosons with rydberg atoms, Science 365, 775 (2019).
  * de Léséleuc  _et al._ [2019] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected topological phase of interacting bosons with rydberg atoms, Science 365, 775 (2019).
  * Kim  _et al._ [2018] H. Kim, Y. Park, K. Kim, H.-S. Sim, and J. Ahn, Detailed balance of thermalization dynamics in rydberg-atom quantum simulators, Phys. Rev. Lett. 120, 180502 (2018).
  * Scholl  _et al._ [2022] P. Scholl, H. J. Williams, G. Bornet, F. Wallner, D. Barredo, L. Henriet, A. Signoles, C. Hainaut, T. Franz, S. Geier, A. Tebben, A. Salzinger, G. Zürn, T. Lahaye, M. Weidemüller, and A. Browaeys, Microwave engineering of programmable $xxz$ hamiltonians in arrays of rydberg atoms, PRX Quantum 3, 020303 (2022).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Bluvstein  _et al._ [2021] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletić, and M. D. Lukin, Controlling quantum many-body dynamics in driven rydberg atom arrays, Science 371, 1355 (2021).
  * Serbyn  _et al._ [2021] M. Serbyn, D. A. Abanin, and Z. Papić, Quantum many-body scars and weak breaking of ergodicity, Nature Physics 17, 675 (2021).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * Jaksch  _et al._ [2000] D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Fast quantum gates for neutral atoms, Phys. Rev. Lett. 85, 2208 (2000).
  * Lukin  _et al._ [2001] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole blockade and quantum information processing in mesoscopic atomic ensembles, Phys. Rev. Lett. 87, 037901 (2001).
  * Urban  _et al._ [2009] E. Urban, T. A. Johnson, T. Henage, L. Isenhower, D. D. Yavuz, T. G. Walker, and M. Saffman, Observation of rydberg blockade between two atoms, Nature Physics 5, 110 (2009).
  * Ho  _et al._ [2019] W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Periodic orbits, entanglement, and quantum many-body scars in constrained models: Matrix product state approach, Phys. Rev. Lett. 122, 040603 (2019).
  * Ljubotina  _et al._ [2022] M. Ljubotina, B. Roos, D. A. Abanin, and M. Serbyn, Optimal steering of matrix product states and quantum many-body scars, PRX Quantum 3, 030343 (2022).
  * Michailidis  _et al._ [2020a] A. A. Michailidis, C. J. Turner, Z. Papić, D. A. Abanin, and M. Serbyn, Slow quantum thermalization and many-body revivals from mixed phase space, Phys. Rev. X 10, 011055 (2020a).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Fendley  _et al._ [2004] P. Fendley, K. Sengupta, and S. Sachdev, Competing density-wave orders in a one-dimensional hard-boson model, Phys. Rev. B 69, 075106 (2004).
  * Samajdar  _et al._ [2020] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Complex density wave orders and quantum phase transitions in a model of square-lattice rydberg atom arrays, Phys. Rev. Lett. 124, 103601 (2020).
  * Note [1] The Hilbert space dimension of $N$ spins $d_{H}(N)$ satisfies the equation $d_{H}(N+1)=d_{H}(N)+d_{H}(N-1)$, which implies $d_{H}\sim\phi^{N}$ for large $N$, where $\phi$ is the golden ratio.
  * Turner  _et al._ [2018a] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Quantum scarred eigenstates in a rydberg atom chain: Entanglement, breakdown of thermalization, and stability to perturbations, Phys. Rev. B 98, 155134 (2018a).
  * Khemani  _et al._ [2019] V. Khemani, C. R. Laumann, and A. Chandran, Signatures of integrability in the dynamics of rydberg-blockaded chains, Phys. Rev. B 99, 161101 (2019).
  * Scholl  _et al._ [2023] P. Scholl, A. L. Shaw, R. B.-S. Tsai, R. Finkelstein, J. Choi, and M. Endres, Erasure conversion in a high-fidelity rydberg quantum simulator, Nature 622, 273 (2023).
  * Samajdar  _et al._ [2021] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Quantum phases of rydberg atoms on a kagome lattice, Proceedings of the National Academy of Sciences 118, e2015785118 (2021).
  * Verresen  _et al._ [2021b] R. Verresen, M. D. Lukin, and A. Vishwanath, Prediction of toric code topological order from rydberg blockade, Phys. Rev. X 11, 031005 (2021b).
  * Chepiga and Mila [2019] N. Chepiga and F. Mila, Floating phase versus chiral transition in a 1d hard-boson model, Phys. Rev. Lett. 122, 017205 (2019).
  * Samajdar  _et al._ [2018] R. Samajdar, S. Choi, H. Pichler, M. D. Lukin, and S. Sachdev, Numerical study of the chiral ${\mathbb{z}}_{3}$ quantum phase transition in one spatial dimension, Phys. Rev. A 98, 023614 (2018).
  * Giudici  _et al._ [2019] G. Giudici, A. Angelone, G. Magnifico, Z. Zeng, G. Giudice, T. Mendes-Santos, and M. Dalmonte, Diagnosing potts criticality and two-stage melting in one-dimensional hard-core boson models, Phys. Rev. B 99, 094434 (2019).
  * Rader and Läuchli [2019] M. Rader and A. M. Läuchli, Floating phases in one-dimensional rydberg ising chains (2019), arXiv:1908.02068 [cond-mat.quant-gas] .
  * Kalinowski  _et al._ [2022] M. Kalinowski, R. Samajdar, R. G. Melko, M. D. Lukin, S. Sachdev, and S. Choi, Bulk and boundary quantum phase transitions in a square rydberg atom array, Phys. Rev. B 105, 174417 (2022).
  * Turner  _et al._ [2018b] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nature Physics 14, 745 (2018b).
  * Michailidis  _et al._ [2020b] A. A. Michailidis, C. J. Turner, Z. Papić, D. A. Abanin, and M. Serbyn, Stabilizing two-dimensional quantum scars by deformation and synchronization, Phys. Rev. Res. 2, 022065 (2020b).
  * Heller [1984] E. J. Heller, Bound-state eigenfunctions of classically chaotic hamiltonian systems: Scars of periodic orbits, Phys. Rev. Lett. 53, 1515 (1984).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014).
  * Sachdev  _et al._ [2002] S. Sachdev, K. Sengupta, and S. M. Girvin, Mott insulators in strong electric fields, Phys. Rev. B 66, 075128 (2002).
  * Verstraete  _et al._ [2006] F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Criticality, the area law, and the computational power of projected entangled pair states, Phys. Rev. Lett. 96, 220601 (2006).
  * Windt and Pichler [2022] B. Windt and H. Pichler, Squeezing quantum many-body scars, Phys. Rev. Lett. 128, 090606 (2022).
  * Vanderstraeten  _et al._ [2019] L. Vanderstraeten, J. Haegeman, and F. Verstraete, Simulating excitation spectra with projected entangled-pair states, Phys. Rev. B 99, 165121 (2019).
  * Ji  _et al._ [2011] S. Ji, C. Ates, and I. Lesanovsky, Two-dimensional rydberg gases and the quantum hard-squares model, Phys. Rev. Lett. 107, 060406 (2011).