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 CuO2 planes. Here, we use a dynamical 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-twodimensional (2D) CuO 2 planes after doping additional carriers into these layers. The undoped parent compounds are chargetransfer insulators due to the large Coulomb repulsion U dd on the Cu 3d orbitals and comparatively smaller charge transfer energy, 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 half-filled 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 hole doping, 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 singleband 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 dynamical 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 holedoped 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 one-lattice system 14 .
Moreover, an analysis of resonant inelastic X-ray scattering data has found that a single-band model fails to describe the highenergy magnetic excitations near optimal doping 15 . For these reasons, numerous numerical studies of the three-band Hubbard model have been carried out to date [16][17][18][19][20][21][22][23][24][25][26][27][28] ; however, the crucial task of studying its effective interaction, and, in particular, determining its orbital structure is currently lacking. Such a study will also provide new insight into the nature of high-temperature superconductivity that is not available from the previous singleband studies. Here, 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.

RESULTS
Pairing structure of the three-band model To study the structure of the pairing interaction, we solved the Bethe-Salpeter equation (BSE) in the particle-particle singlet channel to obtain its leading eigenvalues and eigenvectors (see the online supplementary materials) 3 (see "Methods"). 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 29 and is larger for hole doping (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 30,31 . (Although λ d is largest at half-filling, we expect that it asymptotically approaches one as the temperature decreases but never actually crosses one due to the opening of a Mott gap. We observe such behavior in explicit calculations on smaller threeband clusters (see Supplementary Fig. 1 in the online supplementary materials).
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 Supplementary Fig. 2 in the online supplementary materials). For example, we can reach much lower temperatures on 2 × 2 clusters, where we resolve the superconducting T c explicitly (see Supplementary Fig. 1 in the online supplementary materials). In that case, we observe that while the eigenvalue has a strong temperature dependence near T c , its corresponding eigenvector does not vary much with temperature. 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 Supplementary Figs. 3 and 4 in the online supplementary materials), 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 32 . The real-space structure of ϕ r β ðr α Þ shown in Fig. 2a-f for the hole-and electron-doped cases, respectively, display an extended and rich orbital structure. Here, the size and color of the data points 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 (unit-cell) neighbors. The pairing between the individual O 2p x,y orbitals is much weaker in comparison.
Pairing structure in the molecular basis We now transform the leading eigenvector for the hole-doped case from the O p x and p y basis to the bonding L and anti-bonding L 0 basis (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), while the anti-bonding L 0 state does not. The resulting antiferromagnetic exchange interaction between the Cu and L holes is then argued to bind them into the Zhang-Rice spinsinglet state, which provides the basis for the mapping onto a single-band model.
The orbital structure of the leading eigenvector simplifies considerably after one transforms to the bonding L and antibonding L 0 combinations. Figure 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 that is dominated by the (nearest-neighbor) cosðk x aÞ À cosðk y aÞ harmonic; however, both channels also have indications of additional higher-order harmonics [i.e., cosð2k x aÞ À cosð2k y aÞ and cosð2k x aÞ cosðk y aÞ À cosð2k y aÞ cosðk x aÞ]. Interestingly, the contribution from holes occupying neighboring bonding molecular orbitals exhibits similar behavior (L-L, Fig. 3c). The L 0 -related pairing contributes very little as will be discussed in Fig. 4 and in the supplement (see Supplementary Fig. 5 in the online supplementary materials). Figure 3a-c establishes that the pairing between the different orbital components of the ZRS all possesses the requisite d x 2 Ày 2 symmetry. This observation indicates that the ZRS picture-a singlet state made up of holes in the d and L orbitals-is valid for describing pairing correlations in the three-band Hubbard model of the cuprates. To show the pair structure for the ZRS, we plot in Fig. 3d the sum over the d-d, d-L, L-d, and L-L components (with a factor of 0.5 applied to all components). One sees that the ZRS pair structure has a vanishing cosð2k x aÞ À cosð2k y aÞ component, while higher-order harmonics remain.
To compare this result to the single-band picture, we also computed the real-space structure of the leading particle-particle BSE eigenvector for the single-band Hubbard model. Here, we considered cases with next-nearest-neighbor hopping t 0 =t ¼ 0 (panel e), −0.2 (f), −0.3 (g), which are commonly used in the literature, as well as −0.4 (h). The single-band model reproduces the short-range pairing structure of the three-orbital model (panel d), regardless of the value of t 0 ; however, the longer-ranged pairing in Fig. 3d is only captured correctly for large jt 0 =tj. In particular, we observe that with increasing jt 0 =tj, the relative amplitude of the third-nearest-neighbor ½cosð2k x aÞ À cosð2k y aÞ term is suppressed. For t 0 =t ¼ À0:4 (panel h), the single-band pair structure is very similar to that for the ZRS (panel d), with differences appearing at the longest length scales. This value of t 0 is close to the value t 0 ¼ À0:453t that we obtain by downfolding our three-band model parameters onto the single-band model by diagonalizing small Cu 2 O 7 clusters 33,34 . A sizeable negative t 0 is also consistent with parametrizations of the bandstructure extracted from angle-resolved photoemission spectroscopy 35 . 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 0 . The latter conclusion further underscores the crucial role of t 0 for determining the superconducting properties of the single-band model 8,31,36,37 .
Weights for orbital-resolved pair components Figure 3 shows that the structure of the leading eigenvector ϕ αβ is closely linked to the orbital structure of the ZRS. Figure 4 examines how this internal structure evolves with doping by plotting the orbital-dependent hole density (panel a) and the orbital composition of the eigenvector ϕ αβ (panel b) on a 4 × 4 cluster (adequate to capture the essential pairing structure) at a lower temperature. Figure 4a shows that the single hole per unit cell in the undoped case has~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 doping, there is a significant redistribution of the hole density from the d-to the Lorbital, showing that doped holes mainly occupy the O-L molecular orbital. The hole density on the anti-bonding O-L 0 orbital is negligible. Figure 4b shows that the total weight of the nearest-neighbor 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 nearestneighbor contribution. The partial contributions to the nearestneighbor 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 0 contributions remains negligible over the full doping range (see Supplementary  Fig. 5 in the online supplementary materials).

DISCUSSION
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 showed 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.

Model parameters
The three-band Hubbard model we study can be found in ref. 18 (see the online supplementary materials). We adopted a parameter set appropriate for the cuprates 18,[38][39][40] (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 Supplementary  Fig. 4 in the online supplementary materials) but worsens the sign problem significantly 18 . Therefore, we keep U pp = 0 for this study. Fig. 2 Pairing structure of the three-band model. 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 p x , p y ) reference site and all other orbitals as a function of distance. All panels set the Cu-d orbital at the origin, as labeled.

Dynamical cluster approximation
We study the single-and three-band Hubbard models using DCA with a continuous-time QMC impurity solver [41][42][43][44][45] . We determine the structure of the pairing interaction by solving the BSE in the particle-particle singlet channel to obtain its leading eigenvalues and (symmetrized) eigenvectors 3 (see the online supplementary materials). 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 .

DATA AVAILABILITY
The data that support the findings of this study can be obtained at https://github. com/JohnstonResearchGroup/Mai_etal_3bandPairs_2021. Fig. 3 A comparison of the pairing structure in the three-band and single-band models. Each panel plots the real-space components of the leading particle-particle BSE (symmetrized) eigenvector at 15% hole doping. The first row shows d-d, d-L, L-L, and ZRS-ZRS pairing components for the three-band model at β = 10 eV −1 . The second row shows the pair structure for the single-band model (U = 6t, β = 5t −1 ) at t 0 =t ¼ 0, −0.2, −0.3, and −0.4. Fig. 4 Orbital-resolved density and pair components. a Ratios of the orbital-resolved hole densities to the total density n h . b Weights of different orbital components of the nearest-neighbor pairs, as defined in Fig. 1, and their total weight. All results were obtained on a 4 × 4 cluster at β = 16 eV −1 . P. Mai et al.