## Abstract

Exotic phases of matter can emerge from strong correlations in quantum many-body systems. Quantum gas microscopy affords the opportunity to study these correlations with unprecedented detail. Here, we report site-resolved observations of antiferromagnetic correlations in a two-dimensional, Hubbard-regime optical lattice and demonstrate the ability to measure the spin-correlation function over any distance. We measure the in situ distributions of the particle density and magnetic correlations, extract thermodynamic quantities from comparisons to theory, and observe statistically significant correlations over three lattice sites. The temperatures that we reach approach the limits of available numerical simulations. The direct access to many-body physics at the single-particle level demonstrated by our results will further our understanding of how the interplay of motion and magnetism gives rise to new states of matter.

Quantum many-body systems exhibiting magnetic correlations underlie a wide variety of phenomena. High-temperature superconductivity, for example, can arise from the correlated motion of holes in an antiferromagnetic (AFM) Mott insulator (*1*, *2*). It is possible to probe physical analogs of such systems using ultracold atoms in lattices, which introduce a degree of control that is not available in conventional solid-state systems (*3*, *4*). Recent experiments exploring the Hubbard model with cold atoms are accessing temperatures where AFM correlations form but have only observed these correlations via measurements that were averages over inhomogeneous systems (*5*, *6*). The recent development of site-resolved imaging for fermionic quantum gases (*7*–*13*) provides the ability to directly detect the correlations and fluctuations present in a fermionic quantum many-body state. As demonstrated in experiments with both bosons (*14*, *15*) and fermions (*12*, *13*, *16*), microscopy gives access to the spatial variation in observables that occurs in an inhomogeneous system, yielding precise comparisons with theory. The low energy scales in cold-atom systems also allow for time-resolved observations of many-body dynamics, which typically occur on millisecond time scales. In bosonic systems, temporal and spatial resolution have been combined to observe strongly correlated quantum walks (*17*), to measure entanglement entropy (*18*), and to study the dynamics of magnetic correlations (*19*, *20*).

Here, we report trap-resolved observations of magnetic correlations in a Fermi-lattice system. Fermionic atoms in a two-dimensional (2D) optical lattice can be well described by the Hubbard Hamiltonian, a simple model in which there is a competition between the nearest-neighbor tunneling energy *t* and the on-site interaction energy *U*. Despite the apparent simplicity of the Hubbard model, it has a rich phase diagram, containing, for example, the transition from a metal to a Mott insulator. AFM spin correlations begin to appear near half-filling when the temperature scale becomes comparable to the exchange energy, which in the strongly interacting regime is . In the thermodynamic limit, what happens as the temperature is lowered further depends on the dimensionality of the system: In three dimensions, there is a finite-temperature phase transition to a state with long-range AFM order, whereas in two dimensions, such an order is prohibited by the Mermin-Wagner-Hohenberg theorem (*21*). Nonetheless, AFM correlations do arise, decaying exponentially with a correlation length that diverges as the temperature goes to zero as , where is of order unity (*22*). We use quantum gas microscopy to reveal precisely these correlations, which for our finite-size 2D system are expected to lead to long-range order at a finite temperature, where becomes comparable to the system size.

Our experiments begin with a low-temperature, 2D gas of fermionic atoms in a mixture of the two lowest hyperfine states ( and ), as described in (*12*). By adjusting a magnetic bias field in the vicinity of the Feshbach resonance located at 832G, we set the s-wave scattering length to , where is the Bohr radius (*23*). Using a 30-ms linear ramp, the atoms are adiabatically loaded into an isotropic, square optical lattice with a depth of , where the recoil energy is with Planck constant . We detect magnetic correlations by removing atoms in either spin state and measuring the resulting charge correlations with site-resolved imaging (*24*), as shown in Fig. 1. Because our imaging technique removes doubly occupied sites, both doubly occupied and unoccupied sites show up as empty sites after imaging. However, proper linear combinations of the different particle and hole correlators (measured both with and without spin-dependent removal) will account for the contribution to the signal from imperfect unity filling (*24*). From this, we determine the spin correlator (*24*)
(1)Here, , with denoting the number of particles of spin on the site at . We take an average of over all where to obtain . The nearest-neighbor, diagonal next-nearest-neighbor, straight next-nearest-neighbor, etc., correlators are given by , , and , etc. From images where neither spin was removed, we directly obtain a spatial map of the single-occupation probability, which also corresponds to the local moment .

After loading atoms into the lattice, we observe AFM correlations for nearest neighbors and diagonal next-nearest neighbors. These correlations are strongest in the cloud center, where the local chemical potential is set to approximately half-filling. The spatial maps, , and for colder (top) and hotter (bottom) temperatures are shown in panels A, B, and C, respectively, of Fig. 2. For these data, the interaction is tuned to , with (). The chemical potential is tuned to approximately in the center of the cloud for the colder data by adjusting the atom number to maximize in the center. In Fig. 2A, the suppression of in the center of the cloud is caused by the formation of doubly occupied sites and indicates that the chemical potential in the center of the cloud slightly exceeds . To heat the cloud, we hold the atoms in the optical dipole trap for 4 s before loading the lattice. After heating, the maximum detected occupation decreases from 0.89(1) to 0.84(1), with a slight broadening of the density profile, whereas the largest magnitude of the nearest-neighbor correlator decreases from 0.154(3) to 0.052(6). In this regime, where the exchange energy is much smaller than both *U* and the bandwidth, an increase in temperature quickly saturates the entropy available in the spin degree of freedom while creating little entropy in the charge degree of freedom, making the nearest-neighbor correlator much more sensitive than the density to temperature changes. For the colder data, we observe significant negative correlations in away from half-filling, which requires further theoretical investigation.

We take azimuthal averages along the equipotentials of the underlying harmonic trap to obtain and . The resulting profiles are simultaneously fit to the results of a numerical linked-cluster expansion (NLCE) of the 2D Hubbard model under a local density approximation (LDA) (*24*–*26*). From these fits, we obtain a temperature and chemical potential for the cooler data and and for the hotter data. The excellent agreement with theory provides a strong indication that the local density approximation and the assumption of thermal equilibrium are valid.

By evaporatively cooling further before lattice loading, we are able to prepare samples with even larger nearest-neighbor correlations. However, for this data, the NLCE theory error is too large away from half-filling for the fit to converge, owing to the low temperature. Because the averaged correlator in the center may not be at exactly half-filling, by comparing this value for the coldest data set to a quantum Monte Carlo (QMC) calculation at half-filling (*22*), we can determine an upper bound on the temperature. The correlator value of gives , the lowest temperature reported for a Hubbard-regime cold-atom system. The QMC calculation predicts that the nearest-neighbor correlator settles as to a value of ; our largest measured nearest-neighbor correlation is therefore 53% of the largest predicted value. In Fig. 2D, we plot our largest measured value of the correlator for samples prepared at different temperatures, the temperature derived from the NLCE fits where they converge (*x* axis), and the QMC upper bound. We find very good agreement between our data and theoretical prediction, which is consistent with half-filling at the cloud center.

We see statistically significant antiferromagnetic correlations to distances of three sites, and the sign of every measured correlator value is consistent with antiferromagnetic ordering. Our ability to measure correlations at all length scales allows us to directly extract the correlation length (Fig. 2E). Samples are prepared at three different temperatures, with the atom number optimized to achieve half-filling in the center of the cloud. Values for the correlator are taken by averaging the spatial maps over a region in the center of the cloud with a six-site radius. To determine the correlation length, we perform an exponential fit of in the center of the cloud versus *d*, where *i* = 0 (1) if *d* is such that the two sites are on the same (different) sublattice. The correlation lengths are 0.24(9), 0.39(2), and 0.51(4) sites for temperatures of 1.53(18), 0.54(7), and 0.45(2), respectively. The asterisk in Fig. 2, D and 2 shows the QMC prediction of -0.36 for the nearest-neighbor correlator at half-filling as .

Quantum gas microscopy also allows for a detailed study of the thermalization of the atomic cloud when loading into the lattice. We investigate the formation of spin correlations and the thermalization of the density distribution for different lattice loading times in Fig. 3. For these data the lattice is ramped on linearly with a varying duration . We determine the radius , where is maximized. For a cloud in thermal equilibrium with not in the center of the cloud, corresponds to half-filling (). Figure 3A shows and as a function of loading time. The detected density grows from 0.6 at very short loading times and settles at about 0.9 for ms. The loading time required for the density to settle also corresponds to the maximum absolute values for both the nearest-neighbor and diagonal next-nearest-neighbor spin correlators. The matching time scales suggest that the suppression of magnetic correlations at ms is caused by the low initial density and not by exchange dynamics. The density at short loading times is determined by the confinement of the optical dipole trap preceding the lattice loading (*12*). For loading times larger than 10 ms, both and decay, consistent with heating. The faster decay of is further indication that the spin correlators are much more sensitive than the density to temperature in this regime of parameters.

We also study thermalization by fitting the data for different loading times to the NLCE theory and performing a reduced chi-squared () analysis. Figure 3B shows versus loading time, and Fig. 3C shows individual NLCE fits at of 0.4 ms, 20 ms, and 150 ms from top to bottom. The value of settles to approximately 1 on a 20-ms time scale, which is slightly longer than the settling times for the density and spin correlator. The value of remains near unity up to our largest loading times, showing that the density and spin-correlator distributions remain consistent with thermal equilibrium.

Whereas in Bose-Hubbard systems AFM correlations appear only in the Heisenberg limit , Fermi-Hubbard systems exhibit AFM correlations at all , with a maximum in the correlations occurring near . For large , AFM correlations are suppressed because the exchange energy becomes small compared with the temperature. For , where the interaction energy is smaller than the bandwidth, charge fluctuations destroy the magnetic correlations. We study these effects by varying the scattering length for fixed Hz. In Fig. 4, we plot versus the scattering length, along with the predictions of the NLCE theory for three different temperatures. We show the calculated from Wannier functions in the lowest band; for our parameters, corrections to this single-band approximation could play a role (*27*). The data show the expected dependence on from the simple picture mentioned above. We also compare our data with theoretical isothermal curves at half-filling. In this comparison, additional factors should be considered. First, the atom number is fixed, so the chemical potential in the center of the cloud varies with . Second, we anticipate the loading entropy to be approximately fixed, as opposed to the temperature, so the data are not expected to strictly follow a single isotherm. The comparison of data with the theory reflects differences between the thermodynamic preparation of atomic and conventional solid-state systems.

Our ability to observe the in situ, site-resolved distribution of spin correlations at all distances has enabled high-precision comparison with numerical studies and detailed verification that the atomic sample behaves in a manner consistent with thermal equilibrium. These experimental benchmarks on thermal equilibrium affirm our understanding of the entropy distribution, paving the way for the implementation of entropy redistribution techniques to achieve finite-system-size long-range order (*28*, *29*). Implementation of such techniques would require precise trap-shaping protocols, which have been fruitfully demonstrated in bosonic quantum gas microscopes (*30*). Numerical simulations provide evidence that a pseudogap phase in the hole-doped 2D Hubbard model arises in conjunction with long-range AFM correlations (*31*) and should therefore be accessible in our experiment in the near future. At lower temperatures of a d-wave superconductor is expected (*32*). Further thought is required to understand how the real-space observables that we can measure might shed light on these low-temperature phases. Beyond equilibrium physics, we could also exploit our ability to take temporally resolved snapshots of the correlations in a many-body wave function, allowing for in-depth studies of nonequilibrium physics beyond the capability of existing theoretical tools (*33*).

## Supplementary Materials

www.sciencemag.org/content/353/6305/1253/suppl/DC1

Supplementary Text

Figs. S1 to S6

Table S1

Reference (*36*)

## References and Notes

**Acknowledgments:**We thank E. Khatami and M. Rigol for providing the NLCE calculations, as well as T. Paiva and N. Trivedi for the Quantum Monto Carlo calculations at half-filling. We also thank J. P. F. LeBlanc and E. Gull for providing additional data based on a dynamical cluster expansion, used for theory verification. We thank E. Demler, A. Eberlein, F. Grusdt, J. Hoffman, A. Kaufman, M. Kanász-Nagy, M. Lemeshko, L. Tarruell, L. Cheuk, M. Nichols, K. Lawrence, M. Okan, H. Zhang, and M. Zwierlein for insightful discussions. Recently, antiferromagnetic correlations have been observed in the Munich lithium quantum gas microscope and the \MIT potassium quantum gas microscope (

*34*,

*35*). We acknowledge support from the Air Force Office of Scientific Research, the Multi-University Research Initiative, and NSF. D.G. acknowledges support from the Harvard Quantum Optics Center and the Swiss National Foundation. M.F.P, A.M., and C.S.C. acknowledge support from the NSF. The authors declare no competing financial interests.