A neutral-atom Hubbard quantum simulator in the cryogenic regime

Abstract

Ultracold fermionic atoms in optical lattices offer pristine realizations of Hubbard models¹, which are fundamental to modern condensed-matter physics^2,3. Despite notable advancements^4,5,6, the accessible temperatures in these optical lattice material analogues are still too high to address many open problems^7,8,9,10. Here we demonstrate a several-fold reduction in temperature^6,11,12,13, bringing large-scale quantum simulations of the Hubbard model into an entirely new regime. This is accomplished by transforming a low-entropy product state into strongly correlated states of interest via dynamic control of the model parameters^14,15, which is extremely challenging to simulate classically¹⁰. At half-filling, the long-range antiferromagnetic order is close to saturation, leading to a temperature of $T/t=0.0{5}_{-0.05}^{+0.06}$ based on comparisons with numerically exact simulations. Doped away from half-filling, it is exceedingly challenging to realize systematically accurate and predictive numerical simulations⁹. Importantly, we are able to use quantum simulation to identify a new pathway for achieving similarly low temperatures with doping. This is confirmed by comparing short-range spin correlations to state-of-the-art, but approximate, constrained-path auxiliary-field quantum Monte Carlo simulations^16,17,18. Compared with the cuprates^2,19,20, the reported temperatures correspond to a reduction from far above to below room temperature, at which physics such as the pseudogap and stripe phases may be expected^{3,19,21,22,23,24}. Our work opens the door to quantum simulations that solve open questions in material science, develop synergies with numerical methods and theoretical studies, and lead to discoveries of new physics^8,10.

Main

The Hubbard model is a paradigmatic description of strongly correlated electrons that is central to our understanding of quantum materials³. After decades of concentrated study, it is strongly believed that this model hosts a variety of exotic quantum states^25,26, including the intriguing phases observed in the cuprate high-temperature superconductors^2,3,22. Despite its apparent simplicity, solving the two-dimensional (2D), square lattice Hubbard model has been an outstanding challenge. In recent years, building on decades of algorithmic development, combined use of the most advanced numerical methods has established certain aspects of the properties of the model^{21,23,24,27,28,29,30}. Nevertheless, a full description of the low-temperature phase diagram remains beyond the reach of theoretical tools and computational methods⁹.

Its broad applicability and numerical intractability make the Hubbard model an ideal candidate for quantum simulation^8,10. In particular, ultracold fermionic atoms in optical lattices can natively implement the requisite fermionic degrees of freedom and serve as a pristine model system for strongly correlated electrons moving and interacting in a crystalline lattice¹. The Hubbard model is described by the geometry of the lattice, a tunnelling amplitude t between neighbouring sites and a contact interaction energy U, which can all be widely programmed in cold-atom systems. However, the enlarged length and time scales in cold-atom systems lead to very small energy scales, making it challenging to reach regimes corresponding to real materials at low temperatures⁷. The lowest reported experimental temperatures of T = 0.25(2)t (refs. ^6,11,12,13) (Fig. 1a) in 2D Hubbard quantum simulators with intermediate interaction strengths of U/t ≃ 8 correspond to an effective temperature of approximately 700 K in the cuprates by normalizing temperatures with superexchange energy J = 4t²/U (refs. ^2,19,20) (Methods). One needs to reach well below room temperature to investigate the properties relevant to cuprates, including pseudogap behaviour, stripe order and unconventional superconductivity^{3,19,21,22,23,24}. New cooling schemes are, therefore, needed for cold-atom-based platforms to offer competitive simulations of quantum materials.

Fig. 1: Cooling ultracold atoms by transforming a low-entropy product state to strongly correlated states.

a, Schematic of the 2D Hubbard phase diagram. The lowest temperatures reported in refs. ^6,11,12,13 (grey dashed) and the estimated temperatures in this work (red band) are marked. b, To reach lower entropy, we prepare a gapped BI in contact with a gapless metallic reservoir. Entropy flows from the BI to the reservoir³¹. After isolating the two parts, the BI is transformed into a strongly correlated state of interest by dynamically changing the lattice geometry, Hubbard parameters and density, while preserving the very low entropy of the original BI. c, The BI has a filling of two atoms per site, which appear as empty sites in a parity-projected fluorescence image (left), in contrast to the dilute reservoir. We then halve the lattice filling by doubling the number of lattice sites, thereby making the BI visible (centre). The strongly correlated state at the end of the preparation spans about 340 sites, in which empty sites correspond to coherent doublon–hole pairs indicating low entropy (right). The central state of interest and reservoir are shaped by optical potentials programmed with DMDs (not shown). d, To double the number of sites, we continuously decrease the long-spacing lattice depth V_L and increase the short-spacing lattice depth V_S. This effectively splits doubly occupied sites into dimers, which are then connected to form a square lattice.

In this work, we demonstrate a substantial reduction in the temperatures achievable in cold-atom-based quantum simulations of the Hubbard model on a 2D square lattice. At half-filling, we achieve temperatures of $T=0.0{5}_{-0.05}^{+0.06}\,t$ with U/t ≃ 8 in a large system consisting of about 340 lattice sites, which corresponds to a several-fold improvement of the state of the art^6,11,12,13. We compare the measured spin correlations, which extend out to long range, to exact numerical results, and find excellent agreement.

The behaviour of the Hubbard model is not fully understood at finite doping, which has hindered the development of cooling techniques in this regime, even at the level of theoretical proposals. Here, we experimentally propose and demonstrate a novel pathway to realize low temperatures with finite values of doping δ between 2% and 21%. Detailed comparisons with state-of-the-art approximate numerical computations^16,17,18 indicate that our temperatures are T ≲ 0.1t, comparable to those achieved at half-filling. These temperatures correspond to a reduction to markedly below room temperature in cuprates.

Our scheme hinges on the efficient preparation of product states with extremely low entropies³¹. Seminal proposals suggested that these states can be adiabatically connected to strongly correlated states of interest to reach low temperatures^{14,15,32,33,34}. Despite experimental progress with bosonic systems^35,36,37, practical implementations of these schemes for fermions have been challenging beyond proof-of-principle realizations at half-filling^31,38,39,40. This holds especially true in large Hubbard systems, in which finite-time, out-of-equilibrium dynamics are out of reach of numerical simulations^41,42. Our work, therefore, showcases how quantum simulators can experimentally address these optimization problems.

Experimental scheme

The spirit of our protocol is similar to most cooling cycles. In a first generalized compression step, entropy is extracted from the system of interest to a reservoir in thermal contact by reducing the density of states of the system^14,31,43. After the system is isolated from the reservoir, increasing the density of states by generalized expansion decreases the temperature. In practice, our scheme involves initializing cold atoms in a low-entropy band insulator (BI), and transforming the BI into a strongly correlated state of interest^31,39 (Fig. 1b).

The BI can be prepared with extremely high fidelity because the chemical potential lies within the band gap featuring low density of states, which allows for efficient entropy redistribution to a metallic reservoir^{11,14,31,36,43}. We load a spin-balanced mixture of fermionic ⁶Li atoms in the lowest two hyperfine states into an optical lattice (long-spacing lattice). We set the magnetic bias field at 550 G close to a broad Feshbach resonance to minimize interaction strengths for BI formation. Programmable optical potentials created by two separate digital micromirror devices (DMDs) first confine about 340 atoms into a BI covering 170 sites, then isolate the BI from the reservoir. We estimate a low initial entropy per particle of s = 0.025(4)k_B (ref. ³¹), where k_B is the Boltzmann constant, based on a measured singly occupied site density of n_s = 0.7(2)% in a central disk of around 110 sites, although an accurate estimation is complicated by gradients at the edge of the confining potential due to finite optical resolution (Methods).

Keeping the subsequent transformation slow relative to the relevant many-body timescales is crucial to realizing low final temperatures. This is particularly challenging as this protocol involves a substantial change in lattice filling from two in the BI to around one atom per site in strongly correlated Hubbard systems. This massive change in the state typically involves very slow many-body timescales, which requires excellent local control of density^{31,36,38,39,40}. To address this challenge, we take advantage of an optical lattice whose geometry is dynamically adjustable^12,44,45. In particular, we continuously split single lattice sites into double wells and subsequently connect these double wells together into a square lattice, realizing a scheme similar to the one proposed in ref. ¹⁵. This doubles the number of sites at fixed total atom number (Fig. 1c,d and Extended Data Fig. 1) and naturally converts the BI into half-filled Fermi systems. However, this doubling is not exact in the experiment because of variations of the optical potentials. To precisely reach the target filling, our protocol further relies on tunable interactions and DMD potentials to ensure controlled preparation, isolation and expansion of the system of interest (Methods).

Splitting an insulator

Dynamically changing the lattice geometry leads to rich physics even at very large interactions of U/t = 18.6(8), in which the half-filled state resulting from splitting the BI is a Mott insulator with a large charge gap. Here, the Hubbard model simplifies to a spin model with antiferromagnetic Heisenberg coupling J = 4t²/U. In thermodynamic limit, this Heisenberg model is known to host a quantum phase transition between a disordered phase in the dimerized lattice, described as a product of isolated singlet states within double wells, and a Néel ordered phase in the square lattice with long-range antiferromagnetic correlations^46,47 (Fig. 2a). In our experiment, which is a finite-sized system, the many-body spin gap is expected to decrease monotonically during splitting, with a minimal value of Δ_H ∝ J/L² corresponding to the Anderson tower of states with broken symmetry⁴⁸. These two lattice geometries are continuously connected by tuning the coupling between double wells, which preserves total spin S_tot = 0 and enables efficient low-entropy state preparation¹⁵. We perform parametrization of the transformation of lattice geometry by a single parameter α describing the ratios of tunnelling amplitudes within and between double wells (Fig. 2a and Methods).

Fig. 2: Splitting a BI into a Heisenberg antiferromagnet.

a, Similar to the long-spacing lattice limit, the ground state at half-filling on a disconnected dimerized lattice is a product state of singlets. This state is then adiabatically connected to a long-range antiferromagnet on a square lattice by coupling nearest dimers. b, Spin correlations ${C}_{{\bf{d}}}^{zz}({\bf{r}})$ between nearest neighbours are localized on dimers when they are weakly coupled (α = 0.3, left) and become uniform across all nearest-neighbour bonds in the square lattice limit (α = 1.0, right). c, Spin correlations as a function of bond displacement ${C}_{{\bf{d}}}^{zz}$, averaged over an r = 6 region at half-filling at U/t = 18.6(8). The range of the antiferromagnetic correlations starts to grow around α = 0.7. d, Energy levels of the Heisenberg model on dimerized lattices as coupling strengths are tuned, in units of the spin exchange coupling J = 4t²/U (Methods). e, The measured staggered magnetization increases at couplings α > 0.7 similar to the critical point of the quantum phase transition in the ground-state Heisenberg model⁴⁷. Solid lines indicate simulation of dimerized Heisenberg model at temperature T/J = 0.5.

Source Data

To track how spin order changes with lattice geometry, we experimentally measure the spin correlation function in the z-direction between sites at positions r and r + d (ref. ¹²):

$${C}_{{\bf{d}}}^{zz}({\bf{r}})=\frac{1}{{S}^{2}}(\langle {S}_{{\bf{r}}}^{z}{S}_{{\bf{r+d}}}^{z}\rangle -\langle {S}_{{\bf{r}}}^{z}\rangle \langle {S}_{{\bf{r+d}}}^{z}\rangle ).$$

(1)

In the α → 0 limit of disconnected dimers, which corresponds to the BI, we detect only saturated spin correlations between nearest-neighbour sites, |d| = 1 within the dimers and consistent with 0 everywhere else (Fig. 2b). As these dimers are connected, spin correlations start to grow between the dimers at α ≃ 0.7 and become uniform and isotropic in the square lattice α = 1.0.

To probe spin correlations at longer bond distances |d| > 1, we obtain ${C}_{{\bf{d}}}^{zz}$ by averaging ${C}_{{\bf{d}}}^{zz}({\bf{r}})$ over a Mott insulating region within a central disc of radius r = 6, and account for spatial inversion symmetry (Methods). We confirm that spin correlations are antiferromagnetic and localized on the intra-dimer bonds in the dimerized lattice α = 0.3 and that they extend to long range in the square lattice with α = 1.0 (Fig. 2c).

The nature of magnetic order during the lattice transformation is captured by the staggered magnetization m^z (ref. ¹¹). We compute (m^z)² by averaging the sign-corrected spin correlations ${(-)}^{| {d}_{x}+{d}_{y}| }{C}_{{\bf{d}}}^{zz}$ up to a cutoff distance d_max (Methods and equation (3)). We observe an increase in (m^z)² around α ≃ 0.7 (Fig. 2e), consistent with the range of critical values reported in previous numerical studies of dimerized Heisenberg models^46,47. The experimental data are also consistent with quantum Monte Carlo simulations performed at a fixed temperature of T = 0.5J = 0.12t.

We experimentally check ramp adiabaticity and confirm that the spin correlations have saturated by varying the ramp duration, which could be limited by heating. Surprisingly, we find that the strengths of the DMD potentials confining the Mott insulator, which may lead to atom transport into and out of the central region, can strongly affect the final spin correlations, and therefore temperatures (Methods and Extended Data Fig. 3). These measurements suggest that charge transport may be crucial to temperature optimization^49,50 apart from spin dynamics of the idealized model.

Ultralow temperatures in the strongly correlated regime

At intermediate interaction strengths of U/t ≃ 8 (Fig. 3a), which are most relevant to cuprate physics, the Hubbard model remains an antiferromagnet at half-filling and can thus be reached by the same adiabatic ramp as in the Heisenberg regime. However, owing to the reduced interaction strength, the system has increased compressibility, which can lead to enhanced charge transport. As a result, we find that the harmonic confinement produced by the lattice laser intensity profile results in density variation across the system after splitting. Thus, if the central density is to remain n = 1, the edge of the system must be at n < 1. As the entirety of the system began as a BI, this means atoms at the edge of the trap must be spilled out as the insulating gap is reduced during splitting. We achieve this by dynamically adjusting the strength of the DMD potentials (Fig. 3b). During this step, we empirically find that it is helpful to allow for transport to occur at low interaction strengths, which could be related to the fact that a Fermi liquid is more compressible than a Mott insulator^41,42,50 (Fig. 4a). The timing and shape of the above ramps are extensively optimized based on the final temperature.

Fig. 3: Cold Hubbard antiferromagnets in the strongly correlated regime.

a, Intermediate interaction strengths, here U/t = 8.3(2), lead to density fluctuations apart from spin fluctuations at half-filling. b, To achieve the target filling (here n = 1), we adjust the confining potential created by DMDs to spill excess atoms. A repulsive wall isolates the cold region from the reservoir, whereas a volcano-shaped potential enables a spilling of excess atoms from the centre. c, Spin correlations ${C}_{{\bf{d}}}^{zz}$ as a function of bond displacement d, and averaged over a radius r = 3 region, show nearly saturated long-range antiferromagnetic correlations. The error bars denote 1 s.e.m. d, Experimental data are compatible with numerically exact simulations (lines, DQMC and AFQMC) at temperatures T/t < 0.1 and interactions U/t = 8. e, Low temperatures yield a sharp peak of the spin structure factor S(q) at quasimomentum q = (π, π). f, We estimate temperature by computing the staggered magnetization with cutoff d = 6 on bond distance, for both DQMC–AFQMC data (interpolated to U/t = 8.3(2), line) and experimental data (1σ confidence interval shown as shaded area). Their comparison yields a temperature $T/t=0.0{5}_{-0.05}^{+0.06}$ (Methods).

Source Data

Fig. 4: Hole-doping cold Fermi Hubbard systems.

a,b, To coherently introduce hole dopants into the half-filled state after splitting the BI, we increase the kinetic energy and decrease the interaction strength by reducing the lattice depths out of the tight-binding regime. This prepares cold weakly interacting Fermi liquid that inherits the low initial entropy, which is adiabatically expanded to reach a given target density. Reloading the lattices and ramping up interactions create cold, hole-doped Hubbard systems at U/t = 8.0(3) (III). Schemes of splitting the BI into a half-filled antiferromagnetic Mott insulator at strong interaction of U/t = 18.6(8) (I) and intermediate interaction of U/t = 8.3(2) (II) are also marked. c, Doping scheme in quasimomentum space. Splitting the BI, which completely fills the Brillouin zone, into the square lattice doubles the number of sites and the size of the Brillouin zone. The population of quasimomentum states remain nearly unchanged. Expanding in real space then decreases the size of the Fermi surface and coherently introduces doping. d, Spin correlations as a function of bond displacements ${C}_{{\bf{d}}}^{zz}$ measured in an ROI of r = 3. The range and magnitude of the antiferromagnetic correlation decreases with increased doping. e, Azimuthal average of the sign-corrected spin correlations shown in d. At small dopings of δ ≤ 10%, the antiferromagnetism still remains long-ranged over the ROI. As we increase the number of hole dopants, the strengths of spin correlations at all distances are reduced by similar amounts. f, Short-range spin correlations ${C}_{d}^{zz}$ at different distances $d=1,\sqrt{2},2,\sqrt{5}$ as we dope the system. The nearest-neighbour correlations at d = 1 show quantitative agreement between experimental data and CP-AFQMC simulations at U/t = 8 and T/t ∈ [0, 0.1] (solid line). At longer bond distances d > 1, experimental data show stronger correlations than CP-AFQMC simulations. The error bars denote 1 s.e.m.

Source Data

With the optimized experimental sequence, we observe antiferromagnetic spin correlations ${C}_{{\bf{d}}}^{zz}$ extending across the entire half-filling region with a radius r = 5 for an interaction strength of U/t = 8.3(2) (Fig. 3c and Methods). We compare these experimentally measured correlations after sign-correction and azimuthal average ${(-)}^{| {d}_{x}+{d}_{y}| }{C}_{d}^{zz}$ with the numerical results from simulations using finite-temperature determinant quantum Monte Carlo (DQMC) and ground-state auxiliary-field quantum Monte Carlo (AFQMC) methods, under open boundary conditions. At half-filling, the numerical results are in principle exact, although care must be taken to obtain accurate and unbiased spin correlations, especially at larger bond distances d (Methods). We find the measured ${(-)}^{| {d}_{x}+{d}_{y}| }{C}_{d}^{zz}$ are close to saturation and consistent with numerical data at T/t ∈ [0, 0.1] (Fig. 3d). These strong, long-ranged correlations translate to a narrow peak in the spin structure factor S(q) at quasimomentum q = (π, π) (Fig. 3e), which is obtained as the Fourier transformation of ${C}_{{\bf{d}}}^{zz}$ (Methods). Comparing numerical DQMC data with the experimental value of the squared staggered magnetization (m^z)² computed up to a cutoff d_max = 6 allows us to extract a temperature of $T/t=0.0{5}_{-0.05}^{+0.06}$ (Fig. 3f and Extended Data Fig. 7).

Doping cold antiferromagnets

The most intriguing and poorly understood physics of the Hubbard model resides in the doped regime. At low temperatures of T/t < 0.2, and with intermediate interactions of U/t ≃ 8, the doped Hubbard model is extremely challenging. Unlike for half-filling, no numerically exact results for spin correlations ${C}_{d}^{zz}$ are available in this regime for sufficiently large system sizes, which have systematically explored parameters such as temperature and interaction. Extending our cooling scheme to the above regime is, therefore, very valuable to advance our understanding, as quantum simulation can offer direct measurements, potentially adding a powerful new element in conjunction with approximate numerical simulations and analytic theory.

An outstanding challenge is that dopants must be coherently introduced and delocalized into the nearly half-filled system after splitting to avoid detrimental heating^31,49,50. To address this challenge, we must adapt our cooling scheme to reach a final state with a central density n < 1. We begin with a BI and then double the density of lattice sites, resulting in a central half-filled region. Additional expansion of atoms out of this region is then required to achieve n < 1. In previous work, however, this expansion incurred substantial heating³¹, resulting in final temperatures of T/t > 0.4. In this work, we circumvent these challenges by performing the expansion in a shallow lattice, which increases the tunnelling amplitudes and reduces interaction strengths, facilitating the efficient transport of particles. At the same time, this decreases dissipative effects from the lattice light, reducing the background heating rate.

Our cooling scheme for doped systems is thus as follows (Fig. 4a,b): after forming the BI, we reduce the lattice depth by 80% and perform splitting. The large tunnellings and weak interactions transform the BI into a Fermi liquid. We then expand the atoms by reducing the strength of the confining DMD potential, while ramping the magnetic bias field to its final value of 590 G. The depth of the square lattice is then ramped up to its final value, reaching the strongly correlated regime of U/t = 8.0(3) (Methods). At the same time, the strengths of the DMD potentials are adjusted to minimize transport during lattice loading^49,50 (Methods).

In Fig. 4d, we report the spin correlations ${C}_{{\bf{d}}}^{zz}$ averaged over a region of interest (ROI) of radius r = 3. We observe antiferromagnetic spin correlations extending over this entire region for hole dopings of up to δ = 10%. To understand how spin correlations evolve with increasing doping, we plot the azimuthal average of the sign-corrected spin correlations ${(-)}^{| {d}_{x}+{d}_{y}| }{C}_{d}^{zz}$. Larger doping induces stronger relative suppression of spin correlations at longer bond distances, which reduces the range of antiferromagnetism. The absolute magnitude of correlation reduction seems insensitive to bond distances, which suggests that doping decreases a long-range offset instead of affecting the short-range structure (Fig. 4e).

We compare the experimentally measured spin correlations systematically with the numerical results obtained with the constrained-path auxiliary-field quantum Monte Carlo (CP-AFQMC) method^16,17,18, which is a state-of-the-art approximate numerical technique. The spin correlations predicted by CP-AFQMC grow monotonically as temperature decreases. For a hole-doping range of δ ∈ [2%, 21%], we find good agreement between the experimentally measured nearest-neighbour spin correlations ${C}_{1}^{zz}$ and CP-AFQMC simulations at temperatures of T/t ∈ [0, 0.1]. Simulations at T/t = 0.15 exhibit significantly weaker correlations (Fig. 4f).

For bond distances of d > 1, although the agreement is good on the scale of the variation between modestly low temperatures (T/t = 0.25) and the ground state, the experiment shows stronger correlations than predicted by CP-AFQMC. It is not quite at the level of the agreement found at half-filling between experimental data and numerically exact results for all bond distances. We verify that our observations do not depend on the weak spatial variation of doping due to harmonic confinement by changing the radius r of the experimental ROI and further confirm that the system is at thermal equilibrium (Extended Data Fig. 3). Uncertainty in the calibration of experimental parameters (Methods) is insufficient to explain the deviation, but other details in Hamiltonian difference (for example, higher-order corrections) could potentially contribute, especially at lower temperatures. The discrepancy may also hint at an underestimation of long-range correlations by CP-AFQMC. This would be consistent with data at elevated temperatures of T/t = 0.25, 0.33, at which numerically exact DQMC results seem to be more consistent with experimentally measured spin correlations, and CP-AFQMC results show slightly lower correlations (Methods and Extended Data Fig. 5). These systematic comparisons are an example of the synergies that now exist between quantum simulations and numerics in the regime in which obtaining exact numerical results are very challenging.

Previous attempts³¹ at adiabatic state engineering in the Hubbard model suffered from excessive heating during the transformation from BI to Hubbard systems. Our scheme departs from these attempts in two ways: the use of a geometry tunable optical lattice enables us to halve the filling of the system without macroscopic density redistribution and working at weak interactions in a shallow optical lattice minimizes heating effects during the remaining expansion. To probe the the importance of the latter choice, we prepare a BI in a square lattice and expand it to a half-filled system in the same lattice by shaping the DMD potentials and measure the temperature of the final state as a function of the expansion duration τ and the normalized lattice depth η during expansion (η = 1 indicates full depth).

To characterize the adiabaticity of the expansion, we scan the expansion duration at fixed depth η = 0.2 (Fig. 5a). We observe a rapid increase in the temperature at short expansion durations, indicating τ = 20 ms ≃ 35ħ/t as a critical timescale below which diabatic heating limits the temperature. As the lattice depth is increased, the tunnelling energy decreases, requiring longer expansion durations to remain adiabatic. To probe the effects of lattice depth, we, therefore, measure the final temperature as a function of lattice depth at a fixed normalized expansion duration τ ≃ 60ħ/t, in which t is computed from η using a bandstructure calculation (Methods). As shown in Fig. 5b, the final temperature rapidly rises when the lattice depth is increased beyond a critical value η ≃ 0.6, which cannot be because of purely single-particle adiabaticity effects. Rather, we attribute it to a combination of interaction-induced heating during particle transport, which worsens as U/t is increased, and dissipative effects induced by the lattice light and DMD light, which worsen because of the increase in η and τ.

Fig. 5: Adiabaticity of expansion.

a, Adiabaticity of expansion from a BI to a to half-filled Hubbard system in the same lattice as a function of expansion duration τ at normalized lattice depth η = 0.2. We find the ramp adiabaticity and, therefore, temperature improves quickly as the τ is increased from 1.76ħt (1 ms), which starts to saturate at τ ≃ 35ħ/t (20 ms). b, Spin correlations and estimated temperatures after expansion of fermions for a fixed normalized duration measured by tunnelling times τ ≃ 60ħ/t (taking ħ = 1) at different lattice depths η. Expansion with η ≥ 0.8 show strong heating that cannot be explained exclusively by adiabaticity. At critical lattice depth η ≃ 0.6, the interaction strength U/t ≃ 1.5 is similar to the interaction at which a rapid increase in transport timescales was observed in ref. ⁵². This indicates the heating may be largely dependent on a combination of interaction-induced effects and lattice heating. The error bars denote 1 s.e.m.

Source Data

The interaction-induced heating may be related to theoretical and experimental studies that indicate that transport in strongly correlated systems can experience an exponential slowdown^42,51,52 and induce heating from density redistribution^49,50. We note that the interaction strength U/t ≃ 1.5 at the critical lattice depth indicated in Fig. 5b is similar to the interaction at which a marked increase in transport timescales was observed in ref. ⁵². These measurements highlight the importance of performing the expansion step in a shallow lattice, apart from the use of a superlattice to double the density of sites.

Outlook

In this work, we have demonstrated a substantial reduction in the temperatures achievable in large-scale quantum simulations of the Hubbard model, both with and without doping. Through careful dynamical tuning of parameters in a programmable optical lattice, we are able to transform a trivial low-entropy state into a cold, strongly correlated quantum many-body state in equilibrium. Moreover, the dynamics involved in this transformation are exceedingly difficult to tackle by analytical techniques, or to simulate classically, and so must be optimized empirically using the quantum simulator itself. The above approach is broadly applicable to situations in which one can prepare low-entropy product states with high fidelity, for example, with optical tweezers⁵³. These techniques can also be extended to other lattice geometries, including triangular and kagome lattices, which may host quantum spin liquids²⁶, and square lattices with diagonal tunnelling t′, which may host unconventional superconductivity⁵⁴.

The temperatures and system sizes achieved in this work probably allow one to enter phases of the Hubbard model that have not yet been explored in cold-atom quantum simulators, including phases involving charge order^3,6,21,22,23. Recent advancements involving comparisons of different approximate numerical simulations in equilibrium^27,28 have shed some light on these behaviours, including the interplay between the stripe and pseudogap phases and d-wave superconductivity^24,54. However, a complete microscopic understanding has remained unknown. Quantum simulations in the appropriate regimes offer a unique opportunity to study the features that are challenging to simulate classically, including spectral properties and real-time dynamics. Furthermore, in equilibrium, quantum simulations provide a valuable extra data point that complements approximate numerical simulations by taking different assumptions and approximations.

Although advanced numerical simulations have long informed the design and operation of cold-atom-based quantum simulators, this work opens the possibility of a fruitful exchange, in which the results obtained by quantum simulation can also be used to develop and benchmark more efficient and accurate numerical techniques. The benefits of combining classical and quantum tools go beyond simply serving as mutual benchmarks. The success of the preparation scheme in this work suggests that, in the future, hybrid classical-quantum algorithms that use the results obtained from a quantum simulator to optimize state preparation^55,56 could provide substantial advantages, leading to lower temperatures, and possibly a definitive understanding of the Hubbard model in and out of equilibrium.

Methods

Lattice potential

In this work, the lattice potential is formed by three retro-reflected laser beams. Two beams are overlapped and mode-matched, propagating in the x-direction, and are referred to as the X and $\bar{X}$ beams, respectively. The third beam is propagating along the y-direction, orthogonal to the other two, and is referred to as the Y beam. As described in previous works^12,57, the X and Y beams are phase stabilized to form an interfering lattice. The frequency of the $\bar{X}$ beam is detuned by about 1.6 GHz from X and Y (Extended Data Fig. 1d,e). The frequency offset is converted to a lattice phase shift of π upon reflection from the retro-reflecting mirror, shifting the $\bar{X}$ lattice by half a site relative to the X lattice. This allows the $\bar{X}+Y$ lattice to split each unit cell in the X + Y lattice symmetrically into two sites. All three beams are reflected from a super-polished substrate in the z-direction at an angle of θ = 69.2(1)° before being retro-reflected¹², which forms a three-dimensional lattice potential. The lattice potential in the z-direction provides confinement with a trapping frequency much larger than all relevant energy scales in the xy plane lattices, and the tunnelling along the z-direction is negligible during the experimental sequences. This allows for a near-ideal realization of a 2D system in the xy plane after the atoms are selectively loaded into a single layer of the z lattice. The potential of the 2D lattice can be written as

$$\begin{array}{l}V(x,y)\,=\,-\frac{{V}_{x}{\bar{r}}^{4}(1+\cos (2\theta ))}{4}\cos (2{k}_{x}\,x)+\frac{{V}_{\bar{x}}{\bar{r}}^{4}(1+\cos (2\theta ))}{4}\cos (2{k}_{x}\,x)\\ \,\,\,\,\,-\frac{{V}_{y}{\bar{r}}^{4}(1+\cos (2\theta ))}{4}\cos (2{k}_{y}\,y)\\ \,\,\,\,\,-{\bar{r}}^{2}{\cos }^{2}\theta \frac{\sqrt{{V}_{x}{V}_{y}}}{4}(\cos ({k}_{x}\,x-{k}_{y}\,y+\phi )+{\bar{r}}^{2}\cos ({k}_{x}\,x+{k}_{y}\,y+\phi )\\ \,\,\,\,\,+\,{\bar{r}}^{2}\cos ({k}_{x}\,x+{k}_{y}\,y-\phi )+{\bar{r}}^{4}\cos ({k}_{x}\,x-{k}_{y}\,y-\phi ))\\ \,\,\,\,-\frac{{V}_{x}{\bar{r}}^{2}(1+{\bar{r}}^{4})(1+\cos (2\theta ))}{8}-\frac{{V}_{\bar{x}}{\bar{r}}^{2}(1+{\bar{r}}^{4})(1+\cos (2\theta ))}{8}\\ \,\,\,\,\,-\frac{{V}_{y}{\bar{r}}^{2}(1+{\bar{r}}^{4})(1+\cos (2\theta ))}{8}\end{array}$$

(2)

Here, $\bar{r}=8.27(1) \%$ denotes the Fresnel loss present at the surface of the glass cell, k_x = k_y = 2π sin θ/λ are the lattice vectors in the xy plane, and λ = 1,064 nm. Owing to loss, the four interference terms cannot be combined except when the interference time phase ϕ between X and Y is 0 or π. Note that the lattice depths ${V}_{x,\bar{x},y}$ do not directly correspond to the lattice depths in an ideal retro-reflected square lattice, because of the finite incident angle in the z-direction and losses due to the presence of the glass cell.

Experimental methods

As in previous work¹², we prepare an ultracold, spin-balanced Fermi gas of ⁶Li atoms in the lowest two hyperfine states by evaporative cooling in a crossed optical dipole trap. Spin balance is achieved through a microwave mixing process⁵⁸ with a duration of 300 ms, ensuring the SU(2) symmetry of the Fermi gas. The ultracold ⁶Li atoms are then loaded from the optical dipole trap into the interfering (long-spacing) lattice. The interfering lattice is a square lattice with √2 larger spacing than the non-interfering (short-spacing) lattice and is formed by the actively phase-stabilized X and Y beams¹² using equal intensity. This initial lattice loading is performed in 100 ms, and the final lattice depths after loading are V_X = V_Y = 2.88(1)E_R. Here, E_R = h²/(8ma²) = 25.49(4) kHz denotes the recoil energy and a denotes the lattice spacing. We set the magnetic field to 550 G, resulting in an s-wave scattering length of a_s ≃ 84a₀, and ensure that the applied field settles before lattice loading. At the start of lattice loading, we turn on a blue-detuned confining potential formed by one of two DMDs, which we refer to as DMD0 and DMD1. DMD0 is turned on for 60 ms (see section ‘Characterization of BI’) and, together with the harmonic confinement provided by the lattice, controls the density profile of the atoms. This redistributes the entropy in the system, creating a BI that is in thermal contact with a dilute metallic reservoir (see section ‘Characterization of BI’). Once the loading is complete, we turn on a second blue-detuned DMD potential projected using DMD1 in 5 ms (ref. ³¹), which acts as a wall to separate the BI from the reservoir in subsequent steps in the experiment.

Experimental sequence for the half-filling data

As described in ref. ¹⁵, the phase transition from a BI to a half-filled Mott insulator with long-range antiferromagnetic correlations can be tuned by adjusting the ratio between inter- and intra-dimer tunnelling amplitudes. In the Hubbard model, the intra-dimer tunnelling t_d is proportional to the band gap Δ_g in the BI limit. The large band gap of Δ_g ≃ 100 kHz allows for a fast ramp of the $\bar{X}$ and Y beams in 2 ms, while remaining adiabatic relative to motional energy scales. The final Y and $\bar{X}$ lattice depths are 9.29(3)E_R and 6.28(2)E_R, respectively. These values satisfy the relation ${V}_{X}+{V}_{\bar{X}}\approx {V}_{Y}$, ensuring that the harmonic confinement induced by the Gaussian intensity profile of the lattice beams is approximately rotationally symmetric. Next, we ramp the depth of X to 0.19(1)E_R in 15 ms, and simultaneously increase the depth of $\bar{X}$ to 9.06(3)E_R. This compensates for the decrease of V_X, and maintains ${V}_{X}+{V}_{\bar{X}}\approx {V}_{Y}$. During this ramp, each site in the long-spacing lattice is split into a dimer, as shown in Extended Data Fig. 1d. In quasimomentum space, this corresponds to reducing the band gap Δ_g between the ground band and the first excited band, while keeping the gaps from the ground band to higher bands nearly unchanged. The ramp needs to be slower than the band gap Δ_g between the ground and first excited bands to be adiabatic. This gap decreases with decreasing X lattice depth. We separate the ramp into two segments with a duration of 3 ms before and 10 ms after depths of ${V}_{X}=0.52(1){E}_{{\rm{R}}},{V}_{\bar{X}}=8.74(3){E}_{{\rm{R}}}$. The slower ramp at lower X depths helps to accommodate the requirement of adiabaticity.

The set point of V_X is further reduced to 0.032E_R, and ${V}_{\bar{X}}$ increased to its final value of 9.22(3)E_R. Note that there is a systematic uncertainty in V_X below 0.1E_R because of residual leakage from the retro-reflected $\bar{X}$ lattice onto the X lattice intensity regulation photodiode. This causes an offset as large as 0.01E_R and decreases the actual X lattice depth. During the above ramp, the phase transition from spin singlets to antiferromagnetic order is expected to occur. The ramp speed is 50 ms, which is a balance between adiabaticity and heating from the optical lattice. Despite the presence of long-range antiferromagnetic correlations, the resulting state still exhibits dimerized spin correlations due to the imbalanced inter- and intra-dimer tunnelling along the x-direction (see section ‘Calibration of tunnelling and lattice parameters’). We further decrease the set point of the X lattice over 30 ms to 0.016E_R, which is the lowest depth that allows us to maintain phase stabilization between the X and Y lattices. The depth of V_X = 0.016E_R is not sufficient to remove the dimerization of the tunnelling, and so we further ramp the setpoint of the phase stabilization from ϕ = 0 to π/2. Note that the contribution from phase noise to intensity noise is negligible at these low values of V_X. Although this step is supposed to eliminate the interference lattice⁵⁹, loss in the retro-reflected lattice beam (see section ‘Lattice potential’) results in imperfect suppression of the interference term⁶⁰, and thus a potential offset between the A and B sublattices of the square lattice. We use numerical simulations to confirm that this sublattice offset does not alter the relevant physics at half-filling (see section ‘Lattice potential calibration’).

To carry out the experiments shown in Fig. 4 (see section ‘Experimental sequence for doped data’) on the doped Hubbard model, we upgraded the apparatus with a non-reciprocal attenuator (see section ‘Non-reciprocal attenuator’). This allows us to decrease V_X to a depth of 3.2(3) × 10⁻⁵E_R, corresponding to negligible tunnelling dimerization. With this upgrade, we can avoid ramping the phase stabilization setpoint ϕ and simply trigger the attenuator at the beginning of the dimerized lattice ramp. As a verification, we repeat the experiments at half-filling with the attenuator while keeping ϕ = 0. Under these conditions, we obtain a temperature of T/t ≃ 0.1, as shown in Extended Data Fig. 2. These measurements were performed without careful optimization and obtained over a 12-h period without realignment. The achieved temperature is consistent with those obtained without the attenuator for similar data taking durations without realignment. We therefore attribute the slight increase in temperature compared to the data shown in Fig. 3 to lattice drifts (see section ‘Effects of alignment’).

In the short-spacing square lattice, the harmonic confinement imposed by the Gaussian intensity envelope is not compensated. As a result, the final density profile is not flat except in a radius r = 5 region near the lattice centre at half-filling. If all of the atoms in the initial BI are kept inside the DMD0 confining potential, the harmonic confinement will increase the chemical potential μ. According to DQMC simulations at the Hubbard interaction studied in most parts of this work of U/t ≃ 8, |μ| > t will lead to non-negligible doping δ > 1%. Therefore, harmonic confinement may introduce particle dopants into the system. If the DMD0 potential had arbitrarily fine spatial resolution, allowing for an infinitely sharp wall, then as long as the DMD0 potential V_D remains smaller than the band gap and V_D ≫ U, the system would have a well-defined open boundary. However, the resolution of the DMD0 potential is limited by the numerical aperture (NA = 0.7) of the microscope objective, and the wavelength of 650 nm used to create the DMD potentials. The excess atoms residing on the finite slope of the DMD0 potential will further increase the chemical potential μ of the central system as a function of the strength of the DMD0 potential, causing additional particle doping. This doping leads to charge transport through the system at strong Hubbard interaction, which we suspect leads to substantial heating. To minimize charge transport and heating, as described in the main text, the potential strength of DMD0 is scanned and optimized for each ramp. We find the lowest temperatures are realized when the DMD0 potential is set to maintain μ ≃ 0 in the lattice centre, such that no transport occurs at this location throughout the entire ramp. Close to the boundary of the DMD0 potential, the decreased local chemical potential causes excess atoms to flow out of the central region and form a secondary reservoir. This could be beneficial for reducing the final temperature by realizing a secondary entropy redistribution step. We start to ramp the magnetic field to its final value of 620 G after the dimerized lattice has mostly been formed, and where the density distribution matches that of the short-spacing square lattice. We find this is important to reach low temperatures, which we ascribe to interaction-induced heating in transport at large U/t and to interaction-driven changes in the density profile.

The tunnelling amplitudes for all datasets are summarized in Extended Data Table 1. Once the lattice and DMD0 ramp finishes, we immediately quench the lattice potential to ${V}_{\bar{X},Y} > 60{E}_{{\rm{R}}}$ within 50 μs to freeze the dynamics before imaging. We confirmed that this quench is faster than doublon–hole dynamics and slower than band excitations. The probed singlon density starts to increase once the quench time is slower than 100 μs because of the combination of virtual doublon–hole excitations and starts to decrease once the quench time is faster than 20 μs because of excitation to higher bands.

Experimental sequence for doped data

The transport described in section ‘Experimental sequence for the half-filling data’ is necessary only because of imperfections of the potential, including finite optical resolution and harmonic confinement. Splitting a BI in the long-spacing lattice naturally gives a half-filled state in the short-spacing lattice with μ ≃ 0. To reach a final state with finite doping, coherent redistribution of atom density across the entire system is required, which poses more challenges to suppress heating during atom transport.

After isolating the BI, we first reduce the lattice depth to V_X = V_Y = 0.50(1)E_R in 30 ms, which increases the kinetic energy scale and further reduces the interaction strength before splitting into a Fermi liquid. This depth corresponds to around 20% of the original lattice depth and is chosen to maximize the tunnelling energy in the lowest band to enhance adiabaticity and reduce lattice-induced heating, while still maintaining a bandwidth smaller than the height of the DMD0 potential. The latter constraint is important because otherwise the atoms would have enough kinetic energy to move out of the confining potential, leading to uncontrolled transport. Another heuristic intuition is to maintain a finite band gap such that the Fermi surface before and after splitting are approximately matched. In the non-interacting limit, the quasimomenta of the atoms are conserved in a homogeneous lattice, and deformation of the occupation in quasimomentum space leads to diabatic transfer of atoms to excited states. In an ideal tight-binding model with only nearest-neighbour tunnelling, the full first Brillouin zone of the long-spacing lattice matches the Fermi surface of the half-filled first Brillouin zone of the short-spacing lattice (Fig. 4c). The approximate match between the occupied states before and after splitting may minimize the deformation of the Fermi surface, and therefore the heating process above. A finite band gap also suppresses diabatic excitation into higher bands. The reduced lattice depths before splitting weaken the confinement in the z-direction, potentially leading to tunnelling in the z-direction. To avoid this, we turn on a vertical lattice to allow the system to always be treated as 2D.

Unlike at half-filling, where the Mott insulating state is robust against potential variations, doped Hubbard systems are compressible and sensitive to potential offsets. We, therefore, upgraded the lattices with a non-reciprocal attenuator that allows us to almost completely turn off the X lattice, while maintaining phase stabilization between the X and Y lattices at ϕ = 0 (see section ‘Non-reciprocal attenuator’). An interference phase setpoint of ϕ = 0 ensures no potential offset, and therefore no density offset between the A and B sublattices. We ramp on $\bar{X}$ to 0.50(1)E_R, and ramp X to 3.2(3) × 10⁻⁵E_R in 30 ms. Owing to its low value, the X depth is directly calibrated using a power meter (PM100D and S121C from Thorlabs). We trigger the attenuator at the beginning of this ramp to maintain phase stabilization while ramping X. At the end of this ramp, the lattice has negligible tunnelling dimerization. At this stage, the splitting sequence is complete, and the first Brillouin zone doubled.

Expansion of the atom cloud is achieved by lowering the DMD0 potential in 50 ms, accompanied by a ramp of the magnetic field from 550 G to 590 G. This expansion time is chosen to allow the magnetic field to settle. Decreasing the DMD0 potential allows atoms to flow out of the confining potential, forming a dilute secondary reservoir confined by the DMD1 wall. The DMD0 potential is chosen to adjust the central density to the desired value in the final state, which will minimize transport during lattice reloading. As expansion corresponds to a reduction of the Fermi momentum, a change in quasimomenta is necessary. As a result, it is possible that increased scattering lengths due to the higher magnetic field may facilitate thermalization during expansion.

The prepared state is a weakly interacting Fermi gas in equilibrium, which is similar to the state used to load the lattice in previous works¹¹, but at a much lower temperature. We then load the non-interfering lattice by ramping the $\bar{X}$ and Y lattice depths up to ${V}_{\bar{X}}=11.0(1){E}_{R},$${V}_{Y}=11.0(1){E}_{R}$. Here, the lattice depths are deeper than in the half-filling protocol to allow us to reach an interaction strength of U/t ≃ 8 at a lower magnetic field of 590 G (see section ‘Choice of lattice depths and magnetic field’) to maintain the atom density profile, and minimize transport in the regime in which U/t is large, we lower the DMD0 potential and the DMD1 wall to accommodate reduced tunnelling energies at higher lattice depths during the above ramp. After lattice reloading, we freeze the dynamics within 50 μs before imaging.

Potentials of DMD0

We use an incoherent light source at 650 nm to illuminate the DMDs, which is blue detuned of the D1 and D2 transitions in Li, and forms a repulsive potential. For experiments at half-filling, the volcano-shaped potential created by DMD0 is composed of a paraboloid, which is chosen to compensate the harmonic confinement imposed by the lattice beams^11,31, and a circular central cutout region that we will refer to as the crater (Extended Data Fig. 1). The transition from the paraboloid to the crater is sharp (within 1 pixel) on DMD0, which induces diffraction fringes on the potential due to the finite resolution of the imaging system.

To allow expansion for doped systems, we create the flattened volcano-shaped potential (see section ‘Adiabaticity of expansion after splitting’) by applying a cutoff on the maximum amplitude of the volcano potential (Extended Data Fig. 1g). This forms a flat ring-shaped region surrounding the crater. As described in section ‘Characterization of BI’, the flattened volcano-shaped potential may result in worse BI fidelity than without flattening. We choose to work with the largest flattened region for which we do not detect a reduction in BI fidelity.

Non-reciprocal attenuator

As described in section ‘Experimental sequence for doped data’, even when ϕ = 0, the intensity of the X lattice needs to be reduced by more than that in the half-filled case to avoid sublattice offsets. This is challenging because of the reciprocal nature of most optical attenuators, including acousto-optical modulators or ND filters, which act on the beam both in the forward direction (towards the atoms) and in the reverse direction (on retro-reflection, returning to the detection photodiode). A 50-dB attenuation of the X lattice beam, therefore, leads to a 100-dB attenuation of the laser intensity on the phase stabilization photodiode, resulting in greatly decreased gain and signal-to-noise ratio, making it impossible to perform effective stabilization. We solve this problem by introducing a non-reciprocal attenuator⁶⁰, which applies variable and differing levels of attenuation to the forward and reverse beams. This allows us to reduce the X lattice depth by five orders of magnitude, while keeping the interference phase actively stabilized. Under these conditions, we detect no tunnelling dimerization or potential offset (see section ‘Experimental verification’).

Effects of alignment

We find that the main source of instability in the experiment is the drift of the lattice positions (especially of the $\bar{X}$ lattice relative to the X lattice). This results in drifts in the lattice harmonic confinement and corresponding shifts in the position of the peak atom density in the cloud. The lattice position is affected by temperature and humidity fluctuations in the lab and remains stable for about 0.5 h, which is the typical duration of a contiguous scan. We observe drifts in lattice position on longer timescales, which may lead to heating due to excess transport during the lattice ramps. Empirically, we find the strongest effects of heating due to lattice drifts at half-filling.

For the doped data presented in Fig. 4, we average over ROIs at the centre of the trap with different radii, in which the effects of harmonic confinement are minimized. However, the lattice drifts shift the atom distribution relative to the DMD0 potential, which introduces a potential gradient within the ROI. Therefore, for data presented in Figs. 3 and 4, we re-centre the lattice beam positions about every 1 h by maximizing the light back-coupled into the lattice optical fibres.

The drifts of the DMD potentials along the x- and y-directions are strongly suppressed by the magnification of the high NA objective used to project them. We find that focus drifts along the z-direction occur over the course of weeks, which is convenient to correct.

Choice of lattice depths and magnetic field

The Hamiltonian of ultracold fermionic atoms moving and interacting in optical lattices naturally realizes the Hubbard model, with corrections in the form of beyond-nearest-neighbour tunnelling, density-assisted tunnelling, off-site interactions and higher bands effects⁶¹. Shallow lattices are preferred from an experimental perspective because, at a fixed target value of U/t, larger tunnelling energies make all energy scales higher. Harmonic confinement is reduced in shallower lattices too, which, combined with the increased tunnelling energies, results in a more homogeneous density distribution. However, once the lattices are too shallow, longer-range tunnelling grows appreciably, and the band gap decreases. These effects lead to deviations from the Hubbard Hamiltonian and a breakdown of the tight-binding and single-band approximation. As lattice depth increases the band gap increases, the Wannier function are more localized, a smaller scattering length a_s is needed to achieve the same value of U/t, and the above corrections are exponentially suppressed compared with nearest-neighbour tunnelling t and on-site interaction U. Deep lattices are, therefore, preferred from a theoretical perspective, and a tradeoff between experimental performance and asymptotic realization of the Hubbard model in an exact manner needs to be made. In the non-interfering (short-spacing) lattice, in the half-filled case with ${V}_{\bar{X}}=9.22(3){E}_{{\rm{R}}},{V}_{Y}=9.29(3){E}_{{\rm{R}}}$, which yields a radial band gap of Δ_xy,g ≃ 83 kHz and a vertical band gap Δ_z,g ≃ 42 kHz, we find the next-nearest-neighbour d = (2, 0), (0, 2) tunnelling to be t″ ≃ 0.042t. In the doped case with ${V}_{\bar{X}}={V}_{Y}=11.0(1){E}_{{\rm{R}}}$, which yields a radial band gap of Δ_xy,g ≃ 93 kHz and a vertical band gap Δ_z,g ≃ 47 kHz, we find the next-nearest-neighbour tunnelling to be t″ ≃ 0.029t. Note that the radial lattice is nearly separable in the x- and y-directions, resulting in vanishing tunnelling along directions other than x or y. Therefore, for the experimental parameters in this work, the single-band, 2D and tight-binding approximations are well satisfied. Note that the lattice depths provided here take Fresnel loss and the angle of polarization in the apparatus into account. The quoted depths are, therefore, higher than those of an idealized retro-reflected square lattice with the same band properties.

In practice, we find that the maximum s-wave scattering length a_s available to achieve the targeted Hubbard parameter U/t sets the limit on the magnitude of tunnelling amplitudes. This is because shallower lattices yield larger tunnellings and smaller integrals of the Wannier functions, which requires increased a_s to achieve a targeted U/t. In the vicinity of the Feshbach resonance, the universal scaling of the fermion three-body recombination rate is $\kappa \propto {a}_{{\rm{s}}}^{6}$ (ref. ⁶²). Increasing the s-wave scattering length will eventually lead to excessive three-body loss and thus heating. We find that if the atoms are kept in a lattice with a depth of about 10E_R, as applies to the half-filling data, a_s = 512a₀ at 620 G is the highest scattering length that does not lead to noticeable excess heating. To achieve U/t = 8, we set the lattice depths to 9.2E_R. To obtain the doped data, the lattices are ramped down to 0.5E_R for expansion. Here, the lack of protection against three-body recombination from the lattice leads to significant heating at 620 G magnetic field. We, therefore, choose to work at a_s = 294.5a₀ at 590 G and set the final lattice depth after reloading to 11E_R. Details of the lattice parameters are listed in Extended Data Table 1.

Imaging procedure and fidelities

We perform site-resolved fluorescence imaging in the short-spacing square lattice as described in ref. ⁶³. The fidelity of correctly determining the occupation of a lattice site is F_i = 99.4(6)%.

Imaging in the long-spacing lattice differs from imaging in the short-spacing lattice and is described in ref. ¹². To image the BI with parity projection, we set the frequency detuning between $\bar{X}$ and X to 850 MHz using a radiofrequency synthesizer, which ensures good overlap with the X–Y and imaging lattices. This allows for deterministic atom transfer between the physics and imaging lattices, despite the fact that the imaging lattice contains twice as many sites as the physics lattice. As described in ref. ⁶⁰, doublons are converted into molecules by ramping through a narrow Feshbach resonance at 543 G and are subsequently lost because of light-assisted collisions. We report a combined imaging fidelity, including detection fidelity and physics-imaging transfer fidelity, of 99% (ref. ⁵⁷).

Imaging with full charge resolution in the long-spacing and dimer lattices is described in ref. ⁵⁷. The dimer lattice is adiabatically connected to the long-spacing lattice before the gap between the ground band and the first excited band closes. We set the frequency detuning of $\bar{X}$ relative to X to 1,552 MHz, which sets the position of the potential minimum of the $\bar{X}+Y$ lattice to be symmetric with respect to each unit cell in the long-spacing lattice formed by X + Y, and in the dimer lattice formed by $X+\bar{X}+Y$. Each site in the long-spacing and dimer lattice is symmetrically split into two, with a negligible potential offset between the two minima of $\bar{X}+Y$. Moreover, we ramp the magnetic field to 610 G to generate strong on-site interactions between the atoms on doubly occupied sites. This facilitates adiabatic transfer of doubly occupied sites in the long-spacing and dimer lattice to singly occupied sites in the $\bar{X}+Y$ lattice. We find the doublon detection fidelity to be 98% after image reconstruction.

Data analysis

In Figs. 2 (excluding Fig. 2b), 3 and 4, two-point spin correlation functions are spatially averaged over all pairs of sites within a circular ROI centred on the atomic cloud. Atomic density decreases away from the centre as a result of the confining potential imposed by the lattice beams and the DMDs. A large potential gradient would enhance the effective superexchange interaction J for a site-to-site potential offset of Δ_V < U and suppress magnetic interactions with a kinetic origin⁵⁷. We, therefore, chose an ROI radius of r = 6 sites in Fig. 2, r = 5 sites in Fig. 3, and r = 3, 4 and 5 sites in Fig. 4, to limit the potential variation. This limits the variation of the radially averaged singlon density to about 2% within the ROI. The correlation maps for r = 4 and 5 of the doped data are shown in Extended Data Fig. 3a,b. For large dopings, for which only short-range antiferromagnetic correlations are present, the sign-corrected, radially averaged spin correlations ${(-)}^{| {d}_{x}+{d}_{y}| }{C}_{d}^{zz}$ may become negative at certain long bond distances. This could be induced by harmonic confinement or residual disorder of the lattice potentials and will be explored in future works.

In bipartite dimerized lattices (Fig. 2), we separately average the spin correlation function ${\langle {S}_{{\bf{r}}}^{z}{S}_{{\bf{r}}+{\bf{d}}}^{z}\rangle }_{{\mathcal{S}}}$ over reference sites ${\bf{r}}\in {\mathcal{S}}$ belonging to one of the two sublattices ${\mathcal{S}}={\mathcal{A}},{\mathcal{B}}$ and for all bond vectors d. Assuming inversion symmetry, we obtain and plot the sublattice-averaged spin correlation function $\langle {S}_{{\bf{r}}}^{z}{S}^{z}{\bf{r}}+{\bf{d}}\rangle =({\langle {S}_{{\bf{r}}}^{z}{S}_{{\bf{r}}+{\bf{d}}}^{z}\rangle }_{{\mathcal{A}}}+{\langle {S}_{{\bf{r}}}^{z}{S}_{{\bf{r}}-{\bf{d}}}^{z}\rangle }_{{\mathcal{B}}})/2$. All error bars indicate the 1σ confidence interval obtained using bootstrap sampling across all experimental snapshots of a given dataset with 500 randomly selected samples. The number of experimental realizations for each dataset is reported in Extended Data Table 1.

To convert the experimentally measured singlon density ${n}_{{\rm{s}},det}$ to doping δ in Fig. 4, we first correct for imaging fidelity to extract ${n}_{{\rm{s}},{\rm{corrected}}}={n}_{{\rm{s}},\det }/{F}_{i}$. We then use the doublon density n_d as a function of density n obtained from CP-AFQMC simulations at T/t = 0 to reconstruct the doublon density as a function of singlon density n_d(n_s) = n_d(n − 2n_d). Finally, we compute the density ${n}_{\exp }={n}_{{\rm{s}},{\rm{corrected}}}\,+$ $2{n}_{{\rm{d}}}({n}_{{\rm{s}},{\rm{corrected}}})$ and doping ${\delta }_{\exp }=1-{n}_{\exp }$ of the experimental data.

In Extended Data Fig. 4b, we plot the temperature dependence of the doublon density up to n = 0.85. Similar to half-filling, the temperature dependence is negligible compared with the statistical errors in our detected densities for T/t < 0.15. Although in these CP-AFQMC calculations, the n_d results are not as accurate at higher densities of n ∈ (0.85, 1), interpolating between n = 0.85 and half-filling suggests that the temperature dependence of the doublon density is still negligible for T/t < 0.15. This is supported by the fact that the energy scale associated with doublons is the interaction strength of U ≫ t ≫ T. Therefore, using the ground-state doublon density leads to negligible systematic errors in the reported quantities.

Discrepancy between experimentally measured spin correlations and CP-AFQMC results

To better understand the discrepancy between the experimental measurements of beyond-nearest-neighbour spin correlations in the presence of doping and CP-AFQMC simulations in Fig. 4, we explore a few possible sources of similar deviations.

As described in section ‘Calibration of Hubbard interaction U’, we calibrate U/t using the singlon density at half-filling, where n_s reaches its maximum. We also compare experimentally measured spin correlations as a function of the measured singlon density n_s to T = 0 data from CP-AFQMC. The singlon density is directly measured in the experiment and does not involve converting n_s to doping using numerical data. We confirmed that a miscalibration of U/t seems insufficient to explain the observed discrepancy. However, higher-order corrections to the experimental Hamiltonian, including density-dependent terms, may lead to systematic deviations at finite doping from the calibration performed at half-filling. A complete comparison between experimental and numerical data requires a comprehensive characterization of the experimental Hamiltonian and the development of calibration techniques that directly probe the doped regimes.

The constrained-path approximation used in the numerical simulations may also introduce systematic errors. To explore this possibility, we compare experimental and CP-AFQMC data at elevated temperatures with numerically exact DQMC simulations. When simulating the doped Hubbard model, the sign problem in DQMC leads to an exponential overhead with decreasing temperature. Given accessible computational resources, we first perform DQMC at T/t = 0.25 and U/t = 8 in an 8 × 8 system and compare the results to those obtained with CP-AFQMC with the same parameters but in a 12 × 12 system (Fig. 4). Both the DQMC and CP-AFQMC simulations use periodic boundary condition (PBC). For these DQMC simulations, we used the QUEST package⁶⁴ with a time step of δτ = 0.02, 5,000 thermalization sweeps and 200,000 measurement sweeps, which are repeated with random seeds for up to 4,000 runs for values of the density with severe sign problems. The experimental data are obtained with U/t = 8.2(2) and using the standard loading sequence (without engineering the DMD potential) in refs. ^12,63. The harmonic confinement due to the lattice intensity profile leads to a slow variation of the filling from the centre to the edge of the trap. We tune the atom number such that it reaches half-filling at the centre of a region of radius r = 6. This allows us to perform accurate thermometry of the entire atom cloud by comparing with numerically exact simulation results at half-filling without the need for numerical results in the doped regime. At larger distances from the trap centre of r > 6, the filling slowly decreases. Radially binning the data allows us to probe spin correlations as a function of varying singlon density and, therefore, doping at the calibrated temperatures.

We focus on the short-range spin correlations ${C}_{d}^{zz}$ up to $d=\sqrt{5}$ in Extended Data Fig. 5a,b. Given the short correlation lengths at this elevated temperature, and at finite doping, finite-size effects are expected to be negligible. At a density of n = 0.995 and chemical potential of μ/t = −1.0, in which the correlation lengths are expected to be longer than at larger doping, we find that the DQMC and CP-AFQMC simulations are in good agreement for all short-range correlations. Moreover, at T/t = 0.33, we confirm that there is no difference between DQMC simulations performed in 8 × 8 and 12 × 12 lattices.

We find that good agreement between experimental data, DQMC and CP-AFQMC on the nearest-neighbour correlation ${C}_{1}^{zz}$ holds for all doping values at which we performed measurements or simulations. However, for a doping range of δ ∈ [5%, 10%] and bond distances $d=\sqrt{2},2,\sqrt{5}$, the results obtained from DQMC suggest stronger correlations, with which experimental data show good agreement (Extended Data Fig. 5a). This trend is consistent with the behaviour shown in Fig. 4. The large statistical error bars in the DQMC data make it difficult to draw a definitive conclusion, especially because the estimated errors are themselves unreliable because of the vanishing average signs.

At even higher temperatures of T/t = 0.33, the sign problem is less severe, and we can compare DQMC and CP-AFQMC in a 12 × 12 system. We find smaller, but qualitatively similar, discrepancies in which spin correlations ${C}_{\sqrt{2}}^{zz},{C}_{2}^{zz}$ computed from DQMC are stronger than CP-AFQMC for doping below 15% (Extended Data Fig. 5b). Owing to the reduced magnitude of the spin correlations, experimental data performed at this temperature seem to be statistically consistent with both numerical simulations. Future work should be able to address this by using a better trial density matrix or an improved form of the self-consistent constraint in CP-AFQMC and more systematic studies.

Lattice potential calibration

To measure the lattice trapping potential, we measure the density profile of a non-interacting spin-polarized Fermi gas loaded into the lattice, and, by taking the local density approximation, invert the non-interacting equation of state to obtain the local chemical potential. The spin-polarized Fermi gas is prepared by performing evaporative cooling with states 1 and 2 at 321 G, followed by a magnetic gradient-assisted spill-out of state 1 at 27 G, where state 2 is magnetically insensitive and therefore remains trapped. Pauli exclusion in a spin-polarized Fermi gas prevents double occupancy (and associated parity projection during imaging), meaning that the measured density can be mapped unambiguously to a particular value of the chemical potential. The uncertainty in the above calibration procedure is dominated by statistical errors and is not strongly dependent on the assumed temperature of the Fermi gas. We take a conservative estimate of T/t = 0.5 based on independent calibrations of the spin-polarized Fermi gas, but the results do not change significantly when assuming a temperature in the range T/t ∈ [0, 1].

Using the above procedure, the harmonic confinement of doped data (Fig. 4) is measured as V_H = 0.0152(6)t/(site)²r², where r denotes the radial distance measured in sites from the lattice centre.

The half-filled data at U/t = 8.3(2) (Fig. 3) was taken before the addition of the non-reciprocal attenuator; so an offset between the A and B sublattices is present. Using the procedure described above, we find that the difference in the mean local potential in the A and B sublattices is Δμ = 0.75(3)t. This is in good agreement with the expected offset of 0.8(4)t given the geometry of the lattice. The uncertainty in the expected offset is primarily because of uncertainty in the applied X lattice depth at very low values.

To investigate the effects of sublattice offsets on spin correlations, we perform DQMC simulations at T/t = 0.15, 0.3 on an 8 × 8 system and CP-AFQMC simulations at T/t = 0 on a 12 × 12 system. DQMC simulation at T/t ≥ 0.15 show that an offset as large as Δμ = 2t (defined as a symmetric sublattice offset about mean μ = 0) has little effect on spin correlations (Extended Data Fig. 4a). At T/t = 0, CP-AFQMC shows no effects on the spin correlations for Δμ = t and a small decrease in spin correlations for Δμ = 2t. This decrease is smaller than the statistical error in the experimental data; so we conclude that sublattice offsets are not a concern for the above dataset.

Characterization of BI

BI fidelity

We characterize the fidelity of the BI using the singlon density F_BI = 1 − n_s, assuming that a perfect BI has a doublon on every site. We can image the BI state with parity projection in the long-spacing lattice, in which doublons are lost because of light-assisted collisions⁶³ (Fig. 1c, left). Alternatively, we can directly image the doublon population by splitting each long-spacing site into two after freezing the dynamics to reconstruct the population with full density resolution⁵⁷ (Fig. 1c, middle).

Using density-resolved imaging, we measure the doublon density to be n_d = 98.2(5)%, the singlon density to be n_s = 1.8(5)%, and the hole density n_h to be consistent with zero within a central ROI that does not include the boundary of the crater. The ROI covers a circular region with a radius of r = 6 sites in the long-spacing lattice (about r = 9 in the short-spacing lattice), which is one site (two sites) smaller than the radius of the crater. Note that given the fidelity of doublon detection in the experiment, the above numbers are consistent with a doublon fraction of unity. The total number of atoms detected in the crater is N_r=10 = 342(1).

The above measurements confirm that the empty sites appearing in the parity-projected images are doublons instead of holes. Given this observation, parity-projected imaging offers better sensitivity to the doublon population because it is immune to atom loss during fluorescence imaging. The remaining errors in parity-projected imaging of doublons include the fidelity of removing doublons through light-assisted collisions and atoms hopping to a different 2D layer during imaging. At present, we do not have an independent calibration of these errors; so the measurements of the BI fidelity below constitute a lower bound.

Using parity-projected imaging, we find that most of the singly occupied sites occur at the edge of the crater. We note that the crater is not aligned to the optical lattice in a site-resolved manner. When combined with the finite resolution of the imaging system, this results in imperfectly controlled local potentials on the sites falling along the edge of the crater. The equation of state on these sites may not favour doublons. Moreover, the local density approximation may not hold in the presence of this abrupt potential variation. We minimize the population of singlons on the edge by increasing the lattice depths and DMD0 potential strengths, such that the tunnelling energy t is small compared with the site-to-site potential difference at the crater edge. We also set the magnetic field at 550 G to lower the Hubbard interaction U, which prevents the formation of a wide Mott plateau. The detected singlon density is n_s,r=4 = 0.5(3)% in a central disk with r = 4, and n_s,r=6 = 0.7(2)% in a disk with r = 6. Including the edge, the singlon population in the crater is N_r=10 = 11.9(7). However, the presence of singlons on the edge obscures the estimation of entropy in the BI because these singlons are not necessarily excitations. By contrast, because the density of states along the crater edge can be finite, these singlons can host a significant amount of entropy. Taking this into account, we choose to define the BI fidelity F_BI = 1 − n_s,r=6 in the central region with r = 6 as a figure of merit to optimize the entropy redistribution procedure.

To optimize the initial loading of the BI, the DMD0 potential ramp is split into a linear ramp and a holding phase. Once the lattice depths are high and the tunnelling energy correspondingly small, the atoms may not redistribute for continued changes in the DMD0 potential (that is, the ramp is not adiabatic). We also limit the initial ramp rate to avoid creating band excitations in shallow lattices.

The volcano-shaped DMD0 potential allows the formation of a BI that is in thermal contact with a reservoir when the lattice is shallow, and isolated when the lattice is deep. When using the flattened volcano-shaped potential, we find that the quality of the separation between system and reservoir deteriorates if the plateau surrounding the crater becomes too large, because some atoms may remain in the region with no potential gradient.

Heating in the BI

Heating processes in a BI are strongly suppressed because of the vanishing density of states, and the only active processes must excite atoms into higher bands. Relevant heating mechanisms include incoherent light scattering from the optical lattice, intensity and position noise of the optical lattice, and background gas collision. To quantify the loss, we hold the BI for a variable duration of up to 1,600 ms and measure the total density n using both parity-projected imaging and full charge-resolved imaging as a probe of BI fidelity. In the parity-projected imaging, we measure the defect singlon density in the long-spacing lattice. In the full charge-resolved imaging, we measure the density after splitting the doublons into singlons in the short-spacing lattice. We find that the loss rate of density in the lattice centre is dn_r=6/dt = 5.1(3)% s⁻¹ using full charge-resolved imaging, and dn_r=6/dt = 4.8(4)% s⁻¹ using parity-projected imaging, assuming that the detection of a singlon corresponds to the loss of a single particle (Extended Data Fig. 3c,d). The zero-time intercept extracted from charge-resolved imaging is n = 1.983(3). The intercept from parity-projected imaging is n = 0.008(4), which is consistent with the directly measured singlon density. The difference between the parity-projected measurement and the charge-resolved measurement is consistent with the limitations imposed by imaging fidelity. Given the loading time of 100 ms, the above loss rate is consistent with the measured BI fidelity of F_BI,r=6 = 99.3(2)%.

Calibration of Hubbard interaction U

Owing to light-assisted collisions, in the short-spacing lattice, we can only perform parity-projected density imaging, and therefore detect the density of singly occupied sites ${n}_{{\rm{s}},det}$ (ref. ⁶³). We correct the detected singlon density by the imaging fidelity to estimate the actual singlon density n_s. At half-filling (Figs. 2 and 3), the singlon density is a function of interaction and temperature: n_s(U, T). We rely on numerical simulations using DQMC to obtain the doublon density n_d(U, T) and compute the singlon density as n_s(U, T) = 1 − 2n_d(U, T). We find that n_d(U, T) is sensitive to interaction strengths but has only a weak dependence on temperature for 0.15 < T/t ≤ 0.25 and saturates for T/t ≤ 0.15 (Extended Data Fig. 4c). Using both the singlon density and spin correlations ${C}_{{\bf{d}}}^{zz}$, we estimate the possible ranges for interaction and temperature as 8 < U/t < 9 and T/t ≤ 0.15. The fact that n_d is insensitive to temperature, and therefore effectively becomes a function of only interaction, allows us to invert this equation of state to calibrate U/t using only the measured singlon density n_s. With U/t calibrated, we can then estimate the temperature using spin correlations (see section ‘Calibration of temperature T at half-filling’).

To obtain the doped data shown in Fig. 4, the lattice depths and magnetic fields are slightly different from the ones used to obtain the half-filling data shown in Fig. 3 (see section ‘Experimental methods’). Therefore, interaction strength U/t needs to be recalibrated for the doped systems. We adjust the amount of expansion to prepare a sample with half-filling n = 1 in the centre, which can be identified as the maximum of singlon density n_s as filling n is increased. We then apply the same calibration protocol described above to this half-filled state to obtain the Hubbard interaction U/t given by the lattice and scattering length parameters used for the doped systems.

There is no available numerical simulation for the large U/t Hubbard model at low temperatures, so a different approach is required to calibrate the interaction strength U/t for the data taken in the Heisenberg limit. We first perform measurements at 580 G in the same lattice configuration as the final data, and with T/t ≤ 0.15 such that the temperature dependence of the singlon density is negligible. The resulting measured value of the singlon density is n_s(580 G) = 0.898(7), corresponding to U/t = 8.3(4). Next, using a_s(620 G)/a_s(380 G) = 312/233 ≈ 2.18, we obtain U/t = 18.6(8). If we ignore the temperature dependence of the singlon density and compare with numerical linked cluster expansion data at T/t = 0.2 (ref. ⁶⁵), our measured singlon density of n_s(620 G) = 0.976(3) gives 19.4(1.7), which is consistent with the above measurement.

Estimation of doping

To convert singlon density n_s to density n, we rely on numerical simulations using CP-AFQMC to obtain the doublon density n_d(n, U, T) as a function of density, interaction strength and temperature. From this, we can compute the singlon density n_s(n, U, T) = n − 2n_d(n, U, T). Similar to half-filling, the variation of doublon density reduces with temperature T and becomes negligible for T/t ≤ 0.15 (Extended Data Fig. 4b).

We plot the experimentally measured spin correlations as a function of singlon density n_s together with numerically simulated spin correlations as a function of computed singlon density n_s(n, U, T) in Extended Data Fig. 3f. CP-AFQMC simulation data of n_d is not shown in the density range 0.85 < n < 1 because of poor convergence. Based on comparisons of the nearest-neighbour spin correlations, the data are consistent with temperatures of T/t ≤ 0.15 for all values of the doping for which CP-AFQMC is performed.

Taken in combination, the above observations allow us to invert the equation of state and estimate the doping δ = 1 − n(n_s, U, 0) using numerical simulations performed at T/t = 0, which have been studied systematically in previous works¹⁷.

Calibration of temperature T at half-filling

For the half-filled data taken at U/t = 8.3(2), good agreement between experimental and numerical data (Fig. 3) allows us to measure temperature using a physically motivated observable that aggregates information from spin correlations beyond nearest neighbours. Specifically, we consider the staggered magnetization m^z, which is defined as a sign-corrected average of the spin correlation function up to a cutoff d in bond distance:

$${({m}^{z})}^{2}=\frac{1}{{{\mathcal{N}}}_{{\Omega }_{d}}}\sum _{(i,j),{i}^{2}+{j}^{2}\le {d}^{2}}{(-1)}^{i+j}{C}_{(i,j)}^{zz}\,.$$

(3)

We estimate the value of ${m}_{\text{exp}}^{z}$, as well as the 1σ error $\delta {m}_{\text{exp}}^{z}$ in this estimate, from bootstrapping, and compute the expected value of ${m}_{{\rm{num}}}^{z}$ as a function of temperature using the DQMC.

We check the effect of boundary conditions by obtaining DQMC numerical data on 12 × 12 and 16 × 16 square systems with open boundary conditions and PBCs (Extended Data Fig. 7a). PBCs tend to overestimate spin correlations ${C}_{d}^{zz}$ at long range, whereas the difference between 12 × 12 and 16 × 16 data is not substantial for d ≤ 6. This motivates our choice of performing subsequent calculations using open boundary conditions and a 12 × 12 system.

Spin correlations are computed at U/t = 8 and U/t = 9 (Extended Data Fig. 7b), converted into magnetization ${m}_{{\rm{num}}}^{z}$ and interpolated to produce a continuous function of temperature and interaction strength ${m}_{{\rm{num}}}^{z}={f}_{U}(T)$ (using cubic splines along T and linear functions along U). At a fixed value of U based on experimental calibrations, we convert our knowledge of the true experimental magnetization, which we model as a Gaussian random variable $\widetilde{{m}^{z}} \sim {\mathcal{N}}({m}_{\text{exp}}^{z},{(\delta {m}_{\text{exp}}^{z})}^{2})$, into temperature by computing the probability distribution function of $\widetilde{T}=g(\widetilde{{m}^{z}})$, where g is the inverse of the interpolated numerical data, $g({m}^{z})={f}_{U}^{-1}({m}^{z})$ if m^z ≤ f_U(0), and g(m^z) = 0 if m^z > f_U(0). We report the mean value and the ±1σ confidence interval of the temperature distribution.

Experimental and numerical magnetizations and their associated temperature estimates are shown in Extended Data Fig. 7c,d. At the experimentally calibrated value of U/t = 8.3, increasing the bond cutoffs from d = 3 to d = 6 (Extended Data Fig. 7a, vertical dashed lines) does not change the inferred mean temperature of T = 0.05. This indicates a good agreement between numerical and experimental data at all ranges and confirms the calibrated interaction strength (Extended Data Fig. 7c). With a cutoff d = 3 on bond distance, below which numerical data only weakly depends on boundary conditions (Extended Data Fig. 7a), the inferred mean value of T/t varies by about 0.01 within the reported uncertainty on interaction strength U/t = 8.3 ± 0.2, with again a weak variation of the confidence interval. We report in the main text a temperature and asymmetric errors of $T/t=0.0{5}_{-0.05}^{+0.06}$ obtained with a cutoff d = 6, and a single notable digit capturing the overall magnitude of the systematic errors related to interaction calibration and finite-size simulation effects.

Calibration of tunnelling and lattice parameters

As described in section ‘Experimental methods’, the key step to ramp from a dimerized lattice to the short-spacing lattice is to ramp down the long-spacing lattice by reducing the X lattice depth and to ramp up the short-spacing lattice by increasing the $\bar{X}$ lattice depth. Ramping the interference time phase ϕ from 0 to π/2 eliminates tunnelling dimerization but introduces a potential offset between the A and B sublattices. The non-reciprocal attenuator allows for ϕ to be stabilized over a much larger range of X lattice depths, allowing us to suppress tunnelling dimerization without introducing potential offsets. We use the calibrated lattice depths and interference parameters^12,57 to numerically compute the band structure of the 2D lattices. The absolute amplitudes of tunnellings (in Hz) can then be computed by constructing maximally localized Wannier orbitals⁶⁶ or by fitting the band structure with a tight-binding model. We find that the two methods give consistent results for our lattice geometry (Extended Data Fig. 6a) and report the resulting tunnelling amplitudes in Extended Data Table 1.

The tunnelling dimerization is given by the interference term in the lattice potential, which scales as $\sqrt{{V}_{{x}}}$. It is important to quantitatively determine the dependence of intra- and inter-dimer tunnelling energies t_d, t_i and the perpendicular tunnelling energy t_p on V_X. With lattice depths V_Y and ${V}_{X}+{V}_{\bar{X}}$ fixed, the variation of the different tunnelling energies with V_X is shown in Extended Data Fig. 6a. At the end of the ramp V_X = 3.2(3) × 10⁻⁵E_R, and the tunnelling energies are all balanced within systematic uncertainties. To obtain the doped data, the X is ramped down to V_X = 3.2(3) × 10⁻⁵E_R in the shallow lattice condition and kept the same for the rest of the experimental sequence.

Experimental verification

Adiabaticity of splitting at strong interaction

In the strong interaction limit U/t = 18.6(8), the adiabaticity is decided relative to the many-body spin gap for an ideal finite-size system with sharp open boundaries and homogeneous potentials. The gap separating the ground state and the Anderson tower of states scales as 1/L², whereas the gap for spin wave excitations scales as 1/L, where L is the linear system size^15,48. In the absence of heating, slower ramp speed results in more adiabatic evolution.

Experimentally, as we increase the ramp duration from 20 ms to 80 ms, crossing the critical point predicted by previous numerical works⁴⁷, we find the spin correlations in the final state reach saturation and do not increase (Extended Data Fig. 6g). We suspect that for longer ramp durations heating competes with the improvement of adiabaticity. Holding during or after the ramp results in strong heating, which manifests as a reduction of spin correlations. We also find that heating is related to the magnetic bias field. When holding in the middle of the ramp, the heating is much stronger at 620 G, where we operate to achieve U/t = 18.6(8), than at 550 G, where the BI is prepared. This hints at heating due to three-body processes, which are strongly enhanced as the scattering length a_s is increased.

The saturation of spin correlations with increasing ramp time is measured at the same DMD potential strengths as when the BI is formed. We find that with weaker DMD potentials, 20 ms ramps are no longer adiabatic and result in reduced singlon density and spin correlations (Extended Data Fig. 6g). This suggests that atoms are leaving the crater region. However, we do not have a complete understanding of this dynamical process because the time dependence of the magnetic field ramp also plays an important part. Higher magnetic bias fields lead to stronger interaction strengths, which compete against the confining potential and may alter the physics of the transport process.

Adiabaticity of splitting compared with tunnelling dimerization ramp

In a homogeneous Hubbard system at half-filling, or in a Heisenberg spin system, the adiabaticity criterion is determined by the speed of the ramp compared with the energy gap of spin excitations¹⁵. In our case the relevant ramp is of the tunnelling dimerization, which is controlled by the interference phase. Experimentally, at U/t = 18.6(8), we find that a ramp time longer than 5 ms is sufficiently adiabatic and allows intra- and inter-dimer spin correlations to equalize. Ramps with a duration of less than 1 ms are diabatic and result in clear deviations. This ramp has no detectable effects on longer-range correlations, which is probably decided by the rate of the previous ramp across the quantum critical point. For longer ramp times, we find significant heating when the DMD0 potential is not optimized, leading to atom transport. The coldest reported temperatures are, therefore, obtained by dynamically balancing the confinement of the DMD0 potential, the tunnelling energy and interaction strengths. At U/t = 8.3(2), we find that heating is much less problematic and choose a ramp time of 25 ms to ensure adiabaticity. A systematic study of the phase transition and the resulting requirements on adiabaticity would require proper flattening of the potential and site-resolved preparation of the initial BI.

Characterization of tunnelling dimerization in doped systems

The lattice phase stabilization is kept active between consecutive experimental sequences, and the interference time phase is maintained. The optical lattice position is, therefore, tracked at the single-site level in fluorescence atom images. This allows us to probe spin correlations in a bond-resolved manner (see section ‘Data analysis’). Different y, intra-dimer and inter-dimer tunnelling energies t_y, t_d, t_i lead to different superexchange interactions ${J}_{y},{J}_{{\rm{d}}},{J}_{{\rm{i}}}=4{t}_{y}^{2}/U,4{t}_{{\rm{d}}}^{2}/U,4{t}_{{\rm{i}}}^{2}/U$ on the respective bonds, and thus different spin correlations. In Fig. 2, we show how ramping the interference time phase ϕ = π/2 eliminates the tunnelling dimerization and correlation imbalance in the Heisenberg limit. To obtain the doped data, at the end of the ramp to V_X = 3.2(3) × 10⁻⁵E_R (using the non-reciprocal attenuator), we observe that the y, intra-dimer and inter-dimer spin correlations agree to within statistical errors (Extended Data Fig. 6d). The density profile also shows no offset between the A and B sublattices, confirming the absence of potential offsets even with ϕ = 0 (Extended Data Fig. 6b,c).

Tunnelling dimerization has striking effects on the doped Hubbard systems. In contrast to the half-filled case, in which long-range correlations are negligibly affected by dimerization under the conditions V_X = 0.016E_R, ϕ = 0, the doped systems show only correlations within the dimers. This suggests that dimerization will lead to different physics in the doped Hubbard systems, so it is crucial to remove this dimerization to observe the physics of the isotropic doped Hubbard model.

Adiabaticity of expansion after splitting

Splitting a band insulating state naturally gives rise to a half-filled state. To obtain doped Hubbard systems, we decrease the DMD0 potential such that the atoms can expand out of the crater. The boundary of the expansion is set by the wall formed by the DMD1 potential. Expansion in this ring-shaped region can be strongly affected by the shape of the DMD0 potential. We find that a volcano-shaped DMD0 potential (Extended Data Fig. 1f) with a sharp, outward slope (in which the potential decreases as distance from the centre increases) optimizes the separation between the BI and the reservoir. However, it is more difficult to expand the atoms adiabatically with this potential. The sharp slope acts as a potential barrier separating the crater and the region closer to the DMD1 wall potential, which forms two potential minima in the radial direction. This is similar to the broken symmetry state in a double-well system. The initial BI state is elevated from the true ground state with a small gap given by the tunnelling amplitude between the double well, which is exponentially suppressed with increasing strength of the potential barrier created by DMD0. An alternative, classical picture is that as the atoms flow out of the crater, they will accelerate downhill on the sharp slope, which is not a reversible process and therefore not adiabatic. A bowl-shaped potential that increases until reaching the DMD1 wall as the radial distance from the centre increases alleviates the above issues but will trap additional atoms during the BI preparation leading to increased entropy of the initial state. To balance these considerations, we choose the flattened volcano-shaped potential as a compromise (Extended Data Fig. 1g). We confirmed that the energy gap during expansion is similar to the case in which the potential has an inward slope using exact diagonalization in one dimension.

Larger tunnelling and smaller U/t leads to more adiabatic expansion for fixed expansion duration, suggesting that shallow lattices are favourable. However, we find this does not hold true in the limit for which the lattices are almost fully extinguished. We perform round-trip measurements of the BI fidelity by first decreasing the lattice depths to the values used in expansion and then increasing the lattice depths back to those used to obtain the BI state. Excitations created during the ramp will manifest as singlon defects in the BI. We find a significant increase in the number of these defects when lattice depths are decreased to below 0.5E_R during expansion. As a result, we opt to use a lattice depth of 0.3E_R during expansion, which corresponds to about 20% of the lattice depth used to load the BI. We believe this behaviour is because of band gaps becoming smaller than the tunnelling energy, allowing diabatic excitations to higher bands. After ramping to these conditions, we start the splitting procedure, which merges the two lowest bands in the long-spacing lattice into a single band in the short-spacing lattice.

We experimentally investigate the dynamics during expansion by varying the duration of the expansion step. We measure the final spin correlations after reloading for the expansion times of τ_exp = 1 ms, 15.5 ms and 30 ms. No statistically significant variation is observed in these measurements (Extended Data Fig. 6e), suggesting that 1 ms is sufficient for the expansion to be adiabatic. This is consistent with the large bandwidth and Fermi energy of around 10 kHz during expansion. To further facilitate thermalization during expansion, and to allow the field to settle before reloading happens, we ramp the magnetic field from 550 G to 590 G during the expansion step, which increases the scattering length a_s from around 86a₀ to 295a₀, and set the expansion time to a much longer duration of 50 ms.

Adiabaticity of expansion with fixed duration

In Fig. 5, we show the adiabaticity of expansion in the long-spacing lattice with fixed duration in units of the tunnelling time ħ/t. The temperatures are obtained by comparing the d = 1 and d =  $\sqrt{2}$ spin correlations ${C}_{d}^{zz}$ averaged in the central half-filled region to DQMC simulations. The goal is to probe adiabaticity as a function of interaction strengths, independent of the absolute value of tunnelling amplitudes. Although expansion in deeper lattices may lead to more lattice heating because of the increased laser intensity and smaller tunnelling amplitudes, little heating was observed in previous work⁵⁷ involving measurements performed at similar lattice depths and ramp times as for the deepest lattice explored in Fig. 5 (η = 1). The increased kinetic energy scales at shallower lattices may also reduce the effects of potential disorder, whose contribution to expansion adiabaticity requires future studies. We also measured the adiabaticity of expansion with a fixed duration of τ = 120 ms as shown in Extended Data Fig. 3e, which corresponds to 18.5ħ/t for η = 1. We find similar heating as in Fig. 5a despite the holding time being about four times shorter. This suggests that diabatic heating and dissipative heating are occurring on a similar scale. According to Fig. 5b, however, 120 ms should be slow enough to incur only about 0.2t increase in temperature, even at η = 1 (it would correspond to roughly 4 ms in Fig. 5). This indicates that interactions probably play a part in inhibiting transport at τ = 120 ms by placing more stringent requirements on adiabaticity.

Thermalization of splitting at half-filling

To check whether the half-filled state in Fig. 3 has equilibrated, we added a hold time of τ_h = 0 ms, 10 ms, and 20 ms after ramping into the final antiferromagnetic state. During this measurement, the lattice alignment relative to the DMD potentials was not optimized, which may explain the higher resulting temperatures. We find no change in the singlon density as a function of radial distance r from the trapping potential centre, or in the overall profile of the spin correlations as a function of bond distance d (Extended Data Fig. 6h). However, we observe a slight decrease of the nearest-neighbour correlation, which may be attributed to heating during the hold time. This suggests that the ramp to split a BI into a half-filled system realizes a state in thermal equilibrium.

Thermalization during reloading

After unloading the BI, splitting into a half-filled Fermi liquid, and expanding into a Fermi liquid with filling n < 1, the sample is in a similar state to previous works^11,63, but at lower entropy. We reload the lattice using a linear ramp-up of the lattice depth in 100 ms, which is chosen to balance adiabaticity and lattice heating. Moreover, we experimentally confirmed the final quantum state is in equilibrium by holding for τ_h = 0 ms, 20 ms and 40 ms, corresponding to 0ħ/t, 20ħ/t and 40ħ/t. No statistically significant variation is observed for these hold times (Extended Data Fig. 6f).

Furthermore, in the doped system, we can cross-check the thermalization throughout the sample by comparing temperature estimates obtained at different locations in the cloud. As described in section ‘Data analysis’, the temperatures shown in the main text are obtained in an ROI with radii r = 3, 4 and 5 in the centre of the cloud by including only pairs of sites within the ROI (which we refer to as strict binning). We can also obtain temperature estimates from the radially averaged spin correlations, although these estimates suffer from averaging over larger potential variations and different bin shapes. Nevertheless, in Extended Data Fig. 8, we reanalyse the data from Fig. 4f to plot ${C}_{{\rm{bin}}}^{zz}(d)$ against the measured singlon density n_s,bin. ${C}_{{\rm{bin}}}^{zz}(d)$ denotes the measured zz spin correlation at bond distance d for the corresponding bin. Note that the above analysis does not apply strict binning, so only one of the two sites used to compute a given correlation function must be contained within the bin. This leads to systematic deviations from the results obtained with strict binning. For example, the correlations with longer bond distances will be underestimated for the bin with the shortest distances, as sites at lower densities are involved. Nevertheless, this analysis offers a qualitative view that the atom cloud has thermalized across different densities.

Exact diagonalization of BI splitting

We use the QuSpin package⁶⁷ to compute the evolution of energy levels during splitting with exact diagonalization. For the Heisenberg model (Fig. 2), we choose a system size of 6 × 5, with a PBC in the long direction and twisted boundary condition in the short direction. The twisted boundary condition is defined as connecting (i, 5) → ((i + 1)%6, 1). The spin exchange coupling along the y-direction is fixed to J = 1, and $\alpha =\sqrt{{J}_{{\rm{d}}}/J}=\sqrt{J/{J}_{{\rm{i}}}}$ are varied from 10 to 1.

Numerical simulations at half-filling

Our simulations at half-filling were performed with finite-temperature DQMC and ground-state AFQMC methods. AFQMC and DQMC are numerically exact formulations for studying quantum many-body systems. The approach begins with the Trotter decomposition, which breaks the original imaginary-time propagator into smaller pieces. To handle the on-site interaction, an auxiliary field is introduced using the Hubbard–Stratonovich transformation, transforming the interaction operator into a sum of one-body operators. These steps allow the partition function to be expressed as a sum over the products of two Slater determinants—one for spin-up electrons and one for spin-down electrons—whose matrix elements can be computed analytically. However, the product may take both positive and negative values, leading to a sign problem (except in specific cases, such as the half-filled 2D repulsive Hubbard model, in which particle–hole symmetry eliminates this issue).

In our simulations (Fig. 3), finite-temperature results are obtained using the DQMC method with 7,000 thermalization sweeps, 200,000 measurement sweeps and a time step of Δτ = 0.05. Ground-state data are generated with the AFQMC method, using the Metropolis algorithm with force bias updates⁶⁸, thousands of sweeps and a time step of Δτ = 0.02. We verify that the Trotter error at these values of Δτ is smaller than the statistical error, ensuring the robustness of our results. To mitigate ergodicity issues and fluctuations in spin correlations caused by the SU(2) symmetry breaking in the spin-z Hubbard–Stratonovich transformation, we adopt the charge Hubbard–Stratonovich transformation for both the finite-temperature and ground-state calculations. Additionally, the infinite variance problem⁶⁹ in spin correlation measurements is significantly reduced when using the charge decomposition instead of the spin decomposition.

CP-AFQMC simulation

For doped systems, we use both finite-temperature^16,70 and ground-state⁷¹ CP-AFQMC methods to compute the properties of the system at zero and finite temperatures, respectively. The finite-temperature CP-AFQMC method shares several of the building blocks from the standard DQMC method, including the use of Trotter decomposition and the Hubbard–Stratonovich transformation. However, it introduces two key differences. First, in systems suffering from the sign problem, an exact condition is derived for the paths in auxiliary-field space that cause the sign problem⁷⁰. It is shown that only a small subset of paths in the auxiliary-field space contribute to the partition function. The remaining paths are symmetrically distributed and cancel each other, which adds noise to the computation of observables. To address this, the CP-AFQMC method imposes a constraint, acting as a gauge condition, to isolate and sample the relevant paths. This constraint is exact when the many-body Hamiltonian is used in the constraint. Second, to implement the constraint and sample configurations efficiently, the algorithm uses a branching random walk with importance sampling to construct complete paths. In practical applications, an effective Hartree–Fock Hamiltonian is used as a trial Hamiltonian to define the constraint¹⁶. The parameters U_eff and β_eff are tuned so that the spin correlations obtained by the effective Hamiltonian closely match those of the many-body system¹⁸.

Ground-state CP-AFQMC⁷¹ is a projection-based quantum Monte Carlo method that leverages the principle that the ground state can be accessed by applying an imaginary-time evolution operator to an initial wave function, provided the trial wave function has a nonzero overlap with the ground state. Similar to the finite-temperature algorithm, ground-state CP-AFQMC also combines the Trotter decomposition with the Hubbard–Stratonovich transformation to reformulate the imaginary-time evolution as a stochastic process—a random walk through the space of Slater determinants. Ground-state properties are calculated as statistical averages over the configurations sampled during this random walk. To address the fermion sign problem, an unrestricted Hartree–Fock solution is used as a trial wave function, introducing a constraint that ensures each Slater determinant in the random walk maintains a positive overlap with the trial wave function. Extensive benchmarks in the literature^21,23,27,28 demonstrate that CP-AFQMC delivers high accuracy for studying ground-state and finite-temperature properties of quantum many-body systems. Our ground-state AFQMC calculations followed similar procedures to ref. ¹⁷ (although at half-filling we applied the charge decomposition as discussed earlier). In our simulations (Fig. 4), the Trotter time step is typically 0.05 and 0.01 in finite temperature and ground-state calculations, and we verify that the Trotter errors are limited. The population of random walkers is around 5,000 in these calculations.

Quantum Monte Carlo for Heisenberg model

Finite-temperature numerical simulations of the Heisenberg model in Fig. 2 and Extended Data Fig. 9 are performed with the SpinMonteCarlo.jl package (www.juliapackages.com/p/spinmontecarlo), using a loop quantum Monte Carlo algorithm with 1,024 warmup updates and 8,192 measurement updates, over 16 × 16 dimerized lattice unit cells (512 sites) with PBCs. Individual spin correlators $\langle {S}_{z}({\bf{r}}){S}_{z}({{\bf{r}}}^{{\prime} })\rangle$ are grouped by pairs of sites $({\bf{r}},{{\bf{r}}}^{{\prime} })$ with the same bond vector ${\bf{r}}-{{\bf{r}}}^{{\prime} }$, distinguishing pairs that belong to the same sublattice from those that do not. These correlators are averaged at each measurement update, resulting in a relative statistical error of the order of 10⁻³. The resulting correlation maps are shown in Extended Data Fig. 9a for representative values of the coupling parameter α at temperature T/J = 0.5. We estimate that the systematic errors due to finite-size effects are at most 5 × 10⁻³ in the square lattice at this temperature, based on a comparison to a 32 × 32 system.

Staggered magnetization m^z is obtained by averaging the sign-corrected spin correlations over bond distance (Extended Data Fig. 9b), in this case truncated to a square 13 × 13 region (Extended Data Fig. 9a), and matching the ROI used in the analysis of experimental data. We observe a distinct increase and temperature dependence of m^z at couplings above α = 0.65, which we attribute to the existence of a quantum phase transition to an antiferromagnetic phase at zero temperature. This critical value is comparable to ground-state Heisenberg simulations performed in similar dimerized lattice geometries⁴⁷—note that the parameter α in our work is defined as a ratio of tunnel couplings and is related to ratios of antiferromagnetic spin couplings in the equivalent Heisenberg model by a square root.

After normalizing experimental data by the square spin moment ${n}_{{\rm{s}}}^{2}$, where n_s = 0.95 is the averaged singlon density in this dataset, we observe good agreement between numerical Heisenberg data compared with bond distance and experimental spin correlations $| {\bf{r}}-{{\bf{r}}}^{{\prime} }|$ in the square lattice at U/t = 18.6(8) (Extended Data Fig. 9c). A least-squares fit of correlations up to $| {\bf{r}}-{{\bf{r}}}^{{\prime} }| \le 8$ results in an estimated temperature of T/J = 0.458(3).

Effective temperature in Hubbard systems compared with cuprates

The single-band Hubbard model can be related to the high-T_c superconducting cuprates through a few steps of reduction and approximation. A more realistic approximate model for the cuprates is the three-band Hubbard model, in which tunnelling t and interaction strength U can be correspondingly related to energy scales in cuprates². Therefore, it is not straightforward to relate the temperatures in the one-band Hubbard model with doping to an effective temperature in cuprates.

However, at half-filling, the physics is well understood in the Hubbard model³ and cuprates^2,19,20 and can be described by an antiferromagnetic Heisenberg model. In this regime, the characteristic energy scale is the exchange energy J, which can be readily probed in cuprates and is J = 4t²/U in the one-band Hubbard model. This correspondence allows us to convert the temperature in the half-filled Hubbard model using T/J to an effective temperature in cuprates. Taking J = 125(5) meV in yttrium barium copper oxide^19,20 and U/t = 8, a temperature of T/t = 0.25 corresponds to about 725 K, and a temperature of T/t = 0.05 corresponds to around 145 K.