Experimental sequence
We prepare a spin-balanced sample of ultracold 6Li atoms in the two lowest hyperfine states |F = 1/2, mF = ±1/2⟩ in a single layer of an optical lattice following our previous work36,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 ax = 1.11 μm and ay = 1.14 μm. Their depths in the following are given in units of their respective lattice recoil ER = h2/(8Ma2), 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. 36 (see Extended Data Fig. 1c). We first load into decoupled 1d chains along x by exponentially ramping to Vx = 3ER, Vy = 35ER and a scattering length of 1,160aB, 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 Vx = 9ER and Vy = 7ER. 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,293aB. 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. 56. For detection, we freeze out the system by ramping to Vx/y = 43.5ER within 1.5 ms and perform spin-resolved single-site detection as described in ref. 36.
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 tx = h × 163(10) Hz, interactions U = h × 4.4(1) kHz (thus U/tx = 27(2)) and superexchange Jx = 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 ty is negligible. The superexchange coupling Jy is nonzero in both cases with \({J}_{y}^{{\prime} }=h\times 58(7)\,{\rm{Hz}}\) for the 2d system and Jy = h × 104(23) Hz for Δ = 0.65(5)U. Owing to the strongly anisotropic spin couplings and large U/tx, 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 kBT/tx ≈ 0.3(1) and kBT/tx ≈ 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 Vx = 40ER, Vy = 8ER and \({V}_{y}^{{\rm{SL}}}=21{E}_{{\rm{R}}}\), leading to an intrawell coupling of ty(ϕ = 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. 56 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. 56). 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
$$\begin{array}{l}{C}_{{\rm{hhh}}}^{{\rm{c}}}({{\bf{d}}}^{{\rm{h}}},{\bf{d}})=\frac{1}{{{\mathcal{N}}}_{{{\bf{d}}}^{{\rm{h}}},{\bf{d}}}}\sum _{\begin{array}{c}{\bf{i}}\\ {\bf{j}}={\bf{i}}+{{\bf{d}}}^{{\rm{h}}}/2+{\bf{d}}\end{array}}\frac{1}{\langle {\widehat{n}}_{{\bf{i}}}^{h}\rangle \langle {\widehat{n}}_{{\bf{i}}+{{\bf{d}}}^{{\rm{h}}}}^{{\rm{h}}}\rangle \langle {\widehat{n}}_{{\bf{j}}}^{{\rm{h}}}\rangle }\times (\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{n}}_{{\bf{i}}+{{\bf{d}}}^{{\rm{h}}}}^{{\rm{h}}}{\widehat{n}}_{{\bf{j}}}^{{\rm{h}}}\rangle \\ \,\,\,\,\,\,-\,\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{n}}_{{\bf{i}}+{{\bf{d}}}^{{\rm{h}}}}^{{\rm{h}}}\rangle \langle {\widehat{n}}_{{\bf{j}}}^{{\rm{h}}}\rangle -\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}\rangle \langle {\widehat{n}}_{{\bf{i}}+{{\bf{d}}}^{{\rm{h}}}}^{{\rm{h}}}{\widehat{n}}_{{\bf{j}}}^{{\rm{h}}}\rangle \\ \,\,\,\,\,\,-\,\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{n}}_{{\bf{j}}}^{{\rm{h}}}\rangle \langle {\widehat{n}}_{{\bf{i}}+{{\bf{d}}}^{{\rm{h}}}}^{{\rm{h}}}\rangle +2\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}\rangle \langle {\widehat{n}}_{{\bf{i}}+{{\bf{d}}}^{{\rm{h}}}}^{{\rm{h}}}\rangle \langle {\widehat{n}}_{{\bf{j}}}^{{\rm{h}}}\rangle ).\end{array}$$
(7)
Similarly, we can define a connected hole–spin–spin correlator as
$$\begin{array}{l}{C}_{{\rm{hss}}}^{{\rm{c}}}({{\bf{d}}}^{{\rm{s}}},{{\bf{d}}}^{{\rm{h}}})=\frac{1}{{{\mathcal{N}}}_{{{\bf{d}}}^{{\rm{s}}},{{\bf{d}}}^{{\rm{h}}}}}\sum _{\begin{array}{c}{\bf{i}}\\ {\bf{k}}-{\bf{j}}={{\bf{d}}}^{{\rm{s}}},\\ ({\bf{k}}+{\bf{j}})/2-{\bf{i}}={{\bf{d}}}^{{\rm{h}}}\end{array}}\frac{1}{\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}\rangle \sigma ({\widehat{S}}_{{\bf{j}}}^{z})\sigma ({\widehat{S}}_{{\bf{k}}}^{z})}\times (\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{S}}_{{\bf{j}}}^{z}{\widehat{S}}_{{\bf{k}}}^{z}\rangle \\ \,\,\,\,\,\,-\,\langle {\widehat{n}}_{{\bf{i}}}^{h}\rangle \langle {\widehat{S}}_{{\bf{j}}}^{z}{\widehat{S}}_{{\bf{k}}}^{z}\rangle -\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{S}}_{{\bf{j}}}^{z}\rangle \langle {\widehat{S}}_{{\bf{k}}}^{z}\rangle \\ \,\,\,\,\,\,-\,\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{S}}_{{\bf{k}}}^{z}\rangle \langle {\widehat{S}}_{{\bf{j}}}^{z}\rangle +2\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}\rangle \langle {\widehat{S}}_{{\bf{j}}}^{z}\rangle \langle {\widehat{S}}_{{\bf{k}}}^{z}\rangle ).\end{array}$$
(8)
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 Ly = 4, δ = 0.125 and kBT/tx = 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
$$\sum _{i,j}\varGamma (i,j)=\sum _{i,j}\left(\frac{\langle {\widehat{n}}_{i}{\widehat{n}}_{j}\rangle }{{n}_{i}{n}_{j}}-1\right)\approx \left(\frac{\langle {\widehat{N}}^{2}\rangle }{{n}^{2}}-{V}^{2}\right),$$
(9)
in which we used \({\widehat{N}}^{2}={\sum }_{i,j}{\widehat{n}}_{i}{\widehat{n}}_{j}\) and ni ≈ nj ≈ n in our homogeneous system. This we can relate to the variance as
$$\frac{1}{{V}^{2}}\sum _{i,j}\varGamma (i,j)=\frac{{\rm{Var}}(\widehat{N})}{{N}^{2}}.$$
(10)
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
$$\frac{1}{{V}^{2}}\sum _{i\ne j}\varGamma (i,j)=\frac{{\rm{Var}}(\widehat{N})-N(1-n)}{{N}^{2}}.$$
(11)
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
$${o}_{\delta }=\frac{{\rm{Var}}({\widehat{N}}^{{\rm{h}}})-{N}^{{\rm{h}}}(1-{n}^{{\rm{h}}})}{{({N}^{{\rm{h}}})}^{2}},$$
(12)
and the corrected correlation as equation (2)
$${g}_{{\rm{hh}}}^{(2)}({\bf{d}})-1=\frac{1}{{{\mathcal{N}}}_{{\bf{d}}}}\sum _{{\bf{i}}}\left(\frac{\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{n}}_{{\bf{i}}+{\bf{d}}}^{{\rm{h}}}\rangle }{\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}\rangle \langle {\widehat{n}}_{{\bf{i}}+{\bf{d}}}^{{\rm{h}}}\rangle }-1-{o}_{\delta }\right),$$
(13)
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 features36. 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
$${\widetilde{g}}_{{\rm{hh}}}^{(2)}({\bf{d}})-1=\frac{1}{{{\mathcal{N}}}_{{\bf{d}}}}\sum _{{\bf{i}}}\left(\frac{\langle {\widehat{n}}_{{\bf{i}}}^{{\rm{h}}}{\widehat{n}}_{{\bf{i}}+{\bf{d}}}^{{\rm{h}}}\rangle }{{n}^{{\rm{h}}}{n}^{{\rm{h}}}}-1\right)$$
(14)
in which, compared with equation (2), we replace the normalization by the local densities with the global doping level nh. 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 ≫ tx, ty 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 walks56. 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 dy = 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,
$$\widehat{{\mathcal{H}}}=\widehat{{\mathcal{P}}}\left(-{t}_{x}\sum _{{\langle {\bf{i}},{\bf{j}}\rangle }_{x}}\sum _{\sigma =\uparrow ,\downarrow }{\widehat{c}}_{{\bf{i}},\sigma }^{\dagger }{\widehat{c}}_{{\bf{j}},\sigma }+\text{h.c.}\right)\widehat{{\mathcal{P}}}+{J}_{x}\sum _{{\langle {\bf{i}},{\bf{j}}\rangle }_{x}}\left({\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{j}}}-\frac{{\widehat{n}}_{{\bf{i}}}{\widehat{n}}_{{\bf{j}}}}{4}\right)+{J}_{y}\sum _{{\langle {\bf{i}},{\bf{j}}\rangle }_{y}}\left({\widehat{{\bf{S}}}}_{{\bf{i}}}\cdot {\widehat{{\bf{S}}}}_{{\bf{j}}}-\frac{{\widehat{n}}_{{\bf{i}}}{\widehat{n}}_{{\bf{j}}}}{4}\right),$$
(15)
(see equation (1)) for Jy/tx = 0.5 and Jx/Jy = 0.3 at finite temperature using MPS through mixed state purification schemes57,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,
$${\langle \widehat{O}\rangle }_{{\rm{T}}}=\frac{\langle \varPsi (\beta )| \widehat{O}| \varPsi (\beta )\rangle }{\langle \varPsi (\beta )| \varPsi (\beta )\rangle }.$$
(16)
During the imaginary time evolution, we conserve the particle number in each row Nℓ, ℓ = 1,…, Ly, 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 Ly + 2 symmetries used by the DMRG implementation, leading to marked speed-ups over a single global U(1) conservation in the overall physical system11.
The maximally entangled state needed as a starting point of the imaginary time evolution, \(| \varPsi (\beta =0)\rangle \), is generated using specifically tailored entangler Hamiltonians11,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 Lx × Ly = 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 Ly = 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 Ly = 3 and Ly = 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 txΔτ = 0.01. Weight cut-offs are set to 10−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) method59 with time steps of txΔτ = 0.03, weight and truncation cut-offs of 10−10 and 10−12, respectively, and cutting edge maximum bond dimensions of χTDVP = 30,000. We evolve the system to τtx = 2.0, corresponding to a temperature of kBT/tx = 0.25.
Spin–spin correlations, as well as hole distributions in each leg, are exemplarily shown in Extended Data Fig. 8a for kBT/tx ≈ 0.4 for a system of size Lx × Ly = 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 (Ly = 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 dy = 1 and dy = 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 i1 and i2, 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 tx ≫ Jx, Jy, quantum correlations between strongly fluctuating holes and spins in squeezed space (defined in refs. 51,61) can be neglected62,63,64,65,66. Hence, we make the ansatz
$$| \psi \rangle ={| \psi \rangle }_{{\rm{sq}}}\otimes {| \psi \rangle }_{{\rm{c}}},$$
(17)
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
$${| \psi \rangle }_{{\rm{c}}}=\underset{y=-\infty }{\overset{\infty }{\bigotimes }}{| {\phi }^{(0)}\rangle }_{y},$$
(18)
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,
$${| {\phi }^{(0)}\rangle }_{y}=\mathop{\sum }\limits_{\varSigma =-\infty }^{\infty }{\phi }_{\varSigma }^{(0)}| y,\varSigma \rangle ,$$
(19)
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}}})\),
$$\frac{{\langle \widehat{{\mathcal{H}}}\rangle }_{0}}{{L}_{y}}=\frac{{E}_{0}}{{L}_{y}}-{t}_{x}\sum _{\varSigma }({\phi }_{\varSigma +1}^{(0)* }{\phi }_{\varSigma }^{(0)}\,+\,\text{c.c.})+\sum _{\varSigma ,{\varSigma }^{{\prime} }}| {\phi }_{\varSigma }^{(0)}{| }^{2}| {\phi }_{{\varSigma }^{{\prime} }}^{(0)}{| }^{2}{V}_{{\rm{pot}}}(\varSigma -{\varSigma }^{{\prime} }).$$
(20)
Here Vpot(Σ) is the interchain potential defined by the potential energy of two holes in neighbouring chains separated by the string Σ,
$${V}_{{\rm{pot}}}(\varSigma )={J}_{y}[(| \varSigma | -1+{\delta }_{\varSigma ,0}){C}_{2}-(| \varSigma | +1){C}_{1}^{y}],$$
(21)
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 E0 (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,
$${\widehat{{\mathcal{H}}}}_{{\rm{MF}}}=\frac{{E}_{0}}{{L}_{y}}-{t}_{x}\sum _{\varSigma ,{\varSigma }^{{\prime} }}[{\widehat{h}}_{{\varSigma }^{{\prime} }}^{\dagger }{\widehat{h}}_{\varSigma }\,+\,\text{h.c.}]+\sum _{\varSigma }{\widehat{h}}_{\varSigma }^{\dagger }{\widehat{h}}_{\varSigma }{V}_{{\rm{eff}}}(\varSigma ),$$
(22)
in which \({\widehat{h}}_{\varSigma }^{\dagger }| y,0\rangle =| y,\varSigma \rangle \) and
$${V}_{{\rm{eff}}}(\varSigma )=2\sum _{{\varSigma }^{{\prime} }}| {\phi }_{{\varSigma }^{{\prime} }}^{(0)}{| }^{2}{V}_{{\rm{pot}}}({\varSigma }^{{\prime} }-\varSigma ).$$
(23)
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,
$$\widehat{\rho }={\widehat{\rho }}_{{\rm{sq}}}\otimes \left(\underset{y=-\infty }{\overset{\infty }{\bigotimes }}{\widehat{\rho }}_{{\rm{MF}}}^{(0)}\right),$$
(24)
in which
$${\widehat{\rho }}_{{\rm{MF}}}^{(0)}=\frac{1}{Z}{{\rm{e}}}^{-\beta {\widehat{{\mathcal{H}}}}_{{\rm{MF}}}({\widehat{\rho }}_{{\rm{MF}}}^{(0)},T)}$$
(25)
defines the self-consistency equation through
$$\begin{array}{c}\frac{{\widehat{{\mathcal{H}}}}_{{\rm{MF}}}({\widehat{\rho }}_{{\rm{MF}}}^{(0)})}{{L}_{y}}=\frac{{E}_{0}}{{L}_{y}}-{t}_{x}\sum _{\varSigma ,{\varSigma }^{{\prime} }}[{\widehat{h}}_{{\varSigma }^{{\prime} }}^{\dagger }{\widehat{h}}_{\varSigma }\,+\,\text{h.c.}]+\sum _{\varSigma }{\widehat{h}}_{\varSigma }^{\dagger }{\widehat{h}}_{\varSigma }{V}_{{\rm{eff}}}(\varSigma \,;{\widehat{\rho }}_{{\rm{MF}}}^{(0)},T),\\ {V}_{{\rm{eff}}}(\varSigma \,;{\widehat{\rho }}_{{\rm{MF}}}^{(0)},T)=2\sum _{{\varSigma }^{{\prime} }}\langle \varSigma | {\widehat{\rho }}_{{\rm{MF}}}^{(0)}| \varSigma \rangle {V}_{{\rm{pot}}}({\varSigma }^{{\prime} }-\varSigma \,;T).\end{array}$$
(26)
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 Vpot 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 Jx/Jy = 0.3 on a Lx × Ly = 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 tx/Jy = 2 and temperatures kBT/tx = [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 kBT/tx ≈ 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 Vpot (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 (Lx + 1) × (Ly + 1) reads
$${\widehat{{\mathcal{H}}}}_{{\rm{MHZ}}}=\mathop{\sum }\limits_{{\ell }=1}^{{L}_{y}}{V}_{{\rm{pot}}}(| {x}_{{\ell }}-{x}_{{\ell }+1}| \,;T),$$
(27)
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 Vpot(|xℓ − xℓ+1|; T) in equation (21).
The partition function, Z, decouples when being expressed solely by distances dℓ = xℓ − xℓ+1,
$$\begin{array}{l}Z\,=\,\sum _{\{{x}_{{\ell }}\}}\mathop{\prod }\limits_{{\ell }=1}^{{L}_{y}}\exp [-\beta {V}_{{\rm{pot}}}(| {x}_{{\ell }}-{x}_{{\ell }+1}| \,;T)]\\ \,\,=\,\mathop{\sum }\limits_{{d}_{1}=-{L}_{x}}^{{L}_{x}}\ldots \mathop{\sum }\limits_{{d}_{{L}_{y}}=-{L}_{x}}^{{L}_{x}}\mathop{\prod }\limits_{{\ell }=1}^{{L}_{y}}\exp [-\beta {V}_{{\rm{pot}}}(| {d}_{{\ell }}| \,;T)]\\ \,\,=\,{\left[\mathop{\sum }\limits_{d=-{L}_{x}}^{{L}_{x}}\exp [-\beta {V}_{{\rm{pot}}}(| d| ;T)]\right]}^{{L}_{y}}\,=\,{[{Z}_{1}]}^{{L}_{y}}.\end{array}$$
(28)
The probability of finding two adjacent holes at distance d in chains ℓ, ℓ + 1 is given by
$$p(d)=\exp [-\beta {V}_{{\rm{pot}}}(| d| \,;T)]/{Z}_{1},$$
(29)
shown for various temperatures kBT/Jy 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 kBT/Jy = 0.8, we observe similar features to the experiment and the results from mean-field theory for kBT/tx = 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 Cd → 0 and thus Vpot(Σ) ≡ 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 gas65.
