Observing the space- and time-dependent growth of correlations in dynamically tuned synthetic Ising antiferromagnets

We explore the dynamics of artificial one- and two-dimensional Ising-like quantum antiferromagnets with different lattice geometries by using a Rydberg quantum simulator of up to 36 spins in which we dynamically tune the parameters of the Hamiltonian. We observe a region in parameter space with antiferromagnetic (AF) ordering, albeit with only finite-range correlations. We study systematically the influence of the ramp speeds on the correlations and their growth in time. We observe a delay in their build-up associated to the finite speed of propagation of correlations in a system with short-range interactions. We obtain a good agreement between experimental data and numerical simulations taking into account experimental imperfections measured at the single particle level. Finally, we develop an analytical model, based on a short-time expansion of the evolution operator, which captures the observed spatial structure of the correlations, and their build-up in time.


I. INTRODUCTION
The study of non-equilibrium dynamics is currently one of the most challenging areas of quantum many-body physics. In contrast to the equilibrium case, where statistical physics provides a general theoretical framework and where very powerful numerical methods are available, the out-of-equilibrium behavior of quantum matter presents a wide variety of phenomena and is extremely hard to simulate numerically, especially in dimensions d > 1. When the parameters of a quantum many-body system are quenched abruptly, or more generally ramped with a finite rate into a new quantum phase, correlations and entanglement build up and propagate over the system, which may, at long times, either thermalize or retain memory of its initial state when many-body localization occurs [1,2]. The speed at which the correlations propagate depends on the range of the interaction and is limited by the Lieb-Robinson bounds [3][4][5][6][7][8].
An attractive way to study this physics has emerged in the last years, and consists in using quantum simulators, i.e. wellcontrolled, artificial quantum systems that implement experimentally the Hamiltonian of interest [9]. Spin Hamiltonians that are used in condensed-matter physics to describe e.g. quantum magnets are arguably the simplest quantum many-body systems that can be used to study non-equilibrium dynamics: even though they involve distinguishable particles with only internal degrees of freedom, the interplay between interactions, geometry and dimensionality provides a wealth of distinct quantum phases into which the system can be driven. In recent years, many experimental platforms for the quantum simulation of spin Hamiltonians have been developed using the tools of atomic physics. For example, equilibrium properties of synthetic quantum magnets have been studied using trapped ions [10,11] or ultra-cold atoms in optical lattices [12][13][14][15][16] including e.g. the observation of long-range antiferromagnetic order [17]. Many experiments using these platforms were also devoted to the study of non-equilibrium dynamics, including the investigation of the Lieb-Robinson bound [18][19][20][21][22][23][24].
Recently, a new platform, using arrays of individually resolved atoms excited to Rydberg states, has been shown to provide a versatile way to engineer synthetic quantum Ising magnets [25]. Pioneering experiments in quantum gas microscopes studied quenches [26] or slow sweeps [27] in a regime where the blockade radius R b , i.e. the distance over which interatomic interactions prevent the excitation of two atoms, was much larger than the lattice spacing a, rendering the underlying lattice hardly relevant. In this case, the observed correlations are liquid-like, and observing the crystallike ground state of the system [28] would require exponentially long ramps [29]. More recently, experiments with arrays of optical tweezers allowed exploring the regime R b a, studying non-equilibrium dynamics following quenches [30] or slow sweeps [31].
Here, we use a Rydberg-based platform emulating an Ising antiferromagnet to study the growth of correlations during ramps of the experimental parameters in 1d and 2d arrays of up to 36 single atoms with different geometries. We operate in the regime R b a where the interactions act to a good approximation only between nearest neighbors. We dynamically tune the parameters of the Hamiltonian and observe the build-up of antiferromagnetic order. We also observe the influence of the finite ramp speed on the extent of the correlations, and follow the development in space and time of these correlations during a ramp. Numerical simulations of the dynamics of the system without any adjustable parameters are in very good agreement with the experimental data, and show that single-particle dephasing arising from technical imperfections currently limits the range of the observed correlations. Finally, we observe a characteristic spatial structure in the correlations, which can be understood qualitatively by a short-time expansion of the evolution operator for both, square and triangular lattices.
The paper is organized as follows. In Sec. II, we describe the experimental platform. After recalling the phase diagram for the different array geometries (Sec. III), we explore it for the square array (Sec. IV). In Sec. V A we study the influence of ramp speeds on the correlations. In Sec. V B, we observe a delay in the build-up of the spin-spin correlations, a feature linked to their finite speed of propagation. Finally, in Sec. V C, we analyze the 2d spatial structure of the AF correlations on the square and triangular geometries and show that it is qualitatively captured by an analytical model based on short-time expansion.

II. EXPERIMENTAL PLATFORM
Our experimental platform (see Appendix A) is based on user-defined two-dimensional arrays of optical tweezers, each containing a single 87 Rb atom [30]. Here we use the arrays shown in Fig. 1(a) containing up to N = 36 atoms: a 1d chain with periodic boundary conditions (PBC), a square lattice, and a triangular lattice. We achieve full loading of the arrays using our atom-by-atom assembler [32]. The atoms are prepared in the ground state |↓ = 5S 1/2 , F = 2, m F = 2 by optical pumping, and then coupled coherently to the Rydberg state |↑ = 64D 3/2 , m j = 3/2 with a two-photon transition of Rabi frequency Ω and a detuning δ, while the traps are switched off. The system is described by the Hamiltonian: where n i = |↑ ↑| i is the projector on the Rydberg state for atom i, and σ x = |↑ ↓| + |↓ ↑| is the x-Pauli matrix. The interaction term U ij arises from van der Waals interactions between the atoms, and thus scales as 1/r 6 ij with the distance r ij between atoms i and j. This short-range character allows us to neglect interactions beyond nearest-neighbor (NN) atoms for this work [33], and we thus restrict Eq. (1) to NN terms only. For the D states we use, the van der Waals interaction is anisotropic [25,34], and the lattice spacings in the arrays are tuned such that the NN interactions anisotropy is below 10%. We use typical values Ω max /(2π) ∼ 2 MHz and U/h ∼ 1 − 3 MHz, see Table I. The driving parameters Ω and δ can be considered being constant over the entire array (see Appendix A).
The system thus realizes an Ising-like model with a transverse field ∝ Ω and a longitudinal field ∝ δ, which gives rise to AF order for U > 0 (see below). We probe the system by using time-dependent ramps Ω(t) and δ(t) as shown in Fig. 1(b). The Rabi frequency Ω is switched on and off in t rise , t fall ∼ 250 ns at a constant detuning δ, and, in between, δ is ramped linearly from δ 0 /(2π) = −6 MHz to δ final during the time t sweep . The total duration of the ramp is then t tot = t rise + t sweep + t fall . After this, the trap array is switched on again. Atoms in |↓ are observed by fluorescence, while those in |↑ are lost from their trap. For a given set of parameters, the experiment is repeated a few hundred times to reconstruct quantities of interest such as the Rydberg density (equivalent to the magnetization), spin-spin correlation functions, or the sublattice density histogram. When δ final lies in the AF region (see Fig. 1c), correlation functions show the emergence of short-range order [ Fig. 1(d)] (see more details on the measurements in Sec IV).

III. THEORETICAL PHASE DIAGRAMS AND STATE PREPARATION CONSIDERATIONS
The calculated ground state phase diagrams for the nearestneighbor Ising model are shown in Fig. 1(c) for the three geometries considered in this paper. On the 1d chain, the phase diagram is well-known and features an AF phase and a paramagnetic phase delimited by a second-order quantum phase transition line of the (1+1)d Ising universality class [35]. For ( δ/U ) TFI = z/2, with the coordination number z = 2, the Hamiltonian corresponds exactly to the analytically solvable transverse field Ising (TFI) model without longitudinal field, where in 1d the critical point ( Ω/U ) c = 1/2 is known analytically.
The phase diagram for the NN Hamiltonian on the square lattice (z = 4) is qualitatively similar to the phase diagram on the chain, with again an AF phase and a paramagnetic phase delimited by a second-order quantum phase transition line in the (2+1)d Ising universality class. For the transverse field Ising line ( δ/U ) TFI , the critical point is known to high precision from Monte-Carlo simulations, ( Ω/U ) c = 1.52219(1) [36]. Finally, on the triangular lattice (z = 6) the NN Hamiltonian features a much richer phase diagram. A Rydberg crystal with filling fraction 1/3 (at Ω = 0), where one of the three triangular sublattices is occupied by Rydberg states and the atoms in the other sublattices remain in the ground state, appears in a region within 0 < δ/U < z/2. The conjugate crystals obtained by flipping all the spins (filling fraction 2/3 at Ω = 0) lies in the region within z/2 < δ/U < z. The most interesting feature occurs when δ/U = z/2, where an "order by disorder" process occurs (see Appendix B). When including the open boundary conditions for the square and triangular N = 36 arrays studied here, the "phase diagrams" present the same phases with some modifications (see Ap-pendix C).
In the experiments reported here, we initialize the system in the product state with all atoms in their ground state |↓ . Since the system is not at equilibrium with a thermal bath, we cannot simply cool to the ground state. Therefore we use an experimental protocol involving sweeps of Ω(t) and δ(t) in order to characterize the ground states. Ideally, we would use adiabatic state preparation to reach the targeted ground state. In order for this approach to succeed, the duration of the sweep, which scales as the inverse of the square of the energy gap ∆ above the instantaneous ground states [37], should be smaller than the coherence time of the system. In the 1d chain and the 2d square lattice, the minimal gap at the quantum phase transition scales as the inverse of the linear size of the system: ∆ ∼ 1/N in 1d and ∆ ∼ 1/L = 1/ √ N in 2d [35,38]. The gap also depends on the excitation velocity at the quantum critical point, which can vary substantially along the quantum phase transition lines. For the triangular lattice, we expect a first-order quantum phase transition between the paramagnet to either the 1/3 or 2/3 filling Rydberg crystals, resulting in a minimal gap exponentially small in N [39]. Due to these scalings, the gaps are small for the number of atoms used here. As a consequence adiabatic state preparation would require long pulses. In the absence of any imperfections such as dephasing, pulse durations t tot of a few µs would allow to reach strong correlations extending over the entire system. While such durations are experimentally accessible, state-of-the-art platforms [27,30,31] show significant dephasing over these timescales. For this reason, we approach here the question of state preparation from the opposite side: we ask how much correlations and what kind of structures we can expect being built up in a given amount of time. Answering these questions also informs us about the minimal time required to build up highly correlated states starting from a product state [4], i.e. about a quantum speed limit in our many-body system.

IV. EXPLORING THE SQUARE LATTICE PHASE DIAGRAM
In a first set of experiments we map out the Ω = 0 "phase diagram" on a L × L square array with L = 4 and L = 6. In order to do so, we use the ramps shown in Fig. 1(b) with Ω max /U = 2.3. This value is larger than the critical point ( Ω/U ) c ≈ 1.5, such that the quantum phase transition line is crossed while ramping down Ω, see Fig. 1(c). From the analysis of the final fluorescence images, we reconstruct the Rydberg density n = i n i /N , and the connected spinspin correlation function where the sum runs over atom pairs (i, j) whose separation is r i − r j = (ka, la), and N k,l is the number of such atom pairs in the array [40]. For a perfect antiferromagnetic Néel ordering, g (2) takes the values ±1/4 for |k|+|l| even and odd, respectively. Figure 2(a) shows the spin-spin correlation function at the end of the ramp as a function of x = δ final /U . We observe strong AF correlations, i.e. the sign of g (2) (k, l) changes according to the parity of the Manhattan distance |k| + |l|, in the region 0 < x < 4 where AF order is expected, while the correlations vanish outside of this region. The amplitude of the AF correlations decreases with distance, in a way which is well captured by an exponential decay with a correlation length ξ, defined by g (2) Repeating the same experiment with a 1d chain yields the same correlation length (the particularity of the triangular lattice is discussed in Appendix G). Importantly, even though the correlation length is smaller than two sites for both the 1d chain and the square lattice, we are able to detect finite correlations with the expected sign structure for up to five Manhattan shells, i.e., almost over the whole array.
Another way to highlight AF correlations is to partition the array into the two Néel sublattices A and B and plot a twodimensional histogram P (n A , n B ) with the |↑ populations n A and n B of each sublattice as axes. For a perfect AF ordering, one would observe populations only in the two corners (n A = 0, n B = N/2) and (n A = N/2, n B = 0). The experimental results in Fig. 2(b) show that for x = 2.5 (central plot), the sublattice population histogram is substantially elongated along the anti-diagonal n B = N/2 − n A , which is not observed for x < 0 or x > 4. For comparison, Fig. 2(c) shows the corresponding histograms that would be obtained for an uncorrelated random state with the same average Rydberg density: there the elongation along the anti-diagonal is absent.
In Fig. 2(d,e) we locate the boundaries of the AF phase. Panel (d) shows the mean Rydberg density n. For the Ω = 0 ground state of Eq. (1) in the NN limit with periodic boundary conditions, it should rise in steps, from 0 for x < 0, to 1/2 for 0 < x < 4, and to 1 for x > 4 [41]. For open boundary conditions, additional steps are present in the 0 < x < 4 region, see Appendix C. Experimentally, the curve n(x) varies continuously due to the finite duration of the sweep, as also observed in Refs. [20,27]. This mean density is therefore not a good observable to differentiate the paramagnetic and AF regions. Instead, we introduce the Néel structure factor S Néel to detect antiferromagnetic correlations in the AF region of the phase diagram: This quantity is an estimator for the correlation volume, i.e., the number of spins correlated antiferromagnetically with a given spin. In a situation with true long range order, S Néel diverges linearly with the total "volume" N of the system, while it stays almost constant for short-range ordered correlations when the system sizes are larger than the correlation length. As shown in Fig. 2(e), this quantity indicates the presence of substantial short-range AF correlations for 0 < x < 4, as expected from the phase diagram. The results presented so far demonstrate that by analyzing the correlation functions we can locate phase boundaries in the phase diagrams. In the next section, we will study how the correlations build up in time.

V. EXPLORING THE TIME-AND SPACE-DEPENDENCE OF CORRELATIONS
In the following we first investigate in Sec. V A how S Néel depends on the ramp speed and find a duration t tot that maximizes its value. Then, with these settings, we study in Sec. V B the temporal build-up of correlations during the optimized sweep and observe delays between the growth of correlations at increasing distances. Eventually, in Sec. V C, we analyze the spatial structure of these correlations and give a qualitative explanation of the observed patterns in terms of a short-time expansion.

A. Correlations after ramps of varying durations
We perform experiments on a 36-site square lattice using the sweep shown in Fig. 1(b) with fixed parameters δ final /U = 1.7, Ω max /U = 0.7 < ( Ω/U ) c . In contrast to the previous experiment, the quantum phase transition line is crossed while ramping the detuning δ and we vary the duration of the sweep t sweep (and therefore t tot ). The initial parameters are chosen such that the ground state at t = 0 corresponds to |↓↓ · · · (no Rydberg excitation) and the ground state at the final parameters is an AF [42]. Fig. 3(a) shows two-dimensional plots of the experimental g (2) correlation functions for four different values of t tot . For the shortest sweeps, the correlations remain weak. For intermediate durations (t tot ∼ 1 µs) strong correlations emerge, with a staggered structure extending over many Manhattan shells, as shown in Fig. 3(b). For even larger sweep durations, the correlation signal decreases with increasing sweep duration. In Fig. 3(c) we show the value of S Néel as a function of t tot : starting from small values for short ramps, it exhibits a broad maximum for t tot ∼ 1 µs and slowly decreases for longer ramp durations.
In order to gain a better understanding of the behavior de-  scribed above, we compare the data with numerical simulations of the dynamics of the system governed by the Hamiltonian of Eq. (1) with the full van der Waals interactions for the sweeps used in the experiment [43]. We show in Fig. 3(c) the results of the simulation on a 4 × 4 square lattice. The green line corresponds to a unitary evolution which models well the experiment only for short sweep durations. The dashed yellow line includes the decoherence observed in the experiment on the excitation of one atom through an empirical dephasing model (see Appendix D 2). We observe a remarkable agreement between this local (i.e. single atom) dephasing model with no adjustable parameters and the experiment over a large range of t tot . This indicates that the saturation and the decay of the correlations for longer sweeps is due to (single-particle) decoherence.
We now analyze the spin-spin correlations for different Manhattan distances |k| + |l|. In Fig. 3(d) we observe the build up of correlations up to the fourth shell |k| + |l| = 4, all of them being antiferromagnetically staggered, g (2) (k, l) ∝ (−1) |k|+|l| . For short sweeps the correlations sharply increase with increasing duration, saturate for longer sweeps, and decay, again due to decoherence. The simulation including dephasing (dashed lines) reproduces well this trend.

B. Build-up of correlations along a ramp of fixed duration
In a subsequent experiment, we analyze how the AF correlations build up in time, during the sweep that maximizes S Néel (t tot = 1.0 µs) identified in the previous subsection. Stopping the dynamics after a variable time 0 < t < t tot , Fig. 4 shows the time evolution of the correlations after entering the AF region identified in Fig. 2 at δ(t) > 0 (t > 0.5 µs). We observe that most of the correlations build up for 0.6 < t < 0.8 µs and then freeze at larger times. When comparing again to the numerical simulation of the dynamics including the singleparticle dephasing, we obtain a remarkable agreement with the data. To quantify this effect we normalize the correlations for each distance such that the corresponding dephasing-free simulation reaches a maximum of 1 at large time [ Fig. 4(b)]. Fixing an arbitrary threshold of 0.2, we observe that the nearest neighbor shell (m = 1) reaches this threshold at t ≈ 0.64 µs, the second shell (m = 2) at t ≈ 0.71 µs and finally the third shell (m = 3) at t ≈ 0.79 µs (see gray vertical lines). This delay is a manifestation of the finite propagation speed of correlations as theoretically predicted by Lieb-Robinson bounds [3,4]. Such bounds are well explored in the context of quench dynamics [18,20,23,44,45], but their importance for state preparation protocols is not equally well known (see Appendix E for a brief review on Lieb-Robinson bounds adapted for the present context). Fig. 4(b) also reveals that when the correlations on the first shell reach the threshold, correlations for higher Manhattan distances are suppressed, but still detectable for m = 2. This illustrates the known fact that the correlations outside the Lieb-Robinson cone are finite, albeit exponentially suppressed (see also Appendix E). To recover this fact, we introduce a simple analytical approach based on a short-time expansion method (see Appendix F for details and specific results). We find analytical expressions in powers of the duration T of the ramp for the connected g (2) correlation functions for different Manhattan distance m valid in the limit (U T / , ΩT, δT 1). We obtain nearest-neighbor correlations (m = 1) of order T 6 and next-nearest neighbor correlations (m = 2) of order T 10 . More generally, the leading order of the expansion for two sites separated by a Manhattan distance m appears at order T 2+4m thus suggesting that at a given (short) time the correlation decreases exponentially with m (see more details in Appendix F). This scaling explains qualitatively the observed time dependence in Figs. 4(a) and (b), even though strictly speaking the range of applicability of the short-time expansion is limited to times much shorter than those shown there (for which the correlations would still be extremely small, and thus very challenging to measure experimentally). However, our numerical exact diagonalization results without and with dissipation as well as the experiment show that the qualitative features of the short-time expansion prevail for the considered non-perturbative times and the actual ramp shapes.

C. Build-up of spatial structures on the square and triangular lattices
We finally analyze the spatial structure of the correlations in more detail. Both Fig. 3(d) and Fig. 4(a,b) show that the correlations do not depend only on the Manhattan distance m = |k| + |l|, but also, for fixed m, on k − l. For instance, for second neighbors (m = 2), the correlations for (k, l) = (0, 2) or (2, 0) is about half of those along the diagonal (k, l) = (1, 1). This spatial structure, absent in a 1d setting, has not, to our knowledge, been experimentally observed in 2d.
The observed spatial structures in the correlations within a given Manhattan shell m can also be captured by the shorttime expansion. The leading order coefficient of a given correlator g (2) (k, l) depends on the number of paths on the lattice linking the two sites as detailed in Appendix F 2. Fig. 5 shows a comparison of the spatial structure in experimental data for both a square and a triangular lattice after a ramp of finite duration. In panels (a) and (c), the number of linking paths is given for k ≥ l ≥ 0 by the binomial coefficient C m k−l for both lattices and is shown in brackets. In panels (b) and (d) we analyze the correlations on the square and triangular lattice. The green dotted lines show the k − l dependence of g (2) (k, l) assumed to be given only by the binomial coefficients, normalized to the maximal correlation value for each subplot.
The precision of the experiment even allows us to observe a small asymmetry in each shell m = const., due to a slight residual anisotropy of the interaction. An extension of the short-time expansion to the anisotropic case yields an analytical expression from which we extract the ratio U z /U w of the nearest-neighbor interactions. We use this interaction anisotropy in the analytical expression to obtain the correlation for larger |k| + |l|. The results are shown as red dashed lines.
Interestingly, for triangular lattices, in contrast with the case of square arrays, the spatial structures of the correlations that one observes experimentally, and which are well reproduced by the short-time expansion, do not reflect directly those one would obtain in the ground state (see Appendix G). A detailed study of this specific behavior, which may be a dynamical signature of the frustrated character of the triangular geometry, is however beyond the scope of the present work.

VI. CONCLUSION AND OUTLOOK
In conclusion, we have used a Rydberg-based platform to study antiferromagnetic correlations in a synthetic Ising magnet with different geometries. Using dynamical variations of the parameters, we explored the phase diagram and prepared arrays exhibiting antiferromagnetic order with pronounced Manhattan structures. We also studied the growth of the correlations during the sweeps of the parameters. We observed delays in the build-up of correlations between sites at different distances, a feature linked to the Lieb-Robinson bounds for the propagation of correlations in a system with nearestneighbor interactions. We were able to understand the spatial structure of the correlations after short to intermediate evolution times using an analytical short-time expansion of the evolution operator. Finally, we obtained remarkable agreement between the data on the dynamics of the correlations and numerical simulations using a local (i.e. single particle) dephasing model.
In the future, the time over which coherent simulation can be performed in such artificial quantum magnets could be extended by a better understanding of the dephasing mechanisms at play in the coherent excitation of Rydberg states [46]. We will then be able to explore geometric frustration on triangular or Kagome lattices [47,48]. Another promising avenue is the study of magnets described by the XY model implemented using resonant dipole-dipole interactions [49,50]: there the coherent drive is obtained by using microwaves, and much lower dephasing rates are expected; moreover the interaction decays slowly, as 1/r 3 , and long-range effects should be more prominent [51].
Note: Similar antiferromagnetic correlations in 2d have been observed recently at Princeton University, using a system of Li Rydberg atoms in a quantum gas microscope. Here we give more details about our experimental setup and on the mapping of the Rydberg system onto an Ising antiferromagnet.
FIG. 6. Influence of the sign of C6 on the many-body spectrum.
Eigenenergies of a small 2 × 2 square system as a function of δ for Ω/U 0.5. (a) For a repulsive interaction C6 > 0, the ground state of H is antiferromagnetic for 0 < δ < 2U . It can be reached from |↓ ⊗4 by adiabatically following the ground state starting from negative values of δ (trajectory in dashed arrows). (b) For an attractive interaction C6 < 0, the antiferromagnetic state is the most excited state, and the sign of the detuning needed to reach it from |↓ ⊗4 is reversed.
We use a two-photon excitation scheme to excite the atoms to the Rydberg state. The two lasers at 795 nm and 475 nm, with corresponding Rabi frequencies Ω r and Ω b , and a detuning ∆ 2π×740 MHz from the intermediate 5P 1/2 , F = 2 state result in an effective Rabi frequency Ω = Ω r Ω b /(2∆) for the coupling between |↓ and |↑ , but also introduce lightshifts (Ω 2 r −Ω 2 b )/(4∆) that add up to the two-photon detuning δ. In order to generate the ramps shown in Fig. 1(b), we thus compensate for these additional lightshifts. An AOM is used for changing dynamically the amplitude and frequency of the red beam. Due to the finite size of our excitation beams, the atoms do not all experience exactly the same δ and Ω. For our arrays with a size ∼ 40 µm, the inhomogeneities of Ω are below 15%, while the detuning does not differ from its value on the central atom by more than 150 kHz.
The detection of the state of each atom relies on the loss of atoms in the Rydberg state, which are not recaptured in the optical tweezers. This detection method is thus subject to small detection errors [46]. In particular, an atom actually in |↓ has a small probability ε to be lost and thus incorrectly inferred to be in |↑ . As in [30], we measure ε in a calibration experiment, and then include the effect of detection errors on all the theoretical curves.
An antiferromagnetic phase is expected when interaction between parallel spins are repulsive, corresponding to the case where U > 0 in the Hamiltonian Eq. (1). Fig. 6 (a) shows the energy levels for a 2×2 square matrix of atoms as a function of the detuning δ, with U > 0 and Ω 0.5 U . The lowest energy configuration with two |↑ spins is then an antiferromagnetic state for 0 < δ < 2U . Then, starting from a detuning δ < 0 with all atoms in |↓ , and ramping it up to δ ≈ U/ (dashed arrow), the lowest energy state is evolving adiabatically to the antiferromagnetic state. In our experiment, the interaction between |↑ = 64D 3/2 , m j = 3/2 spin states are in fact attractive (U < 0). Consequently, the antiferromagnetic state is never the lowest energy configuration whatever the value of δ, see Fig. 6 (b). Nevertheless, starting from a detuning δ > 0 and ramping it down to δ ≈ U/ , the system evolves from the initial state to an antiferromagnetic state while following the highest-energy level. In this case, preparing the antiferromagnetic state actually means obtaining the ground state of −H, hence the change of sign of the detunings needed to reach the correlated phase. For our parameters, the van der Waals interaction is not only attractive, but it is also slightly anisotropic [25,34,52], being about three times as small in the horizontal direction as in the vertical one (alonĝ z). We compensated for this difference by slightly distorting the arrays alongẑ by a factor ∼ 3 1/6 . Appendix B: Description of the "order-by-disorder" process As mentioned in Sec. III, the triangular lattice presents an interesting feature for δ/U = z/2: there the model is fully frustrated in the classical limit; the interactions on the bonds of each single triangle of the lattice cannot be simultaneously satisfied under the given density constraint n = 1/2 leading to an extensive ground-state degeneracy [53]. Upon switching on the transverse field, this degeneracy is lifted by a process called "order-by-disorder" (OBD), i.e. the addition of "disorder" (here the quantum fluctuations provided by the transverse field Ω) selects a subset of configurations displaying an ordered pattern. In the case under consideration here, OBD leads to the selection of a 3-sublattice or- dered phase with Rydberg occupations of the three sublattices (n A , n B , n C ) = (λ, 1 − λ, 1/2), λ = 1/2. This phase later undergoes a quantum phase transition in the 3d-XY universality class to the paramagnetic phase for a larger Ω [54,55].

Appendix C: Finite-size effects
The two-dimensional arrays used here present open boundary conditions (OBC). The lattice sites along the boundary have less neighbors than the sites in the bulk and therefore experience less interactions. This modifies the "phase diagram" with respect to the bulk phase diagrams presented in Fig. 1. In Fig. 7 we show the classical (i.e. Ω = 0) ground-state configurations for the N = 36 square and triangular arrays used in the experiment and compare them to the bulk phase diagrams. The sites at the boundary are more easily excited to the Rydberg state when δ increases, leading to a larger number of density plateaus. We have checked that for the parameters used in this paper the OBC and bulk phases feature the same typical short-range correlations.
The Néel structure factor S Néel shows almost no finite-size effects for short-range ordered systems as long as the system sizes are larger than the correlation length. For finite arrays the number of pairs of sites with larger distances |k| + |l| is strongly reduced leading to larger errors for the connected correlations g (2) (k, l) due to reduced statistics. In order to keep reasonable values and errors for the Néel structure factor we have restricted the summation to |k| + |l| ≤ 4 in the experimental evaluation of S Néel .

Appendix D: Numerical Methods
We perform numerical exact diagonalization simulations of the time evolution on medium-size lattices for both Hamiltonian systems without any dissipation and the same systems with additional dephasing terms in the framework of a master equation in Lindblad form. In both cases we construct the complete many-body Hilbert space of the system on the finite lattice and do not perform any truncation to calculate the time evolution. Therefore, during the time evolution the system can explore the entire Hilbert space and all phases in the phase diagram can potentially be obtained.

Unitary time evolution
To simulate the unitary time evolution of a time-dependent Hamiltonian H(t), we approximate H(t) by a piecewise constant Hamiltonian H approx (t) on n steps intervals of length ∆t such that for n∆t ≤ t < (n + 1)∆t. We then compute the evolution operator within each interval with a Krylov-type matrix exponential approach [56] and use the evolved state of the previous interval as starting state. We choose n steps such that the results do not differ from a simulation with n steps /2 intervals within a demanded accuracy. In most of the presented simulations n steps = 200.
We are able to calculate the unitary time evolution of system sizes of up to about 30 lattice sites with a reasonable amount of resources. In the present study most of the numerical results have been obtained on 4×4 and 5×5 square lattices. Due to the relatively short correlation lengths observed in the experiments, the numerical results on the smaller systems can nevertheless be used to model the observed correlation functions on the larger, experimentally realizable lattices.

Time evolution in the presence of dephasing
The experimental system features several sources of imperfections [46] and the description of its evolution as a pure Hamiltonian is not exact. In particular, phase noise of the excitation lasers and Doppler shifts lead to dephasing, already for a single particle. We thus perform simulations with a phenomenological local dephasing model with a rate γ obtained by fitting experimental single atom Rabi oscillations. We obtain dephasing rates γ/U ≈ 1.1 − 1.4 for the data shown in this paper. The time evolution of the density matrix ρ(t) of the many-body system is then described by the following master equation in Lindblad form with a Liouvillian method [57][58][59] where only a single wave-function has to be stored such that one can simulate system sizes of up to around twenty sites. The MCWF method evolves a starting state with an effective non-hermitian Hamiltonian; at each time step a quantum jump corresponding to the given dissipation operator collapses the evolved state with a given probability. Averaging over many such quantum trajectories -we use 1280 for the presented results-allows reconstructing the density matrix and thus computing any observable during the time evolution. In practice, to perform simulations with dissipation, we use the Python toolbox "QuTiP" [60].

Appendix E: Lieb-Robinson bounds
In non-relativistic quantum mechanics, one might believe that information propagation is instantaneous because there is no explicit speed of light limiting the propagation. However, in 1972 Lieb and Robinson [3] proved that in an extended quantum many-body system with a finite-dimensional local Hilbert space and sufficiently local interactions, there is nevertheless a characteristic velocity emerging, which defines an approximate light-cone for the propagation of information and implementing causality.
In Ref. [4] the Lieb-Robinson bound has been generalized to equal-time connected correlation functions of two operators (normalized to unity) acting on spatial regions A and B (of size |A| and |B|) at a distance d apart from each other after some time evolution of duration t and starting from an initial state with exponentially decaying correlations with a correlation length χ. Then the connected equal-time correlation is bounded as follows: The velocity v is particularly important and is called the Lieb-Robinson velocity.
The importance of this result for the present study is that correlations at a distance d are exponentially suppressed in (d − 2vt)/χ as long as the time t < d/(2v). After that time, the correlations are bounded by a number of O(1). This time dependence is visualized in Fig. 8. In closing we note that the theorem proves rigorous theoretical bounds, however the actual dynamics does not necessarily saturate these bounds.

Appendix F: Short-time expansion
In this appendix we provide more details on the short-time expansion of the connected correlation function g (2) .

Theoretical description
In order to simplify the notation we use the nearestneighbor Hamiltonian ( = 1): Here, we keep the possibility that nearest-neighbor interactions could be different from atom to atom. We consider simple ramp shapes characterized by a constant Ω, a duration T and a linear δ(t)-dependence between δ(0) = δ 0 and δ(T ) = δ final . We are interested in the limit of very short duration T and in particular how the correlations at various distances build up as we vary the duration T . Formally, we compute the full many-body propagatorÛ (T ) which solves the time-dependent Schrödinger equation for the time-dependent Hamiltonian Eq. (F1), allowing to determine the many-body wave function at any time T via |ψ(T ) =Û (T )|ψ(0) . We then express the connected correlation function as For the linear ramps considered here, a Magnus expansion [61] of the propagator is appropriate. We rely on the leading Magnus expansion term, which can be written aŝ with For the ramp shape considered here, H avg is independent of T and is of the same form as (F1) with U ij and Ω unchanged, while δ(t) is replaced by δ avg = (δ 0 + δ final )/2. We now calculate symbolic expressions for the power series in T of the connected correlators (F2) relying on the propagator in the leading Magnus expansion form and starting from an initial state with all sites in the atomic ground state. When (U T / , ΩT, δT 1), we find the following expressions for the leading order in T for a single path on a lattice linking sites (0, 0) and (k, l) (for the sake of simplicity the results are only given for U ij = U ): • nearest neighbor correlation |k| + |l| = 1 • next nearest neighbor correlation |k| + |l| = 2 g (2) (T ) = Ω 6 T 10 2419200 • third nearest neighbor correlation |k| + |l| = 3 While the expressions look complicated, a common structure emerges: the leading T dependence for a Manhattan distance m = |k| + |l| is proportional to T 2+4m . Neglecting the actual value of the coefficients for the moment in an admittedly crude approximation, these results suggest that correlations at a next shell are suppressed by a factor T 4 compared to the previous shell, suggesting an exponential spatial decay of the connected correlations at short times with m. These results also provide an insight as to why more distant shells require longer times to develop appreciable correlations: they require larger values of T to overcome the initial high powers in T suppression of the correlations. Importantly, the neglected higher-order terms in the Magnus expansions do not alter the leading order in T for the considered ramp shapes, but only affect the sub-leading Tcoefficients. Finally, we note that at very short times, correlations induced by the direct van der Waals tail of the interactions compete with the high-order nearest-neighbor considerations here. Although strictly speaking this leads to a breakdown of the exponential suppression of correlations beyond the light-cone, the smallness of this effect prevents its observation in the present experiments.

Lattice embedding
The expressions derived above are lattice independent, as long as there is a chain of m (m = |k| + |l| above) successive nearest neighbor interactions linking the two considered sites [62]. These expressions can now be embedded in any lattice (e.g. cubic, kagome, honeycomb, ...), and the actual coefficients of the short-time series can be derived by determining the corresponding embedding coefficients. The embedding counting is illustrated in Fig. 9 for the square and the triangular lattice, yielding binomial coefficients multiplying the symbolic expressions for a single path derived above. Note that for the case U z = U w the coefficients of the single paths can already differ before taking the embedding factors into account.