Stripe order from the perspective of the Hubbard model

A microscopic understanding of the strongly correlated physics of the cuprates must account for the translational and rotational symmetry breaking that is present across all cuprate families, commonly in the form of stripes. Here we investigate emergence of stripes in the Hubbard model, a minimal model believed to be relevant to the cuprate superconductors, using determinant quantum Monte Carlo (DQMC) simulations at finite temperatures and density matrix renormalization group (DMRG) ground state calculations. By varying temperature, doping, and model parameters, we characterize the extent of stripes throughout the phase diagram of the Hubbard model. Our results show that including the often neglected next-nearest-neighbor hopping leads to the absence of spin incommensurability upon electron-doping and nearly half-filled stripes upon hole-doping. The similarities of these findings to experimental results on both electron and hole-doped cuprate families support a unified description across a large portion of the cuprate phase diagram. The phase diagram of the Hubbard model is studied numerically by varying parameters and suggests that spin stripe order can be observable at accessible temperatures. A team led by Thomas P. Devereaux from Stanford University and colleagues from SLAC National Accelerator Laboratory and University of North Dakota investigate emergence of spin stripe orders in the Hubbard model by tuning various parameters in the determinant quantum Monte Carlo simulations and the density matrix renormalization group calculations. They show that including the next-nearest-neighbor hopping term, which was often neglected in previous studies, in the Hubbard model leads to nearly half-filled spin stripes upon hole-doping, while no stripes upon electron-doping. The consistence of these findings with experimental results on both electron and hole-doped cuprate superconductors supports a unified description across a large portion of the cuprate phase diagram.


INTRODUCTION
The lack of an analytic solution to the Hubbard model in twodimensions has led to development of various numerical methods to study its low temperature and ground state properties. Calculations to benchmark these techniques have revealed that different candidate ground states all lie close in energy, 1,2 with small differences possibly associated with specific aspects of each method. Density matrix renormalization group, exact diagonalization/dynamical mean-field theory, constrained path auxiliary field Monte Carlo, infinite projected entagled-pair states, and density matrix embedding theory all find evidence for stripes, [1][2][3][4][5][6][7] having stronger amplitudes and longer correlation lengths than d-wave superconductivity. However, dynamical cluster approximation and cellular dynamical mean-field theory calculations have not shown evidence for stripes, instead finding a finite temperature transition into a d-wave superconductor. [8][9][10][11][12][13] These seemingly different ground states with similar energies reflect a delicate balance, sensitive to the specific nuances and biases of each approach.
Numerically discerning energy differences to ascertain lowtemperature properties requires rigorous effort to eliminate biases, and techniques may or may not reveal true ground states if the treatments are variational. On the other hand, provided that fluctuating orders are appreciable, calculations at higher temperatures provide an alternative perspective and carry the benefit that shortened correlation lengths reduce finite size effects. Here we use determinant quantum Monte Carlo (DQMC), an exact finite temperature technique, for this purpose. Although the fermion sign problem sets a lower bound on the range of temperatures amenable to simulation, we show that fluctuating stripe order is nevertheless observable at accessible temperatures.

RESULTS
We first describe the doping dependence of spin correlations for the Hubbard model with interaction strength U/t = 6 and nextnearest-neighbor hopping t′/t = −0.25, where t is the nearestneighboring hopping. Figure 1 displays the real space, equal-time spin-spin correlation functions obtained from finite temperature DQMC simulations on 16 × 4 clusters with periodic boundary conditions. At half-filling (Fig. 1a), antiferromagnetic spin correlations are evident from the checkerboard pattern, or equivalently from the uniform phase of the staggered spin-spin correlation functions (lower portion of Fig. 1a), defined with a (π, π) phase factor that flips the sign of the correlation function on every other site. Hole-doping (Fig. 1b, p = 0.042) first induces a decrease in correlation length, followed by development of antiphase domains with increasing hole concentration. The size of each domain is inversely proportional to the hole doping level; and for p ≥ 0.125, multiple sets of antiphase domain walls become visible for this cluster geometry and size. This behavior, qualitatively and quantitatively similar to previous findings for the three-band Hubbard model, 14 demonstrates that stripe behavior at finite temperatures emerges in the Hubbard model through incommensurate spin correlations. To ensure that these findings are not artifacts of the anisotropic cluster geometry, we present and discuss results for a square geometry in Fig. S1 of the Supplementary Materials.
Previous finite-temperature calculations of the Hubbard model failed to demonstrate spontaneous development of either spin or charge incommensurability, absent imposing inhomogeneity from external fields not part of the original model (e.g., see refs. 15,16 ). In light of the results presented here, a necessary ingredient appears to be large enough clusters capable of supporting multiple stripe domains. While calculations utilizing small clusters (, N = 8) have been used to demonstrate antiferromagnetism, superconductivity, and pseudogap physics, 8,13,17 their inability to host incommensurate states leaves open important questions regarding the interplay of stripes with the aforementioned orders.
In contrast to these previous finite temperature findings, zerotemperature calculations, from a variety of methods, have indicated striped ground states in the Hubbard model. A recent comparison found close agreement in the ground state energies using four different techniques, 2 providing evidence for period-8 stripes in the ground state of the 1/8-hole-doped Hubbard model with only nearest-neighbor hopping (t′ = 0) and canonical interaction strength (U/t = 8, the non-interacting bandwidth). Instead, our findings show stripes with a period~5 for 1/8 holedoping (Fig. 1b, p = 0.125), in better agreement with experimental results on cuprates. [18][19][20] To understand these differences, we study the impact of varying the model parameters, starting with the next-nearestneighbor hopping t′, which induces a particle-hole asymmetry. For t′ = −0.1 and t′ = 0, we find antiphase domain walls still present, but with an increased period of~8 for 1/8 hole-doping ( Fig. 2), similar to the previously mentioned results of ground state calculations. In contrast, the period~5 stripes from simulations using a negative value of t′ = −0.25 correspond well to neutron scattering experiments where multiple hole-doped compounds show incommensurablity corresponding to period 4-5 spin stripes at 1/8 hole-doping. 20 Also at 1/8 hole-doping, Fig. 2 shows the staggered spin-spin correlation function for t′/t = 0.1 and t′/t = 0.25, equivalent to 1/8 electron-doping for negative t′. In contrast to previous results, no antiphase domains are present and only antiferromagnetism is visible. This is additionally corroborated by our DMRG simulations in Fig. S2 of the Supplementary Materials. As neutron scattering 21,22 on electron-doped compounds similarly finds only commensurate antiferromagnetic excitations at low energy, our simulations show that a negative value of t′ that properly captures the cuprates' Fermi surface topology also correctly describes the spin behavior in both directions of doping.
Variations in the interaction strength U make little direct impact on the presence or periodicity of stripes. We first consider results for the lowest temperature accessible to simulation. For U/t = 5, at a temperature of T/t = 0.20, we again find period-5 stripes at 1/8 hole-doping (Fig. 3, top left). Increasing to U/t = 7 (Fig. 3, top right), the worsened sign problem constrains the lowest accessible temperature to T/t = 0.26. Here, the stripes instead have an increased period of~7. We attribute this to the change in temperature: for the same ratio T/J = 0.45 of temperature to exchange coupling, similar period~7 domains are present for U/t = 5 and U/t = 6 ( Fig. 3, second row). At U/t = 8, the sign problem is too severe to achieve temperatures of T/J = 0.45, but at accessible temperatures we find no indication of different behavior. The similarities between these results for the same T/J imply a marginal role of the value of U, at least in the range of explored values.
Generally, with increasing temperatures we find slight increases in stripe period and substantial reduction in correlation length (Fig. 3). This is consistent with neutron scattering data on La 1.875 Ba 0.125 CuO 4 (LBCO), where spin incommensurability (inversely proportional to the period) decreases with increasing temperature. 19 In our data, reduced correlation lengths at higher temperatures make it increasingly difficult to see π-phase shifted domains in correlation functions. At the temperature T/J = 0.75, nearly all correlations expected from the nearest π-phase shifted  We additionally compare our finite-temperature DQMC results to zero-temperature DMRG calculations 25 for the same model parameters and lattice geometry with cylindrical boundary conditions. Figure 4 displays the spin-spin correlation function and density profile from both techniques. Each shows anti-phase domains characteristic of spin stripes (Fig. 4a, b), with period-4 stripes at T = 0 and period-5 stripes at T/t = 0.22. The correlations at zero temperature (Fig. 4a) show immobile anti-phase domain walls that are pinned by the open boundaries. At finite temperatures, we observe short-ranged and mobile antiphase domain walls, as evidenced by their tendency to follow the reference spin position on the cluster, as expected from fluctuating stripes (Fig. 4b). This contrasting behavior is evident in the density profiles as well (Fig. 4c, d), where zero-temperature DMRG results reveal a static charge density wave, with peaks and troughs in the electron occupation coinciding with the antiphase domains, in agreement with the long-established picture of stripes known from the earliest mean-field studies. 26,27 The finite temperature DQMC results instead show only minor modulation of the occupation due to boundary effects, without any indication of static charge order.
Finally, to understand the dynamical properties associated with spin stripes, we use the maximum entropy method (MEM) for analytic continuation 28,29 to extract the dynamical spin structure factor S(Q, ω) from the finite temperature, unequal-time spin-spin correlation function obtained using DQMC. At zero doping (Fig.  5a), the dispersion is conical around (π, π), as expected for antiferromagnetism. Upon hole doping to p = 0.125 (Fig. 5b), the most drastic change is hardening and loss of spectral weight at (π, π), as expected when antiferromagnetism is no longer dominant. While the spin excitations at the closest wavevectors accessible to our cluster (7π/8, π) and (9π/8, π) exhibit similar behavior, the excitations at (3π/4, π) and (5π/4, π) instead soften. This change is also reflected in zero-energy spin structure factor, which shows two peaks (Fig. 5b, bottom) split from (π,π). In neutron scattering experiments on hole-doped cuprates, low-energy incommensurate peaks at these wavevectors near (π,π) are found and interpreted as evidence for stripes. [30][31][32] Although our data do not resolve sharp peaks as in neutron scattering data, which are taken at far lower temperatures, the agreement with the realspace data of Fig. 1 strongly supports that doping induces not only weakened antiferromagnetism but also incommensurate spin correlations.
As in Fig. 2, we study the effect of t′ by fixing the doping at p = 0.125 and varying t′/t from −0.25 (Fig. 5b) to 0 (Fig. 5c) and 0.25 (Fig. 5d). At t′/t = 0, similar behavior is present as in t′/t = −0.25, but with a significantly smaller separation of the two low-energy peaks. This is consistent with our real-space data in Fig. 2 that shows increased stripe period for t′/t = 0. In contrast, for positive t′/t = 0.25, low-energy spectral weight is evidently centered at (π, π), in agreement our previous real-space results in Fig. 2 that  spin correlations are predominantly antiferromagnetic for t′/t = 0.25. Subject to the limited momentum resolution of our calculation, the low-energy behavior seen here is reminiscent of the commensurate excitations at (π, π) seen in electron-doped cuprates. 21,22

DISCUSSION
We have presented the first unbiased study demonstrating fluctuating spin stripes in the Hubbard model at finite temperatures. Without applying external fields, we observe incommensurate spin correlations for a wide range of hole-doping at temperatures below roughly T/J < 0.7. We study the stripes' dependence on kinetic frustration (t′), finding that a value appropriate for cuprates (t′/t = −0.25) captures the experimentally observed doping asymmetry in the spin incommensurability. This result provides a simple resolution, without non-local Coulomb interactions, for the mismatch between experimental data on cuprates and ground-state calculations of the t′ = 0 Hubbard model.
The role of kinetic frustration for stripes may be understood from a strong-coupling perspective where magnetic moments are localized. Consider first the Hubbard model with only nearest neighbor hopping. When doped, a stripe state is more favorable than a uniformly doped antiferromagnetic state. In the former, virtual hopping of holes on antiphase domain walls preserves the local antiferromagnetism of stripe domains; in the latter, the motion of each individual hole leaves behind a trail of frustrated spins. The contribution of diagonal next-nearest neighbor hopping is the opposite: diagonal motion of individual holes in an antiferromagnetic background does not magnetically frustrate, but diagonal motion of holes in a stripe state breaks up the stripes and also frustrates the local antiferromagnetism. Thus, when diagonal motion becomes more energetically favorable (i.e., larger values of t′/t), the system tends more toward a uniform antiferromagnetic state. This is precisely what our simulations demonstrate in Fig. 2; as t′/t is adjusted from −0.25 to +0.25, the stripes grow in period until only uniform antiferromagnetism is present.
In our DMRG calculations, we find interlocked spin and charge stripes in the ground state of the simulated cluster. In contrast, for DQMC calculations using the same parameters and cluster we see instead only fluctuating and unpinned spin stripes down to the lowest accessible temperatures. Taken together, these results suggest that in the Hubbard model, static charge order sets in at temperatures beyond the scope of DQMC.
While incommensurate spin and charge density waves are present among multiple hole-doped cuprate families, interlocked stripes have been seen only in La-based compounds. In other compounds, while the charge ordering has wavevector (and hence periodicity) in the same range as that of La-based cuprates and is similarly pronounced around 1/8-doping, the doping dependence differs. 33 In La-based cuprates, both charge ordering wavevectors and spin incommensurability (periodicities) increase (decrease) with doping. For other compounds, and possibly LBCO at higher temperatures, 34 the charge ordering wavevector decreases with doping while the spin incommensurability, in observed cases, increases, 20 indicating decoupling of charge and spin ordering. Microscopic model calculations, such as the DMRG results presented here, demonstrate behavior similar to that in Labased cuprates. To date, the diverse and materials-specific behavior of charge ordering in the cuprates has not been captured in such calculations.
The growth of modern computing power combined with recent developments in computational techniques has allowed for significant progress toward understanding static properties of the Hubbard model. 2 However, to connect with experimental results, it has become necessary to determine and benchmark the dynamical properties. Here, our calculation of the dynamical spin structure factor captures both the fall-off of antiferromagnetism and the development of low-energy incommensurate peaks indicative of stripes. As our spectra's resolution is ultimately limited by cluster size and temperature, both of which are severely constrained by the fermion sign problem, development of new techniques and algorithms capable of calculating dynamics in larger-scale simulations will allow for a closer comparison against data from neutron and X-ray scattering experiments on the cuprates.

Hubbard model
The Hubbard model Hamiltonian is where c y iσ (c iσ ) creates (annihilates) an electron with spin σ at site i; n iσ ¼ c y iσ c iσ ; the hopping amplitude t ij is equal to t (the simulation energy scale) if i and j are nearest-neighbors and equal to t′ for next-nearest neighbors; U is the on-site repulsive Coulomb interaction; and the chemical potential μ controls the doping level.

Determinant quantum Monte Carlo
We perform DQMC simulations on the Hubbard model 35,36 with parameters U/t and t′/t and temperature T/t as indicated in the corresponding text and figure captions. The chemical potential is tuned to achieve the desired doping level to within an accuracy of O(10 −4 ). The imaginary time interval [0, β] is partitioned into steps of size 0.07 < Δτ < 0.12, resulting in a negligible Trotter error for our parameters.
To ensure numerical stability in computing the equal-time Green's functions, we use the prepivoting stratification algorithm as described in ref. 37 , allowing up to 8 matrix multiplications before performing a QR decomposition. The unequal-time Green's functions are constructed using the Fast Selected Inversion algorithm described in ref. 38 , with blocks corresponding to the product of matrices from eight time steps.
Equal-time measurements are performed 20 times per full space-time sweep. As computing the unequal-time Green's function is of comparable computational cost as a Monte Carlo sweep, unequal-time quantities are measured every 4th sweep. For each Markov chain, we use 50,000 sweeps for warmup and 1,000,000 sweeps for measurements. Between 2000 and 60,000 independently seeded Markov chains are used for each parameter set. This large amount of data allowed for excellent statistics despite the severe fermion sign problem: standard errors in the plotted spin correlation functions are O(10 −6 ) while the average sign dipped as low as~0.015 for p > 0.1.

Density matrix renormalization group
Using the same model parameters, we perform the standard DMRG simulations 25 with up to 26 sweeps and keep up to m = 10,000 DMRG block states with a typical truncation error of 2 × 10 −6 per step. This leads to excellent convergence for the results that we report here.

Equal-time spin-spin correlation function
The equal-time spin-spin correlation function is S z ðiÞS z ðjÞ h i ; where S z (i) = 1 2n i" Àn i# À Á is the z-component of spin on site i. To reduce statistical errors, we average correlations over pairs of sites equivalent by the translation and reflection symmetries of the cluster.

MEM analytic continuation
The imaginary time susceptibility is obtained by Fourier transforming the imaginary time spin-spin correlation function T τ S z ði; τÞS z ðj; 0Þ h i measured in DQMC. It is related to the real frequency susceptibility by χðQ; τÞ ¼ Z 1 0 dω π e Àτω þ e ÀðβÀτÞω 1 À e Àβω ImχðQ; ωÞ The dynamical spin structure factor is calculated by the fluctuation-dissipation theorem, which simplifies to SðQ; ωÞ ¼ ImχðQ; ωÞ 1 À e Àβω Since inverting Eq. (3) is numerically ill-posed, we use Maximum Entropy analytic continuation 28 to extract Imχ(Q, ω) from the DQMC data. The classic variant of MEM is used with Bryan's method for optimization. A Lorentzian centered at 2J with width 1J is chosen as a minimally informative model function with the expected high-frequency decay.

Error analysis
For our DQMC data, we use jackknife resampling to calculate standard errors. For the analytically continued spectra, we verify the repeatability and reliability of the MEM procedure by repeating the calculation on disjoint sets of data bins.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.