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}$

TL;DR

This paper develops a method to quantify electron-phonon coupling in high-temperature superconductors using resonant inelastic X-ray scattering, and applies it to Nd-based cuprate films to analyze phonon interactions.

Contribution

It introduces a theoretical framework to extract EPC strength directly from RIXS data and demonstrates its application to various cuprate samples, enabling detailed EPC analysis.

Findings

EPC strength varies with momentum transfer in cuprates.

Coulomb screening affects EPC measurements in doped samples.

Framework sets stage for future high-resolution RIXS studies.

Abstract

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.

Tables

Table 1. Table 1: The theoretical HWHM Γ Γ \Gamma for selected elements and absorption edges in units of eV. Adapted from Refs.  Keski-Rahkonen and Krause, 1974 ; Krause and Oliver, 1979 .

$Z$	Elements	$K$ edge	$L_3$ edge	$M_5$ edge
6 - 9	C - F	0.05 - 0.10
11 - 17	Na - Cl	0.15 - 0.34
21 - 30	Sc - Zn	0.43 - 0.85	0.09 - 0.33
39 - 48	Y - Cd	1.7 - 3.6	0.07 - 1.2
57 - 70	La - Yb	7.0 - 16	1.7 - 2.3	0.15 - 0.30
72 - 80	Hf - Hg	18 - 27	2.4 - 2.7	1.0 - 1.5

Equations

$$\mathcal{H}=\sum _i{\omega }_0^{}{b}_i^{†}{b}_i^{}+M\sum _i{d}_i^{†}{d}_i^{}(b_i^{†}+b_i^{}),$$

$$\mathcal{F}=T_{\mathrm{el}}({\mathbf{ϵ}}^{\prime },\mathbf{ϵ})\sum _i{e}^{i\mathbf{Q}\cdot {\mathbf{R}}_i}\sum _{n_i=0}^{\infty }\frac{⟨n_i^{\prime }|e^{-S_i}|n_i⟩⟨n_i|e^{S_i}|n_i^0⟩}{z+M^2/{\omega }_0-n_i{\omega }_0}.$$

$$B_{ab}=\sqrt{e^{-g}a!b!}\sum _{l=0}^b\frac{(-1{)}^a(-g{)}^l{g}^{(a-b)/2}}{(b-l)!l!(a-b+l)!},$$

$$\frac{{\mathrm{d}}^2\sigma }{\mathrm{d}\Omega \mathrm{d}\omega }\propto \sum _f|\mathcal{F}{|}^2\delta (\omega -n^{\prime }{\omega }_0)=N|T_{\mathrm{el}}({\mathbf{ϵ}}^{\prime },\mathbf{ϵ}){|}^2\sum _{n^{\prime }=0}^{\infty }{\left|\sum _{n=0}^{n^{\prime }}\frac{B_{n^{\prime }n}(g)B_{n0}(g)}{z+(g-n){\omega }_0}+\sum _{n=n^{\prime }+1}^{\infty }\frac{B_{n{n}^{\prime }}(g)B_{n0}(g)}{z+(g-n){\omega }_0}\right|}^2\delta (\omega -n^{\prime }{\omega }_0),$$

$$I_1$$

$$I_2$$

$$\frac{W}{{\omega }_0}=1.8+1.3\text{leavevmode}{g}^{\frac{1}{2}}=1.8+1.3\text{leavevmode}\frac{M}{{\omega }_0}.$$

$${\left(\frac{M}{\Gamma }\right)}_{\mathrm{breath}}=0.45\sin (\pi h),$$

$${\left(\frac{M}{\Gamma }\right)}_{\mathrm{buckle}}=0.98\cos (\pi h).$$

$${\mathcal{H}}_{\mathrm{e}-\mathrm{ph}}=\frac{1}{\sqrt{N}}\sum _{\mathbf{k},\mathbf{q},\sigma ,\nu }{M}_{\nu }^{}(\mathbf{k},\mathbf{q})c_{\mathbf{k}-\mathbf{q},\sigma }^{†}{c}_{\mathbf{k},\sigma }^{}(b_{\mathbf{q},\nu }^{†}+b_{-\mathbf{q},\nu }^{}).$$

$$I(\mathbf{Q},\omega )\propto |{\mathcal{F}}_{fg}{|}^2\delta \left(E_g-E_f+\omega \right).$$

$$D_{\mathrm{in}}$$

$$D_{\mathrm{out}}$$

$${\mathcal{F}}_{fg}$$

$${\mathcal{F}}_{fg}=\frac{1}{\sqrt{N}}\sum _{\begin{array}{c}\mathbf{p},\mathbf{k},\mathbf{q}\\ \sigma ,\nu \end{array}}{M}_{\nu }^{}(\mathbf{k},\mathbf{q}){\varphi }_{\mathrm{Cu}}^{\star }(\mathbf{p}+\mathbf{Q}){\varphi }_{\mathrm{Cu}}^{}(\mathbf{p})⟨{\psi }_{\mathrm{el}},n_{\mathbf{q},\nu }=1\left|c_{\mathbf{p},\sigma }^{}{c}_{\mathbf{k}-\mathbf{q},\sigma }^{†}{c}_{\mathbf{k},\sigma }^{}{c}_{\mathbf{p}+\mathbf{Q},\sigma }^{†}\right|{\psi }_{\mathrm{el}},n_{\mathbf{q},\nu }=1⟩.$$

$${\mathcal{F}}_{fg}=\frac{1}{\sqrt{N}}\sum _{\mathbf{k},\sigma ,\nu }{M}_{\nu }^{}(\mathbf{k},\mathbf{Q}){\varphi }_{\mathrm{Cu}}^{\star }(\mathbf{k}){\varphi }_{\mathrm{Cu}}^{}(\mathbf{k}-\mathbf{Q})[1-n_{\mathrm{F}}({ϵ}_{\mathbf{k}-\mathbf{Q}})][1-n_{\mathrm{F}}({ϵ}_{\mathbf{k}})].$$

$${\mathcal{F}}_{fg}$$

$$M_{\pm }(\mathbf{k},\mathbf{Q})\propto \left[\sin \left(\frac{k_{x}a}{2}\right)\sin \left(\frac{p_{x}a}{2}\right)\cos \left(\frac{q_{y}a}{2}\right)\pm \sin \left(\frac{k_{y}a}{2}\right)\sin \left(\frac{p_{y}a}{2}\right)\cos \left(\frac{q_{x}a}{2}\right)\right].$$

$$M_{\mathrm{breath}}(\mathbf{Q})\propto {\left[{\sin }^2\left(\frac{q_{x}a}{2}\right)+{\sin }^2\left(\frac{q_{y}a}{2}\right)\right]}^{1/2}.$$

$$I_{+}(\mathbf{Q})\propto {\cos }^2\left(\frac{q_{x}a}{2}\right){\cos }^2\left(\frac{q_{y}a}{2}\right).$$

Full text

I Introduction

The role of the electron-phonon coupling (EPC) in the high critical temperature (high-$T_{\mathrm{c}}$) superconducting cuprates is still an open problem deserving further research. Indeed, even if pairing is not of the phonon-mediated BCS type Scalapino (2012), the question of a potential role for the EPC remains extremely interesting. For instance, it has been suggested theoretically Savrasov and Andersen (1996); Andersen et al. (1996); Sakai et al. (1997); Jepsen et al. (1998); Ishihara and Nagaosa (2004); Johnston et al. (2010) that a synergy between a suitable phonon and other pair-driving excitations can greatly enhance the critical temperature $T_{\mathrm{c}}$. In particular, in the case of magnetic excitations, even a small amount of EPC (per se irrelevant) should be sufficient to considerably increase $T_{\mathrm{c}}$ Johnston et al. (2010). Another reason for studying the EPC is the recent observation of transient superconductivity induced in cuprates by strong illumination with mid-infrared pulses Fausti et al. (2011); Kaiser et al. (2014); Hu et al. (2014); Liu et al. (2019), i.e., more precisely, the increase of $T_{\mathrm{c}}$ under optical pumping of apical oxygen phonons in the non-linear regime. Moreover, the interplay of phonons with the electronic states is very important in the presence of charge-density-wave (CDW) correlations. Indeed, the low-energy acoustic Le Tacon et al. (2013); Miao et al. (2018) and high-energy Cu-O bond-stretching McQueeney et al. (1999); Pintschovius et al. (2004); Uchiyama et al. (2004); Reznik et al. (2006); Graf et al. (2008); Chaix et al. (2017) phonon branches are strongly modulated in intensity and energy 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 Kastner et al. (1998); Alexandrov (2000); Rösch et al. (2005); Shen et al. (2007).

Clarifying the role of the EPC in the cuprates and 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 Bohnen et al. (2003); Giustino et al. (2008); Reznik et al. (2008); Heid et al. (2008); Veenstra et al. (2010) and experiment Lanzara et al. (2001); Pintschovius et al. (2004); Zhou et al. (2005); Cuk et al. (2005); Reznik et al. (2006); Johnston et al. (2012) but also among different theoretical approaches Savrasov and Andersen (1996); Bohnen et al. (2003); Giustino et al. (2008); Reznik et al. (2008). For example, the Cu-O bond-stretching modes have been proposed to mediate both an effective attractive Shen et al. (2002); Ishihara and Nagaosa (2004) and, more recently, repulsive Bulut and Scalapino (1996); Sandvik et al. (2004); Johnston et al. (2010) interaction.

The experimental verification of the phenomenology of the EPC, possibly with momentum resolution, is crucial in this context. This goal calls for the improvement of the existing techniques and for the introduction of new methods. Among traditional measurements, bulk-sensitive neutron scattering Pintschovius (2005); Reznik (2010) suffers from sensitivity limitations as it requires massive homogeneous samples, whereas surface-sensitive angle-resolved photoemission spectroscopy Damascelli et al. (2003); Cuk et al. (2005) and scanning tunneling microscopy Lee et al. (2006) need high quality surfaces that are not always available. Moreover, the comparison with bulk sensitive results is a nontrivial task. In this field, resonant inelastic soft x-ray scattering (RIXS) is likely to offer an important contribution, made possible by the recent progress in instrumentation Brookes et al. (2018). For instance, state-of-the-art RIXS can be performed with a combined energy resolution of $\sim 30$ – 40 meV at the Cu $L_{3}$ edge, which facilitates studies of high-energy phonons ($\gtrsim 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.

Determining the EPC from RIXS data was first proposed by Ament et al. Ament et al. (2011) (for a detailed account see Ref. Ament, 2010). This idea has also been reinforced by the theoretical treatment by Devereaux et al. Devereaux et al. (2016), based on the lowest order Feynman diagrams. Devereaux et al.’s 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 Geondzhian and Gilmore (2018) with a model conceptually similar to Ref. Ament et al., 2011.

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 Lee et al. (2013); Johnston et al. (2016), having in general more pronounced phonon peaks and van Hove singularities. A recent remarkable work has also reported the first application to $n$SrIrO3/$m$SrTiO3 multilayers Meyers et al. (2018). For the low-Z $3d$ oxides, data on Ti $L$ edge are also available Moser et al. (2015); Fatale et al. (2016).

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 one. Specifically, we introduce innovative ways of using the theory by Ament et al. Ament et al. (2011) and identify some scaling laws with a wide range of application. The proposed methods are then employed on the cuprate Nd1+xBa2-xCu3O7-δ (NBCO). We examine the EPC in the antiferromagnetic (AF) compound and the effect of charge order Ghiringhelli et al. (2012); da Silva Neto et al. (2014, 2015); Comin and Damascelli (2016); Arpaia et al. (2019) on the EPC at different doping levels. From the RIXS data, we quantitatively determine the EPC strength in NBCO-AF and qualitatively 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 perhaps the strongest in this case Comin and Damascelli (2016). 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 an in-depth discussion of our theoretical framework. We first review the main results of the theory as developed previously. We then thoroughly explore the origin of lattice excitations during the RIXS process and discuss the role of the various parameters that govern the phonon intensity. Importantly, we unveil universal scaling laws that involve the main parameters of the theory. We also examine three methods for determining the EPC from RIXS data. Although some of them have been already introduced elsewhere Ament et al. (2011); Rossi et al. (2019), we provide novel and efficient ways to exploit them based on some general scaling laws. We also compare the methods and point out their advantages and limitations. Moreover, we show how to extract the momentum dependence of the EPC by combining experimental and theoretical information. Finally, we conclude Section III by briefly addressing the case of absorption edges other than the Cu $L_{3}$ edge. Section IV presents the experimental RIXS results on the bilayer cuprate NBCO, with emphasis on the identification of the phonon modes, determination of the EPC on a quantitative basis, and an inspection of doping effects. Finally, Section V provides some additional discussion of our results.

II Experimental details

II.1 The NBCO samples

The crystal structure of Nd1+xBa2-xCu3O7-δ (NBCO) is displayed in Fig. 1(a) for the undoped antiferromagnetic parent compound ($x=0$ and $\delta=1$). NBCO belongs to the so-called “123” family of high $T_{\mathrm{c}}$ superconductors and is isostructural to YBa2Cu3O7-δ (YBCO). Of primary importance in the centrosymmetric unit cell is the presence of the CuO2 bilayer where superconductivity occurs. The hole density per CuO2 plane can be changed by modifying the oxygen content $\delta$ (similarly to YBCO Arpaia et al. (2018)) or by introducing excess Nd $x$ at the Ba sites. In the latter case, 0.5 holes are removed from the CuO2 planes for each excess Nd Salluzzo et al. (2002, 2005). At optimal doping, $T_{\mathrm{c}}$ reaches 93 K. NdBa2Cu3O6 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 CuO2 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 Chmaissem et al. (1999).

Epitaxial NBCO films were deposited on a $(001)$-oriented SrTiO3 substrate with an almost perfect in-plane matching of the $a$ and $b$ lattice parameters and the buckling practically unchanged with respect to the corresponding bulk structure Salluzzo et al. (2002, 2005). The thickness of the NBCO films is around 150 nm, hence they can be considered infinitely thick for the scattering angles used in our experiment (where the incidence angle with respect to the surface ranged from $25\text{\,}\mathrm{\SIUnitSymbolDegree}$ to $65\text{\,}\mathrm{\SIUnitSymbolDegree}$).

We investigated three NBCO samples:

• •

NBCO-AF: the sample is antiferromagnetic with $x\approx 0$, $\delta\approx 0.9$, and with a Cu-O-Cu buckling angle of $6.35\text{\,}\mathrm{\SIUnitSymbolDegree}$ Takita et al. (1988); Kramer et al. (1994);

• •

NBCO-UD: the sample is underdoped with $T_{\mathrm{c}}=63$ K, hole doping of 0.11, $x\approx 0.2$, $\delta\approx 0$, and a buckling angle of $6.74\text{\,}\mathrm{\SIUnitSymbolDegree}$ Izumi et al. (1987);

• •

NBCO-OP: the sample is very close to optimal doping with $T_{\mathrm{c}}=90$ K, $x\approx 0$, $\delta\approx 0$, and a buckling angle of $7.75\text{\,}\mathrm{\SIUnitSymbolDegree}$ Takita et al. (1988). The slightly lower $T_{\mathrm{c}}$ compared to optimal doping is due to a slight overdoping, which has been identified by measuring the $c$ lattice parameter.

While the superconducting NBCO films are stable in time, the AF samples are more likely to grow in doping by capturing oxygen from the atmosphere when stored for weeks or months at room temperature. As a result, the actual doping of the sample here named AF was not perfectly zero, and the upper limit for its hole doping level is $p\simeq 0.04$. In RIXS spectra, a fingerprint of the doping level appears in the formation of a continuum in the otherwise empty region around 0.8 eV energy loss. We were thus able to monitor the sample during the experiment to ensure that it stayed in the AF insulating region of the phase diagram. In any case, the exact doping level of this sample is not relevant for the present work. The UD and OP samples are at risk of loosing oxygen if kept in vacuum at room temperature. This risk is why we have performed all the measurements at 20 K, where radiation induced and spontaneous loss of oxygen is minimized.

II.2 Experimental layout and energy resolution

The measurements were performed at beam line ID32 of the European Synchrotron (ESRF, France), which is equipped with the new spectrometer ERIXS (European RIXS) Brookes et al. (2018). 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 ($\mathbf{q}_{\mathrm{in}}$ and $\mathbf{q}_{\mathrm{out}}$, respectively). The scattering angle $2\theta$ was fixed to $149.5\text{\,}\mathrm{\SIUnitSymbolDegree}$. Most of the measurements were performed with a combined (beam line $+$ spectrometer) energy resolution of $\sim 40$ meV at the Cu $L_{3}$ edge ($\sim 930$ eV). The contribution to the resolution from the incident beam is $\sim 25$ meV. One spectrum was collected with high statistics and high energy resolution (32 meV) to be used as benchmark. The sample was kept at a temperature of 20K during the measurements.

The momentum transferred to the sample $\mathbf{Q}=\mathbf{q}_{\mathrm{in}}-\mathbf{q}_{\mathrm{out}}$ and projected on the basal $ab$ plane is the relevant observable in dispersion studies of two-dimensional cuprates. The projected in-plane momentum $\mathbf{Q}_{\parallel}$ is scanned in the $\Gamma$-$X$ direction of the Brillouin zone (BZ) by rotating the sample around the $b$ axis, i.e., $\mathbf{Q}_{\parallel}=\left(\frac{2\pi}{a}h,0\right)$. 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 [$\sigma$ incident polarization, red arrow in Fig. 1(b)]. Thus, the electric field vector lay in the basal plane, simplifying the problem (see Sec. III.1). We did not investigate small momentum transfers $|\mathbf{Q}_{\parallel}|<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. In principle the lattice and magnetic excitations can be separated by measuring the scattered beam polarization. For example, such an analysis has been performed in Ref. Hepting et al., 2018, where the use of the polarimeter allowed the authors to isolate the plasmon signal from the magnetic one. The present case is technically more complex, however, and so we have postponed this part of the work until after the commissioning of the new storage ring at the ESRF is complete. Note that the analysis of the scattered beam polarization in the soft x-ray range is possible at present only at beamline ID32 of the ESRF.

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.

III Tools to determine the electron-phonon coupling from RIXS spectra

III.1 Theoretical background

The present work is based on the theory by Ament et al. Ament (2010); Ament et al. (2011), which was laid down for localized electrons coupled to dispersive phonons via a momentum-dependent EPC $M({\bf k},{\bf q})$ and then simplified to the case of coupling 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 Pintschovius (2005); Reznik (2010). Moreover, the assumption of a localized electrons can be justified by appealing to the strong core-hole potential that will have a tendency to localize the excited core electron. We will come back to this issue in section V. Here, we do not introduce any changes in the theory but, instead, 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. Ament, 2010.

In the case of dispersionless phonons – which is the situation we are interested in – the starting Hamiltonian is

[TABLE]

where $\omega_{0}$ is the phonon energy, $M$ is the EPC matrix element, and $b_{i}^{\dagger}$ ($b^{\phantom{\dagger}}_{i}$) and $d_{i}^{\dagger}$ ($d^{\phantom{\dagger}}_{i}$) are the creation (annihilation) operators for phonons and electrons, respectively, at site $i$. We have neglected the spin index for brevity. Eq. (1) can be diagonalized exactly using a canonical Lang-Firsov transformation $\mathcal{H^{\prime}}=e^{S}\mathcal{H}e^{-S}$, where $S=\sum_{i}S_{i}$ and $S_{i}=\frac{M}{\omega_{0}}d^{\dagger}_{i}d^{\phantom{\dagger}}_{i}\left(b^{\dagger}_{i}-b^{\phantom{\dagger}}_{i}\right)$. This same transformation can be used to derive an exact expression for the scattering amplitude:

[TABLE]

Here, $T_{\mathrm{el}}$ is the polarization-dependent atomic elastic scattering factor; $z=\Omega+i\Gamma$ is a complex number whose real part corresponds to the detuning energy $\Omega$ (i.e., the difference between the incoming photon energy and the resonance energy) and whose imaginary part $\Gamma$ is related to inverse lifetime of the core hole. Note that $\Gamma$ is the Half Width at Half Maximum (HWHM) of the core line. For simplicity we will refer to $1/\Gamma$ as the inverse lifetime. $|n_{i}^{0}\rangle$ denotes the initial ground state of the Hamiltonian, which corresponds to the ground state of a shifted Harmonic oscillator.

By introducing the Franck-Condon (FC) coefficients

[TABLE]

the cross section can be written as:

[TABLE]

where $g=(M/\omega_{0})^{2}$ is a dimensionless coupling constant and the index $n^{\prime}$ identifies the number of phonons in each final state $f$ indexing the first summation. We note that a sum with index $n$, the number of phonons excited in the intermediate state, has to be extended to a potentially very large number (in principle to infinity).

For convenience, we report the exact expressions for the one-phonon intensity $I_{1}$ and two-phonon intensity $I_{2}$, corresponding to the $n^{\prime}=1,2$ final states, respectively:

[TABLE]

Note that a common prefactor has been dropped from the two equations.

III.2 The genesis of the phonon excitations during the RIXS process

The excitation of a phonon during the Cu $L$-edge RIXS process is sketched in Fig. 2. Upon photon absorption at the Cu $L_{3}$ resonance, an excited state is created having a $2p_{3/2}$ core hole and an extra electron in a $3d$ orbital. In layered cuprates, the ground state has a $3d_{x^{2}-y^{2}}$ symmetry of a single local hole, corresponding to an empty antibonding molecular orbital, which is temporarily filled in the intermediate state. The intermediate state is, therefore, not an eigenstate of the electron-lattice coupled system, and the lattice will deform towards the new equilibrium structure where the oxygen atoms are at a different distance from the central Cu ion. This intermediate state can decay back to the original $3d^{9}$ configuration but with one or several excited phonons. The detailed nature of the lattice excitations is specific to the intermediate state (and thus depends on the core and valence states involved in the RIXS process) but the general scheme is always the one we have described. Note that the intermediate state has excitonic character, i.e., it is charge neutral. In a sense, we can consider 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_{\mathrm{t}}$ involved in the transport measurements, as stressed in Ref. Geondzhian and Gilmore, 2018, because of the core hole effect and the symmetry of the intermediate state. However, the strong screening of the core hole by the excited electron reduces the difference between $g$ and $g_{\mathrm{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_{\mathrm{t}}$ close.

Let us first consider the effect of the core hole lifetime, which introduces a time scale into the problem. It is intuitive that the phonon signal increases with the EPC matrix element $M$ and that it decreases with the core hole lifetime ($1/\Gamma$), as shorter lived core hole states provide less time for the lattice to evolve towards the new equilibrium configuration. The logical consequence is that the phonon excitation efficiency by RIXS depends on the ratio $M/\Gamma$, irrespective of the value of $g$, over a broad range of $M/\Gamma$.

The universal dependence of the scattering cross section on $M/\Gamma$ is confirmed by numerical calculations made with Eqs. (5) and (6) and is summarized in Fig. 3. On the left axis, we plot the intensity of the one-phonon excitation $I_{1}$, normalized to $(\omega_{0}/\Gamma)^{2}$, for several values of the dimensionless coupling constant $g=(M/\omega_{0})^{2}$. It is evident that all curves corresponding to different $g$ values collapse on the same “universal” curve when $(M/\Gamma)^{2}\lesssim 2$. Above this value, the curves begin to diverge and the behavior is no longer universal. Similar considerations apply to the ratio of the intensities of the one-phonon and two-phonon excitations $I_{2}/I_{1}$ (right axis). We note that for an Einstein phonon, the two-phonon excitation has exactly 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/\Gamma)^{2}$ is independent of $g$ for $(M/\Gamma)^{2}\lesssim 1$. At the Cu $L_{3}$ resonance, $\Gamma\approx 250$ – $280$ meV Keski-Rahkonen and Krause (1974); Krause and Oliver (1979) so the universality range extends up to about $M\approx 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/\Gamma)^{2}$.

For weak interactions, the intensity $I_{2}$ is directly proportional to $(I_{1})^{2}$ by definition Ament et al. (2011), 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/\Gamma)^{2}\approx 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 to extract the EPC from experimental data. A critical comparison between the different methods is presently lacking in the literature and we summarize a number of key points in the following subsections.

III.3 The use of $I_{2}/I_{1}$

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/\Gamma)^{2}$. This approach has been used in early literature at the O $K$ edge Lee et al. (2013); Johnston et al. (2016); Meyers et al. (2018) and at the Ti $L_{3}$ edge Moser et al. (2015); Fatale et al. (2016) 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 $\omega_{0}$, and the core hole lifetime $1/\Gamma$. The value of $I_{2}/I_{1}$ corresponds to the point P1 on the universal curve (left panel) and thus to a specific value of $M/\Gamma$. The horizontal line through P1 intersects the value of $(\omega_{0}/\Gamma)^{2}$, which is the entrance value of the right panel. This defines the point P2. The straight line from the origin through P2 defines the angle $\alpha$ whose tangent is $(M/\Gamma)^{2}/(\omega_{0}/\Gamma)^{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.

The $I_{2}/I_{1}$ method may seem ideal, but this approach has 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. 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. This method does, however, work particularly well in the case of quasi-1D systems characterized by a single non-dispersing lattice mode, which can give rise to sharp and intense multi-phonon peaks Lee et al. (2013); Johnston et al. (2016); Vale et al. (2019).

III.4 The use of $I_{1}$ on resonance

The universal curve of Fig. 3 can also be used to recover $(M/\Gamma)^{2}$ from the intensity of the single phonon excitation $I_{1}$. The geometrical construction shown in Fig. 4 is formally the same, but the ratio $I_{2}/I_{1}$ is 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 relation between $I_{1}$ and the EPC is determined at least at one point by another method, and if a relative measurement of $I_{1}$ can be made as a function of the relevant parameter (momentum transfer, temperature, doping, etc.). Fig. 3 also highlights the non-proportionality and non-monotonicity of $I_{1}$ as a function of 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 electron-lattice interaction. The existence of a maximum for $I_{1}$ means that the plot equivalent to that of Fig. 4 based on $I_{1}$ would no longer be single valued, creating some ambiguities for strong EPC. One might be able to resolve this issue by considering the nature of the physical problem.

III.5 The use of $I_{1}$ upon detuning below threshold

This method has been recently introduced by some of us in a concise way Rossi et al. (2019). We add here some relevant information, including a systematic overview of the method, a straightforward way to retrieve the EPC from energy-detuned RIXS spectra, and the identification of an approximate scaling law. On this basis we also introduce a remarkably simple rule allowing, in many cases, an immediate approximate estimate of the EPC. The detuning approach 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 $\Omega=\omega-\omega_{\mathrm{res}}$, of $i\Gamma$ and of $M^{2}/\omega_{0}$. Thus, detuning has a larger effect if $M$ is small. Indeed, there is a kind of trade-off between detuning and $M^{2}/\omega_{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 to $-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, normalized 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 $\Gamma/\omega_{0}$ and contains the detuning curves for different values of $g$. Since $\Gamma/\omega_{0}$ is known for a given absorption edge and phonon mode, 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 $\Gamma/\omega_{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(a), which reports the width $W$ (in units of $\omega_{0}$) of the detuning curve, defined as the detuning value of the point at half height. The width is displayed as a function of $\sqrt{g}$, with $\Gamma/\omega_{0}$ as a parameter. The lines are remarkably linear in a wide range and at lower values of $\Gamma/\omega_{0}$ they are very close to each other, suggesting an approximate scaling law.

The procedure based on Fig. 7(a) is exact within the limits of Ament et al.’s model, since no further approximation has been introduced. However, if one accepts a modest inaccuracy, it is possible to further simplify the procedure in many cases. Here, we take advantage of the fact that the lines in Fig 7(a) become close to each other when $\Gamma/\omega_{0}$ decreases. In this parameter region, we approximate the state of the system with a single line. This is done in Fig. 7(b) by introducing the red straight line to represent the system in the interval $2\leq\Gamma/\omega_{0}\leq 8$. By using the red line the scaling becomes universal, i.e., independent of $\Gamma$. The red line has the equation:

[TABLE]

If one enters the diagram with a value of $W$, the working point is the intersection of the horizontal line representing $W$ with the red line in place of the exact point, as illustrated in Fig. 7(b) for the case $W=10$. The inaccuracy is of the order of 10%. Obviously the error is maximum at $\Gamma/\omega_{0}$ equal to 2 and 8. This astonishingly simple rule can be extremely useful during experiments since it offers a shortcut to estimate the EPC from the raw data. We will refer to this rule as the Simplified Detuning Rule (SDR).

The method based on detuning is advantageous because it uses the value of $I_{1}$ normalized to its value at resonance, which is easily measured. 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.

III.6 Access to momentum-dependent EPC

As mentioned above, the direct estimation of the EPC from our spectra is based on the important assumptions made by Ament et al.; in particular, the electronic state involved in the RIXS process is localized (i.e. an excitonic intermediate state) and the phonon is assumed to have a local character. We will address later the effect of the localized vs itinerant nature of the electron states in connection with EPC. Here, we limit ourselves to note that the two assumptions are less stringent than it might seem. Very often the RIXS spectra of correlated electron systems involve an excitonic intermediate state, even when the ground state electronic structure implies delocalized states, such as Hubbard bands in cuprates. Therefore, the RIXS process is non-selective on the wave vector $\mathbf{k}$ of the electron for which one wants to determine the EPC and the result is related to weighted integrals over the Brillouin zone. On the contrary, the RIXS experiment is intrinsically selective on the transferred momentum, which is supplied to the phonons excited in the final state. In particular, in experimental RIXS spectra, the one-phonon peak disperses and changes its intensity as the momentum transfer changes, following the phonon branch dispersion and the evolution of the EPC in the reciprocal space. Therefore, the use of the methods presented here at each momentum point is particularly well suited for optical phonon branches with little energy dispersion which correspond to good approximation to Einstein phonons. The method is as follows: for each transferred momentum, dictated by the energy of the photons and the experimental geometry, one can use one of the methods presented above or a combination of them to estimate $g$ or $M$. The same theoretical treatment is used to go from the measured quantities ($I_{1}$, $I_{2}$) to the EPC parameters. Eventually the momentum dependence of the EPC is thus obtained. We note that the procedure is safe because $\Gamma$ does not depend on the momentum and, $\omega_{0}$ is practically constant for optical oxygen phonons in the cuprates.

In the experimental results, it is important to separate a true variation of the EPC with momentum transfer from spurious effects like the dependence of atomic RIXS cross sections on the scattering geometry and photon polarization Devereaux et al. (2016). In our experiment, we used $\sigma$ incident polarization to accomplish this. As shown by measurements with polarization analysis of the scattered beam Fumagalli et al. (2019), the $\sigma\sigma^{\prime}$ channel (with incident and scattered $\sigma$ polarization) has an overwhelming intensity with respect to the $\sigma\pi^{\prime}$ channel. The ionic cross section for Cu $L_{3}$ resonant elastic scattering in the $\sigma\sigma^{\prime}$ channel is isotropic; moreover the $\sigma\pi^{\prime}$ cross section has negligible variation in our range of experimental momenta Moretti Sala et al. (2011); Fumagalli et al. (2019). For these reasons, we do not expect any sizable angular variation of the phonon intensity due to the resonant form factor. Therefore, 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 their ratio.

III.7 Application to other absorption edges

This work work is focused on the Cu $L_{3}$ edge due to our interest in understanding the role of the EPC in the high-$T_{\mathrm{c}}$ cuprates; however, this goal should not obscure the generality of our results and the great amount of work that remains to be done. The case of other absorption edges is briefly discussed here.

A situation very similar to the Cu $L_{3}$ edge is RIXS at the $L$ edges of other $3d$ transition metals in oxides, a vast and important class of materials. The levels involved and the nature of the cation-ligand bond are the same as in cuprates, so the mechanism governing the phonon excitation by RIXS is also the same. However, by going from Cu to Ti the available momentum decreases, the effects of electronic correlation weaken, and the lifetime of the intermediate state can change case by case. For example, the lifetime of the $L_{3}$ core hole in Ti4+ is longer than Cu by about a factor of 2.5 (tabulated $\Gamma=0.11$ eV for Ti, $\Gamma=0.28$ eV for Cu Keski-Rahkonen and Krause (1974); Krause and Oliver (1979)), resulting in much stronger peaks for single- and multi-phonon excitations.

The $L$ edges of heavier elements, namely $4d$ and $5d$ transition metals and lanthanides, are characterized by much larger values of $\Gamma$, i.e., shorter living intermediate states, whereas the $K$ edges of C, N, and O are narrower. By looking at Table 1, we reckon that phonon peaks can be very strong at C and O $K$ edges Lee et al. (2013); Johnston et al. (2016); Meyers et al. (2018), sizable at $3d$ transition metals $L_{3}$ edges, and hardly detectable at $L_{3}$ edges of heavier elements or at $K$ edges of any element heavier than chlorine ($Z=17$). The phonon signal in RIXS spectra at the $M_{4,5}$ edges of lanthanides is also expected to be very weak despite the relatively narrow line width because the EPC of the very localized $4f$ states is notoriously small.

IV Experimental results

IV.1 Overview and assignment of phonon modes

Figure 8(a) reports an overview of the momentum dependence of the low-energy RIXS features of NBCO-AF. Raw data (black solid lines) are plotted after subtraction of the elastic line. As anticipated, we only consider the momentum range of $|h|\geq 0.1$ r.l.u. and energy loss range $\lesssim 120$ – 140 meV to avoid the superposition of the phonon excitations with the rapidly dispersing magnetic excitations. The data are very expressive at a glance: even without optimizing the fitting, two main features are clearly detected with opposite momentum dependence of their intensity. These features are qualitatively plotted as red and blue Gaussian peaks, while the vertical dotted and dashed lines guide the eye across different spectra. While the present resolution allows us to separate these two features, it is insufficient to further decompose each into potential contributions from individual phonon branches that may be overlapping. Nevertheless, the structure of the phonon excitations observed here allows us to divide the excitations into two distinct energy regions, each containing modes that behave qualitatively different as a function of momentum. This observation enables us to assign the two features while keeping in mind that they may represent “effective branches” (i.e., the superposition of different phonons).

The behavior in high-energy region [50 – 100 meV, blue shaded area in Fig. 8(a)] is typical of the breathing (Cu-O in-plane bond-stretching modes) and apical oxygen modes. The EPC of the breathing mode increases on going from the $\Gamma$ to the $X$ point of the BZ, as theoretically shown in Ref. Johnston et al., 2010. Also, its intensity increases with momentum according to Ref. Devereaux et al., 2016 so that the intensity and the EPC go together. The increase of the intensity with momentum is very clear from Fig. 8(a), thus we assign the higher energy feature to the breathing modes. We note that the apical modes are not expected to influence appreciably the situation within our region of interest: the RIXS signal from the apical branch is expected to be about $1/3$ of that from the breathing branch and, even more importantly, its intensity is concentrated at small momenta outside of our working window Devereaux et al. (2016). In passing, we note that the separation of the apical and breathing modes in reciprocal space can be very useful in future RIXS experiments in nonlinear conditions at Free Electron Lasers, since the apical mode is often used to pump the sample in pump-probe experiments Hu et al. (2014); Kaiser et al. (2014).

The feature at lower energy [0 – 50 meV, red shaded area in Fig. 8(a)] is in the typical energy range of the out-of-plane polarized modes (e.g. the Cu-O bond buckling modes). The EPC of these phonons depends on both phonon and electron momentum so that a variety of situations can occur Johnston et al. (2010) (see Appendix A), including a decrease of the intensity with increasing momentum transfer which is, in a sense, a fingerprint of buckling branches. Thus, it is quite natural to assign the lower energy feature to buckling modes. This assignment is consistent with the strong static buckling in NBCO, which breaks inversion symmetry across the CuO2 plane and results in a larger EPC. We note that in the absence of such static buckling, one should expect a much weaker intensity for the buckling phonon modes; the case of Bi2212 is important with this respect, because the buckling mode intensity in experimental RIXS spectra is very small (below the present detection limit) and the static buckling is much smaller and with different local symmetry than in YBCO and NBCO. This family of phonon modes includes both in-phase (A1) and out-of-phase (B1) vibrations of the planar oxygen ions. In our case, the out-of-phase one is silent due to symmetry of the selection rules, while the in-phase one contributes to the RIXS signal (see the calculations in Appendix A). In the model by Devereaux et al. Devereaux et al. (2016), both branches are visible but the B1 is weaker by a factor $\sim 20$ with respect to A1, in agreement with our present assignment.

Further support for our assignments comes from the consistency between our observed phonon dispersion and the measured phonon dispersions in isostructural YBCO-AF Pintschovius and Reichardt (1998). This analogy is valid for our purposes, as one expects the phonon dispersions of the two compounds differ at mainly at very low energy, where the phonons involving the motion of Y and Nd are found Liu et al. (1988); Cardona et al. (1988); Yoshida et al. (1990); Limonov et al. (1998); Bohnen et al. (2003). The dispersion relations of the buckling and breathing branches of NBCO-AF are plotted in Fig. 8(b) as red circles and blue squares, respectively. The weak dispersion extracted from the fit of the data also provides post hoc justification for our use of a local phonon model.

IV.2 Determining the EPC in NBCO-AF

We now proceed with a detailed 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 $\Gamma/\omega_{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\approx 4$, which corresponds to $M=\omega_{0}\sqrt{g}\approx 0.13$ eV. The value of $\Gamma/\omega_{0}$ is in the range of validity of the SDR (Simplified Detuning Rule). It is remarkable that this rule gives $M=120$ meV, comparable with the above fitting. We note that from the fitted momentum dependence of the EPC of NBCO-AF reported by Rossi et al., one estimates $M\approx 0.16$ eV at $h=-0.4$ r.l.u. Rossi et al. (2019). By taking into account an uncertainty on the fitting results of about 10 – 15% (see Appendix B), the two values agree within error bars. We also note that the two different samples employed in the measurements may differ slightly in oxygen content, for the reasons stated in Sec. II.1. A slightly larger oxygen content in one of the samples may enhance the screening of the photoexcited electron, whose coupling to phonons may be consequently reduced.

Extracting the EPC strength for the buckling branch is more delicate and gives the us 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 squares), which is decomposed into an elastic line (orange solid line) and an inelastic signal (gray circles). The latter is further decomposed in Fig. 10(b) into the buckling (red dotted line) and breathing (blue dashed line) contributions each with their first overtone (orange dotted and light blue dashed lines, respectively).

The step-by-step decomposition of the inelastic spectrum is reported in Appendix B. We note that the robustness of the fit can be improved if we constrain the ratio $I_{2}/I_{1}$ for the breathing mode to the value corresponding to $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 beam time. At small 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 small $\mathbf{Q}_{\parallel}$. This is a typical case in which the direct connection between intensity and EPC does not work or works only in a subset of the parameter space. The final results on the EPC (in units of $\Gamma$) 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 (blue squares) shows the opposite trend.

IV.3 Doping effects

In the doped samples, one expects that the photoexcited electron and the EPC are 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) samples at $h=-0.4$ r.l.u. (Fig. 12). Note that this momentum transfer is well above the CDW wave vector ($\approx 0.31$ r.l.u. for NBCO Ghiringhelli et al. (2012)). 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 qualitatively mimics the doping dependence of the CDW signal, which is stronger in the UD sample and weaker in the OP sample Blanco-Canosa et al. (2014); Arpaia et al. (2019), as sketched in Fig. 12. The trend suggests that the CDW effect dominates over the screening even at $\mathbf{Q}_{\parallel}$ 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, which is evident from Fig. 13, where we compare the behavior at $h=-0.4$ r.l.u. (a, c) and $h=-0.1$ r.l.u. (b, d) of the NBCO-AF (black line) and NBCO-OP (red line) samples. 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 energy region of the phonon excitations. 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 NBCO-OP compound upon detuning by $\Omega=-1$ eV, whereas it has almost disappeared in the NBCO-AF sample (black squares). This fact 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 modes (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.

The momentum dependence of the EPC in the parent compound, plotted in Fig. 11, is fitted with simple trigonometric functions. The coupling to the buckling branch scales as $\cos(\pi h)$ (red dotted line) and the coupling to the breathing branch as $\sin(\pi h)$ (blue dashed line) (see Appendix A). Once again, we stress that these data may contain 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 result is important because the simple dependence of $M/\Gamma$ 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 scaling factor so that:

[TABLE]

These results demonstrate that RIXS provides easy and direct access to the momentum-dependent EPC. 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 useful 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 Devereaux and Einzel (1995); Andersen et al. (1996).

Until now, we have been using a strictly local perspective, where both the phonon and electron degrees of freedom are localized. While this approach produces good agreement with experimental data, as shown in Fig. 11, it is important to also consider what happens when itinerant description of the electron states is used instead. Interestingly, in the RIXS studies of EPC, the theoretical results turn out to be quite similar in the two approaches. The comparison is shown in Fig. 15, where the experimental data are also reproduced (the calculations are presented and compared in Appendix A). At least in this case the sensitivity to the electron states description is rather modest. This observation has profound implications on our work since it legitimates the use of the same theoretical tools in treating AF compounds and doped superconducting cuprates. The limits of applicability to other systems will have to be assessed with further work, since the results depend on a delicate combination of many factors, such as electronic screening, excitonic effects and time scales.

Before concluding, let us consider briefly the effect of CDWs. This is not the object of our research and we have concentrated the work on $\mathbf{Q}$ values far from the critical CDW wave vector. However, in doped cuprates it is known that the CDWs are ubiquitous, at least in the UD samples and with lower intensity in the OP compounds Blanco-Canosa et al. (2014); Comin and Damascelli (2016); Arpaia et al. (2019). 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 Bi2Sr2CaCu2O8+δ reported in Ref. Chaix et al., 2017 in that the increase of the intensity at $h=-0.4$ r.l.u. is analogous to the so-called phonon anomaly Le Tacon et al. (2013); Chaix et al. (2017).

• •

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 here, while the coupling of the breathing mode is close to zero. This situation as a whole is favorable to $d$-wave 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. Johnston et al., 2010 the buckling modes are beneficial to pairing, in particular the out-of-phase B1 mode, whose RIXS intensity vanishes for symmetry reasons (see Appendix A). Since the screening is comparable for the A1 and B1 modes, we argue that the out-of-phase mode couples to the electrons with a strength comparable to the in-phase mode. The effect of the breathing modes, instead, is now considered detrimental, although this has been controversial Bulut and Scalapino (1996); Shen et al. (2002); Ishihara and Nagaosa (2004); Sandvik et al. (2004); Johnston et al. (2010).

• •

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 Ghiringhelli et al. (2012); however, the recent discovery that in UD and OP NBCO and YBCO charge density fluctuations persist at all temperatures Arpaia et al. (2019) is consistent with the phenomenology of phonon intensity at high momenta. What remains as a question of interpretation is whether the enhancement of phonon peaks in RIXS is a signal of an actual increase of the EPC, which in turn causes the CDW, or whether the CDW has a different origin and the phonon enhancement at those momenta is the manifestation of the modulated charge density, at constant EPC Chaix et al. (2017); Lin et al. (2020). In all cases the sensitivity of the EPC to the CDWs is extremely high.

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_{\mathrm{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_{\mathrm{eff}}\sim 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_{\mathrm{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 three ways to obtain quantitative information on the EPC with momentum resolution from high energy resolution RIXS spectra. We applied the methods to antiferromagnetic NBCO. We also discussed the evolution of the EPC in underdoped and optimally doped NBCO and its relation to charge order. The three approaches were critically compared by pointing out their merits and limitations, and we showed that often the optimal way to go is a suitable combination of the three methods. Importantly, these approaches can be applied without numerical simulations of the spectra because we have introduced simple procedures based on rescaled variables and on energy detuning. Thus, they can be applied to a variety of strongly correlated materials. The present results have an impact not only on the “123” family of cuprates but also on material science in a broad sense. With the upcoming progress in the instrumentation the RIXS studies on electron-phonon coupling, at present in their infancy, will reach maturity.

Acknowledgements.

RIXS data were collected at the beam line ID32 of the ESRF using the ERIXS spectrometer jointly designed by the ESRF and the Politecnico di Milano. This work was supported by the ERC-P-ReXS project (2016-0790) of the Fondazione CARIPLO and Regione Lombardia, in Italy. R.A. is supported by the Swedish Research Council (VR) under the project “Evolution of nanoscale charge order in superconducting YBCO nanostructures”. S.J. is supported by the National Science Foundation under Grant No. DMR-1842056. J.v.d.B. acknowledges financial support from the DFG through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 39085490) and through SFB 1143 (project-id 247310070).

Appendix A Theoretical arguments on RIXS intensity and selection rules

We outline here a theoretical approach more general than the one introduced by Ament et al. Ament et al. (2011). We aim at validating the local electron model, which is implicit in our treatment. In fact, we demonstrate that the itinerancy of the electrons does not significantly change the landscape of the EPC in cuprates. We also derive the momentum dependence of the RIXS intensity for the buckling and breathing modes.

The generic form of the EPC is given by the following Hamiltonian:

[TABLE]

Here, $c_{\mathbf{k},\sigma}^{\dagger}$ ($c^{\phantom{\dagger}}_{\mathbf{k},\sigma}$) creates (annihilates) an electron with momentum $\mathbf{k}$ and spin $\sigma$ in the $pd$-antibonding band of the CuO2 plane, $b_{\mathbf{q},\nu}^{\dagger}$ ($b^{\phantom{\dagger}}_{\mathbf{q},\nu}$) creates (annihilates) a phonon in branch $\nu$ with momentum $\mathbf{q}$, and $M_{\nu}(\mathbf{k},\mathbf{q})$ is the EPC constant.

In the general case, the EPC depends on both the electron momentum $\mathbf{k}$ and the phonon momentum $\mathbf{q}$. The latter is related to the momentum transfer $\mathbf{Q}$. To determine how RIXS probes this momentum dependence, we consider the scattering process within the ultra-short core hole lifetime expansion. To first order, the RIXS intensity at the Cu $L_{3}$ edge is then given by:

[TABLE]

Here, the scattering cross section is defined as $\mathcal{F}_{fg}=\langle f|D_{\mathrm{out}}\mathcal{H}_{\mathrm{e-ph}}D_{\mathrm{in}}|g\rangle$, where $|g\rangle$ and $|f\rangle$ are the initial and final states, respectively, with energies $E_{g}$ and $E_{f}$, $\omega$ is the energy loss, and $D_{\mathrm{in}}$ and $D_{\mathrm{out}}$ are the dipole transition operators. The latter are defined as

[TABLE]

where $d_{i,\sigma}^{\dagger}$ creates an electron in the Cu $3d_{x^{2}-y^{2}}$ orbital located at position $\mathbf{R}_{i}$ and $p_{i,\sigma}^{\dagger}$ creates an electron in the relevant Cu $2p$ core level, and $\mathbf{q}_{\mathrm{in}}$ and $\mathbf{q}_{\mathrm{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\rangle=|\psi_{\mathrm{el}},n_{\mathbf{q},\nu}=0\rangle$, 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\rangle=|\psi_{\mathrm{el}},n_{\mathbf{q},\nu}=1\rangle$. 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:

[TABLE]

In the last step we have introduced the momentum transfer ${\bf Q}={\bf q}_{\mathrm{in}}-{\bf q}_{\mathrm{out}}$.

The Cu orbital operator $d_{\mathbf{k},\sigma}^{\dagger}$ is related to the $c_{\mathbf{k},\sigma}^{\dagger}$ band operator by $c_{\mathbf{k},\sigma}^{\dagger}=\phi_{\mathrm{Cu}}^{\star}(\mathbf{k})d_{\mathbf{k},\sigma}^{\dagger}$, where $\phi_{\mathrm{Cu}}^{\star}(\mathbf{k})$ measures the Cu character of the antibonding band. Therefore, the scattering matrix element written in band space is:

[TABLE]

The expectation value appearing in Eq. (13) 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 $\mathbf{Q}$, we require $\mathbf{p}=\mathbf{k-q}$ and $\mathbf{Q}=\mathbf{q}$, and the scattering amplitude simplifies to:

[TABLE]

Here, $n_{\mathrm{F}}(x)$ is the Fermi factor and $\epsilon({\bf k})$ is the band dispersion. Similarly, we obtain the localized limit by inverse Fourier transforming Eq. (13) and then retaining only the local operators inside the expectation value:

[TABLE]

Thus, in both limits, RIXS at the Cu $L_{3}$ edge measures a $\mathbf{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 Lee et al. (2013). 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 Johnston et al., 2010 argued that $\phi_{\mathrm{Cu}}(\mathbf{k})$ can be approximated by a constant $\phi_{\mathrm{Cu}}(\mathbf{k})\sim A_{\mathrm{Cu}}$, which leads to simplified momentum dependencies for the various coupling constants. In this case, $A_{\mathrm{Cu}}$ can be absorbed into the prefactor of the coupling constant and the scattering amplitude for the localized limit reduces to $F_{fq}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k},\sigma,\nu}M^{\phantom{\star}}_{\nu}(\mathbf{k,Q})$. For the in-phase ($+$) and out-of-phase ($-$) Cu-O bond-buckling branches one has ($\mathbf{p=k-Q}$):

[TABLE]

Similarly, for the breathing mode:

[TABLE]

Because the coupling constant for the breathing mode does not depend on the Fermion momentum $\mathbf{k}$, the $\mathbf{k}$-averaging is trivial and $I_{\mathrm{breath}}(\mathbf{Q})\propto|M_{\mathrm{breath}}(\mathbf{Q})|^{2}$. The situation is different for the buckling modes. Due to the sign change in the $\mathbf{k}$-dependence of $M_{-}(\mathbf{k,Q})$, the intensity for the out-of-phase buckling mode $I_{-}(\mathbf{Q})$ vanishes. For the in-phase buckling phonon branch one arrives at:

[TABLE]

Figure 15 of the main text displays the EPC strength $M(\mathbf{Q})$ of the in-phase buckling (red circles) and breathing (blue squares) modes, normalized to the natural width $\Gamma$ 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\pi$). 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. Markiewicz et al., 2005. 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.

Appendix B Recovering the phonon signal from RIXS spectra

We summarize here the procedures used to recover the phonon components from the RIXS spectra. These rather technical issues are collected here to improve the flow and cohesion of the main text. Nevertheless, the content of this Appendix must be considered as an integral part of the work.

We consider the benchmark RIXS spectrum shown in Fig. 16(a) (black solid line) and we describe its decomposition. For completeness, we show in panel (a) the subtraction of the long and weak tail coming from the magnon peak. This tail is dominated by the Lorentzian shape of the magnon (purple dash dotted line), which is broader than the instrumental line width. Strictly speaking, a Voigt function should be used but the Lorentzian shape is a good approximation. The resolution-limited Gaussian lineshape defining the elastic line is shown as orange dashed line. The resulting inelastic spectrum is plotted as gray circles in panel (b).

We must use the appropriate line shape when decomposing the inelastic phonon region. Since the instrumental line width is much larger than the phonon width, a resolution-limited Gaussian seems to be acceptable. However, as discussed in the main text, the phonon features that we observe may be an average of several phonon modes. Thus, a broadening is generated with respect to the instrumental response function. Our numerical analysis shows that this is a small effect (typically 2 – 3 meV). To take this effect into account, we leave the line width as a free parameter between 34 meV and 37 meV. Another issue is the energy of the overtone, which in the Einstein model is exactly twice the phonon energy. However, such an extreme situation is very rare in a real material and we use the Einstein model as an approximate description of a system having a tiny dispersion. The factor of two between the phonon energy and its overtone is an upper limit so that also this parameter is left free.

With the above precautions in mind we address the decomposition of the inelastic spectrum by exploiting first the experimental information that is better observed. This is done iteratively with the following logical sequence. Initially, we fit the low-energy shoulder of the buckling mode as shown in panel (b) (dotted red line). The next step is the addition of the breathing component at an energy that is very well defined already in the raw data [dashed blue line in panel (c)]. Moreover, the tail at high energy suggests the approximate position of the breathing overtone (dashed light blue line). The intensity of the main peak and of the overtone of the breathing mode are chosen so that their sum reproduces roughly the turning point of the inelastic signal at around 100 meV [solid line in panel (d), where the energy axis is zoomed into the range of interest]. This is a delicate issue we will discuss in more detail in the next paragraph. By subtracting the two main phonon peaks from the inelastic spectrum, a contribution coming only from the overtones emerges [black squares in panel (e), where both axes are zoomed into the appropriate region]. The black solid line of panel (e) is the smoothed total signal from the overtones. The difference between the black line and the light blue dashed line gives the orange dotted curve representing the buckling overtone. At this time a decomposition cycle is completed, but it is certainly necessary to repeat the procedure to refine the parameters.

The final result is given in panel (f) as a dark yellow line. The residuals are plotted as a thin green line and are very small. The weak point of this procedure comes from the difficulty in finding the breathing weights in panels (c) and (d). As a matter of fact, if this value is changed by 10%, it is still possible to obtain a good decomposition so that there is some ambiguity in the decomposition. One can improve the procedure, however, if information on $I_{1}$ or $I_{2}/I_{1}$ is exploited, which can be obtained using another method. This is what we have done on the basis of the results from detuning. Once again a suitable mix of the different approaches improves the quality of the decomposition. Note that the energies of the overtones differ from the value of the Einstein model, as expected: for the breathing mode the ratio of the overtone energy and the main peak energy is 1.7 instead of 2, while for the buckling mode it amounts to 1.55.