Coupled electronic and magnetic excitations in the cuprates and their role in the superconducting transition

The formation of Cooper pairs, a bound state of two electrons of opposite spin and momenta by exchange of a phonon, is a defining feature of conventional superconductivity. In the cuprate high temperature superconductors, even though the superconducting state also consists of Cooper pairs, the pairing mechanism remains intensely debated. Here, we investigate superconducting pairing in the Bi2Sr2CaCu2O8+δ (Bi2212) cuprate by employing spectral functions obtained from angle-resolved photoemission as input to the Bethe-Salpeter equation. Assuming Cooper pairing is driven by spin fluctuations, we construct the spin-fluctuation-mediated pairing interaction and use it to compute the eigenfunctions and eigenvalues of the Bethe-Salpeter equation for multiple Bi2212 samples. The leading d-wave eigenvalue increases as the temperature is decreased toward Tc, reaching a value of approximately 1 at the Tc corresponding to each doping value. This suggests that spin fluctuations can approximately account for Tc and mediate pairing in the cuprate superconductors. The pairing mechanism behind the Cooper pairs in unconventional superconductors has been debated since the discovery of the high-Tc cuprates. Here, using angle resolved photoemission spectroscopy data for Bi2212 applied to the Bethe-Salpeter equation, the authors suggest that a spin-mediated pairing interaction can account for the high superconducting temperature as well as play a dominant role in the Cooper pair formation.

S uperconductivity in the cuprates emerges following quantum melting of the antiferromagnetic Mott insulating state of the parent compound, either via electron or hole doping. It has been proposed that the spin fluctuations resulting from the melted Mott state act as the "glue" leading to Cooper pairing in high-temperature superconductors 1,2 . In support of this proposal, neutron scattering experiments on the YBa 2 C 3 O 6.95 cuprate superconductors demonstrate that the change in magnetic exchange energy between the superconducting and normal states can provide sufficient superconducting condensation energy 3 . Moreover, it has been shown that the interaction of electrons with spin fluctuations can account for several anomalies in charge, spin, and optical response measurements in the cuprates 2 . An outstanding question is whether the spin-fluctuation-mediated pairing interaction is strong enough to yield the high transition temperatures observed in these materials.
While several studies have attempted to answer this important question by using a combination of theoretical approaches and various experimentally measured electronic and spin excitation spectra, the results have been inconclusive. Using the electron spectral functions from angle-resolved photoemission spectroscopy (ARPES) experiments in near-optimally doped Bi2212 and solving a generalized gap equation, Mishra et al. 4 did not obtain a superconducting transition at any temperature. In contrast, Dahm et al. 5 , using the spin excitation spectrum measured via inelastic neutron scattering experiments (INS) in YBa 2 C 3 O 6.6 and computing the electron spectral function theoretically rather than utilizing experimental ARPES data, obtained a superconducting transition temperature of about 150 K. In a subsequent theoretical study, Maier et al. 6 ascribed the failure of Mishra et al. to find any superconducting transition to the fact that the experimental ARPES data were taken at a temperature of 140 K, well above the critical temperature of T c = 90 K. Taken together, these studies suggest that, in order to reliably determine whether the spinfermion interaction can be the source of the unconventional superconducting state in the cuprate superconductors, three requirements need to be satisfied: (i) experimental ARPES data need to be taken above, but near T c , (ii) the effective superconducting pairing vertex needs to be computed from these ARPES data, and (iii) the solution of the gap equation, using the ARPES data and the derived pairing vertex as input, should reproduce the transition temperature of the material.
Here, we use experimental ARPES excitation spectra from Bi2212 near and above T c to construct the spin-fluctuationmediated pairing interaction and solve the gap equation. Following previous work 4-7 , we employ the Bethe-Salpeter equation for a generalized order parameter Φðk; ω n ; TÞ: The transition to superconductivity takes place when the eigenvalue λ reaches λ(T c ) = 1. The operatorÔ T ð Þ is given by (we solve the equation in Matsubara frequency space): where G is the electron's Green's function and V the effective electron-electron pairing interaction. Equations (1) and (2) are quite universal in that the specific mechanism leading to the emergence of superconductivity is reflected only in the form of V, as noted previously by Mishra et al. For the three samples considered (see "Results-Bethe-Salpeter Equation"), we find solutions to this equation which approximately reproduce their observed T c values.

Results and discussion
Spectral function. The required quantities for the Bethe-Salpeter equation can be obtained, within certain approximations, from the electron spectral function Aðk; ωÞ obtained from the ARPES signal I k; where I 0 is the matrix element, f ðωÞ the Fermi function, and BðωÞ the inelastic background signal. To extract the spectral function, we first subtract the background signal, while simultaneously normalizing out the matrix element I 0 8 . We then note that, although ARPES only measures the occupied states, the cuprates exhibit particle-hole mixing below 9 and slightly above 10 T c . Therefore, the symmetrized spectral function Aðk; ωÞ ¼ Iðk; ωÞ þ IðÀk þ 2k F ; ÀωÞ also yields the unoccupied part, for k values close to the Fermi momentum k F . The validity of this assumption in this work is justified by (i) the approximate insensitivity of the pairing interaction constructed from ARPES data to particle-hole asymmetry put in by hand 8 and (ii) the robustness of the Bethe-Salpeter eigenvalues (see "Results-Bethe-Salpeter Equation") against the stronger assumption Aðk; ωÞ ¼ Aðk; ÀωÞ. Finally, the spectral function is normalized to unit area within the binding energies probed in our experiments. (For more details, see "Methods-Data analysis" and Supplementary Note 2) The Green's functions in the Matsubara representation then follow from the relation We analyze ARPES data for three Bi2212 samples: two thin films, one underdoped with T c = 67 K (UD67), the other optimally doped with T c = 80 K (OP80), and a single crystal, optimally doped with T c = 91 K (OP91). (The data for the thin films have been presented elsewhere 11 .) Representative data are shown in Fig. 1. Figure 1a shows the experimentally determined phase diagram for our samples 12 , with the circles, diamonds, and squares indicating the spectra used in this work. In Fig. 1b, we plot the spectral function in the Y-quadrant of the Brillouin zone, with a tight-binding fit to the Fermi surface 13 appearing as a continuous yellow line. The measured spectra for the UD67 sample along the high symmetry directions, highlighted by the coloured dots in panel b, are shown in Fig. 1c.
Spin susceptibility. For pairing induced by spin fluctuations, we consider the pairing vertex 14 where χ is the spin susceptibility and U is an effective spinfermion coupling energy. The spectral function allows us to calculate the (real frequency) spin fluctuation propagator in the random phase approximation (RPA) 15 : where χ 0 is the bare susceptibility, and the coupling U q is assumed to have a super-exchange momentum dependence 13,16-18 : (see Supplementary Note 4). Here, U 0 reflects the strength of the antiferromagnetic interaction in the particle-hole channel. As such, it is in general distinct from the spin-fermion coupling U 19 .
However, our estimates of U 0 and U (see Supplementary Note 5) indicate that these two energies differ by <11%, and we will therefore take them to be equal for simplicity. Since U 0 is a high energy scale not provided by the RPA nor the low-energy formalism adopted here, its value was chosen so as to reproduce the salient features of the spin susceptibility measured in INS experiments. In particular, in the superconducting state, the cuprates exhibit a sharp resonance in their spin excitation spectra at the commensurate wave-vector Q ¼ ðπ; πÞ: 20,21 The energy of this resonance throughout the phase diagram of the cuprate YBa 2 Cu 3 O 6+x has been extensively reported; however, due to technical limitations, such information is lacking for underdoped Bi2212. For this reason, we used the empirical relation 21-23 to fix the resonance energies for the three samples (E 67 R ¼ 32 meV, E 80 R ¼ 37 meV, and E 91 R ¼ 41 meV). We then adjusted U 0 such that ImχðQ; ΩÞ exhibited resonance peaks at the energies given by Eq. (8). As the temperature increases, INS data suggest that the resonance broadens and decreases in intensity, but disperses weakly with temperature, even above T c 24-26 . This observation, which is consistent with the reported near-constancy of the ARPES antinodal gap in the temperature range considered here 11 , justifies the use of Eq. (8) even above T c . Near T c , the values of U 0 extracted from this procedure are~730 meV for UD67,~650 meV for OP80 and~580 meV for OP91 (see Supplementary Note 3 for a discussion of the contribution from the anomalous Green's function and Supplementary Note 6 for the effects of varying E R ). Figure 2 shows the calculated spin susceptibilities (Imχ) for the two Bi2212 films (additional results for the single crystal can be found in Supplementary Note 1). Here, we include results well below T c to compare with the more extensive INS data in this regime. The commensurate peak intensity decreases with increasing temperature and doping (Fig. 2a, b), as seen in INS experiments on YBa 2 Cu 3 O 6+x 24,27 . Also, the resonance peak widths deep in the superconducting state are~15-20 meV (blue curves in Fig. 2a, b), in agreement with INS data from optimally doped Bi2212 20 . In addition to the resonance, our calculations also reproduce the ubiquitous "upward branch" (ω>E R ), observed in all cuprate families, but not the "lower branch", which is material-specific 5,28 (Fig. 2c, d). The energy-momentum integral of Imχðq; ΩÞ (multiplied by the matrix element 2μ 2 B 13 ) gives a total fluctuating moment of m 2 $ 0:43μ 2 B for the UD67 sample and m 2 $ 0:39μ 2 B for OP80, both varying weakly with temperature. Moreover, for the UD67 sample at 20 K, the total spectral weight taken up by the resonance (integrated over all momenta and from 0 to 70 meV) is about 0:16μ 2 B , whereas the remaining spectral weight (up to 200 meV), $ 0:27μ 2 B . These values compare reasonably well with the corresponding figures for YBa 2 Cu 3 O 6.6 (with a similar T c ) 29 . While the use of weakcoupling approaches like the RPA may be called into question in the strongly correlated cuprates, the above results suggest that this approximation can semi-quantitatively account for the spin Black circles correspond to the UD67 sample, red diamonds to OP80, and blue squares to OP91. There are two superconducting domes, one dashed for the thin films and one solid for the single crystal, with a higher optimal T c . b Raw Fermi Surface (after interpolation and reflection about symmetry axes) for the UD67 sample at 20 K. The thin yellow line is a tight-binding fit to the Bi2212 dispersion 13 and is used here as a guide to the eye. c Energy Distribution Curves corresponding to (b) along the path indicated by the coloured dots. Throughout this article, we take a = 1 for the lattice constant.
susceptibility probed in INS experiments (as has been reported before 4,8 ), which is a necessary ingredient in the calculation of the Bethe-Salpeter eigenvalues.
Bethe-Salpeter equation. With the experimentally derived values of U 0 and the functions Gðk; ω n Þ and Imχðq; ΩÞ, Eq. (1) can be solved after obtaining the Matsubara frequency representation for the spin response function: For the linearized Bethe-Salpeter equation, we considered only data sets near and above T c , and solved for the leading eigenvalue and eigenvector Φðk; ω n ; TÞ using the power method 30 (see "Methods-Numerical analysis"). Figure 3a ð Þ À1 $ 10 fs, as has been observed in ultrafast optical spectroscopy measurements 32 . Finally, the calculated Φ allows us to recover the superconducting gaps measured in ARPES along the Fermi surface, as shown in Fig. 3e for the UD67 film, where Φ was calculated from the data at 62 K and the ARPES gap values were obtained from the same film at 20 K. Figure 3f shows the Bi22212 Fermi surface indicating the angle plotted in Fig. 3e. Similarly, the pairing eigenfunction calculated from superconducting state (T ¼ 70 K) data from the OP80 film and the ARPES gaps from the same film at 20 K show remarkable agreement, as shown in Fig. 3g. Figure 4 summarizes our results for the leading eigenvalues at temperatures slightly below and above T c , for both thin films and the single crystal. Note the temperature dependence of the eigenvalues (Fig. 4a-c); for the three samples, the eigenvalues increase as T is lowered toward T c , as observed in some Hubbard model calculations 6,7,33 . This temperature dependence is a direct consequence of using temperature-dependent, experimental spectral functions. Our main result is that, to the accuracy of the present calculation (see Supplementary Note 6), the eigenvalues are essentially equal to unity near T c .
It is worth noting that the temperature extracted from the condition λ T ð Þ ¼ 1 in this work is not the temperature scale associated with incoherent pair formation (T pair ) reported recently 34 . Indeed, we solve the Bethe-Salpeter equation in momentum space, which implicitly assumes the existence of a phase-coherent, long-range superconducting order parameter (as the Greens functions and the pairing interaction that enter Eq. (2) are coherent). As such, the temperature at which the eigenvalue of Eq. (1) reaches 1 must be identified with T c .

Conclusions
Our results show that the magnetically mediated pairing interaction, constructed from ARPES data within the RPA, is strong  a Momentum dependence along the ðπ; 0Þ ! ð0; πÞ direction (as shown in the inset of (b)) of the eigenfunction Φ for the UD67 sample at T = 62 K (blue) and T = 80 K (red). b Same as (a) but for the OP80 sample at T = 70 K (purple) and T = 90 K (green). The black, dashed line in (a) and (b) is the d-wave function 1 2 ðcosk x À cosk y Þ along the path k y ¼ π À k x . c Matsubara energy dependence of the eigenfunction Φ at ðπ; 0Þ for the UD67 sample at the temperatures indicated in (a). d Same as (c) but for the OP80 sample. e Comparison of angle-resolved photoemission spectroscopy (ARPES) gaps (Δ 0 , grey triangles) for UD67 at 20 K with the calculated Φðk; k B TÞ at 62 K (blue line) along the Fermi surface shown in (f). g Same as (e), but the ARPES gap values correspond to the OP80 sample at 20 K and the pairing eigenfunction was calculated from the OP80 data at 70 K (purple line). In (e), (g) the pairing eigenfunctions were re-scaled to match the ARPES gaps at θ ¼ 0 and the gap values were adapted from Kanigel et al. 11 . The curve colours are consistent throughout each column of panels. enough to mediate high-temperature, d-wave superconductivity in the cuprates at optimal and sub-optimal doping. We note that the emergence of superconductivity in the underdoped sample (UD67) arises from the increased strength of the pairing interaction which compensates for the loss of low-energy electronic states in the pseudogap region. This result is consistent with the notion that antiferromagnetic correlations promote high-T c superconductivity in a dry Fermi sea 6 .

Methods
Experimental. Bi2212 thin films were epitaxially grown by radio frequency magnetron sputtering 35 up to a thickness of about 2000 Å. The resistivity curves for these films showed transition widths of about 10 K for the OP80 film (optimally doped) and 15 K for the UD67 film (underdoped) 36 . These were measured with 22 eV photons at the U1 undulator beamline at the Synchrotron Radiation Center in Madison, Wisconsin. The spectra were collected with a Scienta R2002 spectrometer, with an energy resolution of 15 meV. The films were mounted with the Cu-O bond parallel to the photon polarization and cleaved in situ at a pressure below 2 10 À11 Torr. The optimally doped single crystal was grown by the travelling solvent, floating zone technique. ARPES measurements for this sample were performed at the University of Illinois at Chicago, using the unpolarized He-I line (21.2 eV) from a Helium discharge lamp. Photoelectrons were analyzed with a Scienta R4000 spectrometer, the energy resolution of the system was set at 25 meV, and the sample was measured at a pressure below our measurement limit of 2 10 À11 Torr. Throughout all the experiments, antinodal spectra were repeatedly measured and compared with the corresponding spectra collected at the beginning of each run. These measurements showed no signs of surface contamination or doping (see Supplementary Note 2).
Data analysis. Starting from the raw ARPES data, we first subtract the contribution due to second order light. In the second step, an EDC far away from the Fermi momentum k F in the "unoccupied" side is used as an energy-dependent background. This background is subtracted from the EDCs at each momentum point after normalizing it to each EDC at a certain binding energy ω 0 (640 meV for the thin films and 510 meV for the single crystal). The reason why we can subtract the same background EDC throughout the Brillouin Zone is that this spectrum is roughly momentum-independent, as we reported elsewhere 37 . We also checked that changing ω 0 within an energy window of 140 meV does not significantly alter the end results. Subsequently, we divide the spectra by a resolution-broadened Fermi function to get an effective ImGðk; ωÞ. Afterwards, we determine the unoccupied part of ImGðk; ωÞ by adopting particle-hole symmetry with respect to the Fermi surface: ImGðk þ k F ; ωÞ ¼ Im Gðk F À k; ÀωÞ where k is normal to the Fermi surface.
We then normalize the derived ImGðk; ωÞ so that the integral of the spectral function Aðk; ωÞ ¼ À1=π ImGðk; ωÞ is equal to unity over the energy range ½Àω 0 ; ω 0 . Even though the above-described steps cannot be claimed to be exact, they collectively minimize the effect of dipole matrix elements. Finally, we interpolate the data on to a rectangular grid and apply reflections appropriate to the symmetries of the CuO 2 planes to obtain the spectral function across most of the BZ. (See Supplementary Note 2 for more details.) In previous work 38, 39 we have used the above data normalization procedure to conduct an autocorrelation analysis of ARPES data and compare quantitatively with Fourier transformed STS (FT-STS) data, with excellent agreement. In addition, our current analysis provides estimates for the spin-fermion coupling constant U 0 and the time scale for the pairing interaction in cuprate HTSCs, which are consistent with reported values in the literature 32 . This agreement with experiments provides confidence in our methodology for extracting Aðk; ωÞ from ARPES data.
Numerical analysis. The method employed in solving the eigenvalue problem in Eqs. (1) and (2) followed very closely that of Monthoux 30 . In assessing whether a hypothetical pairing interaction is strong enough to mediate superconductivity at a given temperature, it is sufficient to focus on the largest eigenvalue ofÔ. To this end, we used the power method, which is an iterative procedure to calculate the maximum eigenvalue of an operator. There are several equivalent ways of carrying out this procedure; in this work, the implementation goes as follows. We start with a guess for the eigenfunction Φ, which is a generalized gap function reflecting the momentum and energy dependence of the pairing interaction. Motivated by the well-established d-wave form of this interaction, we take Φ 0 ðk; ω n Þ ¼ cos k x À cos k y as our starting guess.
With this guess, we perform the first iteration F 1 ¼ÔΦ 0 and factor out the largest entry in F 1 , so that F 1 ¼ λ 1 Φ 1 , where λ 1 ¼ maxðF 1 Þ and Φ 1 ¼ F 1 =λ 1 . We then use Φ 1 as input to the second iteration F 2 ¼ÔΦ 1 to obtain λ 2 ¼ maxðF 2 Þ and Φ 2 ¼ F 2 =λ 2 , as before. The above process is repeated until the iteration (N) is reached where the difference λ N À λ NÀ1 is negligible. At this point, the function Φ N is the normalized eigenfunction ofÔ corresponding to the largest eigenvalue λ N .
Since Eq. (2) is essentially a convolution equation, the operationsÔΦ j were performed using Fourier techniques for increased computational speed. The numerical efficiency afforded by the convolution structure of the Bethe-Salpeter equation was our main reason for studying the pairing problem in Matsubara frequencies.
From the Lehmann representation of the Green's function, it can be seen that Pðk; ω n Þ ¼ Gðk; ω n ÞGðÀk; Àω n Þ is even in Matsubara frequency, real, and positive. For the Matsubara spin response, the property Imχðq; ÀΩÞ ¼ ÀImχðq; ΩÞ implies that χðq; ω m Þ is also even in ω m , real and positive. These properties are not only useful as consistency checks of the numerical procedure, but also enable us to use the power method, which breaks down when the eigenvalues are complex.

Data availability
One of the data sets supporting the findings of this study is available at https:// github.com/Francisco-UIC/Spin_fluctuations/. Other data sets are available upon reasonable request.

Code availability
The codes employed in this study are available at https://github.com/Francisco-UIC/ Spin_fluctuations/.