The emergence of quasiparticles in quantum many-body systems underlies the rich phenomenology in many strongly interacting materials. In the context of doped Mott insulators, magnetic polarons are quasiparticles that usually arise from an interplay between the kinetic energy of doped charge carriers and superexchange spin interactions [1][2][3][4][5][6][7][8] . However, in kinetically frustrated lattices, itinerant spin polaronsbound states of a dopant and a spin flip-have been theoretically predicted even in the absence of superexchange coupling [9][10][11][12][13][14] . Despite their important role in the theory of kinetic magnetism, a microscopic observation of these polarons is lacking. Here we directly image itinerant spin polarons in a triangular-lattice Hubbard system realized with ultracold atoms, revealing enhanced antiferromagnetic correlations in the local environment of a hole dopant. In contrast, around a charge dopant, we find ferromagnetic correlations, a manifestation of the elusive Nagaoka effect 15,16 . We study the evolution of these correlations with interactions and doping, and use higher-order correlation functions to further elucidate the relative contributions of superexchange and kinetic mechanisms. The robustness of itinerant spin polarons at high temperature paves the way for exploring potential mechanisms for hole pairing and superconductivity in frustrated systems 10,11 . Furthermore, our work provides microscopic insights into related phenomena in triangular-lattice moiré materials [17][18][19][20] .

One of the key questions in quantum condensed-matter physics is how doped Mott insulators give rise to exotic metallic and superconducting phases. Understanding this problem is crucial for explaining the emergence of the unusual physical properties of many families of strongly correlated electron systems, including the high-critical-temperature (T c ) cuprates 21 , organic charge transfer salts 22 and moiré materials 23,24 . An important aspect of this problem is the interplay between spin order and the quantum dynamics of mobile dopants. So far, most studies have focused on Mott insulators on a square lattice where the motion of charge carriers disturbs spin correlations, resulting in an adversarial relationship between doping and spin order 1,[3][4][5][6][7][8] . This explains why many theoretical studies of doped high-T c cuprates are usually done from the perspective of Mott states in which spin order has been suppressed by fluctuations [25][26][27] .

Recently, experiments on moiré materials have provided a strong motivation for understanding doped Mott insulators in triangular lattices 24 . Here we explore this problem microscopically using a cold-atom triangular Fermi-Hubbard system 28,29 . One surprise of our experiments is that in contrast to square-lattice systems, there is a symbiotic relation between mobile holes and antiferromagnetism. This manifests in the formation of antiferromagnetic (AFM) itinerant spin polarons in the hole-doped system, which we directly image by measuring spin correlations around mobile holes. In striking contrast, we find that particle doping favours the formation of ferromagnetic (FM) polarons similar to those discussed previously for the square-lattice Fermi-Hubbard model 15,16 .

Some of the most important implications of our results are for systems in which the local interaction U is much larger than the single electron tunnelling t, in which case the magnetic superexchange J is strongly suppressed. Indeed, this is the relevant regime for most moiré systems (neglecting nearest-neighbour interactions). Intuition based on earlier studies would suggest that at temperatures higher than the superexchange scale, the regime we explore here, one can not expect coherent propagation of quasiparticles 1 . Our results demonstrate that this does not have to be the case in triangular lattices. Formation of polarons around mobile dopants facilitates their propagation and makes their dynamics more coherent. This robustness of the quasiparticle can also be understood as the result of effective magnetic interactions with energy scale t induced by the motion of dopants in the frustrated system 9 .

including the triangular one. To release this frustration and lower the kinetic energy of the dopants, the system develops magnetic correlations. The resulting magnetism, known as kinetic magnetism, is closely related to that studied by Nagaoka is his seminal work 15 , but is more robust in that it is predicted to survive for finite interactions and doping.

Indeed, signatures of kinetic magnetism above half-filling have been observed recently in doped van der Waals heterostructures through measurements of the spin susceptibility 17,19 , while separate measurements of magnetization plateaus attributed to kinetic effects below half-filling have also been measured in these materials 20 . Our experimental results provide a microscopic picture underlying these observations. More broadly, our results motivate studying the properties of doped Mott insulating states in triangular lattices, including superconductivity, from the perspective of self-organization of itinerant spin polarons 11,27,34,35 .

Our system consists of a two-dimensional degenerate gas of 6 Li that is an equal mixture of two spin species corresponding to the first ( ↑⟩ | ) and third ( ↓⟩ | ) lowest hyperfine states of the atom. The gas is loaded adiabatically into an optical lattice realizing the triangular-lattice Hubbard model

where  c iσ † (c iσ  ) creates (destroys) a fermion of spin σ at lattice site i, the number operator

measures site occupation and ⟨i, j⟩ denotes nearest-neighbour sites. In the model, particles hop with t > 0. With this sign of the tunnelling, a particle in an empty lattice can lower its energy by delocalizing on each lattice bond in a symmetric spatial orbital (Fig. 1a). The corresponding band structure is particle-hole asymmetric, and the particle attains its minimal energy of -6t at zero quasi-momentum. Kinetic frustration can be understood by considering the opposite scenario of a single hole moving in a spin-polarized background. In this case, the Hamiltonian is better expressed in terms of hole operators with

Crucially, the change in the sign of the tunnelling resulting from the anticommutation of fermionic hole operators in the Hamiltonian favours antisymmetric spatial orbitals for the hole on each bond. Indeed, this is manifested in the band structure of the hole, which is mirrored about zero energy relative to the particle. The hole kinetic energy is thus minimized at a value of -3t, larger than in the unfrustrated system.

A simplified picture explaining the emergence of the itinerant spin polaron in the doped interacting system can be obtained by considering a triangular plaquette with two fermions. In the limit of strong interactions, double occupancies are energetically forbidden and the motion of the hole on a closed loop on the plaquette will exchange the two spins. In the spin singlet sector, this produces a spin Berry phase of π, whereas the phase is zero in the spin-symmetric triplet sector (Fig. 1b) 36 . The phase acquired by the hole in the singlet sector returns the tunnelling to a positive value, thereby releasing the kinetic frustration and allowing the hole to reach a lower ground-state energy.

The resulting object, a singlet bond bound to a hole with a binding energy of order t, is predicted to persist in the many-body setting for light hole doping (Fig. 1c). This corresponds to a polaron with AFM spin correlations in the vicinity of a hole. The situation is reversed for particle doping, favouring FM correlations in the vicinity of a doublon. We directly detect the itinerant spin polaron in our system using a connected three-point correlation function, which probes the spin correlations in the environment of a hole or doublon. Such correlators have been previously used to identify magnetic polarons in the square lattice, although, in that case, the mechanism that leads to the formation of the polaron is different and the binding energy is on the superexchange scale 6,8,37 .

To realize a lattice with triangular connectivity 28,29,[38][39][40][41][42] , we superimpose two non-interfering lattices, a strong one-dimensional (1D) optical lattice with spacing a = 532 nm and depth V 532 = 6.7(2) E R,532 and a weak square optical lattice with larger spacing a = 752 nm and depth V 752 = 2.9(1) E R,752 , where one recoil energy is defined as E R,a ≡ h 2 /8ma 2 with m the mass of the atom and h is Planck's constant 41 (Fig. 1d). The frequency detuning between the two lattices is used to tune their relative alignment to obtain a triangular geometry. Their relative depths are chosen to produce an isotropic triangular lattice by equalizing the tunnelling strength along the original square lattices axes and one diagonal. The gas is prepared at a magnetic field near a Feshbach resonance at 690 G, allowing us to freely tune the scattering length. In this way, we tune the coupling strength U/t to explore the evolution of the correlations from the metallic to the Mott insulating regime.

We use a quantum gas microscope to measure site-resolved correlations associated with the polaron in the many-body system 43 . We further implement a bilayer imaging technique [44][45][46][47] , wherein a magnetic field gradient is first used to separate the two spin states into different layers before imaging them simultaneously (Fig. 1e). From the reconstructed images, we can calculate arbitrary n-point correlation functions involving both spin and density operators averaged over experimental cycles. In the strongly interacting regime, the atoms order in a Mott insulator and show short-range 120° spiral AFM correlations that have been observed in previous experiments 28,29 . We use the two-point spin correlations for thermometry by comparison with determinant quantum Monte Carlo (DQMC) calculations 48 (Methods). The typical peak density of the clouds in the lattice is n = n ↑ + n ↓ = 1.2, allowing us to study a range of dopings δ = n - 1 on either side of half-filling of the Hubbard system in each experimental snapshot due to the harmonic confinement of the lattice beams.

To detect the polaron, we evaluate connected three-point chargespin-spin correlation functions. For a hole dopant, the relevant correlation function,

is computed as:

where we have assumed a spin-balanced system S ⟨ ⟩= 0

For the range of dopings we consider, the two-point spin correlator is always negative due to the dominant superexchange aniferromagnetism. The second term of the connected three-point correlators defined in equation ( 3) removes any uncorrelated charge-spin-spin signal associated with this background AFM signal.

We start by studying the doping dependence of C (3) above and below half-filling, shown in Fig. 2a.

) in the vicinity of holes (doublons) are shown for a strongly interacting sample at U/t =11.8 (4). For δ < 0, C h (3) shows negative correlations for = (1, 0),

indicating an enhancement of AFM order in the immediate vicinity of a hole. In contrast, for δ > 0, the doublon correlation function C d (3) is positive for the same bond, revealing a preference for FM order around doublon dopants. We then examine correlations out to further distances to explore the structure of the polaron at a few different dopings in Fig. 2b. In particular, we emphasize the structure of the correlations in the limit of vanishing doping in the right column of this panel. This doping regime is the closest to capturing the behaviour associated with the idealized case of a single dopant 9,15 where polaron-polaron interactions are absent. Around a hole (Fig. 2b,c), we find that although at the shortest distance C h (3) is negative, the next furthest ring (distances 1.32 and 1.5) shows positive correlations. The structure of the AFM hole polaron can be understood using a picture

e a b c d -0.8 -0.6 -0.4 -0.2 0 0.2 -12 -10 -8 -6 -4 -2 0 2 -0.8 -0.6 -0.4 -0.2 0 0.2 -2 0 2 4 6 = -0.10 = -0.04 = 0.15 = 0.02 | | → 0 Doping Correlation distance = -0.04 = 0.02 -8 -6 -4 -2 0 2 4 6 8 -6 -4 -2 0 2 4 6 0.5 1.0 1.5 2.0 2.5 3.0 Correlation distance -1 0 1 2 0.5 1.0 1.5 2.0 2.5 3.0 -3 -2 -1 0 1 2 -0.8 -0.6 -0.4 -0.2 0 0.2 -0.05 0 0.05 0.10 Doping Doublon C (3) correlation (×10 -3

) Hole C (3) correlation (×10 -3 ) C (3) correlation (×10 -3 ) C (3) correlation (×10 -4 ) Doublon C (3) Doublon C (3) Hole C (3) Hole C (3) Doublon C (3) correlation (×10 -3 ) Hole C (3) correlation (×10 -3 versus distance at dopings of δ = -0.04(2) and δ = 0.02(2), respectively. Correlation distance is defined as the distance from the dopant to the bond midpoint, shown as the inset black arrow. Bonds are unity length, implying the closest correlation is at a distance 3 /2. d, Dopant propagation coherently re-orders the surrounding 120° order, resulting in an alternating pattern of correlations, as imaged in the experiment. e, Nearest-neighbour three-point correlators normalized by the dopant density, where red and green points are measured hole and doublon correlators. DQMC theory is shown for the interacting systems (grey bands) and in the non-interacting limit (dotted red and green lines). Error bars represent 1 s.e.m.

of a mobile hole dopant that coherently modifies the surrounding 120° Néel order to facilitate lowering the kinetic energy, as illustrated in Fig. 2d. Strikingly, this is a different structure than on the doublon side, where predominantly positive correlations exist up a distance d ≈ 1.8 away from the doublon. This indicates an energetic preference towards a locally FM environment, as for a particle dopant the motion is unfrustrated in a background of polarized spins. We interpret these short-range polaronic correlations we observe as the precursors to Haerter-Shastry AFM and Nagaoka FM expected at lower temperatures on the hole-and particle-doped sides, respectively. We also note from Fig. 2a the approximate linear dependence of the correlators with relevant dopant density δ for |δ| ≲ 0.1. This indicates that in this regime, the description of the system in terms of weakly interacting polarons is valid. Polaron interactions become important for larger dopant densities. These observations motivate introducing a normalized version of the correlators by dividing out the relevant dopant density δ, where

(Fig. 2e). While the non-normalized correlations close to zero doping appear reduced in magnitude, the normalized spin correlation emphasizes the fact that the spin correlations per dopant are in fact strongest close to half-filling.

The itinerant spin polaron picture we have presented so far is in the regime of strong interactions, but it is also interesting to explore how the three-point correlations evolve with U/t. Figure 3 shows these correlations in the metallic (U/t = 4.4), Mott insulating (U/t = 11.8) and intermediate (U/t = 8.0) regimes in the temperature range T/t ≈ 0.7-0.9. Surprisingly, many of the qualitative features of the correlations are similar, including the minimum in the AFM correlations around a hole at δ ≈ -0.3. This can again be understood from a single plaquette in the alternative limit of vanishing interactions, which predicts correlations of the same sign as the itinerant spin polaron 49 (Methods). For all interactions, the measured correlations shows reasonable agreement with DQMC calculations with a systematic deviation in C d (3) and C h (3) for larger fillings, possibly owing to an increase in reconstruction errors (Methods). The correlations differ significantly from those expected for the non-interacting gas, especially for the two stronger interactions.

As U/t increases, the onset of the correlation moves closer to half-filling as contributions from virtual doublon-hole fluctuations are increasingly suppressed. The characteristic linear growth of the correlations, expected in the polaronic regime and observed for U/t = 11.8, is, however, present for only the strongest interactions. Additional evidence for an experimental observation of the itinerant spin polaron at the largest interaction strength comes by combining the observed three-point correlators with a measurement of the singles fraction n s = n - n d at half-filling, which ensures that the system is in the strongly interacting regime. For U/t = 4.4, U/t = 8.0 and U/t = 11.8, this is 0.71(1), 0.85(1) and 0.93(1), respectively.

As superexchange-induced AFM correlations are present in the system for any finite interaction strength, it is illuminating to quantify their strength compared with kinetic magnetism around dopants. This can be done using four-point correlation functions. For any two nearestneighbour sites j and k, there are two additional sites i and l (which we call conditioning sites) coupled to both of them (Fig. 4, inset). As a dopant on either conditioning site may affect the correlation strength between j and k, a four-point correlation function is required to determine the influence of the background on the j-k bond. We define a conditional four-point correlator as the spin correlator on the shared bond, conditioned on the occupancy observables n n , a b   on sites i and l

i j z k z l i l ab (4) a b a b

where the labels a, b ∈ {h, s}. Figure 4 shows the four-point correlator at U/t = 11.8 for three different occupancies of the conditioning sites. Below half-filling, we find that

, directly indicating the enhancement of AFM correlations on a given bond by the presence of holes on the conditioning sites. Kinetic magnetism, therefore, strengthens the existing AFM correlations below zero doping that . Data are shown for U/t = 4.4(1), T/t = 0.68(2) (a), U/t = 8.0(2), T/t = 0.84(3) (b) and U/t = 11.8(4), T/t = 0.94(4) (c). For all interactions, we observe the largest negative C h (3) correlator at a doping of around -0.3. For increasing U/t, the C h (3) and C d (3) correlators become more linear near zero doping, indicating a region where there are weakly interacting polarons. The blue dashed line in c, which is fit to the DQMC in the doping range -0.09 < δ < 0.06, illustrates this region. DQMC theory is shown for the interacting systems (grey bands) and in the non-interacting limit (dotted red and green lines). Error bars represent 1 s.e.m. arise due to superexchange, which is the only mechanism at play on the level of the four-site plaquette when the conditioning sites are singly occupied. Above zero doping, where holes are due to virtual fluctuations, there is no enhancement of AFM correlations from kinetic magnetism and all three correlators are similar in value (Extended Data Table 1).

In this work, we have directly imaged itinerant spin polarons in a triangular Hubbard system by measuring three-and four-point correlation functions. We have characterized their evolution with doping and interactions, and compared the strength of correlations induced by superexchange and kinetic effects. In future work, it would be interesting to study the polaron spectroscopically 50 , which would allow direct characterization of its binding energy as well as its dispersion and effective mass. While the polarons we have focused on here are part of the physics of the doped triangular Hubbard model at high temperatures, pushing to lower temperatures would shed light on its rich ground-state phase diagram. Theoretical work suggests that this may include magnetically ordered phases as well as a quantum spin liquid with fractionalized excitations at intermediate interactions 51,52 .

Higher-order connected correlations may be useful in identifying more complex multi-particle bound states in the system 11 , which can lead to hole-pairing mechanisms and superconductivity at high temperatures 10,11,51,[53][54][55] .

Note added in proof: During completion of the paper, we became aware of related work studying three-point correlations in a triangular Fermi-Hubbard system 56 .

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-024-07356-6. Conditional four-point correlators C (4) show the spin-spin correlator on the j-k bond when surrounded by two holes (red), one hole and one singly occupied site (green), and two singly occupied sites (blue). Below zero doping, bonds have increased AFM spin correlations with increasing neighbours that are holes. DQMC theory is shown as coloured bands. Error bars represent 1 s.e.m.

We use a degenerate mixture of hyperfine states 1⟩ | and 3⟩ | , where i⟩ | represents the ith lowest energy level of the ground hyperfine manifold of atomic 6 Li, to simulate the two-component Fermi-Hubbard model on a triangular lattice. State preparation of a degenerate Fermi gas before loading the science lattice largely proceeds as detailed in previous work 57 . After the final stage of evaporation, we are left with a spin-balanced sample of approximately 400 atoms in each spin state. At this stage, the atoms are confined in a single layer of an accordion lattice, created with 532-nm light, with spacing a z = 3.6(3) μm and trap frequency ω z = 2π × 16.4(2) kHz in the vertical (z) direction. The combined in-plane two-dimensional lattices (see below) are then ramped to their final depths following a cubic spline trajectory in 100 ms. An additional 1,070-nm optical dipole trap propagating along z with waist w 0 = 100 μm is used to provide variable confinement in the x-y plane in the final science configuration. In particular, for strongly interacting samples, the reduced compressiblity necessitates greater confinement to achieve comparable densities. We worked at magnetic fields ranging from 587(1) G to 612(1) G, where the scattering length varies between 330 (15) and 945 (30) Bohr radii, respectively.

The triangular lattice is formed as in ref. 41 by combining two noninterfering lattices of different polarizations and detunings (Extended Data Fig. 1, inset). Both lattices are created using light of wavelength 1,064 nm. The first is a square lattice with a spacing of 752 nm, created by retroreflecting a single vertically polarized laser beam in a bowtie geometry. The depth of this lattice is calibrated using amplitude modulation spectroscopy. Both losses of laser power as the beam traverses the vacuum chamber and and non-orthogonal beam alignments can cause a significant tunnelling imbalance along the axes of the square lattice. In our system, we specifically tune the angle between the lattice beams to 90.7(1)° (measured using atomic fluorescence images), which approximately cancels the imbalance due to power losses. AFM correlations along the two axes of the square lattice show a systematic difference of 4(3)%, indicating a difference in the tunnellings of approximately 1(1)%.

The second lattice is a 1D optical lattice with a spacing of 532 nm and wavevector aligned with a diagonal of the 752-nm square lattice. The light for this lattice is horizontally polarized and detuned by about 330 MHz with respect to the square lattice, preventing any electric field interference. Both lattices share a common retroreflecting mirror, avoiding the need for active phase stabilization as in other schemes 29,58 .

The frequency detuning between the two lattices introduces a relative spatial phase between the two potentials at the atoms which is given by ϕ = 4πLΔ/c, where Δ is the relative detuning, L is the distance from the atoms to the retroreflecting mirror and c is the speed of light. The triangular-lattice configuration is obtained for the case of constructive interference, that is, ϕ = 0, which we calibrate using in situ measurements. The superlattice depth is set to a weak value (V 532 = 0.49(1) E R,532 ) relative to the dominant square lattice (V 752 = 40.2(3) E R,752 ) and then modulated at the frequency of the square-lattice p-band resonance. As this is an odd-parity transition, excitation should be maximized when ϕ = π/2 or ϕ = 3π/2, which induces a sloshing motion. The two corresponding prominent resonance peaks versus superlattice detuning at a constant modulation frequency are shown in Extended Data Fig. 1.

We perform two identical measurements separated by 1 week (red and blue data) to assess the long-term stability of the set-up. The agreement of the resonance peaks between the two datasets is at or below the uncertainty (about 1 MHz) of the spectroscopic measurement, indicating phase stability at or below 0.02π radians. Explicit band-structure calculations show that such a phase drift results in a negligible change of the tunnelling values of less than 0.2 Hz on top of a tunnelling strength of 400 Hz.

We note that the ϕ = 0 (triangular) and ϕ = π (honeycomb) conditions are indistinguishable from spectroscopic measurements alone as they both produce an even-parity drive. To distinguish these two phases, a dense Mott insulator is prepared and subsequently allowed to expand in the combined superlattice potential with V 752 = 42.0(3) E R,752 and V 532 = 3.7(1) E R,532 for 1 second. The constructive interference in the triangular lattice results in a deeper potential compared with the destructive interference present in the honeycomb lattice, resulting in a much denser cloud following the same period of expansion. The combination of these in situ measurements uniquely determines the superlattice phase.

In principle, each lattice may be independently calibrated to give a full reconstruction of the potential in the plane of the atoms. However, owing to limited power available in the 1D lattice, independent calibration with modulation spectroscopy is difficult as the band transitions are not truly resolved. Instead, precise knowledge of the depth of the square lattice, V 752 , combined with knowledge of the relative tunnellings (obtained from correlation maps of the system), can be used to obtain the 1D lattice depth. We empirically find the depth of the 1D lattice that equalizes t x , t y with the diagonal tunnelling t d . This is done by experimentally equalizing the nearest-neighbour two-point spin correlations in the triangular lattice. The depth of the square lattice used in the experiment is measured to be V 752 = 2.9(1) E R,752 . At the point where we obtain an isotropic triangular-lattice connectivity, we infer the depth of the 1D lattice using the computed band structure to be V 532 = 6.7(2) E R,532 . This corresponds to absolute tunnelling strengths of t x = t y = t d = h × 400(20) Hz.

Full spin-charge readout in a bilayer imaging scheme Simultaneous imaging of charge and spin information is performed using a bilayer imaging scheme using Raman sideband cooling, similar to the method discussed in a previous publication 47 . Minor differences from the previous scheme are discussed here.

Imaging consists of four steps: (1) Tunnelling is quenched by deepening the 2D lattice depth to 56.3(4) E R,752 in 170 μs. The axial confinement lattice is turned off in 20 ms. Atoms in the ground hyperfine state 1⟩ | are transferred to hyperfine state |2⟩ using a radiofrequency Landau-Zener sweep lasting 50 ms. In this state, the magnetic moment is of opposite sign to state |3⟩. The magnetic Feshbach field is turned off in 10 ms.

(2) A magnetic field gradient of 336 G cm -1 is applied to separate the two spin components 2⟩ | and 3⟩ | along the z axis. For this step, the two-dimensional lattice is increased to 160(1) E R,752 .

(3) Each spin component is trapped by a light sheet potential with vertical waist w z ≈ 5 μm, which is turned on in 20 ms. The two potentials are moved further apart in the z direction using a minimum-jerk trajectory in 4 ms to a final separation of 16 μm. (4) Both layers are imaged simultaneously using Raman sideband cooling over a 2-s duration. The fluorescence photons are collected by a microscope objective and focused on two separate areas of a complementary metal-oxide-semiconductor camera.

During the imaging procedure described above, various types of error can accumulate that will affect imaging fidelities. Radiofrequency spin-flip fidelities exceed 99%, and these errors have a negligible role compared with other infidelities. Transport fidelity encompasses various errors that occur during the vertical motion of the atoms. First, we observe some atom loss that could be due to off-resonant scattering or background gas collisions. A surviving atom may hop to other sites of the same layer, which disturbs magnetic correlations or, if the final site is already occupied, leads to atom loss due to parity imaging. Finally, atoms may be transported into the wrong layer, so that they are assigned to the wrong spin state.

These effects are difficult to isolate and characterize independently. We instead benchmark a related quantity: by preparing an almost unity-filled Mott insulator, we observe the proportion of singly occupied sites with and without the transport step. We prepare Mott insulating states with 97.1(4)% singles fraction, verified by imaging both spin states in a single layer. Any sites with zero atoms or two atoms appear dark from parity imaging. The transport step is tested by adding the Stern-Gerlach and optical transport (steps 2 and 3), and then reversing those steps to transport both spin states back into a single layer. Then, the visible singles fraction drops to 95.4(5)%. We assume that transport hopping errors populate randomly distributed sites and are irreversible. A hopping event will create a hole and a double occupancy. Therefore, this test indicates a transport infidelity of at most 0.9(3)%.

In addition, errors may accrue during the Raman sideband cooling, appearing as loss (3.9%), interlayer hopping (0.5%) and intralayer hopping (negligible).

Finally, errors can be introduced during image processing when we digitize the images into an occupancy matrix. Compared with the bilayer readout of a sparse tweezer array of fewer than 50 atoms 47 , our current bilayer imaging scheme must reliably reconstruct atomic distributions of hundreds of atoms with high filling. Each layer adds an out-of-focus background on the image of the opposite layer, decreasing our signal to noise. We choose a 2-s Raman imaging time as a compromise between increasing the ratio between the desired signal and the background layer noise and minimizing hopping and loss errors. The problem of the out-of-focus background ultimately limits the peak densities that we can reliably probe. Empirically, we find that beyond dopings of 0.2, distinguishing between empty and occupied sites becomes difficult. We include a representative image of the occupation histograms for both imaging layers and the corresponding Gaussian fits to the zero and single atom peaks in Extended Data Fig. 2.

The experimental correlation functions presented in the text are computed as the fully connected three-point correlation function of a three-observable operator:

In particular, we do not a priori assume a perfectly spin-balanced gas, compared with the simplified definition in equation ( 3).

As we equalize the tunnellings in all three directions and the same lattice depths are used for all datasets, we average over all 120° and reflection symmetric higher-order correlators for plots versus doping. In addition, all DQMC calculations are done with the assumption that t x = t y = t d . For completeness, in Extended Data Figs. 5 and 6, the same data as in Fig. 2a are shown without symmetrization. There appear to be no major systematic differences between correlators of different orientations. This feature holds for all datasets and three-and four-point correlators shown in the main text.

We use a bootstrapping analysis technique to obtain vertical error bars for the U/t = 11.8 (1,146 experimental runs), U/t = 8.0 (535 experimental runs) and U/t = 4.4 (360 experimental runs) datasets for all correlators. The experimental runs are randomly separated into 80 groups, and the relevant three-and four-point correlations are calculated for each group. We sample from these 80 groups with replacement 10,000 times to obtain 10,000 bootsamples. We average over all lattice symmetries before taking the standard deviation of the bootsamples.

To create strongly interacting samples with a high central density, additional radial confinement is provided by an external dipole trapping beam at 1,070 nm propagating approximately orthogonal to the atom plane. For the most strongly interacting datasets at U/t ≈ 12, the radial trap frequency ω r is approximately 2π × 370 Hz. Although the gradient near the centre of the trap (regions of highest density) remains small, gradients away from the centre of the lattice have the potential to impact resonant tunnelling and affect correlations, particularly long-range and multi-point correlations. We empirically test for such by comparing two datasets with different global chemical potentials (different total atom number) but otherwise identical science parameters 7 . A global chemical potential shift will displace a bin of given density to a different radial position and hence cause it to sample a different spatial gradient. Disagreement between the two sets, particularly at low densities where the gradient is largest, would therefore indicate an effect due to the spatial gradient. In Extended Data Fig. 3, we show measured three-and four-point correlations versus doping for two datasets with different total atom number. We find no significant systematic deviations between these two datasets, from which we conclude that the gradient at the level present in the experiment does not affect measured correlation functions within experimental error bars.

The four-point correlations shown in Fig. 4 are for the strongest interacting sample. In Extended Data Fig. 4, we compare these measurements with the more weakly interacting dataset at U/t = 8 to probe the evolution with interaction strength. Close to half-filling, we do not measure a significant difference in the strongest correlation C hh (4) , whereas the two correlations with singlon nearest neighbours C hs (4) and C ss (4) are reduced in magnitude with increasing interactions. This is consistent with kinetic magnetism being enhanced relative to superexchange. Nonetheless, the variable temperature between datasets makes quantitative comparison difficult.

In this section, we present the ground-state wavefunctions and correlation functions on a three-site triangular plaquette. This toy example shows many of the important features that are present in the larger lattice system and is therefore instructive to consider in detail. The ground state of the triangular plaquette above zero doping with N = 4 spin has energy E N=4 = -2t + U in all magnetization sectors S z ∈ {-1, 0, 1} for U > 0, with ground-state wavefunctions given by 49 : | -0.6 -0.4 -0.2 0 0.2 -12 -8 -4 0 -0.8 -0.6 -0.4 -0.2 0 0.2 -10 -5 0 5 x10 -3 C (3) Correlation C (4) Correlation Doping Doping hole doublon x10 -2

C (4) ss C (4) hs C (4) hh a b

Extended Data Fig. 3 | Multi-point Correlations at Different Global Chemical Potentials. a, C h (3) (red data) and C d (3) (green data) evaluated at the bond closest to the dopant, and b, C (4) for two different datasets with distinct global chemical potentials. As these two datasets track each other well, particularly at low filling, we conclude spatial gradients do not appreciably suppress resonant tunneling to affect measured correlations. Filled data points have mean atom number 799 (35) while empty data points have mean atom number 622 (27).