## Numerics converging on stripes

The Hubbard model (HM) describes the behavior of interacting particles on a lattice where the particles can hop from one lattice site to the next. Although it appears simple, solving the HM when the interactions are repulsive, the particles are fermions, and the temperature is low—all of which applies in the case of correlated electron systems—is computationally challenging. Two groups have tackled this important problem. Huang *et al.* studied a three-band version of the HM at finite temperature, whereas Zheng *et al.* used five complementary numerical methods that kept each other in check to discern the ground state of the HM. Both groups found evidence for stripes, or one-dimensional charge and/or spin density modulations.

## Abstract

Upon doping, Mott insulators often exhibit symmetry breaking where charge carriers and their spins organize into patterns known as stripes. For high–transition temperature cuprate superconductors, stripes are widely suspected to exist in a fluctuating form. We used numerically exact determinant quantum Monte Carlo calculations to demonstrate dynamical stripe correlations in the three-band Hubbard model, which represents the local electronic structure of the copper-oxygen plane. Our results, which are robust to varying parameters, cluster size, and boundary conditions, support the interpretation of experimental observations such as the hourglass magnetic dispersion and the Yamada plot of incommensurability versus doping in terms of the physics of fluctuating stripes. These findings provide a different perspective on the intertwined orders emerging from the cuprates’ normal state.

Recent experiments have established charge stripes as universal in underdoped cuprate superconductors (*1*, *2*). In contrast, no consensus exists regarding the universality of spin stripes, which are present and intimately tied to charge stripes in many doped Mott insulators (*1*, *3*–*5*) but are absent, at least in the static long-range form, in the majority of cuprates. Whether spin stripes exist in a more subtle fluctuating form in these cuprates remains an open and controversial question, of importance because of theoretical proposals suggesting a link between fluctuating stripes and the mechanism of high-*T*_{c} superconductivity (*6*–*10*). The evidence for fluctuating spin stripes in the cuprates has revolved around ubiquitous observations of an hourglass-shaped magnetic excitation spectrum (*11*, *12*). Its presence both in compounds that exhibit static stripe order (*13*) and in those that do not (*14*, *15*) finds a natural explanation in the concept of fluctuating stripes (*9*, *16*). However, alternative interpretations based on itinerant electrons exist (*17*), and conclusive experimental evidence for fluctuating stripes remains elusive. Characterizing the nature of stripes in microscopic models provides an important alternative lens for investigating the physics of stripes in the cuprates.

Early mean-field studies of the Hubbard model (*18*, *19*) have revealed some essential attributes of stripes: a propensity for doped holes to aggregate into lines of charge that correspond to antiphase boundaries of antiferromagnetic domains. Since then, more sophisticated methods also have substantiated the presence of charge and spin stripes in the ground state of the Hubbard model (*20*–*23*), including recent tensor network studies indicating that charge stripes and uniform superconducting states have very close ground-state energies (*24*). These efforts have investigated only ground-state properties; stripe phenomena in the disordered phase and the role of thermal fluctuations have received less attention. Furthermore, existing numerically exact finite-temperature calculations of the doped Hubbard model show only short-range antiferromagnetism and no sign of incommensurate spin or charge ordering (*25*). However, these studies are limited by cluster sizes smaller than the expected stripe periodicity. Here, we used numerically exact determinant quantum Monte Carlo (DQMC) simulations on substantially larger rectangular clusters than those used in previous studies. The horizontal dimensions are large enough to support multiple stripe domains, mitigating boundary effects that may frustrate striped correlations; the total system size is kept sufficiently small to be computationally tractable and to avoid an unmanageable sign problem (*26*). Comparing our results with density matrix renormalization group (DMRG) simulations on identical systems allowed us to connect zero- and finite-temperature results to fully characterize the presence of stripes.

We chose a three-band Hubbard model of a Cu-O plane, defined by the Hamiltonian(1)where the hopping integral *t _{ij}* for the bond ⟨

*ij*⟩ between orbitals

*i*and

*j*is ±

*t*for nearest-neighbor copper-oxygen bonds and ±

_{pd}*t*for next-nearest-neighbor oxygen-oxygen bonds. The chemical potential μ controls the overall doping, and the charge transfer energy Δ

_{pp}*introduces a site energy on oxygen orbitals*

_{pd}*p*. The local Coulomb interaction

*U*on orbital

_{i}*i*is

*U*for copper orbitals and

_{d}*U*for oxygen orbitals. We used the parameters

_{p}*U*= 6 eV,

_{d}*U*= 0,

_{p}*t*= 1.13 eV,

_{pd}*t*= 0.49 eV, and Δ

_{pp}*= 3 eV. The DQMC simulation temperature was set to*

_{pd}*T*= 1/12 eV ≈ 970 K; lower temperatures led to a prohibitively severe sign problem. See (

*27*) for details, including an exploration of a range of parameters consistent with those found in the literature that yield good agreement with experiments (

*28*).

Figure 1 presents the real-space, equal-time spin correlation function from our finite-temperature DQMC simulations, defined as *S*(**i**, **j**) = 〈*S _{z}*(

**i**)

*S*(

_{z}**j**)〉, where

*S*(

_{z}**i**) = (

*n*

_{i}_{,↑}–

*n*

_{i}_{,↓})/2 is the

*z*component of the spin on the copper orbital of unit cell

**i**. Correlations equivalent by the translational symmetry of the periodic cluster are averaged and plotted as

*S*(

**r**) = (1/

*N*)Σ

_{i}

_{=}

_{j}

_{+}

_{r}

*S*(

**i**,

**j**), where

*N*= 64 is the number of unit cells in the 16 × 4 cluster. In the undoped state with one hole per unit cell (Fig. 1A), as in prior studies, copper spin correlations are dominated by commensurate antiferromagnetism, evident through the checkerboard pattern in the spin correlation function or equivalently the uniform phase of the staggered spin correlation function . At a hole doping of

*p*= 0.042 (Fig. 1B), where the doped holes predominantly reside on oxygen orbitals, antiferromagnetism persists but with decreased correlation length. Further doping reveals copper spin correlations that do not simply decay but exhibit periodic phase inversions. This can be seen in the pattern of the staggered spin correlation functions of Fig. 1B, where regions of uniform signs are separated by distinct antiphase domain walls. The presence of such domain walls is a definitive signature of stripe ordering, and their periodicity of approximately (2

*p*)

^{–1}agrees with results from previous work (

*23*) and presents a direct confirmation of stripe behavior in the disordered phase. To illustrate the stripes’ fluctuating nature, we performed a direct comparison between DQMC and ground-state DMRG simulations with identical model parameters and cluster geometry.

For the comparison, we used a cluster with periodic boundaries in the four–unit cell vertical direction and open boundaries in the horizontal direction to break horizontal translational symmetry and potentially pin any stripe ordering. Figure 2A shows the staggered copper spin correlation function calculated by DMRG for a hole doping of *p* = 0.125. Here, antiphase domains with periodicity similar to that in the *p* = 0.125 panel of Fig. 1B are present. By varying the reference point of the correlation function, it is clear that the locations of the phase inversions are pinned by the open boundaries, corresponding to a picture of static stripes. We note that for some reference points, the nearest domain walls are sometimes shifted by one unit cell with respect to those seen for other reference points, owing to contributions from short-ranged antiferromagnetic correlations; however, the pinned locations of the domain walls are immediately clear by comparing with the panels for the other reference points.

This behavior stands in sharp contrast to the results from our finite-temperature DQMC simulations with the same open boundary conditions and model parameters (Fig. 2B). In every panel of Fig. 2B, the structure and periodicity of the domain walls relative to the reference point are nearly identical to what is seen in the periodic boundary result of Fig. 1B. This qualitative departure from the ground-state behavior seen in the DMRG simulations demonstrates that at sufficiently high temperature, stripes are delocalized and fluctuating rather than pinned by the open boundaries. For the temperature of the DQMC simulation, a lack of static long-range stripes (as in the DMRG results) is not surprising. Seeing vestigial signatures of the ordered phase in Figs. 1B and 2B is far less expected and provides compelling evidence for the widespread nature of fluctuating stripes in the three-band model.

The elevated temperatures where stripe correlations are seen imply surprisingly strong stripe correlations over a broad doping range. As shown in (*27*), stripe order is robust to different choices of Hubbard model parameters (figs. S1 and S2). Moreover, stripe order persists for larger rectangular clusters (figs. S3 and S4), and additional stripes begin to develop as the transverse dimension increases (8 × 8 and 10 × 10; fig. S5). This is consistent with DMRG results showing strong stripe tendencies for larger cluster sizes (*20*, *21*), indicating that our observations are not artifacts of our choice of cluster geometry. The fact that the DMRG and DQMC results are consistent with each other corroborates the usefulness of both methods and confirms the robustness of the measured stripe phenomena.

To draw a closer connection to experimental results, we calculated the dynamical spin structure factor *S*(**Q**, ω) by analytically continuing our DQMC data using the maximum entropy method, which is regarded as a standard procedure for extracting real-frequency spectra from imaginary-time data (*29*). Figure 3 displays the calculated spectra along a horizontal cut through the antiferromagnetic ordering wave vector (π, π), in units where the lattice constant *a* = 1. We first consider the undoped spectra (Fig. 3A) as a reference. Despite the broadening effects of the temperature (*T* = 1/12 eV ≈ 1/4 *J*) and finite cluster size, the structure factor exhibits a clear intensity peak and minimum in dispersion at (π, π), as expected in linear spin wave theory for antiferromagnets. Upon hole doping (Fig. 3, B to F), although the high-energy portions are unaffected, the soft excitations at and nearest to (π, π) lose spectral weight while hardening, corresponding to the increase in the singlet-triplet gap and the destruction of antiferromagnetic ordering. In the intermediate region, at wave vectors with incommensurability corresponding to the real-space periodicity in Fig. 1B, a qualitatively distinct behavior emerges. In particular, at **Q** = (3π/4, π) and **Q** = (5π/4, π) (corresponding to period-4 antiphase domain walls), systematically tracking the evolution of the structure factor with doping in Fig. 3G reveals that until roughly 0.125 hole doping, the spectral weight is maintained while the excitations soften. As the low-energy intensity peak originally at (π, π) continues to bifurcate, further doping beyond 0.125 results in hardening and loss of spectral weight at **Q** = (3π/4, π) and **Q** = (5π/4, π).

This nonmonotonic behavior motivates a comparison to the universal hourglass spectrum seen in inelastic neutron scattering. In Fig. 4A, we plot the center positions of momentum distribution curve (MDC) fits (*27*) to our calculated structure factor for 0.125 hole doping together with experimental data from three compounds, similarly derived from MDC fits but taken at lower temperatures (typically 10 K). The high-energy excitations at ω > 0.6 *J* show good agreement with the neutron scattering data. For lower energies, our MDC fits do not resolve the neck of the hourglass; this is to be expected given the high temperature and limited momentum resolution of the DQMC simulation. The energy distribution curve (EDC) fits, however, correctly resolve the collection of spectral intensity around ω = 0.5 *J*. With this in mind, the low-energy incommensurability agrees reasonably well with the experimental results. Furthermore, the doping dependence of the incommensurability ε [defined as the separation of the ω = 0 peaks in Fig. 4B from **Q** = (π, π), divided by 2π] follows a trend similar to the points of the Yamada plot (*12*, *30*) (Fig. 4C). The correspondence with well-established experimental results implies that the three-band Hubbard model may be capable of capturing the microscopic features necessary to understand essential collective properties of the cuprates.

The idea that thermal and quantum fluctuations cause static stripes to melt into a fluctuating state with dynamic correlations has often been discussed theoretically (*6*, *9*, *16*), but the experimental evidence remains sparse and seldom direct (*31*). Our state-of-the-art numerical calculations have shown that in the disordered phase, stripes maintain their characteristic antiphase behavior and periodicity in a fluctuating form, while being robust to variations in parameters, cluster size, and boundary condition. The fluctuating stripe order observed up to such high temperatures is a strong piece of corroborating evidence that these phenomena are strong enough to affect all electronic properties in the phase diagram. In particular, they may have a bearing on the controversy in previous studies (*24*) over the true ground state of microscopic models for cuprates; a benchmark of dynamical properties determined numerically is highly desired to go beyond comparisons of solely static properties (*32*).

## Supplementary Materials

www.sciencemag.org/content/358/6367/1161/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S10

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**We thank A. Kampf, S. Kivelson, W.-S. Lee, Y. Lee, D. Scalapino, R. Scalettar, J. Tranquada, and J. Zaanen for helpful discussions. Supported by the U.S. Department of Energy (DOE), Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under contract DE-AC02-76SF00515; the Alexander von Humboldt Foundation via a Feodor Lynen fellowship (C.B.M.); and the University of Tennessee’s Science Alliance Joint Directed Research and Development program, a collaboration with Oak Ridge National Laboratory (S.J.). Computational work was performed on the Sherlock cluster at Stanford University and on resources of the National Energy Research Scientific Computing Center, supported by DOE under contract DE-AC02-05CH11231. Data supporting this manuscript are stored on the Sherlock cluster at Stanford University and are available from the corresponding author upon request. Source code for the simulations can be found at https://github.com/cmendl/hubbard-dqmc. Sample input files are included for the three-band Hubbard model as studied in this work.