Orbital structure of the effective pairing interaction in the high-temperature superconducting cuprates

The nature of the effective interaction responsible for pairing in the high-temperature superconducting cuprates remains unsettled. This question has been studied extensively using the simplified single-band Hubbard model, which does not explicitly consider the orbital degrees of freedom of the relevant CuO$_2$ planes. Here, we use a dynamic cluster quantum Monte Carlo approximation to study the orbital structure of the pairing interaction in the three-band Hubbard model, which treats the orbital degrees of freedom explicitly. We find that the interaction predominately acts between neighboring copper orbitals, but with significant additional weight appearing on the surrounding bonding molecular oxygen orbitals. By explicitly comparing these results to those from the simpler single-band Hubbard model, our study provides strong support for the single-band framework for describing superconductivity in the cuprates.

Introduction -Cuprate superconductivity emerges in their quasi-two-dimensional (2D) CuO 2 planes after doping additional carriers into these layers. The undoped parent compounds are charge transfer insulators due to the large Coulomb repulsion U dd on the Cu 3d orbitals, and, to a good approximation, a spin-1 2 hole is located on every Cu 3d x 2 −y 2 orbital. This situation is well described by a 2D square lattice Hubbard model or Heisenberg model in the large U dd limit.
Upon doping the additional holes or electrons primarily occupy the O or Cu orbitals, respectively. The minimal model capturing this asymmetry is the three-band Hubbard model, which explicitly accounts for the Cu 3d x 2 −y 2 , O 2p x , and 2p y orbitals (Fig. 1a) [1]. Even at finite doping, the low energy sector of the three-band model can be mapped approximately onto an effective single-band Hubbard model [2]. One expects this in the case of electron-doping since the additional carriers go directly onto the Cu sublattice, on which the holes of the undoped materials already reside. The case of holedoping, however, is more subtle. Here, the additional carriers predominantly occupy the O sublattice due to the large U dd on the Cu orbital, and the appropriateness of a single-band model is less clear. In their seminal work, Zhang and Rice [2] argued that the doped hole effectively forms a spin-singlet state with a Cu hole, the "Zhang-Rice singlet" (ZRS, Fig. 1b), which then plays the same role as a fully occupied or empty site in an effective single-band model, again facilitating a single-band description.
The nature of the single-band 2D Hubbard model's pairing interaction has been extensively studied [3][4][5][6][7][8]. Detailed calculations of its momentum and frequency structure using dynamic cluster approximation (DCA) quantum Monte Carlo (QMC) [3] find that it is well described by a spin-fluctuation exchange interaction [4]. The single-band model, however, cannot provide any information on the orbital structure of the interaction. For example, in the hole doped case, the spins giving rise to the spin-fluctuation interaction are located on the Cu sublattice, while the paired holes are moving on the O p x/y sublattice. This situation can produce a different physical picture than if the interaction and the pairs both originate from the same orbital on the same lattice [9][10][11][12][13]. And indeed, studies have observed two-particle behavior in a two sublattice system that is not observed in a onelattice system [14]. Moreover, an analysis of resonant inelastic x-ray scattering studies has found that a singleband model fails to describe the high-energy magnetic excitations near optimal doping [15]. Studying the effective interaction in a three-band model, and, in particular, determining its orbital structure is, therefore, critical. Such a study will also provide new insight into the nature of high-temperature superconductivity that is not available from the previous single-band studies. In this letter, we use a QMC-DCA method to explicitly calculate the orbital and spatial structure of the effective interaction in a realistic three-band CuO 2 model, and compare the results with those obtained from a single-band model.
Model and Methods -The three-band Hubbard model we study can be found in Refs. [16,17]. We adopted a parameter set appropriate for the cuprates [16,[18][19][20] (in units of eV): the nearest neighbor Cu-O and O-O hopping integrals t pd = 1.13, t pp = 0.49, on-site interactions U dd = 8.5, U pp = 0, and charge-transfer energy ∆ = ε p − ε d = 3.24, unless otherwise stated. Since we use a hole language, half-filling is defined as hole density n h = 1 and n h > 1 (< 1) corresponds to hole (electron)doping. A finite U pp only leads to small quantitative changes in the results (see Fig. S4 [17]) but worsens the sign problem significantly [16]. Therefore, we keep U pp = 0 for this study.
We study the single-and three-band Hubbard models using DCA with a continuous time QMC impurity solver [21][22][23]. We determine the structure of the pairing interaction by solving the Bethe-Salpeter equation (BSE) in the particle-particle singlet channel to obtain its leading eigenvalues and (symmetrized) eigenvectors [3,17]. A transition to the superconducting state occurs when the leading eigenvalue λ(T = T c ) = 1, and the magnitude of λ < 1 measures the strength of the normal state pairing correlations. The spatial, frequency, and orbital dependence of the corresponding eigenvector, which is the normal state analog of the superconducting gap, reflects the structure of the pairing interaction [3,8].
Results - Figure 1c shows the leading eigenvalue of the BSE for the three-band model as a function of hole concentration n h obtained on a 4 × 4 cluster with β = 1/k B T = 16 eV −1 . We find that it always corresponds to a d-wave superconducting state [24] and is larger for holedoping (n h > 1) compared to electron-doping (n h < 1). The latter observation suggests a particle-hole asymmetry in T c consistent with experiments and prior studies of the single and two-band Hubbard models [25,26]. (Although λ d is largest at half-filling, we expect that it asymptotically approaches one as the temperature decreases but never actually cross one due to the opening of a Mott gap. We observe such behavior in explicit calculations on smaller three-band clusters, see Fig. S1 [17]. ) We now analyze the spatial and orbital structure of the leading eigenvector φ αβ (k) (α and β denote orbitals), by Fourier transforming φ αβ (k) to real space to obtain φ r β (r α ), where r β denotes the position of the orbital taken as the reference site. We employed a 6×6 cluster to allow for long-ranged pairing correlations at T = 0.1 eV. While this relatively high temperature is needed to mitigate the Fermion sign problem, we have found that the leading eigenvector changes very slowly as the system cools (see Fig. S2) [17]. From here on, we focus on results obtained at optimal (15%) hole-or electron-doping. We have obtained similar results for different cluster sizes and for finite U pp (see Figs. S3 and S4) [17], indicating that our conclusions are robust across much of the model phase space.
In the single-band Hubbard model, the pairs are largely comprised of carriers on nearest neighbor sites in a d-wave state, i.e. with a positive (negative) phase along the x-(y)-directions. The internal structure of the pairs in the three-band model seems more complicated [27]. The real-space structure of φ r β (r α ) shown in Figs.2 ac and Figs.2 d-f for the hole-and electron-doped cases, respectively, display an extended and rich orbital structure. Here, the size and color indicate the strength and phase of φ r β (r α ), respectively, on each site after adopting the central Cu 3d x 2 −y 2 or O 2p x,y orbital as a reference site at r β . The form factors φ r β (r α ) are similar for both electron and hole doping, decaying over a length scale of ∼ 2-3 lattice constants. Moreover, while the d-wave pairing between nearest Cu sites dominates, there is also a significant contribution from d-p pairing, with a comparable amplitude for up to the third-nearest neighbors. The pairing between the individual O 2p x,y orbitals is much weaker in comparison.
We now transform the leading eigenfunction from the O-p x and p y basis to the bonding L and anti-bonding L combinations (Fig. 1d). These combinations, formed from the four O orbitals surrounding a Cu cation, are the relevant states for the ZRS, in which the doped holes are argued to reside in. The bonding L state strongly hybridizes with the central Cu 3d x 2 −y 2 orbital (Fig. 1b), The real space components of the leading particle-particle BSE (symmetrized) eigenvector for the three-band model at optimal doping and β = 10 eV −1 on a 6 × 6 cluster. Each column describes the pairing between a Cu d (or O px, py) reference site and all other orbitals as a function of distance. All panels set the Cu-d orbital at the origin, as labelled.
while the anti-bonding L state does not. The resulting antiferromagnetic exchange interaction between the Cu and L holes is then argued to bind them into the Zhang-Rice spin-singlet state, which provides the basis for the mapping onto a single-band model.
The orbital structure of the leading eigenfunction simplifies considerably after one transforms to the bonding L and anti-bonding L combinations. Fig. 3 plots the pairing amplitudes for a hole on Cu paired with another hole on a neighboring Cu (d-d, Fig. 3a) or bonding molecular orbital (d-L, Fig. 3b). Both components exhibit a clear d x 2 −y 2 symmetry; however, both channels also have indications of a higher momentum harmonic [cos(2k x a) − cos(2k y a)]. Interestingly, the contribution from holes occupying neighboring bonding molecular orbitals exhibits similar behavior (L-L, Fig. 3c). The Lrelated pairing contributes very little as will be discussed in Fig. 4 and in the supplement (see Fig. S5) [17].
Figs. 3a-c establishes that the pairing between the different orbital components of the ZRS all possess the requisite d x 2 −y 2 symmetry. This indicates that the ZRS picture is valid for describing pairing correlations in the three-band Hubbard model of the cuprates. To confirm this, we also computed the real-space structure of the leading particle-particle BSE eigenfunction in the singleband Hubbard model. Here, we considered cases with next-nearest-neighbor hopping t /t = 0 (Fig. 3d) and −0.2 (Fig. 3e), which are commonly used in the literature, as well as −0.453 (Fig. 3f), which we obtained by downfolding our three-band model parameters onto the single-band model [28,29]. The single-band model reproduces the short-range pairing structure of the threeorbital model, regardless of the value of t ; however, the second and third neighbor pairing is only captured correctly for t /t = −0.453. These results provide remarkable support for the validity of the ZRS construction but also indicate that single-band models may not capture the correct longer-ranged correlations without a suitable choice of t . The latter conclusion further underscores the crucial role of t for determining the superconducting properties of the single-band model [8,30,31]. Figure 3 shows that the structure of the leading eigenfunction φ αβ is closely linked to the orbital structure of the ZRS. Fig. 4 examines how this internal structure evolves with doping by plotting the orbital-dependent hole density (panel a) and the orbital composition of the eigenfunction φ αβ (panel b) on a 4 × 4 cluster (adequate to capture the essential pairing structure) at a lower temperature. Fig. 4a shows that the single hole per unit cell in the undoped case has approximately 65% Cu-d character, while 35% of the hole is located in the bonding O-L molecular orbital. With electron doping, there is a small decrease of n d /n h indicating that the holes are removed mainly from the Cu-d orbital. In contrast, with hole  Figure 4b shows that the total weight of the nearestneighbor pairing increases from ∼ 70% in the undoped case to almost 100% with either hole or electron doping. Since the BSE eigenvector reflects the momentum structure of the pairing interaction, this dependence can be understood from an interaction that becomes more peaked in momentum space as n h = 1 is approached. This behavior leads to a more delocalized structure of φ r β (r α ) and, therefore, a reduction of the relative weight of the nearest-neighbor contribution. The partial contributions to the nearest-neighbor pairing weight, D d and D dL , have a doping dependence very similar to the corresponding orbital densities n d and n L in panel a, closely linking the orbital structure of the pairing to the orbital makeup of the ZRS. The weight of the L contributions remains negligible over the full doping range [17].
Conclusions -We have determined the orbital structure of the effective pairing interaction in a three-band CuO 2 Hubbard model and shown that it simplifies considerably when viewed in terms of a basis consisting of a central Cu-d orbital and a bonding L combination of the four surrounding O-p orbitals. These states underlie the ZRS singlet construction that enables the reduction of the three-band to an effective single-band model. By explicitly comparing the three-band with single-band results, we show that the effective interaction is correctly described in the single-band model. In summary, these results strongly support the conclusion that a single-band Hubbard model provides an adequate framework to understand high-T c superconductivity in the cuprates.

The Hamiltonian of the three-band Hubbard model is
(S1) Here, d † i,σ (d i,σ ) creates (annihilates) a spin σ (=↑, ↓) hole in the copper d x 2 −y 2 orbital at site i; p † αjσ (p αjσ ) creates (annihilates) a spin σ hole in the oxygen p α (α = x, y) orbital at site j; for nearest neighbor, j = i ±x/2 (orŷ/2); n d iσ = d † iσ d iσ and n pα jσ = p † αjσ p αjσ are the number operators; d and p are the onsite energies of the Cu and O orbitals, respectively; µ is the chemical potential; t i,j is the nearest neighbor Cu-O hopping integral; t j,j is the nearest neighbor O-O hopping integral; and U dd and U pp are the on-site Hubbard repulsions on the Cu and O orbitals, respectively. The hopping integrals are parameterized [5] as t ij = P ij t pd and t jj = Q jj t pp , where P ij = 1 for j = i +ŷ/2 or j = i −x/2, P ij = −1 for j = i −ŷ/2 or j = i +x/2 and Q jj = 1 for j = j −x/2 +ŷ/2 or j = j +x/2 −ŷ/2, Q jj = −1 for j = j +x/2+ŷ/2 or j = j −x/2−ŷ/2. Throughout, we adopted (in units of eV): t pd = 1.13, t pp = 0.49, U dd = 8.5, U pp = 0, and ∆ = ε p − ε d = 3.24 [5][6][7][8], unless otherwise stated. Since we use a hole language, half-filling is defined as hole density n h = 1 and n h > 1 corresponds to hole-doping and n h < 1 corresponds to electron-doping. A finite U pp only leads to small quantitative changes in the pair structure (see Fig. S4), but worsens the sign problem significantly [5]. Therefore, we keep U pp = 0 for this study except for the results presented in Fig. S4.
The downfolded single-band Hubbard model is where c † iσ (c iσ ) creates (annihilates) an electron with spin σ at site i, t i,j = t and t for nearest-and next-nearestneighbor hoping, respectively, and zero otherwise. U is the on-site Hubbard repulsion, and µ is the chemical potential, which is adjusted to fix the electron filling. Throughout, we set t = 1, U = 6t, and vary t as indicated in the text.

SYMMETRIZED EIGENVECTORS OF THE BETHE-SALPETER EQUATION
To determine the structure of the effective pairing interaction, we solve the Bethe-Salpeter equation in the particleparticle singlet channel Here, K = (K, iω n ), andχ α1,α2,α3,α4 (K, ω n ) = (N c /N ) k [G α1α3 (K + k , ω n )G α2α4 (−K − k , −ω n )] is the coarsegained bare particle-particle propagator. The irreducible particle-particle vertex Γ c,pp is extracted from the twoparticle cluster Green's function G 2,c α1,α2,α3,α4 (K, K ) with zero center of mass momentum and frequency by inverting the cluster Bethe-Salpeter equation To remove the ambiguity between left and right eigenvectors of the eigenvalue equation (S3), we symmetrize the pairing kernel entering Eq. (S3). Using matrix notation in (K, α, β), we first diagonalize the bare particle-particle propagator, We use the eigenvectors of the symmetrized BSE, φ ν αβ (K), for the analysis presented in the main text. They are related to the right eigenvectors of the BSE in Eq. (S3) by (S6)

THE BASIS TRANSFORMATION TO THE MOLECULAR L, L ORBITALS
The construction of the Zhang-Rice singlet relies on a transformation from the oxygen p x , p y orbital basis to bonding and anti-bonding molecular orbitals, denoted here as L and L , respectively. The two basis are related by a unitary transformation [1][2][3] and where γ 2 k = sin 2 (k x a/2) + sin 2 (k y a/2), p αkσ = N −1/2 c j p αjσ exp(−i k · R j ), and we have set the lattice constant a = 1. In this basis, only the L state hybridizes with the Cu-d orbital, while the L state only hybridizes with the L state. The Fourier transform of the L and L orbitals to real-space is defined as

SUPERCONDUCTING TRANSITION TEMPERATURE IN THE THREE-BAND HUBBARD MODEL
For the 6 × 6 and 4 × 4 DCA calculations presented in the main text, the QMC Fermion sign problem prevents calculations down to temperatures low enough to determine the superconducting transition temperature T c from the temperature where the leading eigenvalue of the Bethe-Salpeter equation crosses 1, i.e. λ d (T = T c ) = 1. This temperature can be reached on a 2 × 2 cluster, however, and T c (n h ) can be determined as a function of hole density n h in that case. Fig. S1 shows the DCA results for T c (n h ) obtained in a 2 × 2 cluster for the same model parameters as used in the main text. Similar to the electron-hole asymmetry found in Fig. 1 c for the leading eigenvalue λ d (n h ) of the particle-particle Bethe-Salpeter equation, as well as the asymmetry found in experiments, the T c versus n h phase diagram exhibits a higher maximum T c on the hole doped side than on the electron-doped side. Moreover, these results are similar to previous DCA 2 × 2 cluster calculations for a similar two-band model [4], although the critical hole doping where T c vanishes is reduced compared to those earlier calculations. This difference may originate in the difference in model parameters, in particular the neglect of the direct oxygen-oxygen hopping t pp in the earlier two-band model calculations.

DEPENDENCE OF THE LEADING EIGENFUNCTION ON TEMPERATURE, CLUSTER SIZE AND OXYGEN COULOMB REPULSION
While the leading eigenvalue λ d (T ) shows a very strong increase with decreasing temperature, the temperature dependence of the corresponding eigenfunction is found to be rather weak. Fig. S2 shows how the orbital and spatial structure of the leading eigenfunction φ r β (r α ) of the (symmetrized) Bethe-Salpeter equation changes with decreasing temperature between β = 1/T = 10 eV −1 (top panels a-c from Fig. 2 in the main text) and β = 12 eV −1 (bottom panels d-f). We only observe small quantitative changes between these two temperatures.
The cluster size dependence of the leading eigenfunction is studied in Fig. S3, which shows the results of an N Cu = 4 × 4 cluster calculation for a 15% hole doped and a 15% electron doped system. These results should be compared with Fig. S2 (or Fig. 2 in the main text), which displays the same calculation for a larger N Cu = 6 × 6 cluster. From this comparison, one sees that the 4 × 4 cluster is large enough to contain the important components of the eigenfunction. Since the Fermion sign problem is much less severe in the 4 × 4 cluster, it allows for calculations at lower temperatures or with an additional on-site Coulomb repulsion U pp on the oxygen orbitals.
An additional U pp = 4.1 eV term is considered in the data for the pair structure shown in Fig. S4. Other model parameters are unchanged from those considered in Fig. S3. Comparing these images with those in Fig. S3, one sees that the structure of the eigenfunction remains almost unchanged by the additional U pp . Only a very slight suppression of the components that involve the O-p orbitals is observed. This justifies the neglect of U pp in most of our calculations, and provides evidence that our main conclusions reached from those calculations are general and not affected by U pp .

ROLE OF THE ANTI-BONDING MOLECULAR L ORBITAL
Finally, we show the components of the leading eigenfunction that involve the antibonding L orbital in the bottom row of Fig. S5, compared to the bonding L components that were already shown in Fig. 3 of the main text. From the results for the orbital hole densities in Fig. 4 a of the main text, it is clear that the L molecular orbital remains almost completely unoccupied over the entire doping range considered, despite the finite hybridization between the L and L states. As a consequence, and as seen from the bottom row of Fig. S5, the L state is not involved in the pairing. This provides strong support for the Zhang-Rice picture, which only considers the d-and L-states in the mapping to an effective single-band model. The middle column describes pairings with respect to a px orbital reference. The right column describes pairings with respect to a py orbital orbital reference. All panels set the Cu 3d orbital at the origin, as labeled by the open ring. Compared with Fig. 2 in the main text, the 4 × 4 cluster contains the same essential pair structure as the 6 × 6 cluster and makes it possible to explore lower temperatures and stronger interactions.