mygraygray0.6

# Quantum Phases of Matter on a 256-Atom Programmable Quantum Simulator

Sepehr Ebadi1, Tout T. Wang1, Harry Levine1, Alexander Keesling1, Giulia
Semeghini1, Ahmed Omran1,2, Dolev Bluvstein1, Rhine Samajdar1, Hannes
Pichler3,4, Wen Wei Ho1,5, Soonwon Choi6, Subir Sachdev1, Markus Greiner1,
Vladan Vuletić7, and Mikhail D. Lukin1 1Department of Physics, Harvard
University, Cambridge, MA 02138, USA
2QuEra Computing Inc., Boston, MA 02135, USA
3Institute for Theoretical Physics, University of Innsbruck, Innsbruck A-6020,
Austria
4Institute for Quantum Optics and Quantum Information, Austrian Academy of
Sciences, Innsbruck A-6020
5 Department of Physics, Stanford University, Stanford, CA 94305, USA.
6Department of Physics, University of California Berkeley, Berkeley, CA 94720,
USA
7Department of Physics and Research Laboratory of Electronics, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA.

###### Abstract

Motivated by far-reaching applications ranging from quantum simulations of
complex processes in physics and chemistry to quantum information processing
Preskill (2018), a broad effort is currently underway to build large-scale
programmable quantum systems. Such systems provide unique insights into
strongly correlated quantum matter Zhang  _et al._ (2017); Gross and Bloch
(2017); Choi  _et al._ (2016); Mazurenko  _et al._ (2017); Bernien  _et al._
(2017), while at the same time enabling new methods for computation Aleiner
_et al._ (2020); Wang  _et al._ (2020); Arute  _et al._ (2019); Zhong  _et
al._ (2020) and metrology Giovannetti  _et al._ (2011). Here, we demonstrate a
programmable quantum simulator based on deterministically prepared two-
dimensional arrays of neutral atoms, featuring strong interactions controlled
via coherent atomic excitation into Rydberg states Browaeys and Lahaye (2020).
Using this approach, we realize a quantum spin model with tunable interactions
for system sizes ranging from 64 to 256 qubits. We benchmark the system by
creating and characterizing high-fidelity antiferromagnetically ordered
states, and demonstrate the universal properties of an Ising quantum phase
transition in (2+1) dimensions Sachdev (2011). We then create and study
several new quantum phases that arise from the interplay between interactions
and coherent laser excitation Samajdar  _et al._ (2020a), experimentally map
the phase diagram, and investigate the role of quantum fluctuations. Offering
a new lens into the study of complex quantum matter, these observations pave
the way for investigations of exotic quantum phases, non-equilibrium
entanglement dynamics, and hardware-efficient realization of quantum
algorithms.

Recent breakthroughs have demonstrated the potential of programmable quantum
systems, with system sizes reaching around fifty trapped ions Zhang  _et al._
(2017); Pagano  _et al._ (2020); Friis  _et al._ (2018) or superconducting
qubits Aleiner  _et al._ (2020); Wang  _et al._ (2020); Arute  _et al._
(2019), for simulations and computation. Correlation measurements with over
seventy photons have been used to perform boson sampling Zhong  _et al._
(2020), while optical lattices with hundreds of atoms are being used to
explore Hubbard models Gross and Bloch (2017); Choi  _et al._ (2016);
Mazurenko  _et al._ (2017). Larger-scale Ising spin systems have been realized
using superconducting elements Harris  _et al._ (2018), but they lack the
coherence essential for probing quantum matter.

Neutral atom arrays have recently emerged as a promising platform for
realizing programmable quantum systems Browaeys and Lahaye (2020); Labuhn  _et
al._ (2016); Bernien  _et al._ (2017). Based on individually trapped and
detected cold atoms in optical tweezers with strong interactions between
Rydberg states Morgado and Whitlock (2020), atom arrays have been utilized to
explore quantum dynamics in one- and two-dimensional systems Bernien  _et al._
(2017); Lienhard  _et al._ (2018); Guardado-Sanchez  _et al._ (2018); Kim  _et
al._ (2018); de Léséleuc  _et al._ (2019); Keesling  _et al._ (2019), to
create high-fidelity Madjarov  _et al._ (2020) and large-scale Omran  _et al._
(2019) entanglement, to perform parallel quantum logic operations Graham  _et
al._ (2019); Levine  _et al._ (2019), and to realize optical atomic clocks
Madjarov  _et al._ (2019); Young  _et al._ (2020). While large numbers of
atoms have been trapped Young  _et al._ (2020) and rearranged in two and three
dimensions Schymik  _et al._ (2020); Kumar  _et al._ (2018); Ohl de Mello  _et
al._ (2019); Barredo  _et al._ (2018), coherent manipulation of programmable,
strongly interacting systems with more than a hundred individual particles
remains an outstanding challenge. Here, we realize a programmable quantum
simulator using arrays of up to 256 neutral atoms with tunable interactions,
demonstrating several novel quantum phases and quantitatively probing the
associated phase transitions.

Figure 1: Programmable two-dimensional arrays of strongly-interacting Rydberg
atoms. a. Atoms are loaded into a 2D array of optical tweezer traps and
rearranged into defect-free patterns by a second set of moving tweezers.
Lasers at 420 nm and 1013 nm drive a coherent two-photon transition in each
atom between ground state $|g\rangle=|5S_{1/2},F=2,m_{F}=-2\rangle$ and
Rydberg state $|r\rangle=|70S_{1/2},m_{j}=-1/2,m_{I}=-3/2\rangle$. b.
Fluorescence image of initial random loading of atoms, followed by
rearrangement to a defect-free 15$\times$15 (225 atoms) square array. After
this initialization, the atoms evolve coherently under laser excitation with
Rabi frequency $\Omega(t)$ and detuning $\Delta(t)$, and long-range
interactions $V_{ij}$. Finally, the state of each atom is read out, with atoms
excited to $|r\rangle$ detected as loss and marked with red circles. Shown on
the far right is an example measurement following quasi-adiabatic evolution
into the checkerboard phase. c, d. Similar evolution on honeycomb and
triangular lattices result in analogous ordered phases of Rydberg excitations
with filling $1/2$ and $1/3$, respectively.

##  I Programmable Rydberg Arrays in 2D

Our experiments are carried out on the second generation of an experimental
platform described previously Bernien  _et al._ (2017). The new apparatus uses
a spatial light modulator (SLM) to form a large, two-dimensional (2D) array of
optical tweezers in a vacuum cell (Fig. 1a, Methods). This static tweezer
array is loaded with individual 87Rb atoms from a magneto-optical trap (MOT),
with a uniform loading probability of 50–60% across up to 1000 tweezers. We
rearrange the initially loaded atoms into programmable, defect-free patterns
using a second set of moving optical tweezers that are steered by a pair of
crossed acousto-optical deflectors (AODs) to arbitrary positions in 2D (Fig.
1a) Barredo  _et al._ (2016). Our parallel rearrangement protocol (see
Methods) enables rearrangement into a wide variety of geometries including
square, honeycomb, and triangular lattices (left panels in Fig. 1b-d). The
procedure takes a total time of $50$–$100$ ms for arrays of up to a few
hundred atoms and results in filling fractions exceeding $99\%$.

Qubits are encoded in the electronic ground state $|g\rangle$ and the highly-
excited $n=70$ Rydberg state $|r\rangle$ of each atom. We illuminate the
entire array from opposite sides with two counter-propagating laser beams at
420 and 1013 nm, shaped into light sheets (see Methods), to coherently couple
$|g\rangle$ to $|r\rangle$ via a two-photon transition (Fig. 1a).

The resulting many-body dynamics $U(t)$ are governed by a combination of the
laser excitation and long-range van der Waals interactions between Rydberg
states ($V_{ij}=V_{0}/|\mathbf{x_{i}}-\mathbf{x_{j}}|^{6}$), described by the
Hamiltonian

$\frac{H}{\hbar}=\frac{\Omega}{2}\sum_{i}(|g_{i}\rangle\langle
r_{i}|+|r_{i}\rangle\langle
g_{i}|)-\Delta\sum_{i}n_{i}+\sum_{i<j}V_{ij}n_{i}n_{j}$ (1)

where $\hbar$ is the Planck’s constant, $n_{i}=|r_{i}\rangle\langle r_{i}|$,
and $\Omega$ and $\Delta$ are the two-photon Rabi frequency and detuning,
respectively. After evolution under the Hamiltonian (1), the state of each
atomic qubit is read out by fluorescence imaging that detects only atoms in
$|g\rangle$, while atoms in $|r\rangle$ are detected as loss. Detection
fidelities exceed 99% for both states (see Methods).

The Rydberg blockade mechanism Jaksch  _et al._ (2000); Lukin  _et al._ (2001)
is central to understanding the programmable dynamics driven by the
Hamiltonian (1). It originates from the long-range interactions between
Rydberg states, providing an effective constraint that prevents simultaneous
excitation of atoms within a blockade radius
$R_{b}\equiv(V_{0}/\Omega)^{1/6}$. We control the effective blockade range
$R_{b}/a$ by programming the lattice spacing $a$ for the atom array. Using
these control tools, we explore quantum evolution resulting in a wide variety
of quantum phases.

##  II Checkerboard Phase

Figure 2: Benchmarking of quantum simulator using checkerboard ordering. a. A
quasi-adiabatic detuning sweep $\Delta(t)$ at constant Rabi frequency $\Omega$
is used to prepare the checkerboard ground state with high fidelity. b. Two-
site correlation function $G^{(2)}(k,l)$, averaged over all pairs of atoms on
a 12$\times$12 array, showing near-perfect alternating correlations throughout
the entire system. c. Exponential fits of rectified horizontal and vertical
correlations are used to extract correlation lengths in the corresponding
directions $\xi_{H}$ and $\xi_{V}$. d. Histogram of many-body state occurrence
frequency after 6767 repetitions of the experiment on a 12$\times$12 array.
The two most frequently occurring microstates correspond to the two perfect
checkerboard orderings, and the next four most common ones are those with a
single defect in one of the corners of the array. e. Probability of finding a
perfect checkerboard ground state as a function of array size. The slightly
higher probabilities in odd$\times$odd systems is due to commensurate edges on
opposing sides of the array. All data in this figure are conditioned on
defect-free rearrangement of the array.

The smallest value of $R_{b}/a$ that results in an ordered phase for the
quantum many-body ground state of the system corresponds to $R_{b}/a\approx
1$, where only one out of every pair of nearest-neighbor atoms can be excited
to $|r\rangle$. On a square array, this constraint leads to a
$\mathbb{Z}_{2}$-symmetry-broken checkerboard phase with an antiferromagnetic
(AF) ground state. To realize such a state, we initialize the array at
$R_{b}/a=1.15$ ($a=6.7~{}\mu$m, $\Omega=2\pi\times 4.3$ MHz) with all atoms in
$|g\rangle$. We then dynamically sweep the detuning $\Delta$ from negative to
positive values while keeping the Rabi frequency $\Omega$ fixed to bring the
system quasi-adiabatically into the checkerboard phase (Fig. 1b and Fig. 2a).
A similar approach can be used to create analogous ordered phases on other
lattice geometries (Fig. 1c, d).

We quantify the strength of antiferromagnetic correlations in the checkerboard
phase over many experimental repetitions using the connected density-density
correlator $G^{(2)}(k,l)=\frac{1}{N_{(k,l)}}\sum_{i,j}(\langle
n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle n_{j}\rangle)$, where the sum is
over all pairs of atoms $(i,j)$ separated by the same relative lattice
displacement $\mathbf{x}$ $=$ $(k,l)$ sites, normalized by the number of such
pairs $N_{(k,l)}$. Our measurement of $G^{(2)}(k,l)$ on a 12$\times$12 system
(Fig. 2b) yields horizontal and vertical correlation lengths of $\xi_{H}=$
11.1(1) and $\xi_{V}=$ 11.3(1) respectively (Fig. 2c), showing long-range
correlations across the entire 144 atom array. These exceed the values
reported previously for two-dimensional systems Lienhard  _et al._ (2018);
Guardado-Sanchez  _et al._ (2018) by nearly an order of magnitude.

Single-site readout also allows us to study individual many-body states of our
system (Fig. 2d). Out of 6767 repetitions on a 12x12 array, the two perfectly
ordered states $|\text{AF}_{1}\rangle$ and $|\text{AF}_{2}\rangle$ are by far
the most frequently observed microstates, with near-equal probabilities
between the two. We benchmark our state preparation by measuring the
probability of observing perfect checkerboard ordering as a function of system
size (Fig. 2e). We find empirically that the probability scales with the
number of atoms according to an exponential $0.97^{N}$, offering a benchmark
that includes all experimental imperfections such as finite detection
fidelity, non-adiabatic state preparation, spontaneous emission, and residual
quantum fluctuations in the ordered state (see Methods). Remarkably, even for
a system size as large as 15$\times$15 (225 atoms), we still observe the
perfect antiferromagnetic ground state with probability 0.10${{}^{+5}_{-4}}$%
within the exponentially large Hilbert space of dimension $2^{225}\approx
10^{68}$.

##  III (2+1)D Ising Quantum Phase Transition

We now describe quantitative studies of the quantum phase transition into the
checkerboard phase. Quantum phase transitions fall into universality classes
characterized by critical exponents that determine universal behavior near the
quantum critical point, independent of the microscopic details of the
Hamiltonian Sachdev (2011). The transition into the checkerboard phase is
expected to be in the paradigmatic—but never previously observed—quantum Ising
universality class in (2+1) dimensions Samajdar  _et al._ (2020a) (with
expected dynamical critical exponent $z=1$ and correlation length critical
exponent $\nu=0.629$).

To explore universal scaling across this phase transition for a large system,
we study the dynamical build-up of correlations associated with the quantum
Kibble-Zurek mechanism Zurek  _et al._ (2005); Keesling  _et al._ (2019) on a
$16\times 16$ (256 atoms) array, at fixed $R_{b}/a=1.15$. We start at a large
negative detuning with all atoms in $|g\rangle$ and linearly increase
$\Delta/\Omega$, stopping at various points to measure the growth of
correlations across the phase transition (Fig. 3a, b). Slower sweep rates
$s=d\Delta/dt$ result in longer correlation lengths $\xi$, as expected (Fig.
3c).

The quantum Kibble-Zurek mechanism predicts a universal scaling relationship
between the control parameter $\Delta$ and the correlation length $\xi$.
Specifically, when both $\Delta$ and $\xi$ are rescaled with the sweep rate
$s$ (relative to a reference rate $s_{0}$)

$\displaystyle\tilde{\xi}=\xi(s/s_{0})^{\mu}$ (2)
$\displaystyle\tilde{\Delta}=(\Delta-\Delta_{c})(s/s_{0})^{\kappa}$ (3)

with exponents $\mu\equiv\nu/(1+z\nu)$ and $\kappa\equiv-1/(1+z\nu)$, then
universality implies that the rescaled $\tilde{\xi}$ vs. $\tilde{\Delta}$
collapses onto a single curve Keesling  _et al._ (2019) for any sweep rate
$s$. Taking $z=1$ to be fixed (as expected for a Lorentz-invariant theory), we
extract $\nu$ for our system by finding the value that optimizes this
universal collapse.

In order to obtain $\nu$, we first independently determine the position of the
critical point $\Delta_{c}$, which corresponds to the peak of the
susceptibility $\chi=-\partial^{2}\langle H\rangle/\partial\Delta^{2}$ and is
associated with a vanishing gap Sachdev (2011). For adiabiatic evolution under
the Hamiltonian (1), the susceptibility $\chi$ is related to the mean Rydberg
excitation density $\langle n\rangle$ by $\chi=\partial\langle
n\rangle/\partial\Delta$ according to the Hellman-Feynman theorem. We measure
$\langle n\rangle$ vs. $\Delta$ along a slow linear sweep to remain as
adiabatic as possible. We take the numerical derivative of the fitted data to
obtain $\chi$, finding its peak to be at $\Delta_{c}/\Omega=1.12(4)$ (see
Methods).

Having identified the position of the critical point, we now extract the value
of $\nu$ that optimizes data collapse (inset of Fig. 3d and Methods). The
resulting $\nu=0.62(4)$ rescales the experimental data to clearly fall on a
single universal curve (Fig. 3d). This measurement is in good agreement with
the predicted $\nu=0.629$ for the quantum Ising universality class in (2+1)
dimensionsSamajdar  _et al._ (2020a), and distinct from both the mean-field
valueSachdev (2011) of $\nu=1/2$ and the previously verified value in (1+1)
dimensions Keesling  _et al._ (2019) of $\nu=1$. Despite imperfections
associated with non-adiabatic state preparation and decoherence in our system,
this demonstration of universal scaling highlights opportunities for
quantitative studies of quantum critical phenomena on our platform.

Figure 3: Observation of the (2+1)D Ising quantum phase transition on a
16$\times$16 array. a. The transition into the checkerboard phase is explored
using a linear detuning sweep $\Delta(t)$ at constant $\Omega$. The resulting
checkerboard ordering is measured at various endpoints. b. Example of growing
correlations $G^{(2)}$ with increasing $\Delta/\Omega$ along a linear sweep
with sweep rate $s=15$ MHz/$\mu$s. c. Growth of correlation length $\xi$ for
$s$ spanning an order of magnitude from 15 MHz/$\mu$s to 120 MHz/$\mu$s. $\xi$
used here measures correlations between the coarse-grained local staggered
magnetization (see Methods). d. For an optimized value of the critical
exponent $\nu$, all curves collapse onto a single universal curve when
rescaled relative to the quantum critical point $\Delta_{c}$. Inset: distance
$D$ between all pairs of rescaled curves as a function of $\nu$ (see Methods).
The minimum at $\nu=0.62(4)$ (red dashed line) yields the experimental value
for the critical exponent (red and gray shaded regions indicate
uncertainties).  Figure 4: Phase diagram of the two-dimensional square
lattice. a. Example fluorescence image of atoms in the checkerboard phase and
the corresponding Fourier transform averaged over many experimental
repetitions $\langle\mathcal{F}(\mathbf{k})\rangle$, highlighting the peak at
$(\pi,\pi)$ (circled). b. Image of atoms in the striated phase and the
corresponding $\langle\mathcal{F}(\mathbf{k})\rangle$ highlighting peaks at
$(0,\pi)$, $(\pi,0)$ and $(\pi,\pi)$ (circled). c. Image of atoms in the star
phase with corresponding Fourier peaks at $(\pi/2,\pi)$ and $(\pi,0)$
(circled), as well as at symmetric partners $(\pi,\pi/2)$ and $(\pi,0)$. d.
The experimental phase diagram is constructed by measuring order parameters
for each of the three phases for different values of the tunable blockade
range $R_{b}/a$ and detuning $\Delta/\Omega$. Red markers indicate the
numerically calculated phase boundaries (see Methods). e. The order parameters
evaluated numerically using DMRG for a 9$\times$9 array (see Methods).

##  IV Phase Diagram of the Square Lattice

A rich variety of new phases have been recently predicted for the square
lattice when Rydberg blockade is extended beyond nearest neighbors Samajdar
_et al._ (2020a). To map this phase diagram experimentally, we use the Fourier
transform of single-shot measurement outcomes
$\mathcal{F}(\mathbf{k})=\left|\sum_{i}\text{exp}(i\mathbf{k\cdot
x_{i}}/a)n_{i}/\sqrt{N}\right|$, which characterizes long-range order in our
system. For instance, the checkerboard phase shows a prominent peak at
$\mathbf{k}=(\pi,\pi)$, corresponding to the canonical antiferromagnetic order
parameter: the staggered magnetization (Fig. 4a). We construct order
parameters for all observed phases using the symmetrized Fourier transform
$\tilde{\mathcal{F}}(k_{1},k_{2})=\langle\mathcal{F}(k_{1},k_{2})+\mathcal{F}(k_{2},k_{1})\rangle/2$,
averaged over experimental repetitions, which takes into account the
reflection symmetry in our system (see Methods).

When interaction strengths are increased such that next-nearest (diagonal)
neighbor excitations are suppressed by Rydberg interactions
($R_{b}/a\gtrsim\sqrt{2}$), translational symmetry along the diagonal
directions is also broken, leading to the appearance of a new striated phase
(Fig. 4b). In this phase, Rydberg excitations are mostly located two sites
apart and hence appear both on alternating rows and alternating columns. This
ordering is immediately apparent through the observation of prominent peaks at
$\mathbf{k}~{}=~{}(0,\pi)$, $(\pi,0)$, and $(\pi,\pi)$ in the Fourier domain.
As discussed and demonstrated below, quantum fluctuations, appearing as
defects on single shot images (Fig. 4b), play a key role in stabilizing this
phase.

At even larger values of $R_{b}/a\gtrsim 1.7$, the star phase emerges, with
Rydberg excitations placed every four sites along one direction and every two
sites in the perpendicular direction. There are two possible orientations for
the ordering of this phase, so Fourier peaks are observed at $\mathbf{k}$ =
$(\pi,0)$ and $(\pi/2,\pi)$, as well as at their symmetric partners $(0,\pi)$
and $(\pi,\pi/2)$ (Fig. 4c). In the thermodynamic limit, the star ordering
corresponds to the lowest-energy classical configuration of Rydberg
excitations on a square array with a density of 1/4.

We now systematically explore the phase diagram on 13$\times$13 (169 atoms)
arrays, with dimensions chosen to be simultaneously commensurate with
checkerboard, striated, and star orderings (see Methods). For each value of
the blockade range $R_{b}/a$, we linearly sweep $\Delta$ (similar to Fig. 3a
but with a ramp-down time of 200 ns), stopping at evenly spaced endpoints to
raster the full phase diagram. For every endpoint, we extract the order
parameter corresponding to each many-body phase, and plot them separately to
show their prominence in different regions of the phase diagram (Fig. 4d).

We compare our observations with numerical simulations of the phase diagram
using the density-matrix renormalization group (DMRG) on a smaller 9$\times$9
array with open boundary conditions (Fig. 4e and red markers in Fig. 4d). We
find excellent agreement in the extent of the checkerboard phase. For the
striated and star phases, we also find good similarity between experiment and
theory, although due to their larger unit cells and the existence of many
degenerate configurations, these two phases are more sensitive to both edge
effects and experimental imperfections. We emphasize that the numerical
simulations evaluate the order parameter for the exact ground state of the
system at each point, while the experiment quasi-adiabatically prepares these
states via a dynamical process. These results establish the potential of
programmable quantum simulators with tunable, long-range interactions for
studying large quantum many-body systems that are challenging to access with
state-of-the-art computational tools Felser  _et al._ (2020).

##  V Quantum Fluctuations in the Striated Phase

Figure 5: Probing correlations and coherence in the striated phase via quench
dynamics. a. Unit cell of striated ordering (dashed box) with (0,0) and (1,1)
sublattices outlined in red and blue, respectively. The fill shade on each
site reflects the mean Rydberg excitation. b. The variational states for the
(0,0) and (1,1) sublattices are illustrated on the Bloch sphere (see Methods).
The black arrow illustrates the phase $\phi_{q}$ of $\Omega$ during the
quench. c. Probability $P^{(d)}$ of an excitation, conditioned on observing no
nearest-neighbor excitations, and zero (red), three (light blue), or four
(dark blue) diagonal next-nearest neighbor excitations. $P^{(0)}$ is plotted
for $\phi_{q}=\pi/2$, showing resonant de-excitation of the (0,0) sublattice
near the bare-atom resonance (leftmost vertical line). $P^{(3)}$ and $P^{(4)}$
are plotted for $\phi_{q}=-\pi/2$, showing excitation peaks for the (1,1)
sublattice at interaction shifts corresponding to 3 or 4 diagonal neighbors
(two rightmost vertical lines). d, e. $P^{(0)}$ and $P^{(4)}$ vary with quench
phase $\phi_{q}$ at their corresponding resonances ($\Delta_{q}/2\pi$ = 1.4
and 20.4 MHz, respectively), demonstrating coherence on both the (0,0) and
(1,1) sublattices. Solid line fits are used to extract Bloch vector
components.

We now explore the nature of the striated phase. In contrast to the
checkerboard and star phases, which can be understood from a dense-packing
argument Samajdar  _et al._ (2020a), this phase has no counterpart in the
classical limit ($\Omega\to 0$) (see Methods). Striated ordering allows the
atoms to lower their energy by partially aligning with the transverse field,
favoring this phase at finite $\Omega$. This can be seen by considering the
$2\times 2$ unit cell, within which one site has a large Rydberg excitation
probability (designated the (0,0) sublattice) (Fig. 5a). Excitations on its
nearest-neighbor (0,1) and (1,0) sublattices are suppressed due to strong
Rydberg blockade. The remaining atoms on the (1,1) sublattice have no nearest
neighbors in the Rydberg state and experience a much weaker interaction from
four next-nearest (diagonal) neighbors on the (0,0) sublattice, thus allowing
the (1,1) atoms to lower their energy by forming a coherent superposition
between ground and Rydberg states (Fig. 5b).

We experimentally study quantum fluctuations in this phase by observing the
response of the system to short quenches (with quench times
$t_{q}<1/\Omega_{q}$). The dependence on the detuning $\Delta_{q}$ and laser
phase $\phi_{q}$ of the quench contains information about local correlations
and coherence, which allows us to characterize the quantum states on the
different sublattices. The quench resonance for each site depends on the state
of its nearest and next-nearest neighbors. Due to the large difference between
the interaction energies on the (0,0) and (1,1) sublattices, when one
sublattice is resonantly driven, the other is effectively frozen.

The nature of the striated phase is revealed using nine-particle operators to
measure the state of an atom, conditioned on its local environment.
Specifically, we evaluate the conditional Rydberg density $P^{(d)}$, defined
as the excitation probability of an atom if all nearest neighbors are in
$|g\rangle$, and exactly $d$ next-nearest (diagonal) neighbors are in
$|r\rangle$ (see Methods). For $d=0$, we observe a dip in $P^{(0)}$ near the
bare atom resonance (Fig. 5c), corresponding to resonant de-excitation of the
(0,0) sublattice. Meanwhile, $P^{(3)}$ and $P^{(4)}$ have two separate peaks
that correspond to resonant excitation of the (1,1) sublattice with $d=3$ and
$d=4$ next-nearest neighbor excitations, respectively (Fig. 5c). Remarkably,
we find that the quench response of both the (0,0) and (1,1) sublattices
depends on the phase $\phi_{q}$ of the driving field during the quench (Fig.
5d,e). The measured visibilities, together with a simple mean-field model (see
Methods), enable the estimation of unknown Bloch vector components on the two
sublattices, yielding $\langle\sigma_{x}\rangle=-0.82(6)$,
$\langle\sigma_{y}\rangle=0.25(2)$ for the (0,0) sublattice, and
$\langle\sigma_{x}\rangle=-0.45(4)$, $\langle\sigma_{y}\rangle=0.09(1)$ for
the (1,1) sublattice. We emphasize that accurate characterization requires the
use of more sophisticated variational wavefunctions (based on e.g. tensor
networks) and warrants further investigation. This approach can also be
extended through techniques such as shadow tomography Huang  _et al._ (2020).

##  VI Outlook

These experiments demonstrate that two-dimensional Rydberg atom arrays
constitute a powerful platform for programmable quantum simulations with
hundreds of qubits. We expect that system size, quantum control fidelity, and
degree of programmability can all be increased considerably via technical
improvements. In particular, array sizes and rearrangement fidelities, along
with atomic state readout, are currently limited by collisions with background
gas particles, and can be improved with an upgraded vacuum system Madjarov
_et al._ (2020) and increased photon collection efficiency. Quantum coherence
can be enhanced using higher-power Rydberg lasers and by encoding qubits in
hyperfine ground states Levine  _et al._ (2019); Morgado and Whitlock (2020).
Tweezers with different atomic Saskin  _et al._ (2019); Madjarov  _et al._
(2020); Young  _et al._ (2020) and molecular Anderegg  _et al._ (2019); Liu
_et al._ (2018) species can provide additional features and lead to novel
applications in both quantum simulations and metrology. Finally, rapidly
switchable local control beams can be used to perform universal qubit
operations in parallel across the system.

Our experiments realize several new quantum phases and provide unprecedented
insights into quantum phase transitions in two-dimensional systems. These
studies can be extended along several directions, including the exploration of
non-equilibrium entanglement dynamics via rapid quenches across quantum phase
transitions Turner  _et al._ (2018); Surace  _et al._ (2020); Bluvstein  _et
al._ (2020), the investigation of topological quantum states of matter on
frustrated lattices Samajdar  _et al._ (2020b); Verresen  _et al._ (2020), the
simulation of lattice gauge theories Bañuls  _et al._ (2020); Notarnicola  _et
al._ (2020), and the study of broader classes of spin models using hyperfine
encoding Weimer  _et al._ (2010). Quantum information processing can also be
explored with hardware-efficient methods for multi-qubit operations Levine
_et al._ (2019) and protocols for quantum error correction and fault tolerant
control Auger  _et al._ (2017). Finally, our approach is well suited for
efficient implementation of novel algorithms for quantum optimization Farhi
_et al._ (2014); Zhou  _et al._ (2020) and sampling Wild  _et al._ (2020),
enabling experimental tests of their performance with system sizes exceeding
several hundred qubits.

##  VII Methods

2D Optical Tweezer Array
Our 2D tweezer array is generated by a free-running 810-nm Ti:Sapphire laser
(M Squared, 18-W pump). The laser illuminates a phase-control spatial light
modulator (Hamamatsu X13138-02), which imprints a computer generated hologram
on the wavefront of the laser field. The phase hologram is calculated using
the phase-fixed weighted Gerchberg-Saxton (WGS) algorithm Kim  _et al._ (2019)
to produce an arbitrary arrangement of tweezer spots after propagating to the
focus of a microscope objective (Mitutoyo: 3.5 mm glass thickness corrected,
50$\times$, NA=0.5). Using this method, we can create tweezer arrays with
roughly 1000 individual tweezers (Extended Data Fig. 1). When calculating the
phase hologram, we improve trap homogeneity by pre-compensating for the
variation in diffraction efficiency across the tweezer array (roughly given by
$\operatorname{sinc}^{2}(\frac{\pi}{2}(\theta_{\text{trap}}/\theta_{\text{max}}))$
where $\theta$ denotes the deflection angle from zeroth order).

We also use the phase control of our SLM to correct for optical aberrations on
tweezers within the experimentally-used field of view at the plane of atoms
(Extended Data Fig. 2). Aberrations reduce the peak intensity of focal spots
(characterized by the Strehl ratio), and correspondingly reduce the light
shift of our tweezers on the atoms. By measuring these light shifts as we scan
several low-order Zernike polynomials, we quantify and correct for various
aberrations in our optical system. Using this method, we compensate for 70
milliwaves of aberrations, observe a total increase of $18\%$ in our trap
intensity (Extended Data Fig. 2c), and measure a corresponding reduction in
the range of trap frequencies (Extended Data Fig. 2d). Aberration correction
additionally allows us to place tweezers closer together (minimum separation 3
$\mu$m) to reach larger blockade ranges $R_{b}/a$.

Tweezers in the array have waists $\sim$ 900 nm, trap depths of
$\sim~{}2\pi\times 17$ MHz, and radial trap frequencies of
$\sim~{}2\pi\times~{}80$ kHz. In each experimental cycle, the tweezers are
loaded from a magneto-optical trap (MOT) with uniform loading probabilities of
50–60% after 50–100 ms loading time.

Atom Rearrangement
Atoms are rearranged using an additional set of dynamically moving tweezers,
which are overlaid on top of the SLM tweezer array. These movable tweezers are
generated by a separate 809-nm laser source (DBR from Photodigm and tapered
amplifier from MOGLabs), and are steered with a pair of independently-
controlled crossed acousto-optic deflectors (AODs) (AA Opto Electronic
DTSX-400). Both AODs are driven by an arbitrary waveform which is generated in
real time using our home-built waveform generation software and an arbitrary
waveform generator (AWG) (M4i.6631-x8 by Spectrum Instrumentation).
Dynamically changing the RF frequency allows for continuous steering of beam
positions, and multi-frequency waveforms allow for multiple moving tweezers to
be created in parallel Endres  _et al._ (2016).

While many 2D sorting protocols have been described previously Barredo  _et
al._ (2016); Lee  _et al._ (2017); Sheng  _et al._ (2020); Ohl de Mello  _et
al._ (2019); Schymik  _et al._ (2020), we implement a novel protocol which is
designed to leverage parallel movement of multiple atoms simultaneously. More
specifically, we create a row of moving traps which scans upwards along the
SLM tweezer array to move one atom within each column up in parallel. This is
accomplished by scanning a single frequency component on the vertical AOD to
move from the bottom to the top of the SLM array, during which individual
frequency components are turned on and off within the horizontal AOD to create
and remove tweezers at the corresponding columns. This protocol is designed
for SLM tweezer arrays in which traps are grouped into columns and rows. While
this does constrain the possible geometries, most lattice geometries of
interest can still be defined on a subset of points along fixed columns and
rows.

Rearrangement Algorithm
Here we detail the rearrangement algorithm, which is illustrated in Extended
Data Fig. 3. It operates on an underlying rectangular grid of rows and
columns, where the SLM traps correspond to vertices of the grid. We pre-
program a set of ‘target traps’ that we aim to fill.

_Pre-sorting:_ We begin by ensuring that each column contains a sufficient
number of atoms to fill the target traps in that column. In each experimental
cycle, due to the random loading throughout the array, some columns may
contain excess atoms while other columns may lack a sufficient number of
atoms. Accordingly, we apply a ‘pre-sorting’ procedure in which we move atoms
between columns. To fill a deficient column $j$, we take atoms from whichever
side of $j$ has a larger surplus. We identify which atoms to take by finding
the nearest atoms from the surplus side which are in rows for which column $j$
has an empty trap. We then perform parallel horizontal sorting to move these
atoms into the empty traps of $j$ (not all surplus atoms need to be from the
same source column).

If the one-side surplus is insufficient to fill column $j$, then we move as
many surplus atoms as possible from this one side and leave $j$ deficient. We
then proceed to the next deficient column, and cycle through until all columns
have sufficient atoms. In typical randomly loaded arrays, this process takes a
small number of atom moves compared to the total number of moves needed for
sorting. This specific algorithm can fail to properly distribute atoms between
columns due to lack of available atoms, but these failures are rare and do not
limit the experimental capabilities.

_Ejection:_ After pre-sorting, we eject excess atoms in parallel by scanning
the vertical AOD frequency downward, beginning at a row in which we want to
pick up an atom, and ending below the bottom row of the array. In each
downward scan, we eject a single atom from each column containing excess
atoms; we repeat this process until all excess atoms are ejected.

_Parallel sorting within columns:_ After pre-sorting and ejection, each column
has the correct number of atoms to fill all of its target traps by moving
atoms up/down within the column. We now proceed to shuffle the
$i^{\text{th}}$-highest loaded atoms to the $i^{\text{th}}$-highest target
traps. As the atoms cannot move through each other, in a single vertical scan
atoms are moved as close as possible to their target locations, reaching their
targets unless they are blocked by another atom. We repeat upward/downward
scans until all atoms reach their target locations.

Rearrangement Parameters and Results
When using moving tweezers to pick up and drop off atoms in the SLM traps, the
moving tweezers ramp on/off over $15~{}\mu$s while positioned to overlap with
the corresponding SLM trap. The moving tweezers are approximately twice as
deep as the static traps, and move atoms between SLM traps with a speed of 75
$\mu$m/ms. Typical rearrangement protocols take a total of 50-100 ms to
implement in practice, depending on the size of the target array and the
random initial loading. Alignment of the AOD traps onto the SLM array is pre-
calibrated by measuring both trap arrays on a monitor CMOS camera and tuning
the AOD frequencies to match positions with traps from the SLM array.

A single round of rearrangement results in typical filling fractions of
$\sim~{}98.5\%$ across all target traps in the system. This is limited
primarily by the finite vacuum-limited lifetime ($\sim$ 10 s) and the duration
of the rearrangment procedure. To increase filling fractions, we perform a
second round of rearrangement (having skipped ejection in the first round to
keep excess atoms for the second round). Since the second round of
rearrangement only needs to correct for a small number of defects, it requires
far fewer moves and can be performed more quickly, resulting in less
background loss. With this approach, we achieve filling fractions of $\sim
99.2\%$ over more than 200 sites, with a total experimental cycle time of
$400$ ms.

Rydberg Laser System
Our Rydberg laser system is an upgraded version of a previous setup Omran  _et
al._ (2019). The 420-nm laser is a frequency-doubled Ti:Sapphire laser (M
Squared, 15-W pump). We stabilize the laser frequency by locking the
fundamental to an upgraded ultra-low expansion (ULE) reference cavity (notched
cylinder design from Stable Laser Systems), with finesse $\mathcal{F}=30,000$
at 840 nm. The 1013-nm laser source is an external-cavity diode laser (Toptica
DL Pro), which is locked to the same reference cavity ($\mathcal{F}=50,000$ at
1013 nm). To suppress high-frequency phase noise from this diode laser, we use
the transmitted light through the cavity, which is filtered by the narrow
cavity transmission spectrum ($30$ kHz linewidth) Levine  _et al._ (2018).
This filtered light is used to injection-lock another laser diode, whose
output is subsequently amplified to 10 W by a fiber amplifier (Azur Light
Systems).

Using beam shaping optics to homogeneously illuminate the atom array with both
Rydberg lasers, we achieve single-photon Rabi frequencies of
$(\Omega_{\text{420}},\Omega_{\text{1013}})=2\pi\times(160,50)~{}$MHz. We
operate with an intermediate state detuning $\delta=2\pi\times 1$ GHz,
resulting in two-photon Rabi frequency
$\Omega=\Omega_{\text{420}}\Omega_{\text{1013}}/2\delta\sim 2\pi\times 4$ MHz.
Small inhomogeneities in the Rydberg beams result in Rabi frequency variations
of $\sim 2\%$ RMS and $\sim 6\%$ peak-to-peak across the array. With these
conditions, we estimate an off-resonant scattering rate of $1/(20~{}\mu$s) for
atoms in $\ket{g}$ and $1/(150~{}\mu$s) for atoms in $\ket{r}$ at peak power.

Rydberg Beam Shaping
We illuminate our 2D atom array with counter-propagating Rydberg laser beams
from each side. Instead of using elliptical Gaussian beams, we shape both
Rydberg excitation beams into one-dimensional top-hats (light sheets) to
homogeneously illuminate the plane of atoms (Extended Data Fig. 4). To ensure
homogenous illumination over the entire array, we define our target field
profile in the plane of the atoms with both uniform amplitude cross section
and flat phase profile. Using a single phase-only SLM in the Fourier plane to
control both phase and amplitude in the image plane is inherently limited in
efficiency; therefore, in practice, we compromise between optimizing hologram
efficiency and beam homogeneity. We generate these holograms using the
conjugate gradient minimization algorithm (Extended Data Fig. 4c)Bowman  _et
al._ (2017). In all experiments in this work, we use 1D top-hat beams with a
flat-width of $105~{}\mu$m and a perpendicular Gaussian width of $25~{}\mu$m.
The conversion efficiencies into the top-hat modes are 30% for 420 nm and 38%
for 1013 nm.

Since holographic beam shaping relies on the intricate interplay of different
high spatial frequency components in the light field, it is extremely
sensitive to optical aberrations. We correct for all aberrations up to the
window of our vacuum chamber by measuring the amplitude and phase of the
electric field as it propagates through the optical beampath (Extended Data
Fig. 4a,b) Zupancic  _et al._ (2016). We do so by picking off a small portion
of the Rydberg beam and observing it on a camera with small pixel size and
with sensor cover removed for high-fidelity beam characterization (Imaging
Source DMM 27UP031-ML). In this way, we reduce the wavefront error in our beam
down to $\lambda/100$, and also use the measured field profile as the starting
guess in our hologram generation algorithm (Extended Data Fig. 4 a,b).
Furthermore, by imaging the top-hat beams we also correct for remaining
inhomogeneities by updating the input of our optimization algorithm (Extended
Data Fig. 4e,f). Due to aberrations and imperfections of the vacuum windows,
we observe slightly larger intensity variations on the atoms than expected
($\sim 3\%$ RMS, $\sim 10\%$ peak-to-peak).

Rydberg Pulses
After initializing our atoms in the ground state $|g\rangle$, the tweezer
traps are turned off for a short time ($<$5 $\mu$s) during which we apply a
Rydberg pulse. The pulse consists of a time-dependent Rabi frequency
$\Omega(t)$, time-dependent detuning $\Delta(t)$, and a relative instantaneous
phase $\phi(t)$. This is implemented by controlling the amplitude, frequency,
and phase of the 420-nm laser using a tandem AOM system, similar to what is
described previously Omran  _et al._ (2019).

_Quasi-Adiabatic Sweeps:_ To prepare many-body ground states with high
fidelity, we use an optimized quasi-adiabatic pulse shape (Fig. 2a). The
coupling $\Omega(t)$ is initially ramped on linearly at large fixed negative
detuning, held constant during the detuning sweep $\Delta(t)$, and finally
ramped down linearly at large fixed positive detuning. The detuning sweep
$\Delta(t)$ consists of a cubic spline interpolation between five points:
initial detuning, final detuning, an inflection point where the slope reaches
a minimum, and two additional points that define the duration of the slow part
of the sweep. The sweep used for finding perfect checkerboard ground state
probabilities (Fig. 2e) was obtained by optimizing the parameters of the
spline cubic sweep to maximize the correlation length on a 12$\times$12 (144
atoms) array. The sweep used in detection of the star and striated phases was
optimized based on maximizing their respective order parameters. In
particular, the inflection point was chosen to be near the position of the
minimum gap in these sweeps in order to maximize adiabaticity.

_Linear Sweeps:_ To probe the phase transition into the checkerboard phase
(Fig. 3), we use variable-endpoint linear detuning sweeps in which $\Omega$ is
abruptly turned off after reaching the endpoint. This ensures that projective
readout happens immediately after the end of the linear sweep instead of
allowing time for further dynamics, and is essential for keeping the system
within the quantum Kibble-Zurek regime. Linear sweeps are done from $\Delta$ =
$-16$ to 14 MHz ($\Delta/\Omega$ = -3.7 to 3.3) at sweep rates $s$ = 15, 21,
30, 42, 60, 85, and 120 MHz/$\mu$s. Data for locating the quantum critical
point (Extended Data Fig. 7a) is taken from the slowest of these sweeps ($s$ =
15 MHz/$\mu$s) to remain as close as possible to the ground state. For mapping
out the 2D phase diagram (Fig. 4), we use the same variable-endpoint linear
sweeps at fixed sweep rate $s=12~{}$MHz / $\mu$s, except that $\Omega$ is
ramped down over $200~{}$ns after reaching the endpoint.

State Detection
At the end of the Rydberg pulse, we detect the state of atoms by whether or
not they are recaptured in our optical tweezers. Atoms in $\ket{g}$ are
recaptured and detected with fidelity $99\%$, limited by the finite
temperature of the atoms and collisions with background gas particles in the
vacuum chamber.

Atoms excited to the Rydberg state are detected as a loss signal due to the
repulsive potential of the optical tweezers on $|r\rangle$. However, the
finite Rydberg state lifetimeBeterov  _et al._ (2009) ($\sim 80~{}\mu$s for
70S1/2) leads to a probability of $\sim 15\%$ for $|r\rangle$ atoms to decay
to $|g\rangle$ and be recaptured by the optical tweezers. In our previous work
Omran  _et al._ (2019), we increased tweezer trap depths immediately following
the Rydberg pulse to enhance the loss signal for atoms in $\ket{r}$. In 2D,
this approach is less effective because atoms which drift away from their
initial traps can still be recaptured in a large 3D trapping structure created
by out-of-plane interference of tweezers.

Following an approach similar to what has been previously demonstratedGraham
_et al._ (2019), we increase the Rydberg detection fidelity using a strong
microwave (MW) pulse to enhance the loss of atoms in $|r\rangle$ while leaving
atoms in $|g\rangle$ unaffected. The MW source (Stanford Research Systems
SG384) is frequency-tripled to $6.9$ GHz and amplified to 3 W (Minicircuits,
ZVE-3W-183+). The MW pulse, containing both $6.9$ GHz and harmonics, is
applied on the atoms using a microwave horn for $100$ ns. When applying a
Rydberg $\pi$-pulse immediately followed by the MW pulse, we observe loss
probabilities of $98.6(4)\%$. Since this measurement includes both error in
the $\pi$-pulse as well as detection errors, we apply a second Rydberg
$\pi$-pulse after the MW pulse, which transfers most of the remaining ground
state population into the Rydberg state. In this configuration, we observe
$99.1(4)\%$ loss probability, which is our best estimate for our Rydberg
detection fidelity (Extended Data Fig 5). We find that the loss signal is
enhanced by the presence of both MW fundamental and harmonic frequencies.

Coarse-Grained Local Staggered Magnetization
We define the coarse-grained local staggered magnetization for a site $i$ with
column and row indices $a$ and $b$, respectively, as:

$m_{i}=\frac{(-1)^{a+b}}{N_{i}}\sum_{\langle j,i\rangle}(n_{i}-n_{j})$

where $j$ is summed over nearest neighbors of site $i$ and $N_{i}$ is the
number of such nearest neighbors (4 in the bulk, 3 along the edges, or 2 on
the corners). The value of $m_{i}$ ranges from $-1$ to 1, with the extremal
values corresponding to the two possible perfect antiferromagnetic orderings
locally on site $i$ and its nearest neighbors (Extended Data Fig. 6a,b). The
two-site correlation function for $m_{i}$ can then be defined as an average
over experiment repetitions
$G^{(2)}_{m}(k,l)=\frac{1}{N_{(k,l)}}\sum_{i,j}(\langle
m_{i}m_{j}\rangle-\langle m_{i}\rangle\langle m_{j}\rangle)$, where the sum is
over all pairs of sites $i,j$ separated by a relative lattice distance of
$\mathbf{x}=(k,l)$ sites and normalized by the number of such pairs
$N_{(k,l)}$ (Extended Data Fig. 6c). We obtain the correlation length $\xi$ by
fitting an exponential decay to the radially averaged $G^{(2)}_{m}(k,l)$
(Extended Data Fig. 6d). The coarse-grained local staggered magnetization
$m_{i}$ is defined such that the corresponding $G^{(2)}_{m}(k,l)$ is isotropic
(Extended Data Fig. 6c), which makes for natural radial averaging. This radial
average captures correlations across the entire array better than purely
horizontal or vertical correlation lengths $\xi_{H}$ and $\xi_{V}$, which are
more sensitive to edge effects.

Determination of the Quantum Critical Point
To accurately determine the location of the quantum critical point
$\Delta_{c}$ for the transition into the checkerboard phase, we measure mean
Rydberg excitation $\langle n\rangle$ vs. detuning $\Delta$ for a slow linear
sweep with sweep rate $s=15$ MHz/$\mu$s (Extended Data Fig. 7a). To smoothen
the measured curve, we fit a polynomial for $\langle n\rangle$ vs. $\Delta$
and take its numerical derivative to identify the peak of the susceptibility
$\chi$ as the critical pointSachdev (2011) (Extended Data Fig. 7b).

Small oscillations in $\langle n\rangle$ result from the linear sweep not
being perfectly adiabatic. To minimize the effect of this on our fitting, we
use the lowest-degree polynomial (cubic) whose derivative has a peak, and
choose a fit window in which the reduced chi-squared metric indicates a good
fit. Several fit windows around $\Delta/\Omega=0$ to 2 give good cubic fits,
and we average results from each of these windows to obtain
$\Delta_{c}/\Omega$ = 1.12(4).

We also numerically extract the critical point for a system with numerically-
tractable dimensions of 10$\times$10\. Using the density-matrix
renormalization group (DMRG) algorithm, we evaluate $\langle n\rangle$ as a
function of detuning $\Delta$, and then take the derivative to obtain a peak
of the susceptibility at $\Delta_{c}/\Omega=1.18$ (Extended Data Fig. 7c,d).
To corroborate the validity of our experimental fitting procedure, we also fit
cubic polynomials to the DMRG data and find that the extracted critical point
is close to the exact numerical value (Extended Data Fig. 7d). This numerical
estimate of the critical point for a 10$\times$10 array is consistent with the
experimental result on a larger $16\times 16$ array. Moreover, our experiments
on arrays of different sizes show that $\Delta_{c}/\Omega$ does not vary
significantly between $12\times 12$, $14\times 14$, and $16\times 16$ arrays
(Extended Data Fig. 8b).

Data Collapse for Universal Scaling
Optimizing the universal collapse of rescaled correlation length $\tilde{\xi}$
vs. rescaled detuning $\tilde{\Delta}$ requires defining a measure of the
distance between rescaled curves for different sweep rates $s_{i}$. Given
$\tilde{\xi}^{(i)}_{j}$ and $\tilde{\Delta}^{(i)}_{j}$, where the index $i$
corresponds to sweep rate $s_{i}$ and $j$ labels sequential data points along
a given curve, we define a distance Bhattacharjee and Seno (2001)

$D=\sqrt{\frac{1}{N}\sum_{i}\sum_{i^{\prime}\neq
i}\sum_{j}\left\lvert\tilde{\xi}^{(i^{\prime})}_{j}-f^{(i)}\left(\tilde{\Delta}^{(i^{\prime})}_{j}\right)\right\rvert^{2}}.$
(4)

The function $f^{(i)}(\tilde{\Delta})$ is the linear interpolation of
$\tilde{\xi}^{(i)}_{j}$ vs. $\tilde{\Delta}^{(i)}_{j}$, while $N$ is the total
number of terms in the three nested sums. The sum over $j$ only includes
points that fall within the domain of overlap of all data sets, avoiding the
problem of linear interpolation beyond the domain of any single data set.
Defined in this way, the collapse distance $D$ measures all possible
permutations of how far each rescaled correlation growth curve is from curves
corresponding to other sweep rates.

Applied to our experimental data, $D$ is a function of both the location of
the critical point $\Delta_{c}$ and the critical exponent $\nu$ (Extended Data
Fig. 8a). Using the independently measured $\Delta_{c}/\Omega=1.12(4)$, we
obtain $\nu=0.62(4)$ for optimal data collapse, and illustrate in particular
the better collapse for this value than for other values of $\nu$ (Extended
Data Fig. 8c-e). The quoted uncertainty is dominated by the corresponding
uncertainty of the extracted $\Delta_{c}/\Omega$, rather than by the precision
of finding the minimum of $D$ for a given $\Delta_{c}/\Omega$. Our experiments
give consistent values of $\Delta_{c}/\Omega$ and $\nu$ for systems of size
12$\times$12, 14$\times$14, and 16$\times$16 (Extended Data Fig. 8b).

Order Parameters for Many-Body Phases
We construct order parameters to identify each phase using the Fourier
transform to quantify the amplitude of the observed density-wave ordering. We
define the symmetrized Fourier transform
$\tilde{\mathcal{F}}(k_{1},k_{2})=\langle\mathcal{F}(k_{1},k_{2})+\mathcal{F}(k_{2},k_{1})\rangle/2$
to take into account the $\textit{C}_{4}$ rotation symmetry between possible
ground-state orderings for some phases. For the star phase, the Fourier
amplitude $\tilde{\mathcal{F}}(\pi,\pi/2)$ is a good order parameter because
ordering at $\mathbf{k}=(\pi,\pi/2)$ is unique to this phase. The striated
phase, on the other hand, shares its Fourier peaks at $\mathbf{k}$ = $(\pi,0)$
and $(0,\pi)$ with the star phase, and its peak at $\mathbf{k}$ = $(\pi,\pi)$
with the checkerboard phase; hence, none of these peaks alone can serve as an
order parameter. We therefore construct an order parameter for the striated
phase to be $\tilde{\mathcal{F}}(0,\pi)-\tilde{\mathcal{F}}(\pi/2,\pi)$, which
is nonzero in the striated phase and zero in both checkerboard and star.
Similarly, the checkerboard shares its $\mathbf{k}=(\pi,\pi)$ peak with the
striated phase, so we construct
$\tilde{\mathcal{F}}(\pi,\pi)-\tilde{\mathcal{F}}(0,\pi)$ as an order
parameter which is zero in the striated phase and nonzero only in
checkerboard.

Numerical Simulations of the 2D Phase Diagram
We numerically compute the many-body ground states at different points in the
$(\Delta/\Omega,R_{b}/a)$ phase diagram using the density-matrix
renormalization group (DMRG) algorithm White (1992, 1993), which operates in
the space of the so-called matrix product state (MPS) ansätze. While
originally developed for one-dimensional systems, DMRG can also be extended to
two dimensions by representing the 2D system as a winding 1D lattice
Stoudenmire and White (2012), albeit with long-range interactions. A major
limitation to two-dimensional DMRG is that the number of states required to
faithfully represent the ground-state wavefunction has to be increased
exponentially with the width of the system in order to maintain a constant
accuracy. For our calculations, we employ a maximum bond dimension of $1600$,
which allows us to accurately simulate $10\times 10$ square arrays Samajdar
_et al._ (2020a). We also impose open boundary conditions in both directions
and truncate the van der Waals interactions so as to retain up to third-
nearest-neighbor couplings. The numerical convergence criterion is set by the
truncation error, and the system is regarded to be well-converged to its true
ground state once this error drops below a threshold of $10^{-7}$. In
practice, this was typically found to be achieved after $\mathcal{O}(10^{2})$
successive sweeps.

Since the dimensions of the systems studied in Figure 4, (13 $\times$ 13
(experimentally) and 9 $\times$ 9 (numerically), are both of the form $(4n+1)$
$\times$ $(4n+1)$, the two phase diagrams are expected to be similar. In
particular, both these system sizes are compatible with the commensurate
ordering patterns of the crystalline phases observed in this work, and can
host all three phases (at the appropriate $R_{b}$/a) with the same boundary
conditions. Likewise, for extraction of the QCP, we use a 10$\times$10 array
as it is the largest numerically accessible square lattice comparable to the
16$\times$16 array used in our study of the quantum phase transition.

Mean-Field Wavefunction for the Striated Phase
To understand the origin of the striated phase, it is instructive to start
from a simplified model in which we assume that nearest-neighbor sites are
perfectly blockaded. Since we always work in a regime where $R_{b}/a>1$, this
model should also capture the essential physics of the full Rydberg
Hamiltonian.

In the classical limit of $\Omega=0$, the perfect checkerboard state has an
energy per site of $-\Delta/2+V(\sqrt{2}a)+V(2a)$, with $V(x)$ being the
interaction between sites at a distance $x$, whereas the corresponding energy
for the star-ordered state is $-\Delta/4$ (neglecting interactions for
$x>2a$). Accordingly, there is a phase transition between the checkerboard and
star phases when $\Delta=4[V(\sqrt{2}a)+V(2a)]$. On the other hand, for the
same density of Rydberg excitations, the striated phase has a classical energy
per site of $-\Delta/4+V(2a)/2$, which is always greater than that of the star
phase; hence, striated ordering never appears in the classical limit.

At finite $\Omega$, however, the striated phase emerges due to a competition
between the third-nearest-neighbor interactions and the second-order energy
shift upon dressing a ground state atom off-resonantly with the Rydberg state.
We can thus model the ground state of the striated phase as a product state,
where (approximately) $1/2$ of the atoms are in the ground state, $1/4$ of the
atoms are in the Rydberg state, and the remaining $1/4$ are in the ground
state with a weak coherent admixture of the Rydberg state. A general mean-
field ansatz for a many-body wavefunction of this form is given by

$\displaystyle\rvert\Psi_{\textsc{str}}(a_{1},a_{2})\rangle=$
$\displaystyle\bigotimes_{\mathbf{i}\in A_{1}}\left(\cos a_{1}\rvert
g\rangle_{\mathbf{i}}+\sin a_{1}\rvert r\rangle_{\mathbf{i}}\right)$ (5)
$\displaystyle\bigotimes_{\mathbf{i}\in A_{2}}\left(\cos a_{2}\rvert
g\rangle_{\mathbf{i}}+\sin a_{2}\rvert
r\rangle_{\mathbf{i}}\right)\bigotimes_{\mathbf{j}\in B}\rvert
g\rangle_{\mathbf{j}},$

where $A_{1}$ and $A_{2}$ represent the two sublattices of the (bipartite) $A$
sublattice, and $a_{1,2}$ are variational parameters. If $a_{1}=a_{2}$, then
our trial wavefunction simply represents a checkerboard state, but if
$a_{1}\neq a_{2}$, this state is not of the checkerboard type, and leads to
the striated phase.

Based on this ansatz, we can now explicitly see how the striated phase may
become energetically favorable in the presence of a nonzero $\Omega$. Consider
the atoms on the partially excited sublattice to be in the superposition
$\rvert g\rangle+[\Omega/\\{4V(\sqrt{2}a)-\Delta\\}]\rvert r\rangle$; this
describes the state of the atoms on the $(1,1)$ sublattice in the notation of
Fig. 5. The net energy per site of the system is then

$-\frac{\Delta}{4}+\frac{V(2a)}{2}-\frac{\Omega^{2}}{4\,(4V(\sqrt{2}a)-\Delta)}+\frac{\Omega^{2}\,V(\sqrt{2}a)}{2\,(4V(\sqrt{2}a)-\Delta)^{2}}$

where the third and fourth terms are the second-order energy shift and mean-
field interaction shift, respectively. From this expression, we observe that
if the energy gained from the dressing (these last two terms) is larger than
$V(2a)/2$, then the striated phase prevails over the star phase.

Dynamical Probe of the Striated Phase
We prepare striated ordering using an optimized cubic spline sweep along
$R_{b}/a=1.47$, ending at $\Delta/\Omega=2.35$. Immediately after this sweep,
the system is quenched to detuning $\Delta_{q}$ and relative laser phase
$\phi_{q}$. We quench at a lower Rabi frequency $\Omega_{q}=\Omega/4\approx
2\pi\times 1$ MHz to improve the resolution of this interaction spectroscopy.
For the chosen lattice spacing, the interaction energy between diagonal
excitations is $2\pi\times 5.3~{}$MHz. The reference phase for the atoms
$\phi=0$ is set by the instantaneous phase of the Rydberg coupling laser at
the end of the sweep into striated ordering. In the Bloch sphere picture,
$\phi=0$ corresponds to the $+x$ axis, so the wavefunctions on (0,0) and (1,1)
sublattices correspond to vectors pointing mostly up or mostly down with a
small projection of each along the $+x$ axis. In the same Bloch sphere
picture, quenching at $\phi_{q}=\pi/2$ or $-\pi/2$ corresponds to rotations
around the $+y$ or $-y$ axes (Fig. 5a).

To resolve the local response of the system, we use high-order correlators
which are extracted from single-shot site-resolved readout. In particular, we
define an operator $\hat{\mathcal{O}}_{i}^{(d)}$ on the eight atoms
surrounding site $i$. This operator projects the neighboring atoms into
configurations in which all four nearest atoms are in $\ket{g}$ and exactly
$d$ of the diagonal neighbors are in $\ket{r}$. Specifically, the operator
$\hat{\mathcal{O}}_{i}^{(d)}$ decomposes into a projector $\hat{A}_{i}$ on the
four nearest neighboring atoms and $\hat{B}_{i}^{(d)}$ on the four diagonal
neighbors, according to
$\hat{\mathcal{O}}_{i}^{(d)}=\hat{A}_{i}\hat{B}_{i}^{(d)}$. Defining
$\bar{n}_{i}=|g\rangle_{i}\langle g|$ and $n_{i}=|r\rangle_{i}\langle r|$, the
nearest neighbor projector is written as $\hat{A}_{i}=\prod_{\langle
j,i\rangle}\bar{n}_{j}$, where $\langle.\rangle$ denotes nearest neighbors.
The projector $\hat{B}_{i}^{(d)}$ sums over all configurations of the diagonal
neighbors (indexed $k_{1},k_{2},k_{3},k_{4}$) with $d$ excitations:

$\displaystyle\hat{B}_{i}^{(4)}$
$\displaystyle=n_{k_{1}}n_{k_{2}}n_{k_{3}}n_{k_{4}}$ (6)
$\displaystyle\hat{B}_{i}^{(3)}$
$\displaystyle=\bar{n}_{k_{1}}n_{k_{2}}n_{k_{3}}n_{k_{4}}+n_{k_{1}}\bar{n}_{k_{2}}n_{k_{3}}n_{k_{4}}+\ldots$
(7) $\displaystyle\hat{B}_{i}^{(2)}$
$\displaystyle=\bar{n}_{k_{1}}\bar{n}_{k_{2}}n_{k_{3}}n_{k_{4}}+\bar{n}_{k_{1}}{n}_{k_{2}}\bar{n}_{k_{3}}n_{k_{4}}+\ldots$
(8)

These operators are used to construct the conditional Rydberg density

$P^{(d)}=\frac{\sum_{i}\langle
n_{i}\hat{\mathcal{O}}_{i}^{(d)}\rangle}{\sum_{i}\langle\hat{\mathcal{O}}_{i}^{(d)}\rangle}$

which measures the probability of Rydberg excitation on site $i$ surrounded by
neighboring-atom configurations for which $\hat{\mathcal{O}}_{i}^{(d)}=1$.

To quantify coherences, we measure these conditional probabilities on their
corresponding resonances, after a fixed quench with variable quench phase
$\phi_{q}$. For a single particle driven by the Hamiltonian
$H=\Omega(\cos\phi_{q}\sigma_{x}+\sin\phi_{q}\sigma_{y})/2+\Delta\sigma_{z}/2$
for time $\tau$, the resulting Heisenberg evolution is given by
$\sigma_{z}^{\prime}=U^{\dagger}\sigma_{z}U$, where $U=e^{-iH\tau}$. The
resulting operator can be expressed as

$\displaystyle\sigma_{z}^{\prime}$ $\displaystyle=\tilde{\Omega}\sin
2\alpha(-\sigma_{x}\sin\phi_{q}+\sigma_{y}\cos\phi_{q})$ (9)
$\displaystyle+2\tilde{\Delta}\tilde{\Omega}\sin^{2}\alpha(\sigma_{x}\cos\phi_{q}+\sigma_{y}\sin\phi_{q})$
(10)
$\displaystyle+(\cos^{2}\alpha-(1-2\tilde{\Delta}^{2})\sin^{2}\alpha)\sigma_{z}$
(11)

where $\tilde{\Delta}=\Delta/\sqrt{\Delta^{2}+\Omega^{2}}$,
$\tilde{\Omega}=\Omega/\sqrt{\Delta^{2}+\Omega^{2}}$, and
$\alpha=\frac{1}{2}\tau\sqrt{\Delta^{2}+\Omega^{2}}$.

We fit the conditional probabilites $P^{(0)}$ and $P^{(4)}$ as a function of
$\phi_{q}$ (Fig. 5d,e), taking $\Delta$ as the effective detuning from
interaction-shifted resonance, and measuring
$\langle\sigma_{z}^{\prime}\rangle$ as a function of $\phi_{q}$ to extract the
Bloch vector components
$\langle\sigma_{x}\rangle,\langle\sigma_{y}\rangle,\langle\sigma_{z}\rangle$
on the two respective sublattices. For the (1,1) sublattice response, we model
the evolution averaged over random detunings, due to $\sim 15\%$ fluctuations
of the interaction shifts associated with thermal fluctuations in atomic
positions, which broaden and weaken the spectroscopic response. For both
sublattices we also include fluctuations in the calibrated pulse area ($\sim
10\%$ due to low power used). The extracted fit values are
$\sigma_{x,y,z}^{(0,0)}=-0.82(6),0.25(2),-0.32(4)$, and
$\sigma_{x,y,z}^{(1,1)}=-0.46(4),0.01(1),0.91(5)$.

Acknowledgements We thank many members of the Harvard AMO community,
particularly Elana Urbach, Samantha Dakoulas, and John Doyle for their efforts
enabling safe and productive operation of our laboratories during 2020. We
thank Hannes Bernien, Dirk Englund, Manuel Endres, Nate Gemelke, Donggyu Kim,
Peter Stark, and Alexander Zibrov for discussions and experimental help. We
acknowledge financial support from the Center for Ultracold Atoms, the
National Science Foundation, the Vannevar Bush Faculty Fellowship, the U.S.
Department of Energy, the Office of Naval Research, the Army Research Office
MURI, and the DARPA ONISQ program. T.T.W. acknowledges support from Gordon
College. H.L. acknowledges support from the National Defense Science and
Engineering Graduate (NDSEG) fellowship. G.S. acknowledges support from a
fellowship from the Max Planck/Harvard Research Center for Quantum Optics.
D.B. acknowledges support from the NSF Graduate Research Fellowship Program
(grant DGE1745303) and The Fannie and John Hertz Foundation. W.W.H. is
supported by the Moore Foundation’s EPiQS Initiative Grant No. GBMF4306, the
NUS Development Grant AY2019/2020, and the Stanford Institute of Theoretical
Physics. S.C. acknowledges support from the Miller Institute for Basic
Research in Science. R.S. and S.S. were supported by the U.S. Department of
Energy under Grant DE-SC0019030. The DMRG calculations were performed using
the ITensor Library Fishman  _et al._ (2020). The computations in this paper
were run on the FASRC Cannon cluster supported by the FAS Division of Science
Research Computing Group at Harvard University.

Author contributions S.E., T.T.W., H.L., A.K., G.S, A. O. and D.B. contributed
to the building experimental setup, performed the measurements, and analyzed
the data. Theoretical analysis was performed by R.S., H.P., W.W.H., and S.C.
All work was supervised by S.S., M.G., V.V., and M.D.L. All authors discussed
the results and contributed to the manuscript.

Competing interests M.G., V.V., and M.D.L. are co-founders and shareholders of
QuEra Computing. A.O. is a shareholder of QuEra Computing.

Correspondence and requests for materials should be addressed to M.D.L.

EXTENDED DATA FIG. 1: Large arrays of optical tweezers. The experimental
platform produces optical tweezer arrays with up to $\sim 1000$ tweezers and
$\sim 50\%$ loading probability per tweezer after $100~{}$ms of MOT loading
time. a. Camera image of an array of 34$\times$30 tweezers (1020 traps),
including aberration correction. b. Sample image of random loading into this
tweezer array, with 543 loaded atoms. Atoms are detected on an EMCCD camera
with fluorescence imaging. EXTENDED DATA FIG. 2: Correcting for aberrations in
the SLM tweezer array. The aberration correction procedure utilizes the
orthogonality of Zernike polynomials and the fact that correcting aberrations
increases tweezer light shifts on the atoms. To independently measure and
correct each aberration type, Zernike polynomials are added with variable
amplitude to the SLM phase hologram, with values optimized to maximize tweezer
light shifts. a. Two common aberration types: horizontal coma (upper) and
primary spherical (lower), for which $\sim 50$ milliwaves compensation on each
reduces aberrations and results in higher-depth traps. b. Correcting for
aberrations associated with the thirteen lowest order Zernike polynomials. The
sum of all polynomials with their associated coefficients gives the total
aberrated phase profile in the optical system, which is now corrected (total
RMS aberration of $\sim 70$ milliwaves). c. Trap depths across a $26\times 13$
trap array before and after correction for aberrations. Aberration correction
results in tighter focusing (higher trap light shift) and improved
homogeneity. Trap depths are measured by probing the light shift of each trap
on the $\ket{5S_{1/2},F=2}\to\ket{5P_{3/2},F^{\prime}=2}$ transition. d.
Aberration correction also results in higher and more homogeneous trap
frequencies across the array. Trap frequencies are measured by modulating
tweezer depths at variable frequencies, resulting in parametric heating and
atom loss when the modulation frequency is twice the radial trap frequency.
The measurement after correction for aberrations shows a narrower spectrum and
higher trap frequencies (averaged over the whole array). EXTENDED DATA FIG. 3:
Rearrangement protocol. a. Sample sequence of individual rearrangement steps.
There are two pre-sorting moves (1, 2). Move (3) is the single ejection move.
Moves (4-6) consist of parallel vertical sorting within each column, including
both upward and downwards move. The upper panel illustrates the frequency
spectrum of the waveform in the vertical and horizontal AODs during these
moves, with the underlying grid corresponding to the calibrated frequencies
which map to SLM array rows and columns. b. Spectrograms representing the
horizontal and vertical AOD waveforms over the duration of a single vertical
frequency scan during a realistic rearrangement procedure for a 26$\times$13
array. The heat-maps show frequency spectra of the AOD waveforms over small
time intervals during the scan. EXTENDED DATA FIG. 4: Generating homogeneous
Rydberg beams a. Measured Gaussian-beam illumination on the SLM for shaping
the 420-nm Rydberg beam. A Gaussian fit to this data is used as an input for
the hologram optimization algorithm. b. Corrected and measured wavefront error
through our optical system, showing a reduction of aberrations to
$\lambda/100$. c. Computer-generated hologram for creating the 420-nm top-hat
beam. d. Measured light intensity of the 420-nm top-hat beam (top), and the
cross section along where atoms will be positioned (bottom). Vertical lines
denote the 105-$\mu$m region where the beam should be flat. e. Using the
measured top-hat intensity, a phase correction is calculated for adding to the
initial hologram. f. Resulting top-hat beam after feedback shows significantly
improved homogeneity. EXTENDED DATA FIG. 5: Characterizing microwave-enhanced
Rydberg detection fidelity. The effect of strong microwave (MW) pulses on
Rydberg atoms is measured by preparing atoms in $\ket{g}$, exciting to
$\ket{r}$ with a Rydberg $\pi$-pulse, and then applying the MW pulse before
de-exciting residual Rydberg atoms with a final Rydberg $\pi$-pulse. (The
entire sequence occurs while tweezers are briefly turned off.) a. Broad
resonances are observed with varying microwave frequency, corresponding to
transitions from $\ket{r}=\ket{70S}$ to other Rydberg states. Note that the
transition to $\ket{69P}$ and $\ket{70P}$ are in the range of $10-12$ GHz, and
over this entire range there is strong transfer out of $\ket{r}$. Other
resonances might be due to multi-photon effects. b. With fixed $6.9$-GHz MW
frequency and varying pulse time, there is a rapid transfer out of the Rydberg
state on the timescale of several nanoseconds. Over short time-scales, there
may be coherent oscillations which return population back to $\ket{r}$, so a
$100$ ns pulse is used for enhancement of loss signal of $\ket{r}$ in the
experiment. EXTENDED DATA FIG. 6: Coarse-grained local staggered
magnetization. a. Examples of Rydberg populations $n_{i}$ after a faster (top)
and slower (bottom) linear sweep. b. Corresponding coarse-grained local
staggered magnetizations $m_{i}$ clearly show larger extents of
antiferromagnetically ordered domains (dark blue or dark red) for the slower
sweep (bottom) compared to for the faster sweep (top), as expected from the
Kibble-Zurek mechanism. c. Isotropic correlation functions $G^{(2)}_{m}$ for
the corresponding coarse-grained local staggered magnetizations after a faster
(top) or a slower (bottom). d. As a function of radial distance, correlations
$G^{(2)}_{m}$ decay exponentially with a length scale corresponding to the
correlation length $\xi$. The two decay curves correspond to faster (orange)
and slower (blue) sweeps.  EXTENDED DATA FIG. 7: Extracting the quantum
critical point. a. The mean Rydberg excitation density $\langle n\rangle$ vs.
detuning $\Delta/\Omega$ on a 16$\times$16 array. The data is fitted within a
window (dashed lines) to a cubic polynomial (red curve) as a means of
smoothening the data. b. The peak in the numerical derivative of the fitted
data (red curve) corresponds to the critical point $\Delta_{c}/\Omega=1.12(4)$
(red shaded regions show uncertainty ranges, obtained from varying fit
windows). In contrast, the point-by-point slope of the data (gray) is too
noisy to be useful. c. Order parameter $\tilde{\mathcal{F}}(\pi,\pi)$ for the
checkerboard phase vs. $\Delta/\Omega$ measured on a 16$\times$16 array with
the value of the critical point from b. superimposed (red line), showing the
clear growth of the order parameter after the critical point. d. DMRG
simulations of $\langle n\rangle$ vs. $\Delta/\Omega$ on a 10$\times$10 array.
For comparison against the experimental fitting procedure, the data from
numerics is also fitted to a cubic polynomial within the indicated window
(dashed lines). e. The point-by-point slope of the numerical data (blue curve)
has a peak at $\Delta_{c}/\Omega=1.18$ (blue dashed line), in good agreement
with the results (red dashed line) from both the numerical derivative of the
cubic fit on the same data (red curve) and the result of the experiment. f.
DMRG simulation of $\tilde{\mathcal{F}}(\pi,\pi)$ vs. $\Delta/\Omega$, with
the exact quantum critical point from numerics shown (red line). EXTENDED DATA
FIG. 8: Optimization of data collapse. a. Distance $D$ between rescaled
correlation length $\tilde{\xi}$ vs. $\tilde{\Delta}$ curves depends on both
the location of the quantum critical point location $\Delta_{c}/\Omega$ and on
the correlation length critical exponent $\nu$. The independently determined
$\Delta_{c}/\Omega$ (blue line, with uncertainty range in gray) and the
experimentally extracted value of $\nu$ (dashed red line, with uncertainty
range corresponding to the red shaded region) are marked on the plot. b. Our
determination of $\nu$ (red) from data collapse around the independently
determined $\Delta_{c}/\Omega$ (blue) is consistent across arrays of different
sizes. c-e. Data collapse is clearly better at the experimentally determined
value ($\nu=0.62$) as compared to the mean-field ($\nu=0.5$) or the (1+1)D
($\nu=1$) values. The horizontal extent of the data corresponds to the region
of overlap of all rescaled data sets.

## References

  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Zhang  _et al._ (2017) J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z. X. Gong, and C. Monroe, Nature 551, 601 (2017).
  * Gross and Bloch (2017) C. Gross and I. Bloch, Science 357, 995 (2017).
  * Choi  _et al._ (2016) J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Science 352, 1547 (2016).
  * Mazurenko  _et al._ (2017) A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature 545, 462 (2017).
  * Bernien  _et al._ (2017) H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).
  * Aleiner  _et al._ (2020) I. Aleiner _et al._ , arXiv:2012.00921 (2020).
  * Wang  _et al._ (2020) C. S. Wang, J. C. Curtis, B. J. Lester, Y. Zhang, Y. Y. Gao, J. Freeze, V. S. Batista, P. H. Vaccaro, I. L. Chuang, L. Frunzio, L. Jiang, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. X 10, 021060 (2020).
  * Arute  _et al._ (2019) F. Arute _et al._ , Nature 574, 505 (2019).
  * Zhong  _et al._ (2020) H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, Science 370, 1460 (2020).
  * Giovannetti  _et al._ (2011) V. Giovannetti, S. Lloyd, and L. Maccone, Nat. Photonics 5, 222 (2011).
  * Browaeys and Lahaye (2020) A. Browaeys and T. Lahaye, Nat. Phys. 16, 132 (2020).
  * Sachdev (2011) S. Sachdev,  _Quantum Phase Transitions_ , 2nd ed. (Cambridge University Press, Cambridge, UK, 2011).
  * Samajdar  _et al._ (2020a) R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Phys. Rev. Lett. 124, 103601 (2020a).
  * Pagano  _et al._ (2020) G. Pagano, A. Bapat, P. Becker, K. S. Collins, A. De, P. W. Hess, H. B. Kaplan, A. Kyprianidis, W. L. Tan, C. Baldwin, L. T. Brady, A. Deshpande, F. Liu, S. Jordan, A. V. Gorshkov, and C. Monroe, Proc. Natl. Acad. Sci. U.S.A. 117, 25396 (2020).
  * Friis  _et al._ (2018) N. Friis, O. Marty, C. Maier, C. Hempel, M. Holzäpfel, P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, R. Blatt, and B. Lanyon, Phys. Rev. X 8, 021012 (2018).
  * Harris  _et al._ (2018) R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao, Science 361, 162 (2018).
  * Labuhn  _et al._ (2016) H. Labuhn, D. Barredo, S. Ravets, S. de Léséleuc, T. Macrì, T. Lahaye, and A. Browaeys, Nature 534, 667 (2016).
  * Morgado and Whitlock (2020) M. Morgado and S. Whitlock, arXiv:2011.03031 (2020).
  * Lienhard  _et al._ (2018) V. Lienhard, S. de Léséleuc, D. Barredo, T. Lahaye, A. Browaeys, M. Schuler, L.-P. Henry, and A. M. Läuchli, Phys. Rev. X 8, 021070 (2018).
  * Guardado-Sanchez  _et al._ (2018) E. Guardado-Sanchez, P. T. Brown, D. Mitra, T. Devakul, D. A. Huse, P. Schauß, and W. S. Bakr, Phys. Rev. X 8, 021069 (2018).
  * Kim  _et al._ (2018) H. Kim, Y. Park, K. Kim, H.-S. Sim, and J. Ahn, Phys. Rev. Lett. 120, 180502 (2018).
  * de Léséleuc  _et al._ (2019) S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Science 365, 775 (2019).
  * Keesling  _et al._ (2019) A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 568, 207 (2019).
  * Madjarov  _et al._ (2020) I. S. Madjarov, J. P. Covey, A. L. Shaw, J. Choi, A. Kale, A. Cooper, H. Pichler, V. Schkolnik, J. R. Williams, and M. Endres, Nat. Phys. 16, 857 (2020).
  * Omran  _et al._ (2019) A. Omran _et al._ , Science 365, 570 (2019).
  * Graham  _et al._ (2019) T. M. Graham, M. Kwon, B. Grinkemeyer, Z. Marra, X. Jiang, M. T. Lichtman, Y. Sun, M. Ebert, and M. Saffman, Phys. Rev. Lett. 123, 230501 (2019).
  * Levine  _et al._ (2019) H. Levine, A. Keesling, G. Semeghini, A. Omran, T. T. Wang, S. Ebadi, H. Bernien, M. Greiner, V. Vuletić, H. Pichler, and M. D. Lukin, Phys. Rev. Lett. 123, 170503 (2019).
  * Madjarov  _et al._ (2019) I. S. Madjarov, A. Cooper, A. L. Shaw, J. P. Covey, V. Schkolnik, T. H. Yoon, J. R. Williams, and M. Endres, Phys. Rev. X 9, 041052 (2019).
  * Young  _et al._ (2020) A. W. Young, W. J. Eckner, W. R. Milner, D. Kedar, M. A. Norcia, E. Oelker, N. Schine, J. Ye, and A. M. Kaufman, Nature 588, 408 (2020).
  * Schymik  _et al._ (2020) K.-N. Schymik, V. Lienhard, D. Barredo, P. Scholl, H. Williams, A. Browaeys, and T. Lahaye, arXiv:2011.06827 (2020).
  * Kumar  _et al._ (2018) A. Kumar, T.-Y. Wu, F. Giraldo, and D. S. Weiss, Nature 561, 83 (2018).
  * Ohl de Mello  _et al._ (2019) D. Ohl de Mello, D. Schäffner, J. Werkmann, T. Preuschoff, L. Kohfahl, M. Schlosser, and G. Birkl, Phys. Rev. Lett. 122, 203601 (2019).
  * Barredo  _et al._ (2018) D. Barredo, V. Lienhard, S. de Léséleuc, T. Lahaye, and A. Browaeys, Nature 561, 79 (2018).
  * Barredo  _et al._ (2016) D. Barredo, S. de Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, Science 354, 1021 (2016).
  * Jaksch  _et al._ (2000) D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Phys. Rev. Lett. 85, 2208 (2000).
  * Lukin  _et al._ (2001) M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 87, 037901 (2001).
  * Zurek  _et al._ (2005) W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 95, 105701 (2005).
  * Felser  _et al._ (2020) T. Felser, S. Notarnicola, and S. Montangero, arXiv:2011.08200 (2020).
  * Huang  _et al._ (2020) H.-Y. Huang, R. Kueng, and J. Preskill, Nat. Phys. 16, 1050 (2020).
  * Saskin  _et al._ (2019) S. Saskin, J. T. Wilson, B. Grinkemeyer, and J. D. Thompson, Phys. Rev. Lett. 122, 143002 (2019).
  * Anderegg  _et al._ (2019) L. Anderegg, L. W. Cheuk, Y. Bao, S. Burchesky, W. Ketterle, K.-K. Ni, and J. M. Doyle, Science 365, 1156 (2019).
  * Liu  _et al._ (2018) L. R. Liu, J. D. Hood, Y. Yu, J. T. Zhang, N. R. Hutzler, T. Rosenband, and K.-K. Ni, Science 360, 900 (2018).
  * Turner  _et al._ (2018) C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Nat. Phys. 14, 745 (2018).
  * Surace  _et al._ (2020) F. M. Surace, P. P. Mazza, G. Giudici, A. Lerose, A. Gambassi, and M. Dalmonte, Phys. Rev. X 10, 021041 (2020).
  * Bluvstein  _et al._ (2020) D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletić, and M. D. Lukin, submitted for publication (2020).
  * Samajdar  _et al._ (2020b) R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, arXiv:2011.12295 (2020b).
  * Verresen  _et al._ (2020) R. Verresen, M. D. Lukin, and A. Vishwanath, arXiv:2011.12310 (2020).
  * Bañuls  _et al._ (2020) M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. Van Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, Eur. Phys. J. D 74, 165 (2020).
  * Notarnicola  _et al._ (2020) S. Notarnicola, M. Collura, and S. Montangero, Phys. Rev. Research 2, 013288 (2020).
  * Weimer  _et al._ (2010) H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nat. Phys. 6, 382 (2010).
  * Auger  _et al._ (2017) J. M. Auger, S. Bergamini, and D. E. Browne, Phys. Rev. A 96, 052320 (2017).
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, arXiv:1411.4028 (2014).
  * Zhou  _et al._ (2020) L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Phys. Rev. X 10, 021067 (2020).
  * Wild  _et al._ (2020) D. S. Wild, D. Sels, H. Pichler, and M. D. Lukin, arXiv:2005.14059 (2020).
  * Kim  _et al._ (2019) D. Kim, A. Keesling, A. Omran, H. Levine, H. Bernien, M. Greiner, M. D. Lukin, and D. R. Englund, Opt. Lett. 44, 3178 (2019).
  * Endres  _et al._ (2016) M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Science 354, 1024 (2016).
  * Lee  _et al._ (2017) W. Lee, H. Kim, and J. Ahn, Phys. Rev. A 95, 053424 (2017).
  * Sheng  _et al._ (2020) C. Sheng, J. Hou, X. He, P. Xu, K. Wang, J. Zhuang, X. Li, M. Liu, J. Wang, and M. Zhan, arXiv:2011.10390 (2020).
  * Levine  _et al._ (2018) H. Levine, A. Keesling, A. Omran, H. Bernien, S. Schwartz, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Phys. Rev. Lett. 121, 123603 (2018).
  * Bowman  _et al._ (2017) D. Bowman, T. L. Harte, V. Chardonnet, C. D. Groot, S. J. Denny, G. L. Goc, M. Anderson, P. Ireland, D. Cassettari, and G. D. Bruce, Opt. Express 25, 11692 (2017).
  * Zupancic  _et al._ (2016) P. Zupancic, P. M. Preiss, R. Ma, A. Lukin, M. E. Tai, M. Rispoli, R. Islam, and M. Greiner, Opt. Express 24, 13881 (2016).
  * Beterov  _et al._ (2009) I. I. Beterov, I. I. Ryabtsev, D. B. Tretyakov, and V. M. Entin, Phys. Rev. A 79, 052504 (2009).
  * Bhattacharjee and Seno (2001) S. M. Bhattacharjee and F. Seno, J. Phys. A: Math. Gen. 34, 6375 (2001).
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
  * Stoudenmire and White (2012) E. M. Stoudenmire and S. R. White, Annu. Rev. Condens. Matter Phys. 3, 111 (2012).
  * Fishman  _et al._ (2020) M. Fishman, S. R. White, and E. M. Stoudenmire, arXiv:2007.14822 (2020).