Bosonic excitation spectra of superconducting $\mathrm{Bi_2Sr_2CaCu_2O_{8+\delta}}$ and $\mathrm{YBa_2Cu_3O_{6+x}}$ extracted from scanning tunneling spectra

From a scanning tunneling spectrum below T_c we obtain the differential conductance d$I/$d$U(eU)\equiv \sigma^\mathrm{tot}(eU)$. This function consists of the purely elastic part $\sigma^\mathrm{el}$ and the inelastic part $\sigma^\mathrm{inel}$. The elastic part is directly proportional to the superconducting density of states in the sample ν_s . In this step, we determine the functional form of the elastic contribution by fitting a model function to the low bias region of the d$I/$dU spectrum that keeps the complexity as low as possible while still capturing the relevant features of the band structure and pairing strength. We opt for a generalized Dynes function [58] with a momentum dependent gap:

$$\begin{align} \sigma^\mathrm{el}_s\left(\omega\right)& = \frac{\sigma_0}{\mathcal{N}}\sum\limits_\lambda C_\lambda\int^{\pi/2}_0 \mathrm{d}\varphi \nu^\mathrm{F}_\mathrm{n,\lambda}\left(\varphi\right) \nonumber\\ & \quad \times \left|\Re\left(\frac{\omega+{i\mkern1mu} \Gamma_\lambda\left(\varphi\right)}{\sqrt{\left(\omega + {i\mkern1mu} \Gamma_\lambda\left(\varphi\right)\right)^2-\Delta_\lambda\left(\varphi\right)^2}}\right)\right| \,. \end{align} \tag{ 3 }$$

Here, $\mathcal{N} = \sum_\lambda\int \mathrm{d}\varphi \nu^\mathrm{F}_\mathrm{n,\lambda}(k_{\mathrm {F}},\varphi)$ is a normalization factor, $\Delta_\lambda(\varphi) = \Delta_{0,\lambda}\cos(2\varphi)$ is the d-wave pairing potential, $\Gamma_\lambda(\varphi) = \gamma_\lambda\cdot |\Delta_\lambda(\varphi)|$ the quasiparticle scattering rate, $\nu^\mathrm{F}_\mathrm{n}(\varphi)$ is the angle-dependent DOS at the Fermi energy in the normal conducting phase, λ the band index and C_λ the relative tunneling sensitivity for the band. The function $\nu^\mathrm{F}_\mathrm{n}(\varphi)$ weighs gap distributions for different $(k_x,k_y)$ by their abundance along the Fermi surface (FS). It is derived from a microscopic tight-binding approach that models the dispersion relation. In the case of Bi2212 we used a single-band model whereas for Y123 we considered two CuO₂ plane bands and one CuO chain band (see appendix A). It should be noted that, unlike in fully gapped superconductors, the inelastic spectrum can be non-zero down to vanishing bias voltage because $\Delta(\boldsymbol{k})$ has a nodal structure. This prevents us from directly assigning the differential conductance for $e|U|\lesssim \Delta$ to the purely elastic tunneling contribution as was possible for the $s_\pm$ superconductor monolayer FeSe [56, 57]. We will, however, start from here and refine $\sigma^\mathrm{el}$ in the next step.

We use the fact that $\sigma^\mathrm{tot}(\omega) = \sigma^\mathrm{el}(\omega)+\sigma^\mathrm{inel}(\omega)$ and the physical constraints $\sigma^\mathrm{inel} (0) = 0$ and $\sigma^\mathrm{inel}(e|U|\gt0)\gt0$. In order to relate our fitted curve to $\sigma^\mathrm{el}$, we need to make an initial guess for $\sigma^\mathrm{inel}$ in the range $0\lt e|U|\lt\Delta$. As a first guess, we assume $\sigma^\mathrm{inel}$ to be directly proportional to $\nu_s(\omega)$. This is equivalent to assuming that $B(\Omega) \sim \delta(\Omega - 0)$ and we can rewrite for that energy regime:

$$\begin{align} \sigma^\mathrm{tot} = \sigma^\mathrm{el}\left(1+\mathrm{const.}\right) = \sigma^\mathrm{el}\left(1/\eta\right)\,,\quad 0 < \eta < 1. \end{align} \tag{ 4 }$$

In fact, we rather expect $B(\Omega)$ to be a slowly increasing function with a soft spin gap for $0\lt e|U|\lt\Delta$ due to the ungapped parts of the Fermi surface, but our choice of $B(\Omega)$ is only a first guess and will later be refined. Additionally, the delta distribution at zero is far away from the physically interesting energy regime, which avoids a misinterpretation of the final extracted boson spectrum. The introduction of the numerical scaling factor η is also necessary from a computational standpoint as it can ensure that $\sigma^\mathrm{inel}$ and $\sigma^\mathrm{el}$ are always positive, which is a prerequisite for our extraction procedure. We approximated η using a boundary condition for the total number of states (see appendix B) and in order to keep the condition $\sigma^\mathrm{el}(e|U|\rightarrow D) \rightarrow \sigma_0$ (D denotes the energy scale of the band edge), we scale up our experimental curve by $1/\eta$ instead of scaling down our fitted curve. $\sigma^\mathrm{inel}$ is then simply the difference of the scaled up (by $1/\eta$) experimental data points and the Dynes fit ($\sigma^\mathrm{el}$). We took care that the choice of this numerical factor, that simply helps to perform our deconvolution algorithm and paint a more realistic picture of $\sigma^\mathrm{el}$, does not influence the extracted boson spectrum in a qualitative manner (see appendix B).

We compare two methods by which the boson spectrum was determined: The direct deconvolution in Fourier space and the iterative Gold algorithm [59, 60] to perform the deconvolution. The advantage of the Gold algorithm is that for positive kernel and signal, the result of this deconvolution method is always positive. This is in agreement with our physical constraint that the bosonic excitation spectrum is strictly positive. We used the implementation of the one-fold Gold algorithm in the TSpectrum class of ROOT system [61, 62] in C language wrapped in a small python module.

The bosonic function from direct deconvolution in Fourier space is obtained from

$$\begin{align} B\left(\Omega > 0\right) = \mathcal{F}^{-1}\left(\frac{\mathcal{F}\left(\sigma^\mathrm{inel}\right)\left(t\right)}{\sqrt{2\pi}\mathcal{F}\left(\nu_s\cdot \Theta\right)\left(t\right)}\right) \end{align} \tag{ 5 }$$

where $\mathcal{F}$ denotes a Fourier Transform and $\mathcal{F}^{-1}$ the inverse transform. The abrupt change in elastic conductance at zero energy (multiplication with Heaviside distribution in equation (2)) leads to heavy oscillations in the Fourier components. Therefore the result of this deconvolution procedure can contain non-zero contributions for E < 0 and negative contributions for $0\lt E\lt\Delta$. They are exact solutions to equation (2) but from a subset of non-physical solutions that we are not interested in. Because the solutions obtained in this way are highly oscillatory we show the result after Gaussian smoothing. In order to obtain a positive solution to equation (2), the result of the direct deconvolution method is used as a first guess to the Gold algorithm. The results shown in this work are obtained after 2000 000 iterations at which point convergence has been reached.

In order for our extraction method to be applicable, several simplifying assumptions were made:

• 1.  
  Quasiparticles with energy ω couple to bosonic excitations of energy Ω and effective integrated density of states $B(\Omega)$. The k-dependence of the interaction is thus disregarded.
• 2.  
  $\sigma^\mathrm{inel}$ has a simple relation to $\sigma^\mathrm{el}$ for $0\lt e|U|\lt\Delta$ (see section 2.1.2) which is generally oversimplified for d-wave superconductors, especially near Δ.

In general, retrieving the source function from a convolution integral is an ill-posed problem which means that we only obtain one solution from a large set of valid solutions to the convolution equation. Additional aspects that complicate the problem are the following:

• 1.  
  The kernel function $\nu_s(\omega)$ is a guess which is very dependent on the modeling of the superconducting density of states and is further questioned by lack of energetic regions of purely elastic processes in the scanning tunneling spectrum. We are in fact on the verge of a necessity for blind deconvolution algorithms.
• 2.  
  Strong electron-boson coupling leads to spectral features of $\sigma^\mathrm{el}$ outside the gap [52, 63, 64] that are neglected here as the contribution of $\sigma^\mathrm{inel}$ is expected to be much larger. They can in principle be reconstructed within an Eliashberg theory using the extracted boson spectrum. Using this refined $\sigma^\mathrm{el}$ and repeating the procedure until $B(\Omega)$ leads to the correct $\sigma^\mathrm{el}$ and $\sigma^\mathrm{inel}$ could further improve our result.
• 3.  
  Electronic noise in the recorded spectra is ignored and consequently ends up in either ν_s or $B(\Omega)$.
• 4.  
  With $B(\Omega)$ we obtain only an 'effective tunnel Eliashberg function' which includes all bosonic excitations that are accessible to the tunneling quasiparticles. Hence, no disentanglement into lattice and spin degrees of freedom is possible.
• 5.  
  The k-dependence of $B(\Omega)$ is inaccessible.

Figure 1 shows experimental $\mathrm{d}I/\mathrm{d}U$ and $\mathrm{d}^2I/\mathrm{d}U^2$ spectra recorded at 0.7 K and 84 K. In order to remove a linear background that stems from a slope in the DOS due to hole doping [67, 68] and prepare the data for further analysis of particle-hole symmetric features, we follow a previous work [57] and symmetrized the $\mathrm{d}I/\mathrm{d}U$ spectra in figure 1(a). This procedure does not affect the quasi-particle spectrum due to particle-hole symmetry of the superconducting state and the inelastic contribution to the tunneling current. The dI/dU spectrum for superconducting Bi2212 in figure 1(a) shows a single but smeared gap with residual zero-bias conductance due to the nodal $d_{x^2-y^2}$ gap symmetry. Similarly, the coherence peaks are smeared due to the gap symmetry and possibly also due to short quasiparticle lifetimes. This is typical for the underdoped regime and may be caused by its proximity to the insulating phase [66]. Outside the gap, a clear dip of the superconducting spectrum below the normal conducting spectrum, followed by a hump reapproaching it, are visible. The V-shaped conductance in the normal state hints towards strong inelastic contributions to the tunneling current from overdamped electronic excitations discussed in context of cuprates and iron-based superconductors in [56, 69, 70]. A similar V-shape due to inelastic scattering off phonons has even be reported for conventional superconductors like Pb [41]. In contrast to phonons, electronic excitations become partly gapped in the superconducting state causing a redistribution of the spectrum in form of dip-hump features in the differential conductance. The hump shows as a peak in the second derivative of the tunneling current that exceeds the curve of the normal state at $\approx {120}~\textrm{mV}$ in figure 1(b). The relatively round shape of the superconducting gap in Bi2212 is atypical for a classic d-wave superconductor, in which the naive expectation is a V-shaped conductance minimum. As will be shown later on, the round shape of the gap can be generated without admixture of an s-wave pairing term by respecting the anisotropy of the Fermi surface in the normal state. The Fermi surface and gap anisotropy are summarized schematically in figure 2(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We followed the step-by-step extraction procedure outlined in section 2.1 starting from the determination of the superconducting density of states ν_s . The optimal Dynes fit to our experimental spectrum in the superconducting state is shown in figure 2(b) with $\Delta_0 = {63.31}~\textrm{meV}$ ($\Delta_\mathrm{max} = {59.14}~\textrm{meV}$ and $\bar{\Delta} = {49.60}~\textrm{meV}$) and $\gamma = {0.15}{}$. The resulting (in)elastic contribution is shown in red(grey) in figure 2(c). Here, a numerical scaling factor of η = 0.6 was used. The value for Δ lies within the range of previously reported gap values on the Bi2212 surface [66], especially in the slightly underdoped regime, where variations of the local gap from the average gap tend to be larger [72].

The regularized bosonic function from direct deconvolution in Fourier space is shown in figure 2(d) in orange. The contributions at low energies are an artifact from the scaling with factor η. Despite our uncompromising simplifications, the bosonic spectrum recovers well the tendency of the total conductance in the forward convolution (figure 2(c) orange) and shows the expected behavior for coupling to spin degrees of freedom at medium and high energies, i.e. a resonance mode at $\Delta \lt E \lt 2\Delta$ and approach of the normal state bosonic function B for $E \gtrsim 3\Delta$ [56], that, in contrast to the Eliashberg function in the case of phonon-mediated pairing, remains finite for energies well above 2Δ. While the long-lived resonance mode is associated with a spin resonance due to the sign-changing gap function, the broad high-energy tail of the bosonic spectrum is due to the coupling to overdamped spin fluctuations, or paramagnons [18]. Resonant inelastic x-ray scattering (RIXS) studies have shown that these paramagnons dominate the bosonic spectrum for energies larger than ${\approx}{100}~\textrm{meV}$ in several families of cuprates as almost all other contributors, e.g. phonons, lie lower in energy [73–76].

By application of the Gold algorithm we obtained the bosonic spectrum shown in green in figure 2(d). Again, the high value at E = 0 is a consequence of the scaling with factor η. We could get rid of negative contributions and find a bosonic function that recovers well the total conductance (figure 2(c) green), especially the dip-hump structure, and shows a very clear resonance at $\Omega_\mathrm{res} \approx {63}~\textrm{meV}\approx 1.0 \Delta_0 \approx 1.1 \Delta_\mathrm{max}\approx 1.3 \bar{\Delta}$. The sharp peak at ${10}~\textrm{meV}$ is due to inadequacies of our elastic fit in the region of the coherence peak. It is e.g. not present in our analysis of Y123 (see section 3.2) and vanishes once we take an elastic DOS with a larger gap (here ${65}~\textrm{meV}$, not least-square minimum) as we show in the thin dark green line in figure 2(d). This is in contrast to the resonance at $\Omega_\mathrm{res}$, which remains and only decreases slightly in energy for the fit with a larger gap. This behavior is expected according to the hump in the tunneling spectrum at $\sim\Delta+\Omega_\mathrm{res}$.

The resonance mode extracted in this work is higher in energy than reported in INS experiments ($\Omega_\mathrm{res} \approx {43}~\textrm{meV}$ at the antiferromagnetic ordering vector) [12] and closer to the resonance determined by optical scattering ($\Omega_\mathrm{res} \approx {60}~\textrm{meV}$) [32]. Due to the loss of momentum information in tunneling, the center of the resonance is expected to be shifted to higher energies compared to the INS results [56]. Due to the large inhomogeneity of Δ on the surface of Bi2212 [66], it is more instructive to compare the ratio $\Omega_\mathrm{res}/\Delta$ to other works rather than the absolute value of $\Omega_\mathrm{res}$. The ratio $\Omega_\mathrm{res}/\bar{\Delta} = 1.3$ lies within the current range of error of $\Omega_\mathrm{res}/\Delta = 1.28\pm 0.08$ by Yu et al [25]. In most other extraction methods of the bosonic mode energy, the normal state DOS is not respected, which is why, depending on the method, the Δ₀ used there is most similar to what is here called $\Delta_\mathrm{max}$ or $\bar{\Delta}$. $\Delta_\mathrm{max}$ is the largest gap value that contributes to the elastic conductance spectrum and $\bar{\Delta}$ the momentum averaged and DOS weighted gap.

The symmetrized dI/dU spectrum for superconducting Y123 in figure 3(a) is qualitatively in excellent agreement with previous STM measurements [77–82] and shows three low-energy features: (i) a superconducting coherence peak at ${\approx}{25}~\textrm{meV}$ that is sharper than in Bi2212, (ii) a high-energy shoulder of the coherence peak and (iii) a low-energy peak at ${\approx}{10}~\textrm{meV}$. The high-energy shoulder as well as the sub-gap peak are believed to arise from the proximity-induced superconductivity in BaO planes and CuO chains [80, 83–85]. This would certainly account for the fact that these states are missing in the Bi-based compounds and that the sub-gap peak shows a direction-dependent dispersion in ARPES data [86, 87]. At energies larger than Δ, we again find a clear dip-hump feature, similar as in Bi2212. The hump lies at ${\approx}{60}~\textrm{meV}$ as can also be seen from the second derivative of the tunneling current in figure 3(b). The V-shaped background conductance in the high-energy regime of the superconducting spectrum is in agreement with the predicted inelastic contribution by magnetic scattering in the spin-fermion model [56, 69].

Zoom In Zoom Out Reset image size

We proceeded as in the case for Bi2212, but incorporated the one-dimensional band from the CuO chains as well as the bonding and anti-bonding band from the CuO₂ planes into the calculation of the normal state DOS to remodel the sub-gap peak and coherence peak shoulder in the estimated $\sigma^\mathrm{el}$ of the superconducting state. The optimal Dynes fit with gaps $\Delta^\mathrm{AB}_0 = {20.62}~\textrm{meV}$ for the anti-bonding (AB), $\Delta^\mathrm{BB}_0 = {25.77}~\textrm{meV}$ for the bonding (BB) and $\Delta^\mathrm{CHSS}_0 = {5.66}~\textrm{meV}$ for the chain band is shown in red in figure 3(a). Because vacuum-cleaved surfaces favor tunneling into states of the CuO chain plane [77, 80, 88], the sub-gap peak is pronounced and the contribution to the total DOS of the CH$_\mathrm{SS}$ band is, in our analysis, roughly five times higher than for the AB and BB band. The size of $\Delta^\mathrm{BB}$ is in good agreement with other scanning tunneling spectroscopy results spanning around 20 experiments, in which the extracted gap value lies between $\Delta = 18-{30}~\textrm{meV}$ for optimally doped samples [66]. For comparison: From Raman spectra, Δ₀, i.e. the gap in the antinodal direction, is frequently found to be ${34}~\textrm{meV}$ for optimally doped Y123 samples [89–92]. It should be noted that vacuum cleaved surfaces of Y123 tend to be overdoped [77] which goes hand in hand with a steep decline of Δ. The reason for discrepancy between the gap measured in STS and ARPES [93, 94] (also yielding $\Delta_0 \sim {34}~\textrm{meV}$) is expected due to two factors: (i) Although less influenced by a local gap variation than Bi2212, the gap of Y123 is expected to be inhomogeneous on a wider scale of ${\gt}{100}~\textrm{nm}$ [93]. While ARPES yields an average gap over several of these domains, STS yields a more local gap. (ii) The measurement of a k-averaged gap value in STS naturally tends to give smaller values for a d-wave superconductor than the maximum gap size measured in ARPES. We try to eliminate this last effect by respecting the k-dependence of Δ and $\nu_\mathrm{n}^\mathrm{F}$ in our fit. Nevertheless, despite the large T_c , the spectroscopic results on Y123 in this work do not support an effective gap value of ${\gt}{30}~\textrm{meV}$ because the total conductivity is already on the decrease at this energy.

Analogous to the case of Bi2212, the elastic part was, as a first guess, approximated by the Dynes fit to the total conductance times a scalar factor η. Here, η = 0.9 was chosen in order to secure the constraint $\sigma^\mathrm{inel}(e|U|\gt0)\gt0$. The (in)elastic parts to the total conductance are shown in red(grey) in figure 3(c).

We compare the extracted bosonic DOS obtained from direct deconvolution and Gold algorithm for Y123 in figure 3(d). The resonance mode at $\Omega_\mathrm{res}\approx {61}~\textrm{meV}$ is significantly higher in energy than experimentally found by INS in (nearly) optimally doped samples with $\Omega_\mathrm{res}\approx {41}~\textrm{meV}$ [11, 95–97] and even lies at the onset of the spin scattering continuum at $\Omega_c \approx {60}~\textrm{meV}$ [97]. Apart from the k-space integration, which shifts the peak center to higher energies, several other factors can play a key role: (i) The well-studied ${41}~\textrm{meV}$ odd-parity mode is paired with an even-parity mode at $\Omega^\mathrm{e}_\mathrm{res}\approx 53-{55}~\textrm{meV}$ [97–99] which may be of the same origin as it vanishes at T_c . This mode appears with a ${\approx}3-20$ times lower intensity in INS than the odd-parity mode, but this does not necessarily have to hold for a tunneling experiment. (ii) The bosonic spectrum extracted here is essentially poisoned by phononic contributions from every k-space angle. A disentanglement of phononic and electronic contributions to the total bosonic function by non-equilibrium optical spectroscopy showed that for $\Omega\gt{100}~\textrm{meV}$ the bosonic function is purely electronic, yet in the energy range of the spin resonance, the contribution of strong-coupling phonons is almost equal to that of electronic origin [8]. (iii) Apart from physical arguments, there can also be made sceptical remarks on the deconvolution procedure: Evidently, it heavily depends on the guess of the elastic tunneling conductance, which in this case does not contain strong-coupling features from an Eliashberg theory. (iv) The pronounced contribution of the CuO chains to the total conductance essentially causes the resonance mode to appear at roughly $\omega_\mathrm{hump}-\Delta^\mathrm{CHSS}$ instead of $\omega_\mathrm{hump}-\Delta^\mathrm{BB}$. Correcting for the 20 meV difference between the two gaps, it is likely that without sensitivity to the CuO chain gap, our extraction procedure will yield $\Omega_\mathrm{res}\approx {41}~\textrm{meV}\approx 1.6\,\Delta^\mathrm{BB}_0 \approx 1.8\,\Delta^\mathrm{BB}_\mathrm{max}\approx 2.4\,\bar{\Delta}^\mathrm{BB}$.