N ov 2 01 9 Cooperative enhancement of superconducting correlations by electron-electron and electron-phonon interactions in the quarter-filled band

We present the results of Quantum Monte Carlo calculations for a two dimensional frustrated Hubbard model coupled to bond phonons. The model is known to have a d-wave superconducting ground state in the limit of large phonon frequency for sufficiently strong electron-phonon coupling. In the absence of electron-phonon coupling the Hubbard interaction U enhances superconducting pairing in the quarter-filled (density ρ = 0.5) band. We show here that at ρ = 0.5 electron-electron and electron-phonon interactions cooperatively reinforce d-wave pairing, while competing with each other at all other densities. Cooperative degrees of freedom are found in many phase transitions and are essential to understanding superconductivity in strongly correlated materials.

Introduction.-There is strong evidence that any successful theory of superconductivity (SC) in the high T c cuprates and other unconventional superconductors must incorporate the effects of electron-electron (e-e) interactions. It is now often assumed that superconducting pairing in these materials can be mediated solely by e-e interactions, leading to numerous studies of model Hamiltonians with e-e interactions such as the Hubbard model. While the properties of the two dimensional (2D) Hubbard model remain under active debate, calculations with the best unbiased methods available have found no longrange superconducting pairing in the weakly doped Hubbard model on a square lattice [1,2] or the half-filled Hubbard model on the anisotropic triangular lattice [3].
One possible reason that such calculations have failed to find SC is that the superconducting state is of a more complex form than conventionally assumed, for example a pair density wave [4][5][6][7][8][9]. It is also possible that repulsive e-e terms alone are not sufficient for SC and that other interactions must be considered. The most obvious candidate is the electron-phonon (e-p) interaction. In the cuprates there is strong experimental evidence for the importance of phonons in the electronic properties and SC. This includes giant softening of the Cu-O bond stretching frequency in the underdoped cuprates [10], the oxygen isotope effect [11,12], and kinks observed in photoemission experiments [13]. Phonons also appear to play some role in the ubiquitous charge ordering found across the cuprate family [10,14].
Away from the weakly doped half-filled band there have been fewer investigations of SC in interacting models. Calculations on frustrated lattices have shown that e-e interactions enhance superconducting pairing at quarter-filling (corresponding to a carrier density of ρ = 0.5) while suppressing pairing at all other densities [15][16][17]. Clay and collaborators have proposed a valencebond theory of SC at ρ = 0.5 [18]. In this picture SC emerges from a proximate spin-gapped insulating state, similar to the proposal of resonating valence-bond (RVB) SC emerging from a 2D valence-bond solid (VBS) [19]. The RVB theory of SC in the weakly doped Hubbard model remains controversial, and as noted above calculations for ρ 1 have not found SC. In two dimensions antiferromagnetism (AFM) is preferred over singlet formation and calculations at ρ = 1 have also failed to find a VBS state [20].
A valence bond insulator, the Paired Electron Crystal (PEC), does occur at ρ = 0.5 [18,21,22]. The PEC is a density wave of pairs, with nearest-neighbor (n.n.) singlet pairs separated by pairs of vacant sites. The PEC is characterized by a spin gap and coexisting charge and bond order, and becomes favored over AFM at ρ = 0.5 in the presence of lattice frustration [21,22]. The PEC is often adjacent to SC in the organic charge-transfer solids (CTS), which are all ρ = 0.5 materials [18]. Besides the CTS many other ρ = 0.5 superconductors are known [18], and this mechanism of SC has been further generalized to the cuprates [23]. E-p interactions are required to realize the PEC because of the bond distortion that simultaneously occurs with n.n. singlet formation at ρ = 0.5 [21,22]. It is natural to expect that SC evolving from a destabilized PEC must involve both e-e and e-p interactions.
Model.-In many theories e-p and e-e interactions have been considered to be mutually exclusive in mediating SC. This belief is partly due to the failures of the phononbased pairing mechanism of the BCS model in unconventional superconductors [24], but is also reinforced by the study of Hamiltonians that combine e-e and e-p interactions. The simplest is the Hubbard-Holstein model (HHM), where at each site a dispersionless phonon is coupled to the charge density: In Eq. 1 c † i,σ creates an electron of spin σ on site i, n i,σ = c † i,σ c i,σ , t ij is the electron hopping integral, and U is the onsite Hubbard interaction. x i and p i are coordinate and momentum operators for the phonon oscillator at site i with mass M and frequency ω. The e-p coupling constant is g. We give energies below in units of the bare hopping t. The physics of the HHM is governed by the competition between onsite pairing driven by the e-p interaction and opposed by U . This is most clearly seen in the limit ω → ∞ where the effective Hubbard interaction between electrons is U eff = U − 2g 2 /ω.
We consider here a model with dispersionless bondcoupled (Su-Schrieffer-Heeger-type [25]) phonons, In Eq. 2 x (ij) is the phonon coordinate associated with the bond connecting sites i and j and all the other terms have identical meaning to Eq. 1. Compared to the large amount of work on the HHM there are few numerical studies of bond-coupled phonons beyond the classical limit. Most work at finite ω has been in one dimension [26][27][28][29][30]; a 2D multi-orbital lattice was studied recently, although only for U = 0. [31]. The ω → ∞ limit of Eq. 2 on a square lattice results in a more complex effective interaction than in the HHM, with four types of terms [30,32,33]. These include a nearest neighbor repulsion for parallel spins, on-site pair terms (which will be suppressed by U ), and pair hopping of nearest-neighbor singlet pairs. Recent work has shown that this pair hopping interaction leads to strongly bound bipolarons with light effective masses that are stable against strong e-e interactions [30]. This suggests the possibility of e-p mediated SC in the presence of strong e-e interactions at finite carrier densities. The effective ω → ∞ Hamiltonian, H W , was studied using quantum Monte Carlo near half-filling for U = 4 on the square lattice [33][34][35]. In this density region the pair hopping terms d x 2 −y 2 mediate pairing but compete with the Mott insulating state driven by U [33][34][35]. While this shows that d-wave SC can in principle result from e-p interactions, for ρ ≈ 1 d-wave SC is only possible provided the coupling is strong enough to overcome e-e interactions, again reinforcing the belief that e-e and e-p interactions cannot simultaneously drive SC. In this Letter we investigate the combined effects of bond-phonon coupling and e-e interactions on superconducting pairpair correlations over a wide density range including both ρ = 1 and ρ = 0.5. We find that at ρ = 0.5 e-e and e-p interactions cooperatively enhance pairing.
Numerical results.-We define singlet pair creation operators ∆ i , where g(ν) is a relative sign that determines the pairing symmetry. The pair-pair correlation function is P (r) = ∆ † i ∆ i+ r . A theory of correlated electron SC should satisfy two requirements, (i) at zero temperature P (r) must have long-range order, and (ii) P (r) in the presence of e-e interactions should be enhanced over its value for non-interacting fermions. Finite and zero temperature calculations on frustrated Hubbard models of up to 128 sites have shown that U enhances d-wave pairing preferentially at ρ = 0.5 while suppressing pair-pair correlations relative to their U = 0 value at all other ρ [15][16][17]. While enhanced by U the magnitude of P (r) however decreases with distance, consistent with either a zero or possibly small long-range superconducting order parameter [15][16][17].
To solve Eq. 2 we use finite-temperature determinant quantum Monte Carlo (DQMC) which provides unbiased results [36]. The lowest temperatures that can be reached in DQMC are strongly limited by the Fermion sign problem. In the density region of most interest here, ρ ≈ 0.5, the sign problem is considerably less severe than for the more intensively studied ρ ≈ 0.8. Nevertheless, we are not able to reach low enough temperatures to determine whether the ground state has long-range superconducting order; however it is important to recall that the effective model in the ω → ∞ limit does have a superconducting ground state [33][34][35]. Besides the sign problem, Monte Carlo autocorrelation times often increase exponentially with e-p coupling near phase transitions [37]. To help mitigate this we implemented the block phonon updates of Reference [38]. We use an imaginary time discretization of ∆τ = 0.1, which is small enough that this source of systematic error can be neglected. We report results in terms of the dimensionless e-p coupling strength λ = α 2 t x /(M ω 2 ) and set M = 1. Further information on the method is given in the Supplemental information [39].  We performed calculations on 4×4, 6×6, and 10×10 anisotropic triangular lattices [15]. This lattice has a single frustrating bond t ′ across each plaquette; in the limit t ′ = t (t ′ = 0) it is the triangular (square) lattice. Because the PEC requires lattice frustration [21,22], we work in the strongly frustrated limit with t ′ = 0.8. The lattice dimensions are chosen under the constraint that the ρ = 0.5 single-particle state is non-degenerate [15]. We also take t y = 0.9 slightly different from t x = 1.0 in order to increase the number of non-degenerate densities. This lessens the severity of the Fermion sign problem and makes calculations feasible over a range of ρ. The precise choice of t x , t y , and t ′ are however not critical to the results we report. P (r = 0) can be decomposed into combinations of charge and spin correlations [40]. AFM order leads to a trivial increase of the short-range component of P (r), even as the long-range component is strongly suppressed [3,40]. To mitigate such finite-size effects we exclude small r correlations and measure the average long-range value of P (r) [41], In Eq. 4 only correlations over distances greater than two lattice spacings are kept; N p is the number of such terms. We considered four pairing symmetries: s-wave involving n.n. pairs, s xy with next-nearest-neighbor (n.n.n) pairs, d x 2 −y 2 , and d xy .
For the 4×4 and 6×6 lattices we calculatedP over the density range 0.2 ρ ≤ 1.0. In Fig. 1 we plotP for s pairing versus density for the 6×6 lattice (similar data for the 4×4 lattice is in the Supplemental information [39]). We find that at all densities U suppresses s pairing; as U increasesP becomes closer to zero across the entire density range. There is an increase inP for s pairing over its value for non-interacting electrons for strong λ in the low density region ρ 0.4. An s or s xy SC state may exist in the model in the very low density region provided U is not too large; we will not consider this parameter region further here. s xy pairing shows some enhancement by λ for ρ > 0.4, but is also suppressed by U [39].
In the thermodynamic limit on the anisotropic triangular lattice we expect a pairing symmetry that mixes d xy and d x 2 −y 2 ; on finite lattices either d x 2 −y 2 or d xy is favored [15,42]. In Fig. 2 we plotP for the d xy symmetry versus density for the 6×6 lattice. Fig. 2(a) shows the effect of increasing U at fixed e-p coupling strength. Compared to the noninteracting system (solid line), pairing correlations are enhanced by U selectively at ρ = 0.5. At all other densities pairing is suppressed by U . Fig. 2(b) shows the effect of increasing e-p coupling at fixed U . Here phonons enhance the pairing over a wide density range for ρ 0.4, including at both ρ = 0.5 and ρ = 1. However, only at ρ = 0.5 do e-e and e-p interactions both enhanceP ; at other densities the interactions compete.
For the 10×10 lattice we find that the Fermion sign is reasonable for both ρ ≈ 0.5 and ρ = 1. Fig. 3 summarizes our results for ρ = 1 on all of the lattices. At ρ = 1P for either d x 2 −y 2 or d xy pairing increases with λ for λ 1. This is consistent with the d-wave superconducting state found in References [33][34][35] in the ω → ∞ limit of the model. However, U competes with the e-p interaction at ρ = 1 withP decreasing with increasing U . We expect on less frustrated lattices at ρ = 1 a competition between dwave SC mediated by the bond phonons and AFM order mediated by U . On the 4×4 lattice at large λP for s xy pairing becomes comparable toP for d x 2 −y 2 pairing [39]. This may be the reason for the weak decrease of the d x 2 −y 2P with U at large λ in Fig. 3(a).
The behavior ofP at ρ = 0.5 (Fig. 4) is very different from ρ = 1. At ρ ∼ 0.5 zero-temperature calculations find thatP increases with U [15][16][17]. Here, for all the lattices we find thatP increases with increasing U and λ at the same densities where T = 0 calculations find enhancement by U alone. Cooperative interactions however, should not merely both increase the value of an order parameter, but the effect of the first interaction should be strengthened in the presence of the second, and vice versa. This is indeed what we see at ρ = 0.5: First, we see that λ enhances the effect of U : in Fig. 4 the increase inP between U = 0 and U = 3 for all lattices is larger for λ > 0 than at λ = 0. In fact, at the temperatures we can access here, on some latticesP decreases with U at λ = 0. Second, nonzero U enhances the increase inP with λ. As a function of λ,P reaches a broad maximum at λ = λ max . Comparing the value of P at λ = 0 and λ max , there is a larger increase for U > 0 compared to U = 0. These data show that at ρ ≈ 0.5 U and SSH phonon interactions not only both enhance pairing, but their effect is cooperative. This is the central result of our work.
In all lattices we find that as a function of λ,P reaches a broad peak at an intermediate value of λ before beginning to decrease. A similar broad maximum inP is seen as a function of U in zero temperature λ = 0 calculations at ρ = 0.5 [16,17]. We expect the decrease inP at larger λ and/or U is caused by the increasing effective mass of pairs. λ max and the amount of enhancement also depend on ω. In Fig. 5 we show the effect of the phonon frequency ω onP for the 6×6 lattice. As ω decreases, λ max shifts towards stronger coupling and the amount of enhancement increases. This shows that the phonon dispersion relation will play an important role in understanding a superconducting state mediated by the combination of e-e and e-p interactions.
Conclusion.-We have shown that at ρ = 0.5 bondcoupled phonons act cooperatively with onsite Coulomb interactions in enhancing superconducting pairing. Both interactions promote a superconducting state through their effect on short-range (n.n.) singlet pairing; the Hubbard U through antiferromagnetic superexchange, and bond-coupled phonons through effective pair-hopping. These interactions act cooperatively only for density ρ ≈ 0.5 while for all other densities they compete. The presence of cooperative degrees of freedom is essential to understand phase transitions in real materials, for example the non-superconducting phases of the organic CTS [18]. Producing an unconventional superconducting state from e-p interactions alone in the presence of competing e-e interactions would require a careful tuning of parameters that is unrealistic. As noted above, model calculations have suggested that e-e interactions are not sufficient for SC. Cooperative e-e and e-p interactions resolve both of these problems.
The observation of cooperative enhancement of pairing only at ρ ≈ 0.5 supports our proposal of SC emerging from the PEC [18]. The PEC has period four charge and bond order lattice [21,22]. This is only commensurate on the 4×4 lattice, where with classical phonons a metal-PEC transition occurs at a critical e-p coupling strength [21,22]. On the 4×4 lattice we find the peak bond-bond correlation for ρ = 0.5 is at Q = (π/2, π), which is consistent with the PEC [21,22]. However, with the λ and temperatures accessible to DQMC we are not able to reach the metal-PEC transition. Intersite Coulomb interactions V ij (below the critical value V c for Wigner crystal formation) strengthen the PEC and might be necessary to see the metal-PEC transition on large lattices [22]. Further calculations over a range of t ′ and on commensurate lattices are in progress. Studies of the effective interaction in the ω → ∞ limit across a wider density range will also be useful.
We thank S. Mazumdar for useful discussions. Some calculations in this work used the Extreme Science and Engineering Discovery Environment [43] (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. Specifically, we used the Bridges system [44] which is supported by NSF award number ACI-1445606, at the Pittsburgh Supercomputing Center under awards TG-DMR190052 and TG-DMR190068. We also used the Comet system of the San Diego Supercomputer Center and the Stampede2 system of the Texas Advanced Computing Center under award TG-DMR190052.