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—constrained-path (CP) 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 nonsuperconducting in the moderate to strong coupling regime in the vicinity of optimal hole doping.

Popular Summary

Understanding high-temperature superconductivity in copper oxide materials known as cuprates is one of the greatest challenges in theoretical condensed-matter physics. For more than 30 years, the leading mathematical model of this phenomenon has been the so-called 2D Hubbard model, a relatively simple model of particles interacting in a lattice. While the model does a reasonable job of describing many properties of cuprates, superconductivity is itself a bit more delicate. It remains to be seen whether the ground state of this model is superconducting in the relevant regime. Here, we employ advances in state-of-the-art algorithms to explore this question and find the ground state of this model in its simplest form is not superconducting.

In our analysis, we combine two complementary techniques for computing behavior of many-body systems. One, known as the density matrix renormalization group, provides nearly exact calculations for systems best described as narrow, long cylinders. We then use that technique to validate another method, known as the auxiliary-field quantum Monte Carlo method, which we find to be very accurate on larger systems. We find that the ground state of the Hubbard model is nonsuperconducting for the case where particle interactions are moderate to strong and carrier concentrations are intermediate, a regime relevant to superconducting materials.

Further research will need to extend these results to a wider range of model parameters to establish the full superconducting phase diagram of the family of Hubbard models.

Figure 1

Illustration of the approach to compute order parameters. The AFM order parameter at half-filling ($U=4$) is computed by applying staggered magnetic fields to the periodic supercell. (a) The computed magnetic order parameter $M$ is shown as a function of the pinning field strength $h_m$ for different system sizes. (b) A quadratic fit is performed at each $h_m$ to extrapolate $M(h_m)$ to the TDL. The procedure is shown for two values $h_m=0.06$ and 0.1, as marked by the vertical dashed lines in panel (a). (c) An extrapolation of $M_{\infty }(h_m)$ to the $h_m\to 0$ limit is then performed, again using a quadratic fit. The resulting $M_{\infty }(0)$ is shown by the open symbol. The value is in excellent agreement with the order parameter determined from direct computation of spin-spin correlation functions [92], shown by the red star.

Figure 2

Upper panel: Comparison of the ground-state energies computed from AFQMC (blue) and DMRG (red), as a function of the applied pairing field strength $h_p$, at $U=8$. Two cylindrical systems are shown, $16\times 2$ and $16\times 4$, with the chemical potential held fixed in each so that the doping is $1/8$ when $h_p\to 0$. The inset shows the difference of the energy computed from AFQMC with respect to DMRG. Lower panel: Comparison of the computed SC pairing order parameter for the same systems.

Figure 3

Comparison of the computed pairing order parameter at $U=4$, as a function of the applied pairing field strength $h_p$. The system is a $24\times 4$ cylinder, at $h\simeq 1/6$. DMRG results, with $U(1)$ symmetry and a two-site update, with different bond dimensions are shown. The results after extrapolation to zero truncation error are also plotted. The inset illustrates the extrapolation (with respect to the truncation error) at $h_p=0.003\times \sqrt{2}\approx 0.00424$. The red dot represents the extrapolated value from DMRG, while the blue dot with the error bar is the AFQMC result.

Figure 4

Comparison between DMRG and AFQMC of the pair-pair correlation function for a $24\times 4$ cylinder, with $U=4$ at $1/8$ doping. The reference bond is [(12,1), (12,2)]. Vertical and horizontal bonds are along the $y$ and $x$ directions, respectively.

Figure 5

SC order parameters as a function of doping level, with $U=8$. Results from three systems are shown, each at a modest value of the applied $d$-wave pairing field (with strength $h_p$ indicated in the legend). The doping level or particle density is controlled by varying the chemical potential.

Figure 6

SC pairing order parameter on the vertical bonds at a fixed $y$, versus position $x$, for a $48\times 4$ cylinder at $1/8$ hole doping, computed from DMRG with $U(1)$ symmetry ($S_{\mathrm{tot}}^z$) and a two-site update. Pairing fields with $h_p=0.25$ are applied to the vertical bonds at the left edge $x=0$. Results with different bond dimension $m$ are shown. Both exponential (red line) and algebraic (blue curve) fits are shown; the solid (dashed) region indicates the region (not) used in the fits. The oscillation of the pairing order parameter coincides with the stripe period. The inset shows the (negative) pairing order parameters on the vertical ($\parallel$) [horizontal ($\perp$)] bonds for $m=12000$.

Figure 7

SC pair-pair correlation on the vertical bonds at a fixed $y$, versus pair separation, for a $48\times 6$ cylinder at $1/8$ hole doping, computed from DMRG with $U(1)\times SU(2)$ symmetry and single-site update. Results of different bond dimensions $m$, as well as the extrapolation to the infinite bond dimension, are shown. The pair separation $(x^{\prime }-x)$ is measured with respect to the reference bond, a vertical bond at $x=5$. Both the exponential (red line) and the algebraic (blue curve) fits are shown; the solid (dashed) region indicates the region (not) used in the fits. The correlation length is approximately 2.9 from the fitting. The inset shows the (negative) correlations on the vertical ($\parallel$) [horizontal ($\perp$)] bonds for $m=22000$.

Figure 8

Pair-pair correlation from QMC for a $32\times 6$ system with $U=8$, $1/8$ doping and PBC. Different trial wave functions are used. In the inset, the results from trial wave functions themselves are shown. Comparing to the variational result in the BCS wave function, the pair-pair correlation in the QMC calculation, using it as a trial wave function, is largely suppressed.

Figure 9

Finite-size scaling of the SC pairing order parameter, at $U=8$ and $h=1/8$. In panel (a), the computed pairing order parameters $\mathrm{\Delta }$ are shown for a variety of system sizes in the small-$h_p$ region. In panel (b), an extrapolation of the results in panel (a) is performed with respect to the width in the periodic direction (linear in $1/L_y$) for each value of $h_p$. The lines, from top to bottom, are for $h_p=0.064$, 0.049, 0.035, 0.021, 0.014, and 0.0078, respectively. The resulting pairing order parameter ${\mathrm{\Delta }}_{\infty }(h_p)$ is shown on the left. (A slight shift in the horizontal position has been applied to some of the data points for better visibility of the results and error bars.) In panel (c), the result from the fit in panel (b), ${\mathrm{\Delta }}_{\infty }(h_p)$, is plotted versus $h_p$. A quadratic fit is then performed, which yields a value ${\mathrm{\Delta }}_{\infty }(0)=0.003(6)$, as indicated by the black star at $h_p=0$. A linear fit of the last few points is also shown. Weighted least-squares fits are used to account for the statistical errors.

Figure 10

Strengths of the stripe and pairing orders versus the applied pairing field, $h_p$. The upper and lower panels show results from $16\times 4$ and $32\times 8$ cylinders, respectively, both at $U=8$ with $h=1/8$.

Figure 11

Dependence of the pairing order parameter on interaction strength $U$. Calculations were done in the same manner as in Fig. 2, on $16\times 4$ systems, varying only $U$. In the inset, the pairing order parameters are plotted versus $U$ for two values of $h_p$ as indicated by the vertical lines in the main graph, $h_p=0.38$ and 0.58.

Figure 12

SC pair-pair correlation on the vertical bonds for $U=4$ at $h=1/8$ on a $48\times 4$ cylinder, computed by DMRG with $U(1)\times SU(2)$ symmetry and the single-site update. The reference bond is a vertical bond at $x=5$. Both exponential (red line) and algebraic (blue curve) fits are shown; the solid (dashed) region indicates the region (not) used in the fits.

Figure 13

Finite-size scaling of the SC pairing order parameter, at $U=4$ and doping of $h\simeq 1/6$. The layout is the same as in Fig. 9. In panel (a), the pairing order parameter computed from DMRG for $24\times 4$ (see Fig. 3) are also shown for comparison. In panel (b), the width-four data (black symbols) were not included in the linear fits, as they are nonmonotonic with larger $L_y$ results. In panel (c), the extrapolation yields a value ${\mathrm{\Delta }}_{\infty }(0)=0.006(4)$, as indicated by the symbol at $h_p=0$.

Figure 14

Comparison of the ground-state energies from CP-AFQMC (red) and DMRG (blue), as a function of the pairing field strength $h_p$. The system is a $16\times 4$ cylinder. A fixed value of $\mu =1.75$ is used with which the doping is $1/8$ when $h_p=0$. A quadratic fit is applied to each set of energies using only values at large $h_p$ (points to the right of the vertical bar). The triangular symbols at $h_p=0$ show the intercept result from the fit, while the open symbols are those obtained from actual calculations done at $h_p=0$. In the inset, a zoom of the main plot near $h_p=0$ is shown.

Figure 15

Comparison of the local hole densities on different rungs for $16\times 4$ systems with $h_p=0.05$ and with zero pairing field. The system with the pairing field is the same as in Fig. 5 with $\mu =1.72$.

Figure 16

The system is a $64\times 4$ cylinder with local chemical potential $\mu$ linearly changing along the longitudinal direction. (Upper panel) Particle densities along the longitudinal direction. The $x$ axis has been replaced by the local $\mu$ on the corresponding position. (Lower panel) The local SC orders on the vertical (blue) and horizontal (red) bonds along the longitudinal direction. The results are obtained by using DMRG without conserving total particle number (grand canonical ensemble).

Figure 17

The density dependence on $\mu$ for a $16\times 4$ cylinder with pairing field $h_p=0.25$ at the left edge. The numbers in the parentheses are the number of holes in each stripe, from left to right. For example, (4,4) means two filled stripes (with four holes). The red color indicates the mixture of different stripe fillings.

Figure 18

SC pair-pair correlation for the $2/3$-filled stripes with $1/8$ hole doping on $48\times 6$ cylinders. Different curves are for different bond dimensions $m$, as well as the extrapolation to the infinite bond dimension. The pair separation $(x_0-x)$ is measured with respect to the reference bond, a vertical bond at $x=5$. Both exponential (red line) and algebraic (blue curve) fits are shown; solid (dashed) region indicates the region (not) used in the fits. The correlation length is approximately 1.9 from the fitting. The inset shows the (negative) correlations on the vertical ($\parallel$) [horizontal ($\perp$)] bonds for $m=15000$.

Figure 19

Linear extrapolation of ground-state energy with the two-site energy variances for the filled and the $2/3$-filled stripes on $48\times 6$ cylinders. The MPS bond dimensions shown for the filled stripes are from 8500 to 22 000, and for the $2/3$-filled stripes, they are from 7000 to 15 000. The extrapolated energy is $-0.7581(6)$ for the filled stripes and $-0.7574(4)$ for the $2/3$-filled stripes.

Figure 20

Pair-pair correlations (dots) and the cubic extrapolations (curves) using the two-site variance on the vertical bonds (left panel) and the horizontal bonds (right panel). The reference bond is the vertical bond on $x=5$. The system consists of the filled stripes on a $48\times 6$ cylinder as in Fig. 7. The MPS bond dimensions $m$ are from 10 000 to 22 000.

Figure 21

Pair-pair correlations (dots) and the cubic extrapolations (curves) by the two-site variance on the vertical bonds (left panel) and the horizontal bonds (right panel). The reference bond is the vertical bond on $x=5$. The system consists of the $2/3$-filled stripes on a $48\times 6$ cylinder as in Fig. 18. The MPS bond dimensions $m$ are from 7000 to 15 000.