Entanglement in the quantum phases of an unfrustrated Rydberg atom array

Recent experimental advances have stimulated interest in the use of large, two-dimensional arrays of Rydberg atoms as a platform for quantum information processing and to study exotic many-body quantum states. However, the native long-range interactions between the atoms complicate experimental analysis and precise theoretical understanding of these systems. Here we use new tensor network algorithms capable of including all long-range interactions to study the ground state phase diagram of Rydberg atoms in a geometrically unfrustrated square lattice array. We find a greatly altered phase diagram from earlier numerical and experimental studies, revealed by studying the phases on the bulk lattice and their analogs in experiment-sized finite arrays. We further describe a previously unknown region with a nematic phase stabilized by short-range entanglement and an order from disorder mechanism. Broadly our results yield a conceptual guide for future experiments, while our techniques provide a blueprint for converging numerical studies in other lattices.


I. INTRODUCTION
Rydberg atom arrays, where cold atoms are trapped in an optical lattice and interact via excitation into Rydberg states [1,2], have generated interest for quantum information processing and to realize exotic many-body states [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].A recent experiment [18], backed by numerical studies [19,20], has suggested a richness in 2D Rydberg atom array ground states on a square lattice.However, although the observed, non-disordered, phases are not all classical crystals, they contain little entanglement [19].Thus it remains unclear whether such arrays realize non-trivial entangled quantum ground-states on simple lattices.A confounding feature suggested more recently [17], is that the long-range tails of the interactions greatly affect the phases, complicating the accurate numerical determination of bulk behaviour.Here, we describe new numerical techniques that greatly reduce finite size effects, allowing us to confidently converge the bulk phase diagram.We also showcase techniques that address large finite two-dimensional lattices realized in experiments, while incorporating all long-range interactions.Unexpectedly, we derive quite different physics from our simulations compared to both previous theoretical and experimental analyses -including the emergence of a non-trivial, entangled nematic phase, even on the unfrustrated square lattice array.
The Rydberg atom array Hamiltonian is Here σx i = |0 i 1 i |+|1 i 0 i | and ni = |1 i 1 i | ({|0 i , |1 i } denote ground and Rydberg states of atom i). a is lattice spacing, Ω labels Rabi frequency, and δ describes laser detuning.V parameterizes the interaction strength between excitations.This can be re-expressed in terms of the Rydberg blockade radius R b , with V /(R b /a) 6 ≡ Ω.We study the square lattice in units a = Ω = 1 [19], yielding two free parameters δ and R b .
The ground states of this Hamiltonian are simply understood in two limits.For δ/Ω 1, R b = 0, the system is classical and one obtains classical crystals of Rydberg excitations [21][22][23][24] whose spatial density is set by the competition between δ and R b .For δ/Ω 1, R b = 0, Rydberg excitations are disfavored and the solutions are dominated by Rabi oscillations, leading to a trivial "disordered" phase [19,25,26].In between these limits, it is known in 1D that no other density-ordered ground states exist besides the classical-looking crystals, with a Luttinger liquid appearing on the boundary between ordered and disordered phases [26].
In 2D, however, the picture is quite different.An initial study [19] using the density matrix renormalization group (DMRG) [27][28][29][30] found additional quantum crystalline (or "density-ordered") phases, where the local excitation density is not close to 0 or 1.A recent experiment on a 256 programmable atom array has realized such phases [18].However, as also discussed there, the density-ordered phases are unentangled quantum meanfield phases, and thus not very interesting.In addition, more recent numerical results [17] highlight the sensitivity of the physics to the tails of the Rydberg interaction and finite size effects.Thus, whether Rydberg atom arrays on a simple unfrustrated lattice -such as the square lattice -support interesting quantum ground-states, remains an open question.
Here, we resolve these questions through high-fidelity numerical simulations.To do so, we employ variational tensor network methods.Tensor networks have led to breakthroughs in the understanding of 2D quantum many-body problems [32], and we rely on two new techniques that address specific complexities of simulating interactions in Rydberg atom arrays.The first we term Γ-point DMRG, which captures interactions out to infinite range, while employing a traditional finite system two-dimensional DMRG methodology [27], removing interaction truncations and boundary effects present in earlier studies [16,17,19,20].This allows us to controllably converge the bulk phase diagram.The second is a representation of long-range interactions [31] compatible with projected entangled pair states (PEPS) [33][34][35][36].With this, we use PEPS to find the ground states of a Hamiltonian with long-range interactions for the first time, and specifically here, model the states of finite arrays of large widths as used in experiment.Both techniques can be used for more faithful simulations of Rydberg atoms in other settings.We first describe the new numerical meth-FIG.1. Numerical methods and strategy.(a) A schematic representation of Γ-point DMRG.A single infinite bulk configuration is given by periodic images of the central supercell configuration.The wavefunction coefficient for this infinite configuration is given by the contraction of a snake MPS, which is defined only within a single supercell.(b) By widely varying the size of the supercell, Γ-point DMRG obtains many different ground states.Identifying all accessible supercells which give the same ground state order (shown with identically colored points), we can ensure that all competing low-energy states are well converged w.r.t.finite size effects, and thus properly identify the true ground state (inset shows ground-state order (dark green) converged w.r.t.supercell size, separated from other low-energy orders by 10 −4 energy units).(c) A PEPS wavefunction ansatz with bond dimension D for a finite system.Each tensor is a different color because they can all be unique.(d) A simplified diagrammatic representation of the long-range Hamiltonian construction for PEPS in Ref. [31].All terms in the Hamiltonian are accounted for by a sum of Lx comb tensor network operators.Tensors of the same color are identical.
ods, before turning to the bulk and finite-size phase behavior of square lattice Rydberg arrays and the question of entangled quantum phases.

II. NUMERICAL STRATEGY AND TECHNIQUES
A. Bulk simulations and Γ-point DMRG A challenge in simulating Rydberg atom arrays is the long-range tails of the interaction.Because itinerancy only arises indirectly as an effective energy scale [25], the main finite size effects arise from interactions.Many previous studies have employed a cylindrical DMRG geometry common in 2D DMRG studies [27].However, there the interaction is necessarily truncated to the cylinder half-width, while along the open direction, edge atoms experience different interactions than in the bulk; both choices produce strong finite size effects.
To avoid these problems, we perform 2D DMRG calculations in a Bloch basis.The resulting Γ-point 2D DMRG formally models an infinite lattice (Fig. 1a) with a wavefunction constrained by the supercell, and periodic boundary conditions in both directions.This differs from using a periodic matrix product state (MPS) as periodicity is enforced by the Bloch basis rather than the MPS, and the underlying 2D DMRG can be carried out using the conventional snake path.It is also different from a cylindrical/toroidal geometry, which are both finite; here the lattice remains infinite.The Bloch basis interactions enter as an infinite lattice sum over lattice vectors (2) where L x , L y are the supercell side lengths.
The only finite size parameter is the supercell size L x × L y .We thus perform exhaustive scans over L x , L y .However, because no interactions are truncated and there are no edge effects even in the smallest cells, finite size effects converge very rapidly (much more quickly than using a cylinder or torus).Using different supercell sizes with up to 108 sites we converge the energy per site to better than 10 −5 , compared to the smallest energy density difference we observe between competing phases of ∼ 10 −4 (see Fig. 1b and supplementary information [37]).

B. Finite simulations and PEPS with long-range interactions
To simulate ground-states of finite arrays, we consider finite systems (with open boundaries) of sizes 9 × 9 up to 16 × 16 atoms.This resembles capabilities of near-term experiments [10,18].The width of the largest arrays challenges what can be confidently described with MPS and DMRG for more entangled states.Consequently, we employ PEPS wavefunctions which capture area law entanglement in 2D, and can thus be scaled to very wide arrays (Fig. 1c).Together with DMRG calculations on moderate width finite lattices, the two methods provide complementary approaches to competing phases and consistency between the two provides strong confirmation.However, PEPS are usually combined with short-range Hamiltonians.We now discuss a way to combine longrange Hamiltonians efficiently with PEPS without truncations.
For this, we rely on the representation we introduced in Ref. [31].This encodes the long-range Hamiltonian as a sum of "comb" tensor network operators (Fig. 1d).As discussed in Ref. [31], arbitrary isotropic interactions can be efficiently represented in this form, which mimics the desired potential via a sum of Gaussians, i.e.
where k max ∼ 7 for the desired accuracy in this work).The combs can be efficiently contracted much more cheaply than using a general tensor network operator.
While Ref. [31] described the Hamiltonian encoding, here we must also find the ground-state.We variationally minimize Ψ| Ĥ|Ψ using automatic differentiation [38].Combined with the comb-based energy evaluation, this allows for both the PEPS energy and gradient to be evaluated with a cost linear in lattice size.(Stably converging the PEPS optimization involves some challenges.Further details in [37]).

III. BULK PHASES
Summary of the phase diagram.Fig. 2a shows the bulk phase diagram from Γ-point DMRG with infiniterange interactions.We first discuss the orders identified by their density profiles (orders of some phase transitions are briefly discussed in [37]).Where we observe the same phases as in earlier work [19] we use the same names, although there are very substantial differences with earlier phase diagrams.R b < 1.8.With weaker interactions, the ground states progress through densely-packed, density-ordered phases starting from checkerboard (pink) (R b ∼ 1.2), to striated (cyan) (R b ∼ 1.5), to star (blue) (R b ∼ 1.6).While the checkerboard and star phases are classical-like crystals, the striated state is a density-ordered quantum phase, seen previously [19].R b > 1.8.Here, the phases look very different from earlier work, which truncated the interactions [19].Ordered ground states start with the 1 5 -"staggered" phase (red) (R b ∼ 1.95), then progress to a "nematic" phase (dark green) (R b ∼ 2.2) and the 1  8 -"staggered" phase (gold) (R b ∼ 2.4).There is also a small region at larger δ (not shown) where the nematic phase and a "3-star" classicallike crystal appear to be essentially degenerate, with an energy difference per site of ∆e < 3 • 10 −5 (see [37]).Effects of interactions.In Fig. 2b we show the phase diagram computed using Γ-point DMRG with interactions truncated to distance 2. This approximation resembles earlier numerical studies [19], but here bulk boundary conditions are enforced by the Bloch basis, rather than cylindrical DMRG.Comparing Figs.2a,b we see the disordered and striated phases are greatly stabilized using the full interaction, and new longer-range orders are stabilized at larger R b .Comparing Fig. 2b and Ref. [19], we see that having all atoms interact on an equal footing (via the Bloch basis) destroys the quantum ordered phases seen in [19].Classical, mean-field, and entangled phases.Without the Rabi term Ω, one would obtain classical Rydberg crystals without a disordered phase.Fig. 2c shows the classical phase diagram.For the δ values here, the 1D classical phase diagram has sizable regions of stability for all accessible unit fraction densities [22,26].However, the connectivity of the square lattice in 2D changes this.For example, only a tiny part of the phase diagram supports a 1  3 -density crystal, and we do not find a stable All ordered quantum phases in Fig. 2a appear as classical phases except for the striated and nematic phases, while there are small regions of classical phases at densities 1 3 and 2 9 with no quantum counterpart.The striated and nematic phases emerge near the 1  3 and 1 7 density gaps respectively, however the nematic phase also supersedes the large region of the 1  6 density "3-star" crystal.Ref. [18] suggested that quantum density-ordered phases are qualitatively mean-field states of the form 2d shows the mean-field phase diagram.The disordered phase does not appear, as it emerges from defect hopping and cannot be described without some entanglement [25].The mean-field phase diagram contains features of both the classical and quantum phase diagrams.The striated quantum phase indeed appears as a mean-field state, confirmed by the match between the mean-field and exact correlation functions (Fig. 3a).However, the nematic phase does not appear, and in its place is the same 1  6 -density crystal stabilized in the classical phase diagram.The nematic phase thus emerges as an example of the non-trivial entangled ground-state we are seeking.Nature of the entangled nematic phase.Fig. 3b shows the density correlation function of the nematic phase, which does not display mean-field character.The bi-partite entanglement entropy and entanglement spectrum are shown in Fig. 3c.Importantly, the entanglement spectrum carries 3 large Schmidt values across every cut along the DMRG snake MPS, showing the state is fully entangled across the supercell, and well approximated by an MPS of bond dimension 3.
To reveal the phase structure, Fig. 3d shows the lowest energy classical states in the same region of the phase diagram.Due to the Rydberg blockade radius (R b = 2.3), excitations are spaced by 3 units within a column, giving 3 column configurations |a , |b , |c .Column-column interactions, however, prevent adjacent columns from being in the same configuration (which would have excitations separated by 2 units); this is an adjacent-column constraint.The lowest classical state is the crystal |abc . . .(the 3-star phase) and its 6-fold degenerate permutations, while configurations such as |abab . . .and |abac . . ., which satisfy the inter-and intra-column constraints and are more numerous, lie higher in energy.However, the ground-state character qualitatively changes in the quantum nematic phase.Fig. 3d gives the weights of the configurations in the quantum state.The classically lowest |abc . . .configurations are strongly disfavored, with the state mainly composed of |abab . . .configurations, which allow for greater itinerancy between different column states and thus energy lowering via σx .Note that distant long-range interactions are essential in this analysis, as with truncated interactions all the classical configurations become degenerate.
We can construct an effective one-dimensional spin-1 model of the low-energy sector, and in the thermodynamic limit the model ground-state generically reproduces the entanglement structure seen in Fig. 3c (see [37]).The emergence of one-dimensional entangled order due to kinetic energy and interaction competition recalls other famous nematic phases, for example in 2D fermionic models [32].

IV. FINITE PHASE DIAGRAM
Current experiments are limited to lattices with open boundary conditions consisting of a few hundred atoms [10,18].To investigate how this modifies the bulk behavior, we computed the phase diagram of selected finite lattices from size 9 × 9 to 16 × 16, using DMRG for the smaller sizes and our PEPS methodology for the larger ones.
We first focus (in Fig. 4a) on understanding the fate of the ordered phases on the 15 × 15 lattice along three slices: δ = 2.7, 4.0, and 5.0 (16 × 16 lattice phases, as well as other lattice sizes, are discussed in [37]).Here, many finite lattice ground state orders resemble those in the bulk.However, their regions of stability are substantially reduced and their patterns are broken by frustration.Out of the density-ordered quantum phases, the striated mean-field phase remains due to its commensurate boundary-bulk configurations, while in the region of strongest interactions the nematic phase is destabilized.A new region of classical order, called here the square phase (Fig. 4b), emerges across much of the R b = 1.5−1.8region where the star phase was stable in the bulk [20].We distinguish the square order from the striated order in the sense that the former has negligible quantum fluctuations on the (1, 1)-sublattice, although it is unclear if the square and striated orders constitute truly distinct phases (in the bulk phases the square order is not stable, only the striated order appears).
In Fig. 5, we directly compare the experimental results on the 13 × 13 lattice to our calculations on the same lattice.The analysis of the experiments in Ref. [18] was based on simulations on the 9 × 9 lattice using truncated interactions.This assigned only part of the experimental non-zero order parameter space to a square/striated phase (see Fig. 5a left panel, note, the order parameter does not distinguish between square/striated orders).4), while the (b) row is 13 × 13 numerical data computed in this work.The first two columns show the order parameters used in [18] to identify the square/striated and star phases, while the third column shows a new, more sensitive order parameter for the star phase.Red dots in (a) denote the phase boundaries assigned in [18], while the cyan dotted lines in (b) indicate the subset of parameter space that was computed.Our calculations support a re-interpretation of the experimental data with a significantly larger square/striated region and much smaller star phase.
However, our simulations (Fig. 5b left panel) in fact reproduce the full region of the non-zero order parameter, and thus the whole region seen experimentally should be assigned to a square/striated phase, with the square order appearing in the upper part of the region.Similarly, the experimental analysis identified a large region of star order (Fig. 5a middle panel).This assignment is complicated by edge effects, which mean that the order parameter used does not cleanly distinguish the star phase from other phases.However, our simulations suggest that the region of the star phase should be considered to be much smaller, located at the very top of the non-zero order region, and this is confirmed using a different, more sensitive order parameter (Fig. 5b, right panel).Overall, the measured data corresponds more closely to our numerics than earlier simulations, giving confidence in our more precise interpretation (more discussion in [37]).
Stabilizing entangled ground-state order.Generally, the impact of boundary physics can be understood in terms of frustration of the bulk order by the boundary order, where excitations concentrate more densely due to the lower energetic penalty from fewer long-range interactions on the edge.Examples of the effects of this frustration, ranging from modified bulk orders, to defect dominated states, are shown in Fig. 4b-c (see also [37]).
We searched for conditions to stabilize the entangled nematic ground-state on a finite lattice by manipulating boundary effects.We scanned various rectangular sizes and explicitly "removed" patterns of atoms from the edges to induce different bulk orders.We found the best conditions to stabilize the nematic phase occur near (δ, R b ) = (3.4,2.1), on a 15 × 14 lattice, while removing edge atoms to create a spacing of 4 on two edges and 3 on the other two edges (Fig. 4d) [39].Although there are strong finite size effects, the density profile and correlation functions (Fig. 4d-e) reveal qualitative similarities to the bulk nematic phase, in particular, the presence of 4-fold correlation peaks at distance √ 5 and √ 8, which are also a feature of the bulk entangled phase (Fig. 3b).

V. CONCLUSION
Using new tensor network simulation methods, we have obtained a converged understanding of the phase diagram of Rydberg atom arrays in both bulk and finite simple square lattices.Surprisingly, our bulk phase diagram is quite different from that predicted in earlier numerical studies, while on finite lattices, our results support a reinterpretation of previous experimental analysis.Theoretically, this is due to the subtle effects of the long-range interactions that are addressed by our techniques, while experimentally, it brings into focus the challenge of more accurate theoretical models to interpret increasing experimental capabilities in quantum many-body physics.Perhaps most intriguingly, we find that the unfrustrated square lattice supports an entangled quantum nematic phase, brought about the competition between emergent itinerancy and the constraints of the Rydberg interaction.
A primary focus of Rydberg atom array experiments has been to realize well-studied short-range Hamiltonians, for example, on frustrated lattices.However, we find that lattice frustration is not necessary to produce interesting entanglement in Rydberg systems.In fact our work highlights the richness and complexity intrinsic to Rydberg atom arrays, due to the non-trivial effects of their native interactions.

VII. NUMERICAL METHODS
This section gives details for the numerical simulations in this work.Principally, it will focus on algorithmic subtleties and sources of error, as well as the strategies employed to resolve the physics of the Rydberg atom system.

Convergence and physical strategy
The Γ-point DMRG uses a lattice sum construction of the interaction terms in order to approximate the bulk physics of the Rydberg system (see Eq. 2, main text).As explained in the main text, this is quite different from the cylindrical boundary conditions employed in previous studies [16,17,19] as the represented system is truly infinite.In particular, to increase the range of interactions, we do not need to incur the exponential increase in cost that arises from the associated increase in cylinder width in the standard cylindrical approach.
We find that we can converge our calculations to sufficiently high accuracy with reasonable bond dimensions.Even in the very complicated region of the phase diagram near δ = 5.0 − 6.0, R b = 2.3, we can distinguish the ground-state orders using a bond dimension of D = 1200, as shown in Fig. 6.However, although this is enough to identify the ground state order, higher bond dimensions would be needed to capture the phase transitions with high precision; given the large region of phase space explored here, we leave such detailed calculations to future work.
The strategy used to generate the bulk phase diagram in main text Fig. 2a, as well as the truncated interaction phase diagram Fig. 2b, is as follows.
• For a given point in phase space (δ, R b ), run a D max = 1000 simulation for all reasonable supercell sizes between 4 × 4 and 10 × 10, as well as 12 × 9.
• Identify all supercells for which the ground state has an energy per site within 10 −2 of the lowest energy.
• If there are competing orders, ensure these solutions are all sufficiently converged by requiring (i) the largest singular value truncated during the final DMRG sweeps is less than 10 −8 , and (ii) corrections to the energy when increasing supercell size (up to 12 × 9 maximally) are smaller than the energy gap between competing states (Fig. 6).
• The ground state phase is then identified by evaluating simple density-based order parameters on the largest supercell size which hosts the ground state order.
The only time this convergence criteria is not satisfied is for disordered phase solutions near the order-disorder phase transition (largest truncated DMRG singular value is ∼ 10 −6 ), for which all large supercells show a disordered solution.The classification of the phase in this region is supplemented by analyzing the ground state entanglement entropy, which shows a distinctive "drop" when the phase becomes ordered (see Fig. 7).Importantly, this strategy completely neglects possible orders with unit cells larger than 10 × 10 or 12 × 9, as well as non-periodic solutions.Although orders with unit cells of this large size are not expected in the region of the phase diagram under investigation in this work due to the relatively high crystal densities (and thus close spatial packing) [19,26], our study cannot definitively rule out the stability of such solutions.

Finite size errors
There are a two main sources of finite size error in this formulation of the bulk system.The first comes from the lattice sum form of the long-range interaction.Given the Γ-point Hamiltonian in Eq. 2 (main text), there are some interaction terms of the form ni ni /| r i − r i+R l |.These represent the interaction of a Rydberg excitation with its own periodic "image" in a distant supercell.For a classical crystal, this image term is exact, but in a quantum phase, it is an approximation.Due to the idempotency of n, this term simplifies to ni /| r i − r i+R l |, which would not exist in the Hamiltonian if the supercell were large enough to contain both points i and i + R l .The effect of this error on the energy per site can be estimated by the quantity, Here, ni is the expectation of the local Rydberg excitation for a single characteristic excited site, while ρ ex is the density of sites which have the characteristic excitation of ni .Note that ∆e is always positive, can be systematically reduced by increasing the supercell size, and it is always close to 0 for (almost) classical crystals with excitation densities close to 1 or 0, regardless of cell size.
The other source of systematic error comes from the constraint on the wavefunction imposed by approximating a bulk system by a supercell.Most obviously, this means that certain orders cannot appear in smaller supercell, even for classical crystals.In the case of quantum orders, even for a fixed order there are finite size effects on the emergent kinetic energy of defects.

B. PEPS
The PEPS simulations in this work combine recent advances in optimizing PEPS wavefunctions using automatic differentiation [38] and 2D operator representations of long-range interactions [31].This combination illuminated many new challenges for PEPS optimization with respect to complicated Hamiltonians.This section will detail the various challenges and the technical solutions used in this work.The instability of PEPS optimization remains an open problem and it is an area of future research to determine a PEPS optimization pipeline (using automatic differentiation) that is fully robust to problem instance.In this section, D will refer to the PEPS bond dimension and χ will refer to the maximum bond dimension allowed during contraction before approximations (via SVD) are performed.

Operator representation
The method proposed in Ref. [31] to represent Hamiltonians with long-range interactions writes the interaction potential as a sum of Gaussians, 1 Using the methods in Ref. [43], we can obtain a K = 7 fit with error = max , which is used throughout the work.

Essential computational techniques
As originally discussed in Ref. [38], when trying to use automatic differentiation to optimize a PEPS there are a few essential techniques that must be employed, which are not typically "default" in standard automatic differentiation libraries.They are "essential" in the sense that without them the computation of the energy expectation value and its derivative will typically not run to completion due to out-of-memory errors or numerical infinities.These techniques are: • Numerical stabilization of the gradient of SVD, by adding Lorentzian broadening to the inverse singular values.
• Significant usage of "checkpointing" when evaluating the energy to reduce the memory load of computing gradients.
Both of these techniques are explained in significant detail in Ref. [38].

Stabilizing the optimization
A straightforward implementation of the energy expectation value as described in [31], with optimization via automatic differentiation including the above techniques, typically fails to find the ground state PEPS for the Rydberg Hamiltonian (see Fig. 8).This failure can be generally attributed to the fact that in the quantity under optimization E = ψ|H|ψ ψ|ψ , both the numerator and denominator are evaluated approximately and thus the computation is not strictly bound by the variational principle.Consequently, the optimization can find pathological regions of the PEPS parameter values which make the PEPS contractions inaccurate for the chosen χ, even when starting from an accurately contractible PEPS.Unfortunately, in this problem we find that simply raising the value of χ does not prevent this behavior until χ is impractically large.
In order to mitigate this problem we use the following four techniques in tandem: • We employ line search methods that minimize the gradient norm as well as the energy.In this work, we use the BFGS algorithm [44] in conjunction with such a line search, as suggested in [38].
• We use the cost function where E 1 and E 2 are the energies of PEPS on lattices rotated by 180 degrees and λ is a penalty factor.This strongly penalizes the optimization from entering parameter space with large contraction error (where E 1 and E 2 would be very different).
• During the first iterations of the gradient optimization we only update small patches of tensors at a time, which are chosen to break spatial symmetries that may be contained in the initial guess.After this has pushed the optimization towards the symmetries of the true ground state order, then all tensors can be updated at each optimization step.
• We evaluate the numerator and denominator of E in a consistent way by using "local normalization" during the computation of ψ|H|ψ .This means that, writing H as a comb tensor sum H = Lx i=1 h i , then for each comb tensor numerator ψ|h i |ψ , the associated denominator uses the identical contraction, but with h i replaced by the identity (the environments are not recomputed).
Combining all four of these techniques removes the most egregious instabilities in the optimization trajectory (see Fig. 8), at the cost of a slightly larger computational burden.However, as in more standard DMRG calculations with small bond dimension, convergence to the correct ground-state (rather than a local minimum) still requires a reasonable initial guess.

Initial guess
Obtaining an accurate ground state PEPS typically relies on starting with an accurate initial guess.The predominant algorithms to generate such a guess for problems with a local Hamiltonian are simple update [45][46][47] or imaginary time projection of a converged small D solution to a larger D guess.However, in the presence of long-range interactions it becomes challenging to generalize either of these methods in an efficient and/or accurate way.We therefore used the following simple scheme to generate initial guesses in this work.
• Sum n manually constructed D = 1 PEPS to obtain an initial PEPS of bond dimension D = n.The configurations of these D = 1 PEPS were set to reproduce specific low energy Rydberg crystals and defects within them.
• For small R b : truncate the long-range interactions in H to next-nearest, or next-next-nearest, neighbor interactions (distance of √ 2 or 2), and run conventional simple update starting from the above manually summed PEPS.This fails once the ground state excitations are spaced by more than 2.
• For large R b : add positive random noise to the manually summed PEPS, and then run a highly approximate, first-order gradient optimization for ∼ 25 iterations using a large step size when updating the parameters.

Convergence and physical strategy
Despite the simple procedure to generate initial guesses, we were usually able to systematically converge PEPS solutions according to the conventional protocol of increasing D and χ until the energies corresponding to multiple increasing (D, χ) pairs all vary by less than 0.01% relative to each other.(e.g.see Fig. 10).In this study, we used maximal values of D = 5, χ = 100.However, for a small number of phase points (δ, R b ) we encountered inconsistent convergence of PEPS solutions (see Fig. 9), where increasing D and χ did not systematically result in finding a PEPS with a lower energy, instead getting stuck in various local minima.We attribute this to the low quality of the initial guesses for larger D and R b .
In these cases when PEPS energies could not be systematically converged to within 0.01%, the observed order of the various low-energy solutions were nonetheless the same.The differing energies arose due small quantitative differences such as single-site defects and variations in the local density ni .To further increase certainty in the observed order, we also compared the PEPS solutions to the results of 2D DMRG on the same finite lattice, since the convergence properties of DMRG are much more well-understood.In all cases, the low-energy PEPS solutions had similar energies to the approximate DMRG (relative difference < 1%), and they all showed the same generic low-energy ground state order.The energy gap between phases appeared to be sufficiently large to allow for a tentative classification of the order of this small number of phase points, even though the DMRG was not necessarily converged to high precision (due to the wide lattices) and the PEPS convergence could not be definitively confirmed.The uncertainty in convergence highlights remaining challenges in simulating complex large 2D interacting problems with competing phases using tensor network techniques.The relevant points in the finite lattice phase diagram are labelled by triangles in Fig. 9 above, and in Fig. 4 of the main text.

C. Finite 2D DMRG
Standard 2D DMRG calculations with open boundaries were used to study the 9 × 9 system, a lowentanglement region of the 13×13 system, and to supplement convergence of PEPS on the larger 15 × 14, 15 × 15, and 16 × 16 lattices.Like the PEPS calculations, these too included all long-range interactions (according to Eq. 1 in the main text).The maximal bond dimension used for the 9 × 9 and 13 × 13 simulations was D max = 1200, which we found was more than enough to accurately study the regions of interest in Fig. 5 (main text) for these lattices (see Fig. 11).For supplementing PEPS convergence on the larger lattices, we used D max = 750.Although this bond dimension is not large enough to capture the ground state energy or entanglement of such large systems with high precision, we found it sufficient to capture the first 3-4 digits of the ground state energy and to help with distinguishing between the different lowentanglement ordered phases present in the finite phase diagram, which have substantially larger gaps than the FIG.11.Accuracy of 2D DMRG on the 9 × 9 and 13 × 13 finite lattices (open boundaries).The displayed regions of parameter space correspond exactly to the computed regions in Fig. 5 (main text) and Fig. 15 (SI).The reported error is the largest truncated singular value during the final DMRG sweeps (i.e.once converged).Note that in the ordered regions the error is ∼ 10 −9 , and it grows to ∼ 10 −7 as the ground state becomes disordered on the 13 × 13 lattice due to increasing entanglement.
bulk system due to edge effects.

D. Mean field and classical
The mean field phase diagram for the bulk system (including all long-range interactions) in Fig. 2d (main text) was generated by the following procedure.
• Parameterize the single site wavefunction as , where |0 is the atomic ground state and |1 is the excited Rydberg state.
• Construct a completely un-entangled many-body wavefunction as a typical product of these singlesite states according to all reasonable unit cells between size 2 × 2 and 8 × 10 (supercells are not necessary for mean-field convergence).
• Initialize all possibly relevant configurations for each unit cell as initial guesses.
• Minimize the Γ-point energy for all guesses with respect to the {θ i } using gradient descent.Analytic gradients are easily derived, or automatic differentiation can be employed.
• Classify the phase of the lowest energy state using the same density-based order parameters as the Γpoint DMRG calculations.
The phase space was scanned with a δ-resolution of 0.1 and a R b -resolution of 0.025.Importantly, these calculations are subject to the same limitation as the Γ-point DMRG -they do not capture any possible low energy states with a unit cell larger than 8 × 10.Although such states are not expected in the phase space under examination, this study cannot definitively rule them out.The classical phase diagram for the bulk system (including all long-range interactions) in Fig. 2c (main text) was generated by the following procedure.
• Run classical Monte Carlo minimization of the Γpoint energy for every unit cell size between 2 × 2 and 10 × 10 at phase space points spaced by ∆δ = 0.3, ∆R b = 0.1.
• For all low energy configurations obtained at all phase points, derive their continuous functional form E(δ, R b ) by numerically integrating the interactions.
• Analytically solve for the intersection line between each adjacent pair of configurations in phase space that have minimal energy.
These calculations are also subject to the same limitation as above -any states with unit cells larger than 10 × 10 are not captured, and we cannot rule out their possible existence.

VIII. BULK PHASE DIAGRAM DEGENERATE REGION
In the main text it was briefly mentioned that there is a small region of the bulk phase diagram where the nematic phase and 3-star phase become essentially degenerate.By this we mean that their gap becomes too small to resolve within the estimated finite size error in the Γpoint DMRG numerics.Using the ∆e finite size error

IX. BULK PHASE TRANSITIONS
The order-disorder phase transitions that occur throughout the bulk phase diagram have been characterized as continuous phase transitions in previous work [19].Although full, precise characterization of all bulk phase transitions is beyond the scope of this work, we are able to estimate the order of some transitions using straightforward numerical differentiation of the energies.Figure 13 shows the first and second derivatives of the energy as a function of δ, for various values of R b .The clear peaks in the second derivatives near the critical values of δ support previous conclusions that the disorder→star and disorder→striated phase transitions are indeed second-order.

X. 1D MODEL OF THE BULK NEMATIC PHASE
In the main text the character of the nematic phase was discussed in terms of the classical configurations that make up the quantum wavefunction.It was pointed out that all the low-energy (and thus the most relevant) classical configurations can be described in a succinct notation like |abcabc... in terms of compositions of 3 individual column states |a , |b , and |c which are defined in Fig. 3 of the main text.
This notation is very suggestive of the idea that a nice model for the 2D state can be written as a 1D MPS with a local Hilbert space of dimension 3, spanning |a , |b , The same Hamiltonian can also be written in terms of spin-1 operators and is conceptually straightforward.The first summation in Ĥ1D contains the local terms, where δ D is a direct mapping of the δn term in the original 2D Hamiltonian (up to a scalar), and t T encodes local hopping between the 3 different column states.Microscopically, the hopping emerges virtually from the σx term in the 2D Hamiltonian.The second summation is a direct mapping of the interaction terms in 2D to the new basis, where the various V ij matrices contain the different long-range interaction matrix elements between different pairs of columns.For a faithful mapping, we have V aa = V bb = V cc and V ab = V ac = V bc , and all of them scale as V ij ∼ 1 |i−j| 6 .Using infinite DMRG [48], we can obtain the ground state of this 1D model in the thermodynamic limit.For (δ, R b , t) = (5.0,2.3, 0) we find that the ground state is the unentangled crystal |abcabc... , as expected since this corresponds to the the classical ground state in 2D and t = 0 is the classical limit of this model.For t > 0 we find an entangled ground state with equal density in the local |a , |b , and |c basis states, as in the 2D nematic ground state.We also observe that the structure of the entanglement spectrum in the 1D model ground state is very similar to the entanglement spectrum of the 2D nematic state in between the columns, as shown in Fig. 14.We conclude that this 1D effective Hamiltonian provides a useful model of the 2D nematic state.The phase diagram of the 15 × 15 lattice reported in the main text contained many of the ground state orders seen in the bulk phase diagram, but it also revealed the strong finite-size effects induced by the boundary.Due to the long-range van der Waals interactions, Rydberg excitations at the edge of the array incur roughly half of the energetic penalty that excitations in the interior do, but lower the energy by an equal amount (δ).Except at small values of R b , this induces excitations along the edge of the array to be more densely packed than what would be expected from the bulk phase diagram at a given point (δ, R b ).This generic effect causes frustration between the boundary and interior of the finite lattices, which gives rise to the square classical order and many defectdominated states at large R b , as discussed in the main text.In these defect states, the optimal bulk density becomes so small relative to the optimal edge density that the ground states are permeated by edge-induced defects, leaving only small regions of any discernible order and making the precise configuration very sensitive to small changes in R b and δ.
In addition to the 15 × 15 lattice, we also studied two slices (δ = 4.0, 5.0) of the phase diagram of the 16 × 16 lattice to probe for bulk-like ordered phases where the 15 × 15 system is dominated by defects.Specifically, we focused on the R b > 1.8 region, for which the results are shown in Fig. 9b.We find a clear region of the stability for the boundary-bulk frustrated 1  5 -stagger phase (red), for which the density profile is shown in Fig. 4c of the main text.Along with a small region of the 3-star phase (green and black), these regions are unique to the 16 × 16 lattice (i.e. they are not seen in 15 × 15).There are also some common features between the two array sizes, namely regions of the star and 1 8 -stagger (gold) phase as well as many defect states.This suggests that the defect states are an intrinsic part of the physics of medium-sized arrays.
As reflected by the triangular markers in Fig. 9b (which reflect inconsistent convergence) we found it more challenging than the 15×15 lattice to systematically converge the PEPS calculations with respect to (D, χ), especially shows results using interactions truncated to zero beyond distance 2. This is identical to the truncation scheme used in numerics in Ref. [18].The first three columns show all three order parameters used in [18] to distinguish the phase diagram, while the fourth column shows a new, more precise order parameter for the star phase.Red dots in (a) denote the phase boundaries assigned in [18], while the cyan dotted lines in (b)-(d) indicate the subset of parameter space that was computed.
in the star phase (blue).In part, this was due to the boundary itself being frustrated; on an even-sided lattice it is not possible to place excitations in all corners and also along all edges spaced by a distance of 2. Because the corner excitations are strongly pinned due to their reduced interaction penalty, this causes the boundary to be frustrated and makes it more difficult to prepare a good initial guess with our rudimentary strategies.The main text discussed discrepancies between our numerical results on the 13 × 13 lattice and analysis reported in a recent experiment [18], specifically concerning the striated, square and star phases.It was noted that the actual experimental data appears to agree with our numerics, but the interpretation of the data offered in Ref. [18] is inconsistent with ours.This section details the effect of the approximations made in the numerics of Ref. [18] on the interpretation of the data, and how relaxing those approximation leads to the interpretation described in our main text.

A. Context
In Ref. [18], the experimental data on the 13 × 13 square lattice was primarily understood with respect to DMRG calculations performed on the 9 × 9 lattice (all open boundaries), in which interactions were truncated to zero beyond a distance of 2. The experimental results of Ref. [18] are reproduced in Fig. 15a, and they are compared to our numerical results on 9 × 9 and 13 × 13 lattices (Fig. 15b-d).The region of the phase diagram that was studied included domains of stability for the disordered, checkerboard, striated, and star phases.The square "phase" was not separately reported, although it may be considered the classical limit of the striated phase.
We also introduce here a useful order parameter for detecting the star phase, where N = L x • L y .O star detects a symmetry breaking that occurs in the star phase but not in the disordered, checkerboard, striated, or square phases.On a finite lattice, this provides a clean way to define the star phase separate from the other orders in this set.We also recapitulate the definition of the order parameters defined in [18] and used in Fig. 15 In Fig. 15c, we recompute the main 9×9 phase diagram numerical results used in [18], which use truncated interactions.The bright region in F(π, π/2) predicts a large domain of stability for the star phase, which is corroborated by the value of O star .This data was used in [18] to draw the expected phase boundary in the 13 × 13 experimental data seen in Fig 15a .However, Fig. 15d shows the analogous results on the 9 × 9 lattice when including all long-range interactions.Surprisingly, the star phase gets completely destabilized!This illustrates the hazard of interpreting the experimental data from smaller lattice simulations.
Unlike the 9 × 9 lattice, we observe that the 13 × 13 lattice phase diagram has a qualitative difference: it hosts a nonzero domain of star phase even when accounting for all long-range interactions.As pointed out in the main text, F(π, π/2) is not a sensitive order parameter for the star phase as it appears on finite lattices, but O star does reveal the tiny stable region of the star phase (see Fig. 15b).

C. Square and striated phases
The overestimation of the extent of the star phase by using numerics from the 9×9 lattice with truncated interactions also results in an underestimation of the extent of the striated order parameter, F(π, 0) − F(π/2, π), since F(π/2, π) is the star order parameter used in [18] (see Fig 15c).These 9×9 results were used in [18] to interpret the striated phase domain in the experimental data, so the boundary drawn in Fig. 15a is too small.In fact, the extent of the experimental data for F(π, 0) − F(π/2, π) (Fig. 15a) is significantly larger than the drawn boundary, corresponding much more closely to the numerical data on the 13 × 13 including long-range interactions (Fig. 15b), as mentioned in the main text.
In this work, we distinguish a region of classical square order from the striated phase where the square order contains (almost) no quantum fluctuations on the (1, 1)sublattice, which are an essential feature of the striated phase in the bulk.F(π, 0) − F(π/2, π) does not help distinguish between square and striated orders, and no classical square order was reported in Ref. [18].In Fig. 16 we show that a large part of the bright region in F(π, 0) − F(π/2, π) on the 13 × 13 lattice should be interpreted as a classical square order by plotting, x,y n x,y ifx mod 2 = 1, y mod 2 = 1 0 else which detects the deformation of the density on the (1, 1)sublattice.This sublattice is defined in terms of a 2 × 2 cell, as in [18].

D. Numerical accuracy
All numerical results in Figs.15-16 were computed using DMRG.It was possible to study the 13 × 13 lattice using DMRG because we only investigated a lowentanglement region of the phase diagram.The level of accuracy for these calculations is shown in Fig. 11 in terms of the largest truncated singular value during the DMRG sweep.In the ordered regions of the results, the largest truncated singular value is below 10 −9 , which is generally considered accurate.

1 7 -
density crystal within unit cell sizes of up to 10 × 10.

FIG. 2 .
FIG. 2. Phase diagrams of the bulk system under various assumptions.The color of a dot/region identifies the ground state order.The density profiles for each color are given in (e) and shown near each phase domain.(a) The phase diagram given by Γ-point DMRG including all long-range interactions.(b) The phase diagram from Γ-point DMRG when interactions are truncated to 0 beyond a distance of | ri − rj| = 2. (c) The classical phase diagram (when all sites are either fully occupied or empty) including all long-range interactions.(d) The mean-field phase diagram, including all long-range interactions.Error bars display the uncertainty of the computed phase boundaries.(e) Representative density profiles for all phases in (a)-(d), identified by the colored dot in each lower right corner.All profiles have Γ-point boundary conditions on all edges.In (a)-(b) dots denote computed data, while shading is a guide for the eye.(c),(d) are computed with very fine resolution/analytically, thus no dots are shown.

FIG. 3 .
FIG. 3. Mean-field striated versus entangled nematic phase.(a) Density-density correlation functions of the mean-field and exact striated ground state, both at (δ, R b ) = (3.1,1.5); these agree, confirming the mean-field nature of the striated phase.(b) Density-density correlation functions for the entangled nematic phase ground state and two different mean-field ground states (from a 6×3 unit cell and a 3×4 unit cell) at (δ, R b ) = (5.0,2.3).In (a)-(b), 2-fold/4-fold degeneracy of a peak is indicated by 2/4 horizontal dots distributed around the proper distance coordinate.8-fold degeneracy in (a) is shown as two rows of 4 dots.The non-mean-field (entangled) character of the nematic phase is evident.(c) Bipartite entanglement entropy for each possible bipartition of the 12 × 9 supercell nematic ground state.One inset shows the "path" that the "partition location" axis follows through the supercell MPS, while the other shows the entanglement spectrum at a central cut.(d) Structure of the nematic state in terms of classical configurations constructed via compositions of 3 individual column states |a , |b , |c .In the classical limit, there are 4 distinct sets of low-energy configurations, all characterized by the absence of adjacent columns in the same state (e.g.|aa... ) and large degeneracies due to permutational symmetry between |a , |b , and |c .The lowest in energy is 6-fold degenerate, corresponding to the 3-star state.However, in the quantum nematic state the configurations that are slightly higher in energy have much larger wavefunction coefficients.The most relevant classical states in the wavefunction are those with the greatest number of possible single "column hops" (e.g. a → b) without introducing unfavorable states like |aa... , revealing the role of itinerancy in the nematic phase.

FIG. 4 .
FIG. 4. Phase diagram of the 15 × 15 finite system and finite lattice orders.(a) The phase diagram, where colors correspond to the same phase classifications as Fig. 2. Triangles represent tentative classification of points showing inconsistent PEPS convergence, see [37].A new "square" order is specified in (b) and various examples of boundary-bulk frustrated ground states in (c).(d) The density profile for a nematic-like ground state that can be stabilized on a 15 × 14 lattice at (δ, R b ) = (3.4,2.1) with manually tailored edge excitations (see text).(e) Comparing the correlations of the finite nematic phase to the "exact" bulk phase.The degeneracy of the peaks is split by the boundary excitations, but the number of peaks is generally conserved between the two (green ovals), which provides a clear distinction from mean-field states (see Fig.3b).

FIG. 5 .
FIG. 5. Comparison to experiment.The (a) row directly reproduces experimental phase diagram data on the 13 × 13 lattice (data extracted from Ref. [18] Fig.4), while the (b) row is 13 × 13 numerical data computed in this work.The first two columns show the order parameters used in[18] to identify the square/striated and star phases, while the third column shows a new, more sensitive order parameter for the star phase.Red dots in (a) denote the phase boundaries assigned in[18], while the cyan dotted lines in (b) indicate the subset of parameter space that was computed.Our calculations support a re-interpretation of the experimental data with a significantly larger square/striated region and much smaller star phase.
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA (Dated:January 11, 2022)

FIG. 6 .
FIG. 6. Convergence of Γ-point DMRG in the most difficult region of the phase diagram (δ, R b ) = (5.0− 6.0, 2.3).(a) Shows the convergence w.r.t.bond dimension of the largest truncated DMRG singular value (red) and the change in energy per site relative to the energy obtained with bond dimension D − 200 (blue).(b) The energies per site of a large variety of supercell sizes.This is adapted from Fig. 1 of the main text to highlight the relevant points.The connected dark green points are the nematic phase, and lime green points are the low energy 3-star 1 6 -density crystalline phase.The inset shows the convergence of the nematic phase energy w.r.t.supercell size and gaps to the other low energy solutions, whose density profiles are shown in (c).Note that, between (a) and (b), the nematic phase is converged to below 10 −5 accuracy while the competing states differ in energy by at least 10 −4 .

FIG. 7 .
FIG.7.Bipartite entanglement entropy of various crystalline phases as δ increases.Each line is a slice over δ values for a fixed R b value.Black line segments denote when the ground state is in the disordered phase.Solid colored line segments denote when the ground state is an ordered crystalline phase (same color classifications as the phase diagram in main text).Dotted line segments denote the "transition zone" of a given line between the disordered phase and an ordered phase.These are simply a result of the finite resolution used to sample phase space in the phase diagrams.

FIG. 8 .
FIG.8.Examples of typical optimization trajectories for longrange PEPS using different automatic differentiation schemes.The blue line often occurs with a naive implementation of the energy evaluation algorithms and use of a line search which does not minimize gradient norm.The red line can occur even when using a more sophisticated energy evaluation including local norms and/or a multi-evaluation cost function.The stable magenta and green lines result from combining the four techniques discussed in Sec.VII B 3. The difference between the magenta and green curves reflects the quality of the initial guess.

FIG. 9 .
FIG. 9. Phase diagrams of 15 × 15 (a) and 16 × 16 (b) arrays, detailing convergence.Circular points indicate systematic convergence with PEPS up to bond dimension D = 5, while triangles indicate intermittent convergence with PEPS, requiring supplemental convergence checks using 2D DMRG.More details are available in Sections XI and VII B. Colors used in these plots correspond identically to the colors used in the main text to identify phases.

FIG. 12 .
FIG. 12. Expanded view of the large-R b part of the bulk phase diagram, computed with Γ-point DMRG.For δ ≤ 5.0, this data is identical to Fig. 2a in the main text.All colors correspond to the same phases as in the main text Fig. 2. The small lime green region indicates the degenerate zone where the gap between the 3-star and nematic phases becomes very small.

FIG. 13 .
FIG.13.Numerical evidence of second-order phase transitions between the disordered phase and the star (blue) and striated (cyan) phases.Top: First derivative of the energy with respect to δ. Bottom: Second derivative of the energy with respect to δ.Both are estimated using standard finite difference formulas.

FIG. 14 .
FIG.14.The inter-column entanglement spectrum of the 2D nematic phase (left) and the ground state of the 1D effective model presented in section X (right) for hopping parameter t = 0.5.They both contain 3 dominant eigenvalues with one doubly degenerate pair.The subsequent levels consist of two double degeneracies and two single degeneracies in both cases, although the spacing and ordering is slightly different.

FIG. 15 .
FIG. 15.Detailed comparison to experimental phase diagram.The (a) row directly reproduces the experimental phase diagram on the 13 × 13 lattice (data extracted from Ref. [18] Fig. 4).Rows (b)-(d) show analogous numerical data on 9 × 9 and 13 × 13 lattices, where (b) and (d) are results from simulations containing all long-range interactions and (c)shows results using interactions truncated to zero beyond distance 2. This is identical to the truncation scheme used in numerics in Ref.[18].The first three columns show all three order parameters used in[18] to distinguish the phase diagram, while the fourth column shows a new, more precise order parameter for the star phase.Red dots in (a) denote the phase boundaries assigned in[18], while the cyan dotted lines in (b)-(d) indicate the subset of parameter space that was computed.