Formation of individual stripes in a mixed-dimensional cold-atom Fermi–Hubbard system

The phase diagram of high-temperature superconducting materials has so far eluded full understanding despite 40 years of extensive theoretical and experimental studies^4,5,6. Especially, the exact nature of the intricate relationship between superconducting pairs and stripes remains an open question^1,2,3,4,5,6. Although experimentally both phenomena may be found in close proximity, numerical studies have long been investigating whether stripes precede, compete or coexist with superconductivity^15,16,17. Recently, by the discovery of superconductivity in bilayer nickelates⁷, a new class of superconducting materials has been found. These materials are conjectured to be of a mixD nature, in which the dynamics of charge is restricted to two dimensions, whereas the spins order in a bilayer system. Here we study mixD systems, which are predicted to show an enhanced version of stripe order^8,9,10.

The repulsive, two-dimensional (2d), spin-1/2 Fermi–Hubbard model and its natural extensions are widely assumed to provide minimal models to capture the physics of these doped antiferromagnets (AFMs). Ultracold atoms in optical lattices provide natural implementations of the Hubbard model with a high degree of control over system parameters and dimensionality^18,19. Although solid-state experiments mostly focus on spectroscopic and dynamical response measurements, quantum simulation, especially with single-site resolution, opened up access to new sets of microscopic observables and correlation functions^20,21,22,23. Previous studies in 2d systems found AFM correlations^{24,25,26,27,28,29}, investigated the effect of doping on the spin order^{19,30,31,32,33,34,35} and observed pairing of dopants in tailored ladder systems³⁶.

Here we present the first, to our knowledge, observation of hole attraction beyond nearest-neighbouring sites in a repulsive, 2d Hubbard system with mixD coupling. Using higher-order charge and spin correlators, we find signatures of extended charge structures, which we identify as individual stripes.

Stripe phases, characterized by charge density waves in combination with incommensurate AFM order, have been found in measurements on solids^1,2,3,37,38 as well as numerical studies^{12,13,14,17,39,40,41,42,43}. These stripes form out of individual dopants of an AFM background, a process governed by the competition between the kinetic energy favouring delocalization and the magnetic energy of the AFM spin order, which is disrupted by dopant motion. Consequently, the energy scale at which stripe order is expected to occur is only around 5% of the tunnelling energy⁴³ (about 10% of the superexchange energy), placing it out of reach for state-of-the-art quantum simulators. In particular, for temperatures around the superexchange energy, the effective repulsion resulting from the fermionic nature of the holes (Pauli blocking) disfavours tightly bound hole pairs and extended structures such as stripes but favours the formation of magnetic polarons³⁰ (see Fig. 1a).

a, Illustration of the isotropic 2d Fermi–Hubbard model. Holes delocalize within small regions and disturb their respective spin background, forming magnetic polarons. The overall hole density is uniform and holes repel each other owing to their fermionic statistics at experimentally accessible temperatures of k_BT ≈ J. There are no ___domain walls in the spin order. b, By raising the potential on every other lattice site along y by Δ, we suppress tunnelling along this direction, thus removing the Pauli repulsion between holes, while preserving the superexchange coupling J_y. At low temperatures, the holes form collective structures, which also result in a ___domain wall in the AFM correlations of the system, indicated by the AFM parity change across the stripe. c, A single raw experimental shot of spin-up (red) and spin-down (blue) atoms and doubly occupied sites (purple), as well as its reconstructed spin and charge distribution with the main system being inside the black circle, surrounded by a low-density reservoir (see Supplementary Information). The green box indicates a stripe-like structure.

Full size image

MixD systems allow to tilt the balance towards collective charge and spin ordering by restricting the hole motion to one dimension, thus reducing the kinetic energy while keeping spin couplings 2d. This leads to an increase in the characteristic energy scales of collective effects as kinetic and magnetic terms in the Hamiltonian are less frustrated, lifting these effects to experimentally accessible regimes⁴⁴. In nickelates, which are mixD bilayer systems, this causes high critical temperatures for superconductivity^7,8,9,10. Similarly, in ladder systems, in which the spin order is dominated by rung singlets, numerics¹¹ and experiments³⁶ confirmed strong pairing in mixed dimensions. Here we apply this mixD approach to the unexplored 2d system, in which the spin sector is not gapped and classical simulations are limited to very small system sizes, especially at the relevant intermediate temperature scales. By adding a sufficiently large potential offset to every other chain within the lattice, we remove nearest-neighbour hopping along the perpendicular direction while increasing spin couplings^45,46. This biases hole attraction and stripe formation along the direction perpendicular to the hole motion because Pauli repulsion is suppressed in this direction. For the same reason, fully filled stripes are favoured but the key concept of charge ordering associated with the stripe phase is retained^42,43 (see Fig. 1b). This allows us to study the poorly understood temperature regime around the superexchange energy in which individual stripes are expected to form.

In the experiment, we realize the spin-1/2 Fermi–Hubbard model by using ⁶Li atoms in an optical superlattice with a homogeneous, circular system of about 110 sites surrounded by a low-density reservoir (see Fig. 1c). In the limit of strong on-site interactions U, the essential physics of the system can be captured by the t–J Hamiltonian using projections \(\widehat{{\mathcal{P}}}\) onto singly occupied sites,

with tunnel couplings t_ij ∈ {t_x, t_y}, spin exchange couplings J_ij ∈ {J_x, J_y}, spin-up/spin-down fermionic creation (annihilation) operators on site i, \({\widehat{c}}_{{\bf{i}},\uparrow /\downarrow }^{\dagger }\) (\({\widehat{c}}_{{\bf{i}},\uparrow /\downarrow }\)) and on-site spin (density) operators \({\widehat{{\bf{S}}}}_{{\bf{i}}}\) (\({\widehat{n}}_{{\bf{i}}}\)). This model suffers from the fermion sign problem, making it numerically challenging to tackle even in the mixD regime⁴⁷.

Here we work at U/t_x = 27(2), J_y/t_x = 0.6(2), J_x/t_x = 0.15(3) and a filling of n ≈ 0.7–0.9 (hole doping δ = 1 − n) with a temperature of k_BT/t_x = 0.3(1) (see Supplementary Information). We make use of an optical superlattice along y to controllably detune neighbouring sites by \(\varDelta =0.65(5)U\gg {t}_{x},{t}_{y}^{{\prime} }\), thus effectively disabling nearest-neighbour tunnelling along y (t_y ≈ 0), and leading to a spin coupling \({J}_{y}=2{({t}_{y}^{{\prime} })}^{2}\left(\frac{1}{U-\varDelta }+\frac{1}{U+\varDelta }\right)\), in which \({t}_{y}^{{\prime} }\) is the tunnel coupling in the 2d system without potential offsets. Owing to the staggered superlattice potential, there is also a second-order next-nearest-neighbour hopping term along y, which reintroduces a weak Pauli repulsion at distance d_y = 2. This term, however, is smaller than J_y, such that it is still expected to be favourable for stripes to form (see Supplementary Information for more details on preparation and subdominant couplings).

To reveal the charge order within the system, we evaluate the connected, normalized two-point hole–hole correlator

with hole density operator \({\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}\) at position i and normalization \({{\mathcal{N}}}_{{\bf{d}}}\) the number of bonds with distance d. Owing to the finite size and particle number fluctuations in our system, there is a global, doping-dependent offset o_δ ∈ [−0.06, −0.03] on this correlator that we subtract (see Supplementary Information). A positive (negative) value of this correlator indicates attraction (repulsion) between holes at distance d.

We consider hole correlations in a mixD system with a doping of δ = 0.18 in Fig. 2a,b. We observe a positive nearest-neighbour correlation along y, whereas along x, we find antibunching caused by the Pauli repulsion of the holes (see Fig. 2a). Furthermore, at larger distances d_y > 1, there are positive correlations, which indicates that, instead of merely forming isolated, nearest-neighbour hole pairs, there is a finite probability that vertically aligned hole structures are extended through the system. Furthermore, there are significant correlations along the diagonals at d = (1, 1), which we interpret as signs of charge fluctuations along x owing to the finite coupling t_x. The correlations at d_y = 2 are slightly suppressed, which we attribute to next-nearest-neighbour hopping (see Supplementary Information). Finally, the positive signal at d_x = ±5 may be related to the presence of a second, vertically aligned charge structure in the system.

a, Hole–hole correlations in a mixD system revealing bunching (attraction) along y at distances d_y ≥ 1 and antibunching (repulsion) along x with δ = 0.18. The symmetrization is indicated by the dashed lines. b, A cut along y (x) is shown in dark (light) green for mixD systems, with the inset showing the equivalent data for standard 2d systems. c, The dependency of the mixD correlator \({g}_{{\rm{hh}}}^{(2)}\) at distance d = (0, 1), (1, 1), (0, 2), (0, 3) on doping is plotted in red, purple, blue and grey with a doping binning of  ±0.009. In a doping region around 0.2, the correlators for distances d_y > 1 are positive, indicating longer range charge correlations. Error bars are estimated using bootstrapping. d, Results for the renormalized correlator (see Supplementary Information) from DMRG calculations for a system size L_x × L_y = 8 × 3 as a function of doping for k_BT/t_x = 0.41.

Source data

Full size image

By considering one-dimensional (1d) cuts along y and x (Fig. 2b), we corroborate the bunching (antibunching) along y (x) through the system in the mixD setting. By contrast, for the standard 2d system (Δ = 0; Fig. 2b, inset), there is antibunching along both directions. Compared with this, the anticorrelations along x are enhanced in mixD by removing t_y as a result of the absent competition between anticorrelations along x and y.

To identify whether an ideal doping level for the emergence of individual stripes exists in our mixD system, we bin our data by doping and calculate \({g}_{{\rm{hh}}}^{(2)}\) per bin (see Fig. 2c). Both the nearest-neighbour and diagonal correlations decrease with doping above δ ≈ 0.15, indicative of the decrease in pairing probability with doping and compatible with a reduction of the spin correlations responsible for the binding. For d = (0, 2), (0, 3), there is a non-trivial dependence on doping with positive correlator values starting at δ = 0.17. This is indicative of a possible transition from the formation of individual pairs to extended stripe-like structures^41,42,43.

We compare the correlations along y to density matrix renormalization group (DMRG) calculations of equation (1) on 8 × 3 sites, J_y/t_x = 0.5, J_x/t_x = 0.15, k_BT/t_x = 0.41 as a function of doping in Fig. 2d (see also Supplementary Information for normalization). Although the trend in correlations at d = (0, 1) and d = (1, 1) is qualitatively comparable with the experimental data, the same cannot be stated about correlations at d = (0, 2). We attribute these differences to the strong finite size limitations in the DMRG along y as well as the simulation of the t–J model instead of the Hubbard model. Further differences could arise owing to the presence of the aforementioned second-order hopping process that introduces further repulsion between holes, as well as the statistical distribution of holes between different chains in the experiment, whereas calculations feature balanced hole numbers.

The connected two-point correlator \({g}_{{\rm{hh}}}^{(2)}\) only provides limited insights into the physics of extended charge structures and how they interact with each other. We extend the analysis by considering the two-point pair–hole and pair–pair correlators

in which we define the pair operator \({\widehat{n}}_{{\bf{i}}}^{{\rm{p}}}={\widehat{n}}_{({i}_{x},{i}_{y})}^{{\rm{h}}}\,{\widehat{n}}_{({i}_{x},{i}_{y}+1)}^{{\rm{h}}}\). These correlators, as shown in Fig. 3a, assume the existence of nearest-neighbour pairs along y (established using \({g}_{{\rm{hh}}}^{(2)}\)) and consider the attraction or repulsion of these pairs to other dopants or pairs. They may be seen as fully connected and normalized two-point correlators of pairs or ‘partially connected’ three-point/four-point hole correlators. Note that, for simplicity, we neglect diagonal pairs (that is, \({\widehat{n}}_{({i}_{x},{i}_{y})}^{{\rm{h}}}{\widehat{n}}_{({i}_{x}\pm 1,{i}_{y}+1)}^{{\rm{h}}}\)) associated with fluctuations along x and may thus underestimate the amount of order within the system.

a, Illustration of pair–hole, pair–pair and hole–hole–hole correlator, in which a pair is defined as a nearest-neighbour pair of holes along y. b,c, Symmetrized correlation map of the pair–hole (pair–pair) correlator. We find an attraction of the pairs along y, which points towards the formation of larger-scale structures. Above the map, the average over d_y (\(\bar{\,{g}^{(2)}}\)) hints at the existence of another charge structure at d_x = 4. Error bars are estimated using bootstrapping and are smaller than the marker size if not visible. d, In the symmetrized, connected three-point hole–hole–hole correlator with d^h = (0, 1), we observe a positive signal at nearest neighbours along y, which indicates the existence of longer charge structures beyond pairs of two holes (see Supplementary Information for statistical significance). The data are evaluated over the hole doping distribution as given in the Supplementary Information.

Source data

Full size image

We present the pair–hole and pair–pair correlations for the mixD system as a function of distance along x and y in Fig. 3b,c. For improved statistics, we include in our analysis all hole doping levels (see Supplementary Information) for which the offset o_δ of equation (2) becomes negligible. In both cases, we observe positive correlations along y, which extend throughout the system, indicating that individual pairs are not repelled from other holes or each other but instead align along y and tend to form stripe-like structures. Meanwhile, there is a strong anticorrelation along x for |d_y| ≤ 1, which we attribute to the antibunching of individual holes in the same chain. We also compute the average of the correlators over d_y as \(\bar{\,{g}^{(2)}}\) (top of Fig. 3b,c). This reveals a slightly positive signal at a distance of d_x = 4, qualitatively similar to Fig. 2. Although this signal is reminiscent of a charge density wave, future studies are required to confirm this hypothesis.

For further insights into the binding of larger structures, we consider the fully connected three-point hole–hole–hole correlator

in which C_disc removes all lower-order disconnected parts of the correlator (see Supplementary Information). We show the correlator for d^h = (0, 1) in Fig. 3d and find a positive signal at the closest distance along y, whereas all other distances are negative (along x) or vanish within the error bars (see Supplementary Information). This signal directly points to extended charge structures being favoured in excess of just individual hole pairs.

To provide further evidence for extended, fluctuating charge structures, we make use of the full information in our snapshots and count ‘stripes’. To this end, we define a fully filled ‘stripe’ as a connected line of holes along y, for which the pairwise distance along x between holes in neighbouring chains is at most 1 (see Fig. 4a, inset). We designate the length ℓ of this structure by the number of chains involved. We then consider the fraction ζ(ℓ) of experimental realizations containing a ‘stripe’ of at least length ℓ. In Fig. 4a, we compare the mixD case (green) with the 2d system (brown) and randomly distributed holes (grey line; see Supplementary Information) at a doping of δ = 0.111 analysed on a subsystem of 9 × 9 sites. For the mixD case, we find an excess of events for large ℓ, consistent with the tendency to form long fluctuating structures, whereas the results obtained for the standard 2d case are consistent with randomly distributed holes. Full numerical calculations are out of reach at our system size and temperature range, but a mean-field model of stripes shows quantitative agreement in the low-doping regime (see green lines in Fig. 4a and Supplementary Information). We next analyse the difference to the random distribution δζ(ℓ) as a function of doping (Fig. 4b). For all doping levels and lengths, this signal is positive in the mixD system, indicating the inclination of the system to form extended structures. The excess probability at longer lengths grows with doping as structures of increasing lengths form.

a, Normalized histogram of ‘stripes’ (as defined in the text) of at least length ℓ in mixD (green) and 2d (brown) systems at a doping of δ = 0.111 compared with a random distribution of holes (grey line). We compare with a mean-field theory (see Supplementary Information) at k_BT/t_x = 0.355 (light-green line). The inset shows the difference δζ(ℓ) to the random distribution (see Supplementary Information). Error bars are estimated using bootstrapping and are smaller than the marker size if not visible. b, The full doping dependence of δζ(ℓ), in which the excess occurrences tend towards longer lengths with doping.

Source data

Full size image

The AFM correlations in the system and their interplay with charge delocalization are crucial for the formation of stripes and leads to characteristic signatures in the spin sector¹³. Most prominently, one expects a change in the parity of the AFM order in the presence of stripes, manifesting as incommensurate magnetism of the system and splitting of the peak at (π, π) in the spin structure factor⁵ (as also known in 1d systems⁴⁸). Although our anisotropic and strongly interacting parameter regime is not favourable to investigate structure factors (see Supplementary Information), our microscopic resolution in both spin and charge sector allows us to evaluate real-space observables inaccessible in solid-state experiments. Most useful in this context are higher-order spin-charge correlators, such as the normalized, bare three-point hole–spin–spin correlator

in which d^s is the spin bond vector, d^h the distance of the bond from the dopant and we normalize by hole density \(\langle {\widehat{n}}^{{\rm{h}}}\rangle \) and the spin standard deviation \(\sigma ({\widehat{S}}^{z})\).

Previous studies have shown that, in square lattice 2d Fermi–Hubbard systems, a single mobile dopant will be surrounded by a dressing cloud of reduced spin correlations, forming a magnetic polaron^30,31,32. In 1d systems, incommensurate magnetism leads to a change in the parity of the AFM pattern across impurities^49,50. The same feature is predicted to prevail in stripe phases, making this correlator suited to revealing this specific feature in our data. We show the bare hole–spin–spin correlator as defined in equation (5) for the mixD system in Fig. 5a, in which specific spin bonds are shown for varying distances from a hole. We focus on the diagonal and next-nearest-neighbour correlators. The most prominent feature is the strongly negative correlation across the hole along x, which is consistent with a change in the parity of the local AFM pattern across a hole. Similarly, the diagonal bonds in the direct vicinity of the hole also become negative. This is another indication of fluctuations along x within charge structures. Meanwhile, the d^s = (0, 2) correlations along y are largely unaffected by the presence of a hole and retain their positive sign. The slightly negative (positive) d^s = (2, 0) (d^s = (1, 1)) bond in the background further away from the dopant is a result of the overall doping level and vanishes in the connected correlator (see Supplementary Information).

a, Hole–spin–spin correlation map. We show the bare correlator for diagonal and selected next-nearest spin bonds as a function of distance from the hole. The strongest signal is found in the sign change of the next-nearest-neighbour bond across the hole along x pointing towards a ___domain wall in the local AFM pattern. Along y, the correlations keep their expected positive sign from the AFM pattern. b, Similarly, by considering the string spin correlator (dark green) and normal spin–spin correlator (light green) at distance d = 2 (see text), we observe a change in sign, consistent with the change in parity of the AFM pattern. Shaded regions are theory results on L_x × L_y = 8 × 3 at k_BT/t_x = 0.3. Error bars are estimated using bootstrapping and are smaller than the marker size if not visible.

Source data

Full size image

Another way to explain the change in spin order across dopants is by using a spin-string correlator^49,51. This spin–spin correlator has extra sign changes for every hole between two spins in the same chain and is defined as

in which \({\widehat{R}}_{i}={{\rm{e}}}^{{\rm{i}}{\rm{\pi }}{\widehat{n}}_{i}^{{\rm{h}}}}\). Note that, for \({\widehat{R}}_{i}={\mathbb{1}}\), the common spin–spin correlator is recovered. For systems with spin-charge separation, this correlator reveals a hidden spin structure in doped AFM systems^49,50. The changes in the phase of the AFM pattern for stripe phases act in a similar fashion and can be revealed by measuring this string correlator along the direction perpendicular to the stripes (that is, along x). We show both the common spin–spin and the string correlator at distance d = 2 in Fig. 5b as a function of doping. We observe a change to a positive sign following the use of the string correlator that only varies weakly with doping, in agreement with theory predictions. These features can be directly related to the characteristic spin ___domain parity flips present in stripe phases. Note that we observe these features even without long-range AFM correlations—which are only expected at lower temperatures—because stripe-like structures already energetically favour such a local spin arrangement.

We have realized a mixD Fermi–Hubbard model using ultracold atoms and found signatures of hole pairing and extended charge ordering in a temperature regime with short-ranged spin correlations, for which the collective behaviour of charges remains poorly understood. We detect effective hole attraction in density correlations and present further evidence for the onset of fluctuating individual stripes and their interplay with the magnetic background using real-space observables. Also, the spin environment is in qualitative agreement with the formation of an AFM ___domain wall across the dopants in both three-point and string correlators. We interpret these features as signatures for the formation of individual stripes as a precursor to the ordered stripe phase. The favourable energy scales of the mixD setting pave the way for quantum simulators to study this collective phase, including the precise periodicity, fluctuations and filling⁶, and thereby provide valuable comparisons with recent results in theoretical calculations^17,52. The direct connection between mixD and 2d systems provides a possible method to study the adiabatic preparation of stripes using mixD couplings. Through the mapping to attractive interactions⁵³, new insights into the stripe phase also directly relate to the exotic Fulde–Ferrell–Larkin–Ovchinnikov phase⁵⁴. Furthermore, direct extensions to bilayer mixD systems connect our work to recently discovered high-T_c compounds, for which the mixed dimensionality seems essential for the emergence of a superconducting phase at around 80 K in bilayer nickelates^7,8,10.

Experimental sequence

We prepare a spin-balanced sample of ultracold ⁶Li atoms in the two lowest hyperfine states |F = 1/2, m_F = ±1/2⟩ in a single layer of an optical lattice following our previous work^36,55. After magnetic evaporation, we load from a crossed dipole trap into a box potential (surrounded by a reservoir with  approximately h × 2 kHz higher chemical potential) projected with a digital mirror device (see Extended Data Fig. 1a). From this, we load into optical lattices along x and y with a_x = 1.11 μm and a_y = 1.14 μm. Their depths in the following are given in units of their respective lattice recoil E_R = h²/(8Ma²), in which M is the atomic mass.

To prepare the mixD system described in the main text without introducing large density inhomogeneities, we cannot directly load into the final lattice configuration but instead follow a procedure similar to that in ref. ³⁶ (see Extended Data Fig. 1c). We first load into decoupled 1d chains along x by exponentially ramping to V_x = 3E_R, V_y = 35E_R and a scattering length of 1,160a_B, corresponding to our final on-site interactions of U = h × 4.4(1) kHz within 200 ms. At this point, we turn on a superlattice along y (\({a}_{y}^{{\rm{SL}}}=2{a}_{y}=2.28\,\mu {\rm{m}}\)) to a depth of \({V}_{y}^{{\rm{SL}}}=2{E}_{{\rm{R}}}\) within 1 ms. By tuning the relative phase between the lattice and the superlattice to the fully staggered configuration, we ensure that the spin couplings remain the same in even and odd bonds along y. This staggering creates a potential offset of Δ = 0.65(5)U between neighbouring sites to suppress tunnelling along y. We then slowly restore coupling along y by ramping the lattices in 56 ms to their final depths of V_x = 9E_R and V_y = 7E_R. We make sure to keep the interactions constant during this second ramp by adjusting the scattering length accordingly, leading to a final scattering length of 1,293a_B. For any 2d system comparison, we perform the same ramps without turning on the y superlattice. For more details on our superlattice design, see ref. ⁵⁶. For detection, we freeze out the system by ramping to V_x/y = 43.5E_R within 1.5 ms and perform spin-resolved single-site detection as described in ref. ³⁶.

The resulting system can be accurately described by a Fermi–Hubbard-type model with parameters (t, U, Δ), which can be mapped onto the t–J model of equation (1). For all settings, we have tunnelling t_x = h × 163(10) Hz, interactions U = h × 4.4(1) kHz (thus U/t_x = 27(2)) and superexchange J_x = h × 24(4) Hz. For the 2d system, we have \({t}_{y}^{{\prime} }=h\times 253(13)\,{\rm{Hz}}\), however for Δ ≠ 0U, the effective coupling t_y is negligible. The superexchange coupling J_y is nonzero in both cases with \({J}_{y}^{{\prime} }=h\times 58(7)\,{\rm{Hz}}\) for the 2d system and J_y = h × 104(23) Hz for Δ = 0.65(5)U. Owing to the strongly anisotropic spin couplings and large U/t_x, the spin correlations are not sufficiently long-ranged to expect any signal in the spin structure factor.

We estimate the temperature of our system using the spin correlations \({C}_{{\rm{ss}}}({\bf{d}})=\frac{1}{{{\mathcal{N}}}_{{\bf{d}}}}{\sum }_{{\bf{i}}}\frac{\langle {\widehat{S}}_{{\bf{i}}}^{z}{\widehat{S}}_{{\bf{i}}+{\bf{d}}}^{z}\rangle -\langle {\widehat{S}}_{{\bf{i}}}^{z}\rangle \langle {\widehat{S}}_{{\bf{i}}+{\bf{d}}}^{z}\rangle }{\sigma ({\widehat{S}}_{{\bf{i}}}^{z})\sigma ({\widehat{S}}_{{\bf{i}}+{\bf{d}}}^{z})}\) as a function of doping and compare this to matrix product states (MPS) calculations in Extended Data Fig. 2. We fit the individual doping bins to the numerical data and extract their respective temperature. As the short y direction may be subject to finite size effects in the DMRG calculations, we determine individual temperatures along x and y. By extracting temperatures per doping level, we estimate a temperature of k_BT/t_x ≈ 0.3(1) and k_BT/t_x ≈ 0.4(1) from the correlations along x and y, respectively. Owing to the doping and high interactions, spin correlations only extend up to approximately three sites along y and are shorter along x. For that reason, fully correlated ___domain walls cannot form, making structure factors particularly challenging to investigate.

Offset phase calibration

To calibrate the detuning Δ, we first need to precisely determine the relative phase between the lattice and the superlattice. For this, we load a dilute cloud into a system of decoupled double wells along y, in which V_x = 40E_R, V_y = 8E_R and \({V}_{y}^{{\rm{SL}}}=21{E}_{{\rm{R}}}\), leading to an intrawell coupling of t_y(ϕ = 0) = h × 724(80) Hz. We vary the phase between the lattice and the superlattice and measure the normalized imbalance, that is, the difference in occupation between the different parts of the double well, normalized by their summed occupation (see Extended Data Fig. 3). When we prepare symmetric double wells, the imbalance approaches zero. However, when we tune away from this configuration, we reach an imbalance of  ±1 within less than 50 mrad. This sharp transition indicates a high degree of stability and homogeneity of the relative superlattice phase within the system (see also ref. ⁵⁶ for more details). For the measurements presented here, we then work at a phase of ϕ = π/2, at which the offset between neighbouring lattice sites is highest for a given lattice depth and the interwell and intrawell couplings are identical.

We confirm the energy scales associated with a given potential offset by comparing it with our interaction energy. We prepare the system at the lattice parameters stated in the previous section and a phase of ϕ = π/2 and vary the depth of the superlattice (that is, the potential offset) in a slightly hole doped system (see Extended Data Fig. 3c,d). For offsets smaller than the bandwidth, tunnelling between sites is not yet suppressed and we create a strong imbalance. On the other hand, for large offsets around the interaction energy, we enable resonant tunnelling between sites whenever both are occupied, thus creating both an imbalance and doublons within the system. The observed scales are consistent with band-structure calculations based on our lattice parameters. Between these two regimes, the imbalance approaches zero. The maximum in imbalance around U/2 can be explained by a second-order process in which a doublon in a lower chain breaks into two atoms in the two adjacent chains (see also ref. ⁵⁶). To avoid this effect and the associated extra holes and doublons, we perform our experiments at an offset slightly above U/2. The resulting mixD system has a small residual normalized density imbalance (as can also be seen in Extended Data Fig. 1b) of about 0.037. This does not affect the validity of the Hubbard Hamiltonian of equation (1) (as the couplings are unchanged) but only leads to a slightly worse doping resolution.

Data statistics, doping histograms

In total, we collect 11,675 experimental realizations. Of these, 1,254 were taken in a 2d system with Δ = 0U, the remaining 10,421 with Δ = 0.65(5)U. Within these measurements, we slightly vary the doping level (as well as the natural fluctuations inherent to our preparation scheme) which yields a range of 10–30% hole doping. To ensure that there is no overall magnetization \({M}^{z}={\sum }_{i}\langle {S}_{i}^{z}\rangle \) within the system, we check the distribution of magnetization normalized by the system size, which is centred around zero and shows a width below shot noise (see Extended Data Fig. 4).

Connected correlators and offsets

Full connected correlator expressions

We present a variety of correlators to characterize the spin and charge order in our system. We here distinguish between bare, ‘partially connected’ and fully connected correlators. Although the bare correlator does not subtract anything, the fully connected correlator subtracts all possible lower-order contributions between all of its constituents, for example, for a two-point correlator, it removes the product of the mean operator values, whereas for a three-point correlator, it also removes all combinations of two-point correlators. Meanwhile, partially connected correlators only subtract some specific lower-order correlators: in this case, as we consider pairs as new objects, we do not subtract any correlations arising from the individual holes in the pair.

All of these different types of correlator are then helpful to extract slightly different information about the system. Although fully connected correlators are especially useful to extract small signals in higher-order correlators dominated by lower-order contributions, the bare correlator may be more interesting when higher-order correlations are actually larger than lower-order correlators.

For this reason, as well as the partially connected pair–hole and pair–pair correlator (equation (3)) and the bare hole–spin–spin correlator (equation (5)), we used the fully connected hole–hole–hole correlator defined as

Similarly, we can define a connected hole–spin–spin correlator as

For this connected correlator, we observe the same main features also shown in Fig. 5a with a dominant negative bond across the hole (see Extended Data Fig. 5). This signal is strong enough to dominate over the AFM background, changing the correlator sign even in the bare correlator shown in Fig. 5. Meanwhile, the positive diagonal and next-nearest-neighbour bonds along y far from the hole shown in Fig. 5a now vanish, as they are not related to the presence of a hole but just stem from the AFM background. We compare this to DMRG calculations with L_y = 4, δ = 0.125 and k_BT/t_x = 0.4, which shows the same main features of strong anticorrelations across the dopant and at the diagonals in the immediate vicinity.

Offset correction

As well as the subtraction of the disconnected part, we also introduce an offset correction o_δ on the hole–hole correlator. This correction arises owing to the doping fluctuations in our finite-sized system. For each realization, we prepare a system with random but fixed total atom number and magnetization (see Extended Data Fig. 4). The calculated correlations in a finite system then obey a sum rule depending on the particle number and variance.

We start by considering N fermions on V sites with density n = N/V. The local two-point correlator \(\varGamma (i,j)=\frac{\langle {\widehat{n}}_{i}{\widehat{n}}_{j}\rangle }{{n}_{i}{n}_{j}}-1\) (with \({n}_{i}=\langle {\widehat{n}}_{i}\rangle \)) after summing over all possible pairs of sites i, j can be expressed as

in which we used \({\widehat{N}}^{2}={\sum }_{i,j}{\widehat{n}}_{i}{\widehat{n}}_{j}\) and n_i ≈ n_j ≈ n in our homogeneous system. This we can relate to the variance as

If we now separate the on-site fluctuations and use fermionic statistics in which \({\widehat{n}}^{2}=\widehat{n}\) (and thus \(\varGamma (i,i)=\frac{1}{n}-1\)), we obtain

Unless the global fluctuations of N are also fermionic fluctuations (that is, multinomial, in which \({\rm{Var}}(\widehat{N})=N(1-n)\)), the sum rule in equation (11) leads to a nonzero value of Γ(i, j) for i ≠ j even at T = ∞. Note that typically \({\rm{V}}{\rm{a}}{\rm{r}}(\hat{N})\propto N\) or less, such that equation (11) is a 1/N correction, which vanishes in the thermodynamic limit.

Identifying \(\widehat{n}\equiv {\widehat{n}}^{{\rm{h}}}\), we use this result in the calculation of the hole–hole correlations in Fig. 2 and thus define the offset o_δ through

and the corrected correlation as equation (2)

with \({{\mathcal{N}}}_{{\bf{d}}}={\sum }_{{\bf{i}},{\bf{j}}}{\delta }_{{\bf{j}},{\bf{i}}+{\bf{d}}}\) and Kronecker delta δ_i,j. Most importantly, the doping-dependent offset we apply is global on all distances. Therefore, we can understand this offset as a global −1/N correction for a fixed number of particles in the system, whereas for exceedingly large global fluctuations, positive offsets can occur.

This offset correction o_δ only plays a role when selecting specific doping levels in a finite-sized system such that the total atom number is almost fixed (\({\rm{Var}}(\widehat{N})\to 0\)) and thereby leads to strong global offsets, that we hereby compensate (see Extended Data Fig. 6a). We show in Extended Data Fig. 6b the offset as a function of doping together with the nearest-neighbour hole–hole correlator values with and without applied offset. As indicated by the dashed lines, the offset without selection on a density bin is negligible. For this reason, we do not apply any corrections in Fig. 3.

Correlator from theory

When comparing the absolute values of hole–hole correlations with simulations, care needs to be taken because of the differences in doping, fluctuations and boundary conditions. All calculations are performed with open boundary conditions along x and y. Meanwhile, the potential at the edges in the experiment has a finite width, which means that the exact position of any charge feature will be fluctuating and therefore be washed out. As a result, we detect signals in \({g}_{{\rm{hh}}}^{(2)}\) but not in the density, in contrast to theory, in which stripes appear as density features³⁶. When using connected correlators on theory data, this will lead to reduced correlations. To analyse numerical results, we hence use the slightly modified correlator \({\widetilde{g}}_{{\rm{hh}}}^{(2)}({\bf{d}})\) defined as

in which, compared with equation (2), we replace the normalization by the local densities with the global doping level n^h. This effectively assumes that the density is homogeneous throughout the system instead of bunched at the centre, allowing for easier comparison with the experiment.

Statistical significance in correlation maps

The correlation maps shown in Figs. 2 and 3 do not give any indication of which data points in the map are statistically significant or fall below the noise floor of the measurement. To address this, we show in Extended Data Fig. 6c–f the same maps as in Figs. 2 and 3 for which we now set all distances with signals compatible with zero (that is, the signal being less than 1σ away from zero) to grey. All features mentioned in the main text are still clearly visible.

Further coupling terms in the mixD Fermi–Hubbard model

In the experiment, we realize a 2d Fermi–Hubbard model with anisotropic tunnel couplings and energy offset on every second site along y. In the limit of strong interactions U ≫ t_x, t_y used here, this is commonly mapped onto the t–J model. However, this approximation neglects higher-order terms that can arise in the expansion, including a crucial second-order hopping term. Although nearest-neighbour hopping is suppressed owing to the potential offset, next-nearest-neighbour hopping remains resonant in a staggered potential. We experimentally confirmed the presence of this term and its scaling \({\widetilde{t}}_{y}={t}_{y}^{{\prime\prime} }+\frac{{t}_{y}^{2}}{\varDelta }\) with direct next-nearest-neighbour tunnelling \({t}_{y}^{{\prime\prime} }\) (which is, however, negligible for our parameters) by performing single-particle quantum walks⁵⁶. This simple expression neglects interaction effects with atoms in the intermediate lattice site. For Δ = 0.65(5)U, this means that \({\widetilde{t}}_{y}\approx \frac{1.54{t}_{y}^{2}}{U}\). This could, in principle, disfavour stripe formation, as the weak Pauli repulsion associated with \({\widetilde{t}}_{y}\) could inhibit pairs at distance 2 such that only d_y = 1 hole pairs would form. In this experiment, the contribution can mostly be neglected as the principal energy scale is given by \({J}_{y}\approx 3{\widetilde{t}}_{y}\), which dominates in our parameter regime over \({\widetilde{t}}_{y}\).

Stripe-length random data generation

To interpret the stripe-length results of Fig. 4, we compare with random hole distributions with different short-ranged correlations. We first simply randomly sample holes on 9 × 9 sites (see Extended Data Fig. 7a), in which we observe strong positive signals in the mixD case and negative signals for the 2d case. However, the strong Pauli repulsion along x might have an influence on this signal. For this reason, we randomly sampled holes for which we included, in Fig. 4, the experimentally obtained anticorrelations along x (see Fig. 2). Finally, we compare with randomly placed pairs along y within the system in Extended Data Fig. 7b, exhibiting similar features. Thus, we conclude that the observed main qualitative features are relatively insensitive to the exact details of the randomly generated data and that we see a genuine stripe signal that cannot be explained by random or short-ranged correlated holes.

Numerical simulations of the mixD t–J model

We simulate the mixD t–J model,

(see equation (1)) for J_y/t_x = 0.5 and J_x/J_y = 0.3 at finite temperature using MPS through mixed state purification schemes^57,58,59. In particular, we expand the system by introducing one auxiliary site for each physical site, which allows for showing mixed physical states as pure states on an enlarged Hilbert space. A pure state in the enlarged system at finite temperature is calculated by evolving the maximally entangled, infinite temperature state \(| \varPsi (\beta =0)\rangle \) in imaginary time under the physical Hamiltonian, \(| \varPsi (\tau )\rangle ={{\rm{e}}}^{-\tau \widehat{{\mathcal{H}}}}| \varPsi (\beta =0)\rangle \), in which τ = β/2, with β the inverse temperature. Thermal expectation values \({\langle \widehat{O}\rangle }_{{\rm{T}}}\) in the physical subset are computed by tracing out the auxiliary degrees of freedom, that is,

During the imaginary time evolution, we conserve the particle number in each row N_ℓ, ℓ = 1,…, L_y, the total particle number in the auxiliary system \({N}_{{\rm{aux}}.}^{{\rm{tot}}}\) and the total spin \({S}_{{\rm{phys}}.+{\rm{aux}}.}^{z,{\rm{tot}}}\) (the latter allowing for thermal fluctuations of the total magnetization in the physical system). This results in a total of L_y + 2 symmetries used by the DMRG implementation, leading to marked speed-ups over a single global U(1) conservation in the overall physical system¹¹.

The maximally entangled state needed as a starting point of the imaginary time evolution, \(| \varPsi (\beta =0)\rangle \), is generated using specifically tailored entangler Hamiltonians^11,60. Because these states (being projected product states) are of low bond dimension (\(\chi (\tau =0)\approx {\mathcal{O}}(100)\)), local approximations of the Hamiltonian and subsequent exponentiation will suffer from large projection errors. Hence, we start by using global methods for a single step in imaginary time, after which the entanglement in the system (and the bond dimension of the thermal MPS) has sufficiently increased to switch to local methods.

Owing to the mapping of the (enlarged) 2d system to a 1d chain, the bond dimension required for a fixed accuracy scales exponentially with linear system size in the y direction. For doping scans, we limit the system size to L_x × L_y = 8 × 3 with open boundaries and hole configurations N_ℓ = 1, 2, 3 for each ℓ = 1, 2, 3. For a single hole per chain, we simulate systems up to L_y = 4. As this mixD model suffers from the fermion sign problem, these limited system sizes are still state of the art for numerical calculations while mostly allowing general qualitative comparison with the much larger experimental system. Larger system sizes have only been achieved at zero temperature, which is numerically much easier to realize in DMRG. We furthermore checked that our temperature estimations (Extended Data Fig. 2) are not affected by the finite size effects of the DMRG calculation by comparing spin correlations for L_y = 3 and L_y = 4 at δ = 0.125 and finding very similar values.

In particular, we evolve \(| \varPsi (\beta =0)\rangle \) using global Krylov schemes by a single step t_xΔτ = 0.01. Weight cut-offs are set to 10⁻¹⁰, expanding the bond dimension to \(\chi (\tau =\varDelta \tau )\approx {\mathcal{O}}(1,000)\). From here on, we switch to the local two-site time-dependent variational principle (TDVP) method⁵⁹ with time steps of t_xΔτ = 0.03, weight and truncation cut-offs of 10⁻¹⁰ and 10⁻¹², respectively, and cutting edge maximum bond dimensions of χ_TDVP = 30,000. We evolve the system to τt_x = 2.0, corresponding to a temperature of k_BT/t_x = 0.25.

Spin–spin correlations, as well as hole distributions in each leg, are exemplarily shown in Extended Data Fig. 8a for k_BT/t_x ≈ 0.4 for a system of size L_x × L_y = 8 × 4 with periodic boundaries along the short direction and N_ℓ = 1 for all ℓ = 1,…, 4. At the centre of the chains, at which the hole density peaks, an AFM ___domain wall forms, signalling the formation of a single, fully filled stripe. For a higher doping of δ = 0.25 (L_y = 3, open boundaries), we show the hole density as well as hole–hole correlations in Extended Data Fig. 8b,c, in which the two separate stripes are visible. Results as a function of temperature are shown in Extended Data Fig. 8d for d_y = 1 and d_y = 2.

Effective descriptions of stripes in the mixD t–J model

Mean-field theory

In this section, we present a mean-field theory for the stripe phase in the mixD t–J model. We focus on describing an individual stripe in the y  direction with exactly one hole per chain, bound by the magnetically mediated confining potentials. In particular, we neglect the interaction between several stripes at positions i₁ and i₂, that is, we focus on the low-doping regime. To illustrate the concept, we first consider a mean-field description of the ground state, before generalizing to finite temperature.

For t_x ≫ J_x, J_y, quantum correlations between strongly fluctuating holes and spins in squeezed space (defined in refs. ^51,61) can be neglected^{62,63,64,65,66}. Hence, we make the ansatz

in which \({| \psi \rangle }_{{\rm{sq}}}\) is the spin state of the undoped Heisenberg model in squeezed space and \({| \psi \rangle }_{{\rm{c}}}\) is the chargon wavefunction. Our starting point for the description of the single stripe is the variational Gutzwiller wavefunction, given by

that is, we describe the charge sector by the product of identical single-leg wavefunctions \({| {\phi }^{(0)}\rangle }_{y}\) in chain y. Assuming that the stripe is centred around x = 0, we express \({| {\phi }^{(0)}\rangle }_{y}\) within the string basis,

in which Σ can be understood as the length of the string measured relative to the centre of the stripe.

Within this variational ansatz, coefficients \({\phi }_{\varSigma }^{(0)}\) can be found by minimizing the energy of the trial state, \({\langle \widehat{{\mathcal{H}}}\rangle }_{0}=\langle \psi | \widehat{{\mathcal{H}}}| \psi \rangle \,=\,({\langle \psi | }_{{\rm{sq}}}\otimes \)\({\langle \psi | }_{{\rm{c}}})\widehat{{\mathcal{H}}}({| \psi \rangle }_{{\rm{sq}}}\otimes | \psi {\rangle }_{{\rm{c}}})\),

Here V_pot(Σ) is the interchain potential defined by the potential energy of two holes in neighbouring chains separated by the string Σ,

in which \({C}_{1}^{\mu }=\langle {\psi }_{{\rm{s}}}| {\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{i}}+{{\bf{e}}}_{\mu }}| {\psi }_{{\rm{s}}}\rangle \), μ = x, y are nearest neighbours and \({C}_{2}=\langle {\psi }_{{\rm{s}}}| {\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{i}}+{{\bf{e}}}_{x}+{{\bf{e}}}_{y}}| {\psi }_{{\rm{s}}}\rangle \) are diagonal spin–spin correlations in the undoped Heisenberg model in the ground state. Note that there are also intrachain contributions, which, however, are constant and only lead to a trivial energy shift on top of the Heisenberg ground state energy E₀ (see Extended Data Fig. 9a).

By averaging over the upper and lower chains for a given leg, we can reformulate the variational problem, equation (20), as a self-consistent ground state search of the mean-field Hamiltonian per chain,

in which \({\widehat{h}}_{\varSigma }^{\dagger }| y,0\rangle =| y,\varSigma \rangle \) and

Note the factor of 2 in the potential energy, arising from energy contributions between chains y ± 1 with chain y. When considering the total energy of the variational wavefunction, equation (20), however, there is no extra factor to not double count interchain energy contributions.

In practice, we set a maximal cut-off for the string length, here chosen as \(| {\varSigma }_{\max }| \approx 15\). By exact diagonalization and self-consistently solving equation (22), the string-length distribution \(| {\phi }_{\varSigma }^{(0)}{| }^{2}\) within the mean-field picture can be calculated.

At finite temperature, we generalize the ansatz to a product of density matrices,

in which

defines the self-consistency equation through

Here \({C}_{1}^{\mu }(T)=\langle {\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{i}}+{{\bf{e}}}_{\mu }}{\rangle }_{T}\), μ = x, y and \({C}_{2}(T)=\langle {\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{i}}+{{\bf{e}}}_{x}+{{\bf{e}}}_{y}}{\rangle }_{T}\) entering V_pot in equation (21) are thermally averaged two-point correlators of the 2d Heisenberg model. Given the self-consistent solution of \({\widehat{\rho }}_{{\rm{MF}}}^{(0)}\), the mean-field string-length distribution is determined by the diagonal elements of \({\widehat{\rho }}_{{\rm{MF}}}^{(0)}\), that is, \({p}_{\varSigma }=\langle \varSigma | {\widehat{\rho }}_{{\rm{MF}}}^{(0)}| \varSigma \rangle \).

We use finite-temperature DMRG methods (see previous section) to calculate thermally averaged nearest-neighbour and diagonal correlations of the undoped Heisenberg model with J_x/J_y = 0.3 on a L_x × L_y = 12 × 4 lattice with periodic boundaries along y; see Extended Data Fig. 9b. Results for the corresponding mean-field estimates of the string-length distributions in the stripe phase are shown in Extended Data Fig. 9c for t_x/J_y = 2 and temperatures k_BT/t_x = [0.2, 0.625].

Using the mean-field theory string-length distributions, we sample snapshots and compare the resulting stripe-length distributions to the experiment (see Extended Data Fig. 9f). At the expected temperature of k_BT/t_x ≈ 0.3, the effective description matches the experiment rather well, with only a slight overestimation of the order in the mean-field description.

Müller–Hartmann–Zittartz estimate

To make further comparisons with statistical models, we reduce the mixD system to a 1d, purely classical model of fluctuating holes bound together by the effective potential V_pot (equation (21); Müller–Hartmann–Zittartz (MHZ) approach). Denoting with x_ℓ the x position of the doped hole in leg ℓ (we again consider one hole per chain, that is, a single fluctuating ___domain wall), the effective Hamiltonian (excluding quantum fluctuations from the hopping of the holes) for a system of size (L_x + 1) × (L_y + 1) reads

in which again the temperature-dependent correlators \({C}_{1}^{\mu }(T)\,=\) \(\langle {\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{i}}+{{\bf{e}}}_{\mu }}{\rangle }_{T}\), μ = x, y and \({C}_{2}(T)=\langle {\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{i}}+{{\bf{e}}}_{x}+{{\bf{e}}}_{y}}{\rangle }_{T}\) enter the effective potential V_pot(|x_ℓ − x_ℓ+1|; T) in equation (21).

The partition function, Z, decouples when being expressed solely by distances d_ℓ = x_ℓ − x_ℓ+1,

The probability of finding two adjacent holes at distance d in chains ℓ, ℓ + 1 is given by

shown for various temperatures k_BT/J_y in Extended Data Fig. 9d.

Fixing the first hole in the centre and sampling distances according to equation (29), we again generate snapshots of the hole configurations. Note that, although in the mean-field theory fluctuating stripes pinned to the centre were described, the classical formulation as given above captures stripes that are not pinned to the boundary and hence naturally form extended hole configurations (see also Fig. 4a). We compare the results with the experimental data in Extended Data Fig. 9g, in which, for k_BT/J_y = 0.8, we observe similar features to the experiment and the results from mean-field theory for k_BT/t_x = 0.36. Finally, we investigated the mean length of the excess stripes in the MHZ approach as a function of temperature (see Extended Data Fig. 9e). We observe an increase of this length below J, marking the onset of stripe-like structures.

Limitations of the effective descriptions

Both descriptions of the fluctuating stripe presented above are approximate, as they rely on assuming an effective confining potential (equation (21)) between holes in neighbouring chains. The latter description is derived within the geometric string approach, that is, assuming that the fluctuating charges merely displace spins in the background when they move around, without affecting their spin correlations. At low temperatures, this is a valid assumption, whereas at higher temperatures—when longer strings play an increasingly important role—we expect corrections to the confining potential. In particular, the question remains whether an unbinding transition of holes out of the stripe can take place. Latest for T → ∞ in the effective model this is expected to happen, for which spin correlations in the background C_d → 0 and thus V_pot(Σ) ≡ 0 in equation (21).

Another limitation of the effective models of an individual stripe is its limitation to one hole per chain. On one hand, extended interactions between neighbouring stripes at low temperatures can lead to ordering and the formation of a stripe phase with long-range charge and (incommensurate) spin order. On the other hand, at higher temperatures, the spatially extended nature of the individual stripes can cause interaction effects between them to play a role in the expected thermal unbinding transitions into a deconfined chargon gas⁶⁵.