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 [1], is a defining feature of conventional superconductivity. In the cuprate high temperature superconductors, even though it has been established that the superconducting state also consists of Cooper pairs, the pairing mechanism remains intensely debated. Here we investigate superconducting pairing in the Bi2Sr2CaCu2O8+\delta (Bi2212) cuprate by employing spectral functions obtained directly from angle-resolved photoemission (ARPES) experiments as input to the Bethe-Salpeter gap equation. Assuming that Cooper pairing is driven solely by spin fluctuations, we construct the single-loop spin-fluctuation-mediated pairing interaction, and use it to compute the eigenfunctions and eigenvalues of the Bethe-Salpeter equation in the particle-particle channel for multiple Bi2212 samples. The key point of our results is that, as the temperature is reduced, the leading eigenvalue increases upon approaching Tc, reaching a value of approximately 1 at the Tc corresponding to each doping value, indicating a superconducting transition with d-wave eigenfunctions. This suggests that spin fluctuations can approximately account for Tc and, consequently, mediate pairing in the cuprate high temperature superconductors.

response measurements in the cuprates [3]. However, an outstanding question is whether the spinfermion 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 ARPES experiments in near-optimally doped Bi2212, Mishra et al. [5] solved a generalized gap equation, but did not find a superconducting transition for any temperature. In contrast, Dahm et al. [6] employed the spin excitation spectrum measured via inelastic neutron scattering experiments (INS) in YBa2Cu3O6.6 to obtain a superconducting transition temperature of about 150 K. However, in this study, the electron spectral function was theoretically computed, and not taken from experimental ARPES data. In a subsequent theoretical study, Maier et al. [7] ascribed the failure of Mishra et al. [5] 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 TC = 90 K. Taken together, these studies have suggested that, in order to determine whether the spin-fermion 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 Tc, (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 ARPES data from three underdoped and optimally doped Bi2212 samples, taken at temperatures near and above Tc, to construct the spin-fluctuation-mediated pairing vertex and solve the generalized gap equation. This yields a superconducting order parameter with 2− 2 −wave symmetry, and a critical temperature for all samples which is close to the ones observed experimentally.
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. [5].
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 ( , ) = 0 ( ) ( , ) + ( ), where I0 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 I0 [9]. We then note that, although ARPES only measures the occupied states, the cuprates exhibit particle-hole mixing below [10] and above [11] Tc. Therefore, the symmetrized spectral function ( , ) = ( , ) + (− + 2 , − ) also yields the unoccupied part, for k values close to the Fermi momentum kF.
The Green's functions in the Matsubara representation then follow from the relation We analyze ARPES data from three Bi2212 samples: two thin films, one underdoped with Tc = 67 K (UD67), the other optimally doped with Tc = 80 K (OP80), and a single crystal, optimally doped with Tc = 91 K (OP91). 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 For pairing induced by spin fluctuations, we consider the pairing vertex [14] ( , ) = where χ is the spin susceptibility and U0 is an effective spin-fermion 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 spin susceptibility,   With the experimentally derived values of U0 and the functions G(k,ωn) and Imχ(q,ω), equation (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 Tc, and    These results show that the magnetically mediated pairing interaction, constructed from ARPES data within the RPA, is sufficiently strong 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 AFM correlations promote high-Tc superconductivity in a "dry Fermi sea" [7]. 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.    Using these data, we followed the same procedure as in the main text to calculate the spin susceptibility and solve the BSE. Supplementary Figure 2 summarizes the results for the OP91 crystal, emphasizing those obtained at 90 K (∼ Tc). In real frequencies, Imχ again broadens and loses intensity as the temperature is increased (Supplementary Figure 2a)

The anomalous Green's function
In the superconducting state, the particle-hole bubble χ0 should include a term quadratic in the anomalous Green's function, which we call FF. We ignored this term in this work for three reasons: As pointed out by Chatterjee et al. [3], in the superconducting state, the calculated dynamic susceptibilities vary only slightly when one includes the anomalous Green's function. Inclusion of the FF term will certainly influence the bare susceptibility, but an appropriate readjustment of the coupling constant U0 will keep the resonance in the interacting susceptibility at the correct energy.
While we considered several data sets below Tc to illustrate the performance of the RPA in constructing the spin fluctuation propagator, for the calculation of the BetheSalpeter eigenvalue we restricted our analysis to temperatures slightly below and above Tc, where the FF term is not important.
Similar results for the eigenvalues were obtained from an independent method of obtaining the prefactor U 2 in the pairing interaction V, as shown in Supplementary Section 4.
However, it is instructive to examine the effect of this term on the RPA energies U0 in the superconducting state. To this end, we first expressed the real part of the bare susceptibility as Re(χ 0 ) = Re(χ 0 G ) + Re(χ 0 ) (where the second term is the anomalous contribution coming from the FF term) and defined the ratio α = Re(χ 0 F /χ 0 G ). Since 0 is unknown, so is α. However, in a d-wave superconductor, Re(χ 0 F ) and Re(χ 0 G ) add constructively at the commensurate wave-vector, so we expect this ratio to be positive at (π,π). From the pole condition for the resonance, 1 − U0Re(χ0) = 0, we see that an enhancement of the real part of the bare susceptibility reduces U0. To estimate this reduction, we computed Re 0 as in the main text (using ARPES data below Tc) and expressed the real part of the total bare susceptibility as where the ratio α was calculated using BCS normal and anomalous spectral functions [4] and a tight-binding fit for the Bi2212 dispersion [1]. We considered superconducting state spectral functions with a Lorentzian lineshape, sampled over a momentum grid of 128×128 points and energies ranging from −1.5 to 1.5 eV. We also allowed for a finite width (Γ) of the spectral function and studied the effect of varying this width. The convergence factor δ was set to 10 meV, and superconductivity was included through the d-wave gap , (2) with ∆0 = 35 meV. Since we were only interested in the effect on the coupling constant U0, we compared results only at the commensurate wave-vector (π,π). For Γ = 40 meV, α = 0.05, leading to a 5% decrease in U0. For Γ = 20 meV, α = 0.10, leading to a 10% decrease in U0.
Finally, for Γ = 0 meV, α = 0.2, leading to a 20% decrease in U0. The values of U0 reported in the main text do not take this reduction (∼ 10% for the more physically relevant case) into account because the data employed were collected near Tc.
3 Obtaining the spin susceptibility from ARPES data: the choice of RPA coupling U q Supplementary Figure 3 shows the ARPES-derived bare susceptibility χ0(q,ω = 0) for an optimally doped Bi2212 sample in the normal state 5 . The results agree qualitatively with those in ref. [1]. There is a pronounced, square-like signal enclosing the q = 0 point, as well as an incommensurate response around Q = (π,π). However, there is also significant response at q = 0, which is not seen in model calculations (Figure 1 of ref. [1]). This response was found to be an artifact of the Fermi cutoff on the spectral function.
To demonstrate the fictitious origin of the q = 0 peak, we calculated χ0 using Lorentzian spectral functions with tight binding fits for the dispersion of Bi2212. The procedure was essentially the same as in Supplementary Section 2, but setting ∆0 = 0 (for the normal state) and Γ = 20 meV. Supplementary Figures 4a-b  presented in ref. [1]. Note the absence of spectral weight at q = (0,0). The second row shows similar calculations but now the spectral function has been multiplied by the Fermi function and symmetrized, to simulate the analysis we did with experimental data. The dispersion does not cover the entire BZ because particle-hole symmetry is only assumed to hold for small binding energies (in this case, we took 80 meV about the Fermi energy). We see that the resulting bare susceptibility (panel f) has significant spectral weight at the zone center (note also the difference in color bar scales). The artifact is thus a direct consequence of not having access to the complete spectral function (including the unoccupied states), a shortcoming which can only partially be mitigated through symmetrization.
This artifact is also apparent in the interacting susceptibility, approximated via RPA as [7]: . ( If Uq is momentum-independent, χ(q,ω) possesses an anomalously strong peak near q = 0, which is not seen in INS measurements. In Supplementary Figure 5 we compare the magnetic dispersion for Imχ calculated with a momentum-independent coupling Uq = U0 from q = (π/2,π/2) to q = (3π/2,3π/2) (as is frequently reported in the literature) with the same quantity plotted from q = (0,0) to q = (2π,2π). The ARPES data employed correspond to an optimally doped Bi2212 sample in the superconducting state 6 . In the restricted momentum view, U0 yields a spin response with rich momentum structure, resembling the famous "hourglass" shape reported in the literature on other cuprates. However, the wider momentum view reveals two peaks at low energy near q = (0,0) and q = (2π,2π) three times larger than the rest of the signal (indicated by white arrows in Supplementary Figure 5b).
The spin response thus constructed displays an inordinate amount of spectral weight at momenta where experiments do not detect it. This artifact was removed by the choice of coupling function Uq defined in the main text, which suppresses the response around q = n(2π,2π). We emphasize that such a choice of Uq was not made to favor spin-fluctuationinduced superconductivity models, but was rather motivated by the fact that it rectifies the small-q intensity problem and yields a spin susceptibility in better agreement with experimental observations.

Alternative estimate of the spin-fermion coupling energy
It has been pointed in ref. [6] that the the effective electron-electron interaction should be written as , where the coupling energy U may differ from U0 due to vertex corrections. Following ref. [6], here we show how this energy can be obtained from the electron self-energy. In this section, we only report results in the important region near Tc.
Assuming the single-particle renormalizations are due to V, we can estimate the selfenergy via where n(ω) and f(ω) are the Bose and Fermi distributions, respectively, and the integral is cut off at Ω = 0.4 eV. By Kramers-Kronig relations, we can then obtain the real part of the selfenergy and the nodal Fermi velocity renormalization, where kn is the nodal Fermi momentum and 0 and are the bare and interacting nodal  [9]. Using the measured spectral functions and the spin susceptibilities calculated in the main text, we computed Z while adjusting U until we obtained the values above. For all samples, U and U0 differ by at most 11%, which is somewhat lower than the deviation reported in [6]. Using these coupling energies, the eigenvalues are still essentially equal to unity near Tc, as shown in Supplementary Table 1.
For comparison, we also included the eigenvalues obtained by assuming U = U0, as in the main text (labeled λ0).

The power method for the Bethe-Salpeter Equation
The method employed in solving the eigenvalue problem with (8) followed very closely that of Monthoux [10]. 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(kx) − cos(ky) as our starting guess.
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 Supplementary Equation (8) is essentially a convolution equation, the operations ̂Φ
were performed using Fourier techniques for increased computational speed. The numerical efficiency afforded by the convolution structure of the BSE 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.
We now present two examples of convolution equations which are one-dimensional analogues of the BSE considered in this work. We begin with an equation frequently found in signal processing applications: .
The convolution in Supplementary Equation (9) can be expressed as ]. This situation corresponds to a band pass filter which is broader than the input signal, in which case the signal is passed unaffected. Therefore, we expect ψ(x) to be an eigenfunction of the sinc-convolution operator with λ = 1. To test this, we choose a Gaussian test function ψ0(x) = e −0.1x2 for our initial guess and set λ0 = 0.5 initially. Because in this simple case the initial guess is the eigenfunction, the power method immediately converges to the eigenfunction ψ = ψ0 with λ = 1, as shown in Supplementary Figure 6.
As a second example, consider another convolution product f(ω) * ψ(ω), where ψ(ω) = e αω (α is a complex number). It can be shown that , with (12) Thus, ψ(ω) is the eigenfunction of a convolution operator (with f(ω) as a kernel) with an eigenvalue given by Supplementary Equation (12). To test our procedure, we take a gaussian "test" kernel f(x) = e −ω2 , and assume α = 0.01i (this small value was chosen so that the imaginary components in Supplementary Equation (11)

Uncertainties in the eigenvalues
The leading eigenvalue in the BSE depends on the spin-fermion coupling through the factor and the denominator of the RPA expression for the spin susceptibility. This coupling, in turn, was fixed from the empirical relation ER = 5.4kBTc.
Therefore, we estimate the uncertainty in the eigenvalue by relaxing the above relation and introducing the uncertainty ∆ER = 10 meV. This rather large variation was chosen to account not only for the INS energy resolution (< 5 meV), but also for the fact that the energy of the (π,π) peak may vary somewhat with temperature (as seen, for example, in ref. [11]).
Supplementary Figure 8 shows plots of the leading eigenvalue as a function of the coupling U0 for both thin films above (but close to) Tc. The dashed vertical lines represent the values of U0 used in the main text and for which the (π,π) peak energy was exactly given by Supplementary Equation 13. The horizontal dashed lines are the resulting eigenvalues (reported in the main text). The horizontal dimension of the shaded rectangles represents the range of U0 values for which the peak energy moves from ER +∆ER meV (lowest U0) to ER