Determining the Electron-Phonon Coupling in Superconducting Cuprates by Resonant Inelastic X-ray Scattering: Methods and Results on Nd$_{1+x}$Ba$_{2-x}$Cu$_3$O$_{7-\delta}$

The coupling between lattice vibration quanta and valence electrons can induce charge density modulations and decisively influence the transport properties of materials, e.g. leading to conventional superconductivity. In high critical temperature superconductors, where electronic correlation is the main actor, the actual role of electron-phonon coupling (EPC) is being intensely debated theoretically and investigated experimentally. We present an in-depth study of how the EPC strength can be obtained directly from resonant inelastic x-ray scattering (RIXS) data through the theoretical approach derived by Ament et al. [EPL 95, 27008 (2011)]. The role of the model parameters (e.g. phonon energy $\omega_0$, intermediate state lifetime $1/\Gamma$, EPC matrix element $M$, and detuning energy $\Omega$) is thoroughly analyzed, providing general relations among them that can be used to make quantitative estimates of the dimensionless EPC $g = (M/\omega_0)^2$ without detailed microscopic modeling. We then apply these methods to very high resolution Cu $L_3$ edge RIXS spectra of three Nd$_{1+x}$Ba$_{2-x}$Cu$_3$O$_{7-\delta}$ films. For the insulating antiferromagnetic parent compound the value of $M$ as a function of the in-plane momentum transfer is obtained for Cu-O bond-stretching (breathing) and bond-bending (buckling) phonon branches. For the underdoped and the nearly optimally doped samples, the effects of Coulomb screening and of charge-density-wave correlations on $M$ are assessed. In light of the anticipated further improvements of the RIXS experimental resolution, this work provides a solid framework for an exhaustive investigation of the EPC in cuprates and other quantum materials.

The coupling between lattice vibration quanta and valence electrons can induce charge density modulations and decisively influence the transport properties of materials, e.g. leading to conventional superconductivity. In high critical temperature superconductors, where electronic correlation is the main actor, the actual role of electron-phonon coupling (EPC) is being intensely debated theoretically and investigated experimentally. We present an in-depth study of how the EPC strength can be obtained directly from resonant inelastic x-ray scattering (RIXS) data through the theoretical approach derived by Ament et al. [EPL 95, 27008 (2011)]. The role of the model parameters (e.g. phonon energy ω0, intermediate state lifetime 1/Γ, EPC matrix element M , and detuning energy Ω) is thoroughly analyzed, providing general relations among them that can be used to make quantitative estimates of the dimensionless EPC g = (M/ω0) 2 without detailed microscopic modeling. We then apply these methods to very high resolution Cu L3 edge RIXS spectra of three Nd1+xBa2−xCu3O 7−δ films. For the insulating antiferromagnetic parent compound the value of M as a function of the in-plane momentum transfer is obtained for Cu-O bond-stretching (breathing) and bond-bending (buckling) phonon branches. For the underdoped and the nearly optimally doped samples, the effects of Coulomb screening and of charge-density-wave correlations on M are assessed. In light of the anticipated further improvements of the RIXS experimental resolution, this work provides a solid framework for an exhaustive investigation of the EPC in cuprates and other quantum materials.

I. INTRODUCTION
The role of the electron-phonon coupling (EPC) in the high critical temperature superconducting cuprates is still an open problem deserving further research. Indeed, even if pairing is not of the phonon-mediated BCS type [1], the question of a potential role for the EPC remains extremely interesting. For instance, it has been around the CDW wave vector. These examples, including results from the recent literature, provide new clues towards understanding the role of the EPC in the cuprates, on top of the long debated polaronic behavior of carriers in this class of materials [20][21][22][23].
Clarifying the role of the EPC in cuprates and in other strongly correlated systems is a formidable task because the EPC acts on top of strong electron-electron interactions, giving rise to a correlated electron liquid where phononic and electronic degrees of freedom are highly entangled. This difficulty is reflected in broadly discrepant estimates of the EPC, not only between theory [24][25][26][27][28] and experiment [15,17,[29][30][31][32] but also among different theoretical approaches [2, [24][25][26]. For example, the Cu-O bond-stretching modes have been proposed to mediate both an effective attractive [6,33] and, more recently, repulsive [7,34,35] interaction. In this context, experimental verification of the phenomenology, possibly with momentum resolution, is crucial. This goal calls for the improvement of the existing techniques and for the introduction of new methods. Among traditional measurements, bulk-sensitive neutron scattering [36,37] suffers from sensitivity limitations as it requires massive homogeneous samples, whereas surface-sensitive angleresolved photoemission spectroscopy [31,38] and scanning tunneling microscopy [39] need high quality surfaces that are not always available. Moreover, the comparison with bulk sensitive results is a nontrivial task. In this field, soft x-ray resonant inelastic x-ray scattering (RIXS) is likely to offer an important contribution, made possible by the recent progress in instrumentation [40]. For instance, state-of-the-art RIXS can be performed with a combined energy resolution of ∼ 30 -40 meV at the Cu L 3 edge, which facilitates studies of high-energy phonons ( > ∼ 30 meV). For detailed investigations, a further increase of the resolution is needed, which is a realistic prospect in the next three to five years. What is needed in this context is a well defined procedure for extracting quantitative information about the EPC from the available RIXS data, preferably in a model independent way.
The proposal of determining the EPC from RIXS data was first introduced by Ament et al. [41] (for a detailed account see Ref. 42). This idea has also been reinforced by the theoretical treatment by Devereaux et al. [43], based on the lowest order Feynman diagrams. These calculations on single layer Bi-based cuprates show very distinct RIXS signals for the main phonons excited at the Cu L 3 edge. The relationship between the EPC and the RIXS signal has been expanded recently by Geondzhian and Gilmore [44] with a model conceptually similar to Ref. 41. Due to the limited energy resolution, RIXS experiments on lattice excitations were mainly carried out at low incident energies, i.e. the O K edge in cuprates and other functional oxides and L edges of systems containing elements at the beginning of the 3d series. A convenient option are the quasi-one-dimensional compounds [45,46], having in general more pronounced phonon peaks and van Hove singularities. A recent remarkable work has also reported the first application to nSrIrO 3 /mSrTiO 3 multilayers [47]. For the low-Z 3d oxides, data on Ti L edge are also available [48,49].
From a conceptual point of view, state-of-the-art RIXS research on the EPC has proceeded along two main lines: i) the search for a simplified universal scheme allowing experimentalists to address a variety of cases without detailed calculations; ii) the development of more specific treatments that exploit numerical computation. While both approaches are useful and cross fertilize each other, in the following we focus on the first approach. Specifically, we introduce innovative ways of using the theory by Ament et al. [41] and identify some scaling laws with a wide range of application. The proposed methods are then employed on the cuprate Nd 1+x Ba 2−x Cu 3 O 7−δ (NBCO). We examine the EPC in the antiferromagnetic compound and the effect of charge order [50][51][52][53][54] on the EPC at different doping levels. From the RIXS data, we extract the EPC strength and discuss its doping dependence in relation to the presence of CDW correlations in the system. We selected a compound of the so-called "123" family of cuprates since the CDW signal is one of the strongest in this case [53]. Our work also paves the way for more advanced experiments in the future.
The article is organized as follows: Section II describes the samples and the experimental setup. Section III provides the theoretical tools used to extract the EPC strength from RIXS data, with subsections on three different methods. Section IV presents the experimental RIXS results for both the antiferromagnetic parent compound and the doped samples. Finally, Section V provides a discussion of our results.

A. The NBCO samples
The crystal structure of Nd 1+x Ba 2−x Cu 3 O 7−δ (NBCO) is displayed in Fig. 1(a) for the undoped antiferromagnetic parent compound (x = 0 and δ = 1). NBCO belongs to the so-called "123" family of high T c superconductors and is isostructural to YBa 2 Cu 3 O 7−δ (YBCO). Of primary importance in the centrosymmetric unit cell is the presence of the CuO 2 bilayer where superconductivity occurs. The hole density per CuO 2 plane can be changed by modifying the oxygen content δ (similarly to YBCO [55]) or by introducing excess Nd at the Ba sites (x). In the latter case, 0.5 holes are removed from the CuO 2 planes for each excess Nd [56,57]. At optimal doping, T c reaches 93 K. NdBa 2 Cu 3 O 6 has a tetragonal crystal structure; however, the system becomes orthorhombic upon doping. For the purpose of the present study, we neglect orthorhombicity and adopt a tetragonal description. The CuO 2 planes are characterized by a buckling of the Cu-O bonds, which depends on the overall stoichiometry. As a general rule, the buckling angle increases from underdoped to overdoped NBCO [58].
Epitaxial NBCO films were deposited on a (001)oriented SrTiO 3 substrate with an almost perfect inplane matching of the a and b lattice parameters and the buckling practically unchanged with respect to the corresponding bulk structure [56,57]. The thickness of the NBCO films is around 150 nm, hence they can be considered infinite for the scattering angles used in our experiment (where the incidence angle with respect to the surface ranged from 25 • to 65 • ).
Three NBCO samples were investigated in the present work: • NBCO-AF: the sample is antiferromagnetic with x ≈ 0, δ ≈ 0.9, and with a Cu-O-Cu buckling angle of 6.35 • [59, 60]; • NBCO-UD: the sample is underdoped with T c = 63 K, hole doping of 0.11, x ≈ 0.2, δ ≈ 0, and a buckling angle of 6.74 • [61]; • NBCO-OP: the sample is very close to optimal doping and has T c = 90 K, x ≈ 0, δ ≈ 0, and a buckling angle of 7.75 • [59]. The slightly lower T c compared to optimal doping is due to a slight overdoping, which has been identified by measuring the c lattice parameter.

B. Experimental layout and energy resolution
The measurements were performed at beamline ID32 of the European Synchrotron (ESRF, France), which is equipped with the new spectrometer ERIXS (European RIXS) [40]. The layout of the experiment is sketched in Fig. 1(b). The sample is mounted on the in-vacuum diffractometer with the ab planes perpendicular to the horizontal scattering plane, defined by the incoming and outgoing photon wave vectors (q in and q out , respectively). The scattering angle 2θ was fixed to 149.5 • . Most of the measurements were performed with a combined (beamline + spectrometer) energy resolution of ≈ 40 meV at the Cu L 3 edge (∼ 930 eV). The contribution to the resolution from the incident beam is ≈ 25 meV. One spectrum was collected with high statistics and high energy resolution (32 meV) to be used as benchmark.
The momentum transferred to the sample Q = q in − q out and projected on the basal ab plane is the relevant observable in dispersion studies of two-dimensional cuprates. The projected in-plane momentum Q is scanned in the Γ-X direction of the Brillouin zone (BZ) by rotating the sample around the b axis, i.e. Q = 2π a h, 0 . The rotation direction was chosen in order to minimize the incidence angle to the surface, which reduces self-absorption effects, one of the experimental difficulties of RIXS. The incident beam polarization was perpendicular to the scattering plane [σ incident polarization, red arrow in Fig. 1(b)]. Thus, the electric field vector lay in the basal plane, simplifying the problem (see Sec. III A). The sample was kept at a temperature of 20 K during the measurements. We did not investigate small momentum transfers |Q | < 0.1 r.l.u. because the phonon and magnon signals begin to overlap in this region, making it difficult to identify the phonon contribution. Figure 1(c) displays an example of RIXS spectrum of NBCO-AF over an extended energy range. The inset provides a close-up of the region of interest to the present work. Note that a phonon peak around 70 meV is clearly visible in the data.

A. Theoretical background
The present work is based on the theory by Ament et al. [41,42], which was derived for the case of dispersionless electrons coupled to local phonons. The approximation of Einstein phonons is appropriate for most of the medium-and high-energy optical branches in cuprates, which are generally found to be weakly dispersing [36,37]. We do not introduce any change in the theory here, instead we propose an innovative way of using it. For completeness, we now summarize the main theoretical results of this approach. The following equations are rather intuitive in spite of their complexity at first sight. For a rigorous derivation, the reader is referred to Ref. 42. In the case of non-dispersing phonons -which is the situation we are interested in -the starting Hamiltonian is where ω 0 is the phonon energy, M is the EPC matrix element, and b † i (b i ) and d † i (d i ) are the creation (annihilation) operators for phonons and electrons, respectively, at site i. We have neglected the spin index for brevity. Equation 1 can be diagonalized exactly using a convenient canonical transformation H = e S He −S , This same transformation can be used to derive an exact expression for the scattering amplitude: (2) Here, T el is the polarization-dependent atomic elastic scattering factor; z = Ω+iΓ is a complex value with a real term corresponding to the detuning energy Ω (i.e. the difference between the incoming photon energy and the resonance energy) and an imaginary term Γ, which is the inverse core hole lifetime. Here, initial state |n 0 i denotes the ground state of the Hamiltonian, which corresponds to the ground state of a shifted Harmonic oscillator.
By introducing the Franck-Condon (FC) coefficients the cross section can be written as: where g = (M/ω 0 ) 2 is a dimensionless coupling constant and the index n identifies the number of excited phonons in the final state f .

B. The genesis of the phonon during the RIXS process
A phonon is excited in the intermediate state of the RIXS scattering via the process sketched in Fig. 2. Upon photon absorption at the Cu L 3 resonance, an excited state is created having a 2p core hole and an extra electron in the lowest unoccupied level compatible with the selection rules. This intermediate state has excitonic character as already stated in Ref. 41; however, from the point of view of the surrounding oxygen ions, the extra electron localized in the excited site strongly screens the core hole. The additional 3d electron also locally perturbs the charge density and repels the nearby oxygen ions. The interaction between the new charge distribution and the oxygen ions is mainly of quadrupolar nature, since the monopole is almost unchanged due to the excitonic character of the intermediate state (neutral excitation) and the dipole vanishes by symmetry of the 3d states and of the crystal field. Note that this mechanism is typical of resonant x-ray scattering and it is not present in nonresonant x-ray scattering. In a sense, we can con-sider RIXS as a way of introducing a probing charge to measure the EPC (i.e. the value of g) while maintaining charge neutrality. In principle, the value of g measured with RIXS is different from the value g t involved in the transport measurements, as stressed in Ref. 44, because of the core hole effect and the symmetry of the intermediate state. Thus, the strong screening helps to reduce the difference between g and g t . The core hole's role in localizing the excited electron is crucial to this process; it is the excitonic nature of the intermediate state that makes g and g t close. The present accuracy of RIXS measurements does not allow a detailed discussion on the relationship between g and g t .
We consider now the effect of the time scale dictated by the core hole lifetime. It is intuitive that the phonon signal increases with the EPC matrix element M , and this observation allows one to interpret a RIXS intensity map as an EPC landscape (although this is not always the case as discussed above). The possibility of mapping the EPC in momentum space is the main reason for studying phonons with RIXS. It is also intuitive that the phonon signal decreases with the core hole lifetime (1/Γ), as shorter lived core hole states cannot provide sufficent time for the lattice to respond. The logical consequence is that the phonon excitation process should exclusively depend on the parameter M/Γ. Upon absorption, an electron is promoted into an empty valence state and a core inner hole (in green) is created in the Cu ion. Due to the attraction between the core hole and the electron in a screening orbital (in blue) the state has excitonic character. Seen from the oxygen ion, the core hole is well screened while the electron pushes/pulls the oxygen ions (red spheres). The black solid line schematically shows the radial part of the wave function of the photoexcited electron. For simplicity, we assumed spherical symmetry in the drawing. The correct symmetry is discussed in the text.
is confirmed by numerical calculations summarized in Fig. 3. On the left axis, we plot the intensity of the one-phonon excitation I 1 , normalized to (ω 0 /Γ) 2 , for several values of the dimensionless coupling constant g = (M/ω 0 ) 2 . It is evident that all curves corresponding to different g values collapse on the same "universal" curve when (M/Γ) 2 < ∼ 2. Above this value, the curves begin diverging so that the behavior is no longer universal. Similar considerations apply to the ratio I 2 /I 1 of the intensities of the one-phonon excitation and its first overtone I 2 (right axis). We note that for an Einstein phonon, the two-phonon excitation has twice the energy of a single phonon, which would not be the case for dispersing phonon branches. In this simplified picture, it is interesting to note that I 2 /I 1 vs. (M/Γ) 2 is independent of g for (M/Γ) 2 < ∼ 1.5. At the Cu L 3 resonance, Γ ≈ 250 -280 meV [62] so the universality range extends up to about M ≈ 350 meV. This value is much larger than most estimates for the EPC constants appearing in the literature. Thus we can safely say that the intensity of the high-energy phonon excitations in cuprates universally scales with (M/Γ) 2 .
For small interactions, the intensity I 2 is directly proportional to (I 1 ) 2 by definition [41], so that I 2 /I 1 coincides with I 1 apart from an overall factor. As a matter of fact, the inset of Fig. 3 shows that the two curves are superimposed at low interactions after rescaling. The region of small interaction extends up to (M/Γ) 2 ≈ 0.12. Note that the condition of linearity is even more restrictive, as it is evident from the Figure. The universal plots can be exploited in several ways  Figure 3. The universal plots of the phonon excitation intensities based on rescaled dimensionless variables. M is the absolute value of the EPC matrix element, Γ is the intrinsic width of the Cu L3 resonance so that 1/Γ is the core hole lifetime, and g = (M/ω0) 2 is a dimensionless coupling constant. Left axis: The intensity of the one-phonon excitation I1 rescaled by (ω0/Γ) 2 and plotted as a function of (M/Γ) 2 . The different lines corresponding to different values of g superimpose almost perfectly, so that the curve is universal. Right axis: The ratio between the two-and one-phonon excitation intensities as a function of (M/Γ) 2 . The behavior is also universal for (M/Γ) 2 < ∼ 1.5. The inset shows that in the limit of small coupling, the universal curves of I1 (solid black line) and I2/I1 (dashed red line) share the same behavior as a function of (M/Γ) 2 apart from an overall factor. to extract the EPC from experimental data. A critical comparison between the different methods is lacking in the literature and is given briefly in the following.

C. The use of I2/I1
When the phonon peak and its overtone can be detected in the experimental RIXS spectra, one can directly determine I 2 /I 1 and use the universal plot of Fig. 3 to obtain (M/Γ) 2 . This approach has been used in early literature at the O K edge [45][46][47] and at the Ti L 3 edge [48,49] with some uncertainty due to the modest resolving power. In order to make the application more transparent, we present a simple geometrical construction allowing at a glance to capture the interplay between the different parameters (see Fig. 4). Let us suppose we can determine experimentally the intensity ratio I 2 /I 1 , the phonon energy ω 0 , and the core hole lifetime 1/Γ. The value of I 2 /I 1 corresponds to the point P 1 on the universal curve (left panel) and thus to a specific value of M/Γ. The horizontal line through P 1 intersects the value of (ω 0 /Γ) 2 , which is the entrance value of the right panel. This defines the point P 2 . The straight line from the origin through P 2 defines the angle α whose tangent is (M/Γ) 2 /(ω 0 /Γ) 2 = g. The I 2 /I 1 approach does not require the knowledge of the absolute efficiency of the instrument since it is based on a ratio of intensities. An important issue is the momentum dependence of the EPC. In our experiment we used σ incident polarization. As shown by measurements with polarization analysis of the scattered beam [63], the σσ channel (with incident and scattered σ polarization) has an overwhelming intensity with respect to the σπ channel. We note that the ionic cross section for Cu L 3 resonant elastic scattering in the σσ channel is isotropic; moreover the σπ cross section has negligible angular variation in our range [63,64]. For these reasons, we do not expect any sizable angular variation of the phonon intensity due to the resonant form factor. Thus the measured change in the phonon intensity as a function of the incidence angle is due to a genuine variation of the EPC influencing I 1 , I 2 and, as a consequence, their ratio. In this way we establish a link between Q and the EPC.
Given the reasons above, the I 2 /I 1 method may seem ideal. Unfortunately, this approach comes with a very serious drawback: It is very difficult and often impossible to experimentally identify the two-phonon spectral feature and to determine its intensity. As a matter of fact, at the Cu L 3 edge this is only possible with great difficulty in parent compounds, as demonstrated below. On the contrary, in the doped cuprates the broadening of the features and the presence of a continuum prevent the detection of the overtones in the majority of cases.

D. The use of I1 on resonance
The universal curve of Fig. 3 can also be used to recover (M/Γ) 2 from the intensity of the single phonon excitation I 1 . The geometrical construction shown in Fig. 4 is formally the same, but with the ratio I 2 /I 1 substituted by I 1 in the left panel. However, the use of this method requires the measurement of the absolute value of I 1 and RIXS cross sections are not typically measured in an absolute way. Hence, this approach can be exploited only if the intensity is calibrated at a minimum of one point by determining M with another method. On the other hand, the advantage is that the relative intensity I 1 at different momentum transfers can be measured very accurately. The diagram in Fig. 3 highlights also the nonproportionality of the I 1 intensity to the EPC. In particular, the broad maximum of I 1 is due to the increase of I 2 spoiling the intensity I 1 at large values of the electronlattice interaction. The existence of a maximum means that the function in the left panel of Fig. 4 is no longer single valued, creating some ambiguities. One might be able to resolve this issue by considering the nature of the physical problem, but this may not always be possible.

E. The use of I1 upon detuning below threshold
This method was introduced recently by some of the authors in a concise way [65]. We add here some relevant information. The method is based on two observations: i) upon detuning below threshold, i.e. using an incident photon energy slightly lower than the absorption resonance, the phonon signal evolves differently with respect to the other features in the same RIXS spectrum; ii) the difference in the phonon behavior depends on the strength of the EPC. Therefore, it is possible to recover the EPC from a suitable set of spectra measured as a function of the detuning. An inspection of the scattering amplitude (Eq. 2) helps one understand why this is the case. The denominator is the sum of the detuning energy Ω = ω − ω res , of iΓ and of M 2 /ω 0 . Thus, detuning has a larger effect if M is small. Indeed, there is a kind of trade-off between detuning and M 2 /ω 0 .
An example of the effect of detuning on the RIXS spectrum of NBCO-AF is shown in Fig. 5, where we compare the spectrum measured at h = −0.4 r.l.u. and with an incident photon energy tuned at the Cu L 3 resonance (black squares) and detuned by −0.5 eV (red circles). It is evident that the phonon and the bimagnon become weaker upon detuning; their processes of excitation are "slow" and the intensity depends on the time duration set by the core hole decay. It is sufficient to fit the measured detuning effect to the cross section to obtain g. More precisely, we define the detuning curve as the phonon intensity as a function of the excitation energy with normalization to the value at resonance. An overview of the theoretical detuning curves is given in Fig. 6, where each panel is identified by the value of Γ/ω 0 and contains the detuning curves for different values of g. Since Γ/ω 0 is known for a given absorption edge, by comparing the experimental data to this set of theoretical curves one can identify the appropriate detuning curve, and thus the corresponding value of g. The sensitivity of the method is clearly decreasing at higher values of Γ/ω 0 , where the approach is useful only in the cases of large g. A very convenient way of handling the data is presented in Fig. 7, which reports the width (in units of ω 0 ) of the detuning curve, defined as the detuning value of the point at half height. The width is displayed as a function of √ g, with Γ/ω 0 as a parameter. The lines are remarkably linear in a wide range and at lower values of Γ/ω 0 they are very close to each other, suggesting an approximate scaling law.
The method based on detuning is advantageous because it uses the value of I 1 , which is easily measurable, normalized to its value at resonance. In essence, the detuning method retains the advantages of the two approaches discussed above without their limitations. The primary drawback of this method is that the detuning curves require long acquisition times due to the loss of intensity below threshold.

IV. RESULTS
A. Overview and assignment of phonon modes features are displayed as resolution-limited Gaussian lineshapes in Fig. 8(a) (shaded areas). They show distinct momentum dependences of their intensities. The intensity of the excitation at higher energy (≈ 70 meV, black dashed vertical bar) increases towards the X point of the BZ. This behavior is typical of the Cu-O bond stretching phonon branch, hereafter referred to simply as breathing branch. The intensity of the excitation at smaller energies (≈ 30 meV, red solid vertical bar) decreases towards X and is compatible with Cu-O bond-bending (buckling) modes. The difference in energy between the two phonon features is sufficient to resolve them, despite an overlap of the regions where their superposition is appreciable. The dispersion relations of the two phonon modes are plotted in Fig. 8(b), as obtained from the fit of the data. The weak dispersion extracted from the data supports our use of a local phonon model. To unambiguously identify the phonon branches seen in the spectra, one should consider that each RIXS feature is an average of the branches weighted over the response of the spectrometer and the probability of being seen with RIXS. As a consequence, only the phonon modes having a substantial interaction with the electrons contribute significantly to these averages. We consider the analogy with the isostructural YBCO, which is better known. This analogy is valid for our purposes as one expects the phonon dispersions of the two compounds to mainly differ at low energy, where the phonons involving the motion of Nd and Y are found [24,[67][68][69][70]. Figure 8(c) reports for convenience the theoretical results on the phonon dispersion relation of YBCO-AF from Ref. 66 (solid lines), limited to the branches whose EPC is particularly strong. The size of the vertical bars is proportional to the momentum-dependent EPC calculated for YBCO-OP [36], assumed to be qualitatively similar in the AF case. The major contributions arise from two branches above 50 meV (black lines) and two branches below 50 meV (red lines). In the upper region this confirms our initial conjecture on the relevance of the breathing branch, with a smaller contribution from the apical oxygen branch. Among the buckling phonons, there are branches involving the in-phase (A 1 ) and outof-phase (B 1 ) oscillations of the planar oxygen atoms. In our model, the out-of-phase one is silent due to symmetry selection rules, while the in-phase one is detected (see the calculations in Appendix A). In the model by Devereaux et al. [43], both branches are visible but the B 1 is weaker by a factor ∼ 20 with respect to A 1 , in agreement with the present result. As can be seen from Fig. 8(c), the EPC strength of the buckling mode is stronger at small momentum, while at large momentum the breathing mode prevails, in qualitative agreement with the phonon behavior observed in Fig. 8(a).

B. Determining the EPC in NBCO-AF
We start with an analysis of the breathing branch, which is clearly resolved in the spectra [ Fig. 9(a)]. From energy-detuned RIXS spectra, one can extract the detuning dependence of the phonon intensity. Figure 9(b) reports a set of theoretical detuning curves for the appropriate value of Γ/ω 0 = 4 corresponding to the Cu L 3 resonance. The comparison of the detuning curves with the phonon intensity (light blue squares) shows that we can achieve a good agreement with the data with a value of g ≈ 4, which corresponds to M = ω 0 √ g ≈ 0.13 eV. We note that a previous estimate of the EPC in NBCO-AF by means of the detuning approach yielded M ≈ 0.18 eV at h = −0.4 r.l.u. [65]. Both samples used in the measurements belong to the AF region of the phase diagram, but differ in aging and oxygen content. Extracting the EPC strength for the buckling branch is more delicate and gives the opportunity to employ the other methods presented in Sec. III. Figure 10(a) represents the benchmark RIXS spectrum of NBCO-AF measured at h = −0.4 r.l.u. (black circles), which is decomposed into an elastic line (red solid line) and an inelastic signal (blue squares). The latter is further decomposed in Fig. 10(b) into the buckling (red solid line) and breathing (black dashed line) contributions each with their first overtone (dotted lines). The decomposition is ambiguous if a naive fitting is applied; however, it can be made more robust if we constrain the ratio I 2 /I 1 for the breathing mode to the value of g found using the detuning method. Once this is done, the ratio I 2 /I 1 of the buckling mode is well determined. Moreover, knowledge of I 2 /I 1 at h = −0.4 r.l.u. allows us to calibrate the intensity of the buckling mode excitations I 1 such that its momentum dependence can be measured directly with I 1 on resonance. We use this approach since we could not measure the detuning curves at all momenta because of the limited availability of beamtime. At low momentum transfer the intensity I 1 tends to saturate because of the saturation of the universal curve, while the results of detuning are free of this problem at the Cu L 3 edge. Thus we used the detuning result at low Q . This is a typical case in which the direct connection between in- tensity and EPC does not work or works only in a subset of the parameter space. The final results on the EPC (in units of Γ) are summarized in Fig. 11. While the EPC of the buckling phonon branch (red circles) decreases going from the BZ center towards X, the EPC of the breathing branch (black squares) shows the opposite trend.

C. Doping effects
In the doped samples one expects that the photoexcited electron is screened by the free carriers. Thus, the phonon intensity is expected to monotonically decrease, as qualitatively shown in the cartoon of Fig. 12. However, a different scenario arises from the comparison between the measured RIXS spectra for the antiferromagnetic NBCO-AF (black circles), the underdoped NBCO-UD (red squares) and the (nearly) optimally doped NBCO-OP (blue diamonds) at h = −0.4 r.l.u. (Fig. 12). Note that this momentum transfer is well above the CDW wave vector (≈ 0.31 r.l.u. for NBCO [50]). The comparison between the intensities (normalized to the photon flux) among the three samples is reliable as can be seen in the spectral region between 0.6 -0.8 eV, where only the electron-hole continuum contributes to the spectra and scales with the doping. The behavior of the phonon intensity mimics the doping dependence of the CDW signal, which is stronger in the UD sample and weaker in the OP sample [54,71], as sketched in Fig. 12. The trend suggests that the CDW effect dominates over the screening even at Q considerably higher than the CDW wave vector.
As a matter of fact, there is an interplay between CDWs and phonons that strongly depends on the momentum transfer. This is evident from Fig. 13 Figure 12. The upper panel qualitatively shows two situations that modify the RIXS phonon signal: in the absence of charge order a smooth decrease of the EPC with doping is expected due to the increased screening from the free carriers. In the presence of CDWs, instead, CDW-induced effects may prevail so that a non-monotonic trend is found as a function of doping. Indeed, the signal from the CDWs is absent in the AF sample, is the strongest in the UD sample, and is weaker in the OP sample. The non-monotonic trend is not a necessary condition but it is sufficient to demonstrate the prevalence of the CDWs. This is indeed the case, as shown by the experimental results of the bottom panel, where spectra of NBCO-AF (black circles), NBCO-UD (red squares) and NBCO-OP (blue diamonds) are compared. The spectra were collected at h = −0.4 r.l.u. and are measured with an energy resolution of 60 meV.
NBCO-OP (red line). After subtraction of the elastic line and the continuum, we find that the phonon intensity is stronger in the doped system at large momentum transfer [ Fig. 13(c)]. The difference spectrum (blue line) does not vanish in the phonon region. In contrast, at small momentum transfer, the two dopings are equivalent within our sensitivity [ Fig. 13(d)]. Moreover, the response to detuning is very different between NBCO-AF and NBCO-OP, as shown in Fig. 14. At h = −0.4 r.l.u., the phonon intensity is still sizable in the OP compound upon detuning by Ω = −1 eV, whereas it has almost disappeared in the AF sample (black squares). This is evident from Fig. 14, in which the spectra have been normalized to the photon flux. The difference spectrum (blue diamonds) shows that both phonon branches (breathing and buckling) are more robust upon detuning in the doped system.

V. DISCUSSION
Here, we focus on the experimental results for NBCO since the methods used to recover the EPC from RIXS spectra have been already discussed. Nevertheless, these two aspects are connected so that the discussion of the results on NBCO will also clarify the limits of the methods.
In the parent compound the momentum dependence of the EPC plotted in Fig. 11 is fitted with simple trigonometric functions. The coupling to the buckling branch scales as cos(πh) (red solid line) and the coupling to the breathing branch as sin(πh) (black dashed line) (see Appendix A). Once again, we stress that these data are averages of unresolved branches. It is remarkable that such simple rules fit the data well, especially for the breathing mode whose error bars are smaller. This is important  Figure 14. The different behavior of the AF and the OP samples upon detuning. The doped sample (red circles) is definitely more robust than the undoped (black squares). The difference spectrum (blue diamonds) shows that both phonon branches detected by RIXS are more robust in the OP sample.
because the simple dependence of M/Γ on the momentum is theoretically obtained from models that consider only nearest neighbor interactions. The amplitude of the trigonometric functions has been adjusted by tuning the multiplicative factor so that: M Γ buckle = 0.98 cos(πh).
These results demonstrate that RIXS provides easy and direct access to the momentum-dependent EPC. Note that our results are in qualitative agreement with the expectations based on the analogy with YBCO [see Fig. 8(c)]. Moreover, the analysis of the RIXS data of NBCO-AF and the internal consistency of the results obtained with different methods demonstrate that our model involving the starting Hamiltonian (Eq. 1) and Einstein phonons is adequate to study high energy phonons in undoped cuprates. It is remarkable that the intensity of the in-phase buckling branch is high and has a strong EPC; this is consistent with the large static buckling in the system, which is known to enhance the coupling to the bond-buckling phonons [3,72].
In the doped cuprates it is known that the CDWs are ubiquitous, at least in the UD samples and with lower intensity in the OP compounds [53,54,71]. We have already shown that in NBCO-OP the role of the CDWs is crucial. More specifically, we stress the following points: • The modification of the spectral weight in the phonon region of NBCO follows the general pattern of Bi 2 Sr 2 CaCu 2 O 8+δ reported in Ref. 19 in that the increase of the intensity at h = −0.4 r.l.u. is analogous to the so-called phonon anomaly [12,19].
• At h = −0.1 r.l.u., after background subtraction, the phonon intensity is basically the same in the AF and OP samples. This means that our model (especially in the detuning formulation) can be used also in OP systems, provided the momentum transfer is sufficiently small. A more precise assessment is left for future investigations, but it is already clear that the model is meaningful around h = −0.1 r.l.u.. This region of the parameter space is of great interest since the coupling to the buckling modes are close to the maximum and the coupling of the breathing mode is close to zero. This situation as a whole is favorable to pairing and makes buckling phonons good candidates for a synergistic action together with magnetic fluctuations. Basically, our approach works well in the region of the parameter space that is most relevant to superconductivity in the cuprates. Note that according to Ref. 7 the buckling modes are beneficial to pairing, in particular the out-of-phase B 1 mode, whose RIXS intensity vanishes for symmetry reasons (see Appendix A).
Since the screening is comparable for the A 1 and B 1 modes, we can argue that the out-of-phase mode has comparable coupling to the electron than the in-phase mode. The effect of the breathing modes, instead, is now considered detrimental, although this has been controversial [6,7,[33][34][35].
• At large momentum transfer, the situation is completely different: the effect of the CDWs cannot be neglected in the OP sample and, a fortiori, in UD compounds. This is remarkable because the CDW signal in OP cuprates is very small and sometimes considered to be zero [50]. Thus the sensitivity of the EPC to the CDWs is extremely high. This prevents the use of our model and calls for a more sophisticated treatment.
The above arguments clarify the limits of validity of the model itself, i.e. far from CDW-induced effects. It is interesting, however, to try to use the model also to evaluate an effective EPC, M eff , in the region where the CDWs are important. This very crude approach can be applied to the energy-detuned RIXS spectra plotted in Fig. 14 and measured at h = −0.4 r.l.u.. For the breathing mode we obtain M eff ∼ 3 -4 times the value in the parent compound, which is a coupling strength more typical of polaronic systems. Note also that, on the basis of the ratio M eff /M , it would be possible to introduce a scale to characterize how far we are from the situation free of CDWs.

VI. CONCLUSIONS
We have presented a way to obtain quantitative information on the electron-phonon coupling in undoped and superconducting NBCO from high resolution Cu L 3 RIXS spectra. Importantly, this approach can be applied without numerical simulations of the spectra be-cause we have introduced simple procedures based on rescaled variables and on energy detuning. These methods, together with the further improvements of the experimental resolution in the coming years, can be extremely relevant in the study of electron-lattice interactions in strongly correlated materials. A support of this view is given by the original results on the EPC in undoped and doped NBCO. These results are likely to become a paradigm for the "123" family of cuprates as a whole. intensity at the Cu L 3 edge is then given by: Here, the scattering cross section is defined as F f g = f |D out H e−ph D in |g , where |g and |f are the initial and final states, respectively, with energies E g and E f , ω is the energy loss, and D in and D out are the dipole transition operators. The latter are defined as where d † i,σ creates an electron in the Cu 3d x 2 −y 2 orbital located at position R i and p † i,σ creates an electron in the relevant Cu 2p core level, and q in and q out are the incident and scattered photon wave vectors, respectively. Note that we have neglected the polarization-dependent prefactors in the dipole matrix elements for simplicity.
We assume that the ground state can be written in the form |g = |ψ el , n q,ν = 0 , which describes the electronic subsystem with no phonon quanta excited. Here, we are interested in quasi-elastic processes, where the energy transferred into the system excites a phonon. Therefore, we can restrict the final states to only those where one phonon has been excited, i.e. |f = |ψ el , n q,ν = 1 . This assumption is equivalent to the view that the phonon and electron subsystems are not deeply entangled such that there are zero phonons present in the ground state. Under these simplifying assumptions we have: In the last step we have introduced the momentum transfer Q = q in − q out . The Cu orbital operator d † k,σ is related to the c † k,σ band operator by c † k,σ = φ Cu (k)d † k,σ , where φ Cu (k) measures the Cu character of the antibonding band. Therefore, the scattering matrix element written in band space is: M ν (k, q)φ Cu (p + Q)φ Cu (p) ψ el , n q,ν = 1 c p,σ c † k−q,σ c k,σ c † p+Q,σ ψ el , n q,ν = 1 .
The expectation value appearing in Eq. (A4) must be evaluated using the correlated many-body wave function. We can, however, simplify the problem by considering either the fully localized or fully itinerant cases. First, we consider the itinerant case and approximate the many-body wave function as the non-interacting Fermi sea. Then, for non-zero values of Q, we require p = k − q and Q = q, and the scattering amplitude simplifies to: Here, n F (x) is the Fermi factor and (k) is the band dispersion. Similarly, we obtain the localized limit by inverse Fourier transforming Eq. (A4) and then retaining only the local operators inside the expectation value: p,k,q σ,ν,l M ν (k, q)φ Cu (p + Q)φ Cu (p)e −i(Q−q)·R l ψ el , n q,ν = 1 |(1 − n l )| ψ el , n q,ν = 1 = 1 N 3 2 p,k σ,ν M ν (k, Q)φ Cu (p + Q)φ Cu (p).
Thus, in both limits, RIXS at the Cu L 3 edge measures a k-integrated coupling constant weighted by the Cu orbital character of the band and additional phase space factors. The involvement of the Cu orbital character is in agreement with conclusions drawn in a previous cluster ED study of a quasi-1D cuprate [45]. In the main text we have adopted the localized limit, since we expect strong electron correlations and the sizable core hole potential to localize the excited core electron to the Cu site where it is created. This approach also has the advantage that simple analytic expressions for the intensity can be obtained.
Reference 7 argued that φ Cu (k) can be approximated by a constant φ Cu (k) ∼ A Cu , which leads to simplified momentum dependencies for the various coupling constants. In this case, A Cu can be absorbed into the prefactor of the coupling constant and the scattering amplitude for the localized limit reduces to F f q = 1 √ N k,σ,ν M ν (k, Q). For the in-phase (+) and out-of-phase (−) Cu-O bond-buckling branches one has (p = k − Q): M ± (k, Q) ∝ sin k x a 2 sin p x a 2 cos q y a 2 ± sin k y a 2 sin p y a 2 cos q x a 2 . (A7) Similarly, for the breathing mode: M breath (Q) ∝ sin 2 q x a 2 + sin 2 q y a 2 1/2 .
Because the coupling constant for the breathing mode does not depend on the Fermion momentum k, the k-averaging is trivial and I breath (Q) ∝ |M breath (Q)| 2 . The situation is different for the buckling modes. Due to the sign change in the k-dependence of M − (k, Q), the intensity for the out-of-phase buckling mode I − (Q) vanishes. For the in-phase buckling phonon branch one arrives at: I + (Q) ∝ cos 2 q x a 2 cos 2 q y a 2 .
(A9) Figure 15 displays the EPC strength M (Q) of the in-phase buckling (red circles) and breathing (black squares) modes (see Fig. 11 and related discussion), normalized to the natural width Γ of the Cu L 3 resonance. The EPC strength is plotted as a function of momentum transfer along the (1, 0) direction of the reciprocal space (h = q x a/2π). Note that the EPC strength is proportional to the square root of the RIXS intensity. The best fits to the data obtained in the two limiting cases of local (thick lines) and itinerant (thin lines) electron models are shown for comparison. In the latter case, the band dispersion is taken from Ref. 73. The two approximations yield similar results and describe the the data fairly well, with a small advantage to the local electron picture close to the X point of the BZ. Liechtenstein, "Out-of-plane instability and electronphonon contribution to s-and d-wave pairing in hightemperature superconductors; LDA linear-response cal-