Temperature Dependence of Spin and Charge Orders in the Doped Two-Dimensional Hubbard Model

Competing and intertwined orders including inhomogeneous patterns of spin and charge are observed in many correlated electron materials, such as high-temperature superconductors. Introducing a new development of the constrained-path auxiliary-field quantum Monte Carlo (AFQMC) method, we study the interplay between thermal and quantum fluctuations in the two-dimensional Hubbard model. We obtain an accurate and systematic characterization of the evolution of the spin and charge correlations as a function of temperature $T$ and how it connects to the ground state, at three representative doping levels $\delta = 1/5$, $1/8$, and $1/10$. We find increasing short-range commensurate antiferromagnetic correlations as $T$ is lowered. As the correlation length grows sufficiently large, a modulated spin-density-wave (SDW) appears. At $\delta = 1/5$, the SDW saturates and remains short-ranged as $T \rightarrow 0$. In contrast, at $\delta = 1/8$ and $1/10$ this evolves into a ground-state stripe phase. We study the relation between spin and charge orders and find that formation of charge order appears to be driven by that of the spin order. We identify a finite-temperature phase transition below which charge ordering sets in and discuss the implications of our results for the nature of this transition.


INTRODUCTION
The complex interplay of spin, charge and pairing correlations is a common feature of many strongly correlated materials, from manganites to cuprates [1][2][3][4][5].In the latter case, inhomogeneous orders such as stripes are a central concept [6].Spin and charge, rather than being uniformly distributed, form intricate collective ordering patterns which often have characteristic wavelengths spanning multiple lattice spacings.For example, stripe orders exhibit antiferromagnetic (AFM) regions separated by πphase shifts that reverse the sublattice occupations across lines where the holes accumulate.Related spin-density wave (SDW) states also show modulated AFM correlation, but with more uniform charge distribution in which holes spread into regions away from the spin phase shift lines.These states tend to be collinear in the spin order, and the wavelengths of the spin and charge modulation are related: λ spin = 2λ charge [7,8].Such states are the outcome of a compromise and balance among the AFM interactions, Coulomb interactions and kinetic energies [9].Some families of cuprates such as Nd-or Ba-substituted La 2−x Sr x CuO 4 (LSCO) indeed display long-range stripe ordering [6,[10][11][12][13].Taking a broader perspec- * bxiao@flatironinstitute.org † heyuanyao@nwu.edu.cn‡ ageorges@flatironinstitute.org§ szhang@flatironinstitute.org tive, charge ordering is a commonly observed feature of the low-doping ground-state of cuprate superconductors when a magnetic field is applied to suppress superconductivity [8,14,15].
The one-band Hubbard model and related models have played a key role in the study of stripe and SDW orders [16][17][18].Indications of these orders were first seen in Hartree-Fock calculations in the two-dimensional (2D) Hubbard model [19][20][21][22][23]. Calculations in the related t-J model with a hole doping δ = 1/8 using the density matrix renormalization group (DMRG) showed evidence for π phase-shifted AFM regions separated by domain walls [24].Further DMRG studies in the doped Hubbard model on cylinders found the existence of a stripe state [25].Ground-state constrained-path (CP) auxiliary-field quantum Monte Carlo (AFQMC) calculations in the Hubbard model found inhomogeneous spin and charge orders [26], which were shown to be SDW states at U/t = 4 and stripe states at larger U [27].At the thermodynamic limit, these states were shown [27] to have modulation along the x-or y-direction, with wavelength of 2/δ for the spin and 1/δ for the charge.A recent multi-method study [28] using four different computational methods confirmed that the ground state of the 2D Hubbard model at δ = 1/8 doping and strong U is indeed a vertical (x or y) stripe state with spin modulation wavelength of 16 and charge modulation wavelength of 8, i.e., filled stripes.Note that this stripe pattern is somewhat different from the one observed in real materials [6,10,11,29,30].
Determining and understanding the full phase diagram arXiv:2202.11741v3[cond-mat.str-el]22 Feb 2023 of the Hubbard model has remained an outstanding challenge.Computational methods need to reach high accuracy and sufficiently large system sizes to accommodate collective states whose modulating wavelengths tend to be large.Ground-state properties sensitively depend on the delicate balance of small energy differences.Often the signal is comparable to or smaller than the uncertainties in the numerical methods (due, e.g., to finite supercell size effects, finite bond dimensions, limitation on the temperatures and statistical error bars, approximations in the computational scheme, etc.).Recent years have seen major advances in computational methods which have led to considerable progress in the determination of the ground-state of the Hubbard model.For example, studies by DMRG [30,31], infinite project-entangled pair states (iPEPS) [32], density matrix embedding theory (DMET) [33], inhomogeneous (unrestricted) dynamical mean field theory (DMFT) [34] and its cluster extensions [35], variational Monte Carlo [36][37][38], AFQMC [39] etc. have all indicated that stripe or modulated spin orders feature prominently throughout a large portion of the ground state phase diagram in the 2D Hubbard model, either as the true ground state or energetically very close by.More subtle and controversial has been the existence of superconducting order in the model [38,[40][41][42][43][44].It is clear that a strong interplay exists between superconductivity and the stripe and SDW orders, and recent results indicate that there is a competition which appears to favor the latter in the ground state of the pure 2D Hubbard model -namely in its original form without hopping beyond nearest neighbors -at intermediate to strong coupling and near optimal doping [45].
Understanding the temperature dependence of the magnetic and charge orders is of fundamental importance.Information about the temperature evolution of stripe and SDW correlations is crucial for making connections with experimental observations.At finite temperatures, new fascinating phenomena arise such as the pseudogap regime whose understanding requires reliable characterization of spin and charge properties.However, the computational challenges are outstanding.In addition to all the requirements for reliable and predictive calculations at T = 0 K, we must now not only retain thermal fluctuations but reach sufficiently low temperatures.As a result, much less is known about the properties of the magnetic and, especially charge correlations in the doped Hubbard model at finite temperatures.
Recently there has been a flurry of activities to study the finite-temperature properties of stripes as well as spin and charge correlations in the doped Hubbard model, including calculations using determinant Quantum Monte Carlo (DQMC) [46,47], minimally entangled typical thermal state (METTS) [48], the dynamical cluster approximation (DCA) [49], and connected determinant diagrammatic Monte Carlo (CDet) [50].This has created an exciting synergy, from which a consistent picture is emerging on the temperature evolution of the antiferromagnetic spin correlations and charge inhomogeneities in the doped Hubbard model.The picture is far from complete, however.Despite impressive advances in the computational methods, technical hurdles in each study limited the scope of questions that could be addressed.These hurdles included the sizes of the simulation cells or clusters that could be accessed (which limit the detection of long-range correlations), the use of relatively narrow (width-4) cylindrical cells in some studies (whose quasione-dimensional nature may not represent the physical behavior in two-dimensions), restriction to relatively high temperatures, etc. Important questions remain on how the spin and charge correlations connect to the ground state, whether and how the system develops into longrange-ordered state, and the interplay and connection between spin and charge ordering as a function of temperature.
In this paper, we address these questions by employing finite-temperature AFQMC [51,52] with a self-consistent constraint [53] and by formulating a more flexible and powerful self-consistency scheme at finite temperatures which introduces an effective temperature in the constraint.Our new approach allows for accurate computations with much larger simulation cells at much lower temperatures than previously possible.We employ two different ways to probe the spin and charge correlations, and combine them with a careful investigation of the size dependence to extract the properties in the twodimensional thermodynamic limit.We focus on the pure Hubbard model (t = 0) and investigate three representative doping levels in the parameter regime most relevant to cuprates: δ = 1/5 (U/t = 6), δ = 1/8 and 1/10 (U/t = 8).In each case we systematically determine the temperature dependence and how the system evolves into its corresponding ground state.We find that the δ = 1/5 case shows qualitatively and fundamentally different behavior from the other two, which develop longrange ground-state order.Notably, we answer a longstanding open question in the field by identifying the onset temperature at which long-range charge order forms.We discuss possible scenarios regarding the interplay and nature of the charge-and spin-ordering transitions.
The paper is organized as follow.In Sec.II, we introduce the Hubbard Hamiltonian and give a brief overview of our computational approach, followed by a discussion and illustration of the new finite-temperature approach that optimizes the trial density matrix self-consistently using an effective temperature.Our results are presented and discussed in Sec.III.In Sec.III A, we examine the temperature dependence of the spin order in real space.This is complemented by momentum space information in Sec.III B. In Sec.III C, we focus on the overdoped regime at δ = 1/5 and U = 6, showing the temperature evolution of a system which has only short-range modulated spin correlations in the ground state.Our findings of the interplay between onsets of spin and charge ordering and the critical T c for charge long-range order at δ = 1/8 and U = 8 are then presented in Sec.III D. We conclude in Sec.IV.

II. MODEL AND METHOD
The Hubbard Hamiltonian [16,17,54,55], which was originally introduced to explain the itinerant ferromagnetism of transition-metal monoxides, serves as a minimal microscopic model to study several key features in strongly interacting quantum systems.It takes the form We consider a 2D lattice with t denoting the nearestneighbor hopping amplitude and i, j denoting a pair of nearest-neighbor lattice sites i and j.We choose t = 1 to set the unit of energy.c † i,σ (c i,σ ) is the creation (annihilation) operator of electrons with spin σ at site i = (i x , i y ), and n i,σ = c † i,σ c i,σ denotes the number operator.U is an on-site interaction between spin-up and spin-down electrons.The chemical potential µ tunes the electron density and µ = 0 corresponds to half-filling due to the choice of the particle-hole symmetric form for the interaction.The hole doping is defined as δ = 1 − ρ where ρ = N e /N denotes the electron density, N e is the particle number, and N = L x × L y is the number of sites in the simulation cell.
We solve this model by a state-of-the-art AFQMC method which controls the fermion sign problem with a constraint on the paths in auxiliary-field space [51,52].In the standard DQMC algorithm [56], the interaction part of the Hamiltonian is formulated in terms of single-particle operators, after applying the Hubbard-Stratonovich (HS) transformation [57,58] to the twobody terms and replacing them by one-body interactions with a set of auxiliary fields.The partition function, e −βH , where β = 1/T is inverse temperature, can be treated as a many-dimensional path integral over random auxiliary fields, which is computed using Monte Carlo (MC) techniques.The sign problem occurs because different configurations of the auxiliary-fields, i.e., different paths, lead to contributions to the partition function with different signs.As the temperature is decreased (i.e., the length of the path in imaginary-time increased), the average sign of the contributions approaches zero exponentially [59][60][61][62][63].The sign problem fundamentally limits the accessible temperatures and lattice sizes in simulations of strongly correlated systems by standard DQMC algorithms.
The AFQMC approach we employ shares with DQMC the decoupling of the two-body interactions by HS transformation and replacing them with one-body interactions in auxiliary-fields.However, it reformulates the partition function as a path integral over a constrained portion of the paths in auxiliary-field space.The full pathintegral gives the many-body partition function via an over-complete space.There exists an exact sign or gauge condition for a reduced subset of paths which preserve the exactness of the partition function [51,52].In practice this condition is realized approximately by a trial density matrix.In the limit T = 0 K, this approach reduces to the ground-state CP method (or more generally, the phaseless AFQMC [52,64] for arbitrary two-body Hamiltonians such as those in electronic structure calculations).There exists a large body of work which shows that this approach leads to highly accurate results in the ground state [26,51,[65][66][67].Recently a self-consistent constraint has been developed for ground-state [68] and finite-T [53] calculations.By coupling the AFQMC calculation to an independent particle (IP) Hamiltonian, an optimized constraint from single-particle wave function or density matrix can be achieved, which further reduces the systematic error.In this paper we propose an improved self-consistent scheme to further enhance the accuracy of the method in the most challenging cases at finite temperature, which we discuss further in Sec.II B below.

A. Two different approaches to characterize spin and charge orders
Orders in the Hubbard model are in general quite delicate, because of competitions and the small energy scales which separate them.At finite-temperatures, the spin and charge orders have collective modes which all decay with distance and are sensitive to finite sizes and boundary conditions and which thus require exquisite accuracy to determine.In this work, we use two different approaches to characterize these orders, and compute several quantities to probe them under each approach.
We perform AFQMC calculations in 2D using either supercells under periodic boundary condition (PBC) in both x and y directions or cylindrical cells which are periodic along y but has open boundary condition (OBC) along the x-direction.In the first case, we measure various correlation functions, for instance the equal-time spin-spin correlation in real space, where Ŝz (i) = 1 2 (n i↑ − n i↓ ) is the z-component of the spin operator at site i.This type of measurements are employed in Sec.III A. Imaginary-time-dependent correlation functions can also be measured.For example the static spin susceptibility discussed in Sec.III B, χ s (q), is which is the Fourier transform of the unequal-time spin correlation function integrated over imaginary time.The charge correlation functions and susceptibilities can be calculated in a similar way.
When the order to be detected is small, it is especially challenging to extract signals from correlation functions (and quantities derived from them such as structure factors and susceptibilities), because one is essentially measuring the square of the (small) order parameter.To ameliorate this problem, we take a complementary approach to the fully PBC calculations discussed above.We introduce "pinning fields" to break translational symmetry in one direction (x) (or in both x and y directions in certain cases as needed), so that we can measure a local "order parameter" in these systems as a proxy to the two-body correlation functions in translationally-invariant systems in the large size limit.For example, we can apply AFM pinning fields to the left edge of the simulation cell by adding the following term to the original Hamiltonian in Eq. (1): where v ix↑ = −v ix↓ = (−1) iy v for all sites with i x = 0.These spin pinning fields break the translational and SU(2) symmetries, allowing the direct detection of the order parameter [69,70]; however, care must be taken to remove any bias they introduce, as we further discuss below.
The approach of pinning fields is used to obtain the results in Sec.III C and some of the results in Sec.III D. We can similarly add charge "pinning fields", in the form H c = i P i n i .For example, we can apply a fully periodic potential to the whole simulation cell: where κ = 2π/λ is the wave number corresponding to the expected wavelength of the charge order, in order to probe the response of the system to the perturbation.This particular form is used to study the interplay between charge and spin orders in Sec.III D. In the presence of pinning fields, we preserve symmetry in the y direction by using PBC.In these cases, we can average over i y the staggered spin density and hole density to help improve statistics in detecting the spin and charge orders.The effect of the edge pinning field can be minimized or removed by separate calculations using smaller values of v.A robust order shows no variation with v except locally at the edge [39]; this can be used in combination of additional calculations with increasing system sizes, especially for determining the long-range behavior.When a full pinning field is applied to the simulation cell, such as in Eq. ( 5), an extrapolation must be performed to P → 0 to obtain the correct response.
B. Controlling the sign problem at finite temperature: self-consistent constraint In AFQMC the simplest trial density matrix we use is the non-interacting type, defined by the following trial Hamiltonian [51,53] where µ i,T needs to be tuned such that the electron density for this single-particle Hamiltonian is equal to that for the original many-body Hamiltonian.In this paper, we set µ i,T = µ T for simplicity and preserve translational invariance.This type of trial density matrix is employed in our first approach discussed above in Sec.II A, namely in fully periodic calculations.
To improve the CP approximation systematically, a self-consistent constraint can be employed, by coupling the AFQMC calculation with an effective independentparticle (IP) calculation [53,68].If we choose unrestricted Hartree-Fock (UHF) as the form for our IP calculation, the IP Hamiltonian that corresponds to the manybody one in Eq.( 1) takes the form where σ denotes the opposite spin of σ.The SU(2) spin rotational symmetry is reduced to U(1) symmetry because we assume the z-axis as the quantization direction.
The effective chemical potential µ eff tunes the particle number in the grand canonical ensemble.The effective interaction strength U eff , as well as the mean-field spin densities, { n iσ }, are to be determined through the selfconsistent iteration with the AFQMC calculations, as discussed further below.The self-consistent constraint is applied in our second approach discussed in Sec.II A, namely with simulation cells in which a symmetry-breaking pinning field is applied.The same external potential H s or H c , in Eq. ( 4) or (5), is applied to both the many-body Hamiltonian for AFQMC and the effective Hamiltonian in the IP calculation.With a pinning field applied at the edge, as mentioned earlier, the symmetry-breaking allows us to use the behavior of the many-body local "order parameters" (such as the one-body spin or density) versus distance as a probe to characterize the two-body correlations in the corresponding translationally invariant system.When a small external potential, e.g., the one in Eq. ( 5), is applied to the whole simulation cell, it allows us to do linear response study of the order built into the potential.
We start the self-consistent procedure with a trial density matrix e −τ H T , where τ ∈ (0, β) is a variable denoting the imaginary time where the constraint is applied [51,53], and H T can be taken as a simple mean-field Hamiltonian, e.g., the non-interacting or UHF Hamiltonian.Using this trial density matrix as a constraint, an AFQMC calculation is carried out (which always uses the physical β and U values), and spin-resolved electron densities { n iσ QMC } are computed.We then search for the best IP density matrix e −β eff HIP(U eff ) which produces spin densities { n iσ IP } closest to the AFQMC result.In other words, we minimize as a function of the effective interaction and effective temperature in the IP calculation.The IP solution with optimal U eff and β eff determined from Eq. ( 10), including the value of U eff and the resulting spin densities, is then fed as the input trial Hamiltonian to the new AFQMC calculation.We perform this procedure iteratively until (U eff , β eff ) have converged.
Before illustrating the self-consistent procedure with an example, we comment on two technical aspects.First, instead of using the IP spin densities in the trial Hamiltonian, we could also feed directly the AFQMC results from the previous iteration [68]; however, the statistical error of the QMC results could affect their accuracy as a constraint, and reducing the error bars during the self-consistency would incur additional computational cost.The introduction of the effective inverse temperature β eff allows the IP calculation to reproduce the QMC densities better.Second, we find it convenient and efficient to perform extensive UHF calculations prior to the self-consistent procedure, to set up a "look-up" table of the spin densities on a two-dimensional grid of (U eff , β eff ) values.In the UHF calculations, we use multiple initial conditions, including randomized spin densities, to facilitate convergence to the global minimum [71].
We illustrate the self-consistent scheme in Fig. 1.A 32 × 4 lattice is used, with spin pinning fields applied to its left edge, at ρ = 0.875 and U = 8, targeting an inverse temperature β = 18.We start the iteration with the UHF trial density matrix corresponding to the physical parameters.At U eff = 8 and β eff = 18, the UHF solution at the thermodynamic limit is a diagonal stripe state [71].This can sill be seen in the first step AFQMC result, with some frustration, as shown in the upper color sketch in panel (a).The UHF solution also greatly over-estimates the strength of the spin order, as seen in the curve with blue circles.It takes approximately four iterations between AFQMC and IP calculations for the crossover to take place from diagonal to vertical stripes in the trial density matrix, after which the self-consistent process continues to systematically remove the overestimation of spin order in the trial Hamiltonian, until convergence at the solution given by the curve with diamonds (spin pattern illustrated in lower color sketch in panel (a)).In panel (b), we show the trajectory of the self-consistency iteration for the IP calculations on a heatmap of χ 2 (U eff , β eff ) as defined in Eq. (10).Such a self-consistency procedure is performed for each system at each temperature, in all our calculations of the second type (see Sec. II A).

III. RESULTS
In this section we present our main findings, first providing an overview of the property of spin correlations in real space (Sec.III A), followed by presenting the picture in momentum space with spin susceptibilities and connecting it to real-space signatures (Sec.III B), across a wide range of temperatures.As mentioned, we systematically investigate three doping levels δ = 1/5, 1/8 and 1/10, i.e., electron densities of ρ = 0.8, 0.875, and 0.9, respectively.Then in Sec.III C we focus on the case of δ = 1/5, which is shown to have only short-range order, to examine the temperature evolution of the spin and charge correlations in this physical regime, as well as to illustrate the delicate requirement computationally to determine short versus long-range correlations.Finally in Sec.III D we discuss in detail the case of δ = 1/8, to probe the relation between spin and charge orders and the transition temperatures for long-range ordering.

A. FINITE-TEMPERATURE SPIN CORRELATIONS IN REAL SPACE
We first investigate the development of spin-spin correlation as a function of temperature in supercells under fully periodic conditions (i.e., without pinning fields).The supercell size is chosen to depend on δ to ensure commensuration with the predicted/expected wavelength, 2/δ, of the stripe or SDW state, should such an order exist in the ground state [28,39].In Fig. 2, we use a 20 × 8 supercell for δ = 1/5, which can accommodate two complete SDWs, a 16 × 8 supercell for δ = 1/8 and a 20 × 8 supercell for δ = 1/10, which can accommodate one complete SDW.We have performed calculations by varying the supercell size and shape, some of which are discussed below or presented in the appendices, to exam-ine the effect of breaking C 4 symmetry and ensure the robustness of the results.In all three cases, we see that short-range AFM correlations are immediately present.As T is lowered, the correlation length increases, and the magnitude of C s ( ) saturates quickly for small separation , and it grows more slowly at larger towards its lowtemperature limit.A modulation on the spin-spin correlation, signaled by the π-phase shift in the staggered spin correlation function (i.e., parallel spins across the domain boundaries), is seen almost as soon as the magnitude of the correlation near the domain boundary exceeds the statistical resolution.
The behaviors at δ = 1/8 and 1/10, as shown in the middle and right columns of in Fig. 2, are qualitatively similar.At an inverse temperature of β = 6, the latter displays a larger AFM correlation length, as further shown in Sec.III B. This is consistent with a stronger effective interaction due to smaller doping (U = 8 in both systems) and a stronger tendency for AFM correlation.As we lower the temperature to β = 9, a clear nodal line develops in δ = 1/8, showing two domains of AFM correlation, while the δ = 1/10 system remains in a single domain.At the very low temperature of β = 30, both systems show two well formed domains in the spin correlation, consistent with the expected filled stripe state at zero temperature.The size of the central AF domain in δ = 1/8 fluctuates between 7 ∼ 9 lattice spacings, consistent with the ground-state stripe state [28].In δ = 1/10, it is somewhat larger than the expected value of 10 from ground-state results [39], likely due to either residual finite-size effects or the use of the simplest, noninteracting uniform trial density matrix as our constraint in these particular calculations.Comparing the last two temperatures, the central domain appears to "breathe" at β = 9, especially in the case of δ = 1/10, being a slightly larger prior to the full formation of the outer domains at very low temperatures.We will examine the spin-spin correlation as a function of temperature more quantitatively with δ = 1/8 in Sec.III D.
Turning to the case of δ = 1/5 in the left column, we see that the short-range spin-spin correlation is weaker, but a clear modulated AFM pattern is visible even at β = 6, in contrast with the other two cases.The main reason for this is that the higher hole concentration significantly reduces the modulation wavelength, such that it can be comparable to the correlation length even at rather high temperatures.As we lower the temperature to β = 9, the size of the central AFM domain does not change but the adjacent domain grows.At β = 30, the spin correlation shows a strong tendency to restore C 4 symmetry (within the confines of the rectangular shape of the supercell).In fact the correlation seems to show, rather than a superposition of two linear waves, more of an isotropic pattern which would be compatible with short-range order.In order to reliably distinguish shortversus long-range orders in the 2-D limit, it is important to have a more accurate and systematic characterization of the correlation, using large and wide simulation cells, as we illustrate in the appendix and in Sec.III C, where it is shown that δ = 1/5 indeed exhibits only short-range order and is fundamentally different from the cases with smaller doping.

B. SPIN SUSCEPTIBILITY, AND COMPARISON WITH REAL-SPACE SIGNATURES
The transition from commensurate AF order to incommensurate (IC) SDW order can be characterized by the spin susceptibility χ s (k x , k y ) in momentum space, in a complementary manner to the real-space correlation functions.A shift of the peak value away from (π, π) to ((1 − δ)π, π) and (π, (1 − δ)π) is a signature of the onset of the SDW.However, finite-size effects can affect the detection, either by shifting the incommensurate peak position due to size resolution, or removing some of the peaks when rectangular supercells are used.Here we analyze the development of spin order as T is lowered using the spin susceptibility χ s (k x , k y ), and relate and contrast it with the real-space correlation function.
For δ = 1/8 both the real-space and momentum-space observables provide a clear signal of the SDW evolution.As shown in the inset of Fig. 3(b), the staggered spin cor-   relation function (−1) x C s ( x , 0) is zero within error bars for all sites x ≥ 5 at β = 6.When T is lowered to β = 8, the staggered spin correlation develops a small but clearly resolved signal at larger distances, with a change of sign at x ∼ 5.This creates a node separating two AFM domains and a π-phase shift across the domain boundaries.As we continue lowering T down to β = 30, the magnitude of the correlation function at long distance increases.However, the position of the node only changes slightly, from site-centered to bond-centered.From the main graph in Fig. 3(b), we see that the peak position of χ s (k x , k y ) shifts from (π, π) to (7π/8, π) starting from β = 11.At β = 30, the peak at (7π/8, π) is dominant, consistent with a long-range IC SDW order in the ground state.
The results at δ = 1/10 provide an instructive comparison.At β = 5, the system only displays short-range AFM order, with larger correlation length than in δ = 1/8, as can be seen from the inset in Fig. 3(c) and comparison with the inset in panel b.In momentum space, the peak values of χ s (k x , k y ) are substantially larger than in δ = 1/8, again consistent with stronger AFM correlations.When we lower the temperature to β = 10, a node appears between x = 7 and 8 in the staggered spin-spin correlation.The spin susceptibility χ s (k x , k y ), on the other hand, splits its peak at (π, π) and shifts to (9π/10, π) between β = 20 and 30.We argue that the lack of precise correspondence between these different signatures is a feature of such calculations, both because of the correlation being algebraic or weaker, and because of finite-size effects which become more severe for lower doping.At very low temperature β = 30, the position of the node in the real-space correlation shifts back to between x = 6 and 7. Note that the spin correlations in real-space at high T could have two subtly different behaviors: a modulation is always present even when the correlation is decaying exponentially; or it is a purely AFM correlation with no nodes until a crossover occurs to an SDW state [48].In momentum space both scenarios would show a signal of the spin susceptibity shifting from (π, π) to incommensurate peaks.Despite the small but discernible shrinking of the central AFM domain with T mentioned above, our results are more consistent with the former physical picture, namely of the modulation appearing as soon as the correlation length is such that a node can be accommodated.
We next contrast these results with the case of δ = 1/5, in Fig. 3(a).The magnitudes of the peak in the spin susceptibility are much lower than for δ = 1/8 and 1/10, indicating weaker AFM correlations as mentioned earlier.At high temperature β 6, the peak location starts to shift from (π, π) to (4π/5, π), as nodes in the real-space staggered correlation develop roughly in sync.These signatures both point to the development of modulated AFM correlation.However, without careful finitesize scaling analysis, neither can exclude the possibility of only short-range order, which turns out to be the case at δ = 1/5 and U = 6, as we show in the next section.

C. SHORT-RANGE SPIN AND CHARGE ORDERING IN THE OVERDOPED REGIME
Compared to lower doping levels, the δ = 1/5 system has the clearest signal of the modulated AFM correlation, which shows up at rather high T and requires relatively small system size to detect.At low temperatures the magnitude of the spin correlation function C s ( ) remains substantial up to | | ∼ 10.It is then tempting to conclude that this system exhibits prototypical SDW or even stripe order.However, a more detailed characterization of the spin and charge correlations reveals a different picture.
We apply our second approach, by introducing a local symmetry-breaking field in a cylindrical simulation cell and measuring the induced local order parameter.The spin and hole densities are used as proxies for the corresponding correlation functions, which magnifies the signal.Furthermore, this approach allows us to employ the FIG. 4. The nature of spin and charge orders with δ = 1/5 and U = 6.The results are obtained with our second type of calculations.A spin pinning field is applied to the left edge ( x = 0) of a 20×8 supercell.The top panel main graph shows the rung-averaged staggered spin density Sz( x), at five temperatures.The bottom panel shows the corresponding hole density (note vertical shift for clarity).The inset in (a) shows an illustration of the spin and hole densities in the supercell at the lowest temperature T = 1/20.The color of the arrows indicates the direction, while the length is proportional to the magnitude of Sz( x, y ).The diameter of black circles is proportional to the hole density h( x, y ).
self-consistent scheme described in Sec.II B, which provides a better constraint in AFQMC and leads to more accurate and robust results.We use an N = 20×8 lattice with cylindrical symmetry (OBC in x direction and PBC in y direction).The use of OBC in one direction can help probe correlations at larger distances than PBC, if the boundary effect is more localized.Spin pinning fields with amplitude v = 0.1 are added on the left edge of the cell, as described in Sec.II A. Both the spin pinning fields and boundary condition break symmetry and favor the unidirectional spin-and charge-density waves to grow along the cylinder.As mentioned, the combined use of this and our first approach, together with finitesize checks varying the shape and size of the cell, helps to remove potential bias in the results.We note that even a width-6 cylindrical cell shows significant finite-size effects which would suggest different physics from widths 8 and 10, as illustrated in Appendix C, so it is crucial to reach sufficiently large simulation cells.Figure 4 shows how the spin and hole densities evolve as temperature is lowered from β = 6 to β = 20.A modulated AFM spin pattern is seen in response to the local spin pinning fields.The magnitude grows with inverse temperature β at short range, and displays a non-zero signal at x ∼ 10, consistent with the results from C s ( ) discussed earlier.At long range, the magnitude of the spin pattern saturates quickly with β, and the rung-averaged staggered spin density Sz ( x ) remains zero within statistical errors.As we will see in the next section, this is in sharp contrast with the behavior in δ = 1/8, where a long-range stripe order develops at T = 0.
In the charge channel, when a long-range charge stripe exists, holes accumulate on the edges of AFM domains, which leads to peaks in hole density at the nodal points of the staggered spin density.In Fig. 4 (b), we observe a peak at x = 2 which corresponds to the first node of the modulated AFM spin pattern.However, no evidence is seen for further peaks in the hole density at other nodal positions.(A peak also appears at the opposite edge of the cylinder, which is induced by the OBC.We have performed systematic checks with separate calculations using longer cylinders to gauge and remove its effect.)This indicates that the charge order at δ = 1/5 is extremely short-ranged.We conclude that the spin and charge correlations remain short-ranged at δ = 1/5, U = 6 at all temperatures, consistent with the phase diagram from ground-state AFQMC calculations [39], an inhomogeneous DMFT study [34] and the recent higher-T CDet diagrammatic Monte Carlo results [50].It should be emphasized that our results are for U/t = 6, to allow direct comparisons with other work [49].Increasing the interaction strength will eventually drive the shortrange spin and charge orders to become long-ranged in the ground state, with U/t ∼ 8 at the border separating the two regimes [39].

D. INTERPLAY BETWEEN LONG-RANGE SPIN AND CHARGE ORDERING
We now investigate systematically the spin and charge orders, and their interplay, in the case of δ = 1/8 at U = 8, applying our second approach of symmetrybreaking pinning fields.In order to gauge finite-size effects and deduce the properties at the thermodynamic limit, we perform large-scale calculations on a sequence of supercell sizes, N = 32 × 4, 32 × 6, and 32 × 8, reaching very low temperatures.As shown in Fig. 5(a), at high temperatures (β = 3 and β = 6), the staggered spin density displays short-range antiferromagnetic correlations without an SDW modulation.When T is lowered, the correlation length grows along with the size of the AFM domain.Correspondingly, in Fig. 5(b), the distribution of hole density is uniform at these temperatures except at the ends of the lattice, adjacent to the pinning fields or to the boundary.
When T is lowered to β = 10 and β = 12, multiple AFM domains form and the staggered spin density displays a π-phase shift at domain wall boundaries.The SDW wavelength is 2/δ = 16 and that of the CDW is 1/δ = 8.At intermediate temperature, the SDW modulation starts from both open ends and extends to the central part of the simulation cell.As shown in the inset of Fig. 5(a), as T is lowered, the spin correlation length grows and the SDW modulation develops as soon as the correlation length becomes larger than the size of an single AFM domain.This is consistent with the observation from spin-spin correlation functions, as discussed in Sec.III B.
A distinguishable AFM domain in the center starts forming at β = 12, but has a size smaller than the expected periodicity (half a wavelength).Meanwhile, the wave-like modulation in the charge channel is very weak.Peaks in the hole density can only be observed at the spin nodal positions closest to the two ends ( x = 4 and x = 27), and have very small amplitudes.At β = 25, two complete incommensurate spin density waves form with the expected stripe periodicity.Furthermore, a modulation in hole density appears which leads to a charge density wave, with the hole density peak locations coinciding precisely with the spin nodes.The results at β = 25 are fully consistent with the ground state properties in this system [39], and provide a smooth "handshake" to T = 0.
Charge order is seen to follow spin in their temperature evolution in the calculations shown in panels (a) and (b).There is no charge order before the appearance of the first spin node at β ∼ 10.However, the signature of the charge order, namely the oscillation of the hole densities, is very small at intermediate temperatures where the transition takes place, which makes it challenging to extract more quantitative information.In order to further investigate the interplay between spin and charge orders, we next apply a periodic charge potential as a perturbation in a periodic supercell.We perform two kinds of calculations: (i) using this perturbation in the presence of the spin pinning fields on one edge, i,e., applying H s + H c ; (ii) using only H c without the finite local spin pinning field.
The results of these calculations are summarized in Fig. 5(c) and (d).We define the strength of the induced charge-density-wave as (∆n , where max x and min x are an expected location for the maximum and minimum density, respectively, according to the minimum and maximum of the applied potential in H c (at max x = 8, 16, 24 and min x = 12, 20 in the Fig. 5(b)).We exclude sites that are near the edges, which eliminates the effect from the spin pinning field in calculations of (i) (while having no effect in calculations of (ii)).Calculations are performed for a sequence of field strengths, P , to obtain (∆n) P .The results are shown in Fig. 5(c) for a set of 8 temperatures using (i), and 2 temperatures using (ii).Extrapolations are performed to the limit of vanishing amplitude of the periodic potential, P → 0, using a linear fit.In (i), we find that the resulting (∆n) P →0 is approximately zero above T /t ≈ 0.1, but becomes nonzero and increases rapidly at lower T .On the other hand,   in (ii) when there is no spin pinning, (∆n) P →0 remains zero at a temperature as low as β = 50.(Note that at high temperature (β = 3), the presence of spin pinning has little effect, and the two sets of calculations yield similar results, given by the blue lines.)These calculations allow us to disentangle the charge order from spin order.The results show, rather unambiguously, that the charge order is driven by the spin correlations.
We next examine the transition temperature of the charge ordering.Fig. 5(d) plots the extrapolated charge order (∆n) P →0 from calculations in (i), as a function of temperature, for three supercell sizes.We vary the width of the supercell from 4 to 8, while keeping the length fixed at 32, which is sufficient to accommodate four (two) full waves of the charge (spin) order in the ground state [39].The extrapolated charge response remains zero at high T for all system sizes, but turns finite abruptly at T ∼ 1/10.(In contrast, we have verified that the corresponding calculation for spin, applying a periodic spin perturbation for inducing an SDW in a 64 × 4 periodic supercell, leads to no extrapolated spin order down to T = 1/20; see Appendix F.) The response at low T becomes stronger as the system size is increased.The results from the three sizes fit well a linear form in 1/N (i.e., 1/4, 1/6, and 1/8), and we indicate the estimate for L y → ∞ from this fit by the dashed line in Fig. 5(d) (See Appendix C for more details on the extrapolation).These results suggest a finite-temperature phase transition into a phase with long-range (or quasi long-range) charge stripe ordering, with critical temperature T c ∼ 1/10.A more systematic finite-size scaling analysis, however, is required in order to determine T c more accurately, which is extremely challenging given the long-wavelength nature of these orders.The precise nature of how the spin and charge correlations evolve from high temperature to the ground-state long-range order is subtle.Our results show that the spin order is the driving force, yet a finite T c is seen for a transition to charge ordering.We comment on the possible nature of this transition in the next section.

IV. DISCUSSION AND CONCLUSION
Determining the properties of the doped 2D Hubbard model is faced with enormous challenges, that continue to motivate the development of ever more powerful and accurate computational methods.At zero temperature, intertwined competing orders are only separated by tiny energy scales, of the order of 10 −2 t or less, which can now be resolved by advanced wave-function based methods [28].At finite temperature, progress has been achieved using a variety of methods such as cluster extensions of DMFT [41,49,[72][73][74][75], diagrammatic Monte Carlo [50,76], DQMC [46,47], as well as METTS [48].Nonetheless, limitations in the range of temperature that can be accessed and/or in the system sizes that can be studied have prevented a comprehensive physical picture of the crossovers or transitions in the spin and charge correlations down to T = 0.
In this work, we have overcome these limitations to realize a full 'handshake' between finite-T and groundstate calculations, by deploying the latest advances in the AFQMC approach.This includes in particular the development of a self-consistent constraint with an effective finite-temperature mean-field performed by using both an effective interaction U eff and an effective inverse temperature β eff .This new approach allows us to reach supercell sizes and temperatures which were previously inaccessible, namely lattices of many hundred sites and temperatures as low as 1/50.The ability to perform computations which can span the full range of temperatures and connect with ground-state methods adds an important dimension to the accurate treatment of strongly correlated systems, and will be useful in many other systems besides the Hubbard model.
Systematic and meticulous checks were carried out by testing different forms of trial density matrices, using supercells with different geometries and boundary conditions and performing calculations both for translationally invariant systems and under symmetry-breaking external fields.As a result, we have obtained an accurate and systematic description of the physical nature and Tdependence of the spin and charge correlations over an extended range of temperature at three representative values of doping and interaction strength: δ = 1/5 at U/t = 6 and δ = 1/8, δ = 1/10 at U/t = 8.
We showed that, at the largest doping δ = 1/5 and U/t = 6, SDW correlations start to develop at rather high T , but remain short-ranged as T is lowered.Correspondingly, no charge ordering is found except for very shortrange correlations.These results are consistent with the very recent findings of Ref. [50], but much lower T was reached in the present study.At the two smaller doping levels and U/t = 8, spin correlations first develop at the antiferromagnetic wave-vector (π, π) at high-T and an incommensurate SDW is found with lowering temperature, as soon as the correlation length is sufficient to accommodate a node.These results confirm the qualitative expectations from mean-field theory [19][20][21][22][23] and are broadly consistent with recent computational results [48,50].We find that, at low temperatures the spin and charge correlation have modulation wavelengths compatible with the ground-state results at these parameters, i.e., filled stripes.In each case our results extend to very low temperatures and accomplish a full handshake with T = 0 properties.
Crucially, our work establishes how this growth of SDW correlations eventually connects to the groundstate stripe order at T = 0 [28,39].We identify a finite-temperature phase transition below which charge ordering sets in.This transition takes place at a fairly low temperature T c /t < 0.1 and was therefore inaccessible to previous studies.It will be valuable to confirm our result with more systematic finite-size scaling studies, which will require further algorithmic advances.By different calculations that impose bulk periodic potentials and/or boundary pinning fields, we were able to probe the different charge and spin responses to study their interplay.We establish that charge order is driven by the development of spin stripe correlations.
The nature of the critical behavior at the charge ordering transition poses questions of broad theoretical interest.This is extremely challenging to resolve numerically, because of the enormous length scale needed.However our results do place some bounds on possible scenarios.At a given commensurate doping δ = 1/q, there is a discrete symmetry Z q corresponding to the translations of the stripe pattern with wavelength q, e.g.Z 8 for 1/8-doping.It is interesting to note that the critical behaviour of the two-dimensional Z q -symmetric q-state clock model has a non-trivial dependence on q.While for q < 5 the phase transition corresponds to the breaking of the discrete Z q symmetry with long-range order for T < T c , an emerging U (1) symmetry has been suggested to occur for q ≥ 5, leading to a low-temperature phase with algebraic quasi long-range order and a BKT transition [77].Hence, one could imagine two different scenarios.The most straightforward one is that longrange charge order takes place below T c , while the spin correlation length remains finite at any non-zero T and long-range spin ordering takes place only at T = 0, in accordance with the Mermin-Wagner theorem [78,79].Another more intriguing scenario is that T c corresponds to a BKT transition below which the charge has algebraic quasi long-range order, and that the coupling between the charge and the spin degrees of freedom also leads to algebraically decaying spin correlations in the low-T phase.These questions are left open for a future study and in all likelihood will require and trigger significant new developments in computational methods.Another obviously crucial and still rather open question is the interplay between the charge and spin orders discussed here and superconductivity.For a mean-field Landau theory of the interplay between charge and spin ordering, see Ref. [80].
From an experimental standpoint, our work has implications in several directions.Experiments with ultracold fermionic atoms in optical lattices combined with siteresolved microscopy are able to image spin and charge correlations at progressively lower temperatures in a setup which directly emulates the Hubbard model with nearest-neighbor hopping only, as studied here [81][82][83].Indeed, in these experiments, growing antiferromagnetic correlations as temperature is lowered have been demonstrated in both the half-filled and doped systems [84], consistent with our finding that the onset of spin correlations takes place at a higher temperature than for charge correlations.Our determination of the temperature below which charge order occurs provides guidance for future experiments, which may be able to address this issue soon thanks to progress in cooling techniques.
Our work also has connections to experimental results on the cuprate materials.Indeed, for cuprates that display long-range stripe ordering such as Nd-LSCO [85] and LBCO [86], spin long-range order develops at a temperature below that of the charge ordering [87] in agreement with our finding.Note however that the stripes seen in our results are 'filled', while in La-based cuprates, they are found to be 'half-filled' [88].This difference, as well as other qualitative differences, points to the importance of the next-nearest neighbor hopping t .Extending our work to include this important parameter, as well as considering more realistic models such as a three-band Hubbard model including oxygen states are extremely valuable directions in which to extend the work presented here.Each UHF solution is used to form the input trial density matrix HT for the next step AFQMC calculation in the iteration.The trial state starts from an incorrect diagonal stripe order and converges to the correct order from the self-consistent coupling to AFQMC.The system is 32 × 4 with δ = 1/8 and the amplitude of edge pinning fields is v = 0.1.
width-4 cylinder can give qualitatively different physics [71,90].Width-6 cylinders are often much better; however in the case of δ = 1/5 it turns out that even a width-6 cylinder is insufficient for predicting the correct order in the TDL.In Fig. 8, we show the temperature evolution of spin order at δ = 1/5, U = 6 on a 30×6 lattice.In contrast to δ = 1/8, antiphase regions form even at temperatures as high as β = 6.The width of each AFM region fluctuates from ∆ x = 3 to ∆ x = 7 sites in different rows.As T is lowered to β = 12 and 15, more AFM regions appear but the amplitude of spin correlation at long distance is small.At T = 0.05, the sizes of the third AFM regions (symmetric about the center) tends to saturate to that of the central one.These behaviors are similar to the results in Ref. [49], in which a 16 × 4 cluster is used in DCA, and would suggest that there is no apparent distinguishing features between them and those in δ = 1/8 or 1/10 (except that the signature for AFM modulation or stripe appears more easily, at higher T ).
However, a different picture emerges as we increase the width of the simulation cells to reach the TDL.In Fig. 9, we use supercells with fixed length L x = 20 but different widths L y = 6, 8 and 10.To determine the spin order more accurately, we add AFM spin pinning fields to the simulation cell and measure the one-point function (the rung-averaged spin density).For each temperature, we employ the self-consistent scheme and plot the final converged results.On the width-6 cylinder, the result is compatible with that of the correlation function shown in Fig. 8.
In sharp contrast, when we increase the width of the cylinder to L y = 8, 10, we find that the spin order is short-ranged.The high-temperature behavior of the spin density on width-8 and width-10 cylinders is the same as that on the width-6 cylinder.As T is lowered, the spin density at large distance ( x ≥ 8) is essentially zero with error bars as shown in Fig. 9 (b) and (c).This is consistent with the conclusion from ground-state calculations [39].
In Fig. 10, we show the size dependence of the hole density in response to an external periodic perturbation which couples to the density.Two representative amplitudes of the periodic charge potential are shown, P = 0.02 and P = 0.1.We select two temperatures (β = 9 and 10) near the estimated phase transition and one temperature (β = 18) deep in the charge ordered phase, based on the analysis in Fig. 5.At higher temperature (β = 9), where charge order is absent, we see a negligible size dependence, as the charge response is essentially an independent-electron response to the external potential.As we lower T and approach the transition, the finite-size dependence becomes different for different perturbation strengths, showing a larger correction at P = 0.02.This trend becomes even stronger at T = 1/18 (well below the transition), when the size dependence remains small at P = 0.1 while a large increase as T is lowered, at δ = 1/5, on a 30 × 6 lattice with U = 6.PBC is used in both directions.AFM regions are observed, and the results in this finite system do not distinguish shortvs.long-range order.is seen in the amplitude of the induced CDW with increasing supercell size at P = 0.02.In Fig. 11, we show the finite-size extrapolation of the charge order (amplitude of the charge response to an external charge perturbation), under the effect of a local spin pinning field at the edge of the supercell.As mentioned, the dimension L x of the supercell is fixed at L x = 32.The finite-size data at each T fit very well a linear function of 1/L y .The response at low T becomes stronger as the width of the lattice increases.The behavior of the response as a function of lattice size is consistent with a transition to long-range charge order at a finite temperature.A more systematic finite-size scaling analysis to determine T c precisely would require further knowledge of the system and additional computations in significantly larger system sizes.We have discussed in the main text the two different approaches to probe spin and charge correlations, namely fully PBC calculations versus pinning field calculations.We have illustrated how they complement each other.It is reassuring that they lead to consistent results on the spin and charge order, even with very different settings and under very different constraining density matrices.However, as we discussed, they can have different behaviors in finite-size systems and give different signal strengths.In this appendix, we show some comparison of the structure factors to supplement the discussions in the main text.
The equal-time structure factors for the spin and charge correlations discussed in the main text are defined by: In Fig. 12 (a), we present the spin structure factor at U = 8 and δ = 1/8.Lattices with two different sizes (16 × 4 and 32 × 4) are used.For both lattice sizes, the peak position of S s (k x , π) is seen to shift from (π, π) to (7π/8, π) and (9π/8, π).Similar to the spin susceptibility χ s in Sec.III B, it indicates the evolution of an incommensurate SDW order.
The signature for charge order is subtle and more challenging to detect, as shown in Fig. 12 (b).On a 16 × 4 lattice, the amplitude of charge structure factor peak at (±π/4, 0) remains extremely small, down to T = 0.This reflects strong finite-size effects, with the signal for long-range correlation being submerged in a large background.To detect the peak that characterizes the charge stripe, we need to use larger lattices, e.g 32 × 4. Then distinguishable but small peaks at (±π/4, 0) appear at T = 0.033.As T is lowered to 0.02, the peaks become more clear and eventually converge to the amplitude at ground state (T = 0).As a comparison, the second approach of applying a pinning field to turn the spin and charge density into a proxy for correlation functions shows much clearer signals.We apply Fourier transform to the one-point func- where the spin and charge densities are shown in the main text, for example in Fig. 5 (a) and (b).We perform the Fourier transform on a subset of the lattice, with sizes of 16 × 4 and 24 × 4, respectively, for the spin and charge densities.Sites on both ends are discarded to remove effects induced by spin pinning fields and the OBC.
As shown in Fig. 12 (c), as T is lowered, the peak of Ss (k) gradually shifts from (π, π) to (7π/8, π) and (9π/8, π), which signals the development of the SDW order.In contrast to the charge structure factor, the shoulder peaks in Sc (k) at (±π/4, 0) appear much more clearly, starting approximately at β = 15, consistent with the result from directly applying charge-inducing field discussed in the main text.Its height is much larger (≥ 10 times) than that from the charge structure factor.As mentioned, recent work (see, e.g., Ref. [17]) indicates that there is no d-wave superconducting long-range order in the ground state of the first two systems that we study, at δ = 1/8 or δ = 1/10 (U/t = 8).The third system, δ = 1/5 with U/t = 6, is near the boundary of the parameter regime explicitly scanned in Ref. [17]); indeed there has been a recent suggestion that there may be superconducting order in the ground state [38].We have computed in this system the d-wave pairing susceptibility χ d , which is defined as In Fig. 13(a), we show the result for three system sizes N = 20 × 4, N = 20 × 6 and N = 20 × 8.As T is lowered, the d-wave pairing susceptibility first grows and then decreases.At low temperature T < 0.1, as the system size increases, the pairing susceptibility decreases, which is incompatible with the development of long-range order.(This, of course, cannot rule out the possibility of a superconducting order developing below the temperature scale we investigate.)For reference, the corresponding spin susceptibility is shown in panel (b).As T decreases, χ spin (4π/5, π) increases (except for a small downward bend in 20 × 8).Recall that this system does not have long-range spin order as T → 0, as we showed in the main text.

Appendix F: Response to periodic spin potential
To complement our studies of the charge response in Sec.III D, we apply periodic spin pinning fields and examine the amplitude of the induced spin order in the limit of zero external perturbation in the linear response regime.Similar to what was done for the charge channel, we apply sinusoidal spin pinning fields to a 64 × 4 lattice at U/t = 8 and δ = 1/8.PBC is used in both directions of the simulation cell.We add the following to the original Hamiltonian where B is the amplitude and κ is the characteristic wavevector.In this system the wavelength of the SDW is twice that of the charge stripe, so a longer cylinder (64 × 4) is used to accommodate the same number of SDWs as for charge stripes in Fig. 5.In Fig. 14, we plot the row-averaged staggered spin density Sz ( x ) as defined in Eq. ( 6).We observe four complete spin dnesity wave periods, with positions of peaks and valleys consistent with the applied sinusoidal potential.As we increase the amplitude of the input spin perturbation, the amplitude of Sz ( x ) increases approximately linearly as a response.Note that the magnitude of the spin response is much larger than that of the charge, as seen in comparison with As shown in Fig. 15 (a), we then perform extrapolations to vanishing amplitude of the periodic spin potential, B → 0, using a linear fit.statistical error bars down to the lowest temperature of T = 0.05, consistent with Mermin-Wagner theorem.The contrasting behaviors between the spin and charge responses provide yet another strong consistency check on our computational approach.

2 FIG. 1 .
FIG. 1. Systematic improvement of the AFQMC accuracy using the self-consistent constraint.(a) The main plot shows the staggered spin density (−1) x +0 Sz( x, y = 0) vs. site x, obtained from AFQMC at five representative steps in the selfconsistent iteration.Statistical error bars are much smaller than symbol size.The insets show visualization of the realspace spin density obtained from AFQMC in the first and last iterations (corresponding to the circle and diamond in the second panel, respectively) .Blue and orange colors mark regions of AFM spin densities with a phase change between them.(White color means the magnitude of the spin density is smaller than twice the AFQMC statistical error bars.)(b) The trajectory of how (U eff , β eff ) evolve in the self-consistent procedure is shown on the heatmap of χ 2 (U eff , β eff ).Each solid symbol denotes the effective interaction and temperature parameters of the IP Hamiltonian at one step of the iteration, and the light blue arrows indication the direction in the iteration.The five selected steps shown in panel (a) are indicated here by their corresponding symbol shapes, while all other steps are shown as stars.The system is 32 × 4 with δ = 1/8 and U = 8, at inverse temperature β = 18.

FIG. 2 .
FIG.2.Evolution of the spin correlation function Cs( x, y ) as T is lowered (T = 1/6, 1/9, and 1/30 from top to bottom panels), for doping δ = 1/5, U = 6 (left column), δ = 1/8, U = 8 (middle), and δ = 1/10, U = 8 (right).These are obtained from our first type of calculations, using fully periodic supercells along both directions, with the correlation averaged over all reference sites i, indicated by the black cross.(The amplitude of Cs(0, 0) is not shown.)The color of the arrows indicates the sign (dark blue = + and purple = −), while the length of the arrows is proportional to the amplitude of Cs( x, y ).The color of the square at each site represents the sign of the staggered spin correlation (−1) x+ y Cs( x, y ).A missing arrow and missing colored square indicate that |Cs( x, y )| is less than twice its Monte Carlo error bar, which is considered approximately zero within our statistical resolution.

FIG. 5 .
FIG. 5.The nature of spin and charge orders with δ = 1/8 and U = 8.The top two panels have similar layout as in Fig. 4. The inset in (a) shows a zoom in the middle region with β = 10 (triangle) and β = 12 (square), when the first nodes of spin modulation develop.Panel (c) shows calculations to probe the response of the system to an applied external periodic charge potential with amplitude P .Solid lines are obtained in the presence of an additional finite edge pinning field for spin.In contrast, dashed lines show otherwise identical calculations, but without the spin pinning field.Panel (d) shows the development of charge order as T is lowered, based on the extrapolation procedure shown by the solid lines in panel (c).Results in panels (a), (b) and (c) are for a 32 × 4 supercell, while panel (d) shows three different sizes, 32 × 4, 32 × 6, and 32 × 8 as labelled, as well as an extrapolation to infinite transverse size.

FIG. 6 .
FIG.6.(a)-(f) Convergence of the trial state in the self-consistent process which couples it to the AFQMC calculation.The spin density in real space is shown for the UHF solution at selected (U eff , β eff ) values on the optimization path shown in Fig.1.Each UHF solution is used to form the input trial density matrix HT for the next step AFQMC calculation in the iteration.The trial state starts from an incorrect diagonal stripe order and converges to the correct order from the self-consistent coupling to AFQMC.The system is 32 × 4 with δ = 1/8 and the amplitude of edge pinning fields is v = 0.1.

25 FIG. 7 .
FIG. 7. Evolution of the spin correlation function Cs( x, y ) as T is lowered.The results are for a 16 × 16 square lattice with U = 8 and δ = 1/8 doping, at four T values as labeled (in β).

20 FIG. 8 .
FIG.8.Evolution of the spin correlation function Cs( x, y ) as T is lowered, at δ = 1/5, on a 30 × 6 lattice with U = 6.PBC is used in both directions.AFM regions are observed, and the results in this finite system do not distinguish shortvs.long-range order.

FIG. 9 .
FIG. 9. Investigating the behavior in the TDL at δ = 1/5 (U = 6).Rung-averaged staggered spin density Sz( x) is plotted vs. site x on cylinders with fixed length Lx = 20 and varying width.Staggered pinning fields with amplitude v = 0.1 are applied to the left edge, at ( x = 0).The three panels show widths (a) Ly = 6, (b) Ly = 8, and (c) Ly = 10, respectively.Converged results obtained from the selfconsistent procedure are plotted at each temperature.Error bars are smaller than the symbols size.

15 FIG. 10 . 18 FIG. 11 .
FIG.10.Finite-size analysis of the charge response to an external periodic charge perturbation.Rung-averaged hole densities h( x) are shown for two strengths (P ) of the external charge potential, at three temperatures spanning the transition Tc (see discussions in Fig.5).The systems are at δ = 1/8 and U = 8, with Lx = 32 and staggered spin pinning fields of amplitude v = 0.1 applied to the left edge at x = 0. Three lattice sizes are studied, with width Ly = 4, 6 and 8. Error bars are smaller than symbol size.
Appendix D: Spin and charge structure factors under the two different approaches

FIG. 12 .
FIG. 12. Spin and charge structure factors under the two different approaches, fully PBC versus pinning field calculations.The left column shows (a) spin structure factor Ss(kx, π) and (b) charge structure factor Sc(kx, 0) vs. momentum kx on 16 × 4 and 32 × 4 lattices with PBC applied in both directions, for two finite (low) temperatures, compared to the corresponding ground-state result.The right column shows the Fourier transforms of the (c) spin density Ss(kx, π) and (d) charge density Sc(kx, 0) vs. kx on segments of a 32 × 4 lattice, for a sequence of temperatures.U = 8 and δ = 1/8 are used in all calculations.

FIG. 15
FIG. 15.(a) Linear extrapolations for the induced SDW order ∆Sz(B) versus the external potential strength B for six different temperatures, for the same system as in Fig. 14.(b) The absence of spin order as T is lowered, based on the extrapolation procedure shown in panel (a).