Introduction

The discovery of charge–density wave (CDW) order1 in underdoped cuprates raised the question of whether it is intimately related to pseudogap physics2,3,4,5,6,7, and thereby yet another signature of strong electronic correlations. Charge modulations with a moderate correlation length are detected by nuclear magnetic resonance8,9,10, scanning tunneling microscopy (STM)11,12,13, and X-ray techniques11,12,14,15,16, with the latter finding incommensurate wave-vectors near qCO = 0.3 reciprocal lattice units, oriented along the crystalline axes with either uniaxial or biaxial character. The CDW wavevector qCO continuously drops with increasing hole doping in YBa2Cu3O6+x (YBCO) and Bi2Sr2−xLaxCuO6+x (Bi-2201)1,11,16,17, whereas the completely opposite trend is observed in charge-stripe ordered La-based 214 cuprates such as La2−xSrxCuO4 (LSCO) and La2−xBaxCuO4 (LBCO)18. Also the concomitant spin-stripe order marks 214 cuprates to be different from YBCO or Bi-2201. Due to these apparent differences we focus on the charge ordering in YBCO.

Stripe formation and competition with d-wave superconductivity or pair-density waves was investigated for multiband Hubbard or t–J models with accurate computational tools19,20,21. Charge-stripe ground states were found with a near-degeneracy of competing states. So far, the axial orientation and the incommensurate character of qCO, as well as its doping variation in YBCO and Bi-2201, proved difficult to reconcile with purely electronic model calculations20,22. Weak-coupling theories typically predict that qCO lies on the Brillouin zone (BZ) diagonal5,23,24 unless the CDW instability is preceded by a Fermi surface reconstruction to form hole pockets7. Based on this evidence, we conjecture that some missing ingredient may still be required to explain charge-order in YBCO and Bi-2201.

In the continuing search for this ingredient, we employ further experimental facts. First, the discovery of charge order in overdoped Bi-2201, with an unreconstructed Fermi surface25, poses a question as to whether the CDW is actually tied to the pseudogap phenomena. Second, the softening or broadening of phonon modes at the CDW wavevector26,27 along with giant phonon anomalies near the CDW instability, all point to a strong electron–phonon (el–ph) coupling28. While the phonon dispersion will necessarly react to a charge–density modulation, phonons may also play a key role in CDW formation, for example selecting the ordering wavevector in CDW-susceptible materials29,30. In cuprates, the frequencies of the phonons decrease only weakly, and a continuous softening to zero frequency is likely precluded by quenched disorder10,31.

A third hint comes from the atomic displacement pattern that accompanies the CDW. A complete set of X-ray diffraction data, by Forgan et al.32, found the by-far largest displacements for planar oxygen atoms; they shift out-of-plane with an out-of-phase pattern (see Fig. 1) that closely resembles the normal mode of the B1g bond-buckling phonon33. In stoichiometric, overdoped YBCO the dispersion of the B1g phonon was monitored upon cooling26; even though the stoichiometric composition of YBCO does not support CDW order, the B1g phonon frequency softens by about 6% at the charge ordering wavevector for underdoped YBCO. The X-ray diffraction data in Forgan et al.32 reveal that also the heavier ions in YBCO slightly move away from the positions they take in the charge homogeneous phase in equilibrium. This naturally suggests that also the low-energy phonons are involved in the CDW formation, and indeed the softening of a low-energy optical phonon was observed by Le Tacon et al.28 and Kim et al.27. While the softening of this mode indeed indicates a coupling to the CDW (as must occur in the absence of special symmetries), the small static displacements of the Ba and Y atoms observed in X-ray diffraction suggests ultimately that there is a relatively weak participation of this mode in the CDW formation.

Fig. 1: Oxygen vibrations with B1g symmetry.
figure 1

A representation of the out-of-plane oxygen vibrations (B1g pattern) in a copper–oxide bilayer. The arrows represent the out-of-phase motions of the oxygen atoms.

Full size image

There is, furthermore, a bond-stretching phonon that develops a Kohn-like anomaly as temperature is reduced34,35,36,37,38,39,40,41. This temperature-dependent softening was initially understood to indicate strong el–ph coupling, and indeed the size of the anomaly correlates with the superconducting Tc40. However, trends with doping do not show any clear correlation with spin-charge stripes39, static charge density waves42, or kinks in the electronic dispersion40, which seems to rule out conventional el–ph coupling mechanisms of the type contempleted here. Rather, it is proposed that the softening reflects coupling to a collective charge mode of the electron liquid40,43. In Forgan et al.32, oxygen displacements do indeed have in-plane components, indicating that the bond-stretching modes contribute alongside the buckling mode, but the out-of-plane components are dominant, and we make the simplifying approximation that the in-plane component can be neglected.

Based on these empirical grounds, we select the B1g bond-buckling phonon, reanalyze its coupling to the electrons in the copper–oxygen planes, and obtain the Landau free energy for a B1g-phonon mediated CDW. We start from a microscopic model for the CuO2 planes, and find that the structure of the el–ph coupling matrix element g(q; k) depends strongly on the orbital content of the Fermi surface. In particular, the Cu-4s orbital is crucial; with it, g(q; k) is maximum at phonon momenta q* that are axial (rather than diagonal) and track the doping dependence of the CDW wavevector. The special role of q* is imperceptible in the free energy when calculated with bare electron dispersions.

Here, we show that with the inclusion of electronic correlations, modeled by a renormalization of the band dispersion and quasiparticle spectral weight, an axial wavevector close to q* emerges as the dominant wavevector in the CDW susceptibility. While the B1g mode is, by itself, too weak (i.e., softens by only a few percent) to induce true long-range order, the inevitable presence of disorder will necessarily slow down and pin the CDW fluctuations, creating local patches of static or quasi-static charge order with finite correlation lengths. This leads us to propose a scenario, in which a phonon-based mechanism is enabled by strong electronic correlations, and the incommensurate wavevector of the concomitant charge correlations is dictated by the momentum–space structure of the el–ph coupling matrix element.

Results

Three-orbital electronic model

To set the stage for the essential features of the electronic structure we select YBCO as the target material. The unit cell of YBCO contains a bilayer of CuO2 planes. The coupling between these two layers split the individual-layer derived bands. With respect to CDW formation it was demonstrated in earlier theoretical work44 that the bilayer splitting determines primarily the relative orientation and phase of the CDWs in the two separate layers but does not have much effect on the structure of the CDW in the individual layers themselves. In YBCO, the bilayer splitting even collapses in the near-nodal region at low doping45. We therefore ascribe no vital role to the bilayer structure for the CDW itself and henceforth focus on a single copper–oxygen plane.

We start by modeling a single CuO2 plane in YBCO including copper 4s and \(3{d}_{{x}^{2}-{y}^{2}}\) as well as oxygen px and py orbitals. An effective three-band model is obtained in terms of the d- and p-orbitals by downfolding the original four-band model to refs. 7,46

$${H}_{{\rm{kin}}}=\sum _{{\bf{k}},\sigma }{\Psi }_{{\bf{k}},\sigma }^{\dagger }\left(\begin{array}{rcl}{\varepsilon }_{d}&2{t}_{pd}{s}_{x}&-2{t}_{pd}{s}_{y}\\ 2{t}_{pd}{s}_{x}&{\tilde{\varepsilon }}_{x}({\bf{k}})&4{\tilde{t}}_{pp}{s}_{x}{s}_{y}\\ -2{t}_{pd}{s}_{y}&4{\tilde{t}}_{pp}{s}_{x}{s}_{y}&{\tilde{\varepsilon }}_{y}({\bf{k}})\end{array}\right){\Psi }_{{\bf{k}},\sigma },$$
(1)

with the three-spinor \({\Psi }_{{\bf{k}},\sigma }^{\dagger }=\left({d}_{{\bf{k}},\sigma }^{\dagger },\ {p}_{x{\bf{k}},\sigma }^{\dagger },\ {p}_{y{\bf{k}},\sigma }^{\dagger }\right)\) and \({s}_{x,y}=\sin ({k}_{x,y}/2)\). εd denotes the onsite energy of the d orbital, tpd and tps are the hopping amplitudes between p- and d- and p- and s-orbitals, respectively. In the downfolding procedure, the hopping processes via the copper 4s orbital renormalize the oxygen energies εp and generate indirect hopping \(4{t}_{pp}^{i}\) between oxygen orbitals:

$${\tilde{\varepsilon }}_{x,y}={\varepsilon }_{p}+4{t}_{pp}^{i}{s}_{x,y}^{2};\ {\tilde{t}}_{pp}={t}_{pp}^{i}+{t}_{pp}^{d};\ {t}_{pp}^{i}=\frac{{t}_{ps}^{2}}{{\varepsilon }_{{\rm{F}}}-{\varepsilon }_{s}},$$
(2)

where εF is the Fermi energy and \({t}_{pp}^{d}\) a small direct hopping amplitude. We adopt all the parameters entering Eqs. (1) and (2) from Andersen et al.46, specifically tpd = 1.6 eV, εd − εp = 0.9 eV, \({t}_{pp}^{d}=0\) and \({t}_{pp}^{i}=-1.0\) eV. We diagonalize Hkin and focus only on the partially filled anti-bonding band; the irrelevant spin index is subsequently suppressed.

Bond-buckling phonon coupling

The out-of-plane B1g vibrations of the oxygen atoms, as depicted in Fig. 1, naturally couple linearly to their local electric field Ez47. At this point, we neglect the motion of the almost four times heavier copper atoms. Consequently, we start from the ansatz

$${{\mathcal{H}}}_{{\rm{ep}}}=e{E}_{z}\sum _{{\bf{n}}}\left[\delta {u}_{x{\bf{n}}}{p}_{x{\bf{n}}}^{\dagger }{p}_{x{\bf{n}}}+\delta {u}_{y{\bf{n}}}{p}_{y{\bf{n}}}^{\dagger }{p}_{y{\bf{n}}}\right],$$
(3)

where δux,yn are the out-of-plane displacements of the oxygen atoms in unit cell n. We project the el–ph Hamiltonian in Eq. (3), onto the anti-bonding band and obtain the effective Hamiltonian \({\mathcal{H}}={{\mathcal{H}}}_{{\rm{el-ph}}}+{{\mathcal{H}}}_{{\rm{ph}}}\) with

$${{\mathcal{H}}}_{{\rm{el-ph}}}=\sum _{{\bf{k}}}{\varepsilon }_{{\bf{k}}}{c}_{{\bf{k}}}^{\dagger }{c}_{{\bf{k}}}+\sum_{{\bf{q}},{\bf{k}}}{\tt{g}}({\bf{q}};{\bf{k}}){c}_{{\bf{k}}+{\bf{q}}}^{\dagger }{c}_{{\bf{k}}}\left({a}_{{\bf{q}}}+{a}_{-{\bf{q}}}^{\dagger }\right),$$
(4)

where \({c}_{{\bf{k}}}^{\dagger }\) is the creation operator for the anti-bonding electrons with dispersion εk, aq annihilates a B1g phonon mode, and \({{\mathcal{H}}}_{{\rm{ph}}}=\hslash {\Omega }_{P}{\sum }_{{\bf{q}}}{a}_{{\bf{q}}}^{\dagger }{a}_{{\bf{q}}}\) (see Devereaux et al.47 and Supplementary Notes 1 and 2). The momentum-dependent el–ph coupling is written \({\tt{g}}({\bf{q}};{\bf{k}})=\gamma \tilde{{\tt{g}}}({\bf{q}};{\bf{k}})\)48, where γ is the coupling strength and

$$\tilde{{\tt{g}}}({\bf{q}};{\bf{k}})=\left[{e}_{{\bf{q}}}^{x}{\phi }_{x}({\bf{k}}^{\prime} ){\phi }_{x}({\bf{k}})+{e}_{{\bf{q}}}^{y}{\phi }_{y}({\bf{k}}^{\prime} ){\phi }_{y}({\bf{k}})\right],$$
(5)

with \({\bf{k}}^{\prime} ={\bf{k}}+{\bf{q}}\). The eigenfunctions ϕx,y signify the orbital contents of the oxygen px,y orbitals in the anti-bonding band. Similarly, the eigenvectors \({e}_{{\bf{q}}}^{x,y}\) correspond to the normal mode of the out-of-plane displacements of the two oxygen atoms in the CuO2 unit cell. The overall strength of the el–ph coupling is \(\gamma =e{E}_{z}\sqrt{\hslash /2m{\Omega }_{{\rm{P}}}}\), where m is the mass of the oxygen atom and ΩP ~40 meV is the frequency of the dispersionless B1g mode. Adopting the electric-field value eEz = 3.56 eV/Å from Johnston et al.49 leads to γ = 0.22 eV. The eigenvectors for the B1g mode are \({e}_{{\bf{q}}}^{x,y}=\mp \frac{\cos ({q}_{y,x}/2)}{{M}_{{\bf{q}}}}\) (see Supplementary Fig. 1), where the normalization factor \({M}_{{\bf{q}}}=\sqrt{{\cos }^{2}({q}_{x}/2)+{\cos }^{2}({q}_{y}/2)}\).

Earlier theoretical work33,50 on the B1g phonon in optimally doped Bi-2212, argued that for an antinodal fermion state, ∣g(q; k)∣2 is strongest for an axial scattering wavevector q to the nearby antinodal final state. This value of q is considerably smaller than qCO. Instead we find, when the Cu-4s orbital is properly included via the finite indirect hopping \({t}_{pp}^{i}\) in Eq. (2), the anisotropic structure of g(q; k) changes qualitatively. The maximum of the coupling ∣g(q; k)∣2 now occurs for larger axial wavevectors q*, that connect initial (\({{\bf{k}}}_{{i}}^{* }\)) and final (\({{\bf{k}}}_{{i}}^{* }\) + q*) Fermi surface states near the nodal points (see Fig. 2a, b). The resulting q* is quantitatively close to experimental values of qCO.

Fig. 2: Intrinsic features of electron-phonon (el–ph) coupling.
figure 2

a Plot of the dimensionless el–ph coupling \(| \tilde{{\tt{g}}}({\bf{q}};{{\bf{k}}}^{* }){| }^{2}\) (in units of 10−2) for phonon wavevectors q that connect the initial \({{\bf{k}}}_{{i}}^{* }\) and the scattered\({{\bf{k}}}_{{i}}^{* }\) + q state, both on the Fermi surface at 10% hole doping. The strongest coupling occurs at the axial wavevector q*. b The maximum value of \(| \tilde{{\tt{g}}}({\bf{q}};{\bf{k}}){| }^{2}\) with respect to all initial momenta ki for axial wavevectors q = (qx, 0). The global maximum is achieved at q* for the initial state at \({{\bf{k}}}_{{i}}^{* }\). c Spectral function Ay(k, εF) (in units of 10−1 eV−1) for the oxygen py-orbital electron on the Fermi surface. d The variation of Ay(k, εF) along the Fermi surface parametrized by the angle ψk indicated in panel (c).

Full size image

Fixing an initial state ki on the Fermi surface, we evaluate the maxima of ∣g(q; ki)∣2 with respect to q where ki + q is the final state on the Fermi surface, using the Nelder–Mead51 gradient approximation. We repeat this procedure as we vary the initial state ki along the Fermi surface branch in the first BZ quadrant and thereby identify the global maximum which we denote as \(| {\tt{g}}({{\bf{q}}}^{* };{{\bf{k}}}_{{i}}^{*}){| }^{2}\). A plot of \(| \tilde{{\tt{g}}}({\bf{q}};{{\bf{k}}}_{{i}}^{* }){| }^{2}\) vs. phonon wavevector q is shown in Fig. 2a. The strongest scattering occurs at the axial q* indicated by the white arrow.

In order to determine what is special about \({{\bf{k}}}_{{i}}^{* }\) and q*, we show the oxygen py-orbital resolved spectral function Ay(k, εF) in Fig. 2c. The highest spectral weight is obtained for \({\bf{k}}={{\bf{k}}}_{{i}}^{* }\). In Fig. 2d, we parametrize the position on the Fermi surface by the angle \({\psi }_{{\bf{k}}}={\tan }^{-1}(| {k}_{y}/{k}_{x}| )\). This panel indicates that it is the oxygen orbital content on the Fermi surface that determines \({{\bf{k}}}_{{i}}^{* }\) and q*. The variation of the oxygen content is specifically controlled by the indirect hopping processes via the Cu-4s orbital.

Landau free energy and linked cluster expansion

We next calculate the doping evolution of q* and \({{\bf{k}}}_{{i}}^{* }\) and collect the results in Fig. 3. The magnitude of q* comes close to the observed CDW ordering wavevector 0.3 (r.l.u.) and decreases with hole doping in a similar fashion as detected in X-ray experiments1,11,12,15,16,25. It is therefore tempting to suspect a close connection. To pursue this idea, we return to the el–ph Hamiltonian Eq. (4) and assume a static mean-field lattice distortion

$$\sqrt{\frac{\hslash }{2m{\Omega }_{{\rm{P}}}}}\left\langle {a}_{{\bf{q}}}+{a}_{-{\bf{q}}}^{\dagger }\right\rangle ={\xi }_{{\bf{q}}},$$
(6)

with the four possible axial wavevectors \({\pm} {\!}{{\bf{q}}}^{* },{\pm} {\!}{\overline{{\bf{q}}}}^{* }\) oriented either along the x- or the equivalent y-direction. We perform a linked-cluster expansion for the free energy to low orders in Δq = eEzξq (see Supplementary Note 3, Fig. 2 and Melikyan et al.52,53,54):

$${\mathcal{F}}= \, {{\mathcal{F}}}_{0}+\sum _{{\bf{q}} = {{\bf{q}}}^{* },{\overline{{\bf{q}}}}^{* }}\left[\frac{| {\Delta }_{{\bf{q}}}{| }^{2}}{4{\gamma }^{2}}\left(\hslash {\Omega }_{{\rm{P}}}-2{\chi }_{{\bf{q}}}^{({\tt{g}})}\right)+| {\Delta }_{{\bf{q}}}{| }^{4}{\chi }_{{\bf{q}}}^{(4)}\right]\\ + \, | {\Delta }_{{{\bf{q}}}^{* }}{| }^{2}| {\Delta }_{{\overline{{\bf{q}}}}^{* }}{| }^{2}{\chi }_{{{\bf{q}}}^{* };{\overline{{\bf{q}}}}^{* }}^{(4)}.$$
(7)

The static susceptibility \({\chi }_{{\bf{q}}}^{({\tt{g}})}\) is defined as

$${\chi }_{{\bf{q}}}^{({\tt{g}})}=-2{\gamma }^{2}\sum _{{\bf{k}}}| \tilde{{\tt{g}}}({\bf{q}};{\bf{k}}){| }^{2}\left[\frac{f({\varepsilon }_{{\bf{k}}})-f({\varepsilon }_{{\bf{k}}+{\bf{q}}})}{{\varepsilon }_{{\bf{k}}}-{\varepsilon }_{{\bf{k}}+{\bf{q}}}}\right],$$
(8)

where f(ε) denotes the Fermi distribution function. Upon cooling from high temperature, a lattice instability occurs at the q value for which the coefficient \((\hslash {\Omega }_{{\rm{P}}}-2{\chi }_{{\bf{q}}}^{({\tt{g}})})\) in the quadratic term first vanishes. This instability will necessarily produce an incommensurate charge modulation at the same q. (In previous theoretical work, such a criterion was successfully employed to identify the correct CDW wavevectors for weakly correlated tellurides53.) Upon further cooling, the fourth-order coefficients decide between uniaxial and biaxial charge order. Because the magnitudes of the coefficients, and even the sign of \({\chi }_{{{\bf{q}}}^{* };{\overline{{\bf{q}}}}^{* }}^{(4)}\), depend sensitively on specific parameter choices, it is difficult to make universal statements about the behavior of the quartic terms.

Fig. 3: Doping variation of phonon momentum.
figure 3

Hole doping (p) dependence of the axial wavevector q* (circles) and the angle \({\psi }_{{{\bf{k}}}_{{i}}^{* }}\) (diamonds) for the initial electronic momentum \({{\bf{k}}}_{{i}}^{* }\) at which the global maximum of ∣g(q; k)∣2 is achieved. The error bars represent the k-space resolution.

Full size image

Returning to the quadratic term, we define for comparison the Lindhard function, \({\chi }_{{\bf{q}}}^{({\rm{L}})}\), obtained by setting \(\tilde{{\tt{g}}}({\bf{q}};{\bf{k}})=1\) in Eq. (8); as shown in Fig. 4a, its weight is concentrated near the (π, π) point without any prominent wavevector related to nesting. Yet, as demonstrated in previous works29,30, the momentum dependence of the el–ph coupling can by itself select the wavevector for a lattice instability and the concomitant charge order. Indeed, if the structure of \(\tilde{{\tt{g}}}({\bf{q}};{\bf{k}})\) is incorporated as in Eq. (8), the momentum–space weight of \({\chi }_{{\bf{q}}}^{({\tt{g}})}\) redistributes (Fig. 4b). But, a clear instability wavevector still cannot be seen, and the largest values of \({\chi }_{{\bf{q}}}^{({\tt{g}})}\) near the 2 meV are an order of magnitude smaller to meet the instability condition \(2{\chi }_{{\bf{q}}}^{({\tt{g}})}=\hslash {\Omega }_{{\rm{P}}}=40\) meV. In particular, the axial wavevector q* for the global maximum of \(| \tilde{{\tt{g}}}({\bf{q}};{\bf{k}}){| }^{2}\) remains invisible.

Fig. 4: Momentum space structure of free susceptibility.
figure 4

a Lindhard function \({\chi }_{{\bf{q}}}^{({\rm{L}})}/{\gamma }^{2}\) and b \({\chi }_{{\bf{q}}}^{({\tt{g}})}/{\gamma }^{2}\) (in units of (eV)−1) in the first quadrant of the Brillouin zone at 110 K and 10% hole doping. Color scheme—the absolute values with the minimum (blue) and the maximum (black).

Full size image

The above analysis for the el–ph Hamiltonian Eq. (4) apparently fails to find an instability at the axial candidate wavevector q*. This naturally compromises the attempt to carry over a strategy, which so far has relied only on the electronic input parameters from density–functional theory (DFT), to strongly correlated cuprates. Seeking for the role the B1g phonon in the CDW formation therefore requires to include quasiparticle renormalization effects beyond the DFT framework. Yet, including phonons in a multi-orbital model with strong electronic correlations is a demanding task. We therefore proceed with a trial ansatz to include correlation physics on a phenomenological level.

Phenomenlogical electronic correlation

The prevailing issue in the physics of underdoped cuprates is the conundrum of the pseudogap55. The Fermi surface appears truncated to Fermi arcs centered around the BZ diagonals (the nodal regions) and it is obliterated near the BZ boundaries (the antinodal regions). In essence, well defined quasiparticles exist for near-nodal directions, while they are wiped out towards the antinodes. This motivates a simple ansatz for the Green function of the anti-bonding band quasiparticles, similar in spirit to Melikyan et al.52

$$G({\bf{k}},\omega )=\frac{{Z}_{{\bf{k}}}}{\omega -{\varepsilon }_{{\bf{k}}}+i\delta }+{G}_{{\rm{incoh}}},$$
(9a)
$${Z}_{{\bf{k}}}=\left\{\begin{array}{l}\frac{1-\cos \left({\psi }_{{\bf{k}}}-{\psi }_{\max }\right)}{1-\cos \left(4{5}^{\circ }-{\psi }_{\max }\right)};\quad {\psi }_{{\bf{k}}}\ge 4{5}^{\circ }\\ \frac{1-\cos \left({\psi }_{{\bf{k}}}-{\psi }_{\min }\right)}{1-\cos \left(4{5}^{\circ }-{\psi }_{\min }\right)};\quad {\psi }_{{\bf{k}}} \, < \, 4{5}^{\circ }\end{array},\right.$$
(9b)

for k and ψk defined in the first quadrant of the BZ (see Fig. 2c). Zk is the quasiparticle weight tailored to continuously decrease from 1 at the nodal point to zero at the BZ faces (see Supplementary Fig. 3 for a pictorial illustration of Zk). \({\psi }_{{\rm{max/min}}}\) denotes the largest/smallest angle ψk for the Fermi surface momenta (kF, π) and (π, kF), respectively. The k-dependence in Eq. (9b) is analogously carried over to the other segments of the Fermi surfaces. We emphasize that Eq. (9b) is an ad hoc, phenomenological ansatz. The desired quantitative information—in particular with sufficient momentum–space resolution—for the anisotropic quasiparticle-weight renormalization at different doping levels is not available. The latter is a hard task by itself for strongly correlated electron theory.

The ansatz in Eqs. (9a) and (9b) leads to the renormalized static susceptibility

$${\chi }_{{\bf{q}}}^{({\tt{g}}Z)}\approx -2{\gamma }^{2}\sum _{{\bf{k}}}{Z}_{{\bf{k}}}{Z}_{{\bf{k}}+{\bf{q}}}{\left|\tilde{{\tt{g}}}({\bf{q}};{\bf{k}})\right|}^{2}\left[\frac{f({\varepsilon }_{{\bf{k}}})-f({\varepsilon }_{{\bf{k}}+{\bf{q}}})}{{\varepsilon }_{{\bf{k}}}-{\varepsilon }_{{\bf{k}}+{\bf{q}}}}\right],$$
(10)

where the contributions from the incoherent part of G(k, ω) are neglected (see Supplementary Note 4).

We examine the impact of the quasiparticle weight factors on the susceptibility \({\chi }_{{\bf{q}}}^{({\tt{g}}Z)}\) in Fig. 5a. The highly anisotropic variation of Zk reduces the susceptibility around the BZ diagonal and creates new structures along the axes including the global maximum at axial wavevectors \({{\bf{q}}}_{{\rm{CO}}}^{(1)}\), and \({\overline{{\bf{q}}}}_{{\rm{CO}}}^{(1)}\), which are larger than, but close to, q* or \({\overline{{\bf{q}}}}^{* }\) for which the el–ph coupling is strongest.

Fig. 5: Momentum space structure of renormalized susceptibility.
figure 5

a Renormalized static susceptibility \({\chi }_{{\bf{q}}}^{({\tt{g}}Z)}/{\gamma }^{2}\) and b \({\chi }_{{\bf{q}}}^{({\tt{g}}Z)}/{\gamma }^{2}\) (in units of (eV)−1) modified further by a renormalized dispersion in the square brackets of Eq. (10), in the first quadrant of the Brillouin zone at 110 K and 10% hole doping. Color scheme—the absolute values with the minimum (blue) and the maximum (black).

Full size image

In a second step we incorporate—besides the anisotropic weight factors Zk—also the correlation induced renormalization of the band dispersion by replacing \({\varepsilon }_{{\bf{k}}}\Rightarrow {\varepsilon }_{{\bf{k}}}^{{\rm{R}}}\) in Eq. (9a) and hence also in the square bracket of Eq. (10). For this purpose, we adopt a phenomenological fit to the measured ARPES dispersion applied previously to data for optimally doped Bi-221256. Compared to the bare dispersion, \({\varepsilon }_{{\bf{k}}}^{{\rm{R}}}\) has an almost three times narrower bandwidth and the nodal Fermi velocity vF is similarly reduced, while the shape of the Fermi surface remains almost preserved. The susceptibility \({\chi }_{{\bf{q}}}^{({\tt{g}}Z)}\) with the renormalized \({\varepsilon }_{{\bf{k}}}^{{\rm{R}}}\) is shown in Fig. 5b. We observe that the axial peaks are retained at wavevectors \({{\bf{q}}}_{{\rm{CO}}}^{(2)}\) with \({q}^{* } \, < \, {q}_{{\rm{CO}}}^{(2)} \, < \, {q}_{{\rm{CO}}}^{(1)}\), whereas the peak at the (π, π) point loses its strength. \({q}_{{\rm{CO}}}^{(2)}\) is about 16% larger than q* and follows the same doping dependence as q*. Furthermore, the strength of the peak at \({{\bf{q}}}_{{\rm{CO}}}^{(2)}\) has tripled. This increase is naturally tied to the downward renormalization of vF. So, we conclude that it apparently requires both a correlation-induced quasiparticle renormalization, and a specific momentum dependence of the el–ph coupling to generate the axial and incommensurate candidate wavevector \({{\bf{q}}}_{{\rm{CO}}}^{(2)}\) for the anticipated lattice and concomitant CDW instability.

In reaching our conclusion, we have manifestly neglected the life-time broadening effects on the nodal quasiparticles; material difference matter. While the majority of ARPES experiments, based on Bi-based cuprates57, report a quantitatively larger quasiparticle broadening (~10 meV near the gap node58), state-of-the-art transport experiments provide much smaller broadening in YBCO. Notably, microwave measurements in YBa2Cu3O6.5 find transport scattering rates of order 0.1 meV59, which is 100 times smaller than typical values quoted for Bi-2212, which is most likely connected to the inhomogeneity of the samples. Indeed, STM experiments60,61 show that Bi based cuprates (Bi-2201 and Bi-2212) host a relatively large material inhomogeneity, while YBCO is believed to be homogeneous. It seems to us likely that scattering rates obtained from ARPES in Bi-2212 include significant disorder broadening at low temperature and energy, and are not relevant to YBCO.

Numerical estimate of B1g phonon softening

Without aiming at a quantitative description we nevertheless translate the obtained result into an estimate. The absolute units in Fig. 5b indicate that \({\chi }_{{{\bf{q}}}_{{\rm{CO}}}^{(2)}}^{{\tt{g}}Z}\ll \hslash {\Omega }_{{\rm{P}}}\), and we are therefore far from meeting the mean-field instability criterion. Nonetheless, the el–ph coupling inevitably also alters the frequency of the participating phonons. Based on Grüner et al.62, the renormalized B1g phonon frequency \({\tilde{\omega }}_{{\bf{q}}}\) follows from

$${\tilde{\omega }}_{{\bf{q}}}^{2}={\Omega }_{{\rm{P}}}^{2}(1-{\lambda }_{{\bf{q}}}),\quad {\lambda }_{{\bf{q}}}=\frac{2{\chi }_{{\bf{q}},{\Omega }_{{\rm{P}}}}^{({\tt{g}}Z)}}{\hslash {\Omega }_{{\rm{P}}}},$$
(11)

where \({\chi }_{{\bf{q}},{\Omega }_{{\rm{P}}}}^{({\tt{g}}Z)}\) is the real part of the susceptibility at the bare phonon frequency ΩP. At 110 K and \({\bf{q}}={{\bf{q}}}_{{\rm{CO}}}^{(2)}\), the dynamical susceptibility is slightly smaller than its static zero-frequency value, from which Eq. (11) gives \({\tilde{\omega }}_{{\bf{q}}}\simeq 0.987{\Omega }_{{\rm{P}}}\), i.e., a softening of 1.3%. Upon cooling \({\chi }_{{\bf{q}},{\Omega }_{{\rm{P}}}}^{({\tt{g}}Z)}\) increases slightly and the softening reaches  ~1.5% at around 3 K. For reference we mention that in Raichle et al.26, the measured softening of the B1g phonon frequency in YBa2Cu3O7 was about 6% at 3 K for q ~(0, 0.3).

Several factors may act to enhance the el–ph coupling in underdoped cuprates. We recall that a dispersion \({\varepsilon }_{{\bf{k}}}^{{\rm{R}}}\) for an optimally doped material was used for the evaluation of \({\chi }_{{\bf{q}}}^{({\tt{g}}Z)}\). On the underdoped side, the nodal Fermi velocity in Bi-2212 drops by as much as 50%63. Such a drop necessarily enhances \({\chi }_{{{\bf{q}}}_{{\rm{CO}}}^{(2)}}^{({\tt{g}}Z)}\) and the corresponding phonon softening is estimated to rise to about 2.12% reflecting somewhat enhanced CDW correlations at wavevectors \({{\bf{q}}}_{{\rm{CO}}}^{(2)}\) for underdoped materials. Furthermore, strong correlations in the underdoped cuprates decrease the copper and increase the oxygen orbital content on the Fermi surface, thereby enlarging the susceptibility at \({{\bf{q}}}_{{\rm{CO}}}^{(2)}\). One may argue that a considerably larger charge susceptibility—as an outcome of our calculation—would strengthen the case for the proposed mechanism. Yet, any theory for the CDW in Y- or Bi-based cuprates has to be reconciled with the prevailing fact that ARPES experiments fail to identify signatures of spectral changes or electronic reconstruction in the charge ordered phase. Even, in charge-stripe ordered Nd-doped LSCO, the interpretation of line-shape changes remains subtle64,65. For these reasons, the anticipated CDW can be considered weak.

The rough estimates presented above for an only weak phonon-based tendency towards charge order is—in this sense —not contradicting, but it is clearly too weak to generate a long range ordered CDW state. Still, the charge order in underdoped cuprates is in fact short ranged with moderate planar correlation lengths14,15,16, and is most likely nucleated by defects10,31,66. A strong magnetic field67,68 or uniaxial pressure27 is needed to enhance the planar correlation length and to even achieve 3D CDW order. The magnetic field suppresses superconductivity, which otherwise stops the charge correlations from growing upon cooling.

Discussion

Our proposed mechanism therefore appears compatible with experimental observations. The momentum-dependence of the el–ph coupling to the B1g bond-buckling phonon and the specific variation of the oxygen orbital content on the Fermi surface select an incommensurate, axial wavevector q* for which the lattice is most susceptible to deform. But only in conjunction with the strong and anisotropic renormalization of the correlated electrons in the copper oxygen planes does the corresponding susceptibility develop the required enhancement near q* to move the el–ph systems at least towards a charge ordered state with an axial wavevector \({{\bf{q}}}_{{\rm{CO}}}^{(2)}\) near q*.

Compared to bilayer materials like YBCO the case of single-layer cuprates is more complicated. In pristine Hg-1201 and Bi-2201, the CuO2 planes have mirror symmetries that prohibit a linear coupling between the B1g phonon and the electrons. However, both of these materials are doped by large concentrations of interstitial oxygen atoms that reside above the CuO2 plane. These unscreened dopants are the source of large electric fields in the CuO2 plane and act as nucleation sites for CDW patches69. In bilayer cuprates viz. Bi-2212, such interstitial oxygen atoms generate an electric field of the order of few eV/Å which in turn strongly amplifies (up to a factor of 5) the strength of the B1g el–ph coupling70. These field strengths are comparable to that obtained for YBCO, and must (by symmetry) produce a linear coupling between the B1g phonon and the CDW. Based on these empirical facts, we expect a similar strong enhancement of the B1g coupling in the single-layer cuprates, such as Bi-2201 and Hg-1201.

An important question is the extent to which such inhomogeneous el–ph coupling is visible in the phonon dispersion. Provided the phonons are harmonic, the lattice distortion associated with the local electric fields should not shift the phonon frequencies; rather, the el–ph matrix element will be distributed across a range of values leading to a broadening of the phonon dispersion near the ordering wavevector. However, the scale of the broadening, and in particular, how it compares with other sources of apparent broadening, as discussed in Ahmadova et al.42, remains unsettled.

A determination of boundaries for the CDW phase in the temperature vs. doping phase diagram is beyond the scope of our current work, but we offer arguments for the relevant ingredients which determine the variation with respect to hole doping. The magnitude of the charge susceptibility at the anticipated ordering wavevector is controlled by (i) the magnitude of the el–ph coupling for initial and scattered electron momenta on the Fermi surface, (ii) the nodal Fermi velocity, and (iii) the quasiparticle weight at those Fermi surface points where g(q; k) achieves its largest values in the near nodal regions.

First, the global maximum value of g(q; k) on the doping dependent Fermi surfaces drops with hole doping, i.e. upon leaving the underdoped region toward the overdoped side (see Supplementary Note 5, Fig. 4). Simultaneously, the nodal Fermi velocity almost doubles63. Both of these trends reduce the charge susceptibility.

The ansatz for the quasiparticle weight along the Fermi surface reflects the Fermi arc formation in underdoped cuprates; the quasiparticles are assumed intact in the near-nodal region with only a weak quasiparticle weight reduction. Necessarily we have to expect that the quasiparticle-weight—even in the near-nodal regions—will shrink upon lowering the hole doping concentration towards the insulator. Taken together, these trends naturally suggest a dome shaped CDW region in the temperature vs. doping phase diagram centered around an underdoped composition. This qualitative reasoning complies with the experimental findings.

For the intra-unit cell symmetry of the charge redistribution in the CDW state, our proposed scenario leads to a predominant s-symmetry form factor. Since, Cu-d-orbital character is admixed to the near-nodal Fermi surface points connected by q*, also the charge on the Cu ion will be sizably modulated, as indeed is detected by nuclear magnetic resonance8. Although a prominent d-wave character for the charge modulations on the oxygen p-orbitals has been reported in early resonant X-ray scattering71 and STM experiments72, this initial conclusion has recently been disputed. Instead, the X-ray data in McMahon et al.73 rather support a dominant s-wave form factor and are therefore compatible with the prediction from the phonon scenario.

At the core of our proposal is that it requires both, correlated electron physics and the coupling to the lattice degrees of freedom to address the CDW in cuprates. The special momentum–space structure of the el–ph coupling matrix element for the B1g bond-buckling phonon has revealed an important ingredient which was not appreciated in electronic theories before. This is the variation of the oxygen orbital content on the Fermi surface, which dictates for which momenta the coupling, here to the B1g phonon, is strongest. This may prove as a relevant step forward to elucidate the true complexity of the CDW phenomenon in cuprates. Strong electronic interactions by themselves develop charge correlations, but we infer that these may lock into an incommensurate charge-ordering pattern in Y- and Bi-based cuprates only in conjunction with a specific momentum-dependent microscopic coupling to the lattice.