Doping-dependence of the electron-phonon coupling in two families of bilayer superconducting cuprates

We have performed O $K$-edge RIXS studies on NBCO and Bi2212 samples with different doping levels. The NBCO samples were 100 nm films grown on SrTiO₃, with two different doping levels: a heavily underdoped sample with $T_{\mathrm{c}}=38$ K (UD38, $p\simeq 0.07$) and an optimally-doped sample with $T_{c}=92$ K (OP92, $p\simeq 0.16$). Lattice constants were $a=b=3.9~{}\textrm{\AA}$, while $c=11.714~{}\textrm{\AA}$ in the UD38 sample and $c=11.74~{}\textrm{\AA}$ in the OP one. RIXS spectra were collected at the I21 beam line of the Diamond Light Source. Energy resolution was determined to be $34$ meV by measuring the width of the elastic line on a carbon tape. We used a fixed temperature of $20$ K.

The Bi2212 samples were single crystals with three doping levels: optimally doped with $T_{c}=91$ K (OP91, $p\simeq 0.16$), overdoped with $T_{c}=82$ K (OD82, $p\simeq 0.19$), and overdoped with $T_{c}=66$ K (OD66, $p\simeq 0.21$). The RIXS measurements on Bi2212 were performed using the soft x-ray inelastic beamline (SIX) $2$-ID at the National Synchrotron Light Source II, Brookhaven National Laboratory Jarrige et al. (2018). The total experimental energy resolution was $\sim 33$ meV. We measured at $35$ K to avoid sample damage. We define our reciprocal lattice units (r.l.u.) in terms of the tetragonal unit cell with $a=b=3.8~{}\textrm{\AA}$ and $c=30.89~{}\textrm{\AA}$.

The RIXS spectra were collected with $\sigma$ incident polarization (perpendicular to the scattering plane) to maximize the charge signal. The scattering angle was fixed at 150^∘, and we used a grazing-in geometry, which by convention means negative parallel transferred momentum $\mathbf{q_{\parallel}}$, as shown in Fig. 1(g).

The first model developed to describe the EPC in RIXS experiments was the single-site model introduced in Ref. Ament et al. (2011) to describe lattice excitations generated at transition metal $L$-edges. The model was derived starting from the local Hamiltonian

where $d^{\dagger}_{\sigma}$ and $b^{\dagger}$ represent the creation operators of a spin $\sigma$ electron and local phonon mode, respectively, $n_{d}=\sum_{\sigma}d^{\dagger}_{\sigma}d^{\phantom{\dagger}}_{\sigma}$, and $M$ is the EPC constant, which is related to a dimensionless coupling parameter $g=(M/\omega_{0})^{2}$. This model is based on the simplifying assumption that both the electrons and phonon are local in the sense that they cannot travel through the lattice at any stage of the RIXS process. In this limit, the RIXS intensity can be computed exactly by use of a Lang-Firsov transformation, yielding

where $\Delta$ is the detuning of the incident energy from the resonant peak, $\Gamma$ is the half-lifetime of the intermediate state of RIXS process, and

with

being the standard Franck-Condon factor ($B_{n,m}$ is understood to mean $B_{\mathrm{max}(n,m),\mathrm{min}(n,m)}$).

Reference Geondzhian and Gilmore (2020) recently generalized the single-site model to the case where the localized electron couples to more than one type of phonon. This generalization allows each phonon to be coherently excited and is relevant to the cuprates where we naturally expect coupling to different phonon branches. Here, we consider the case with two types of phonon modes (e.g. the bond-buckling and bond-stretching modes), again described by the local Hamiltonian Geondzhian and Gilmore (2020)

The exact RIXS intensity can still be computed, yielding

where

Here, $C_{\mathrm{scale}}$ is an overall amplitude factor, independent of momentum and detuning and applies to all phonon branches, which we use to scale overall intensity to match experiment. The coefficients $D$ in Eq. (7) are given by a product of Franck-Condon factors

where $g_{\lambda}=(M_{\lambda}/\omega_{\lambda})^{2}$.

We stress that both models neglect all electron mobility involved in the RIXS scattering process. Nevertheless, relatively simple expressions for the RIXS intensity [Eqs. (2) and (6)] are obtained from these approximations that are useful for fitting experimental data. The parameters setting the intensity of phonon excitations in both models are the EPC strength of the different branches $g_{\lambda}$, the detuning energy $\Delta$, and the core-hole lifetime parameter $\Gamma$. Throughout, we treat $g_{\lambda}$ as a fitting parameter and set $\Delta$ based on experimental conditions. We will also fix $\Gamma=0.15$ eV Lee et al. (2013); Johnston et al. (2016); Braicovich et al. (2020) for all the samples investigated, since the intrinsic lifetime of the state is mostly determined by Auger recombination rate and so independent of the material studied.

A recent theoretical study in the single carrier limit Bieniasz et al. (2021) found that the local model for the RIXS cross-section produces stronger phonon features in comparison to a fully itinerant model. Moreover, the discrepancy between the models becomes more pronounced as the strength of the localizing core-hole potential is reduced. This result suggests that fitting a localized model to phonon spectra may underestimate the strength of the EPC. In the single carrier limit, electron mobility in the intermediate state can also introduce a momentum dependence in the intensity of the phonon excitations that is unrelated to the momentum dependence of the EPC constant $M({\bf k},{\bf q})$. With these caveats in mind, we proceed by fitting the experimental data with the localized models with the assumption that doing so will allow us to extract general trends in the data. For these reasons, the absolute values of the extracted couplings obtained here should not be equated to the value entering a microscopic Hamiltonian.

The X-ray absorption spectra (XAS) of our Bi2212 and NBCO samples are displayed in Fig. 1(a) and Fig. 1(b), respectively. It is easy to identify the Zhang-Rice singlet (ZRS) pre-peak in the O $K$-edge XAS of Bi2212, which originates from the hybridization of the in-plane O $2p_{\sigma}$ and Cu $3d_{x^{2}-y^{2}}$ orbitals Hawthorn et al. (2011). Both the energy and amplitude of the ZRS peak change with doping as expected. On the other hand, we observe two different structures in NBCO whose intensity increases with doping: a main peak at 528.2 eV, and another peak at a lower energy ($\sim 527.6$ eV). The former one is associated with doped holes in the CuO₂ planes, while the latter originates from the CuO₃ chains, which act as a charge reservoir in NBCO Hawthorn et al. (2011). Our UD38 sample has short chains in both in-plane directions due to the tetragonal structure, while the OP92 sample has fully developed long chains, and due to the twin structure, they exist in both in-plane directions. The comparison between UD38 and OP92 sample also nicely highlights the decrease in spectral weight of the upper Hubbard band (UHB) peak around 529.5 eV.

To investigate the electron-phonon interaction, we have acquired RIXS spectra at different incident X-ray energies across the doping-related resonances. Figs. 1(c) and 1(d) present an overview of the RIXS spectra for Bi2212 OP91 and NBCO OP92, respectively. The incident energy of each scan is indicated by the tick marks with the same color in panels (a) and (b) for the two compounds. As highlighted in Fig. 1(e), we recognize four features in the spectra: the elastic peak, the phonons (and their overtones) below 0.2 eV, the broad charge excitations peaked around 0.8 eV, and some higher-energy excitations above 1.5 eV from $dd$ and charge-transfer (CT) excitations. Moreover, we also identify some spectral weight around 0.5 eV, which could be due to broad bimagnon excitations Bisogni et al. (2012). To highlight the phonon contribution to the spectra and compare intensities between different samples at resonant energy, we removed the elastic peak and normalized the spectra to their respective dd+CT excitations as summarized in Fig. 1(f). The phonon excitations have a larger spectral weight in NBCO and the exciations in Bi2212 OD82 ($p_{c}$=0.19) have a stronger intensity than those observed in the OP91 and OD66 samples.

To precisely determine the resonant energy of the phonon signal (i.e., the energy at which the phonon intensity is maximum), we have extracted their integrated spectral weight. This value was determined by removing the elastic peak and charge excitations from RIXS spectra, and then integrating the infrared region of the RIXS spectra ($0-300$ meV). As a reference, we have also extracted the spectral weight of CT excitations by integrating the spectra in a 3 eV window centered around their maximum value. The integrated spectral weight of the phonon and CT excitations are plotted in Figs. 2(a) and 2(b) for Bi2212 OP91 and NBCO OP92, respectively, normalized to their maximum value and superimposed over the corresponding XAS spectra. We notice that the one-dimensional breathing mode of the chains has a stronger EPC strength than its two-dimensional counterpart, which is evident from the higher intensity of its first overtone.

The intensity of CT excitations follows closely the XAS spectra in both samples, confirming the stability of our incident energy and the independence of the lifetime of the intermediate state. However, the phonon has a large shift of $200-250$ meV relative to the maximum of XAS, both in Bi2212 and NBCO. The presence and magnitude of this energy shift are puzzling. A small displacement, of the order of the phonon energy ($\approx 50-100$ meV) between the maximum of the XAS and the maximum of phonon intensity has been observed before Feng et al. (2020); Peng et al. (2020) and can be accounted for by the localized Lang-Firsov model with a strong EPC and a long core-hole lifetime. The large value observed here, however, requires a different explanation. As suggested in Ref. Geondzhian and Gilmore (2018), a possible origin could be a change in the curvature of the potential energy surface between the ground and intermediate state of the RIXS process.

We initially attempted to use the energy detuning method to extract the value of the EPC for the buckling and breathing branches, which has been described in detail in Ref. Braicovich et al. (2020) and has been successfully applied at the Cu $L_{3}$ edge Braicovich et al. (2020); Rossi et al. (2019). This approach is based on the single-site model presented in Sec. II.2 and consists in acquiring scans at different incident energies, at and below the resonant peak in the XAS. The EPC of the different phonon branches can then be extracted by fitting the intensities of the different phonon branches (and of their higher harmonics) with Eq. (6). To reduce the number of parameters, we fit all of the spectra for a given sample at the same time, using a single overall scale factor and fixing $\Gamma=0.15$ eV Lee et al. (2013); Johnston et al. (2016). In this way, the only free parameters are the EPC strengths $M_{\lambda}$ of the two branches and their energies $\Omega_{\lambda}$. We have also subtracted the elastic peaks and considered only the inelastic parts, since the elastic peaks can carry additional contribution from other types of scattering.

As highlighted in Sec. III.1, one difficulty is that the maximum of the phonon signal is displaced from the maximum of XAS spectra by $\sim 200-250$ meV. While the employed model can account for a small displacement between the two energies (of the order of $\omega_{\lambda}$) in the strong coupling limit, the value of $M_{\lambda}$ that is needed is too large and would produce phonon harmonics inconsistent with the RIXS spectra. Therefore, we have decided to artificially put the zero of the detuning in the model at the phonon resonance.

To test the validity of the detuning approach, we have used two different datasets: the chain resonance of NBCO OP92 and the plane resonance of Bi2212 OP91. We chose the chain resonance for NBCO for two reasons. First, it selects the one-dimensional CuO₃ chains as opposed to the two-dimesional CuO₂ planes of Bi2212, which allows us verify whether there are differences related to dimensionality. Secondly, the chain resonance energy ($527.6$ eV) is lower than that of the plane ($528.2$ eV). This difference means that the detuning energies of the chain are not obscured by additional contributions from the plane resonance, while the detuning energies of the plane overlap with the chains.

The spectra with the corresponding fittings are presented in Fig. 3. The solid orange and purple lines are the bond-buckling and bond-stretching phonons including the overtones, respectively, while the green shaded area represents the mixed coherent multi-phonons. However, it is evident from the spectra that the detuning method does not reproduce the measured intensities for all incident photon energies. While the agreement is quite good close to the resonance energy for both the cases, the calculated spectra far from the resonance significantly underestimate (by more than $\sim 50$%) the intensity of the phonon branches. This failure of the detuning method at the ligand $K$ edge is not completely unexpected given the very different nature of the intermediate state between the Cu $L_{3}$- and the O $K$-edges, and of the localized assumptions underpinning the Lang-Firsov models. Indeed, at the Cu $L_{3}$-edge, the excited electron is more strongly bound to the core hole in comparison to the O $K$-edge. This means that the excited electron is more free to propagate around the lattice during its lifetime in the latter case. Reference Bieniasz et al. (2021) recently developed a model for an intermediate state consisting of an itinerant electron interacting with phonons and the core hole in an otherwise empty band. One of the main results derived there is that the detuning curves are quantitatively different from the localized model; they have a weaker dependence on the EPC when the core-hole potentials are weak and they generally decay faster as a function of the detuning energy in comparison to the localized model. We observe a small difference between the decay of the intensity of the two phonon branches with detuning, as evident from Fig. 3. This observation suggests that the real situation close to the Fermi energy in cuprates is intermediate between the two approaches.

Determining the EPC as a function of doping is crucial for understanding its potential contribution to $d$-wave pairing. For this purpose, we fit the RIXS spectrum at the resonant energy using Eq. (6), after subtracting the elastic, charge, and CT excitations from the data, as shown in Fig. 4. As discussed above, the use of Eq. (6) is likely to underestimate the overall strength of the EPC, particularly at the O $K$-edge, where electron itinerancy in the RIXS intermediate state can be more pronounced than at the Cu $L$-edge Bieniasz et al. (2021). We nevertheless proceed with Eq. (6) in a comparative analysis to identify trends in the data, both as a function of doping and between Bi2212 and NBCO. Presumably, the degree of itinerancy in these systems is comparable for similar doping levels. Thus, the analysis for both samples should be subject to the same systematic error when applying the local model to the data. One should keep in mind that the absolute magnitudes of EPC constants extracted in this analysis cannot be mapped directly onto the $M_{\lambda}({\bf q})$ entering in any microscopic Hamiltonian.

To fit the data, we consider the coherent contribution from two phonon modes – one from the bond-buckling branch ($\omega_{1}\simeq 40$ meV) and one from the bond-stretching branch ($\omega_{2}\simeq 70$ meV) – while treating the coupling strengths $M_{\lambda}$ as fitting parameters. We also introduce a prefactor $C_{\mathrm{scale}}$, which multiplies the entire spectrum by a constant amount. We are able to obtain a consistent fit to the doping-dependent RIXS spectra using this approach. The doping evolution of the extracted amplitude $C_{\mathrm{scale}}$ for the Bi2212 sample are plotted in Figs. 4(c-e) for different momentum transfers, while the corresponding $M_{\lambda}$ values are displayed in Figs. 4(f-h). For comparison, we have also overlaid the doping dependence of the extracted $M_{\lambda}$ in NBCO at $H=-0.25$ r.l.u. in Fig. 4(f). Notably, the values of the extracted coupling in Bi2212 and NBCO are almost identical at optimal doping, which provides additional post-hoc justification for our comparative approach.

For momentum transfer $H=-0.25$ r.l.u., we find that the EPC of the lower-energy buckling mode is stronger than the EPC of the breathing mode in both samples and for all the measured doping levels. This observation is consistent with the previous Cu $L$-edge RIXS measurements on NBCO UD38 sample for small momentum transfers Braicovich et al. (2020) and prior ARPES experiments that infer the strongest coupling to the bond-buckling phonon modes Cuk et al. (2004, 2005); Johnston et al. (2010, 2012). As the hole concentration increases, we find that the extracted $M$ values for both branches decrease monotonically at $H=-0.25$ r.l.u.. This behavior can be understood as being a consequence of the improved electronic screening as the system becomes more metallic Johnston et al. (2010). The situation is quite different for smaller momenta. For $H=-0.20$ r.l.u., the inferred EPC constants are very weakly dependent on doping, while at $H=-0.15$ r.l.u. the doping dependence even becomes non-monotonic, with both branches becoming more strongly coupled near $p_{c}=0.19$.

The doping dependence summarized here in both compounds is consistent with the poor screening scenario discussed in Ref. Johnston et al. (2010). There, it was shown that the poor conductivity perpendicular to the CuO₂ plane results in a nonuniform screening of the electron-phonon interaction whereby it was well screened for large momentum transfers but poorly screened at small ${\bf q}$. The suppression of large ${\bf q}$ coupling relative to small ${\bf q}$ will increase the EPC in the $d$-wave channel, suggesting that the contribution of these modes to pairing increases as the doping increases and may become largest at optimal doping.

Both optimally doped NBCO Ghiringhelli et al. (2012); Blanco-Canosa et al. (2014) and Bi2212 Hashimoto et al. (2014); da Silva Neto et al. (2014); Lee et al. (2021) present a CDW ordering with $q_{\text{CDW}}\simeq 0.29$ r.l.u., though the CDW signals become weaker with respect to underdoped samples Blanco-Canosa et al. (2014). Previous RIXS studies at Cu $L_{3}$-edge and O $K$-edge on doped cuprates have found an anomalous enhancement of phonon intensity near the CDW wavevector, which is ascribed to the presence of dispersive CDW excitations that overlap and interfere with the phonon excitations Chaix et al. (2017); Li et al. (2020). We have, therefore, searched for possible enhancements of EPC of buckling and breathing branches due to the presence of charge order.

Employing again the generalized Lang-Firsov model [Eq. (7)] that accounts for coupling to both modes, we have extracted the evolution of the EPC values and the overall scaling prefactor along the $(H,0)$ direction. Figure LABEL:fig:momentumdep shows the momentum dependence of the phonon amplitudes ($C_{\mathrm{scale}}$) and $M_{\lambda}$ values in optimally doped Bi2212 and NBCO. Evidently, the phonon amplitudes increase monotonically on approaching the $q_{\text{CDW}}$. At the same time, the EPC does not show a clear trend with $M_{1}$ larger than $M_{2}$ for all our measured momenta. These results indicates that while the phonon excitations increase in intensity in proximity of the CDW, there is no clear corresponding enhancement in the relative intensity of the phonon overtones. This behavior can be reconciled if there is an overlapping contribution to the elastic line intensity due to the CDW excitations and no corresponding increase in the EPC close to ${\bf q}_{\mathrm{CDW}}$. One should, therefore, be careful in interpreting directly the phonon intensity as the EPC strength at O $K$-edge due to the longer intermediate state. The explanation of this discrepancy remains a topic for further investigation.