Stripe and superconducting order competing in the Hubbard model on a square lattice studied by a combined variational Monte Carlo and tensor network method

The long-studied Hubbard model is one of the simplest models of copper-oxide superconductors. However, the connection between the model and the experimental phase diagram is still under debate, in particular regarding the existence and extent of the $d$-wave superconducting phase. Recent rapid progress in improving the accuracy of numerical solvers has opened a way to answer this question reliably. Here, we study the hole-doping concentration ($\delta$) dependence of the Hubbard model in the ground states on a square lattice at strong coupling $U/t=10$, for the on-site interaction $U$ and the transfer $t$, using a variational Monte Carlo method. The method, which combines tensor network and Lanczos methods on top of Pfaffian wave functions, reveals a rich phase diagram, in which many orders compete severely and degenerate within the energy range of $0.01t$. We have identified distinct phases including a uniform $d$-wave superconducting phase for $0.17\lesssim \delta \lesssim 0.22$ and a stripe charge/spin ordered phase for $\delta \lesssim 0.17$ with the stripe period depending on $\delta$, together with presumable spatially coexisting antiferromagnetic and stripe order for $\delta \lesssim 0.07$ and coexisting stripe and $d$-wave superconductivity for $0.07\lesssim \delta \lesssim 0.17$. The present, improved method revealed a wider region of a charge uniform superconducting phase than the previous studies and shows a qualitative similarity to the phase diagram of the cuprate superconductors. The superconducting order parameter is largest at doping of around $\delta =0.17$ in the ground state, which undergoes phase transitions from an inhomogeneous to a uniform state.

Figure 1

(a) Examples of hole density and $\langle S_z\rangle$ in a $2\lambda \times 2$ cell for an odd and even period stripe obtained by optimization at $U/t=10$. Hole density is proportional to the diameter of the circle, while $\langle S_z\rangle$ is proportional to the length of the arrow. (b) Energies of stripe and uniform states per site in units of $t$ as a function of doping for $L_y=8$ systems ($U/t=10$). $L_x$ is indicated in the legend. As finite-size effects are found to be small (see Appendix pp2) we note that the result here essentially represents the thermodynamic limit. A linear function $f\left(\delta \right)=1.835\delta +0.4304$ has been added to the energy in (b) to improve visibility. (c) Ground-state phase diagram with phase separation (PS), stripe spin and charge order coexisting with weak $d$-wave superconductivity (SO $+$ weak SC), charge uniform $d$-wave SC (SC) and charge uniform paramagnetic normal (PM) regions indicated. Checkerboard region indicates a possible region of phase separation between SO and uniform SC phases. (d) Thermodynamic limit extrapolations of order parameters of the ground state in different regions: peak of charge structure factor is plotted for the SO region, and ${\mathrm{\Delta }}_{\mathrm{SC}}$ is plotted in the PM region. Charge structure of several different stripe periods is plotted simultaneously as different stripe periods are near degenerate.

Figure 2

(a) Peak values of spin structure factor and (b) correlation ratio as a function of doping $\delta$. In (a) dots are from finite-size ($L_y=16$) calculations, while cross symbols are obtained as values after the TL extrapolation. (a) shows results with the tensor network correlation factor, while (b) is calculated without it due to computational cost.

Figure 3

Variance extrapolation of the energy per site in units of $t$ at $\delta =0.0625$ on a $16\times 8$ system at $U/t=10$. Points on this plot from highest to lowest energy are the VMC wave function without TN, VMC wave function with TN and $D=2$, VMC wave function without TN with first-step Lanczos applied and VMC wave function with $D=2$ TN and first-step Lanczos.

Figure 4

Comparison of energy per site (in units of $t$) vs doping plot at $U/t=10$ for two different system sizes: (a) with $L_y=8$ and (b) with $L_y=16$. Essential features of the diagram appear unaffected by increasing size beyond $L_y=8$. A linear function $f\left(\delta \right)=1.835\delta +0.4304$ has been added to the energy to improve visibility.

Figure 5

System-size extrapolation for (a) spin structure factor peak and (b) charge structure factor peak for a $\lambda =7$ stripe at $\delta =0.142$ and $U/t=10$. Error bars are smaller than marker sizes. Extrapolated values for infinitely long $L_y=8$ systems and $L_y=16$ systems are essentially indistinguishable. For these system sizes, the one and two dimensional extrapolations with respect to $N_s$ are nearly identical, therefore, we have taken the one dimensional extrapolation with fixed $L_y=8$ as the TL extrapolated values of $S_s$ and $S_c$.

Figure 6

Size extrapolation of ${\mathrm{\Delta }}_{\mathrm{SC}}$ in charge uniform states. Calculations were performed at $U/t=10$ and the system size is varied. At $\delta =0.125$, only square systems with $L_x=L_y$ were used, however, at other dopings, nonsquare lattices were needed (since square systems could not support the specified doping for an integer number of electron pairs). The system dimension is labeled for $\delta =0.1875$. The same dimensions were used for $\delta =0.2188$. For $\delta =1875$ and 0.2188, values obtained with antiperiodic-periodic and periodic-periodic boundary conditions were averaged to reduce finite-size effects. The tensor-network bond dimension was fixed to $D=2$, as we observed little change with increasing $D$ beyond this.

Figure 7

Size dependence of ${\mathrm{\Delta }}_{\mathrm{SC}}$ in stripe states at $U/t=10$. Although ${\mathrm{\Delta }}_{\mathrm{SC}}$ does not vary smoothly with system size, a linear extrapolation is plotted, providing rough approximations for the TL values of ${\mathrm{\Delta }}_{\mathrm{SC}}$. The tensor network factor is not used in these calculations due to computational cost. Although the TL value of ${\mathrm{\Delta }}_{\mathrm{SC}}$ is not precisely determined, it is clearly less than ${\mathrm{\Delta }}_{\mathrm{SC}}$ in uniform states in the region $0.125\le \delta \le 0.1875$ (see Fig. 6).

Figure 8

Comparison of extrapolated energies for $U/t=8$ at half-filling with numerically exact QMC energies for various system sizes in units of $t$. The green line is without quantum number projections applied. The red line is obtained with ${\mathcal{L}}^K$ and ${\mathcal{L}}^S$ applied to ${\psi }_{\mathrm{pair}}\left(x\right)$ and ${\mathcal{L}}^K$ to $\mathcal{M}\left(x\right)$. $2\times 2$ unit cell is used in the calculations without quantum number projections, while the full unit cell is used when they are applied.

Figure 9

Comparisons of different physical quantities with and without quantum number projections ${\mathcal{L}}^{\mathcal{K}}$ and ${\mathcal{L}}^{\mathcal{S}}$ applied to ${\psi }_{\mathrm{pair}}\left(x\right)$, which are denoted by “$S=0$, $k=0$” and “no proj” in the graph legends, respectively. Structure factor for (a) charge, (b) spin and (c) superconducting correlations in the $x$ direction for a period 7 stripe at $U/t=8$ at $\delta =0.125$ on a $14\times 8$ system. (d) Spin structure factor at half filling on a $8\times 8$ system at $U/t=8$.