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 the 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 the attractive V pp increases further, the system gives way to an unconventional d -wave superconductivity.

The Emery model, also known as the three-band Hubbard model, has long been proposed as one of the minimal models to understand the electronic properties of cuprate high-temperature superconductors Zaanen et al. (1985); Emery (1987); Scalettar (1989); Scalettar et al. (1991); White and Scalapino (2015); Huang et al. (2017); Jiang et al. (2023). In this model, a square lattice of copper (Cu) and oxygen (O) atoms in the CuO 2 plane (see Figure 1) is considered, where the Copper sites are represented by a single 3 d x 2 − y 2 orbital, while each oxygen site has one active 2 p orbital (2 p x or 2 p 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. n ̂ i σ d = d ̂ i σ + d ̂ i σ and n ̂ i σ p = p ̂ i σ + p ̂ i σ are the number operators for spin- σ at the Cu and O sites, respectively, with the total number operators are defined as n ̂ i d = ∑ σ n ̂ i σ d and n ̂ i p = ∑ σ n ̂ i σ p. Δ pd is the energy difference between having a hole on the Cu and oxygen sites. t p d i j and t p p i j are the hole hopping matrix elements between nearest-neighbor (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 frameworks for studying the cuprates superconductors, which captures phenomena like superconductivity, charge, and spin density wave orders Zaanen et al. (1985); Emery (1987); Scalettar (1989); White and Scalapino (2015); Huang et al. (2017); Jiang et al. (2023), it has more recently been extended to investigate the emergence of novel pair density wave (PDW) states Jiang (2023). In a PDW state, the superconducting (SC) order parameter carries finite center-of-mass momentum and varies spatially so that its spatial average vanishes Fulde and Ferrell (1964); Larkin and Ovchinnikov (1965); Berg et al. (2009); Fradkin et al. (2015); Agterberg et al. (2020); Lee (2014); Jian et al. (2020); Lozano et al. (2022). The PDW state has been considered a promising candidate state to understand the physics of cuprates high-temperature 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 Lee (2014); Fradkin et al. (2015); Agterberg et al. (2020); Himeda et al. (2002). 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 Hamidian et al. (2016); Ruan et al. (2018); Edkins et al. (2019); Liu et al. (2021) and La 1.875 Ba 0.125 CuO 4 Agterberg and Tsunetsugu (2008); Li et al. (2007); Berg et al. (2007); Tranquada et al. (2008); Tranquada (2020), Tranquada (2021), kagome superconductor CsV 3 Sb 5 Chen H. et al. (2021) and iron-based superconductor Liu et al. (2023).

The realization of a PDW state in microscopic lattice models remains 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 modulations Berg et al. (2010); Jaefari and Fradkin (2012); Venderley and Kim (2019); Xu et al. (2019); Peng et al. (2021a,b); Han et al. (2020); Huang et al. (2022); Wu et al. (2023); Jiang and Yao (2023). These include the Kondo-Heisenberg model Berg et al. (2010), the extended Hubbard-Heisenberg model Jaefari and Fradkin (2012), the strong coupling limit of the Holstein-Hubbard model Han et al. (2020); Huang et al. (2022) and generalized t - J and Hubbard models Xu et al. (2019); Venderley and Kim (2019); Peng et al. (2021a), Peng et al. (2021b). More recently, it has also been shown by one of the authors that the PDW ground state can also be realized in the three-band Hubbard model on a two-leg square cylinder Jiang (2023), where the SC correlations are dominant between neighboring Cu sites with d x 2 − y 2 -wave pairing symmetry.

In the present work, we consider the Emery model on the square lattice as defined in Figure 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) White (1992), White (1993). The signs of the hopping matrix elements in the related orbital configuration, i.e., Cu 3 d x 2 − y 2 orbital and O x / y 2 p x /2 p y orbitals, of an elementary plaquette centered at a generic Cu site is shown in Figure 1B. It is noted that the sign of t p p ≡ t p p σ − t p p π is taken to be negative, opposite to that is usually chosen Eskes et al. (1990). 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. White and Scalapino (2015); Jiang (2023), 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 White and Scalapino (2015); Armitage et al. (2010); Haule et al. (2014) 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 Figure 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 = 3 L x L y + 2 L y = 3 N u + 2 L y, where N u is the number of unit cells. 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 = 20,000 states with a typical truncation error ϵ ∼ 10 –10.

3.1 Phase diagram

Our main results are summarized in the ground state phase diagram in Figure 2. When V pp between oxygen sites are not strongly attractive, we find that the ground state of the system is consistent with that of a PDW state with power-law and mutually commensurate SC, charge-density wave (CDW), and spin-density-wave (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 White and Scalapino (2015); Jiang (2023), 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 Chen Z. et al. (2021); Qu et al. (2022); Peng et al. (2023), 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 increasing the attractive V pp, the system gives way to d -wave superconductivity where the PDW signatures are much suppressed. Interestingly, similar to the PDW phase, the Cooper pairing between adjacent oxygen sites also dominates in the pairing channel and the corresponding pairing symmetry is consistent with the d xy -wave. For even stronger attractive V pp, the ground state of the system becomes phase separated.

3.2 Pair density wave phase

As shown in Figure 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 Figure 3 and Figure 4 for several characteristic sets of parameters. Our conclusions hold for all parameters in the PDW phase in Figure 2.

3.2.1 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) = 1 2 [ c ̂ (x, y), ↑ † c ̂ (x, y) + α, ↓ † − c ̂ (x, y), ↓ † c ̂ (x, y) + α, ↑ † ] is spin-singlet pair creation operator on the bond α = a, b, d, d ̄, h and u defined in Figure 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 Jiang (2023), our results differ significantly for the negative t pp case. As illustrated in Figure 3A, we observe that the strongest SC correlations are prominently exhibited in Φ dd (r) and Φ d ̄ d ̄ (r) (not shown). 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: Φ d d (r) ∼ Φ d ̄ d ̄ (r) ∼ − Φ d d ̄ (r).

We have closely examined the spatial distribution of SC correlations, specifically targeting Φ dd (r). Our findings, based on three representative parameter choices, are depicted in Figure 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 power-law decay f (r) = A ∗ r − K s c Kühner et al. (2000). 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 Chen Z. et al. (2021); Qu et al. (2022); Peng et al. (2023) and the positive t pp Emery model Jiang (2023), 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 Figure 3C. More comprehensive results of K sc for Φ dd at δ = 1/8 are shown in Figure 3C. These point towards the divergence of the static PDW susceptibility, characterized as χ pdw ∼ T − (2 − K s c) 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 Figure 3D and align well with the fitting function ϕ dd (r) ∼ sin(Qr + θ). This pattern resonates with characteristics observed in the PDW state with a vanishing spatial average of ϕ (r) Agterberg et al. (2020). 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 Figure 3D, whereas λ sc ≈ 12 for δ = 1/12.

3.2.2 Charge density wave

We have calculated the charge density profile n α (x, y) = ⟨ n ̂ α (x, y) ⟩ and its rung average n (x) = ∑ y = 1 L y n a (x, y) / L y (e.g., Figure 4A) to describe the charge density properties of the system, where α = Cu/O x /O y site. Similar with the positive t pp case Jiang (2023), the spatial oscillation of n α (x) is also characterized by two ordering wavevectors at Q and 2 Q, corresponding to wavelengths λ Q ≈ 1/ δ and λ 2 Q ≈ 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 White et al. (2002). Here A Q and A 2 Q are amplitudes, ϕ 1 and ϕ 2 are phase shifts and n 0 is the mean density. Examples of the extracted exponents for δ = 1/8 are K c (Cu) ≈ 1.6 and K c (O y) ≈ 1.7 for V pd = 1 and V pp = 0.75, and K c (Cu) ≈ 1.7 and K c (O y) = 1.7 for V pd = V pp = 0.

3.2.3 Spin-spin correlations

To elucidate the magnetic properties of the ground state, we examine the spin-spin correlation functions F α (r) = ⟨ S ⃗ x 0, y 0 ⋅ S ⃗ x 0 + r, y 0 ⟩ where α = Cu/O x /O y site. Figure 4C illustrates examples of F (r), based on two representative parameter choices at δ = 1/8. Unlike the dominant spin-spin correlations between Cu sites Jiang (2023), our findings highlight a dominant correlation between oxygen sites in the kinetically frustrated case. Notably, this decays as a power-law, described by F (r) ∼ r − K s over extended distances. The associated Luttinger exponent is K s (O y) ≈ 1.1 for V pd = V pp = 0 and K s (O y) ≈ 1.3 for V pd = −0.8 and V pp = −0.2. In line with the characteristics of a PDW state, F (r) exhibits pronounced spatial oscillations (as seen in the inset of Figure 4C) with a characteristic wavelength, λ s = 1/ δ. This aligns closely with λ sc, yielding an ordering wavevector Q ≈ 2 πδ akin to the SC correlation.

3.2.4 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 derivable from the von Neumann entropy, formulated as: S (x) = −Tr ρ x ln ρ x 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 Calabrese and Cardy (2004); Fagotti and Calabrese (2011) that for an open system of length L x, where A ̃ 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 Figure 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.

3.3 d -wave superconductivity

When further increasing the attractive V pp, the system evolves into a d -wave SC phase. Similar to the PDW phase, we find that Φ dd (r) in Figure 5A and Φ d ̄ d ̄ (r) (not shown) exhibit the strongest SC correlations, i.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 Φ d d (r) ∼ Φ d ̄ d ̄ (r) ∼ − Φ d d ̄ (r).

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 Figure 5A, and the SC order parameter does not possess finite momentum, (2) the spin-spin 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 Figure 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 Emery (1987). This is reminiscent of the single-band Hubbard model on four-leg square cylinders as discussed in prior studies Jiang et al. (2018); Jiang and Devereaux (2019); Jiang Y.-F. et al. (2020); Chung et al. (2020); Jiang H.-C. et al. (2020); Gong et al. (2021); Peng et al. (2023).

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 Zhang and Rice (1988). 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 Lee et al. (2006); Fradkin et al. (2015).

To understand these results, we return to cluster calculations that determine the relevant parameters t and J of an effective single band model Eskes et al. (1989), Eskes et al. (1990). 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 1. 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 Figure 3.

Calculations for three holes on the same cluster yield the hopping parameter t, defined as the energy difference between the ground and first excited state. For V pp = V pd = 0, 2 t = −0.673 for t pp = 0.5 while 2 t = −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, the binding of the holes can be examined for the ground state of four holes on the same cluster (see Figure 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.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 Figure 2 may be related to this tendency to bind neighboring O holes.