Article
https://doi.org/10.103
8/s41467-023-41166-0
Entanglement in the quantum phases of an
unfrustrated Rydberg atom array
Matthew J. O
’
Rourke
1
& Garnet Kin-Lic Chan
1
Recent experimental advances have stimulated interest in the use of large,
two-dimensional arrays of Rydberg atoms as a platform for quantum infor-
mation processing and to study exotic
many-body quantum states. However,
the native long-range interactions between the atoms complicate experi-
mental 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
fi
nd a greatly altered phase
diagram from earlier numerical and experimental studies, revealed by study-
ing the phases on the bulk lattice and t
heir analogs in experiment-sized
fi
nite
arrays. We further describe a previousl
y unknown region with a nematic phase
stabilized by short-range entanglement and an order from disorder mechan-
ism. Broadly our results yield a conceptu
al guide for future ex
periments, while
our techniques provide a blueprint for converging numerical studies in other
lattices.
Rydberg atom arrays consist of a set of cold neutral atoms that are
trapped in an optical lattice, interacting strongly with each other via
excitation into Rydberg states
1
,
2
. Advances in experimental control
over a large number of atoms, arranged in two-dimensional arrays,
have generated signi
fi
cant interest in using these systems for a variety
of applications, including quantum information processing and stabi-
lizing quantum states with long-range entanglement
3
–
17
.Arecent
seminal experiment
18
, backed by numerical studies
19
,
20
,hassuggesteda
richness in the ground states of Rydberg atom arrays on a 2D square
lattice. However, although the observed, non-disordered, phases are
not all classical crystals, they contain little entanglement
19
.Thusit
remains unclear whether such arrays realize non-trivial entangled
quantum ground-states on simple lattices.
The Rydberg atom array Hamiltonian is
^
H
=
X
N
i
=1
Ω
2
^
σ
x
i
δ
^
n
i
+
1
2
X
i
≠
j
V
ðj
r
i
r
j
j
=
a
Þ
6
^
n
i
^
n
j
:
ð
1
Þ
Here
^
σ
x
i
=0
i
1
i
+1
i
0
i
and
^
n
i
=1
i
1
i
(
f
0
i
,1
i
g
denote ground
and Rydberg states of atom
i
).
a
is lattice spacing,
Ω
labels Rabi fre-
quency, and
δ
describes laser detuning.
V
parameterizes the interac-
tion strength between excitations. This can be re-expressed in terms of
the Rydberg blockade radius
R
b
,with
V
=
ð
R
b
=
a
Þ
6
Ω
.Westudythe
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
–
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
.Inbetween
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
–
30
found
additional quantum crystalline (or so-called 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 mean-
fi
eldphases,andthusnotveryinter-
esting. In addition, more recent numerical results
17
highlight the sen-
sitivity of the physics to the tails of the Rydberg interaction and
fi
nite
size effects. Thus, whether Rydberg atom arrays on a simple
Received: 17 February 2022
Accepted: 23 August 2023
Check for updates
1
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA.
e-mail:
gkc1000@gmail.com
Nature Communications
| (2023) 14:5397
1
1234567890():,;
1234567890():,;
unfrustrated lattice
—
such as the square lattice
—
support interesting
quantum ground-states, remains an open question.
Here, we resolve these questions through high-
fi
delity numerical
simulations. To do so, we develop and employ new numerical techni-
ques based on variational tensor network methods. Tensor networks
have led to breakthroughs in the understanding of 2D quantum many-
body problems
31
, and our two new techniques address speci
fi
c com-
plexities of simulating interactions in Rydberg atom arrays. The
fi
rst we
term
Γ
-point DMRG, which utilizes a computational spin basis with
periodic boundary conditions, and which can also be viewed as a type
of DMRG that is deployed on a torus with interactions wrapping
around to in
fi
nite range, while employing a traditional
fi
nite system
two-dimensional DMRG methodology
27
. This removes interaction
truncations and boundary effects present in earlier studies
16
,
17
,
19
,
20
,and
allows us to controllably converge the bulk phase diagram. The second
is a representation of long-range interactions
32
compatible with pro-
jected entangled pair states (PEPS)
33
–
36
. With this, we use PEPS to
fi
nd
the ground states of a Hamiltonian with long-range interactions for the
fi
rst time, and speci
fi
cally here, model the states of
fi
nite Rydberg
arrays of large widths as used in experiment. We show that, unex-
pectedly, the faithful inclusion of all long-range terms in our simula-
tions yields quite different phy
sics compared to both previous
theoretical and experimental analyses. Some previously predicted
ground state phases are destabilized, while other unanticipated phases
emerge
–
including evidence of a non-trivial nematic phase stabilized
by short-range entanglement in an order from disorder mechanism
37
,
even on the geometrically unfrustrated square lattice array. In the
following, we
fi
rst describe the new numerical techniques, before
turning to the bulk and
fi
nite-size phase behavior of square lattice
Rydberg arrays and the question of entangled quantum phases.
Results
Bulk simulation strategy 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
,themain
fi
nite size effects arise from interactions. Many
previous studies have employed a cylindrical DMRG geometry com-
monin2DDMRGstudies
27
. However, there the interaction is truncated
to the cylinder half-width, while along the open direction, edge atoms
experience different interactions than in the bulk; both choices pro-
duce strong
fi
nite size effects.
To avoid these problems, we perform 2D DMRG calculations in a
site Bloch basis. Given the site basis
j
n
x
,
y
i
,
n
∈
{0, 1}, the Bloch basis
states are periodic combinations,
j
~
n
x
,
y
i
=
P
l
j
n
ð
x
,
y
Þ
+
R
l
i
where
R
l
=(
n
⋅
L
x
,
m
⋅
L
y
),
n
,
m
2
Z
, are lattice vectors,
L
x
,
L
y
are the supercell
side lengths, and
n
x
,
y
=
n
ð
x
,
y
Þ
+
R
l
, i.e., the occupancies at the transla-
tionally related sites are the same. The above are analogous to Bloch
states at the
Γ
-point of the supercell Brillouin zone. The
fi
nite many-
body Hilbert space under the
Γ
-point restriction is
Q
x
,
y
j
~
n
x
,
y
i
; this Hil-
bert space should be interpreted as a model of the Hilbert space of the
in
fi
nite system, rather than a true subspace of it. The corresponding
matrix product state (MPS) is expressed as
Ψ
ji
=
P
f
e
g
Q
x
,
y
A
~
n
x
,
y
f
e
x
,
y
g
j
~
n
x
,
y
i
where
A
~
n
x
,
y
is the MPS tensor associated with Bloch function
~
n
x
,
y
,
e
x
,
y
denote its bonds, and a 2D snake ordering has been chosen through
the lattice. In the above picture,
Γ
-point 2D DMRG formally models an
in
fi
nite lattice (Fig.
1
a) with a wavefunction constrained by the super-
cell shape. This differs from using a periodic MPS as periodicity is
enforced by the Bloch basis rather than the MPS. The
Γ
-point DMRG
calculation can also be viewed as working on a
fi
nite toroidal geometry
(i.e., the supercell) with the typical pure site basis
j
n
x
,
y
i
,butwherethe
interactions are allowed to wrap in
fi
nitely around the torus, rather than
being cut off. In either interpretati
on, the Hamiltonian per supercell
contains interactions expressed as an in
fi
nite lattice sum,
^
H
=
X
i
1
2
^
σ
x
i
δ
^
n
i
+
1
2
X
i
≠
j
+
R
l
,
R
l
R
6
b
j
r
i
r
j
+
R
l
j
6
^
n
i
^
n
j
,
ð
2
Þ
Further details of this approach and its two interpretations are dis-
cussed in the Methods section.
The only
fi
nite size parameter is the supercell size
L
x
×
L
y
.Wethus
perform exhaustive scans over
L
x
,
L
y
to identify competing ground
state orders. We systematically converge the energy per site of low-
a
b
c
d
68
69
77
87
88
89
97
99
10 8
10 9
10 10
12 9c
12 9q
Supercell size
0.0
0.01
0.02
Energy per site
68
89
12 9q
-10
3
-10
5
0
+10
5
+10
4
ground state order
0 1 0 1
0 1 0 1
1 0 1 0
1 0 1 0
Fig. 1 | Numerical methods and strategy. a
A schematic representation of
Γ
-point
DMRG. A single in
fi
nite bulk con
fi
guration is given by periodic images of the central
supercell con
fi
guration. The wavefunction coef
fi
cient for this in
fi
nite con
fi
guration
is given by the contraction of a snake MPS, which is de
fi
ned 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.
fi
nite 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
fi
nite
system. Each tensor is a different color because they can all be unique.
d
Asim-
pli
fi
ed diagrammatic representation of the long-range Hamiltonian construction
for PEPS in ref.
32
. All terms in the Hamiltonian are accounted for by a sum of
L
x
comb tensor network operators. Tensors of the same color are identical.
Article
https://doi.org/10.103
8/s41467-023-41166-0
Nature Communications
| (2023) 14:5397
2
energy orders by increasing the commensurate supercell sizes to
contain many copies of the order (up to 108 sites). The
fi
nite size
effects converge rapidly because no interactions are truncated and
there are no edge effects even in the smallest cells, allowing us to
converge the energy per site to better than 10
−
5
,comparedtothe
smallest energy density difference we observe between competing
phases of ~10
−
4
(see Fig.
1
b and Supplementary Methods).
Finite simulations and PEPS with long-range interactions
To simulate ground-states of
fi
nite arrays, we consider
fi
nite 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 con
fi
dently 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.
1
c). Together with DMRG cal-
culations on moderate width
fi
nite lattices, the two methods provide
complementary approaches to competing phases and consistency
between the two provides strong con
fi
rmation. However, PEPS are
usually combined with short-range Hamiltonians. We now discuss a
way to combine long-range Hamiltonians ef
fi
ciently with PEPS without
truncations.
For this, we rely on the representation we introduced in ref.
32
.
This encodes the long-range Hamiltonian as a sum of comb tensor
network operators (Fig.
1
d). As discussed in ref.
32
, arbitrary isotropic
interactions can be ef
fi
ciently represented in this form, which mimics
the desired potential via a sum of Gaussians, i.e.
1
r
6
=
P
k
max
k
=1
c
k
e
b
k
r
2
(where
k
max
~ 7 for the desired accuracy in this work). The combs can
be ef
fi
ciently contracted much more cheaply than using a general
tensor network operator.
While ref.
32
described the Hamiltonian encoding, here we must
also
fi
nd the ground-state. We variationally minimize
h
Ψ
j
^
H
j
Ψ
i
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. Further details are discussed
in the Methods section, including some challenges in stably conver-
ging the PEPS optimization.
Summary of the bulk phase diagram
Figure
2
a shows the bulk phase diagram from
Γ
-point DMRG with
in
fi
nite-range interactions. We
fi
rst discuss the orders identi
fi
ed by
their density pro
fi
les (orders of some phase transitions are brie
fl
y
discussed in Supplementary Note 3). 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.
With weaker interactions (
R
b
< 1.8), 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
.
With stronger interactions (
R
b
> 1.8), the phases look very differ-
ent 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)andthe
1
8
-staggered
phase (gold,
R
b
~ 2.4). There is also a small region at larger
δ
(not
1
2
3
4
5
1.2
1.4
1.6
1.8
2.0
2.2
2.4
1
2
3
4
5
1.2
1.4
1.6
1.8
2.0
2.2
2.4
1
2
3
4
5
1.2
1.4
1.6
1.8
2.0
2.2
2.4
1
2
3
4
5
1.2
1.4
1.6
1.8
2.0
2.2
2.4
ab
e
c
d
full interactions
truncated interactions
mean-field
classical
1/2
1/4
1/5
1/6
1/8
2/9
1/3
disordered
checkerboard
striated
star
1/5-stagger
nematic
1/8-stagger
modulated
square
1/3-stagger
diamond
3-star
1.0
0.0
Fig. 2 | Phase diagrams of the bulk system under various assumptions.
The color
of a dot/region identi
fi
es the ground state order. The density pro
fi
les 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
∣
r
i
−
r
j
∣
=2.
c
The classical phase diagram (when all sites are either fully occupied or empty)
including all long-range interactions.
d
The mean-
fi
eld phase diagram, including all
long-range interactions. Error bars display the uncertainty of the computed phase
boundaries.
e
Representative density pro
fi
les for all phases in (
a
–
d
), identi
fi
ed by
the colored dot in each lower right corner. All pro
fi
les 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
fi
ne resolution/analytically, thus no
dots are shown.
Article
https://doi.org/10.103
8/s41467-023-41166-0
Nature Communications
| (2023) 14:5397
3
shown) where the nematic phase and a classical-like crystal (which we
call 3-star) appear to be essentially degenerate, with an energy differ-
ence per site of
Δ
e
<3
⋅
10
−
5
(see Supplementary Note 2).
Effects of interactions on the bulk phases
In Figure
2
b, we show the phase diagram computed using
Γ
-point
DMRG with interactions truncated to distance 2. This approxima-
tion resembles earlier numerical studies
19
, but here bulk boundary
conditions are enforced by the Bloch basis, rather than cylindrical
DMRG. Comparing Fig.
2
a, 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.
2
b
and ref.
19
, we see that having all atoms interact on an equal footing
(via the Bloch basis) destroys some quantum ordered phases seen in
ref.
19
at larger
R
b
.
Classical, mean-
fi
eld, and entangled bulk phases
Without the Rabi term
Ω
, one would obtain classical Rydberg
crystals without a disordered phase. Figure
2
cshowstheclassical
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
fi
nd a stable
1
7
-density
crystal within unit cell sizes of up to 10 × 10. All ordered quantum
phases in Fig.
2
a 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 stri-
ated 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 qua-
litatively mean-
fi
eld states of the form
Q
i
α
i
0
i
+
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1
j
α
i
j
2
q
1
i
.
Figure
2
dshowsthemean-
fi
eld phase diagram. The disordered phase
does not appear, as it emerges from defect hopping and cannot be
described without some entanglement
25
.Themean-
fi
eld phase dia-
gram contains features of both the classical and quantum phase dia-
grams. The striated quantum phase indeed appears as a mean-
fi
eld
state, con
fi
rmed by the match between the mean-
fi
eld and exact cor-
relation functions (Fig.
3
a). However, the nematic phase does not
appear, and in its place is the same
1
6
-density crystal stabilized in the
classical phase diagram. This shows that a treatment of entanglement
is required to describe the nematic phase.
10
10
10
10
large gap
10
10
10
10
0
10
energy per site
exact gs
mean field gs
mean field 3×4
correlation
lattice distance
0.0
0.25
0.5
0.75
123
4
5
exact gs
mean field gs
correlation
lattice distance
0.0
0.25
0.5
0.75
123
4
5
partition location
von Neumann entropy
1.5
1.0
0.5
0.0
020406080100
Schmidt eigenvalue
10
10
10
10
10
classical spectrum
exact gs coefficients
a
b
d
c
striated
mean field
and exact gs
exact gs
mean field 3×4
mean field gs
Fig. 3 | Mean-
fi
eld striated versus entangled nematic phase. a
Density-density
correlation functions of the mean-
fi
eld and exact striated ground state, both at
(
δ
,
R
b
) = (3.1, 1.5); these agree, con
fi
rming the mean-
fi
eld nature of the striated
phase. (b) Density-density correlation functions for the entangled nematic phase
ground state and two different mean-
fi
eld 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 coordi-
nate. 8-fold degeneracy in (
a
)isshownastworowsof4dots.Thenon-mean-
fi
eld
(entangled) character of the nematic phase is evident.
c
Structure of the nematic
state in terms of classical con
fi
gurations constructed via compositions of 3 indi-
vidual column states
a
ji
,
b
,
c
ji
. In the classical limit, there are 4 distinct sets of low-
energy con
fi
gurations, all characterized by the absence of adjacent columns in the
same state (e.g.,
aa
:::
ji
) and large degeneracies due to permutational symmetry
between
a
ji
,
b
,and
c
ji
. The lowest in energy is 6-fold degenerate, corresponding
to the 3-star state. However, in the quantum nematic state the con
fi
gurations that
are slightly higher in energy have much larger wavefunction coef
fi
cients. The most
relevant classical states in the wavefunction are those with the greatest number of
possible single full column hops (e.g.,
a
→
b
) without introducing unfavorable states
like
aa
:::
ji
, revealing the role of itinerancy in the nematic phase.
d
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.
Article
https://doi.org/10.103
8/s41467-023-41166-0
Nature Communications
| (2023) 14:5397
4
Nature of the bulk nematic phase
Figure
3
b shows the density correlation function of the nematic phase,
which does not display mean-
fi
eld character. To reveal the phase
structure, Fig.
3
c 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 con
fi
gurations
a
ji
,
b
,
c
ji
. Column-column interactions,
however, prevent adjacent columns from being in the same con
fi
g-
uration (with excitations separated by 2 units); thus, states such as
abcb
...
are allowed, but
accb
...
are not. Without long-range
interactions, these column constraints give rise to an exponential
classical degeneracy. Long-range interactions partially lift the classical
degeneracy, yielding the
abc
...
crystal (3-star phase) and its 6-fold
degenerate permutations. However, after including quantum
fl
uctua-
tions and entanglement through a 4th order perturbative treatment of
σ
x
(giving rise to defect itinerancy),
abab
...
and related con
fi
gura-
tion energies are lowered below those of the
abc
...
con
fi
gurations;
the
fl
uctuations stabilize non-classical crystal con
fi
gurations (see
Supplementary Note 1). Figure
3
c gives the weights of the con
fi
gura-
tions in the computed quantum ground-state, which are distributed
across the exponentially numerous non-classical
abab
...
,
abcbab
etc., con
fi
gurations, with the classical crystal
abc
...
con
fi
gurations
strongly disfavored. The bi-partite entanglement entropy and entan-
glement spectrum are further shown in Fig.
3
d. Although the
fl
uctua-
tions are presumably of
fi
nite range, the entanglement spectrum
carries 3 large Schmidt values across every cut along the DMRG snake
MPS, showing the state is entangled across the entire supercell, and
well approximated by an MPS of bond dimension 3. The entanglement
structure emerges from the combination of defect itinerancy and the
constraints on adjacent columns. Thus, it is clear that quantum
fl
uc-
tuations are much stronger in this phase than in any of the surrounding
ordered phases. Assuming the entanglement is ultimately short-
ranged (i.e., on scales beyond the supercell sizes we can treat here),
this phase can be identi
fi
ed as containing strong
fl
uctuations around a
non-classical crystal, stabilized by an order from disorder mechanism
37
(further discussion in Supplementary Note 1).
Finite phase diagram
Current experiments are limited to lattices with open boundary con-
ditions consisting of a few hundred atoms
10
,
18
. To investigate how this
modi
fi
es the bulk behavior, we computed the phase diagram of
selected
fi
nite lattices from size 9 × 9 to 16 × 16, using DMRG for the
smaller sizes and our PEPS methodology for the larger ones.
We
fi
rst focus (in Fig.
4
a) on understanding the fate of the ordered
phases on the 15 × 15 lattice along three slices:
δ
=2.7,4.0,and5.0
(16 × 16 lattice phases, as well as other lattice sizes, are discussed in
Supplementary Notes 4
–
5). Here, many
fi
nite 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-
fi
eld
phase remains due to its commensurate boundary-bulk con
fi
gura-
tions, while in the region of strongest interactions, the nematic phase is
destabilized. A new region of classical order, called here the square
phase (Fig.
4
c), emerges across much of the
R
b
=1.5
−
1.8 region 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
fl
uctuations 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 experi-
mental non-zero order parameter space to a square/striated phase (see
Fig.
5
a, note, the order parameter does not distinguish between
square/striated orders). However, our simulations (Fig.
5
c) 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 identi
fi
ed a large region of
star order (Fig.
5
b). This assignment is complicated by edge effects,
which mean that the order parameter used does not cleanly distin-
guish 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 con
fi
rmed using a different, more sensitive order parameter
defect states
2.4
2.2
2.0
1.8
1.6
1.4
1.2
3.0
3.5
4.0
4.5
5.0
16 × 16
15 × 14
manual edge
defects
square
b
c
exact gs
finite 15×14
0.0
0.1
123
4
5
0.2
0.3
lattice distance
correlation
a
e
d
Fig. 4 | Phase diagram of the 15 × 15
fi
nite system and
fi
nite lattice orders. a
The
phase diagram, where colors correspond to the same phase classi
fi
cations as Fig.
2
.
Triangles represent tentative classi
fi
cation of points showing inconsistent PEPS
convergence, see Supplementary Methods. A new order, which we call square, is
speci
fi
ed in (
c
) and various examples of boundary-bulk frustrated ground states in
(
b
).
d
The density pro
fi
le 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
fi
nite nematic phase to the converged 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-
fi
eld states (see Fig.
3
b).
Article
https://doi.org/10.103
8/s41467-023-41166-0
Nature Communications
| (2023) 14:5397
5