Pair density wave and superconductivity in a kinetically frustrated doped Emery model on a square lattice

The quest to understand the nature of superconductivity in cuprates has spotlighted the pair density wave (PDW) -- a superconducting state characterized by a spatially modulated order parameter. Despite significant advances in understanding PDW properties, conclusively demonstrating its presence in systems pertinent to cuprate superconductors remains elusive. In this study, we present a systematic density-matrix renormalization group study to investigate the Emery model (or the three-band Hubbard model) on two-leg square cylinders with negative electron hopping term $t_{pp}$ between adjacent oxygen sites. Kinetic frustration - introduced by changing the sign of oxygen-oxygen hopping - leads to a much reduced Cu-Cu antiferromagnetic exchange along with an enlarged charge transfer energy that changes the local properties of the model. At light doping levels, our findings reveal a ground state remarkably consistent with a PDW, exhibiting mutually commensurate superconducting (SC), charge and spin density wave correlations. Intriguingly, the dominant SC pairing is observed between neighboring oxygen sites, diverging from the expected Cu sites in the positive $t_{pp}$ case. When the system incorporates moderate near-neighbor interactions, particularly an attractive $V_{pp}$ between adjacent oxygen sites, the SC correlations become quasi-long-ranged, accompanied by a pronounced divergence in the PDW susceptibility. When further increase the attractive $V_{pp}$, the system gives ways to an unconventional $d$-wave superconductivity.

Introduction -The Emery model, also knows as the three-band Hubbard model, has long been proposed to as one of the minimal models to understand the electronic properties of cuprate high-temperature superconductors [1][2][3][4][5][6][7].In this model, a square lattice of copper (Cu) and oxygen (O) atoms in the CuO 2 plane (see Fig. 1) is considered, where the Copper sites are represented by a single 3d x 2 −y 2 orbital, while each oxygen site has two active 2p orbitals (2p x and 2p y ).In the hole representation, the model Hamiltonian is defined as Here d+ iσ and p+ jσ create holes with spin-σ on the i th Cu and j th oxygen sites, and ⟨ij⟩ denotes NN sites.nd iσ = d+ iσ diσ and np iσ = p+ iσ piσ are the number operators for spin-σ at the Cu and O sites, respectively, with the total number operators are defined as nd i = σ nd iσ and np i = σ np iσ .∆ pd is the energy difference between having a hole on the Cu and oxygen sites.t ij pd and t ij pp are the hole hopping matrix elements between nearestneighbor (NN) Cu and oxygen sites and the NN oxygen sites, respectively.U d and U p are the on-site Cu and oxygen Coulomb repulsion, and V pd and V pp are the NN Cu-O and O-O Coulomb interactions, respectively.While the Emery model has been proposed as one of the critical framework for studying the cuprates superconductors, which captures phenomena like superconductivity, charge and spin density wave orders [1][2][3][5][6][7], it has more recently been extended to investigate the emergence of novel pair density wave (PDW) states [8].In a PDW state, the superconducting (SC) order parameter carries finite center-of-mass momentum and varies spatially so that its spatial average vanishes [9][10][11][12][13][14][15][16].The PDW state has been considered as a promising candidate state to understand the physics of cuprates hightemperature superconductors and other strongly correlated systems, where it has been proposed that various phases, including the superconductivity, charge and spin density wave orders, can emerge by partially melting the PDW state [12][13][14]17].Recently, intense interest in the PDW state has emerged due to experimental observations in cuprate superconductors Bi 2 Sr 2 CaCu 2 O 8+x [18][19][20][21] and La 1.875 Ba 0.125 CuO 4 [22][23][24][25][26][27], kagome superconductor CsV 3 Sb 5 [28] and ion-based superconductor [29].
The realization of PDW state in microscopic lattice models remain highly nontrivial and usually involves modifying or extending existing frameworks like the Hubbard models to include competing interactions or inhomogeneities that can give rise to spatial modulation [30][31][32][33][34][35][36][37][38][39].These include the Kondo-Heisenberg model [30], the extended Hubbard-Heisenberg model [31], the strong coupling limit of the Holstein-Hubbard model [36,37] and generalized t-J and Hubbard models [32][33][34][35].More recently, it has also been shown by one of us that the PDW ground state can also be realized in the three-band Hubbard model on a two-leg square cylinder [8], where the SC correlations are dominant between neighboring Cu sites Model and method -In the present work, we consider the Emery model on the square lattice as defined in Fig. 1 and Eq.( 2) to study whether the same PDW state or distinct SC state emerges upon doping as well as the associated pairing symmetry using the density-matrix renormalization group (DMRG) [40,41].The signs of the hopping matrix elements in the related orbital configuration, i.e., Cu 3d x 2 −y 2 orbital and O x/y 2p x /2p y orbitals, of an elementary plaquette centered at a generic Cu site is shown in Fig. 1B.It is noted that the sign of t pp ≡ t σ pp −t π pp is taken to be negative, opposite to that is usually chosen [42].The resulting increase in the delocalization energy involving ligand L oxygen orbitals raises the level of the effective charge transfer energy, alters the magnetic exchange among Cu and O, and affects the local symmetry of the ground state orbital configuration.This has a profound impact on the local physics and ground state properties of the system.Following Ref. [5,8], we set t pd = 1 as the energy unit and take a canonical set of parameters U d = 8, U p = 3, ∆ pd = 3 for cuprates [5,43,44] but negative t pp = −0.5, and study the ground state properties of this system as a function of V pd and V pp .We focus on two-leg cylinders as shown in Fig. 1 with width L y = 2 and length up to L x =96, where L x and L y are the number of unit cells along the e 1 and e 2 directions, respectively.The total number of sites is N = 3L x L y + 2L y = 3N u + 2L y , where N u is the number of unit cell.The overall hole density of the system is defined as ρ = 1 + δ, where δ = N h /N u and N h denote the hole doping concentration and number of doped holes away from half-filling, respectively.We consider δ = 1/12 and 1/8, and keep up to m = 20000 states with a typical truncation error ϵ ∼ 10 −10 .
Phase diagram -Our main results are summarized in the ground state phase diagram in Fig. 2. When V pp between oxygen sites are not strongly attractive, we find the the ground state of the system is consistent with that of a PDW state with power-law and mutually commen- surate SC, charge-density wave (CDW) and spin-densitywave (SDW) correlations.The SC correlations oscillate periodically in real space in such a way that its spatial average vanishes and the PDW ordering wavevector Q ≈ 2πδ is incommensurate.Contrary to the positive t pp case [5,8], our results show that the SC pairing is dominant between adjacent oxygen sites instead of Cu sites.Accordingly, the SC pairing symmetry is consistent with d xy rather than d x 2 −y 2 wave.Similar to the single-band Hubbard model on the square lattice [45][46][47], the finite electronic attractions V pd and V pp , especially V pp between oxygen sites, can notably enhance the SC correlations while simultaneously suppress the CDW correlations.For modestly strong V pp interaction, including both repulsion and attraction, the SC correlations become strong enough so that a quasi-long-range PDW order emerges with K sc < 2 and divergent static PDW susceptibility.When further increase the attractive V pp , the system gives ways to the d-wave superconductivity where the PDW signatures get much suppressed.Interestingly, similar with the PDW phase, the Cooper pairing between the adjacent oxygen sites also dominates in the pairing channel and the corresponding pairing symmetry is consistent with the d xy -wave as well.For even stronger attractive V pp , the ground state of the system becomes phase separated.
Pair density wave phase -As shown in Fig. 2, the majority of the ground state phase diagram is occupied by the PDW phase,where the SC correlations decay as a power-law at long distances and oscillate periodically in real space in such a way that its spatial average vanishes.We provide detailed examples in Fig. 3 and Fig. 4 for several characteristic sets of parameters.Our conclusions hold for all parameters in the PDW phase in Fig. 2.
Superconducting correlations -In order to explore the potential for superconductivity, we have calculated the equal-time spin-singlet SC pair-pair correlations defined as Here, ∆ † ĉ † (x,y)+α,↑ ] is spin-singlet pair creation operator on the bond α = a, b, d, d, h and u defined in Fig. 1A.(x 0 , y 0 ) is a reference bond with x 0 ∼ L x /4, r is the distance between two bonds in the e 1 direction.We have comprehensively analyzed the various components of the SC correlations.This includes calculations of Φ aa , Φ ab , Φ bb , Φ dd (r), Φ d d(r), Φ d d(r), Φ hh , Φ uu and Φ uh .While the positive t pp case primarily showcases dominant correlations in Φ hh and Φ uu as discussed in [8], our results differ significantly for the negative t pp case.As illustrated in Fig. 3A, we observe that the strongest SC correlations are prominently exhibited in Φ dd (r) and Φ d d(r).This suggests that the pairing is more dominant between neighboring oxygen sites, rather than Cu sites.Furthermore, even though the pairing symmetry aligns with the d-wave, our findings indicate a shift to the d xy -wave symmetry in this particular scenario, diverging from the previously understood d x 2 −y 2 -wave.This distinction in d xy -wave symmetry is characterized by the relationship: We've closely examined the spatial distribution of SC correlations, specifically targeting Φ dd (r).Our findings, based on three representative parameter choices, are depicted in Fig. 3D.Here, Φ dd (r) exhibits clear spatial oscillations as Φ dd (r) ∼ f (r)ϕ dd (r) over a vast region of r.In this context, f (r) acts as the envelope, while ϕ dd (r) gives rise to the spatial oscillation.As we move to longer distances, the envelope function f (r) adheres to a powerlaw decay f (r) = A * r −Ksc [48].For instance, we derived an exponent K sc ≈ 2.4 at V pd = 1.0 and V pp = 0.75, and K sc ≈ 1.9 at V pd = 0 and V pp = 0, and K sc ≈ 1.7 at V pd = −0.8 and V pp = −0.2.Drawing connections with the established single-band Hubbard model [45][46][47] and the positive t pp Emery model [8], we find that diminishing the NN repulsion or amplifying the NN attraction, especially V pp , can notably enhance SC correlations.This observation is further validated by the K sc values presented in Fig. 3C.More comprehensive results of K sc for Φ dd at δ = 1/8 are shown in Fig. 3C.These point towards the divergence of the static PDW susceptibility, characterized as χ pdw ∼ T −(2−Ksc) when T → 0. We have also calculated the spin-triplet SC correlations, which are markedly weaker than their spin-singlet counterparts.
The spatial oscillation of the SC correlations Φ(r) is captured by the normalized function ϕ(r), as previously defined.Depictions of ϕ dd (r) are presented in Fig. 3D and align well with the fitting function ϕ dd (r) ∼ sin(Qr + θ).This pattern resonates with characteristics observed in the PDW state with vanishing spatial average of ϕ(r) [13].The PDW ordering wavevector appears to be incommensurate as Q ≈ 2πδ with a corresponding wavelength λ sc ≈ 1/δ.For instance, λ sc ≈ 8 for δ = 1/8 as evidenced in Fig. 3D, whereas λ sc ≈ 12 for δ = 1/12.
Charge density wave -We have calculated the charge density profile n α (x, y) = ⟨n α (x, y)⟩ and its rung average n(x) = Ly y=1 n a (x, y)/L y (e.g., Fig. 4A) to describe the charge density properties of the system, where α=Cu/O x /O y site.Similar with the positive t pp case [8], the spatial oscillation of n α (x) is also characterized by two ordering wavevectors at Q and 2Q, corresponding to wavelengths λ Q ≈ 1/δ and λ 2Q ≈ 1/2δ, respectively.
At long distance, the spatial decay of the CDW correlation is dominated by a power-law with an exponent K c , which can be obtained by fitting the charge density oscillations induced by the cylinder boundaries [49] (3) Here A Q and A 2Q are amplitudes, ϕ 1 and ϕ 2 are phase shifts and n 0 is the mean density.Examples of the ex-  [8], our findings highlight a dominant correlation between oxygen sites in the kinetically frustrated case.Notably, this decays as a powerlaw, described by F (r) ∼ r −Ks over extended distances.The associated Luttinger exponent is In line with the characteristics of a PDW state, F (r) exhibits pronounced spatial oscillations (as seen in the inset of Fig. 4C) with a characteristic wavelength, λ s = 1/δ.This aligns closely with λ sc , yielding an ordering wavevector Q ≈ 2πδ akin to the SC correlation.
Entanglement entropy -Our findings indicate the presence of multiple gapless modes, encompassing both charge and spin degrees of freedom.These can be characterized by the central charge, c.This charge is deriv- able from the von Neumann entropy, formulated as: S(x) = −Trρ x lnρ x where where ρ x represents the reduced density matrix for a subsystem of length x.For critical systems in 1+1 dimensions, described through a conformal field theory, it has been established [50,51] that for an open system of length L x , where Ã and S are model dependent fitting parameters, and k F is the Fermi momentum.Our findings reveal a central charge approximately given by c ≈ 2, with illustrative examples provided in Fig. 4D.Specifically, we observed c ≈ 1.95 at V pd = V pp = 0 and c ≈ 2.0 at V pd = −0.8 and V pp = −0.2 at δ = 1/8.These results point to the presence of both a gapless charge mode and a gapless spin mode.d-wave superconductivity -When further increase the attractive V pp , the system envolves into a d-wave SC phase.Similar with the PDW phase, we find that Φ dd (r) and Φ d d(r) exhibit the strongest SC correlations as shown in Fig. 5A, i While there are similarities, several significant distinctions can be drawn between the PDW phase and the d-wave SC phase: (1) in the d-wave SC phase, the SC correlations Φ αβ (r) maintain a consistent sign in real space, as depicted in Fig. 5A, and the SC order parameter does not possess finite momentum, (2) the spinspin correlation functions in this phase are short-ranged and undergo exponential decay, (3) a singular gapless mode with c ≈ 1 is evident, as illustrated in Fig. 5B.
For instance, the extracted central charge c ≈ 1.04 for V pd = −0.8 and V pp = −0.8, and c ≈ 1.2 for V pd = −0.8 and V pp = −0.4.Given these observations, our results affirm that the ground state of the d-wave SC phase aligns with the characteristics of a Luther-Emery liquid [2].This is reminiscent of the single-band Hubbard model on four-leg square cylinders as discussed in prior studies [47,[52][53][54][55][56][57].
Summary and discussion -In conclusion, we have extensively investigated the ground state properties of the lightly doped three-band Hubbard model on two-leg square cylinders, specifically focusing on near-neighbor Cu-O and O-O interactions.Our results strongly suggest that the system's ground state is aligned with the characteristics of a PDW state, showcasing quasi-long-range PDW order and pronounced susceptibility.Several aspects of our findings are unexpected.Within the doped negative t pp Emery model, Cooper pairing prominently emerges between adjacent oxygen sites rather than between neighboring Cu sites.This stands in stark contrast to the prevailing understanding, where Cooper pairing is believed to be dominant between neighboring Cu sites.It has been postulated that cuprate physics can be encapsulated by a single-band effective Hamiltonian, exclusively encompassing the Cu d holes [58].While the pairing symmetry aligns with the d-wave symmetry, it manifests as d xy rather than d x 2 −y 2 .This diverges from the d x 2 −y 2 pairing symmetry characteristic of cuprates [12,59].To understand these results, we return to cluster calculations that determine the relevant parameters t and J of an effective single band model [42,60].Specifically we consider Cu 2 O 7 clusters with the same parameters used for DMRG, and compare the results for positive and negative t pp .Our results are summarized in Table I.Determining the singlet-triplet energy difference for two holes on the cluster yields an exchange energy J = 0.0179 for t pp = −0.5, compared to 0.165 if the sign of t pp is reversed.Largely interpreted as due to an increase in the effective charge transfer energy when ligand delocalization is considered, the effective spin exchange between Cu spins is greatly reduced for negative t pp .While the dependence on V pp is negligible, a negative V pd increases the magnetic exchange for the two-hole ground state configuration.The increase of J for negative V pd may help to favor stronger hole singlet bonding, promoting stronger SC susceptibilities in Fig. 3. Calculations for three holes on the same cluster yields the hopping parameter t, defined as the energy difference between the ground and first excited state.For V pp = V pd = 0, 2t = −0.673for t pp = 0.5 while 2t = −8 × 10 −4 for t pp = −0.5.These numbers increase slightly for V pp < 0, but overall a reversed sign of t pp dramatically affects the magnitude of the NN hopping.
Lastly, binding of doped holes can be examined for the ground state of four doped holes on the same cluster (see Fig. 6).For both t pp positive and negative, the ground state is a spin-singlet, with Cu spins and O spins both forming singlets.However, the spatial orientation of bound holes on O is different: for t pp = 0.5, O holes primarily bind to Cu at the ends of the cluster, without substantial hole occupation on the central, bridging oxygen, while for t pp = −0.

FIG. 1 :
FIG. 1: Emery model on the square lattice.(A) The squares represent Cu d x 2 −y 2 orbitals and circles represent O x/y 2px/2py orbitals.Periodic (open) boundary condition is imposed in the direction specified by the lattice basis vector e2 = (0, 1) (e1 = (1, 0)).The dashed loops represent bonds a, b, d, d, h and u. (B) Signs of hopping matrix elements in an elementary plaquette in the CuO2 plane.

FIG. 2 :
FIG. 2: Ground state phase diagram of the Emery model on two-leg square cylinders at δ = 1/8.The solid symbols are numerical data points and PS denotes phase separation.The shaded regions are guides for eyes.

FIG. 3 :
FIG. 3: Superconducting correlations in the pair density wave phase at δ = 1/8.The magnitude of SC correlations are shown in (A) for Φ(r) at V pd = Vpp = 0, and in (B) for Φ dd (r) with different V pd and Vpp.The dashed lines represent fits to a power-law function f (r) ∼ r −Ksc .Data points far from the envelope and those at short distances are discarded in gray color in the fitting process.(C) Luttinger exponent Ksc as a function of V pd and Vpp.(D) The normalized functions ϕ dd (r) = Φ dd (r)/f dd (r) reflect the spatial oscillation of Φ dd (r).

FIG. 4 :
FIG. 4: Charge density profile, spin-spin correlation and entanglement entropy in the pair density wave phase at δ = 1/8.(A) Charge density profiles nCu(x) on the Cu site and nOy(x) on the Oy site for V pd = Vpp = 0. (B) The magnitude of the spin-spin correlation |F (r)| where the dashed lines represent power-law fits f (r) ∼ r −Ks .(C) The normalized function F (r)/f (r) reflects the spatial oscillation of F (r) in (B).(D) Von Neumann entanglement entropy S(x).Note that a few data points in gray color close to the open ends are excluded to minimize boundary effects.

FIG. 5 :
FIG. 5: Superconducting correlations and entanglement entropy in the d-wave SC phase at δ = 1/8.(A) SC correlations for V pd = −0.8 and Vpp = −0.8.Dashed lines represent fits to a power-law function f (r) ∼ r −Ksc .(B) Von Neumann entanglement entropy S(x).Note that a few data points in gray color close to open ends are excluded to minimize the boundary effect.
.e. the pairing is dominant between neighboring oxygen sites instead of Cu sites.The pairing symmetry is also consistent with the d xy -wave symmetry characterized by the fact Φ dd (r) ∼ Φ d d(r) ∼ −Φ d d(r).

FIG. 6 :
FIG.6: Predominant hole distributions for tpp positive or negative, as indicated.Solid red spheres denote approximately a full hole charge, while shaded red spheres denote an approximate quarter charge.

5 ,
O holes primarily bind to the central oxygen in a local d xy configuration with neighboring oxygens.As this might be expected when antiferromagnetic exchange among oxygen becomes dominant over Cu, the predominance of d xy pairing observed in the phase diagram of Fig. 2 may be related to this tendency to bind neighboring O holes.

TABLE I :
Effective single band exchange parameter J and NN hopping t determined from Cu2O7 clusters for different values of tpp as indicated.All other parameters are the same as used in the main text.