Pairing in the two-dimensional Hubbard model from weak to strong coupling

The Hubbard model is the simplest model that is believed to exhibit superconductivity arising from purely repulsive interactions and has been extensively applied to explore a variety of unconventional superconducting systems. Here we study the evolution of the leading superconducting instabilities of the single-orbital Hubbard model on a two-dimensional square lattice as a function of onsite Coulomb repulsion $U$ and band filling by calculating the irreducible particle-particle scattering vertex obtained from dynamical cluster approximation (DCA) calculations, and compare the results to both perturbative Kohn-Luttinger (KL) theory as well as the widely used random phase approximation (RPA) spin-fluctuation pairing scheme. Near half-filling, we find remarkable agreement of the hierarchy of the leading pairing states among these three methods, implying adiabatic continuity between weak- and strong-coupling pairing solutions of the Hubbard model. The $d_{x^2-y^2}$-wave instability is robust to increasing $U$ near half-filling as expected. Away from half-filling, the predictions of KL and RPA at small $U$ for transitions to other pair states agree with DCA at intermediate $U$ as well as recent diagrammatic Monte Carlo calculations. RPA results fail only in the very dilute limit, where it yields a $d_{xy}$ ground state instead of a $p$-wave state established by diagrammatic Monte Carlo and low-order perturbative methods, as well as our DCA calculations. We discuss the origins of this discrepancy, highlighting the crucial role of the vertex corrections neglected in the RPA approach. Overall, a comparison of the various methods over the entire phase diagram strongly suggests a smooth crossover of the superconducting interaction generated by local Hubbard interactions between weak and strong coupling.

Figure 1

Second-order screening (bubble) and exchange (ladder) diagrams. Note that each interaction line $U$ depicted by a wiggly line connects opposite spins only. The bubble diagram contributes to same spin triplet pairing only, and singlet and opposite spin-triplet solutions arise from the ladder diagram after symmetrization and antisymmetrization, respectively, as stated in Eq. (6).

Figure 2

(a) Fermi surface and (b) real part of bare spin susceptibility ${\chi }_0\left(\mathbf{q}\right)$ for electron filling $n=0.90$ and $t^{\prime }=-0.15$.

Figure 3

Relative eigenvalues ${\lambda }_i/{\lambda }_1$ for subleading instabilities $i\in \{d\left(12\right),g\left(8\right),p^{\prime }\}$ as a function of interaction $U$ at density $n=0.90$. The eigenvalues are derived from Kohn-Luttinger (dashed-dotted lines) and RPA (full lines) in the regime $U=0-1.36$ (blue area) and from DCA for $U=2,4,6,\text{and}8$ (pink area). ${\lambda }_1$ refers to the leading eigenvalue, which always belongs to the $d\left(4\right)$ solution. For clarity, we have omitted subleading singlet instabilities appearing between $g\left(8\right)$ and $p^{\prime }$ for $U\ge 2$.

Figure 4

[(a)–(c)] Three leading solutions of the RPA model ($U=1$) with the leading $d_{x^2-y^2}$ solution, the sub-leading higher order $d_{x^2-y^2}$ solution with 12 nodes and third leading $g$-wave solution. [(d)–(f)] The results of the DCA calculation, ${\varphi }_{\alpha }(\mathbf{K},i{\omega }_1)$ at the lowest frequency are shown as insets. Comparison of RPA and DCA results by angle dependent gap plots in the first quadrant of the Brillouin zone. The angle $\theta$ is defined in the inset in $\left(\mathrm{d}\right)$. The three leading solutions of the RPA calculation ($U=1$) are shown in color and the DCA results ($U=2$) are shown by black symbols. The DCA results are interpolated to the Fermi surface of the noninteracting system to allow a direct comparison although the DCA calculation is only done for $8\times 8$ discrete values.

Figure 5

(a) Symmetry of the leading superconducting instability as a function of $U$ and filling $n$ with $t^{\prime }=0$ and $T=0.015$. RPA results are indicated by filled colors superimposed in the ground-state phase diagram obtained in Ref. [24] shown by open symbols and lines. The blue line displays the phase boundary between a simple harmonic $p$ wave $[sin\left(k_x\right)/sin\left(k_y\right)]$ and $d_{xy}$, and the dark green line indicates the phase boundary between $d_{xy}$ and $d_{x^2-y^2}$. The yellow region indicates an $s^{\prime }$ phase [40]. A region of higher-order $p^{\prime }$ wave (sketch as inset [24]) is visible in both approaches for small values of $U$ around a filling of $n=0.55$. The triplet $p^{\prime }$ state features six or ten nodes depending on temperature. In RPA, the region of leading triplet instability at very low filling is sensitive to resolution and temperature; it essentially vanishes as $T$ is decreased to $T=0.0002$ (small blue arrows and blue dashed line almost coinciding with the $y$ axis). (b) Zoom in of the phase diagram close to $n=0.55$ and comparison to Ref. [24]. The insets show solutions $g\left(\mathbf{k}\right)$ as obtained from RPA at the parameters marked by white points. The dashed red line (and the small red arrows) shows how the boundary between the $p^{\prime }$ and $d_{xy}$ changes within RPA when $T$ is lowered from $T=0.015$ to 0.0002.