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

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][5] 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 [6][7][8][9], as well as spin liquids with topological order [10][11][12][13][14][15] and symmetry protected topological phases [16,17].They can also be used to study thermalization dynamics [18,19] and were instrumental in the experimental discovery [20,21] of the out-of-equilibrium phenomenon now known as quantum many-body scars (QMBS) [22], which have sparked interest as an example of non-thermalizing behaviour in quantum many-body systems, violating the eigenstate thermalization hypothesis [23,24].Remarkably, often the physics of these systems simply emerges from the Rydberg blockade mechanism [25][26][27] 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.
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 constraintsatisfying states that can be elegantly represented as min-imally entangled tensor networks.Our ansatz can be seen as a generalization of minimal matrix-product wavefunctions that have been employed for 1D chains [28,29] and the related tree networks in Ref. [30].
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 [31,32].In 1D, the phase diagram has been studied analytically [33], and in 1D and 2D, these phenomena have been observed experimentally [7,20,21] and studied FIG. 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 θA and θB, corresponding to the two sub-lattices.
numerically with density matrix renormalization group techniques [34].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 A. 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 |↓⟩ and |↑⟩, 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 |↑⟩ state.In physical realizations this constraint arises dynamically from strong statedependent 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 + √ 5)/2 ≃ 1.618 [35].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 P. It can be composed from local projectors P ij that act on pairs of neighboring spins i and j, i.e., P = ⟨ij⟩ P ij .Here P ij = |↓↓⟩ ⟨↓↓| + |↓↑⟩ ⟨↓↑| + |↑↓⟩ ⟨↑↓| simply annihilates components that do not satisfy the constraint on the link between i and j.

B. 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 where σ x = |↑⟩ ⟨↓| + |↓⟩ ⟨↑|.The Hamiltonian (1) can be interpreted as the projection of a non-interacting model of independently rotating spins into the constraintsatisfying subspace.It induces a conditional dynamics, where a spin only precesses if all its neighbours are in the |↓⟩ state.Crucially this projected Hamiltonian is interacting, and in general non-integrable [36,37].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: where n = |↑⟩ ⟨↑|.This Hamiltonian may be interpreted as a toy model for an array of coherently driven Rydberg atoms.In these systems the state |↓⟩ represents an internal electronic ground state of an atom and the state |↑⟩ 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 Here Ω and ∆ denote the Rabi frequency and the laser detuning, and 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/Ω) 1/6 .In the so-called blockade approximation one considers the limit V ij → ∞ if two atoms are less than one blockade radius apart, and V ij → 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 ∆ = 0.If one sets V ij = V for next-nearest-neighbour sites, i.e. for |i − j| = 2a in 1D or |i − j| = √ 2a 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 [38] interesting to consider the case V < 0, i.e. with attractive next nearest neighbour interactions.In what follows, we will set Ω = 1.

C. 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 P, the ground state of (2) can host a plethora of symmetry broken ordered phases [33,34,39] as well as topologically ordered phases [40].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.

Equilibrium
The ground state phase diagram of the Hamiltonian in Eq. ( 2), at fixed finite values of V , depends only on ∆.For ∆ → −∞, the ground state is unique and given by |↓⟩ ⊗N independent of the lattice geometry.For ∆ → +∞ and any finite V , the ground state maximizes the number of spins in the |↑⟩ 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 |↑⟩ state and all spins on sublattice B are in the |↓⟩ state, or viceversa.We denote these states by ⊗B , respectively.In the limit ∆, V → +∞, with constant ∆/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 |Z 2 ⟩ , |Z ′ 2 ⟩ 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 [34].The latter can stabilize these non-Z 2 -ordered phases even at finite V .In fact, in 1D, the interplay between Z 2 and Z 3 order gives rise to several interesting phenomena such as floating phases and chiral transitions [41][42][43][44].In D > 1, the precise nature of the phase transition along all phase boundaries of ( 2) is less clear [34].Here, we focus on the regime where V , when positive, is sufficiently small, such that as ∆ is tuned from −∞ to +∞ the system undergoes a direct quantum phase transition from the gapped, disordered phase to the gapped, Z 2 ordered phase.When D = 1 and for V ≥ 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. [41].The location of the tricritical point is known analytically, since it lies on the integrable line ∆ = V − 1/V for V = −( √ 5 + 1)/2) 5/2 [33].A sketch of the (∆, 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 ∆, V , in 2D and for large enough ∆ the system undergoes a direct second order transition from the Z 2 phase to the striated phase, with V as driving parameter [34,45].Since our ansatz only captures Z 2 order, we will restrict our discussion to the regime where the effect of the 1/4-density ordered phase is negligible.

Quantum many-body scars
When ∆ = 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 [23,24] and have been proposed in Ref. [46] as the explanation of the anomalous dynamics observed in a Rydberg atom simulator initialized in one of the two |Z 2 ⟩ states [20].The anomaly consists in the fact that despite the Hamiltonian (1) being non-integrable, and the initial |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 |Z 2 ⟩ states.It is now believed that this phenomenon arises in the PXP model in any dimension whenever the lattice is bipartite [47], although numerical studies in D > 2 are absent in the literature.The name quantum scars comes from the analogy, put forward in Ref. [46], with single-particle physics, where the quantum system possesses special eigenstates that are "scars" of unstable periodic orbits of their, otherwise chaotic, classical counterpart [48].In this context, the one-dimensional version of the varia- 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 |Z2⟩ states.The dotted line represents the true time evolved state e −iH PXP t |Z2⟩.
tional ansatz which we will introduce in Sec.III was employed in Ref. [28] to map the constrained dynamics of the PXP model into classical equations of motion where the |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.

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 |ψ( ⃗ θ, ⃗ ϕ)⟩ can be represented as tensor network states with bond dimension χ = 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 |↑⟩ state.The variational manifold is parameterized by two angles, θ i ∈ [−π, π) and ϕ i ∈ [0, 2π), 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 (θ A , θ B , ϕ A , ϕ B ), although in the discussion here we will refer to the general case ( ⃗ θ, ⃗ ϕ) when it presents no additional difficulty.
The variational state |ψ( ⃗ θ, ⃗ ϕ)⟩ is a projected entangled pair state (PEPS) [49] with a PEPS tensor M (θ i , ϕ i ) defined at each lattice site i.In 1D, the state is a matrix where |↓⟩ = (1 0) T and |↑⟩ = (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 = M (θ, ϕ) 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.
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 = M (θ, ϕ) where a thick line collectively represents all incoming or outgoing legs.Introducing a bra-ket notation on the virtual indices, the tensor M reads Let I ∈ {0, 1} D be a binary string of length D, let I be the sum over all 2 D such strings, and denote by |I⟩ the corresponding product state in the virtual Hilbert space.We further use the short-hand notation ⃗ 0 = (0, 0, 0, . . . It is straightforward to see that this reproduces Eqs. ( 4) and ( 5) for D = 1, 2, respectively.The variational ansatz |ψ⟩ has several interesting properties.As already noted, it satisfies the blockade constraint.Remarkably, it is always normalized, i.e. ⟨ψ( ⃗ θ, ⃗ ϕ)|ψ( ⃗ θ, ⃗ ϕ)⟩ = 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 |Z 2 ⟩ = |↓⟩ ⊗A |↑⟩ ⊗B ; this can be checked by setting θ i = ±π, θ j = 0 for all sites i ∈ A, j ∈ B. This is in fact a special case of a more general property: product states with all spins on one sublattice in the state |↓⟩, 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 [28].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. (8) 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. ( 8) and provide a more in-depth discussion about its relation to the projected product state manifold in Appendix D.

(θA, θB) space for systems on bipartite lattices
For the applications that follow, we only consider two-site unit-cell translational invariant states |ψ(θ A , θ B , ϕ A , ϕ B )⟩.The ϕ variables are a priori relevant, but we will show that for the variational dynamics, the trajectories contain only states with ϕ A , ϕ B = 0, and that the variational ground states have ϕ A , ϕ B = ±π/2.Thus, the value of ϕ A and ϕ B can be considered fixed, and we have a two dimensional parameter space θ A , θ B ∈ [−π, π), 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 |Z 2 ⟩ states and the product states along the θ A and θ B axes, as well the non-perturbative (For the applications here ϕ is irrelevant and can be fixed to 0 or ±π/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.

A. Calculation of expectation values
The structure of the variational state |ψ⟩ = |ψ( ⃗ θ, ⃗ ϕ)⟩ allows us to efficiently calculate quantities of the form ⟨ψ| Ô|ψ⟩, for product operators Ô 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.
In general, for a bond-dimension χ tensor network state |Ψ⟩, evaluating quantities like ⟨Ψ| Ô|Ψ⟩ requires the contraction of a tensor network of bond dimension χ 2 .However, even though the tensor networks states introduced above have bond dimension χ = 2, it turns out that the evaluation of quantities of the form ⟨ψ| Ô|ψ⟩ 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 ⟨ψ| Ô|ψ⟩ 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.

Perturbative expansion
For D ≥ 2, we introduce a method to express the tensor network contraction as a perturbative expansion.This perturbative method, while powerful, is quite technical and does not have a simple intuitive interpretation.In this section, we give a high-level overview and then simply state the result, while a full derivation of the method can be found in Appendix A 3.
The main idea is to decompose each tensor in the network into a sum of two tensors.This seems counterproductive at first, since it expresses the quantity ⟨ψ| Ô |ψ⟩ 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 sin 2 (θ A /2) and sin 2 (θ B /2) as the expansion parameter.
Let us illustrate this.For simplicity, we focus on the case of a local operator Ô acting on a single site in a 2D lattice.An important object in the tensor network representing the quantity ⟨ψ| Ô |ψ⟩ is the tensor T formed by contracting the PEPS tensor M with its conjugate: As mentioned above, it is useful to express T as a sum of two tensors: The black dot is a constant tensor (i.e.independent of the variational parameters), while the blue dot carries a factor of sin 2 (θ/2).In the notation introduced in Appendix A 1, the black dot is p and the blue dot is − sin 2 (θ/2)q.These tensors have the crucial property that the following contraction is zero: Using the decomposition (10) and the property (11), it is easy to see ⟨ψ| Ô |ψ⟩ can be expressed as a sum of tensor network contractions of the following form: Here the red tensor is obtained by inserting the local operator Ô between the tensors M and M * in eq. ( 9).Eq. ( 12) shows a subset of diagrams that give non-zero contribution at various order of sin(θ/2).To obtain the perturbative expansion, we collect all terms that contribute in the same order of sin(θ/2).For the unit cell translation invariant states, |ψ(θ A , θ B , ϕ A , ϕ B )⟩, that we consider below, the quantity ⟨ψ| Ô |ψ⟩ can thus be expressed as where the prefactor h(θ A , θ B , ϕ A , ϕ B ) is a function that depends on the operator Ô, while f is a simple matrix of "counting factors", that also depends on Ô, but is independent of the variational parameters.We note that Eq. ( 13) 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 ∈ R such that f n,m < a n b m for all n, m ∈ N), causing the series to eventually diverge at high enough order.However, we find that in regions of the (θ A , θ B ) plane where the effective expansion parameters sin 2 (θ A /2) and sin 2 (θ 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 A ⟩GS − ⟨σ z B ⟩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.
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.
Finally, we note that the perturbative expansion can be interpreted as an expansion around product states.For states in the variational manifold that are product states (i.e.θ A = 0 or θ B = 0), the sum (13) contains only a finite number of terms.As seen in Fig. 5, the perturbative regime is located in the region around the set of product states, with the non-perturbative region corresponding to the states in the manifold with higher entanglement.We also note that states with increasing correlation length require higher order expansion terms.For example, if one considers the equivalent of Eq. ( 12) for a two-point function ⟨ψ| Ôi Ôj |ψ⟩, one would have to go to at least order |i − j| to obtain nontrivial results.

IV. GROUND STATE PHASE DIAGRAM
As a first application of the variational manifold, we calculate the energy E = ⟨ψ| H |ψ⟩ 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 |ψ(θ A , θ B , ϕ A , ϕ B )⟩.However, we can restrict the manifold to real states by setting ϕ A = ϕ B = π/2, as the Hamiltonian is real.We consider the energy function E(θ A , θ B ) on the variational space θ A , θ B ∈ [−π, π) and find the variational ground state by minimizing E(θ A , θ B ). Fig. 6 shows the resulting phase diagrams.
Let us first qualitatively describe the landscape of the optimal variational parameters θ GS A , θ GS B as the Hamiltonian is tuned.Suppose we fix V and allow ∆ to vary: when ∆ = −∞, the ground state is the | ⃗ 0⟩ state, and when ∆ remains large and negative, we find that θ GS A = θ GS B at the minimum energy state of the variational manifold.As ∆ reaches a critical value ∆ c , the variational ground state changes from the symmetric phase to a symmetry-broken phase where θ A ̸ = θ B .When ∆ > ∆ c , the variational ground state becomes twofold degenerate, with nonzero θ GS A − θ GS B .For V = 0, we find that the transition occurs at ∆ 1D c ≃ 0.77 in 1D, at ∆ 2D c ≃ −0.45 in 2D, and at ∆ 3D c ≃ −1.0 in 3D.In 1D it is known from exact diagonalization studies that ∆ c ≃ 1.31 [50].In 2D, exact diagonalization on small systems sizes combined with finite-size scaling indicate that ∆ c ≃ −0.2 (see Appendix G).In 3D, the critical point for 64 sites is at ∆ c ≃ −0.61.Thus we see that in all cases the variational method approximates, but underestimates the transition point.
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 (θ GS A , θ GS B ) evolves smoothly with ∆, while for V < V c it changes abruptly at the transition point.We find that V c ≃ −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 ≃ −3.33 [33], although it is still notable that the effect is captured by the variational method.For ∆, we find that ∆ tricritical ≃ −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 (∆, V ) phase diagrams, using the order parameter ⟨σ z A ⟩ GS − ⟨σ z B ⟩ GS .In Fig. 7 we present the one-dimensional plots of minimum energies along lines of constant , for various values of ∆ 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 β = 1/2 (see Fig. 8(c)).

A. 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 npoint function.Here, we focus on the density-density connected two-point functions f (|i − j|) = ⟨n i n j ⟩ − ⟨n i ⟩⟨n j ⟩, which we compute on the optimal variational ground state for various values of ∆ 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|/ξ) we can extract the correlation length ξ, which we plot as a function of ∆ (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 [51], 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) [31,32] to study the time evolution of the PXP model Eq. ( 1) with Ω = 1, and show that it captures the periodic revivals characteristic of quantum scarred systems.The 1D results were already obtained in Ref. [28] and are presented here for comparison with the D > 1 case.

A. TDVP equations of motion
In general, the variational manifold is parameterized by the variables ⃗ θ and ⃗ ϕ.However, since we are mainly concerned with time evolution from the initial states ⊗B , it suffices to consider unit cell translationally invariant states.Also, since the TDVP is energy conserving [31,32], it can be shown that d ⃗ ϕ/dt = 0 when starting from the |Z 2 ⟩ states, and we can set all ⃗ ϕ = 0 (see Ref. [28] and Appendix C 3 b).Thus our variational manifold is two dimensional, parameterized by the position vector ⃗ θ = (θ A , θ B ).We calculate the dynamics within this variational space by minimizing the leakage rate, or "quantum leakage", given by Γ 2 = ⟨δ|δ⟩, where |δ⟩ 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 ⃗ θ = ( θA , θB ) that, due to symmetry, may be expressed in the form θA = f (θ A , θ B ) and θB = f (θ B , θ A ), where f (θ A , θ B ) can be calculated by exact tensor network contraction in D ≤ 2 and via perturbative expansion in any D. The heatmap is the leakage rate γ = ⟨δ|δ⟩ /N ; note that the leakage rate is exactly zero along the x and y axes.
Along the axes of the (θ A , θ B ) space, the equations of motion take on the following, particularly simple form (derived in Appendix F) where θA and θB are the unit vectors.In 1D there also exists a closed-form expression for the equations of motion in general; it was shown [28] that these tensor network calculations can be performed analytically to show that f

B. Periodic trajectories
In the TDVP framework, the interacting spins problem is approximated by a two dimensional nonlinear dynamical system with equations of motion d ⃗ θ/dt = ⃗ θ(θ A , θ B ).These equations can be integrated numerically to obtain trajectories starting from any initial condition, and the vector field ⃗ θ(θ A , θ B ) can be straightforwardly analyzed to infer the stability of these trajectories.In Figs.9-10 we plot ⃗ θ(θ A , θ B ) and the leakage rate in the (θ A , θ B ) plane of the exact equations of motion obtained analytically in Ref. [28] 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 |Z 2 ⟩ and |Z ′ 2 ⟩ states; the all spin down state |↓⟩ ⊗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 In each case (1D, 2D, and 3D), we observe a periodic path connecting the |Z 2 ⟩ and |Z ′ 2 ⟩ states; these nonergodic trajectories correspond to the revival behaviour observed in the quantum spin system and it is easy to see from the vector field plot of ⃗ θ(θ A , θ B ) that they are unstable orbits alike the one-dimensional case.The path periods, i.e. the time it takes to go |Z 2 ⟩ → |Z ′ 2 ⟩ → |Z 2 ⟩, are T TDVP = 4.820, 5.168, 5.345 for 1D, 2D, 3D, respectively.These have to be compared to the revival times T 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 TDVP − T Exact = 0.035, 0.014, 0.005 for D = 1, 2, 3).
We also plot the leakage rate γ = ⟨δ|δ⟩ /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 γdt ≃ 0.17 in 1D, γdt ≃ 0.16 in 2D, and γdt ≃ 0.13 in 3D.The leakage rate is exactly zero along the lines θ A = 0 and θ 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 lowentanglement states.A simple example is the path starting at the origin (i.e. with initial state |↓⟩ ⊗N ) and ending at (π, π), which has leakage γdt ≃ 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π 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 |Z 2 ⟩ = |↑⟩ ⊗A |↓⟩ ⊗B .Ini-tially 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 π.This is corroborated by the variational equations of motion Eq. ( 14): in the limit D → ∞, the velocity is 2 θA along the x axis and 2 θ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 [52].

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 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π.
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 [53].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 [21].

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 ⟨ψ| Ô|ψ⟩, where |ψ⟩ = |ψ( ⃗ θ, ⃗ ϕ)⟩ is the variational state and Ô is a local operator.We will show that it is possible to calculate ⟨ψ| Ô|ψ⟩ 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 ≥ 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. ⟨ψ|ψ⟩ = 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: where parallel in and out legs have been combined, so the legs of T have dimension 4. The quantity ⟨ψ|ψ⟩ 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 ⟨ψ|ψ⟩ or ⟨ψ| Ô|ψ⟩, we can perform a "bond dimension reduction" operation where t abcd = T (aa)(bb)(cc)(dd) for a, b, c, d ∈ {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 Using the matrix elements of M from Eq. ( 8), we can find the matrix elements of T = ⟨M |M ⟩, which is given by Like in the 2D case, we can perform a bond dimension reduction, and then ⟨ψ|ψ⟩ is equal to the contraction of a network of bond-dimension-2 tensors t(θ), which we express diagrammatically as and are given by We can write t as t(θ) = p − sin 2 (θ/2)q, where the constant tensors p and q are p = |α⟩ ⟨ ⃗ 0| (A7a) The decomposition t = p − sin 2 (θ/2)q turns out to be very useful in performing calculations.Observe that in the network of ts, 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. (A7)) 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, = 0.
(A8) 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.Note that Eq. ( 11) is the 2D version of Eq. (A8) here.

a. Normalization
An immediate corollary is that the variational state is normalized in the thermodynamic limit, i.e. that ⟨ψ|ψ⟩ = 1.There is one additional condition that we must impose: we must assume that at most a finite number of the θ i s in ⃗ θ are integer multiples of π.That is, we assume that sin(θ i /2) ̸ = ±1 for all but a finite number of θ i s.
In the expansion of ⟨ψ|ψ⟩ using t = p − sin 2 (θ/2)q, there are no terms with a finite number of qs due to Eq. (A8): there are either no qs, 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 qs give zero contribution, because each q comes with a factor of sin 2 (θ/2).Thus ⟨ψ|ψ⟩ = (network of all ps) = 1.

b. Calculation by explicit contraction
In D ≤ 2, the most straightforward way to calculate ⟨ψ| Ô|ψ⟩ 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 × 2 L transfer matrix (rather than 4 L × 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.

Perturbative expansion
Using the results of the previous section, we show how one can express ⟨ψ| Ô |ψ⟩ as an infinite series.The exposition here is quite formal: for a practical introduction, the reader may skip to Sec.A 3 a below, which is a natural continuation of the high-level overview in Sec.III A 1.
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(• • •) means a contraction over the infinite network, and when tensors inside a Tr(• • •) are "multiplied", they are being placed on lattice sites and then contracted.In the one dimensional case, the tensors are matrices, Tr(• • •) is indeed a trace, and the equations here may be taken literally.
Recall that ⟨ψ|ψ⟩ is the contraction of a network of t tensors.Let Θ be the tensor that results when Ô is sandwiched with the PEPS tensors: then ⟨ψ| Ô |ψ⟩ is the contraction of a network with Θ on a local region and t on all other sites.That is, where I is the region (i.e.set of sites) that Θ covers.Note that due to the mechanics of the bond dimension reduction operation, I may be slightly larger than the region that Ô itself acts on.Using t(θ i ) = p − sin 2 (θ i /2)q, we can expand: We have expressed ⟨ψ| Ô |ψ⟩ 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 Θ 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(Θ q i p j ) will be zero, due to Eq. (A8), unless the region J is "anchored" by Θ.When Tr(Θ q i p j ) is nonzero, it can usually be computed quickly by reading off the tensor elements of Θ.
The perturbative expansion is thus a sum over allowed regions J.This is the general idea; the numerical implementation is described below.

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 |ψ(θ A , θ B , ϕ A , ϕ B )⟩.Using Eq. (A13), we can express ⟨ψ| Ô |ψ⟩ as where the prefactor h(θ A , θ B , ϕ A , ϕ B ) is some function that comes from the aforementioned Tr(Θ q i p j ), and is usually not hard to evaluate, and where Σ(f, θ A , θ B ) is the perturbative expansion: where f is a matrix of "counting factors" that depends on Ô and is obtained by counting the ways that the q i s may be arranged in Tr(Θ q i p j ).To see how exactly Eq.(A15) arises, it is instructive to work out the calculation for a specific example.

a. Example: ⟨ψ| ni |ψ⟩
Consider the case where Ô is the number operator n i , for some site i ∈ A. Note that since the state |ψ⟩ lies inside the constrained subspace, this is equivalent to calculating ⟨ψ| Pn i P |ψ⟩.In this case, the tensor Θ occupies only a single site, and the prefactor h is simply 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. (A15), 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 Θ on one site, and t on all other sites.We use the fact that t = p − sin 2 (θ/2)q, and then expand in powers of sin 2 (θ/2).That is, we want to express the result of the contraction of the tensor network as sum = (Θ with p on all other sites) + (configurations with one q) + (configurations with two qs) + . . . .The allowed configurations are constrained by where the qs can be placed: each q must be above or to the left of another q, or Θ.Also, note that each q comes with a factor of − sin 2 (θ/2).We can write the expansion as in Eq. (A15), where, in each term in the sum, n is the number of qs on sublattice B and m is the number of qs 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 Θ 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 Θ) 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 ): The perturbative expansion only is valid in the region where θ A or θ B or both are small enough that f n,m (sin 2 (θ A /2)) n (sin 2 (θ B /2)) m decreases with larger values of n and m.

b. Two point functions
The method described here also allows for the calculation of two point functions, which in practice corresponds to calculating ⟨ψ| Ô |ψ⟩ where Ô is a non-local operator consisting of two parts.The simplest example would be ⟨ψ| n i n j |ψ⟩, which may be calculated by a generalization of the procedure described above for ⟨ψ| n i |ψ⟩.

Asymptotic nature of the expansion
Here we continue the discussion on the perturbative and non-perturbative regimes from Sec. III A. 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 Ô = n i ; that is, we use f n,m = f (2D) n,m from Eq. (A16).However, the results are valid for any local operator Ô 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. (A15).We define the order of a truncated expansion to be N ≡ max(n + m).Let s n,m refer to a term in the sum, and let S N ≡ n+m≤N s n,m be the partial sum up to order N .To probe the accuracy of the expansion at a certain order, we consider the quantity |S N − S N −1 |.In Fig. 11

a. Proof of non-convergence
Here we prove that the series expansion Eq. (A15) divertges.In Section III A it was claimed that the counting factors f n,m scale super exponentially, i.e. that there do not exist any a, b ∈ R such that f n,m < a n b m for all n, m ∈ N.This would imply that the sum (A15) diverges everywhere.Here, we prove that for any N ∈ N there exists n, m > N such that f n,m > (n + m − 1)!/(n + m), which is a strong enough condition to imply superexponential 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 ≡ 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 "perime-FIG.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.
ter" 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  From these calculations we obtain the energy E(θ A , θ B ), plotted in Fig. 13 and Fig. 14, and find the variational ground state by finding the θ A , θ B that minimizes E(θ A , θ B ).

PXP term
Let us first consider E PXP = ⟨ψ| H PXP |ψ⟩.Using P |ψ⟩ = |ψ⟩, we have 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 ⟨ψ| σ x i |ψ⟩, we must consider the following tensor: 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 ⃗ 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 ⟨ψ| σ x a |ψ⟩ is the contraction of the tensor network with S(θ A , ϕ A ) on a single site, and with T (θ A ) or T (θ B ) on all other sites: the calculation may be performed via either explicit contraction (D ≤ 2) or the perturbative method (D ≥ 2).Upon performing this calculation, we find that the θ, ϕ dependence decouples in a way that we may write ⟨ψ| σ x a |ψ⟩ = sin ϕ A F(θ A , θ B ), and thus for the energy we have It is easy to check that due to the symmetry of F(θ A , θ B ), one can without loss of generality set ϕ A = ϕ B = π/2.

Detuning and next-nearest-neighbour terms
For the detuning and NNN terms, we have, similar to the above: FIG. 13.Plots of the energy E(θA, θB).Top row: 1D, at ∆ = 0 (left) and ∆ = 1 (right).Bottom row: 2D, at ∆ = −1 (left) and ∆ = 0 (right).Thus the plots on the left are for ∆ < ∆c and on the right for ∆ > ∆c.In each plot, the red dot(s) indicate the minimum of E(θA, θB).The line θA = θB is included as a visual aid.The red dot is at the origin at ∆ = −∞.As ∆ is increased, it travels down the line, until it splits into two at ∆ = ∆c when the variational ground state becomes degenerate.
We thus need to consider the tensor: It follows that the quantity ⟨ψ| n i |ψ⟩ is the contraction of a tensor network with J(θ i ) on a single site and T on all other sites, and that ⟨ψ| n i n j |ψ⟩ would be the contraction of a tensor network with J(θ i ) on site i and J(θ j ) on site j and T on all other sites.

Appendix C: TDVP calculations
Here we present the technical details of the variational dynamics calculations.We first give a brief but selfcontained review of the TDVP.

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 ⃗ ϕ = 0, so that ultimately our variational manifold is parameterized by only two variables, θ A and θ B .
Let the variational state be |ψ(θ A , θ B )⟩. Define the vector |δ⟩ 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 d|ψ⟩ dt = H |ψ⟩. In the following we use Greek letters as placeholders for the sub-lattices A and B, and use the notation |∂ µ ψ⟩ ≡ ∂|ψ⟩ ∂θµ .The leakage rate vector is given by The leakage rate is defined as Γ 2 = ⟨δ|δ⟩.To derive the TDVP equations of motion, we minimize the leakage rate with respect to θA and θB .That is, we have ∂Γ 2 /∂ θA = ∂Γ 2 /∂ θB = 0, which gives us the equations of motion In deriving these equations, we used the fact that ⟨∂ A ψ|∂ B ψ⟩ = 0, i.e. the Gram matrix G µν ≡ ⟨∂ µ ψ|∂ ν ψ⟩ is diagonal.This is nontrivial, but we will derive this below.
After performing the calculations for the numerators and denominators of Eq. (C2), it turns out that we can write the equations of motion as θA = f (θ A , θ B ) and θB = f (θ B , θ A ), where where D is the dimension of space, and K, S, and G are each the contraction of a single (infinite) tensor network, the technical details of which are described below.

Gram matrix calculation
Here we describe the calculation of the Gram matrix, i.e. quantities of the form ⟨∂ µ ψ|∂ ν ψ⟩ ≡ G µν .We have, from the product rule of calculus, that |∂ A ψ⟩ = a∈A |∂ a ψ⟩, so that ⟨∂ µ ψ|∂ ν ψ⟩ = i∈µ,j∈ν ⟨∂ i ψ|∂ j ψ⟩, where |∂ i ψ⟩ means that the derivative is only applied on the site i.In other words, |∂ i ψ⟩ is a tensor network state with the tensor ∂M (θ i ) on site i, and the usual M (θ j ) tensor on all other sites j ̸ = i.The tensor ∂M is defined as Define the tensor It turns out that, after performing the bond-dimension reduction, that ∂T ∝ q.This allows us to conclude that ⟨∂ i ψ|∂ j ψ⟩ = 0 for all pairs of sites i ̸ = j.This leads to the conclusion that the Gram matrix is diagonal, and that G AA = (N/2)⟨∂ a ψ|∂ a ψ⟩, 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 performing a bond dimension reduction as before.The elements of g are given by g(θ) = (1/4) sin Thus ⟨∂ a ψ|∂ a ψ⟩ is the contraction of a tensor network with g(θ A ) on a single site and T on all other sites.We define the aforementioned G as G(θ A , θ B ) = ⟨∂ a ψ|∂ a ψ⟩ = G AA /(N/2).By symmetry, the contraction of a single g(θ B ) tensor on a B lattice site is G(θ B , θ A ) = ⟨∂ b ψ|∂ b ψ⟩ = 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.

Calculation of ⟨∂µψ|H|ψ⟩
Here we discuss the calculation of ⟨∂ µ ψ|H|ψ⟩, the numerator of the equations of motion, Eq. (C2).Since |ψ⟩ is inside the Rydberg-blockaded subspace we have We will see that ⟨∂ j ψ|σ x i |ψ⟩ 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).
a. Calculation of K Let us first consider the i = j contribution, which we call K. Define the tensor K as where we find that the elements are given by K(θ) = (i/2) cos 2 (θ/2)| ⃗ 0⟩⟨ ⃗ 10| + (i/2) sin 2 (θ/2)| ⃗ 0⟩⟨ ⃗ 01|.The quantity K(θ A , θ B ) = ⟨∂ a ψ|σ x a |ψ⟩ is the contraction of a network with K(θ A ) on one site a ∈ A and T on all other sites.A ↔ B symmetry gives the analogous quantity K(θ B , θ A ) = ⟨∂ b ψ|σ x b |ψ⟩.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 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 K(θ A , θ B ).

b. Calculation of S
Next we consider ⟨∂ j ψ|σ x i |ψ⟩ when i ̸ = j.To start, consider the tensor S as defined in Eq. (B2), but we now have ϕ A = ϕ B = 0: In this case, we find that = 0. (C11) Note that this implies that ⟨ψ| σ x i |ψ⟩ = 0 which implies that E = ⟨H⟩ = 0 when |ψ⟩ is a variational state with ⃗ ϕ = 0. (Eq.(C11) 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 ̸ = j contributions are the cases where i is next to j.By symmetry, we only need to consider the diagram multiplied by a factor of D. We contract this object with the T tensor on all other sites, to get ⟨∂ j ψ|σ x i |ψ⟩.We define S(θ A , θ B ) = ⟨∂ a ψ|σ x b |ψ⟩ and S(θ B , θ A ) = ⟨∂ b ψ|σ x a |ψ⟩.

Appendix D: Projected product states
Consider the set of projected (i.e.Rydberg blockaded) product states: where the state of each site is a spin coherent state |ϑ i , φ i ⟩ = cos(ϑ i /2) |↓⟩ − ie iφi sin(ϑ i /2) |↑⟩.These states resemble the variational ansatz, but are not normalized.
As noted earlier, P is the application of the twosite operator P ij = |↓↓⟩ ⟨↓↓| + |↓↑⟩ ⟨↓↑| + |↑↓⟩ ⟨↑↓| = 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 = , where the internal line has bond dimension 2. Specifically, P R and P L are given by: P 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 (A4), the PEPS tensor at each site is given by Comparing this expression to Eq. ( 8) one can see that the only difference lies in the coefficient of the component |β⟩ ⟨ ⃗ 0|.In 1D, a gauge transformation exists that maps the PEPS state resulting from the tensors contraction into each other [28].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. (D4) and Eq. ( 8).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 θ A = θ B [54], whereas this phase transition is absent in our ansatz (8), 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 ⟨ψ| Ôlocal |ψ⟩, 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 |δ⟩: 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 ⟨∂ µ ψ|H|ψ⟩ is pure imaginary, so that ⟨∂ µ ψ|H|ψ⟩ = −⟨ψ|H|∂ µ ψ⟩.
We can plug our derived expressions for θA and θB into the leakage rate expression to get The new work that we have to do is to calculate ⟨ψ|H 2 |ψ⟩.First, note that there are two equivalent ways of writing down the PXP Hamiltonian, which we call the "global" and "local" versions: Here, P = |↓⟩ ⟨↓|, and P is the whole-system Rydbergblockade projection, i.e.P = ⟨ij⟩ P ij , where where in the last step we used the fact that we know from the explicit tensor network calculation that ⟨ψ|σ x i σ x j |ψ⟩ is nonzero only for adjacent sites.Specifically, by nonadjacent, we mean neither directly adjacent nor diagonally adjacent.The quantity ⟨ψ|σ x a σ x a ′ |ψ⟩ is nonzero if a and a ′ 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 P operators can be commuted out, except for the one acting on i and j, so we have ⟨ψ| σ x i Pσ x j |ψ⟩ = ⟨ψ| σ x i P ij σ x j |ψ⟩ Using symmetry, we can write the sum as 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 ⟨∂ A ψ|∂ A ψ⟩ = N 2 ⟨∂ a ψ|∂ a ψ⟩ (and the analogous equation for B) for the quantity µ θ2 µ ⟨∂ µ ψ|∂ µ ψ⟩ , we now have the full leakage rate expression: where a and a ′ have a specific relationship as mentioned above, and where the (a ↔ 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 θ A and θ B axes in Fig. 9 and Fig. 10.We derive the equations of motion Eq. ( 14) and demonstrate that the leakage rate is zero, i.e. that γ(θ A , 0) = γ(0, θ B ) = 0.
Let us consider the case θ B = 0; by symmetry, the case θ A = 0 will be the same but with θ A ↔ θ B .When doing the TDVP calculations, we work with the variational state |ψ(θ A , θ B )⟩, and the states |∂ a ψ(θ A , θ B )⟩ and |∂ b ψ(θ A , θ B )⟩. Recall that the latter states are defined as applying a derivative ∂/∂θ A or ∂/∂θ B to |ψ(θ A , θ B )⟩ at any one particular site a ∈ A or b ∈ B. In general, these are tensor network states; however, when θ B = 0, these become product states, allowing one to perform the TDVP and leakage rate calculations by hand.
The state |ψ(θ A , 0)⟩ is the following product state: 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 |Z 2 ⟩ state.To extract the phase boundaries depicted in Fig. 6 we computed the ground state fidelity susceptibility which exhibits a peak at the finite-size phase transition point.We show in Fig. 15 the finite size scaling of F for V = 0 in 2D, from which we extrapolate the thermodynamic limit critical value ∆ c ≃ −0.19, and for V = 0 in 3D for the 4 × 4 × 4 cube of 64 sites.
For the computation of the revival periods reported in Sec.V, we need to compute the fidelity |⟨Z 2 |e −iHt |Z 2 ⟩|.
Although the |Z 2 ⟩ state breaks the lattice symmetry of the Hamiltonian, the revival time can be extracted from the time evolution of the cat state |ψ 0 ⟩ = (|Z 2 ⟩ + |Z ′ 2 ⟩)/ √ 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 ≥ 16 (4×4 square).For the 3D revival time a cube with 4 atoms in each direction (N = 64) is required, otherwise the system is

FIG. 2 .
FIG.2.Sketch of the phase diagram for the generalized PXP model with next-nearest-neighbour interaction Eq. (2).A disordered phase and a Z2 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.

FIG. 5 .
FIG. 5. Depiction of the (θA, θB) space that our variational state |ψ(θA, θB)⟩ is in.Points along the axes are all product states, where we have defined |θϕ⟩ = cos(θ/2) |↓⟩ − ie iϕ sin(θ/2) |↑⟩, and |Z2⟩ = |↑⟩ ⊗A |↓⟩ ⊗B and |Z ′ 2 ⟩ = |↓⟩ ⊗A |↑⟩ ⊗B .(For the applications here ϕ is irrelevant and can be fixed to 0 or ±π/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.

FIG. 6 .
FIG. 6. (∆, V ) phase diagrams for the generalized PXP model (2) with the magnetization ⟨σ zA ⟩GS − ⟨σ z B ⟩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.

FIG. 7 .
FIG. 7. Plots of the minimum energy (per site) along lines of constant θA − θB for D = 1, 2, 3, showing the second order transition (top row) and first order transition (bottom row).

FIG. 9 .
FIG. 9. Flow diagrams and paths for 1D (left) and 2D (right).The heatmap is the leakage rate γ = ⟨δ|δ⟩ /N ; note that the leakage rate is exactly zero along the x and y axes.

FIG. 10 .
FIG.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 (N = 12 (2D) and (N = 9 (3D) at which the calculation is performed.The heatmap is the leakage rate γ = ⟨δ|δ⟩/N .

ACKNOWLEDGMENTS
This research was funded in whole or in part by the Austrian Science Fund (FWF) [grant DOI 10.55776/COE1].For open access purposes, the author has applied a CC BY public copyright license to any author accepted manuscript version arising from this submission.This work is also supported by the European Union's Horizon Europe research and innovation program under Grant Agreement No. 101113690 (PASQuanS2.1), the ERC Starting grant QARA (Grant No. 101041435), 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).
) where the |α⟩ , |β⟩ notation is the same as originally introduced for Eq.(8).The exact meaning of Eq. (A4) 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|⃗ a⟩⟨ ⃗ b| ⊗ |⃗ c⟩⟨ ⃗ d| shows up in the right hand side of Eq. (A4), it means that the element of T with incoming indices set to (a 1 c 1 ), (a 2 c 2 ), . . ., (a D c D ) and outgoing indices set to (b 1 d 1 ), (b 2 d 2 ), . . ., (b D d D ) has the value x.

FIG. 11 .
FIG. 11.Plot of |SN −SN −1| in log scale, across the parameter space, for N = 12 in 2D (left) and N = 9 in 3D (right).Note that the colormap scale has been truncated in the sense that numbers < 10 −5 have been neglected.
we show |S N − S N −1 | for the highest orders of the expansion as used in the calculations in this paper, N = 12 in 2D and 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 N − S N −1 | < 10 −3 at order N = 12 (2D) and 9 (3D).The boundaries of the perturbative regimes in Fig. 5 in the main text are the points |S N − S N −1 | = 10 −3 from the same data as in Fig. 11.
and by induction, f k > √ k!.Now let us recast the result in terms of n and m.Here let fn,m be the largest f n,m such that n + m = k.Clearly k fn,m > f k .Then f k > √ k! implies that (n + m) fn,m > (n + m)!, which is the claimed result.
Appendix B: Energy calculationHere we elaborate on the specifics of the tensor network calculations for the variational energy ⟨ψ| H |ψ⟩. Let us consider the parts of the Hamiltonian (2), H = H PXP + H ∆ + H NNN one by one.

FIG. 14 .
FIG.14.plots of the energy E(θA, θB) for 3D. at ∆ = −1 (left) and ∆ = −0.6 (right).Note that we zoomed in on a smaller window near the origin.The red dot is the minimum.
horizontal leg.We can use this to write the state |ψ( ⃗ ϑ, ⃗ φ)⟩ as a tensor network.In one dimension, the MPS tensor is = .(D3)

19 FIG. 15 .
FIG.15.Ground state fidelity susceptibility for V = 0 as function of ∆ 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 Z2-ordered phase.The inset shows the extrapolation of the thermodynamic limit critical value of ∆ resulting in ∆c ≃ −0.19.

= 1 √ 2 (|Z 2 + |Z 2 FIG. 16 .
FIG.16.Revival fidelity employed to extract the first period of the approximately-periodic many-body dynamics generated by the |Z2⟩ initial state.The upper panel shows a comparison between the symmetric and non-symmetric dynamics starting from the cat state (|Z2⟩ + |Z ′ 2 ⟩)/ √ 2 and from the symmetry-breaking |Z2⟩ state, in 1D for N = 18.The former case produces a revival at a time t ≃ T /2 due to the fact that the periodic orbit from the |Z2⟩ state back to itself goes through the |Z ′ 2 ⟩ 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.
1, 1, 0).In 1D, the local version can also be written as H = i P i−1 σ x i P i+1 .Working with the global version, let us expand ⟨ψ|H 2 |ψ⟩: ⟨ψ|H 2 |ψ⟩ = At this point, it is convenient to use the local H for the i = j terms, and the global H for the i ̸ = j terms.For the i = j terms, we have If i and j are not next to each other, all of the P operators in P can be commuted out (this is easiest to see by drawing 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 |ψ⟩).So, for non-nearest-neighbour sites,