2201.03189.pdf

Entanglement in the quantum phases of an unfrustrated Rydberg atom array

Matthew J. O’Rourke and Garnet Kin-Lic Chan

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA

(Dated: January 11, 2022)

We report on the ground state phase diagram of interacting Rydberg atoms in the unfrustrated

square lattice array. Using new tensor network algorithms, we scale to large systems in two dimen-

sions while including all long-range interactions, revealing the phases in the bulk and their analogs

in accessible finite arrays. We find a greatly altered phase diagram from earlier numerical and ex-

perimental studies, and in particular, we uncover an emergent entangled quantum nematic phase

that appears in the absence of frustration. 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 Ry-

dberg states [1, 2], have generated interest for quantum

information processing and to realize exotic many-body

states [3–17]. A recent experiment [18], backed by nu-

merical 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 entangle-

ment [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 interac-

tions greatly affect the phases, complicating the accurate

numerical determination of bulk behaviour. Here, we de-

scribe new numerical techniques that greatly reduce fi-

nite 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 interac-

tions. Unexpectedly, we derive quite different physics

from our simulations compared to both previous theoret-

ical and experimental analyses – including the emergence

of a non-trivial, entangled nematic phase, even on the un-

frustrated square lattice array.

The Rydberg atom array Hamiltonian is

∑

[

Ω

−

]

∑

(

−

)

(1)

Here ˆ

〉〈

and ˆ

〉〈

(

〉

〉}

denote ground and Rydberg states of atom

is lat-

tice spacing, Ω labels Rabi frequency, and

describes

laser detuning.

parameterizes the interaction strength

between excitations. This can be re-expressed in terms of

the Rydberg blockade radius

, with

(

)

≡

Ω.

We study the square lattice in units

= Ω = 1 [19],

yielding two free parameters

and

The ground states of this Hamiltonian are simply un-

derstood in two limits. For

δ/

Ω

= 0, the system

is classical and one obtains classical crystals of Rydberg

excitations [21–24] whose spatial density is set by the

competition between

and

. For

δ/

Ω

= 0,

Rydberg excitations are disfavored and the solutions are

dominated by Rabi oscillations, leading to a trivial “dis-

ordered” 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 Lut-

tinger liquid appearing on the boundary between ordered

and disordered phases [26].

In 2D, however, the picture is quite different. An ini-

tial study [19] using the density matrix renormalization

group (DMRG) [27–30] found additional quantum crys-

talline (or “density-ordered”) phases, where the local ex-

citation density is not close to 0 or 1. A recent exper-

iment on a 256 programmable atom array has realized

such phases [18]. However, as also discussed there, the

density-ordered phases are unentangled quantum mean-

field phases, and thus not very interesting. In addition,

more recent numerical results [17] highlight the sensitiv-

ity of the physics to the tails of the Rydberg interaction

and finite size effects. Thus, whether Rydberg atom ar-

rays on a simple unfrustrated lattice – such as the square

lattice – support interesting quantum ground-states, re-

mains an open question.

Here, we resolve these questions through high-fidelity

numerical simulations.

To do so, we employ varia-

tional tensor network methods. Tensor networks have

led to breakthroughs in the understanding of 2D quan-

tum many-body problems [32], and we rely on two new

techniques that address specific complexities of simulat-

ing 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 in-

teraction truncations and boundary effects present in ear-

lier studies [16, 17, 19, 20]. This allows us to controllably

converge the bulk phase diagram. The second is a repre-

sentation of long-range interactions [31] compatible with

projected entangled pair states (PEPS) [33–36]. With

this, we use PEPS to find the ground states of a Hamil-

tonian 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-

arXiv:2201.03189v1 [cond-mat.str-el] 10 Jan 2022

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

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

comb tensor network operators. Tensors of the same color are identical.

ods, before turning to the bulk and finite-size phase be-

havior 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 pre-

vious studies have employed a cylindrical DMRG geom-

etry 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 cal-

culations in a Bloch basis. The resulting Γ-point 2D

DMRG formally models an infinite lattice (Fig. 1a) with

a wavefunction constrained by the supercell, and peri-

odic boundary conditions in both directions. This differs

from using a periodic matrix product state (MPS) as pe-

riodicity 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 fi-

nite; here the lattice remains infinite. The Bloch basis

interactions enter as an infinite lattice sum over lattice

vectors

= (

);

n,m

∈

Z

;

∑

[

−

]

∑

−

(2)

where

are the supercell side lengths.

The only finite size parameter is the supercell size

. We thus perform exhaustive scans over

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

−

, compared to the smallest energy den-

sity difference we observe between competing phases of

∼

−

(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 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 en-

tanglement 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 con-

sistency between the two provides strong confirmation.

However, PEPS are usually combined with short-range

Hamiltonians. We now discuss a way to combine long-

range Hamiltonians efficiently with PEPS without trun-

cations.

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.

∑

max

−

(where

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 eval-

uated 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 infinite-

range 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, al-

though there are very substantial differences with earlier

phase diagrams.

. With weaker interactions, the ground states

progress through densely-packed, density-ordered phases

starting from checkerboard (pink) (

∼

2), to striated

(cyan) (

∼

5), to star (blue) (

∼

6). While the

checkerboard and star phases are classical-like crystals,

the striated state is a density-ordered quantum phase,

seen previously [19].

. Here, the phases look very different from ear-

lier work, which truncated the interactions [19]. Ordered

ground states start with the

-“staggered” phase (red)

(

∼

95), then progress to a “nematic” phase (dark

green) (

∼

2) and the

-“staggered” phase (gold)

(

∼

4). There is also a small region at larger

(not

shown) where the nematic phase and a “3-star” classical-

like crystal appear to be essentially degenerate, with an

energy difference per site of ∆

e <

−

(see [37]).

Effects of interactions

. In Fig. 2b we show the phase

diagram computed using Γ-point DMRG with interac-

tions truncated to distance 2. This approximation resem-

bles earlier numerical studies [19], but here bulk bound-

ary 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 us-

ing the full interaction, and new longer-range orders are

stabilized at larger

. Comparing Fig. 2b and Ref. [19],

we see that having all atoms interact on an equal foot-

ing (via the Bloch basis) destroys the quantum ordered

phases seen in [19].

Classical, mean-field, and entangled phases

. With-

out 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]. How-

ever, the connectivity of the square lattice in 2D changes

this. For example, only a tiny part of the phase diagram

supports a

-density crystal, and we do not find a stable

-density crystal within unit cell sizes of up to 10

10.

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

and

with no quantum counterpart. The striated and

nematic phases emerge near the

and

density gaps

respectively, however the nematic phase also supersedes

the large region of the

density “3-star” crystal.

Ref. [18] suggested that quantum density-ordered

phases are qualitatively mean-field states of the form

∏

〉

√

−|

〉

. Fig. 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 quan-

tum 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

ap-

pear, and in its place is the same

-density crystal stabi-

lized 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 spec-

trum are shown in Fig. 3c. Importantly, the entangle-

ment spectrum carries 3 large Schmidt values across ev-

ery cut along the DMRG snake MPS, showing the state

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

−

= 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.

is fully entangled across the supercell, and well approxi-

mated 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 di-

agram. Due to the Rydberg blockade radius (

= 2

3),

excitations are spaced by 3 units within a column, giving

3 column configurations

〉

. Column-column in-

teractions, however, prevent adjacent columns from being

in the same configuration (which would have excitations

separated by 2 units); this is an adjacent-column con-

straint. 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 quan-

tum nematic phase. Fig. 3d gives the weights of the con-

figurations in the quantum state. The classically low-

est

abc...

〉

configurations are strongly disfavored, with

the state mainly composed of

abab...

〉

configurations,

which allow for greater itinerancy between different col-

umn states and thus energy lowering via ˆ

. Note that

distant long-range interactions are essential in this anal-

ysis, as with truncated interactions all the classical con-

figurations become degenerate.

We can construct an effective one-dimensional spin-

1 model of the low-energy sector, and in the thermo-

dynamic limit the model ground-state generically re-

produces 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