Stripes, Antiferromagnetism, and the Pseudogap in the Doped Hubbard Model at Finite Temperature

The interplay between thermal and quantum fluctuations controls the competition between phases of matter in strongly correlated electron systems. We study finite-temperature properties of the strongly coupled two-dimensional doped Hubbard model using the minimally-entangled typical thermal states (METTS) method on width $4$ cylinders. We discover that a phase characterized by commensurate short-range antiferromagnetic correlations and no charge ordering occurs at temperatures above the half-filled stripe phase extending to zero temperature. The transition from the antiferromagnetic phase to the stripe phase takes place at temperature $T/t \approx 0.05$ and is accompanied by a step-like feature of the specific heat. We find the single-particle gap to be smallest close to the nodal point at $\mathbf{k}=(\pi/2, \pi/2)$ and detect a maximum in the magnetic susceptibility. These features bear a strong resemblance to the pseudogap phase of high-temperature cuprate superconductors. The simulations are verified using a variety of different unbiased numerical methods in the three limiting cases of zero temperature, small lattice sizes, and half-filling. Moreover, we compare to and confirm previous determinantal quantum Monte Carlo results on incommensurate spin-density waves at finite doping and temperature.


I. INTRODUCTION
Understanding the physics and phase diagram of copper-oxide high-temperature superconductors is arguably one of the most fundamental and challenging problems in modern condensed matter physics [1][2][3][4][5]. Especially fascinating is the normal state of the hole-doped materials, which displays highly unusual but rather universal features such as the formation of a 'pseudogap' at an elevated temperature [6][7][8][9][10][11] and, upon further cooling, the formation of charge-density wave order for a range of doping levels [12,13].
Early on in the field, this fundamental problem was phrased in the theoretical framework of the twodimensional Hubbard model [14][15][16]. Even though the degree of realism of the single-band version of this model for cuprates can be debated, it has become a paradigmatic model embodying the complexity of the 'strong correlation problem'. Establishing the phase diagram and physical properties of this model is a major theoretical challenge and a topic of intense current effort.
Controlled and accurate computational methods which avoid the biases associated with specific approximation schemes are invaluable to address this challenge, both because of the strongly interacting nature of the problem and because it is essential to establish the physical properties of this basic model beyond preconceptions. Recently, the community has embarked on a major effort to combine and critically compare different computational * awietek@flatironinstitute.org methods, with the aim of establishing some properties beyond doubt and delineating which questions remain open [17][18][19][20].
Among the many methods that have been developed and applied to approach the problem, it is useful to emphasize two broad classes. On the one hand, wavefunction based methods using tensor network representations and extensions of the density-matrix renormalisation group algorithm (DMRG) [21][22][23][24] have been successful at establishing the ground-state properties of systems with a cylinder geometry of infinite length but finite transverse size. On those systems, the ground-state has been established to display inhomogeneous 'stripe' charge and spin ordering [18,[25][26][27][28][29][30]. Substantial support for this picture has come from ground state quantum Monte Carlo methods with approximations to control the sign problem [31][32][33], density matrix embedding methods [34][35][36], and finite-temperature determinantal quantum Monte Carlo simulations [37,38]. For an early discussion and a review of stripes in the Hubbard model and cuprates see [39][40][41][42]. Furthermore, recent studies have established that these states are favored over a possible uniform superconducting ground-state, by a tiny energy difference, for the unfrustrated model with zero next-nearest neighbor hopping [19].
Cluster extensions of dynamical mean-field theory (DMFT) [43][44][45][46] on the other hand, address the problem from a different perspective. These methods use the degree of spatial locality as a control parameter and, starting from the high-temperature limit in which the physics is highly local, follow the gradual emergence of non-local physics as spatial correlations grow upon reducing temperature. These methods have established the occurrence FIG. 1. Hole-densities and spin correlations of a typical METTS state |ψi for U/t = 10 at hole-doping p = 1/16 on a 32 × 4 cylinder. The diameter of black circles is proportional to the hole-density 1 − n l = 1 − ψi|n l |ψi , the length of the red/blue arrows is proportional to the amplitude of the spin correlation S0 · S l = ψi| S0 · S l |ψi . The black cross indicates the reference site of the spin correlation. Red/blue squares indicate the sign of the staggered spin correlation (−1) x+y S0 · S l (a) T /t = 0.025. We observe antiferromagnetic domain walls of size ∼ 6 − 8 bounded by maxima in the hole-density. This indicates a fluctuating stripe phase realized. (b) T /t = 0.100. We observe extended antiferromagnetic domains. No regular stripe patterns are formed.
These two ways of looking at the problem (groundstate T = 0 versus finite-T starting from high-T ), leave the intermediate-T regime as uncharted territory and, crucially, leave a major question unanswered: how does the emergence of the PG at high temperature eventually evolve into the charge inhomogeneous ground states revealed by DMRG and other recent studies?
In this article, we bridge this gap and answer this question for the weakly doped Hubbard model at strong coupling, on a specific lattice geometry. We consider long cylinders of width 4 and length up to 32 sites: this is close to the current limit of ground-state DMRG studies, such as the extensive study recently presented in Ref. [30]. To study this challenging system at finite temperature, we have refined and further developed the minimally entangled typical thermal state (METTS) method [71,72]. This allows us to follow the full evolution of the system as the temperature is increased from the ground state exhibiting stripes. We monitor this evolution through the calculation of several observables, such as the spin and charge structure factors, the momentum distribution function, as well as thermodynamic observables. We determine the onset temperature of the ground state stripe phase. We find that as the system is heated above this onset temperature, a new phase with short-range antiferromagnetic correlations is found, which shares many features with the experimentally observed pseudogap phase. In contrast to the incommensurate correlations observed in the stripe phase, the magnetic correlations are com-mensurate in this regime. We identify the pseudogap onset temperature, which is distinctly higher than the temperature at which the stripes form.
Finite-temperature simulations using DMRG or tensor network techniques beyond ground-state physics were developed early on, but seemed practical mostly for one dimensional systems. The purification method [73,74] is particularly attractive for its simplicity [75][76][77], but in the limit of low temperatures, its representation of the mixed state carries twice the entanglement entropy of the ground state, making it unsuitable for wide ladders, although techniques have been developed to reduce the entanglement growth [78]. The METTS method [71,72] was developed to overcome this obstacle. To simulate finite temperatures several matrix product states (MPS) are randomly sampled in a specified way. The computation of a single such MPS (also called a METTS), which has entanglement entropy similar to or less than the ground state, can be performed with the same computational scaling with system size as ground-state DMRG [71,72,79]. However, the METTS algorithm requires imaginary time evolution rather than DMRG's more efficient sweeping to find the ground state, and also a certain amount of random sampling, both of which increase the calculation time significantly. It has been argued that for certain one-dimensional systems, where entropy scaling is less important than sampling error, the METTS algorithm does not outperform the purification approach [80]. In this manuscript, we develop and refine the METTS methodology and demonstrate that it can be successfully applied to study the finite temperature properties of the doped Hubbard model in the strong coupling limit on a width-4 cylinder.
The manuscript is organised in two major parts. In sections II to VI we discuss our physical results. We define the model, establish basic notations, and give a minimal explanation of the METTS method in section II. The METTS method allows investigating typical states at a given temperature, which we show in section III. Results on charge and magnetic ordering are then presented in section IV. We discuss the nature of the spin and charge gaps in section V and present results on the specific heat and magnetic susceptibility in section VI. In the second part in sections VII to XI we explain and demonstrate the more technical aspects of our simulations. We explain the details of the algorithms we employ for our simulations in section VII. We discuss the accuracy of the imaginary-time evolution algorithms we employ in section VIII and investigate the METTS entanglement section IX. We demonstrate that we can achieve agreement with different numerical methods within their respective limits in sections X and XI. Finally, we discuss the physical implications of our results and the perspectives opened by our work in sections XII and XIII. The statistical properties of the METTS time series are discussed in appendix A. There, we give the reason why these simulations can actually be performed at a reasonable computational cost. We demonstrate that the the variance of several estimators quickly decreases when lowering temperatures as well as increasing system sizes.

II. MODEL AND METHOD
We consider the single-band Hubbard model on the square lattice, where σ =↑, ↓ denotes the fermion spin, c † iσ and c iσ denote the fermionic creation and annihilation operators, and n iσ = c † iσ c iσ denotes the fermion number operator. The system is studied on a cylinder with open boundary conditions in the long direction and periodic boundary conditions in the short direction. We denote the length of the cylinder L, the width W , with the number of sites N = L × W . The cylindrical geometry is adopted since our computations are performed using MPS techniques. When studying the system using METTS, we employ the canonical ensemble with fixed particle number N p . Thermal expectation values of an observable O are given by, Here, β = 1/k B T denotes the inverse temperature and Z = Tr(e −βH ) denotes the partition function. We henceforth set k B = 1. Results are presented depending on the hole-doping, where n = N p /N denotes the particle number density.
We study the finite-temperature properties of the Hubbard model Eq. (1) on a L = 32 and W = 4 square cylinder in the strong-coupling regime at U/t = 10 and focus on the hole-doped case p = 1/16. We also performed simulations at half-filling for comparison. Simulations have been performed for a temperature range, To evaluate thermal expectation values, the METTS algorithm averages over expectation values of pure states Here, · · · denotes averaging over random realizations of |ψ i , which are called minimally-entangled typical thermal states (METTS). We discuss the METTS algorithm in detail in section VII A. While in principle, the basic METTS algorithm is unbiased and exact, the computation of the states |ψ i amounts to an imaginary-time evolution of product states, which is performed using tensor network techniques. As shown in section VIII, this operation can be performed to a high precision using modern MPS time-evolution algorithms. The accuracy is controlled by the maximal bond dimension D max of the MPS representation of each METTS |ψ i and by the number of METTS sampled. Our implementation is based on the ITensor library [81].

III. METTS SNAPSHOTS AT FINITE TEMPERATURE
The METTS |ψ i can be considered "snapshots" of the thermal state at a given temperature. It is, therefore, informative to study properties of individual METTS. The METTS naturally separate the fluctuations of the system into mostly-quantum and mostly-thermal. The mostly-quantum, higher-energy fluctuations are contained within individual METTS, while long-distance low-energy fluctuations tend to appear via the sampling over different METTS. At zero temperature, all METTS are identically the ground state. At infinite temperature, they are classical product states with each site randomly in one of the four basis states of a single site. We illustrate some typical METTS on the 32 × 4 square cylinder at U/t = 10 and p = 1/16 in Fig. 1. Hole-densities at site l, 1 − ψ i |n l |ψ i , where n l = n l↑ + n l↓ , are shown as black circles. Spin correlations, ψ i | S l · S m |ψ i , are shown as red and blue arrows. We choose a reference site l = 0, indicated by the black cross in the middle of the lattice. For states at temperature T /t = 0.025, like the state shown in Fig. 1(a), we observe regular charge density wave patterns with a wavelength of 6 ∼ 8 lattice sites. This charge modulation is accompanied by the staggered spin correlations changing their sign, as indicated by the blue and red squares. We, thus, observe antiferromagnetic domains of the size of the charge wavelength. The Sc(k), of the 32 × 4 square cylinder at U/t = 10 for p = 0 (left) and p = 1/16 (right). We compare different temperatures from METTS with Dmax = 2000 and ground-state DMRG with Dmax = 5000. (a) Magnetic structure factor for p = 0 and ky = π. The peak at k = (π, π) indicates the antiferromagnetism. (b) Magnetic structure factor for p = 1/16 and ky = π. The peak at k = (7π/8, π) (gray dashed line) indicates the stripe order illustrated in Fig. 1(a). (c) Charge structure factor for p = 0 and ky = 0. The quadratic behavior at kx = 0 indicates a gap to charged excitations. (d) Charge structure factor for p = 1/16 and ky = 0. We observe a peak at k = (π/4, 0) (gray dotted line). This indicates a half-filled stripe phase at low temperatures. The approximately linear behavior at kx = 0 indicates a small or vanishing charge gap. In all cases we find the METTS results converging towards the DMRG results in the limit T → 0.
number of observed charge stripes on the 32 × 4 cylinder is 4. Since hole-doping p = 1/16 corresponds to 8 holes on this lattice, we have two holes per stripe and thus observe a half-filled stripe on the width 4 cylinder. We show a typical METTS state at temperature T /t = 0.100 in Fig. 1(b). At this temperature, we do not find regular charge density wave patterns. However, we observe enhanced antiferromagnetic correlations with larger antiferromagnetic domains.

IV. MAGNETIC AND CHARGE ORDERING
To quantify these observations, we investigate the magnetic structure factor, Here, r l denotes the coordinate of the l-th lattice point, and k = (k x , k y ) denotes the quasi-momentum in reciprocal space. The width W = 4 cylinders we focus on, resolve four momenta in the y-direction, k y = 0, ±π/2, π. The ordering vector of the half-filled stripes k = (7π/8, π) is shown in red. (b) Charge structure factor Sc(k) at stripe ordering vector k = (π/4, 0). A transition from stripe-order to short-range antiferromagnetic order takes place at T /t ≈ 0.05. We find agreement between simulations at different bond dimensions.
Furthermore, we investigate the charge structure factor, We show the magnetic structure factors with ymomentum k y = π for p = 0 and p = 1/16 for several select temperatures in Fig. 2(a,b). The METTS simulations shown here have been performed with D max = 2000. The convergence of these results as a function of bond dimension for specific values of k is shown in Fig. 3. Ground-state DMRG calculations have been performed to obtain results for T = 0, shown as the gray line. The DMRG results have been performed using 30 sweeps at maximal bond dimensions up to D max = 5000. We observe convergence in the displayed quantities. At halffilling, p = 0, in Fig. 2(a) we observe a prominent peak at wave vector k = (π, π), which corresponds to the antiferromagnetic ordering vector. We observe convergence of the structure factor at temperatures T /t = 0.050, 0.025 for T → 0 towards the DMRG results.
At hole-doping p = 1/16 in Fig. 2(b) we observe that at temperatures below T /t ≈ 0.050 the peak in the magnetic structure factor is shifted by δ = π/8 away from k = (π, π) towards k = (7π/8, π). This feature is a signature of the formation of stripes, where the shift in the antiferromagnetic ordering vector is induced by the antiferromagnetic domains of opposite polarization, as observed in Fig. 1(a). The accompanying charge modulation is quantified by the charge structure factor shown in Fig. 2(d). We observe a peak at wave vector k = (π/4, 0) = (2δ, 0). Hence, the charge modulations occur at half the wavelength of the antiferromagnetic modulations. We again find that the METTS results tend towards the T = 0 DMRG results as T → 0.
Above T /t ≈ 0.050, we notice that the peak in the magnetic structure factor is shifted back towards the antiferromagnetic ordering vector k = (π, π), as shown for T /t = 0.100. We show the temperature dependence of the magnetic structure factor for both ordering vectors in Fig. 3(a). We find that above T /t ≈ 0.050, k = (π, π) is the dominant ordering vector with a maximum in the structure factor attained at T /t ≈ 0.100. Below T /t ≈ 0.050, the stripe ordering vector k = (7π/8, π) is dominant. This suggests a transition or crossover from the stripe phase, to a new phase with antiferromagnetic correlations. The realization of the stripe phase at low temperatures is also apparent in the behavior of the charge structure factor at ordering vector k = (π/4, 0) in Fig. 3(b). We observe a sharp increase below, T /t ≈ 0.050. We also compare the structure factors from METTS simulations at different bond dimensions D max in Fig. 3 and observe agreement within errorbars for most parameters. The peak height of the magnetic structure factor at k = (π, π) is slightly decreased for D max = 2000.

V. GAPS AND CORRELATIONS IN THE DOPED SYSTEM
We now focus on the charge and spin excitations in the hole-doped case, p = 1/16. Generally, the behavior of a particular type of correlation function is tied to the presence or absence of an associated gap. Consider the density structure factor in the vicinity of k = (0, 0) in Fig. 2(c,d). Its behavior differs significantly between the half-filled case in Fig. 2(c) and the p = 1/16 hole-doped case in Fig. 2(d). We observe that for p = 0, This behavior is expected for a system with a charge gap [33,[82][83][84][85]. The charge gap in question is defined (at T = 0) as where E 0 (N ↑ , N ↓ ) is the ground state energy with N ↑ and N ↓ up and down particles. A charge gap, in the thermodynamic limit, implies that the system favors that particular doping compared to states differing by two particles. Typically, this only happens at commensurate fillings (such as half-filling) where the filling is a smalldenominator simple fraction. We do not expect a doping of p = 1/16 to have a charge gap. For charge-gapless systems, the behavior of the charge structure factor close to k = (0, 0) is expected to be [33,83,84], We observe this behavior for hole-doping p = 1/16 in  (12) and (13). We find that for large cylinder lengths, the single-particle gap ∆c and the spin gap ∆s approach finite values ∆ We also define the single particle gap, and the spin gap, In Fig. 4 we show these three types of gaps for p = 1/16, evaluated using DMRG with fixed particle and spin quantum numbers on cylinders of length 8, 16, 24, and 32. We checked that the computed gaps correspond to bulk excitations, by investigating the density distribution upon adding and removing particles. This showed that additional particles or holes are indeed inserted in the bulk of the system. The single particle gap ∆ (1) c approaches a finite value ∆ (1) c /t ≈ 0.25, and the spin gap ∆ s approaches a finite value of ∆ s /t ≈ 0.07. In contrast, the charge gap decreases approximately as 1/L for L = 16, 24, 32, consistent with a vanishing charge gap in the L → ∞ limit, in agreement with previous DMRG results on the width 4 cylinder [30]. Our observations on the stripe phase agree well with the Luther-Emery 1 (LE1) phase in Ref. [30]. In particular, a finite single particle gap has analogously been reported.
The single particle gap is closely tied to the behavior of the momentum distribution function, as one varies k. Of course, our single particle gap computed from Eq. (12) is the minimum gap as one varies the momentum. Gaps at other momenta could be obtained from spectral functions, but these are beyond the scope of this work. However, we can observe very different behavior in n σ (k) for different k. Results for different k and also several temperatures are shown in Fig. 5. For a system with a Fermi surface, the momentum distribution function is expected to be discontinuous at the Fermi momentum [86]. Although for DMRG at finite length and bond dimension this discontinuity is usually not directly observed, the slope of the momentum distribution function as a function of bond dimension can be investigated to diagnose a Fermi surface [85]. We find the maximal slope of the momentum distribution function n σ (k) to remain finite for all values of k y , consistent with a single particle gap. The steepest slope of n σ (k) is observed at y-momentum k y = π/2. There, we observe that n σ (k) at different temperatures intersect at a specific value close to the nodal point k = (π/2, π/2). We do not observe such an intersection for k y = 0. The small slope for k y = 0 and k y = π suggests a large single-particle gap at these lines in the Brillouin zone, while the larger slope close to the nodal point suggests a smaller single-particle gap.
A spectral gap implies exponentially decaying groundstate correlation functions in real space [87], where a slow exponential decay is a signature of a small gap. We display the single-particle correlation function at p = 1/16 from DMRG after Fourier transform in the y-direction in Fig. 5

(b). The quantity shown is
where W = 4 denotes the width of the cylinder. Results are shown for k y = 0, π/2, π and as a function of bond dimension. We find the slowest decay is found for k y = π/2, with fast exponential decay for k y = 0 and k y = π. This demonstrates that the gap to charged excitations is smallest at k y = π/2.
It is interesting to compare the magnitude of the single particle gap and the spin gap, both with each other and with the temperatures associated with the onset of shortrange antiferromagnetic order and the onset of stripes.
(a) Specific heat, p = 0. Results from all maximal bond dimensions agree. We find a C ∝ T 2 behavior at low temperature, predicted from spin-wave theory of an antiferromagnet. We also show exact AFQMC data on the same system, which agrees within errorbars (b) Specific heat, p = 1/16. We observe a step-like feature around T /t = 0.05, where we locate the transition to the stripe phase in Fig. 2. In the range T /t ≈ 0.075 to T /t ≈ 0.175 we observe an approximately linear behavior. (c) Magnetic susceptibility at p = 0. We observe a maximum at T × /t ≈ 0.29. (d) Magnetic susceptibility at p = 1/16. A maximum is located at T * /t ≈ 0.25.
The onset of local antiferromagnetic correlations is associated with the peak in the specific heat and maximum of the uniform susceptibility, discussed in the next section, which is slightly above 0.2t. This is close in magnitude to the single particle gap, 0.25t, and so it is tempting to tie them together. One could consider a linkage of these energy scales through pair-binding. A simple picture of pair binding is that two holes together disrupt the local antiferromagnetic correlations less than two separate holes. This mechanism could only occur below the onset temperature of local antiferromagnetism.
It is also tempting to tie the spin gap (about 0.07t) to the onset temperature of stripes, 0.05t. However, this linkage is not so clear. The spin gap is probably strongly influenced by the finite width of the cylinder. At halffilling, in the Heisenberg limit, even-width cylinders have a spin gap which vanishes exponentially with the width, and the 2D limit is gapless. It is also not clear that a striped state should have a spin gap, so the similarity in the spin-gap and stripe energy scales for the N y = 4 system may be a coincidence.

VI. THERMODYNAMICS
Finally, we discuss basic thermodynamics quantities for half-filling and p = 1/16. The specific heat is given by, When using the METTS algorithm, we compute the numerical derivative of measurements of the total energy E = H . We perform a total-variation regularization of the differentiation as explained in appendix B. We find this approach advantageous to evaluating the fluctuation of energy as in the second expression of Eq. (16). We observe that evaluating the energy fluctuation requires larger MPS bond dimensions to achieve convergence. The magnetic susceptibility is given by where M = S z tot denotes the total magnetization, and H an applied magnetic field. In contrast to the specific heat, we find that the fluctuation in the second expression of Eq. (17) can be evaluated efficiently.
Results for the specific heat and magnetic susceptibility at half filling and p = 1/16 are shown in Fig. 6. At half filling, shown in Fig. 6(a) we find that the specific heat is well converged at D max = 2000. At low temperature, we observe a behavior C ∝ T 2 , which is expected for an antiferromagnetic insulator from spinwave theory [88][89][90]. The thermodynamics at half filling closely resembles the thermodynamics of the square lattice Heisenberg model [91][92][93][94]. Previous quantum Monte Carlo studies [91][92][93] have shown that the twodimensional antiferromagnetic square lattice exhibits a maximum in the specific heat at T /J ≈ 0.5. The energy scale exchange interaction in the Hubbard model is given by J = 4t 2 /U , which for t = 1 and U = 10 evaluates to J = 0.4. Hence, the observed maximum in the specific heat at T /t ≈ 0.2 = J/2 agrees well with the Heisenberg case. The magnetic susceptibility at half-filling shown in Fig. 6(c) is also converged at a bond dimension D max = 2000. It exhibits a maximum at T × /t ≈ 0.29, which indicates the onset of antiferromagnetic correlations. In the Heisenberg antiferromagnet, the magnetic susceptibility exhibits a maximum at T /J ≈ 1.0. We find, that in our case this maximum is shifted to slightly lower temperatures, as can analogously be observed in previous quantum Monte Carlo results of the Heisenberg model on finite width ladders [95].
Turning to the thermodynamics at hole doping p = 1/16 we also find the specific heat Fig. 6(b) to be well converged at D max = 2000. It exhibits a broad maximum around T /t ≈ 0.2, at a similar temperature as the half-filled case. But at temperatures around T /t ≈ 0.07 we also observe a small step-like feature. This temperature corresponds well with the onset temperature of stripe order from Fig. 3 as well as the spin gap shown in Fig. 4. Between the step-like feature and the maximum from T /t ≈ 0.08 to T /t ≈ 0.175 we observe a regime where the specific heat is approximately linear in temperature, C ∝ T . Also, below the step-like feature at T /t = 0.07 our data suggests a linear temperature regime, although this observation is only based on few data points. The magnetic susceptibility in Fig. 6(d) is well converged at D max = 3000, where we observe a slight uptick at D max = 2000 at lower temperatures. The magnetic susceptibility attains a maximum at T * /t ≈ 0.25. This compares well to previous results from the finitetemperature Lanczos method [96] on small lattice sizes, that also detected a maximum in the magnetic susceptibility at a comparable temperature and doping. It is also in agreement with calculations of the uniform susceptibility with cluster extensions of DMFT [51,97]. In underdoped cuprate superconductors, the pseudogap was indeed first identified experimentally as a suppression of the magnetic susceptibility below a temperature T * larger than the superconducting T c [98,99].

VII. METTS SIMULATIONS OF THE TWO-DIMENSIONAL HUBBARD MODEL
In this second part, we demonstrate that the METTS method can indeed be successfully performed for the Hubbard model approaching two-dimensional geometries. We show that several quantities of interest, like specific heat, magnetic susceptibility, structure factors, and momentum distribution functions can be reliably computed with reasonable computational effort over a wide range of temperatures. We discover several interesting facts about the METTS algorithm. A key practical question is how many METTS have to be random sampled. This depends on the variance of the estimator, where a small variance implies that less METTS are needed to achieve a certain statistical error. Interestingly, we find for several quantities, that both decreasing temperature as well as increasing the system size decreases the variance significantly. The METTS algorithm involves computing an imaginary-time evolution of product states. Modern MPS algorithms allow for performing this time evolution with high accuracy. Our mapping of the linear MPS chain onto the two-dimensional square cylinder geometry is shown in Fig. 8. We employ a combination of the time-dependent variational principle (TDVP) [100,101], and the time-evolving block decimation (TEBD) [102,103] algorithms and demonstrate the accuracy of this approach by comparing to numerically exact Lanczos time evolution on a smaller system size. The algorithms we choose, come with several control parameters. We find some of them can be chosen highly accurate without impacting performance or, otherwise, allow for an optimal choice. We identify the maximal bond dimension for performing the TDVP algorithm to be the main control parameter for the accuracy of the imaginary time evolution. We assess the accuracy by comparing results to two other state-of-the-art methods. Firstly, we compare to the method of thermal pure quantum (TPQ) states [104][105][106], which allows simulating smaller systems in a statistically exact way, also at finite-doping. Secondly, we compare to auxiliary-field quantum Monte Carlo (AFQMC). At half-filling, this method is also statistically exact and allows for simulating larger system sizes. Away from halffilling the method cannot be applied without encountering a sign problem, although they can be performed using the constrained-path approximation [107]. At finitedoping we study our results as a function of the maximal bond dimension and find that several quantities can be converged. The bond dimensions required to do so, are significantly smaller than the bond dimensions reported to converge ground state DMRG calculations.
We focus on the technical and algorithmic aspects of these simulations. The METTS algorithm [71,72,79] combines the advantages of Monte Carlo simulation and tensor network algorithms to simulate finite-temperature quantum many-body systems.

A. Basic METTS algorithm
We briefly review the basic steps of the METTS algorithm. For more details we refer to Refs. [71,72,79]. Given an orthonormal basis {|σ i } of the Hilbert space, the thermal average in Eq. (2) can be written as, where we introduce, The so-called typical thermal states |ψ i are normalized ( ψ i |ψ i = 1) and the weight p i defines a probability dis-tribution, i.e.
The thermal average as defined in Eq. (19) is, thus, amenable to Monte Carlo sampling. Notice, that all weights p i are manifestly real and non-negative. Therefore one does not encounter a sign problem when using such an ensemble. The tradeoff is that the states |ψ i are entangled and one must find a way to represent and manipulate them efficiently.
To construct a Markov chain with stationary distribution p i , Refs. [71,72] introduced the transition probability, which fulfills the detailed balance equations, According to Markov chain theory, the average of the sequence of measurements, converges to the thermal expectation value, i.e.
almost surely, provided ergodicity. R denotes the length of the measurement sequence, i.e., the number of METTS computed.
The algorithm as described above does not make any reference to tensor networks yet. However, its individual steps can be efficiently performed within the framework of matrix product states (MPS), if the basis states |σ i are only weakly entangled. We take the |σ i to be product states of the form which are states with zero entanglement. With this choice, the typical thermal states |ψ i are in a certain sense minimally entangled, hence referred to as METTS.
We illustrate the METTS algorithm in Fig. 7. It consists of the following basic steps: (i) Initial state: Choose a suitable first state of the Markov chain |σ 1 , set i ← 1.
(ii) Time evolution: evolve the state |σ i in imaginary time by an amount β/2 and normalize to compute the state (iv) Collapse: choose a new basis state |σ i+1 according to the probability distribution | ψ i |σ i+1 | 2 . Set i ← i + 1 and iterate from step (ii).
The sequence of measurements { ψ i |O|ψ i } is then analysed using standard time-series analysis techniques to compute an (error) estimate for the thermal average

B. Initial State
The choice of the initial product state in the METTS algorithm is rather important and can result in numerical difficulties if improperly chosen. We begin with a random product state. It is chosen according to a uniform distribution on the space of product configurations with a given particle number. Such a state is in general not related to the low-energy physics of the system. Directly starting the METTS procedure from this state can result in long thermalization times, especially at higher temperatures. Also, such unphysical random states can result in large bond dimensions if time evolved with given accuracy and can, therefore, complicate computations. To circumvent this problem, we perform several initial DMRG sweeps on the uniform random initial state to obtain a more physical state. This state is then collapsed to a new product state |σ 1 , which will then be taken as the initial state of the METTS algorithm. These initial DMRG sweeps are not performed until convergence. We typically choose to do 5 sweeps at maximum bond dimension D = 100. We also add a noise term [108] of 10 −4 , which further randomizes the starting state.

C. Imaginary-time evolution
The computation of the imaginary-time evolution, poses the key algorithmic challenge in the METTS algorithm. Time evolution algorithms for matrix product states are subject of current research and several accurate methods have been proposed [100][101][102][103][109][110][111][112]. The recent review article by Paeckel et al. [113] summarizes and compares a variety of these methods. Among those, the time dependent variational principle (TDVP) [100,101] method has been shown to have favorable properties in various scenarios, including imaginary-time evolution. The TDVP algorithm is closely related to DMRG [101]. Instead of solving an effective local eigenvalue problem, a time evolution of the effective Hamiltonian is performed, followed by a backward time evolution on one less site. We refer the reader to Refs. [101,113] where N denotes the number of sites, D the MPS bond dimension, and d the site-local dimension. In the case of the Hubbard model, d = 4. We, therefore, expect the single-site TDVP algorithm to be approximately four times faster than the two-site variant, given a fixed maximal bond dimension D max . We choose a time evolution strategy, where we first increase the MPS bond dimension using two-site TDVP up to a maximal bond dimension D max . Afterward, we switch to a single-site TDVP algorithm at fixed bond dimension D max for speed. Our implementation is based on the ITensor library and is available online [114]. There is, however, a caveat to directly applying the two-site TDVP algorithm in the context of METTS. TDVP suffers from a projection error onto the manifold of MPS at given bond dimension [100,101,113]. This error becomes non-negligible for MPS of small bond dimension. In particular, time evolution of product states as in Eq. (26) will suffer from a substantial projection error.
We circumvent this problem by starting the time evolution with a different method, the time-evolving block decimation (TEBD) [102,103] algorithm. This method has been widely used in a variety of numerical studies, including applications in the context of METTS [71,72,79]. The TEBD algorithm we apply uses a second-order Suzuki-Trotter decomposition, e −τ H = e −τ /2h1 e −τ /2h2 · · · e −τ /2h2 e −τ /2h1 +O(τ 3 ), (31) where H = ordering, swap gates are applied [72]. In the twodimensional cylinder geometry shown in Fig. 8 this is done for all hopping terms in the long direction and hopping terms "wrapping" the cylinder in the short direction. For details on implementing swap gates in for twodimensional systems, we refer the reader to Ref. [72]. We summarize our time evolution strategy in Fig. 9. Initially, we apply the TEBD with high accuracy up to an imaginary time τ TEBD to obtain an MPS with a suitably large bond dimension. Thereafter, we employ the two-site TDVP algorithm to further increase the bond dimension as we go lower in temperature. The TDVP time evolution can decrease the bond dimension obtained after TEBD at intermediate times since it is usually performed with a lower cutoff ε. Once we encounter maximum bond dimension D max , we switch to single-site TDVP for computational efficiency.
Several parameters control the accuracy of our time evolution strategy. We have performed an extensive investigation of their behavior, which we discuss in section VIII. Summarily we find, that most parameters can be kept at a fixed value. The initial TEBD parameters can be chosen highly accurately to not yield any substantial error. For the TDVP time step ∆τ there appears to be an optimal choice. We find that a choice of ∆τ = 0.02, ε = 10 −12 , τ TEBD = 0.1, (TEBD), (32) yields a negligible time-evolution, as well as TDVP projection error. The two-site TDVP algorithm comes with two control parameters. The time step size ∆τ and the SVD cutoff ε. We analyze the accuracy of the time evolution when varying these two parameters over several orders of magnitude in Fig. 10. We find that the cutoff parameter ε is directly related to the accuracy. Remarkably, we find that for most choices of the cutoff ε, a time step size of ∆τ = 0.5 (TDVP), yields optimal accuracy and is also favorable in terms of computational efficiency. Different studies have analogously reported that TDVP allows for rather large time steps [113,115]. The final control parameter is given by the maximal bond dimension D max used in the single-site TDVP.

D. Collapse
After performing measurements on the state |ψ i , we choose a new product state |σ i+1 in step (iv) of the METTS algorithm. |σ i+1 is chosen according to the probability distribution | ψ i |σ i+1 | 2 . The algorithm we use for sampling a random product state from an MPS is described in detail in Ref. [72]. The computational cost of the collapse step is where N denotes the number of sites, D the MPS bond dimension, and d the site-local dimension. This is a subleading computational cost as compared to the time-evolution, whose computational cost scales as O(N D 3 βd).
There is a freedom of choice of the product state basis. Here, we consider two bases. The local S z -basis is given by the states, A collapse into this basis will conserve particle number and total magnetization. Since also the imaginary-time evolution conserves all quantum numbers, the METTS will always stay in the same particle number and magnetization sector if the METTS |ψ i is collapsed into this basis. For simulating the canonical ensemble we have to allow for fluctuations in the magnetization. This can be achieved by projecting into the local S x -basis, where, Since S x and S z do not commute, the magnetization in S x direction can fluctuate for a state with fixed S z magnetization. Since the full Hubbard Hamiltonian is SU(2) invariant, the state projected into the S x can be rotated into a state in the S z basis. Therefore, we can reinterpret the S x basis as the S z basis with a fixed total magnetization. This allows us to employ S z conservation in the time-evolution algorithm again. Apart from allowing fluctuations in magnetization, the S x updates yield favorable mixing properties for the Markov chain. We discuss this in detail in appendix A. We find that the collapse into the S x -basis not only allows for fluctuation of the magnetization but also reduces autocorrelation effects considerably.

VIII. TIME EVOLUTION ACCURACY
Several parameters set the accuracy of the imaginarytime evolution method we describe in section VII C. We denote the METTS state computed using the MPS techniques by |ψ MPS . In order to assess their accuracy we compute overlaps with METTS states |ψ ED , which have been computed using Exact Diagonalization (ED). We choose a cylindrical system with L = 3 and W = 4. Although this lattice is rather small, it already introduces longer-range interactions, that are present in W = 4 cylinders and, thus, already poses a non-trivial problem for MPS time-evolution techniques. This system is also already too large to easily perform full ED. Therefore, we use a Lanczos method [116,117] to compute the exact reference state |ψ ED . We apply the algorithm and convergence criterion suggested in Ref. [118]. We compute the state |ψ ED up to a precision of 10 −12 , in the sense that, Hence, the state |ψ ED can be considered as quasi-exact. We choose the defect ∆, as a figure of merit to assess the accuracy of the MPS time evolution. For a given MPS state |ψ MPS we compute this overlap exactly, by computing the coefficients in the computational basis of the ED code. We study two different initial product states. We considered the Néel antiferromagnetic state at half-filling, p = 0, and an antiferromagnetic state, with two holes in the center of the system, p = 1/6, We focus on METTS |ψ MPS at temperatures T /t = 0.50, 0.10, 0.02 and choose U/t = 10. In order to avoid the projection error of TDVP when directly applied to product states, we start with a precise TEBD time evolution up to time β = 0.1. For doing so, we choose a time-step τ TEBD = 0.02 and a cutoff ε TEBD = 10 −12 . These parameters are chosen, such that the time evolved state at β = 0.1 is highly accurate, albeit with a potentially large bond dimension. The large bond dimension, however, is beneficial for decreasing the TDVP projection error [113].
The total defect ∆ in Eq. (39) comprises the error from the initial TEBD evolution, ∆ TEBD , the TDVP projection error, ∆ Proj , and the error of the bulk TDVP evolution, ∆ TDVP . Hence, As shown in Fig. 10, we achieve total defects smaller than ∆ < 10 −8 . The total defect is directly related to the TDVP parameters τ and ε. Therefore, we conclude that both the initial TEBD evolution error, ∆ TEBD , and the TDVP projection error, ∆ Proj are negligible as compared to the bulk TDVP evolution error.
We now focus on the two remaining control parameters of the two-site TDVP algorithm. We consider the step-size τ for one TDVP sweep and the cutoff parameter ε, which is used both as the magnitude of the discarded weight and the accuracy of the local effective timeevolution, cf. section VII C. Results for the obtained defect ∆ over a broad range of step-sizes and cutoffs are presented in Fig. 10. We observe that for all temperatures and hole-dopings choosing small time steps does not improve the error. In fact, in several cases, a time step of ∆τ = 0.01 yields the largest defect. This behavior is explained by the fact, that smaller time steps require more time steps to be performed. Since at every time step, an additional truncation of the MPS is performed, more time steps accumulate more truncation errors. Interestingly, we find that a choice of ∆τ = 0. yields optimal results for most temperatures, cutoffs, and dopings. Larger time steps then again decrease accuracy. Remarkably, we do not observe that the accuracy deteriorates strongly for longer time evolutions of the states we have chosen. The defect is directly related to the cutoff ε. In several cases for ∆τ = 0.2 and ∆τ = 0.5 we find, that the defect is of the same order of magnitude as ε.
These observations let us conclude, that a ∆τ = 0.5 for TDVP is optimal in the present context and that the accuracy of the time-evolved state can be precisely controlled by the cutoff. When using single-site TDVP in the final step of our time evolution strategy, the truncation error is controlled by a maximal bond dimension D max . A priori, we cannot predict the required value of D max to achieve converged results. Therefore, we always compare results from different values of D max to show convergence.

IX. METTS ENTANGLEMENT
The accuracy of MPS techniques is determined by the entanglement of the wave functions which are represented as an MPS. In our case, the METTS states |ψ i are MPS and it is thus interesting to investigate their entanglement. We use the Von Neumann entanglement entropy, as an entanglement measure, where ρ A = Tr B |ψ AB ψ AB | denotes the reduced density matrix in the subregion A of the lattice. At infinite temperature, the METTS are product states and, hence, their entanglement entropy vanishes. At lowest temperature, the METTS states approach the ground state, which implies that the METTS entanglement becomes comparable to the ground state entanglement. At intermediate temperatures, the behavior of the METTS entanglement is not a priori clear. The entanglement entropy upon imaginary-time evolving five different initial product states for p = 1/16 of the 32 × 4 cylinder at U/t = 10 is shown in Fig. 12(a). We observe a smooth increase of the entanglement entropy when lowering the temperature. Interestingly, for several initial product states the entanglement entropy decreases below a temperature of T /t = 0.05 again, while it remains growing for others. This is likely the result of the product states defining the METTS having varying overlap with the ground state. Some product states, with lower overlap, could be sampling the slightly larger entanglement of some low lying excited states. As one goes to zero temperature, this overlap factor become unimportant.
We show the average entanglement entropy of the METTS used in our simulations of the 32 × 4 cylinder at U/t = 10 in Fig. 12(b). The subregion A is chosen to be left half of the cylinder. As expected, we find the half-filled case, p = 0, to be less entangled at any temperature. For both values of doping we find the entanglement entropy increasing with decreasing temperature. The entanglement entropy has been computed from runs with different maximal bond dimensions D max = 2000, 3000, 4000. We find across the full range of temperatures that the entanglement entropy is converged T /t = 0.05, which is the approximate temperature where the transition from the antiferromagnetic regime to the stripe regime at low temperatures takes place.
While the maximal value of S vN 2.5 at p = 1/16 constitutes considerable entanglement, it is very well within reach of current MPS techniques. The behaviour in Fig. 12 also demonstrates that the METTS states |ψ i are very different from Hamiltonian eigenstates at higher temperatures, for which we would expect an increase in entanglement when increasing energy.

X. VALIDATION WITH TPQ AND AFQMC
To assess the validity of the METTS calculations across a broad range of temperatures and different dopings, we compare to current state-of-the-art methods. First, we focus on computing the thermodynamic energy E, specific heat C, and magnetic susceptibility χ m . The method of thermal pure quantum (TPQ) states [104,105] has been proven effective to extend the range of system sizes accessible via Exact Diagonalization techniques [106,[119][120][121]. It is also closely related to the finite-temperature Lanczos method [122].
The main idea of the TPQ method is that the trace of any operator O can be evaluated by, Here, D denotes the dimension of the Hilbert space and · · · denotes averaging over normalized random vectors |r . The coefficients of |r are independent and normally distributed. This is used to evaluate thermal averages as, where the TPQ state |β is given by, We notice that this state closely resembles the METTS |ψ i in Eq. (20). The random states |r , however, are not product states and are highly entangled in general. The TPQ state |β can be evaluated using Lanczos techniques [116,117]. Interestingly, the statistical error when random sampling over a finite number R of states is related to the free-energy density at a given temperature [123,124] and can be shown to become exponentially small when increasing the system size [104,105]. We refer the reader to Ref. [106,125] for an in-depth explanation of the method. Apart from a statistical error, the TPQ method has no systematic error and does not apply approximations. It is, however, limited to system sizes for which the Lanczos algorithm can be currently applied. Here, we compare data from METTS and TPQ obtained on a W = 4 and L = 4 cylinder. This system size is beyond the reach of numerical full diagonalization of the Hamiltonian matrix. The choice of a W = 4 cylinder poses a challenge to the METTS calculations due to the long-range nature of the interactions. We think that this benchmark addresses the key difficulties to be overcome in order to simulate longer length cylinders. Also, the case of p = 1/8 holedoping is afflicted by a sign-problem when investigated using conventional QMC methods. We show results in Fig. 11. The METTS calculations have been performed with a cutoff ε = 10 −6 and maximum bond dimension D max = 2000. We find agreement within errorbars between TPQ and METTS. For TPQ we have used R = 200 random vectors.
When comparing to exact results from TPQ, we are limited in system size. At half-filling, however, the auxiliary-field quantum Monte Carlo (AFQMC) method [126][127][128][129][130][131] can be applied without encountering a sign problem that allows for investigating larger lattice sizes. AFQMC at half-filling also does perform any approximation. We compare spin-correlation functions from METTS and AFQMC at various temperatures on t 32 × 4 cylinder in Fig. 13. Again, we employed a cutoff ε = 10 −6 and a maximum bond dimension D max = 2000 to perform the time-evolution in METTS.We find agreement within errorbars.
The comparisons in Figs. 11 and 13 show, that METTS agrees with current state-of-the-art unbiased numerical methods whenever they are applicable. We also find, that a maximum bond dimension D max = 2000 and a cutoff ε = 10 −6 in the TDVP time evolution are sufficient to obtain consistent results. The comparison between METTS and TPQ shows, that METTS is reliable at finite doping. The comparison to AFQMC, on the other hand, shows that our implementation of METTS yields consistent results when considering larger lattices. pling regime have been reported for the single-band [38] and three-band [37] Fig. 14. Indeed, we confirm the presence of incommensurate spin correlations consistent with the DQMC. The sign changes which are observed here are in the tails of the rapidly decaying correlation function, and thus they required substantial statistical averaging of the METTS to observe-here, we used 15000 samples. We find that the width of the antiferromagnetic "domains" (in the correlation function) is 5, which is exactly what has been found in [38]. Moreover, we show the density profile along the cylinder in Fig. 14(d). Here, we observe the same boundary density fluctuations as reported in Fig. 4 of Ref. [38]. In agreement with their findings, the density profile does not exhibit a charge-density wave pattern. We would like to point out, that Refs. [37,38] referred to the incommensurate spin correlations as "stripes", while our definition of a stripe refers to intertwined spin and charge ordering. Incommensurate charge correlations are not observed for  Fig. 4 of Ref. [38].
this particular set of parameters.

XII. DISCUSSION
The results from the METTS simulations presented in the previous sections yield several interesting new insights into the physics of Hubbard model in the strong coupling regime at U/t = 10 at finite temperature.
For the hole-doped case at p = 1/16 we find three different regimes as a function of temperature. At temperatures T /t 0.05 we find a stripe ordered phase. A typical METTS state in this phase is shown in Fig. 1(a). We observe antiferromagnetic domains whose domain walls coincide with peaks in the wave-like hole density modulation. In the present situation the stripes are halffilled. On a width-4 cylinder there are two holes per stripe wavelength. This leads to a magnetic ordering vector of k = (7π/8, π) = (π − π/8, π) and a charge ordering vector of k = (π/4, 0) (hence δk c = 2δk s ). This observation agrees well with recent DMRG results on the width-4 cylinders [30]. Even though the exact value of U/t = 10 and p = 1/16 is not included in the results of these authors, their phase diagrams suggest that this point realizes what is referred to as the Luther-Emery 1 (LE1) phase at T = 0. The LE1 phase is shown to exhibit half-filled stripes, in agreement with our findings. The phase diagrams in Ref. [30] have been established using DMRG with bond dimensions D = 5000, which is also the maximal bond dimension we have used in our DMRG calculations. Both the magnetic and the charge structure factors shown in Figs. 2 and 3 clearly indicate the occurrence of stripe order below T /t 0.05. We would like to mention that a similar crossover from an incommensurate to an antiferromagnetic regime has been found to occur in the three dimensional Hubbard model [132].
At temperatures above the stripe order, we find enhanced antiferromagnetic correlations. As shown in Fig. 4, the spin gap is estimated as ∆ s /t = 0.07 which approximately coincides with the onset temperature of the stripe order. Around this temperature the peak in the magnetic structure factor shifts from k = (7π/8, π) to k = (π, π). The transition or crossover between the two phases is further signalled by a step-like feature of the specific heat in Fig. 6(b). Hence, we identified several quantities that indicate a transition or crossover between the stripe ordered phase and a phase with enhanced antiferromagnetic correlations but no charge ordering.
The system at p = 1/16 differs significantly from the antiferromagnetic insulator at half-filling. The density structure factor at zero temperature in Fig. 2(c) indicates a small or vanishing charge gap by exhibiting approximately linear behavior at k = (0, 0) [33,[82][83][84][85]. We distinguish between the single-particle gap ∆ (1) c and the charge gap ∆ (2) c . While the single-particle gap is shown in Fig. 4 to attain a finite value of ∆ (1) c /t ≈ 0.25 for infinite cylinder length, the charge gap is shown to vanish as ∆ (2) c /t ∝ 1/L. Due to the finite single-particle gap, the system does not exhibit a Fermi surface. This is also evident in the momentum distribution function n σ (k) shown in Fig. 5(a), whose slope at T = 0 remains finite for p = 1/16, as expected for a system without a Fermi surface [86]. We find that for momentum k y = π/2 the slope of n σ (k) becomes largest close to the nodal point k = (π/2, π/2), indicating that the single-particle gap is smallest in this region. To further corroborate this finding we found that the real space electronic correlations Fourier transformed along the y-direction, F y (x l , x m , k y ) in Fig. 5(b), exhibit fast exponential decay at k y = 0, π. At k y = π/2 a slower exponential decay is observed. Therefore, we conclude that while the single-particle gap is still sizeable, the smallest gap to single particle excitations is realized in the nodal region around k = (π/2, π/2).
The novel short-range antiferromagnetic phase shares many features with the experimentally observed pseudogap region, which is characterized by a partial suppression of low-energy excitations [6,7,10]. In particular, angle-resolved photo emission spectroscopy measurements [133] typically detect a suppression of singleparticle excitations close to the antinode, k = (π, 0). Similarly, we find in Fig. 5 that the single-particle gap in our case is largest for y-momenta k y = 0 and k y = π, while it is smallest at k y = π/2.
The onset temperature T * of the pseudogap phase was originally identified experimentally as a decrease of the uniform magnetic susceptibility upon cooling below T * [98,99]. In our simulations at U/t = 10 and p = 1/16 we indeed detect a maximum in the magnetic susceptibility at T * /t ≈ 0.25 in Fig. 6. Below this temperature the onset of antiferromagnetic correlations with (π, π) wave-vector in Fig. 3 agrees well with previous results from cluster extensions of DMFT [62-64, 66, 97] and diagrammatic Monte Carlo [65], which indicate emerging short-range antiferromagnetic correlations in the pseudogap phase.
However, these approaches also point at metallic behaviour in the pseudogap regime, with gapless quasiparticles in the nodal region of the Brillouin zone. Gapless quasiparticles at the node are also expected for superconducting order with d x 2 −y 2 symmetry. In contrast, we find a non-vanishing single-particle gap which, although smallest close to the node, has a rather large magnitude ∆ (1) c /t ≈ 0.25. This is also the temperature scale below which we observe the onset of antiferromagnetism. Hence, the short-range antiferromagnetic correlations in the system that we simulate are not associated with a metallic state hosting gapless single-particle excitations. However, we would like to emphasize that our method is limited by only being able to simulate small width cylinders. It would be very interesting to learn how the single-particle gap evolves with cylinder circumference. Currently, ground state DMRG studies of the Hubbard model are possible on width 6 cylinders [19,28], so some additional information on the evolution of the gap with width should be available soon.
An investigation of pairing (correlations) in both the stripe phase and the pseudogap regime would give additional insights. Since this study requires further technical improvements of our method, this question will be addressed in future work. While the single-particle gap we detect might vanish in the full two-dimensional limit, there are also scenarios where the single-particle gap remains finite while the charge gap vanishes. This behavior is expected for superconductors with a nodeless gap function. Another intriguing possibility is an exotic metallic state, where the electrons fractionalize into chargons and spinons [134,135]. In this case, the single-particle gap would correspond to the spinon gap while chargons are still gapless.
We have also investigated the half-filled case at U/t = 10. We find that the magnetic structure factor converges smoothly towards the DMRG result when lowering the temperature in Fig. 2(a). It exhibits a clear peak at k = (π, π) indicating antiferromagnetism. We use the maximum of the static magnetic susceptibility in Fig. 6 at T × /t ≈ 0.29 to define an onset temperature of antiferromagnetic correlations. The quadratic behavior we observe for the specific heat at low temperatures is in agreement with thermodynamics derived from antiferromagnetic spin-wave theory [88][89][90]. We have also found that the thermodynamics at half-filling matches closely with quantum Monte Carlo results of the square lattice Heisenberg antiferromagnet [91,93]. In particular, the temperature of the maximum in the specific heat in Fig. 6 agrees well with those QMC results [91,93].
To the best of our knowledge, this work constitutes the first application of the METTS algorithm to study the Fermi-Hubbard model on geometries approaching the two-dimensional limit. Therefore, we have investigated the behavior of the algorithm in proper detail. We presented the statistical properties of the measurement time series in appendix A. As METTS is a Markov chain Monte Carlo algorithm, it is crucial to understand autocorrelation and thermalization properties. As we have shown in Fig. 15 the time series can indeed take rather long times to thermalize. At half-filling we found that there can be metastable states with multiple antiferromagnetic domains that are realized before the system becomes thermalized. Interestingly, we find that upon holedoping the system, the measurement time series equilibrate faster. We also investigated autocorrelation effects. We find that applying the S x instead of the S z updates described in section VII D reduces the autocorrelation times significantly, as shown in Fig. 16. For the temperatures and observables O investigated in this manuscript using the S x update, we find autocorraletion time τ [O] ≈ 1. We find that for T /t ≤ 0.5 for U/t = 10 and a 32 × 4 cylinder autocorrelation effects are negligible for most observables. A noticeable exception is the charge density structure factor at temperatures above T /t 0.20 shown in Fig. 3(b) where moderate autocorrelation effects have to be taken into account to compute proper error estimates.
The statistical error estimate, apart from the number of samples and the autocorrelation time, also depends on the variance of the time series. We have investigated the behavior of the variance of several observables as a function of temperature in Figs. 17 and 18 and found rather remarkable behavior. The variance of several estimators decreases rapidly when lowering temperatures. Also, increasing the system size decreases the variance of the energy density estimators. This observation is likely related to the phenomenon called quantum typicality [124]. For the related thermal pure quantum states (TPQ), which we have used to validate the METTS results in section X, a theory of their statistical properties has been developed in Refs. [104,105]. The TPQ and METTS states are obtained by imaginary-time evolving states drawn from some class of random initial states. Whereas for a METTS the random initial states are product states, the initial state of a TPQ state is state with (normally distributed) random coefficients in an arbitrary orthonormal basis of the Hilbert space. For the latter, it has been found that the statistical error when averaging over several TPQ states becomes exponentially small in the system size, under mild assump-tions on the system and observables. It would be interesting to arrive at a similar understanding of the variance of METTS states to explain our observations. Several methods have been proposed to further improve upon the variance of the estimators [136][137][138]. Refs. [136,137] proposed a hybrid approach interpolating between purification and METTS sampling. While using additional auxiliary sites for purification are likely to increase the necessary bond dimensions of the MPS to achieve convergence, they are favorable in the sense that fewer random states have to be sampled. Especially at higher temperatures, these approaches might yield a significant computational speedup. Also, approaches to incorporate additional symmetries in the METTS algorithm have been proposed [139,140]. We would also like to mention that several other tensor network approaches to finite temperature simulations have recently been applied successfully to a variety of problems [141][142][143][144][145].
Apart from the statistical analysis, the METTS method relies on accurate imaginary-time evolution algorithms. We find that the TDVP algorithm for timeevolving matrix product states is an appropriate choice. However, the projection error when applying TDVP directly to MPS of small bond dimension, especially product states, is not negligible. We circumvent this problem by applying an initial TEBD time-evolution up to a certain imaginary-time, thus increasing the bond dimension before applying the TDVP algorithm. We find that this approach yields very accurate results compared to numerically exact Lanczos time evolutions. Also, when fixing a maximal bond dimension, we apply single-site TDVP, which has favorable computational costs. We would like to point out that recently a different approach to solving the projection error problem of TDVP has been proposed [112], which employs a subspace expansion using a global Krylov basis.
Simulating the doped Fermi-Hubbard model at finitetemperature poses a challenging problem for numerical methods. To demonstrate the accuracy of our simulations we have performed comparisons to four different well-established exact numerical methods in their respective limits. First, the limit T /t → 0 is compared to ground state DMRG calculations in Fig. 2. We find that the METTS simulations of the structure factors at T /t = 0.025 very closely resemble the ground state result from DMRG. This also shows, that for the investigated system our METTS simulations can essentially cover the temperature range down to temperatures which can be considered to realize ground state physics. Second, we have compared thermodynamics at half-filling and finite hole-doping on a 4 × 4 cylinder to the TPQ method [104][105][106] over a large range of temperatures. The TPQ method is an extension of exact diagonalization, where traces over statistical density matrices are replaced by random averages of random vectors. The method is considered to be statistically exact. Also here, we find perfect agreement to METTS shown in Fig. 11. Together with the comparison to DMRG this shows, that METTS yields accurate results at finite hole-doping over the full range of temperatures investigate. Finally, we also demonstrated in Fig. 13 that at half-filling METTS simulations of spin-correlations agree within error bars to statistically exact AFQMC simulations on larger system sizes. Finally, we confirm previous results [38] by determinantal quantum Monte Carlo obtained in the challenging parameter set U/t = 6, p = 1/8, t /t = −0.25, and T /t = 0.22 on a 16 × 4 cylinder. We correctly reproduce the incommensurate spin correlations and boundary density fluctuations found earlier.

XIII. CONCLUSION
Using the numerical METTS method, we have investigated finite temperature properties of the doped Hubbard model in the strong coupling regime on a width-4 cylinder. We focused on hole-doping p = 1/16 and half-filling. In the doped case we find that the ground state half-filled stripe phase extends up to a temperature T /t 0.05. Above this temperature, a phase with strong antiferromagnetic correlations is realized. In this regime, the specific heat exhibits behavior linear in T . The onset temperature of the short-range stripe order coincides with the spin gap, which has been shown to attain a finite value. A closer inspection of electronic correlations revealed that no full electron-like Fermi surface is realized, however. Instead, we find a vanishingly small gap to paired charge excitations, while the single particle gap of size ∆ (1) c /t ≈ 0.26 is still sizeable. By investigating the momentum distribution function and electron correlation functions we establish that the single-particle gap is smallest in the nodal region close to k = (π/2, π/2). The magnetic susceptibility realizes a maximum at temperature T * = 0.25t. These features are strongly reminiscent of the pseudogap phase realized in the cuprates. These findings have been made possible by combining recent tensor network techniques to simulate finite-temperature quantum systems and perform imaginary-time evolution of matrix product states. We found that both increasing the system size as well as lowering the temperature significantly reduces the variance of the METTS estimators, hence reducing the need to average over many random samples. Apart from benchmarking the timeevolution with numerical exact diagonalization, we validated our simulations by comparing to four different numerical methods, DMRG, TPQ, AFQMC, and DQMC. These comparisons demonstrate that systems of stronglycorrelated fermions at low temperatures, finite-hole doping, and on large lattice sizes can now be reliably simulated using the METTS method.

ACKNOWLEDGMENTS
We would like to thank Sebastian Paeckel, Salvatore Manmana, Michel Ferrero, Tarun Grover, Thomas FIG. 15. Time series of energy and magnetic structure factor measurements from 5 different random intial states. Data shown for a 32 × 4 cylinder with parameters T /t = 0.300, U/t = 10. We used time evolution parameters Dmax = 3000 and ε = 10 −6 and S x collapses. After initial thermalization no autocorrelation effects are apparent. (a) energy measurements at half-filling, p = 0. We observe several initial plateaux in the energy before thermalization. These plateaux are metastable states that correspond to antiferromagnetic domains. (b) energy measurements at hole-doping p = 1/16. (c,d) measurements of the magnetic structure factor S(k) evaluated at ordering vector k = (π, π) at p = 0 and p = 1/16. We observe a skewed distribution, fast thermalization and no apparent autocorrelation effects.
Köhler, Olivier Parcollet, Thomas Schäfer, Subir Sachdev and André-Marie Tremblay for fruitful discussions. We also thank Matthew Fishman for support in using the ITensor library (C++ version 3.1) [81]. SRW was funded by the NSF through Grant DMR-1812558. The Flatiron Institute is a division of the Simons Foundation.

Appendix A: Time series analysis
In this section, we discuss the statistical properties of the METTS sampling. As described in section VII A, the METTS algorithm is an example of a Markov chain Monte Carlo simulation. As such, thermalization and autocorrelation properties of the resulting time series have to be investigated to derive (error) estimates of physical observables. We find that the behavior of the METTS algorithm as a function of temperature is quite opposite to usual quantum Monte Carlo simulations. Thermalization and autocorrelation times decrease when lowering the temperature. Eventually, using the S x -basis collapse instead of the S z -basis collapse, as discussed in section VII D, we find that autocorrelation effects can be neglected for the temperature range we study. We also show that the variance of several quantities decreases quickly for lower temperatures. This allows reaching a higher statistical precision at lower temperatures, where imaginary-time evolution is more computationally expensive.
For an observable O the METTS algorithm yields a time series of measurements, The thermal average O in Eq. (2) is then estimated by  Error bars indicate a 95% confidence interval. (a) Energy measurements. We observe a fast decrease of the variance as lowering temperature. (b) Magnetic structure factor S(k) at ordering vector k = (π, π). (c) Magnetic structure factor S(k) at k = (0, 0). (d) momentum distribution function n ↑ (k) at k = (π, 0). and the estimated autocorrelation function ρ[O](l) is given by [147], The autocorrelation function ρ[O](l) is typically exponentially decaying in l. The cutoff value M ∝ τ [O] is chosen such to assure convergence of the sum in Eq. (A5), but small enough not to include numerical noise. For a proper choice of M , see Ref. [146].
We show examples of the measurement time series for different quantities in Fig. 15. We choose a temperature of T /t = 0.300 and Hubbard interaction U/t = 10 on a 32 × 4 cylinder. The imaginary-time evolution has been performed with TDVP cutoff ε = 10 −6 and maximum bond dimension D max = 3000. We used the S x basis collapse scheme. For energy measurements in the half-filled case in Fig. 15(a) we see that several hundred steps can be required to thermalize the system. Before reaching equilibrium, the METTS algorithm is temporarily stuck in metastable states. By a closer inspection of these states, we find that they typically contain several antiferromagnetic domains instead of a uniform antiferromagnetic state, which is realized once the system is thermalized at this temperature T /t = 0.300. We consistently find this behavior for temperatures T /t ≥ 0.200. We observe that the time for thermalization decreases when lowering the temperature. This is explained by the fact that a longer imaginary-time evolution constitutes a more Similarly, the measurements of the magnetic structure factor S(k) evaluated at ordering vector k = (π, π) in Fig. 15(c) exhibit prolonged thermalization at half-filling, although less pronounced than for the energy measurements. Fig. 15(b) shows the energy measurement time series for a 32 × 4 cylinder with T /t = 0.300, U/t = 10, L = 32 at hole-doping p = 1/16. Interestingly, no initial metastable plateaux are observed and the system thermalizes more rapidly than in the half-filled case. This could be explained by the additional charge fluctuations "smoothing" the energy landscape around the metastable antiferromagnetic domain states. Fig. 15(d) shows measurements of S(k) at k = M for p = 1/16. Also for this observable, the system is quickly thermalized. The distributions of the magnetic structure factors in Fig. 15(c) and (d) are skewed.
After an initial period of thermalization, we do not observe apparent autocorrelation effects in all time series shown in Fig. 15. This is owed to the fact, that the S x -basis collapses reduce the autocorrelation time substantially as compared to S z -basis updates. In Fig. 16 we compare the two collapse strategies for T /t = 0.400, U/t = 10, L = 16, and W = 4 at half-filling. As shown in Fig. 16(a) and (b), the energy time series shows significant autocorrelation effects when applying the S zbasis collapse. This is characterized by the autocorrelation time, estimated here as τ [H] = 9.46. Using S xbasis collapses instead, we find that subsequent samples are essentially uncorrelated, as shown in Fig. 16(c) and (d). There, we estimate an autocorrelation time of τ [H] ≈ 1.08.
Apart from the autocorrelation time and the number of random samples R, the statistical error estimate in Eq. (A3) is determined by the variance of the measurements, Var [O]. While the number of random samples R can be adjusted to achieve a fixed statistical error σ[O], both the autocorrelation time and the variance are intrinsic properties of the observable and the METTS algorithm. We investigate the behavior of the variance for several quantities in Figs. 17 and 18. We estimated a 95% confidence interval for the variance estimator Eq. (A4) using bootstrap resampling [148]. Fig. 17 compares variances for different system sizes and temperatures at holedoping p = 1/8. Remarkably, we find that both decreasing temperature as well as increasing the system size decreases the variance of the energy density considerably, cf. Fig. 17(a). Low temperatures and large system sizes are of course the more challenging regimes to perform the MPS time evolution. Hence, in order to achieve comparable statistical error estimates, less MPS time evolutions have to be performed. Similarly, we observe in Fig. 17(b,c) that the variance of the magnetic structure factor evaluated at k = (0, 0) and k = (π, π) decreases as a function of temperature. Also for the momentum distribution function at k = (π, 0) in Fig. 17(d) we observe that increasing the system size decreases the variance.
We compare estimators of variances at different holedopings in Fig. 18. The energy and magnetic structure factor of k = (0, 0) variances in Fig. 18(a,c) do not depend strongly on the filling fraction. We observe a rapid decrease of the variance when lowering the temperature. In contrast, the variance of the magnetic structure factor evaluated at M = (π, π) in Fig. 18(b) does not rapidly decrease at half-filling. This can be attributed to the fact, that S(k) develops a peak at k = (π, π) at low temperatures indicating the development of strong antiferromagnetic correlations. The variance of the momentum n(k) distribution function evaluated at k = (π, 0) also shows interesting behavior in Fig. 18(d). The variance is largest for hole-doping p = 1/8, contrary to the magnetic structure factor. Appendix B: Numerical differentiation using total variation regularization In the main text, we compute the specific heat as the derivative of the total energy w.r.t. the temperature, The energy is computed using METTS sampling and is hence afflicted by a statistical error δ. When performing a numerical finite-difference with spacing h, the straightforward error estimate is proportional to δ/h and hence diverges as h → 0.
An established way to estimate derivatives of noisy data is to perform total variation regularization [149]. Crudely speaking, the total variation of a function is a measure of how rapidly it oscillates. More precisely, the total variation of a differentiable function of one variable u is given by, where the prime denotes its derivative. We would like to point out that this corresponds to the L 1 -norm of the derivative of u and not the more conventional L 2 -norm. Thermodynamic observables, like the specific heat, are not expected to rapidly oscillate. Therefore, it can be assumed that their total variation is of small magnitude.
We describe the method proposed by Ref. [149]. Given a function f , we seek to find a regularized derivative u α of f , minimizing the functional, Here, Au(x) = x 0 u denotes the operator of antidifferentiation (assuming f (0) = 0). The prefactor α is the regularization parameter. A minimizing function u α of F is thus in balance between closely integrating to f and having small total variation. Ref [149] proposes an efficient algorithm to solve this minimization problem, based on the lagged-diffusivity algorithm [150]. In our application, f=E and u α= C.
We are given data points {f (x i )} N i=1 with corresponding error estimates {δf (x i )} N i=1 . To derive an estimate for the regularized derivative {u α (x i )} N i=1 and its standard deviation {δu(x i )} N i=1 , we apply a simple bootstrap idea. We randomly sample datapoints according to normal distributions with mean f (x i ) and standard deviation δf (x i ). For every random sample we compute the regularized derivative by minimizing Eq. (B3). From those samples of regularized derivatives we then estimate the mean and standard deviation. We use R = 100 random samples to perform these estimates.