Absence of superconductivity in the pure two-dimensional Hubbard model

We study the superconducting pairing correlations in the ground state of the doped Hubbard model -- in its original form without hopping beyond nearest neighbor or other perturbing parameters -- in two dimensions at intermediate to strong coupling and near optimal doping. The nature of such correlations has been a central question ever since the discovery of cuprate high-temperature superconductors. Despite unprecedented effort and tremendous progress in understanding the properties of this fundamental model, a definitive answer to whether the ground state is superconducting in the parameter regime most relevant to cuprates has proved exceedingly difficult to establish. In this work, we employ two complementary, state-of-the-art many-body computational methods, auxiliary-field quantum Monte Carlo (AFQMC) and density matrix renormalization group (DMRG) methods, deploying the most recent algorithmic advances in each. Systematic and detailed comparisons between the two methods are performed. The DMRG is extremely reliable on small width cylinders, where we use it to validate the AFQMC. The AFQMC is then used to study wide systems as well as fully periodic systems, to establish that we have reached the thermodynamic limit. The ground state is found to be non-superconducting in the moderate to strong coupling regime in the vicinity of optimal hole doping.


I. INTRODUCTION
Understanding high-temperature superconductivity in the cuprates [1] has been a long-standing mystery and one of the greatest challenges in theoretical condensed matter physics [2]. Very early on the single-band twodimensional (2D) Hubbard model [3], along with its cousin, the t-J model, were argued to be the paradigmatic models for this problem [4,5], and in many ways, this suggestion has proven to be accurate. Many of the properties of the cuprates seem to be reasonably well described -or at least mirrored -in the Hubbard model, such as antiferromagnetism [6][7][8] and its abrupt disappearance upon doping, pairing and stripe formation, and pseudogap physics [9]. Pairing, when it occurs, can be seen as a consequence of a sort of frustration, between the hopping/kinetic energy of the holes and antiferromagnetic correlations, which are disrupted by the hopping.
Superconducting long-range order itself, however, is one of the most delicate properties in these systems. Superconductivity appears to have a subtle competition and * These two authors contributed equally to this work. coexistence with stripe formation [10][11][12]. In terms of the models, this means that accurate answers about a possible superconducting phase require simulations which are able to describe all of the possible phases in an unbiased fashion, so that their competition can be resolved. One also needs a systematic approach to the zero-temperature as well as the thermodynamic limit, particularly since stripes can introduce a new length scale somewhat larger than the size of a pair. Numerous studies over the years have addressed pairing order in the Hubbard model. They have often been driven by remarkable methodological advances, and have led to a great deal of insight in the physics of the model (see, for example, Refs. . However, given the competing energy scales and intertwined states, it can reasonably be argued that none has satisfied these rigorous criteria for establishing the nature of superconductivity in the physically relevant parameter regime. Both positive and negative results have been found for d-wave pairing order, reflecting the extreme sensitivity of the ground state and low-lying excitations in the model, and the competition between d-wave and other states [33,[39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54]. The competition between superconductivity and stripes [33,[39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54] or other orders is also strongly affected by modifications of the model, such as the next nearest neighbor hopping t . Given the existence of superconductivity in the cuprates with an apparent elec-tronic mechanism, it seems likely that some modification of the pure model exhibits superconductivity. For example, recent studies on width-four cylinders-where DMRG can be pushed to resolve the competing phases to high accuracy-found a non-superconducting filledstripe state in the pure model, but quasi-long-range pairing correlations with the addition of a t term, coexisting with half-filled stripes [38,[55][56][57]. While superconductivity arising from the addition of a t is encouraging, one clearly needs to go beyond width four. (Note that a width four cylinder is equivalent to a stack of plaquettes, and there is no difference between a pair on a plaquette and a half-filled stripe. Larger systems are needed to properly allow stripes and superconductivity to compete or coexist.) Hopping parameters t and third neighbor (diagonal) t have been predicted using electronic structure methods [58], but even small differences in these parameters can alter the ground state phase and it is difficult to establish whether additional terms, such as hopping mediated by a second hole, are important. It is also not clear whether one needs to study a three-band model in order to connect directly with the cuprates.
Here, we choose to focus on the pure Hubbard model, with parameters U and t only. The existence or absence of superconducting order in this fundamental model at moderate to strong coupling is an outstanding theoretical question. This question has presented a 30-year challenge, magnified by the quest to understand high-T c superconductivity. An intense experimental effort is on-going with ultracold atoms in optical lattices to realize "quantum simulations" of this model [59][60][61][62]. The model has also served as a barometer for the capacity of the computational physics and chemistry community to perform reliable computations in interacting quantum systems. We study pairing correlations and superconductivity using two complementary methods, the density matrix renormalization group (DMRG) and auxiliary field quantum Monte Carlo (AFQMC). Our work follows up on a previous study involving four different methods which determined that the ground state of the Hubbard model has stripe order at 1/8 doping [63]. Although stripes may tend to compete with superconductivity [64], it may be possible for them to coexist [12,34,[65][66][67][68].
The constrained path AFQMC method [69,70] we use treats the fermion sign problem approximately, so validation is important. Here we use DMRG [71] on width four and six cylinders to validate an approach to predict pairing orders in AFQMC. The DMRG calculations involve multiple independent DMRG programs pushed to the limit of current capabilities. We find excellent agreement between the DMRG and AFQMC. The AFQMC does not have DMRG's width restrictions, and we then use the AFQMC to study system as large as 16 × 16, including periodic boundary conditions. In the AFQMC calculations, we devise new techniques to probe the superconducting order, both through a linear response measure of the order parameter, and through the use of a BCS trial wave function to directly measure the pairing correlation function. These simulations allow us to conclude that only short-range pairing occurs in the regime of interest (U/t around 6-8 and dopings 0.1 < h < 0.2), and the system is not superconducting.
In the small U/t limit, controlled results from perturbation theory have shown that the Hubbard model has a superconducting ground state [72][73][74][75]. Diagrammatic Monte Carlo studies [76] indicate that a BCS superconducting state of d-wave symmetry can emerge at weak coupling (U/t < 4) for doping h ≥∼ 0.3. Given the sensitive and delicate nature of the ground state of the model, and, in particular, given that stripe formation is believed not to occur at weak coupling [75], it is a very interesting question how this part of the phase diagram connects with the other parameter regimes. We emphasize that our work does not imply a general statement that there is no superconducting order anywhere in the pure Hubbard model. Rather our focus is on the nature of the pairing order in the pure Hubbard model in the physically important parameter regime as a model for cuprate superconductors.
The rest of this paper is organized as follows. Section II discusses the two different methods we employ, and two different ways in each to probe pairing and superconducting order. Our results are presented in Sec. III: first a general scan of the doping dependence of the superconducting order at U = 8, then a detailed study of the case of U = 8 and h = 1/8, followed by an analysis of the relation between pairing and stripe order, and then the dependence on the interaction strength. We conclude in Sec. IV. Further technical details as well as additional results are included in the Appendix.

II. APPROACH
We study the pure Hubbard Hamiltonian with nearestneighbor hopping and on-site interaction: whereĉ iσ is the fermionic annihilation operator, σ denotes spin (=↑ or ↓),n iσ =ĉ † iσĉ iσ is the particle number operator on site i, and ij denotes the nearest neighboring sites. We study rectangular lattices of size N = L x × L y , typically with periodic boundary conditions (PBC) along the y direction and open boundary conditions along the x direction (i.e., cylinder geometry). We vary the aspect ratios of the cylinders (e.g., 32 × 8, 24 × 14) to ensure that the rectangular cells do not impact our results [77]. Finite-size extrapolations are performed. Additionally, complementary calculations of the pairing correlation functions are performed with PBCs along both directions. We set t as the energy unit, i.e., t = 1.
We denote the number of electrons in the simulation cell by N e , with N e = N ↑ + N ↓ . The electron density or filling factor is n = N e /N , and the hole doping level is then h = 1 − n. These quantities are specified in an average sense, as N e is controlled by the chemical potential µ and will fluctuate in most of our calculations.

A. Two complementary methods
In this work we employ two state-of-the-art methods, constrained path (CP) AFQMC and DMRG. These methods are representative of the leading edge of computational capabilities for interacting quantum manyfermion systems. They involve very different approximations in obtaining ground-state properties in the thermodynamic limit. To quantify the CP error in AFQMC, we benchmark the results in finite systems of narrow cylinders, where DMRG is highly accurate. The AFQMC does not have size or boundary condition restrictions and can reliably approach the thermodynamic limit. The systematic, detailed, and complementary use of these leading computational techniques is a unique and distinguishing feature of the present study. The excellent agreement between the two methods allow us to draw conclusions with confidence.

Constrained-path auxiliary field quantum Monte Carlo
In AFQMC, the interaction part of the Hamiltonian is re-cast into a summation (or an integral) of noninteracting terms through a Hubbard-Stratonovich transformation. As a result, physical quantities are represented as a path integral in many-dimensional auxiliary field space. The high-dimensional summation or integral can be evaluated with Monte Carlo techniques [6]. However, with few exceptions, a minus sign problem is present [78] which causes an exponential growth of the statistical errors with system size. The CP approach overcomes this difficulty by imposing a boundary condition in auxiliary field space, which is derived from an exact property of the path integral [69] but whose practical implementation involves a trial wave function. The use of CP introduces a systematic error, which can be improved with better trial wave functions. Usually simple trial wave functions such as the Hartree-Fock solution have been used as trial wave-functions and previous results [79] show the systematic error is typically small. Recently we introduced an approach [70] to optimize the trial wave-function selfconsistently, further reducing the systematic error. As mentioned, a key feature of this work is the combined use of CP-AFQMC with DMRG, which allows us to systematically gauge the accuracy of CP in cylindrical systems.
Ground-state AFQMC is typically formulated in a sector of the Hilbert space with fixed number of particles, N e , and fixed S z (although a corresponding approach in Hartree-Fock-Bogoliubov space exists [80]). Our computation of the pair-pair correlation function is done in this manner, by separate AFQMC calculations on the original Hubbard Hamiltonian in Eq. (1), using back-propagation [69] and BCS trial wave functions [81]. In this work the order parameter is computed in AFQMC after a particlehole transformation has been applied to Eq. (1), which results in a modified Hamiltonian that conserves the total particle number [36] but breaks total S z (further details in Appendix). As described in the next section, in this formulation the order parameter can be computed from total energy calculations, which leads to very accurate results.

Density matrix renormalization group
DMRG is a variational method [71,82] which can be understood in the language of matrix product states (MPS) [83]. The MPS matrix dimensions, or the so-called bond dimensions, indicate the number of states kept in the reduced Hilbert space, and play central roles in the approximation. A general many-body state can be represented by an MPS with exponential growth of the bond dimension from the edges. In practice one restricts the maximum value of the bond dimension, thus limiting the maximum entanglement allowed in the variational state. Ground states of local Hamiltonians of physical interest generally have low entanglement. DMRG minimizes the energy in this low-entanglement Hilbert space. The accuracy of DMRG can be systematically improved by increasing the bond dimension. Although DMRG is naturally formulated and most powerful for one-dimensional systems, it is now widely applied to 2D systems [84], and remains one of the most-accurate numerical methods in 2D.
In this work, we employ two DMRG schemes with different conserved quantum numbers and using different update schemes. The first scheme conserves only the S z tot with U (1) symmetry, and uses the two-site update in the optimization. This scheme is used when a pairing field is applied to the system, breaking the particle number conservation. In such systems the particle numbers are controlled by the chemical potential. This scheme efficiently enables fluctuations between different quantum numbers in the optimization and is less likely to be stuck in a local minimum. The second scheme [85] conserves both the U (1) total particle number and SU (2) spin symmetries, and uses the single-site update [86,87]. The single-site update is faster than the two-site update and thus allows us to achieve large bond dimension. This scheme is used for systems without pairing fields, which thus conserve total particle number. Since the truncation error is ill-defined in the single-site update, we use the two-site energy variance in the standard extrapolations [88]. The number of states kept in these systems is up to 30000 SU (2) states, which corresponds to ∼ 90000 U (1) states, providing the best accuracy attained to date to our knowledge.

B. Two different ways to characterize superconducting correlation
To study the superconducting properties in the ground state, we use two different probes: pair-pair correlation functions and the pairing order parameter. These are both defined in terms of the pairing operator of a pair of nearest-neighbor sites, i and j: We will compute the pair-pair correlation function and the pairing order parameter where · · · denotes expectation with respect to the many-body ground state. The pair-pair correlation function in Eq. (3) can be obtained directly in a calculation working in a sector with fixed particle numbers. From it the d-wave pairing correlation function, P d (i − i ), can be constructed as a function of pair separation (i − i), by considering all j in ij and all j in i j .
The pairing order parameter in Eq. (4), on the other hand, requires a different approach. We add a term in the Hamiltonian describing SC pairing fields [77,90] applied to the system:Ĥ where the amplitude of h ij p is given by the parameter h p , and the sign of h ij p is positive if the bond (i, j) is vertical (alongŷ-direction) and negative otherwise (alonĝ x-direction), in order to probe pairing order of structure d x 2 −y 2 [91,92].
In AFQMC, we can obtain the superconducting pairing order parameter ∆ from total energy calculations, using the Hellmann-Feynman theorem: where |Ψ 0 (h p ) and E(h p ) are the ground-state wave function and energy of the Hamiltonian (Ĥ +Ĥ p ). We compute the derivative in Eq. (6) by finite difference , where δ is chosen to be sufficiently small to ensure that the error is smaller than our statistical error bar or targeted resolution. As h p → 0, the order parameter in the unperturbed ground state is obtained. This approach allows us to directly compute the pairing order parameter in AFQMC, which had not been possible before. We next use an example to illustrate the above approach to compute order parameters. We consider the antiferromagnetic (AFM) Neel order at half-filling. A staggered inducing field is applied to the periodic supercell of size L x = L y = L, with magnitude h m and alternating signs on the two sublattices. Because of the absence of the sign problem at half-filling, no constraint is needed in the AFQMC calculation, and the results are exact numerically. In Fig. 1(a), we show the computed staggered AFM order parameter M L (h m ) as a function of the applied field strength h m for different lattice sizes. Extrapolation to the thermodynamic limit (TDL) is then performed at each fixed h m , as illustrated in panel (b). The resulting TDL values are plotted in panel (c) versus h m , and extrapolated to the h m → 0 limit to obtain the order parameter. The result of 0.236(3) is in excellent agreement with the previous result of 0.236(1) computed from spin-spin correlation functions [89]. This test provides a validation of our approach for computing superconducting order parameters, which follows identical procedures. (We note that the staggered AFM magnetization in the repulsive Hubbard model at half-filling can be mapped to the s-wave on-site pairing order parameter in the attractive Hubbard model, through a partial particle hole transformation as discussed in the Appendix.) In the following two subsections, we show benchmark results on the two ways to compute the pairing order, respectively. Careful and detailed comparisons are made between AFQMC and DMRG, first for computing the pairing order parameter and then for pair-pair correlation functions, by using cylindrical geometries. Our results applying these approaches to address the physical properties of the Hubbard model are presented in Sec. III.

Pairing order parameter
The ground-state energies of 16 × 2 and 16 × 4 systems computed from AFQMC and DMRG are shown in the top panel of Fig. 2, as a function of the applied pairing field strength h p . Uniform "d-wave" pairing fields are applied to the entire system. A fixed value of µ is used which gives a doping of 1/8 at h p = 0 (µ = 1.75 for 16 × 4 and µ = 1.55 for 16 × 2). In AFQMC the trial wavefunctions are optimized self-consistently by coupling to natural orbitals [70]. The inset shows difference between the energies computed from AFQMC and DMRG. The relative error of the AFQMC energy is less than 0.5% for all h p in Fig. 2 which means the CP error is very small.
In the bottom panel of Fig. 2, we plot the pairing order parameters from AFQMC and DMRG, for the same system. In DMRG the order parameter is directly computed as a ground-state expectation value for each h p , while in AFQMC it is computed with the approach involving Hellman-Feynman theorem described above. Agreement between the two methods is excellent throughout the entire range. The general behavior of the order parameter is similar to that of the AFM order in Fig. 1    larger h p . The behavior is also manifested in the energy results as we show in the appendix: a fit of the energies at h p > h th p (where h th p is a threshold whose precise value does not affect the result), gives a E SC (h p = 0) which lies above the true ground state energy of the system.
We also show the comparison for the pairing order parameter at U = 4 and h = 1/6 in a 24 × 4 cylinder, in Fig. 3. Note that the pinning field range here is much smaller than in Fig. 2, focusing on the weak fields and a very fine scale of the pairing order parameter for comparison. With the lower value of U , this system requires larger bond dimensions in DMRG to converge, and we illustrate the extrapolation with truncation error. Good agreement is seen in the order parameters computed from AFMQC and the extrapolated results from DMRG.

Pair-pair correlation function
The pair-pair correlation function is computed with fixed number of particles (canonical ensemble). Results from CP-AFQMC have been obtained earlier in supercells with PBC using free-electron trial wave functions [35] and also using a BCS type of trial wave function after a particle-hole transformation [36]. Here we employ a more direct and general approach to apply projected BCS trial wave functions [81], and are able to access much larger systems because of algorithmic improvements and especially increased computing power. More unique to this work is the detailed and direct comparison with DMRG to quantify the accuracy.
In Fig. 4, we show a comparison of the pair-pair correlation function for 24 × 4, at doping of 1/8, with U = 4 in a cylindrical geometry between DMRG and AFQMC. The pair-pair correlations have a much smaller signal (roughly ∆ 2 ), as can be seen even in these small system sizes. Agreement is reasonable but the accuracy does not reach the level seen with the order parameter calculations. Hence the order parameter will be the primary tool on which we rely to accurately determine the nature of superconducting orders. This is why the development in this work of direct computation of the order parameter is crucial.

III. RESULTS
This section contains the following four parts. We first scan the pairing susceptibility versus doping h at a representative interaction strength of U = 8. We then carry out a detailed and systematic study of the pairing properties at 1/8 and U = 8 in Sec. III B. This is followed by an examination of the relation between stripe and SC orders in Sec. III C, and then an investigation of the dependence on U in Sec. III D. We first probe the SC response as a function of electron density, by computing the pairing order parameter in the presence of a d-wave pairing field, which is applied to the entire system. The pairing field amplitude h p is chosen to be of modest strength, so that a sizable SC order is induced but the field term does not drive the system far away from its ground state. The electron density n is controlled by the chemical potential µ. For the 16 × 4 cylinder, the µ value was varied in the range of 1.4 to 2.0 which yielded an electron density from ∼ 0.79 to 1. Fig. 5 shows the SC pairing order parameter as a function of density, or doping level. It can be seen from the 16 × 4 scan that the SC order has stronger response between n ∼ 0.81 and 0.92, with the maximum close to n = 0.885. Results in wider systems remain consistent, with the SC order showing slow variations in the vicinity of h = 1/8 doping. At density near the maximum SC order, the system displays charge and spin orders consistent with the ground state at 1/8 doping (n = 0.875), namely a stripe order [63]. This indicates that the system shows a SC order in response to the applied pairing field but remains in a similar ground state as the one when the pairing field is absent.
We have also investigated the doping dependence of the SC response in a 64 × 4 system using a different and complementary approach to the one in Fig. 5. A linearly varying chemical potential, µ(x), is applied along the cylinder, and the SC order and local density are computed without pairing field by allowing particle numbers to fluctuate in the DMRG calculation. The dependence of the SC order as a function of local density is found to be consistent with that in Fig. 5, as shown in the Appendix.
The fact that 1/8 doping is near the maximum response of the SC order for U = 8, and that the SC order shows rather weak dependence on the precise density, leads us to focus on the system of h = 1/8, U = 8, for which there is also detailed data on the spin and charge order, as well as ground-state energy, to compare with. The interaction strength of U = 8 is chosen as representative of the physically relevant regime. The results are presented in the next section. We begin this section by considering a 48 × 4 cylinder at h = 1/8, with a pairing pinning field applied on ver- 10 20 30 40 x − x tical bonds at the left edge only. We measure how the SC pairing order parameter ∆ ij decays as a function of distance from the left edge, which gives an indication of the behavior of the pairing correlation function in the bulk. We expect at least algebraic decay of the pairing order parameter if the system exhibits long range SC order, and exponential decay if there is no such order. The calculations are done with DMRG, without conserving particle number but with U (1) symmetry for S z tot . In Fig. 6 we show the SC pairing order parameter on the vertical (ŷ) bonds along thex-direction. The SC pairing order parameter is well converged when the bond dimension reaches m = 12000, so no extrapolation is needed. We perform both exponential and algebraic fits. The SC pairing order parameter clearly decays exponentially versus distance from the pinning field. As shown in the inset in Fig. 6, the pairing order on the horizontal (x) bonds are perfectly symmetric with the negative values of the vertical bonds. Although the pairing pinning fields are applied only on the vertical bonds at the left edge, the whole system spontaneously builds a d-wave pairing structure. This further confirms the tendency for shortrange d-wave pairing, however the exponential decay of the pairing order indicates that there is no long-range SC order in the filled stripes.
We next investigate the pair-pair correlation directly on a 48 × 6 cylinder. At h = 1/8 doping, the ground state has filled stripes with λ = 8 [63]. Previous study on width-2 ladders found that the SC pair-pair correlation decays algebraically [93], while further study on width-4 cylinder showed that the correlation decays exponentially [55]. Here we study a width-6 cylinder, employing the U (1) × SU (2) symmetry adapted DMRG with single-site update, and keeping bond dimension up to 22000 SU (2), to our knowledge the largest bond dimension to date in studying the SC pairing on width-6 cylinders. The computed SC pair-pair correlation is then extrapolated with respect to the two-site energy variance [88]. (The extrapolation details can be found in the Appendix.) Again we perform both exponential and the algebraic fits. As shown in Fig.7, the SC pair-pair correlations follow clearly an exponential decay with pair separation, showing no long-range order. The inset shows the (negative) values of the correlations on the vertical (horizontal) bonds. The correlations on the horizontal bonds are perfectly symmetric with the vertical bonds but with opposite sign, again confirming the d-wave structure.
On the width-6 cylinders, in addition to the filledstripes with λ = 8, the 2/3-filled stripes with λ ≈ 5 can also be stabilized, and has a slightly higher energy (∼ 0.001t) [63] than the ground state studied above. This state has wavelength closer to the stripes (λ = 4) observed experimentally [94], so it is interesting to also investigate the pairing in this meta-stable state.
We computed with DMRG the pair-pair correlation function in 48 × 6 cylinders at 1/8 hope doping in this state. As in the ground state with filled stripes, the results are shown in Appendix C. The pair-pair correlation is found to decay exponentially, even faster than in the ground state with filled stripes.
It is worth emphasizing that we have used two different DMRG schemes above, one under U (1) symmetry (S z tot ) with two-site updates, and the other under SU (2) × U (1) (spin and particle number) with strictly single-site updates. The consistency between the two approaches is a further confirmation of their reliability. The width of the systems which can be studied efficiently with DMRG is limited due to the linear increase of entanglement entropy with the width. To reach the TDL in two dimensions, we complement DMRG with two kinds of AFQMC calculations, computing both the pair-pair correlation function and the pairing order parameter, as described next.
The pair-pair correlation function in an 32 × 6 supercell with PBC along both directions is shown in Fig. 8. This calculation is performed with fixed number of electrons, at h = 1/8, with U = 8. The calculation with a free-electron trial wave function is consistent with earlier results from CP-AFQMC on square lattices. The calculation with a projected BCS trial wave function, as mentioned, employed a new method [81] which allows direct computation and back-propagation in the Hubbard model working in canonical ensemble. For both free electron and number-projected BCS trial wave-functions, the pair-pair correlations from AFQMC are seen to decay to 0 within the statistical resolution beyond a few lattice spacings.
We next employ AFQMC to calculate the pairing order parameter, by applyingĤ p as in Eq. (5), with the pairing fields chosen to match the d x 2 −y 2 structure and applied throughout the supercell. The pairing order parameter ∆ (averaged over all bonds) is calculated as a function of h p . The chemical potential µ = 1.75 is chosen such that the hole density is h = 1/8 in the ground state when h p = 0, and held fixed for all h p . To detect possible long-range SC pairing order in the pure Hubbard model (h p = 0), we need to reach the TDL at each h p first, then extrapolate h p to zero. This procedure is parallel to what is illustrated in Fig. 1, and is shown in Fig. 9. Following the procedure discussed in Sec. II B, AFMQC calculations are performed on various system sizes up to 32 × 8. We focus on the small h p region where the behavior determines whether long-range SC pairing order exists. In this region the self-consistent trial wave function gives the same results as the non-interacting trial wave-function; the latter form is used here, obtained by setting U = 0 in Eq. 1 and tuning the chemical potential to give a doping of h = 1/8. The computed pairing order parameters as a function of h p are shown in Fig. 9 (a). For the lengths studied here (L x > 16), the results are not sensitive to L x , as can be seen by comparing the 16 × 6 and 24 × 6 results, and also the 16 × 8 and 32 × 8 results. On the other hand, the order increases when the system becomes wider, although the dependence on system size becomes weak beyond L y = 6. To obtain the order parameters at the TDL, ∆ ∞ (h p ), a linear extrapolation with 1/L y is performed for each value of h p , as shown in Fig. 9(b). The resulting ∆ ∞ (h p ) and the statistical uncertainties, are shown in Fig. 9(c). A quadratic fit is then performed, which yields a value ∆ ∞ (0) = 0.003 (6), as indicated by the symbol at h p = 0. A linear fit for h p < 0.05 is also shown, which gives a statistically consistent result. We thus conclude that there is no long-range SC pairing in this system in the TDL.

C. The competition between stripe and superconducting orders
In Fig. 10, we examine the trend of the stripe and SC pairing orders when the strength of the applied pairing field is varied. The stripe order amplitude is defined as the intensity of hole modulation in the stripe state (i.e., the maximum value minus the minimum value of hole density along the longer direction of cylinder). Results are presented for two systems, a 16 × 4 cylinder computed by DMRG and a 32 × 8 cylinder by AFQMC. Both results are for U = 8, with h = 1/8. We can see that, when h p is decreased, the pairing order parameter becomes smaller, while the stripe order increases. The two orders thus compete against each other in the Hubbard model. The results of the previous subsection show that, at the zero pairing field limit, the stripe order dominates and no pairing order survives in this parameter region. In this section we study the pairing order at different interaction strengths U . In Fig. 11, we plot the pairing order parameter for U = 4, 6 and 8 at 1/8 doping. These calculations follow the same procedure as in Fig. 2. performed on a 16 × 4 system, with pairing fields of d-wave structure applied to the entire system. In the large h p region, we find that the SC pairing order parameter increases as U decreases. In the small h p region, the pairing order parameter varies little with the decrease of U .
Given the tendency for the pairing susceptibility to increase as U is reduced, we next focus on a lower but still physically relevant value of U = 4. In Fig. 12, we show the pair-pair correlation function for 1/8 doping from DMRG. This study at U = 4 is similar to the one in Fig. 7 at U = 8. Consistent with the result from Sec. III C, the pairing correlation amplitude is substantially larger than for U = 8. Both exponential and algebraic fits are performed to the correlation function computed with the largest bond dimension kept, m = 30, 000. Here the results are less definitive. The exponential is a slightly better fit to the data, but the algebraic fit cannot be ruled out conclusively. We next study the pairing order parameter at U = 4, using AFQMC to approach the TDL. We target h = 1/6, near optimal doping. In this parameter regime, no stripe or spin-density wave state is observed in the ground state of the Hubbard model [51]. This is the system where one of the extensive cross-checks between DMRG and AFQMC was performed in Sec. II B (Fig. 3). The computed pairing order parameter for a variety of system sizes are shown in Fig. 13. The same procedure as in Fig. 9 is performed. We first carry out an extrapolation with the width L y at each h p , to reach the TDL. In contrast to the U = 8 case, the pairing order parameter here is seen to either decrease or saturate very quickly with system width as L y grows. The results for width-4 cylinders (24×4) are seen to be non-monotonic with wider systems, so they are not included in the fit. The final result extrapolated to the h p = 0 limit is ∆ ∞ (0) = 0.006(4), statistically compatible with a vanishing order parameter.
Of course, based on this and the pairing correlation results above, we can not fully rule out the possibility of a finite pairing order in the ground state at U = 4. (See also discussion on weak coupling in Sec. I.) Our results do put a rather stringent bound on the strength of the pairing order, which is considerably smaller than indicated by the best previous calculations with affirmative results on superconductivity. The small magnitude of this bound suggests that, even if the pure Hubbard model is superconducting in some regime of the parameter space further from the most relevant physical parameterrs, it is likely missing key ingredients as a fundamental model for cuprate superconductors.

IV. SUMMARY AND PERSPECTIVE
In summary, we have carried out a detailed study of the superconducting pairing properties in the ground state of the 2D pure Hubbard model, using two of the most accurate ground-state many-body computational methods available at present. With both methods, we have presented technical advances which enabled new capabilities in probing the superconducting order. The DMRG calculations of pairing correlation functions were performed on up to width-6 cylinders, with unprecedented accuracy. The AFQMC computations were, for the first time, able to compute pairing order parameters relying on total energy calculations. Meticulous comparisons were made between the two methods. Their complementary application allowed us to maintain high accuracy and reach the thermodynamic limit.
In the parameter regime relevant to the cuprates (U ∼ 6-8) we find that the pure Hubbard model does not have a superconducting ground state. We also find that the lack of superconductivity is due to a competition with stripe order, with stripes dominating. At smaller U ∼ 4, the tendency for striped ground states is much weaker. In this case, we still find a pairing response consistent with zero. While we cannot rule out a small nonzero pairing order parameter, our results place an upper bound to its strength which is very small.
In the early days of high temperature cuprate superconductivity, when no numerical approach was adequate to accurately probe the low temperature, doped regime in large system sizes, it was natural to expect that if one could get past the fermion sign problem, one would quickly have a clear picture of the physics involved. Now that we can study this regime, new obstacles have been revealed. A key obstacle is the close competition of a number of different phases, with small Hamiltonian terms mediating which phase is favored. This situation makes simulations more difficult, but equally important is that it is very difficult to know reliably which sets of small terms and parameters (such as next nearest neighbor hopping t ) describe the actual materials. Our work can be viewed as a key initial step, where the iconic simplestto-define model with only U and t is found not to exhibit superconductivity. To go beyond this, one will want to study phases and superconductivity in generalized models including a broad range of parameters. Simultaneously, it is important to improve our techniques for both deriving accurate models, and for simulating real systems with very strong correlation without introducing models.
It is also important to note that the pure Hubbard model does get much of the physics right, including antiferromagnetism and its destruction upon doping, and a tendency for stripes to occur and to compete with d-wave superconductivity. It also produces the crucial physics that there are many intertwined states separated by tiny energy scales, a key part of the reason that the complete nature of superconductivity in the cuprates remains to be resolved.  Fig. 9. In (a) the pairing order parameter computed from DMRG for 24 × 4 (see Fig. 3) are also shown for comparison. In (b), the width-4 data were excluded from the fits to the TDL. In (c), the extrapolation yields a value ∆∞(0) = 0.006(4), as indicated by the symbol at hp = 0.
In Fig. 14 we show the comparison of energies between pairing state and stripe state. The systems are 16 × 4 cylinders with d-wave pairing field applied on the whole systems. The energy for pairing state at h p = 0 (denoted by triangles in Fig. 14) is obtained from a quadratic fit with energies at large h p . The stripe energy (denoted by open square and cycle in Fig. 14)) is the actual value calculated at h p = 0. We can find the energy of pairing state is slightly higher than that of stripe state, by about ∼ 0.01 per site. chemical potential µ(x) linearly changing from 1.4 to 2.1 along the longitudinal (x) direction. Since the local density will vary with x, we obtain information about different dopings in a single calculation. We use DMRG with-out conserving total particle number to allow SC orders to develop. The local densities and the local SC orders along the longitudinal direction are shown in Fig. 15. The SC orders on the horizontal and vertical bonds are symmetric with opposite signs, showing the d-wave symmetry. The maximum SC order appears around the density n ≈ 0.875 (1/8 hole doping), and decays for both higher and lower densities. This result is further confirmation of the doping dependence of the SC order observed in Sec. III A, and validates the choice of doping h ∼ 1/8 as a representative case in studying the SC response. The optimal doping µ opt ≈ 1.73 is actually at the boundary between two different stripe fillings. For µ < µ opt the ground states are filled stripes (4 holes per stripe in the width-4 cylinder) and for µ > µ opt the ground states are half filled stripes (2 holes per stripe). This can be seen from the more abrupt change in density at µ opt (top panel in Fig. 15), and was further confirmed by our calculations with uniform chemical potentials (not shown here). This picture is consistent with the idea that fluctuations between different stripe fillings can help induce SC orders.
Note that, strictly speaking, the SC order should be zero here in a finite-size system, since the Hamiltonian does not break (total) particle number conservation without an applied pairing field. However in the DMRG calculation the variational ground states often break the symmetry due to the finite bond dimensions kept. This feature has been used in the past to study the magnetization and now the SC pairing order.
Appendix D: Filled and 2/3-filled stripes on 6-leg cylinders Besides the filled stripes, we also considered the 2/3filled stripes on width-6 cylinders. The filled and 2/3filled stripes are the only two striped states that can be stabilized on width-6 cylinders in DMRG. In Fig. 16 we show the pair-pair correlations for the 2/3-filled stripes on a 48 × 6 system. The correlations for both the finite bond dimension m as well as the infinite m are shown. The detail of the extrapolation will be shown in the next section. Both the power-law and the exponential fittings are shown. As in the filled stripes, the correlations decay exponentially with distance. The inset shows the absolute values of the correlations on both the vertical bonds and the horizontal bonds. The correlations on the horizontal bonds again are perfectly symmetric with the vertical bonds at the same location but with opposite sign, consistent with the d-wave symmetry.
In Fig. 17 we show the linear extrapolation of the energies with the two-site variance for both the filled and the 2/3-filled stripes. The clearly linear behaviors typically indicate the stability of the MPS toward the zerovariance limit. In other words, the MPS basically stays in the same state for the considered bond dimensions.
The crossing of the lines shows that the filled stripes is lower in energy only when the bond dimension is sufficiently large. We also point our that the filled stripes need larger bond dimension than the 2/3-filled stripes to achieve the same variance, because it contains higher entanglement. Here we discuss some details of the DMRG simulations. In the simulations preserving SU (2) symmetry, i.e., the simulations on the systems without any pairing field, temporary local chemical potentials are applied in the first few sweeps on the expected locations of the stripes to stabilize the states and improve convergence. One often also applies the magnetic pinning field to stabilize the stripes; however in our cases the magnetic field will break the SU (2) symmetry, so no magnetic pinning field is used. For the filled and the 2/3-filled stripes, the temporary chemical potentials are applied up to m = 4000 (m = 1400), and then are switched off for the further sweeps of larger m.
For the same systems, the single-site update is used in DMRG. To eliminate the finite bond-dimension effect, we employ the extrapolations of physical quantities (for example energy and pair-pair correlations) with the two-site energy variances [88]. In two-site DMRG, one usually extrapolates the physical quantities with the socalled truncation error (or alternatively called discarded weight). However in single-site DMRG, the truncation errors are not well defined. We thus extrapolate the physical quantities with the two-site energy variance, which is an approximation of the full variance (Ĥ − E) 2 /N 2 . Physically the variance is a perfect quantity to extrapolate with since it measures the distance of the variational state to an eigenstate. In practice this extrapolation scheme was demonstrated to be as reliable as the extrapolation by the truncation error [88].
At the largest bond dimension we can achieve, which to our knowledge is also the largest bond dimension that has been done to date, the pair-pair correlation vs. energy variance is not yet reach the linear region. We thus perform cubic extrapolations to best fit the data, as shown in Fig. 18 and Fig. 19. Since the correlations on different distance have quite different scales, we show only for part of the distance which is long enough but away enough from the boundary. Although here we show only the cubic extrapolations, we also tested other extrapolations, for example the linear extrapolation of the last three data points, and the results are very similar to what is presented and lead to the same conclusion (exponential decays). x'-x=20 x'-x=21 x'-x=22 x'-x=23 x'-x=24 x'-x=25 x'-x=26 x'-x=27 x'-x=28 x'-x=29