Understanding repulsively mediated superconductivity of correlated electrons via massively parallel density matrix renormalization group

The so-called minimal models of unconventional superconductivity are lattice models of interacting electrons derived from materials in which electron pairing arises from purely repulsive interactions. Showing unambiguously that a minimal model actually can have a superconducting ground state remains a challenge at nonperturbative interactions. We make a significant step in this direction by computing ground states of the 2D U-V Hubbard model—the minimal model of the quasi-1D superconductors—by parallelized DMRG, which allows for systematic control of any bias and that is sign-problem-free. Using distributed-memory supercomputers and leveraging the advantages of the U-V model, we can treat unprecedented sizes of 2D strips and extrapolate their spin gap both to zero approximation error and the thermodynamic limit. Our results for the spin gap are shown to be compatible with a spin excitation spectrum that is either fully gapped or has zeros only in discrete points, and conversely that a Fermi liquid or magnetically ordered ground state is incompatible with them. Coupled with the enhancement to short-range correlations that we find exclusively in the $d_{xy}$ pairing channel, this allows us to build an indirect case for the ground state of this model having superconducting order in the full 2D limit, and ruling out the other main possible phases, magnetic orders, and Fermi liquids.

Figure 1

(a) Overview of the 2D U-V model. Conduction electrons tunnel along effective 1D chains with amplitude $t$, where electrons are subject to strong on-site repulsion $U$ and a weaker nearest-neighbor repulsion $V$. Interchain tunneling $t_{\perp }$ is much weaker than $t$. (b) Schematic parameter space of the 2D U-V model at half filling, for fixed $U=4t$. The only limit in which all possible ground states and the transition between them is understood is the 1D limit, i.e., at $t_{\perp }=0$, when the chains decouple, and are either in a TLL or MI phase, depending on whether $V<V_c$ or $V>V_c$. The red dots show the parameters at which we work in this present paper.

Figure 2

U-V parameter space. Green shaded area: single chain enters MI regime. Red shaded area: two coupled chains can enter MI regime [50]. To be certain we are outside the regime where this counter-intuitive increase in parameter space for the MI insulator takes place in the few-chain regime, we work within the blue-shaded area of parameter space. For reference, the yellow shaded area shows the parameter regime realized in the Bechgaard and Fabre salts. Black lines show lines of constant $K_{\rho }$ for the single U-V chain (the data and figure reprinted from Ref. [51]).

Figure 3

(a) Mapping the physical 2D lattice to the 1D DMRG chain, along the $t$ direction. Many more long-range tunneling terms appear than would be present when mapping along the $t_{\perp }$ direction. (b) How the two different mappings perform: mapping along $t$ (blue) rapidly converges to the exact solution, while $t_p$ mapping is trapped in a metastable state, even at $\chi =5000$ and many sweeps. (c) The primary parallelization layer of pDMRG used for the present work, using the left boundary $L$ as an example. Each element $L\left[b_l\right]$ of $L$ is dispatched to a MPI process, and each such MPI process may control a number of physical nodes. The process then locally handles contractions of $L\left[b_l\right]$ with the wave-function tensors to be optimized, requesting the results of other MPI processes for different $b_l$ as necessary.

Figure 4

[(a) and (b)] Averaging eliminates oscillations in the spin gap. (a) Comparing ${\mathrm{\Delta }}_s[L,N_{\mathrm{ch}}=2]$ (blue line) against ${\overline{\mathrm{\Delta }}}_s[L,N_{\mathrm{ch}}=2]$ (red line) for $U/t=4,V/t=1,t_{\perp }/t=0.1$. Using the averaged version of the spin gap, ${\overline{\mathrm{\Delta }}}_s$, nearly eliminates the oscillations, which are a strong impediment to a clean extrapolation to $1/L\to 0$. (b) Illustrating the cause of oscillatory behavior in ${\mathrm{\Delta }}_s$ is straightforward for $U=V=0$. Plotting the single-particle bands for $L=60$ (left) and 64 (right) ($N_{\mathrm{ch}}=6$—only four bands can be seen as two are doubly degenerate each) around the Fermi level (red line), one sees how the distance between highest occupied/lowest unoccupied levels (pairs marked with blue ellipsis) shifts with $L$, actually increasing here as $L$ increases. The oscillations vanish almost completely with the averaging over the $N_{\mathrm{ch}}$ bands (red boxes). [(c) and (d)] Scaling ground-state energies to zero truncation error ${\epsilon }_{\max }$, for the example of $L=40,N_{\mathrm{ch}}=4, N_{\uparrow }=N_{\downarrow }=40$. (c) Plotting $E_{\mathrm{GS}}$ against ${\epsilon }_{\max }$, computed for $\chi =12000,14000,16000,\text{and}18000$. (d) Plotting the resulting ${\overline{\mathrm{\Delta }}}_s[40,4]$ at $1/\chi =0$, and the corresponding values at the four finite values of $1/\chi$.

Figure 5

Scaling ${\overline{\mathrm{\Delta }}}_s[L,N_{\mathrm{ch}}]$ to $L\to \infty$ using the straight extrapolation protocol (cf. Sec. 4b). (a) Scaling at $V/t=0.5$, for $N_{\mathrm{ch}}=2$ (blue, crosses), 4 (red, squares), 6 (black, diamonds), and 8 (green, circles). (b) Scaling at $V/t=1$, symbols and colors for $N_{\mathrm{ch}}=2,4,$ and 6 as in (a). In both figures, ${\mathrm{\Delta }}_s[L,N_{\mathrm{ch}}=1]$ has been inserted for comparison purposes (grey stars). The raw data are summarized in Table 3 of Appendix pp1.

Figure 6

Plotting ${\overline{\mathrm{\Delta }}}_s[L\to \infty ,N_{\mathrm{ch}}]$ vs $N_{\mathrm{ch}}$ for (a) $V/t=0.5$. Red crosses denote outcomes of the straight extrapolation protocol, while grey lines with uncertainty ranges show outcomes of the extrapolation plus estimates protocol (cf. Sec. 4b and Table 1). (b) $V/t=1$. Symbols as in (a).

Figure 7

[(a) and (b)] Testing the viability of a MO/FL state in the 2D limit as an alternative explanation of the spin gap data, by evaluating the stability of the scaling Ansatz $B\left[N_{\mathrm{ch}}\right]/L+C\left[N_{\mathrm{ch}}\right]/L^2$ with $1/N_{\mathrm{ch}}$ for both $V/t=1$ (blue crosses) and $V/t=0.5$ (red crosses). The strong growth of $C\left[N_{\mathrm{ch}}\right]$ indicates the instability of the quadratic term, making the MO/FL hypothesis an unsuitable one, except potentially at $V/t=0.5$ (cf. text for details). [(c) and (d)] Testing the viability of a nodal USC state in the 2D limit as an alternative explanation of the spin gap data, by evaluating the stability of the scaling Ansatz $A\left[N_{\mathrm{ch}}\right]/L^{\alpha \left[N_{\mathrm{ch}}\right]}$ with $1/N_{\mathrm{ch}}$ [symbols matching (a) and (b)]. For $V/t=1$, there is a quick saturation to seemingly stationary values, indicating that nodal USC cannot definitely be ruled out as an alternative in the 2D limit (cf. text for details).

Figure 8

(a) Correlation functions for $N_{\mathrm{ch}}=8,V/t=0.5$, normalized to their value at $r=1$. Dash-dotted lines show relative decay at first maximum within $r\lesssim r_s$ for $d_{xy}$, (red), $d_{x^2-y^2}$ (green), $s$ (blue), and $C$ (black). This behavior of the normalized correlation functions in representative for all $N_{\mathrm{ch}}\ge 4$. (b) How correlation functions $d_{x^2-y^2},d_{xy}$ and $s\left(r\right)$, defined in Eqs. (6, 7, 8), are evaluated on the two central legs of the lattice.

Figure 9

Fermi surface of the noninteracting lattice electrons for the 2D U-V model (solid black lines) and possible symmetry of the order parameter, shown as a generic mixture between $d_{xy}$ and $g_{xy(x^2-y^2)}$ (see text for details). This symmetry minimizes, but does not eliminate intersections of nodal lines and Fermi surface (red circles), where the gap of USC order would exhibit zero nodes in standard mean-field treatments. From this, these theories predict a nongapped spectrum for spin excitation, i.e., zero spin gap.

Figure 10

[(a) and (b)] Decomposing ${\overline{\mathrm{\Delta }}}_s[L,N_{\mathrm{ch}}]$ into the linear part and oscillatory remnant, and ruling out a quadratic component to the linear part in the noninteracting system (example for $N_{\mathrm{ch}}=8$). (a) ${\overline{\mathrm{\Delta }}}_s[L,8]$ at $t_{\perp }/t=0.1$ and half filling, split into a seemingly linear component (blue crosses, linear fit as dash-dotted line) and weak oscillatory remnant (red crosses). (b) Confirming that the seemingly linear component in (a) is actually linear by plotting the $1/L^2$ coefficient of a second-order polynomial fitting as a function of the upper fitting range $L_{\max }$ (see text for discussion). [(c) and (d)] Instability and stability of fitting procedures respectively for the interacting systems, for example of ${\overline{\mathrm{\Delta }}}_s[L,4]$ at $U/t=4,V/t=1,t_{\perp }/t=0.1$ (magenta diamonds in both figures). (c) Fitting with a second-order polynomial for four largest ${\overline{\mathrm{\Delta }}}_s[L,4]$ values (red dash-dotted line), five largest values (green dashed line), and all six values (blue line). (d) Fitting with a first-order polynomial, line styles match the same fitting ranges as in (c).

Figure 11

Overview of pDMRG: (a) The physical system, e.g., a $M$-site 2D lattice, is mapped to the 1D DMRG chain. Every one of the possible $M+1$ bipartitions, carries left and right boundaries $L_n,R_n$ (cf. Sec. 3), typically by far the largest objects in any simulation that has sizable long-range terms in the Hamiltonian. (b) For every bipartition, the Hamiltonian consists of the (small) local MPOs $W_n,W_{n+1}$ on sites to the left/right of the bipartition, together with the boundaries $L_{n-1},R_{n+1}$. (c) The three layers of parallelization our pDMRG offers. In layer 1, each element $L_n\left[b\right]$ of $L_n$ is dispatched to a MPI group (each of which may be split into $G$ subgroups). Every $L\left[b_l\right]$ is a block-sparse matrix. In layer 2, every one of the dense matrices making up $L\left[b\right]$ is distributed round-robin over the $G$ split-off subgroups held by a global MPI group. In layer 3, the tiles making up a dense block can are distributed among the physical compute nodes held by each split subgroup.