Correlation between optical phonon softening and superconducting $T_\mathrm{c}$ in YBa₂Cu₃O_x within d-wave Eliashberg theory

High-temperature cuprates superconductors were discovered by Bednorz and Müller [1] under the guiding assumption that strong Jahn-Teller type lattice distortions, akin to those in ferroelectric materials, would lead to large polarons and very strong electron–phonon couplings. Over the subsequent decades much new physics in these materials has been unveiled, where non-trivial magnetic effects accompany the superconducting transition, but a microscopic theory is still lacking [2]. Several arguments have been put forward supporting an unconventional pairing mechanism in the cuprates, among which the strongest arguments are perhaps the fluctuating antiferromagnetic ordering, suggesting electron-electron interactions [3], and the strongly 2D physics within the layers [4].

In parallel with research on possible unconventional mechanisms, some experimental evidence still points towards a large electron–phonon coupling [4–6] and the isotope effect [7, 8], possibly enhanced by electronic and magnetic correlations [9, 10], such that suitably modified Bardeen–Cooper–Schrieffer scenarios may be able to account for the observations in doped cuprates, without the need of exotic physics. For an up-to-date discussion of these issues, see e.g. [11, 12].

In this work, we add further evidence in favour of an underlying dominant electron-lattice coupling, by proving a description of the superconducting $T_\mathrm{c}$ that, crucially, incorporates the optical phonon softening as measured experimentally. In particular, earlier experimental results based on Raman scattering by Liarokapis, Müller, Kaldis and co-workers showed that the in-plane A_g oxygen mode displays a dramatic softening upon increasing the oxygen doping x in the model cuprate system YBa₂Cu₃O_x [13, 14]. Concomitantly, the superconducting $T_\mathrm{c}$ is experimentally observed to grow in a correlated fashion with the increase of doping x up to a maximum or dome. This correlation has remained thus far as a mere observational experimental fact, in spite of the crucial importance of oxygen doping for high-T superconductors performance [15]. Furthermore, very recent results based on x-ray absorption spectra support a striking correlation between anharmonic lattice displacements of Sr atoms (via their coupling to apical oxygens) and the superconducting critical temperature in YSr₂Cu_2.75Mo_0.25O_7.54 and Sr₂CuO_3.3 [16]. Again, this kind of correlations cannot be explained within unconventional (e.g. spin fluctuations or electron repulsion) mechanisms.

We build on the recent observation that sharp optical (or acoustic) phonon softening implemented within Migdal–Eliashberg theory can lead to huge electron–phonon couplings and consequent order-of-magnitude increase in $T_\mathrm{c}$ [17]. This generic theoretical prediction has been recently verified by means of first-principles computations on the example of a new compound, of Zr₂S₂C [18]. It was observed that, as predicted by the theory [17], increasing the depth of the acoustic phonon softening can increase the electron–phonon coupling λ by a factor 2 [18].

Implementing the experimentally observed softening of the A_g mode as a function of doping x within the d-wave Migdal–Eliashberg scheme, we mathematically describe the correlation between the oxygen doping, the increase of the softening depth (and sharpening of its width) and the concomitant non-trivial variation of the $T_\mathrm{c}$ with doping, in good qualitative agreement with experimental data. In the comparison between theory and experiments, three free fitting parameters are used: one is the electron–phonon matrix element C which controls the quantitative value of $T_\mathrm{c}$. The other fitting parameter is the shape of the optical phonon softening. Both parameters vary with x. The evolution of the softening with the doping level x is found to be physically meaningful and in qualitative agreement with typical momentum-resolved experimental observations of optical softening in the cuprates [19]. The dependence of the electron–phonon matrix elements on the doping x will be important to capture the correct order of magnitude of $T_\mathrm{c}$ and will be discussed at length in the following. Finally, the third free parameter is the softening wave-vector, which, as we will see, does not affect significantly the results and will be therefore taken in a first approximation as doping-independent (see below for a more detailed discussion on this point).

Numerous experimental studies have established that the plateau in the $T_\mathrm{c} - x$ (oxygen doping) phase diagram corresponds to the suppression of $T_\mathrm{c}$ induced by charge order, charge density wave (CDW) and charge density fluctuations (CDFs) [20]. Also, there are many works confirming that the CDW can give rise to phonon softening at a CDW wavevector q_CDW [21, 22]. For YBCO systems, it has been confirmed that the intensity of CDF reaches its maximum at near optimal doping x = 7, where the phonon is softened to the lowest frequency [23]. There is a strong correlation between the in-phase oxygen mode frequency and the characteristic doping levels of YBa₂Cu₃O_x, shown in [24]. This correlation was discussed in terms of the superstructures developed from the chain oxygen atoms, but it is, in fact, related to the charge ordering. As a matter of fact, a description of these softened-phonon effects on the $T_\mathrm{c}$, within the conventional phonon mediated superconducting framework, is still missing. The simultaneous interplay between charge order, anharmonic phonons and superconductivity still poses an enormous challenge to microscopic theories [25–27]. Both older theoretical and experimental works [8, 28], as well as very recent experimental evidence [16], suggest that phonons cannot be overlooked for explaining the superconductivity in doped cuprates. Our results provide further evidence to the role of electron–phonon coupling for the physical origin of high-temperature superconductivity in oxygen doped rare-earth cuprate oxides.

We choose an Occam's razor strategy and work with a minimal mathematical model in order to: (i) reduce the number of unknown adjustable parameters and (ii) unveil direct connections between experimentally observable trends. Also, working with an effective model allows us to accommodate the experimental input from Raman scattering experiments, which would not be possible within a fully ab − initio approach. In turn, this allows us to focus on physical trends rather than on absolute values of physical quantities, such as $T_\mathrm{c}$, which are anyway out of reach even for the most detailed microscopic calculations. It should be noted that the Raman scattering data that we used [13, 14, 29] were obtained at room temperature, for many values of oxygen content, several (7–8) micro-crystallites and long acquisition times to secure good statistics. Data obtained at low liquid nitrogen temperatures (T = 77 K) show the same characteristics as the room temperature ones, though for high oxygen concentrations the data were collected below $T_\mathrm{c}$. We show also the Raman measurements taken at low temperature (T = 77 K) in the appendix. Since these data present exactly the same features of phonon softening observed at room temperatures, we use the room temperature data as input for the theoretical model, given their better statistics.

The strategy behind our theoretical modeling is as follows:

• (i)  
  The CDW softens the phonon dispersion relation at $q_{\text{soft}} = q_\textrm{CDW}$. There are three parameters controlling the properties (shape) of the softened dispersion relation: the depth of softening δ, which can be determined from experimental data as shown in figure 1(b), the width of the softening region, β, which will be discussed more in detail as an adjustable parameter, and the softening wave vector q_soft which is, usually, fixed. The linewidth of the phonon is also taken directly from experimental data as shown in figure 1(a). The effective dispersion relation and the effect of depth, width, and position of the softened dispersion relation can be seen in the inset of figure 2(a).
• (ii)  
  Once that the softened phonon dispersion relation has been determined, the Eliashberg function $\alpha^2 F(\omega)$ can be obtained by using the method described in [17].
• (iii)  
  By implementing the Eliashberg function $\alpha^2 F(\omega)$ into the simplified Eliashberg equation for the electron–phonon coupling λ [30], $T_\mathrm{c}$ can be obtained via the Allen–Dynes formula, in an approximate way.
• (iv)  
  By direct numerical solution of the d-wave Eliashberg equations, we can verify the T_c trend versus doping x obtained in the point (iii) and produce a physical interpretation of the results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, our strategy will be to build an effective phonon-softening dispersion relation with input from experimental data, and to use this effective dispersion relation to compute the superconducting $T_\mathrm{c}$ within the conventional phonon-mediated framework. The charge order results in phonon softening at $q_{\text{soft}} = q_\mathrm{CDW}$, hence, the charge order is implemented indirectly through the dispersion relation for the softened phonon. Other effects of charge order will not be taken into account explicitly to avoid increasing the number of parameters in the model in the absence of suitable experimental input, although this extension may be considered in future work, e.g. along the lines of [25]. In the following sections, the computational methods will be described with more details.

3.1. Theoretical method and computational details

Let us denote with the symbol q the wave-vector and with the symbol x the doping level. We start by considering an optical phonon mode whose dispersion relation $\omega(q,x)$ is given by:

$$\begin{equation} \omega^2-\Omega\left(q,x\right)^2+\mathrm{i} \omega \Gamma\left(q,x\right) = 0\,.\end{equation} \tag{ 1 }$$

Here, $\Omega(q,x)$ is the renormalized energy of the phonon mode and $\Gamma(q,x)$ the corresponding linewidth. We emphasize that the energy $\Omega(q,x)$ contains already the corrections from the real part of the self-energy arising from phonon-phonon scattering, i.e. anharmonicity, and it is not the bare phonon energy. On the other hand, the linewidth $\Gamma(q,x)$ is related to the imaginary part of the self-energy. The spectral function for this mode is given by the standard Lorentzian form

$$\begin{equation} \mathcal{B}\left(q,\omega,x\right) = \frac{\omega\Gamma\left(q,x\right)}{\pi\left[\left(\omega^2-\Omega^2\left(q,x\right)\right)^2+\omega^2\Gamma^2\left(q,x\right)\right]}.\end{equation} \tag{ 2 }$$

By assuming a constant electronic density of states (DOS) at the Fermi momentum and wave-vector independent electron–phonon matrix elements, the Eliashberg function is given by (see the detailed derivation in [17])

$$\begin{equation} \alpha^2 F\left(\omega,x\right) = C\left(x\right)\int_0^4 \mathcal{B}\left(k_\mathrm{F} \sqrt{\zeta},\omega,x\right) \mathrm{d}\zeta,\end{equation} \tag{ 3 }$$

where $\zeta = 2 k_\mathrm{F}^{2}(1+\cos \theta)$, with θ the polar angle of the Fermi sphere, $k_\mathrm{F}$ the Fermi momentum, and C(x) is an effective parameter incorporating the microscopic information about the screening and the electron–phonon matrix elements. While the electron–phonon matrix in the cuprates is arguably anisotropic, here C(x) is intended to represent an angular average over the different directions [32]. Finally, the Fermi wave-vector is fixed to $k_\mathrm{F} = 0.5 q_{\text{max}}$ (see below for the definition of q_max) to make sure that the ζ variable spans the whole reciprocal space, $\zeta\in [0,4]$. We can then define the electron–phonon coupling parameter λ as customary in Eliashberg's theory [33]

$$\begin{equation} \lambda\left(x\right) = 2 \int_{0}^{\infty}\frac{\alpha^2 F\left(\omega,x\right)}{\omega}\mathrm{d}\omega\,.\end{equation} \tag{ 4 }$$

Finally, we will compute the critical temperature using the Allen–Dynes formula [30]:

$$\begin{equation} T_\mathrm{c}\left(x\right)\, = \,\frac{f_1\left(x\right)\,f_2\left(x\right)\,\omega_\mathrm{log}\left(x\right)}{1.2}\,\exp\left(-\frac{1.04\,\left(1+\lambda\left(x\right)\right)}{\lambda\left(x\right)-\mu^\star\,-\,0.62\,\lambda\left(x\right)\,\mu^\star}\right),\end{equation} \tag{ 5 }$$

where the logarithmic average frequency of the spectrum is defined as:

$$\begin{equation} \omega_\mathrm{log}\left(x\right) = \exp \left(\frac{2}{\lambda\left(x\right)}\int_{0}^{\infty} \mathrm{d}\omega \frac{\alpha^{2}F\left(\omega,x\right)}{\omega} \ln \omega \right).\end{equation} \tag{ 6 }$$

Following [30], we have also defined

$$\begin{align} & f_1\left(x\right) = \left[1+\left(\lambda\left(x\right)/\Lambda_1\right)^{3/2}\right]^{1/3},\nonumber\\ & f_2\left(x\right) = 1+\dfrac{\left(\bar{\omega}_2\left(x\right)/\omega_\mathrm{log}\left(x\right)-1\right)\lambda^2\left(x\right)}{\lambda^2\left(x\right)+\Lambda_2^2\left(x\right)},\nonumber\\ & \bar{\omega_2}\left(x\right) = \left[\dfrac{2}{\lambda\left(x\right)}\int_0^\infty \alpha^{2}F\left(\omega,x\right) \omega\, \mathrm{d}\omega \right]^{1/2},\nonumber\\ & \Lambda_1 = 2.46\left(1+3.8 \mu^\star\right),\nonumber\\ & \Lambda_2\left(x\right) = 1.82\left(1+6.3 \mu^\star\right)\left(\bar{\omega}_2\left(x\right)/\omega_\mathrm{log}\left(x\right)\right),\nonumber\\ & \mu^\star = 0.\end{align} \tag{ 7 }$$

Here, $\mu^\star = 0$ is chosen as typical value for strong-coupled superconductors [34].

At the heart of the model lies the effective description of the soft optical phonon, in this case the A_g in-plane oxygen mode which is experimentally known to exhibit a large softening upon increasing the doping x [13, 14]. We adopt the effective description of [17] where the softened optical phonon is given in terms of a flat dispersion with a Gaussian-shaped dip:

$$\begin{align} & \Gamma\left(q,x\right) = \gamma\left(x\right)\,,\end{align} \tag{ 8 }$$

$$\begin{align} &\Omega\left(q,x\right) = h \left(1-\delta\left(x\right)\, \mathrm{e}^{-\left({q/q_{\text{max}}}- q_{\text{soft}}/q_{\text{max}}\right)^2/\beta^2\left(x\right)}\right).\end{align} \tag{ 9 }$$

We notice that, within this description, the softening of the optical phonon is parameterized at a specific wave-vector q_soft. In other words, the experimental values reported in figure 1(b) for the phonon energy will be identified with $\Omega(q = q_{\text{soft}},x) = h(1-\delta(x))$. On the other hand, the corresponding linewidth, shown in figure 1(a), will be identified with the wave-vector independent quantity $\gamma(x)$. Finally, the constant h, which is doping level and wave-vector independent, corresponds to the frequency of the non-softened optical phonon measured experimentally, that is, $h = 135.6 \,\cdot 10^{11} \text s^{-1}$.

The physical parameter $\delta(x) \in [0,1]$ is the dimensionless depth of the softening, which can be obtained by fitting the experimental data as shown in figure 1(b). Meanwhile, $\beta(x)$ is a free parameter corresponding to the dimensionless width of the softening in units of the maximum wave-vector, for which we keep it a constant $\beta(x) = 0.2$ in our calculation. Moreover, q_soft is the softening wave-vector in units of the maximum wave-vector value q_max, which also represents a free parameter in the theoretical model. As explicitly demonstrated below, the choice of its value does not affect qualitatively the predictions of the model. Finally, q_max is used to set the characteristic length scale of the system and its value is therefore unimportant for the analysis below. Physically, it can be thought of as the size of the Brillouin zone.

The Raman measurements of the in-phase A_g mode at varying doping level x [13, 29] thus fix the values of $\mathrm{Re}[\omega(q_{\text{soft}})]$ upon varying x. This, in turn, fixes h and δ as a function x. Let us emphasize that since $\Gamma \ll \Omega$ (see figure 1), then we can safely assume $\mathrm{Re}[\omega(q_{\text{soft}})] = \Omega(q_{\text{soft}})$ to simplify our analysis. This leads immediately to mapping the Raman experimental data for the phonon energy onto the combination $h(1-\delta(x))$ in our mathematical model, where all the doping dependence goes into the parameter $\delta(x)$. We are thus left with a free parameter $\beta(x)$ corresponding to the width of the softening, which is the only non-trivial free parameter not constrained by the experimental data, along with the overall re-scaling factor C(x) in equation (3).

To reproduce the experimental $T_\mathrm{c}$, it is necessary that the value of C(x) depends on the values of doping x, as well as $\beta(x)$, and that it displays a physically-reasonable trend compatible with typical momentum-resolved measurements of [19]. In particular, it is expected that the width of the softening increases when the softening depth decreases [19, 35]. Let us also emphasize that we assume for C(x) a linear dependence from x as the simplest functional dependence that allows one to reproduce, in an approximate way, the experimental $T_\mathrm{c}$. This procedure gives the linear dependence shown in figure 2. Finally, the value of the softened optical phonon linewidth $\Gamma(q,x)$ is set equal to $\gamma(x)$, independent of q as per Klemens' approximation for optical phonon damping [36], and its value is fixed by the experimental Raman measurements [13, 29] (see panel (b) in figure 1).

3.2. Results

Since the Raman measurement is performed near the Γ-point but does not specify the precise value of q, we consider two different scenarios for the location of the softening. In the first scenario, we assume the softening to be located in the middle of the Brillouin zone, while in the second scenario we consider the softening right at the Γ-point, i.e. $q_{\text{soft}} = 0$.

3.2.1. Softening at finite q

We thus proceed to an empirical fitting (interpolation) of the softened A_g mode data of frequency and linewidth from experimental data [13, 29], as shown in figure 1, as a function of doping level x. This will completely determine the form of $\gamma(x)$ and $\delta(x)$ entering the theoretical model. We are then left only with C(x) (a linear function of x), q_soft and $\beta(x)$ (allowed to vary with x) as fitting parameters, to match the experimentally measured $T_\mathrm{c}$ vs x curve. The linear dependence of C(x) from x is fixed also by imposing that the predicted and experimental $T_\mathrm{c}$ values agree for the lowest doping x ≈ 6.4. By fixing $q_\textrm{soft}/q_\textrm{max} = 0.5$ (half of the maximum wave-vector), we obtain the dispersion relation of the phonon mode for different doping levels, which is shown in figure 2(a) and the corresponding Eliashberg function which is obtained using the phenomenological parameter $\beta(x) = 0.2$. We have explicitly verified that a variation in $q_{\text{soft}}/q_\textrm{max}$ does not strongly affect the final results (see below). All in all, this whole procedure leads to the predicted $T_\mathrm{c}$ vs doping x trend shown in figure 2(c). The corresponding electron phonon coupling $\lambda(x)$ as a function of the doping level x is shown in figure 2(d).

Based on the results shown in figure 2, it can be seen that both the electron–phonon coupling and $T_\mathrm{c}$ are dominated by the effective electron–phonon matrix element C(x), while the effect of other parameters like q_soft, $\beta(x)$, $\delta(x)$ are weaker. We then attempt to implement C(x) to be softening-dependent, i.e. $C(x) \propto \delta(x)$. The results are shown in figure 3, where it can be seen that the plateau of $T_\mathrm{c}$ in the range $x = 6.6-6.8$ can be qualitatively reproduced.

Zoom In Zoom Out Reset image size

The above analysis provides an interesting link between in-phase oxygen optical mode softening due to doping and the high-temperature $T_\mathrm{c}$ in the slightly overdoped/nearly optimally doped cuprates. In particular, it suggests that the (much debated) flattening of $T_\mathrm{c}$ upon increasing doping in the regime $x = 6.5-6.8$ is directly linked to the slower decrease of phonon energy with x in that same doping regime: cfr. the horizontal oval shaded regions in figures 1(b) and 3(b). Less trivial, and even more interesting, is the explanation of the further peak in $T_\mathrm{c}$ vs x in the region $x = 6.8-7.0$. This effect is directly linked to the sharpening of the softening dip as reflected by the decrease of the width parameter $\beta(x)$, as is evident by comparing the vertical shaded ovals in figure 3(b).

3.2.2. Softening at the Γ-point, $q_{\textrm{soft}} = 0$

As emphasized above, the theoretical model contains, by construction, three unknown physical parameters $C(x),\beta(x), q_{\text{soft}}$ which cannot be determined from the experimental data. The first two have been discussed in detail above. The third parameter, which describes the position of the softening within the Brillouin zone, needs further clarifications. In particular, in our first attempt, we have imposed the phonon softening to occur somewhere between the Γ-point and the Brillouin zone boundary. Nevertheless, the experimental data from [13, 29] reported in figure 1 are measured at low wave-vector. In order to prove that the choice of q_soft does not affect our results qualitatively, we have performed the same computation by assuming that the softening region appears right at the Γ point in the Brillouin zone, coinciding with $q_{\text{soft}} = 0$. The corresponding dispersion relation for the softened phonon is shown in figure 4 together with the theoretical predictions for $T_\mathrm{c}$ (compared to the experimental data) obtained with the two schemes of C(x) linear in x as in figure 2 and $C(x)\propto \delta(x)$ as in figure 3. Evidently, the results remain qualitatively the same, proving the subleading effect of the choice of q_soft, as already advertised above.

Zoom In Zoom Out Reset image size

3.3. Physical justification for the x-dependence of C(x)

We have seen that the experimental critical temperature can be reproduced, as a first approximation, with C(x) as a linear function of x. A reasonable physical justification which supports that C(x) grows linearly with the oxygen doping x goes as follows. In [17] the effective parameter C is defined as directly proportional to the electron–phonon matrix element $g_{\boldsymbol{k}\boldsymbol{k^{^{\prime}}}}$, that is, $C \sim g_{\boldsymbol{k}\boldsymbol{k^{^{\prime}}}}$. The latter, according to Bardeen–Pines theory accounting for Thomas–Fermi screening, reads as:

$$\begin{equation} g_{\boldsymbol{k}\boldsymbol{k^{^{\prime}}}}^2 = |v_{q}^{Z}/\epsilon\left(0,q\right)|^2 = \left|\frac{4\pi Ze^2}{q}\left(\frac{N}{M}\right)^{1/2}\frac{q^2}{q^2+k_{s}^{2}}\right|^2,\end{equation} \tag{ 10 }$$

where $\epsilon(0,q) = 1+\frac{k_{s}^{2}}{q^{2}}$ is the static dielectric constant, computed by solving the electrostatic Poisson-equation problem for an electron gas in presence of the positive background charge of the lattice [38], and k_s is the screening wave-vector. For a quantum Fermi gas at T = 0, the screening wave-vector is given by the Thomas–Fermi wave-vector as [39]: $k_\textrm{TF}^{2} = 4 \left(\frac{3 n}{\pi}\right)^{1/3}$, which is the squared inverse Thomas–Fermi screening length, with n the number density of electrons. At higher temperatures, for a non-degenerate electron gas (the Debye–Hückel screening), $k_{s} = \sqrt{\frac{4 \pi e^2 n}{k_\textrm{B}T}}$, with e the electron charge and $k_\textrm{B}$ the Boltzmann constant. For high-temperature SCs such as YBCO we can assume that we are in the non-degenerate Debye–Hückel regime. Hence, upon inserting the Debye–Hückel screening k_s into equation (10), we obtain, in the regime $k_s^2 \gg q^2$ (low phonon wavevectors),

$$\begin{equation} C \sim g_{\boldsymbol{k}\boldsymbol{k^{^{\prime}}}} \sim 1/n.\end{equation} \tag{ 11 }$$

At chemical equilibrium, the concentration of electrons, n, and the concentration of holes, p, are related via [38]:

$$\begin{equation} n\,p = \text{const}\left(T\right).\end{equation} \tag{ 12 }$$

Hence, at a fixed temperature, $n \sim 1/p$. Since in YBCO the concentration of holes increases proportionally to the oxygen doping x [40], upon inserting back into equation (11), we finally obtain

$$\begin{equation} C\left(x\right) \propto x,\end{equation} \tag{ 13 }$$

which thus justifies the ansatz used above in our model.

In the above approach we neglected the well-document d-wave character of Cooper pairs in the cuprates [41]. In what follows, we implement the phonon softening model and the doping-dependent effective electron–phonon coupling C(x) into the d-wave formulation of Eliashberg theory. We wish to reproduce, in a approximate way, the experimental critical temperatures versus oxygen doping with a minimal model so we neglect a lot of previous parameters as q_soft, $\beta(x)$, $\delta(x)$ and we use just the experimental data: the phonon representative energy $\Omega_{0}(x)$ and the Raman linewith $\Gamma(x)$. We do this only to simplify the calculation as much as possible and to minimize the number of free parameters. In the end we will be able to reproduce the trend of T_c as a function of x quite well using a single free parameter for which we will give a physical justification. The first step is to determine the input parameters of the theory as a function of the oxygen doping x. We fit the experimental data presented in figure 1 with a third-degree polynomial for the Raman frequency, $\Omega_{0}(x) = 8496.79564 - 3799.76085x + 570.27049x^{2} - 28.537x^{3}$, and for the linewidth, $\Gamma(x) = 6894.76021 - 3072.10979x + 455.86155x^{2} - 22.5184x^{3}$. The fit is equally good and indistinguishable from the fifth-order polynomial fit shown in figure 1. In this way, we can have the input parameters for calculating, in an analytical way, the spectral function $\alpha^{2}F_{s}(\Omega,x)$.

In what follows, we denote the peak phonon energy with Ω₀, the generic phonon energy with Ω, and the electron energy with ω.

The electron–phonon spectral function now is

$$\begin{equation} \begin{aligned} \alpha^{2}F_{s}\left(\Omega,x\right) = \frac{C\left(x\right)\Omega\Gamma\left(x\right)}{\pi\left[\left(\Omega^{2}-\Omega_{0}\left(x\right)^{2}\right)^{2}+\Omega^{2}\Gamma^{2}\left(x\right)\right]}\\=\frac{\lambda_{s}\left(x\right)\Omega\Omega_{0}\left(x\right)\Gamma\left(x\right)}{\pi\left[\left(\Omega^{2}-\Omega_{0}\left(x\right)^{2}\right)^{2}+\Omega^{2}\Gamma^{2}\left(x\right)\right]} \end{aligned}\end{equation} \tag{ 14 }$$

where $C(x) = \lambda_{s}(x)\Omega_{0}(x)$ and $\lambda_{s}(x)$ is the value of the electron–phonon coupling constant. We calculated the critical temperatures by solving, in a numerical way, the one band d-wave Eliashberg equations [34, 42–48]. In this case two coupled equations for the gap $\Delta(i\omega_{n},\phi)$ and renormalization functions $Z(i\omega_{n},\phi)$ have to be solved (ω_n denotes the Matsubara frequencies and φ and $\phi^{^{\prime}}$ are the azimuthal angles of k and k' in the ab plane). The d-wave one-band Eliashberg equations in the imaginary axis representation, when the Migdal theorem is valid [49], are:

$$\begin{align} \omega_{n}Z\left(\omega_{n},\phi\right)& = \omega_{n}+ \pi T \sum_{m}\int_{0}^{2\pi}\frac{\textrm{d}\phi^{^{\prime}}}{2\pi}\Lambda\left(\omega_{n},\omega_{m},\phi,\phi^{^{\prime}}\right)N_{Z}\left(\omega_{m},\phi^{^{\prime}}\right)\end{align} \tag{ 15 }$$

$$\begin{equation} \begin{aligned} Z\left(\omega_{n},\phi\right)\Delta\left(\omega_{n},\phi\right) = \pi T \sum_{m}\int_{0}^{2\pi}\frac{\textrm{d}\phi^{^{\prime}}}{2\pi}\left[\Lambda\left(\omega_{n},\omega_{m},\phi,\phi^{^{\prime}}\right)- \mu^{*}\left(\phi,\phi^{^{\prime}}\right)\Theta\left(\omega_{c}-|\omega_{m}|\right)\right]N_{\Delta}\left(\omega_{m},\phi^{^{\prime}}\right) \end{aligned}\end{equation} \tag{ 16 }$$

where $\Theta(\omega_{c}-|\omega_{m}|)$ is the Heaviside function, ω_c is a cut-off energy and

$$\begin{equation} \begin{aligned} &\Lambda\left(\omega_{n},\omega_{m},\phi,\phi^{^{\prime}}\right)\\ & = 2\int_{0}^{+\infty}\Omega \textrm{d}\Omega \alpha^{2}F\left(\Omega,\phi,\phi^{^{\prime}}\right)\Big/\left[\left(\omega_{n}-\omega_{m}\right)^{2}+\Omega^{2}\right], \end{aligned}\end{equation} \tag{ 17 }$$

$$\begin{equation} N_{Z}\left(\omega_{m},\phi\right) = \frac{\omega_{m}}{\sqrt{\omega^{2}_{m}+\Delta\left(\omega_{m},\phi\right)^{2}}},\end{equation} \tag{ 18 }$$

$$\begin{equation} N_{\Delta}\left(\omega_{m},\phi\right) = \frac{\Delta\left(\omega_{m},\phi\right)}{\sqrt{\omega^{2}_{m}+\Delta\left(\omega_{m},\phi\right)^{2}}}.\end{equation} \tag{ 19 }$$

We assume [34, 42–46] that the electron–phonon spectral function $\alpha^{2}(\Omega)F(\Omega,\phi,\phi^{^{\prime}})$ and the Coulomb pseudopotential $\mu^{*}(\phi,\phi^{^{\prime}})$ at the lowest order contain separated s- and d-wave contributions $\alpha^{2}(\Omega)F_{s,d}(\Omega,x)$ and $\mu^{*}_{s,d}(\omega_{c})$,

$$\begin{equation} \begin{aligned} \alpha^{2}F\left(\Omega,\phi,\phi^{^{\prime}}\right) = \lambda_{s}\alpha^{2}F_{s}\left(\Omega\right) +\lambda_{d}\alpha^{2}F_{d}\left(\Omega\right)\sqrt{2}\cos\left(2\phi\right)\sqrt{2}\cos\left(2\phi^{^{\prime}}\right), \end{aligned}\end{equation} \tag{ 20 }$$

$$\begin{equation} \mu^{*}\left(\phi,\phi^{^{\prime}}\right) = \mu^{*}_{s}\left(\omega_{c}\right) +\mu^{*}_{d}\left(\omega_{c}\right)\sqrt{2}\cos\left(2\phi\right)\sqrt{2}\cos\left(2\phi^{^{\prime}}\right)\end{equation} \tag{ 21 }$$

as well as the self energy functions:

$$\begin{equation} Z\left(\omega_{n},\phi\right) = Z_{s}\left(\omega_{n}\right)+Z_{d}\left(\omega_{n}\right)\cos\left(2\phi\right),\end{equation} \tag{ 22 }$$

$$\begin{equation} \Delta\left(\omega_{n},\phi\right) = \Delta_{s}\left(\omega_{n}\right)+\Delta_{d}\left(\omega_{n}\right)\cos\left(2\phi\right).\end{equation} \tag{ 23 }$$

The spectral functions $\alpha^{2}F_{s,d}(\Omega)$ are normalized in a way that $2\int_{0}^{+\infty}\frac{\alpha^{2}F_{s,d}(\Omega)}{\Omega}\textrm{d}\Omega = \lambda_{s,d}$ and, we choose, for simplicity, $\lambda_{s} = \lambda_{d}$. Also, we assume that $\alpha^{2}F_{s}(\Omega) = \alpha^{2}F_{d}(\Omega)$ i.e. also the shape of the spectral functions are the same. In this case we search for solutions of the Eliashberg equations in a pure d-wave form, as indicated by the experimental data, for the gap function $\Delta(\omega,\phi^{^{\prime}}) = \Delta_{d}(\omega)\cos(2\phi)$ (the s component is zero and this happens for example when [34] $\mu^{*}_{s} \gg \mu^{*}_{d}$). In the more general case $\Delta(\omega_{n},\phi)$ has d and s components. The renormalization function $Z(\omega,\phi^{^{\prime}}) = Z_{s}(\omega)$ has just the s component because the equation for $Z_{d}(\omega)$ is a homogeneous integral equation whose only solution in the weak-coupling regime is $Z_{d}(\omega) = 0$ [50]. The cut-off energy is $\omega_{c} = 450$ meV and the maximum quasiparticle energy is $\omega_\textrm{max} = 500$ meV. In first approximation, we put $\mu^{*}_{d} = 0$ (if the s component of the gap is zero, the value of $\mu^{*}_{s}$ is irrelevant). If the s component of the gap function is absent, this means that $\mu^{*}_{s} \gg \mu^{*}_{d}\simeq0$ [34] In this model, the experimental input parameters are $\alpha^{2}F_{s,d}(\Omega)$ and $\mu^{*}_{s,d}$ but, with our previous assumptions, we have just one free parameter: $\lambda_{s}(x)$. To start with, we solve the imaginary axis d-wave Eliashberg equations and we fix the value of $\lambda_{s}(x)$ in such a way to obtain the exact experimental critical temperature $T_\mathrm{c}$. The corresponding obtained values of $\lambda_{s}(x)$ are shown in figure 5(a). We can see that the value of the coupling constant increases by a factor three over the examined oxygen doping x range, but remains in a moderate strong-coupling regime, i.e. less than lead (Pb), for example ($\lambda_\textrm{Pb} = 1.55$). The fact that the coupling constant increases by a factor three is connected with the modest variation of the representative energy of the system, which actually decreases as x increases. As we know from theory, T_c is an increasing function of the representative energy of the phonon spectrum (in our case, roughly, Ω₀, the value of the peak of the phonon Lorentzian spectrum) and of the electron–phonon coupling constant.

Zoom In Zoom Out Reset image size

At this point we recall that $\lambda_{s}(x) \propto \frac{1}{N(E_\textrm{F},x)}$ and we use the first value of $\lambda_{s}(x = 6.41) = 0.463$ that produces the exact experimental $T_{c}(x = 6.41) = 21.2$ K for calculating the other values according to the following law: $\lambda_{s}(x) = \frac{N(E_\textrm{F},x \;= \; 6.41)}{N(E_\textrm{F},x)}\lambda_{s}(x = 6.41)\simeq [1+\frac{1}{N(E_\textrm{F},x \;=\; 6.41)}\frac{\partial N(E_\textrm{F},x \;=\; 6.41)}{\partial x}(x-6.41)]\lambda_{s}(x \;= \; 6.41) = [1+k(x-6.41)]\lambda_{s}(x = 6.41)$ where $k = \frac{1}{N(E_\textrm{F},x \; = \;6.41)}\frac{\partial N(E_\textrm{F},x \; = \;6.41)}{\partial x}$. The fit is shown in figure 5(a). We find $k = -1/4$. Now we can calculate the critical temperature without free parameters. The final results are shown in figure 5(b). Our model is too simple to explain the plateau of T_c for x close to 7, but in its extreme simplicity it manages to account for the trend of T_c in almost the entire oxygen doping range. Studies on the dependence of the DOS at the Fermi level, as a function of oxygen doping, would be needed to verify our hypothesis in greater detail.

A qualitative argument to explain the decreasing trend of $N(E_\textrm{F})$ vs x goes as follows. Assuming an electronic DOS for the free electrons, $N(E) \sim E^{1/2}$, at the Fermi level we have $N(E_\textrm{F}) \sim E_\textrm{F}^{1/2}$. For the free electrons, $\sqrt{E_\textrm{F}} \sim n^{1/3}$. At chemical equilibrium between electrons and holes, $n\,p = const$, hence $n \sim 1/p$. Hence, upon further taking $p \sim x$ [40], we get $N(E_\textrm{F}) \sim x^{-1/3}$. Hence, even within the simplest assumptions, it is physically meaningful that the DOS at Fermi level decreases with increasing the level of doping x.

In conclusion, we have developed a simple theoretical framework to rationalize the correlation between optical phonon softening and superconducting $T_\mathrm{c}$ in YBa₂Cu₃O_x, in excellent qualitative agreement with experimental data. The theoretical model is based on the Migdal–Eliashberg theory of strong-coupling superconductivity with the pairing being mediated by the A_g oxygen mode that undergoes strong softening upon increasing the doping x as observed experimentally [13] and couples strongly to the nodal electronic states [51].

We presented two versions of this model. In the first version, we used a simplified Eliashberg-type argument involving three free parameters. One is the effective electron–phonon matrix element C, the other is the width $\beta(x)$ of the softened optical phonon of the in-phase A_g oxygen Raman-active mode, and the last one is the location of the softening region within the Brillouin zone, q_soft. We emphasize that, out of these three parameters, only two, the width of the softening region in the phonon dispersion, and the parameter C(x), depend on the doping level x. In particular, we found that the experimental data can be reproduced if C(x) grows linearly with x. A physical justification is provided in terms of the Bardeen-Pines electron–phonon matrix element, which is approximately inversely proportional to the electron concentration n. Assuming chemical equilibrium between electron and holes, and since the hole concentration p is proportional to the oxygen doping x, we showed that indeed C(x) increases linearly with x. The analysis using this model reveals a connection between the doping and charge order-induced softening of the in-phase oxygen vibrational mode and the superconducting $T_\mathrm{c}$, which is otherwise difficult, if not impossible, to explain with alternative theoretical approaches (e.g. spin fluctuations [52]). The proposed model provides a link between the two phenomena and suggests a moderate role of the phonon linewidth $\gamma(x)$ in such a regime of doping. Additionally, the model predicts that the additional plateau feature observed at optimally doping, x > 6.9, is a direct manifestation of the doping dependence of the softening width $\beta(x)$, which cannot be neglected anymore for optimally doping levels. In other words, the extra plateau in $T_\mathrm{c}$ vs doping x past the flattening region, in the optimally doped range x > 6.9, may be ascribed to the strong sharpening of the softening dip of the oxygen phonon mode in that range. An otherwise constant-with-doping softening width would lead to an unphysical monotonically growing $T_\mathrm{c}$ as a function of x, which is not observed experimentally. Interestingly, the theoretical model returns, as one of its output, the momentum-resolved dispersion relation of the in-phase oxygen mode with a width of the softening that decreases with the depth of the softening, in qualitative agreement with typical observations of optical phonon softening in superconductors.

The second version of the model is based on the same input in terms of A_g phonon mode softening, from the experimental data, which is implemented in the d-wave Eliashberg theory formalism. This is a more appropriate and exact method, which does not rely on the Allen–Dynes formula (which works, in an approximate way, just for intermediate coupling, while it fails for the strong coupling regime). This second version of the model by calibrating the doping-dependence of the electron–phonon coupling λ on a single value of oxygen doping, is able to reproduce the $T_\mathrm{c}$ values as a function of oxygen doping level x in excellent agreement with the experimental data. The calibration implies that the DOS of electrons at the Fermi level is a monotonic decreasing function of x, which can be rationalized, once again, using a simple argument which assumes a Fermi DOS and the chemical equilibrium between electrons and holes.

The present model provides a simple solution to a few puzzles in linking doping-induced anomalous softening of oxygen modes and superconducting critical temperature in the cuprates and may serve as a useful tool for engineering and discovering new high-temperature materials in future work in conjunction with atomistic approaches [53]. Our findings also emphasize that the role of lattice distortions in high-$T_\mathrm{c}$ superconductors might be more important than currently thought and should probably be carefully revisited.

Finally, the role of charge order has been simplified in our model. The important role of charge order leading to anomalous softening and suppression of $T_\mathrm{c}$ has been shown by several studies [20, 22, 54]. Hence, further work on the interplay between phonon, charge order, and superconductivity, is needed.

We would like to thank Chandan Setty and Antonio Sanna for collaboration on related topics and useful comments on a preliminary version of this draft. C J and M B acknowledge the support of the Shanghai Municipal Science and Technology Major Project (Grant No. 2019SHZDZX01) and the sponsorship from the Yangyang Development Fund. A Z gratefully acknowledges funding from the European Union through Horizon Europe ERC Consolidator Grant number: 101043968 'Multimech', and from US Army Research Office through contract nr.W911NF-22-2-0256.

All data that support the findings of this study are included within the article (and any supplementary files).

In this appendix we show, in figure 6 below, the Raman scattering data of the in-phase oxygen mode measured at liquid nitrogen temperatures (T = 77 K) in YBa₂Cu₃O_x. As explained in the main article, these data have intrinsically lower statistics compared to the room temperature data.

Zoom In Zoom Out Reset image size

Nevertheless, these data clearly show the same characteristics of in-phase oxygen optical phonon softening that are retained in the Raman scattering measurements at room temperature (figure 1). The latter were used as input to our theoretical model in the main article.