Emergent Charge Order from Correlated Electron-Phonon Physics in Cuprates

TL;DR

This paper proposes that charge-density wave order in underdoped cuprates arises from electron-phonon interactions, specifically a softening of a bond-buckling phonon, influenced by strong electronic correlations, leading to local charge order.

Contribution

It introduces a mechanism linking phonon softening and electronic correlations to charge order in cuprates, challenging purely electronic explanations.

Findings

Charge-density wave is linked to bond-buckling phonon softening.

Strong correlations enhance charge susceptibility near specific phonon wavevectors.

Lattice instability and charge order can nucleate locally due to defects or disorder.

Abstract

Charge-density wave order is now understood to be a widespread feature of underdoped cuprate high-temperature superconductors, although its origins remain unclear. While experiments suggest that the charge-ordering wavevector is determined by Fermi-surface nesting, the relevant sections of the Fermi surface are featureless and provide no clue as to the underlying mechanism. Here, focusing on underdoped YBa$_2$Cu$_3$O$_{6+x}$, we propose a scenario that traces the charge-density wave formation to the incipient softening of a bond-buckling phonon. The momentum dependence of its coupling to the electrons in the copper-oxygen planes favourably selects the incommensurate and axial ordering wavevector found in experiments. But, it requires strong electronic correlations via their cuprate-specific renormalization of the weight and the dispersion of quasiparticles to enable a unique enhancement of the charge susceptibility near the B$_{1g}$-phonon selected wavevector. The frequency of the B$_{1g}$ phonon softens by a few percent, and a lattice instability with concomitant finite-range charge-density wave correlations will form locally, if nucleated by defects or dopant disorder. These results offer the perspective that the complex phase diagram of underdoped cuprates cannot be understood in the context of strong electronic correlations alone.

Equations

$$H_{kin}=\sum _{\mathbf{k},\sigma }{\mathbf{\Psi }}_{\mathbf{k},\mathbf{\sigma }}^{†}\left(\begin{array}{ccc}{\epsilon }_d& 2t_{pd}{s}_x& -2t_{pd}{s}_y\\ 2t_{pd}{s}_x& {\tilde{\epsilon }}_x(\mathbf{k})& 4{\tilde{t}}_{pp}{s}_x{s}_y\\ -2t_{pd}{s}_y& 4{\tilde{t}}_{pp}{s}_x{s}_y& {\tilde{\epsilon }}_y(\mathbf{k})\end{array}\right){\mathbf{\Psi }}_{\mathbf{k},\sigma },$$

$${\tilde{\epsilon }}_{x,y}={\epsilon }_p+4t_{pp}^i{s}_{x,y}^2;{\tilde{t}}_{pp}=t_{pp}^i+t_{pp}^d;t_{pp}^i=\frac{t_{ps}^2}{{\epsilon }_F-{\epsilon }_s},$$

$$\mathcal{H}_{e\-p}=eE_{z}\sum_{\mathbf{n}}\Big{[}\delta u_{x\mathbf{n}}p^{\dagger}_{x\mathbf{n}}p_{x\mathbf{n}}+\delta u_{y\mathbf{n}}p^{\dagger}_{y\mathbf{n}}p_{y\mathbf{n}}\Big{]},$$

$$\scalebox{0.9}[1.0]{$\mathcal{H}_{\textit{el-ph}}=\sum_{\mathbf{k}}\varepsilon_{\mathbf{k}}c^{\dagger}_{\mathbf{k}}c_{\mathbf{k}}+\sum_{\mathbf{q},\mathbf{k}}\mathsf{g}(\mathbf{q};\mathbf{k})c^{\dagger}_{\mathbf{k}+\mathbf{q}}c_{\mathbf{k}}\left(a_{\mathbf{q}}+a^{\dagger}_{-\mathbf{q}}\right)$},$$

$$\tilde{\mathsf{g}}(\mathbf{q};\mathbf{k})=\Big{[}e^{x}_{\mathbf{q}}\phi_{x}(\mathbf{k}^{\prime})\phi_{x}(\mathbf{k})+e^{y}_{\mathbf{q}}\phi_{y}(\mathbf{k}^{\prime})\phi_{y}(\mathbf{k})\Big{]},$$

$$\sqrt{\frac{\hslash }{2m{\Omega }_P}}⟨a_{\mathbf{q}}+a_{-\mathbf{q}}^{†}⟩={\xi }_{\mathbf{q}},$$

$$\mathcal{F}={\mathcal{F}}_0+$$

$$+|{\Delta }_{{\mathbf{q}}^{\ast }}{|}^2|{\Delta }_{{\overline{\mathbf{q}}}^{\ast }}{|}^2{\chi }_{{\mathbf{q}}^{\ast };{\overline{\mathbf{q}}}^{\ast }}^{(4)}.$$

$$\chi^{(\mathsf{g})}_{\mathbf{q}}=-2\gamma^{2}\sum_{\mathbf{k}}|\tilde{\mathsf{g}}(\mathbf{q};\mathbf{k})|^{2}\Bigg{[}\frac{f(\varepsilon_{\mathbf{k}})-f(\varepsilon_{\mathbf{k}+\mathbf{q}})}{\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{k}+\mathbf{q}}}\Bigg{]},$$

$$G(\mathbf{k},\omega )=\frac{Z_{\mathbf{k}}}{\omega -{\epsilon }_{\mathbf{k}}+i\delta }+G_{\text{incoh}},$$

$$Z_{\mathbf{k}}$$

$$\scalebox{0.9}[1.0]{$\chi^{(\mathsf{g}Z)}_{\mathbf{q}}\approx-2\gamma^{2}\sum_{\mathbf{k}}Z_{\mathbf{k}}Z_{\mathbf{k}+\mathbf{q}}|\tilde{\mathsf{g}}(\mathbf{q};\mathbf{k})|^{2}$}\Bigg{[}\frac{f(\varepsilon_{\mathbf{k}})-f(\varepsilon_{\mathbf{k}+\mathbf{q}})}{\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{k}+\mathbf{q}}}\Bigg{]},$$

$${\tilde{\omega }}_{\mathbf{q}}^2={\Omega }_P^2(1-{\lambda }_{\mathbf{q}}),{\lambda }_{\mathbf{q}}=\frac{2{\chi }_{\mathbf{q},{\Omega }_P}^{(\mathsf{g}Z)}}{\hslash {\Omega }_P},$$

Full text

I Introduction

The discovery of charge-density wave (CDW) order Comin and Damascelli (2016) in underdoped cuprates raised the question of whether it is intimately related to pseudogap physics Caprara et al. (2017); Verret et al. (2017); Chatterjee and Sachdev (2016); Efetov et al. (2013); Tsvelik and Chubukov (2014); Atkinson et al. (2015), and thereby yet another signature of strong electronic correlations. Charge modulations with a moderate correlation length are detected by nuclear magnetic resonance Wu et al. (2011, 2013, 2015), scanning tunneling microscopy Comin et al. (2014); da Silva Neto et al. (2014); Kohsaka et al. (2007), and x-ray techniques Ghiringhelli et al. (2012); Comin et al. (2014); da Silva Neto et al. (2014); Chang et al. (2012); Blackburn et al. (2013), with the latter finding incommensurate wave-vectors near $\mathbf{q}_{CO}=0.3$ reciprocal lattice units, oriented along the crystalline axes with either uniaxial or biaxial character. The CDW wavevector $\mathbf{q}_{CO}$ continuously drops with increasing hole doping in YBa2Cu3O6+x (YBCO) and Bi2Sr2-xLaxCuO6+x (Bi-2201) Comin and Damascelli (2016); Comin et al. (2014); Blackburn et al. (2013); Peng et al. (2016), whereas a completely opposite trend is observed in charge-stripe ordered La-based 214 cuprates such as La2-xSrxCuO4 (LSCO) and La2-xBaxCuO4 (LBCO) Hücker et al. (2011). Not only the concomitant spin-stripe order marks 214 cuprates to be distinctly different from YBCO or Bi-2201; also the tilt-pattern of oxygen octahedra, specifically in the low-tempertaure tetragonal (LTT) phase, points to subtle structural differences, since the LTT tilt pattern breaks the four-fold symmetry in the copper-oxygen planes and thereby offers a structure related source for a uniaxial character of the density waves. Additionally, the orbital symmetry of CDW order in YBCO and La1.875Ba0.125CuO4 was found to be dissimilar Achkar et al. (2016). Due to these apparent differences we focus subsequently on the charge ordering in YBCO.

Stripe formation and the competition of stripe states with d-wave superconductivity or pair-density waves was investigated for multiband Hubbard or $t-J$ models with accurate computational tools Corboz et al. (2014); Zheng et al. (2017); Jiang and Devereaux (2019). Yet, finding the charge-stripe ordered ground states proved elusive due to the near degeneracy of the competing states. So far the axial orientation of $\mathbf{q}_{CO}$ and its doping variation in YBCO and Bi-2201 proved difficult to reconcile with purely electronic model calculations Huang et al. (2017); Zheng et al. (2017). Weak-coupling theories typically predict that $\mathbf{q}_{CO}$ lies on the Brillouin zone (BZ) diagonal Efetov et al. (2013); Sau and Sachdev (2014); Bulut et al. (2013) unless the CDW instability is preceded by a Fermi surface reconstruction to form hole pockets Atkinson et al. (2015). Based on these evidences, we conjecture that some physics ingredient may still be missing in these model calculations.

In the continuing search for the origin of charge order in cuprates, we employ further experimental facts. One is the discovery of charge order in overdoped Bi-2201 with an unreconstructed Fermi surface Peng et al. (2018). This result poses a question as to whether the CDW is actually tied to the pseudogap phenomena. A second experimental hint relates to the softening or broadening of phonon modes at the CDW wavevector Raichle et al. (2011); Kim et al. (2018) and to giant phonon anomalies near the CDW instability, which all point to a strong electron-phonon (el-ph) coupling Le Tacon et al. (2013). 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 materials Eiter et al. (2013); Johannes and Mazin (2008). In cuprates, the frequencies of the phonons decrease only weakly, and a continuous softening to zero frequency is likely precluded by quenched disorder Campi et al. (2015); Wu et al. (2015).

A third hint comes from the atomic displacement pattern that accompanies the CDW. There appears to be only one complete data set for the static lattice distortion pattern in charge-ordered underdoped YBCO. X-ray diffraction data in Ref. Forgan et al. (2015) 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 phonon Cuk et al. (2004). In stoichiometric, overdoped YBCO the dispersion of the B1g phonon was monitored upon cooling Raichle et al. (2011); 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 Ref. Forgan et al. (2015) 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. Le Tacon et al. (2013) and Kim et al. Kim et al. (2018). 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. Besides small additional in-plane displacements of the oxygen atoms, the x-ray experiments Forgan et al. (2015) clearly signal that the lattice distortions in the CDW phase have predominant B1g character.

Based on these empirical grounds, we select the B1g bond-buckling phonon exclusively, reanalyze its coupling to the electrons in the copper-oxygen planes, and pursue first 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 $\mathsf{g}(\mathbf{q};\mathbf{k})$ depends strongly on the orbital content of the Fermi surface. In particular, the Cu-4s orbital is crucial; with it, $\mathsf{g}(\mathbf{q};\mathbf{k})$ is maximum at phonon momenta $\mathbf{q}^{*}$ that are axial (rather than diagonal) and track the doping dependence of the CDW wavevector. The special role of $\mathbf{q}^{*}$ is imperceptible in the free energy when calculated with bare electron dispersions; however, with electronic correlations, modeled by a renormalization of the band dispersion and quasiparticle spectral weight, an axial wavevector close to $\mathbf{q}^{*}$ emerges as the dominant wavevector in the CDW susceptibility. The $B_{1g}$ mode is, by itself, too weak to induce true long-range order; but 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.

II 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 work Maharaj et al. (2014) 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 doping Fournier et al. (2010). 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 $3d_{x^{2}-y^{2}}$ as well as oxygen $p_{x}$ and $p_{y}$ orbitals. An effective three-band model is obtained in terms of the $d$- and $p$- orbitals by downfolding the original four-band model to Andersen et al. (1995); Atkinson et al. (2015)

[TABLE]

with the three-spinor $\bm{\Psi^{\dagger}_{\mathbf{k},\sigma}}=\left(d^{\dagger}_{\mathbf{k},\sigma},\;p^{\dagger}_{x\mathbf{k},\sigma},\;p^{\dagger}_{y\mathbf{k},\sigma}\right)$ and $s_{x,y}=\sin(k_{x,y}/2)$. $\varepsilon_{d}$ denotes the onsite energy of the $d$ orbital, $t_{pd}$ and $t_{ps}$ 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 $\varepsilon_{p}$ and generate indirect hopping $4t^{i}_{pp}$ between oxygen orbitals:

[TABLE]

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

The out-of-plane $B_{1g}$ vibrations of the oxygen atoms, as depicted in Fig. 1, naturally couple linearly to their local electric field $E_{z}$ Devereaux et al. (1995). At this point, we neglect the motion of the almost four times heavier copper atoms. Consequently, we start from the ansatz

[TABLE]

where $\delta u_{x,y\mathbf{n}}$ are the out-of-plane displacements of the oxygen atoms in unit cell $\bm{n}$. We project the el-ph Hamiltonian in Eq. 3, onto the anti-bonding band and obtain the effective Hamiltonian $\mathcal{H}=\mathcal{H}_{\textit{el-ph}}+\mathcal{H}_{\textit{ph}}$ with

[TABLE]

where $c^{\dagger}_{\mathbf{k}}$ is the creation operator for the anti-bonding electrons with dispersion $\varepsilon_{\mathbf{k}}$, $a_{\mathbf{q}}$ annihilates a $B_{1g}$ phonon mode, and $\mathcal{H}_{\textit{ph}}=\hbar\Omega_{P}\sum_{\mathbf{q}}a^{\dagger}_{\mathbf{q}}a_{\mathbf{q}}$ (see Ref. Devereaux et al. (1995) and Supplementary Materials). The momentum-dependent el-ph coupling is written $\mathsf{g}(\mathbf{q};\mathbf{k})=\gamma\tilde{\mathsf{g}}(\mathbf{q};\mathbf{k})$ Devereaux et al. (1999), where $\gamma$ is the coupling strength and

[TABLE]

with $\mathbf{k}^{\prime}=\mathbf{k}+\mathbf{q}$. The eigenfunctions $\phi_{x,y}$ signify the orbital contents of the oxygen $p_{x,y}$ orbitals in the anti-bonding band. Similarly, the eigenvectors $e^{x,y}_{\mathbf{q}}$ 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=eE_{z}\sqrt{\hbar/2m\Omega_{P}}$, where $m$ is the mass of the oxygen atom and $\Omega_{P}\sim 40$ meV is the frequency of the dispersionless $B_{1g}$ mode. Adopting the electric-field value $eE_{z}=3.56$ eV/Å from Ref. Johnston et al. (2010) leads to $\gamma=0.22$ eV. The eigenvectors for the $B_{1g}$ mode are $e^{x,y}_{\mathbf{q}}=\mp\frac{\cos(q_{y,x}/2)}{M_{\mathbf{q}}}$, where the normalization factor $M_{\mathbf{q}}=\sqrt{\cos^{2}{(q_{x}/2)}+\cos^{2}{(q_{y}/2)}}$.

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

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

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

III Analysis and Results

We next calculate the doping evolution of $\mathbf{q}^{*}$ and $\mathbf{k}^{*}_{i}$ and collect the results in Fig. 3. The magnitude of $\mathbf{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 experiments Comin and Damascelli (2016); Comin et al. (2014); da Silva Neto et al. (2014); Chang et al. (2012); Blackburn et al. (2013); Peng et al. (2018). 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

[TABLE]

with the four possible axial wavevectors $\pm\mathbf{q}^{*},\pm\overline{\mathbf{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 $\Delta_{\mathbf{q}}=eE_{z}\xi_{\mathbf{q}}$ (see Supplementary Material and Refs. Melikyan and Norman (2014); Yao et al. (2006); Mahan (2011)):

[TABLE]

The static susceptibility $\chi^{(\mathsf{g})}_{\mathbf{q}}$ is defined as

[TABLE]

where $f(\varepsilon)$ denotes the Fermi distribution function. Upon cooling from high temperature, a lattice instability occurs at the $\mathbf{q}$ value for which the coefficient $(\hbar\Omega_{P}-2\chi^{(\mathsf{g})}_{\mathbf{q}})$ in the quadratic term first vanishes. This instability will necessarily produce an incommensurate charge modulation at the same $\mathbf{q}$. (In previous theoretical work, such a criterion was successfully employed to identify the correct CDW wavevectors for weakly correlated tellurides Yao et al. (2006).) 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^{(4)}_{\mathbf{q}^{*};\overline{\mathbf{q}}^{*}}$, depend sensitively on specific parameter choices, it is difficult to make universal statements about the behavior of the quartic terms.

Returning to the quadratic term, we define for comparison the Lindhard function, $\chi^{(\mathrm{L})}_{\mathbf{q}}$, obtained by setting $\tilde{\mathsf{g}}(\mathbf{q};\mathbf{k})=1$ in Eq. 8; as shown in Fig. 4(a), its weight is concentrated near the $(\pi,\pi)$ point without any prominent wavevector related to nesting. Yet, as demonstrated in Refs. Johannes and Mazin (2008); Eiter et al. (2013), 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{\mathsf{g}}(\mathbf{q};\mathbf{k})$ is incorporated as in Eq. 8, the momentum-space weight of $\chi^{(\mathsf{g})}_{\mathbf{q}}$ redistributes [Fig. 4(b)]. But, a clear instability wavevector still cannot be seen, and the largest values of $\chi_{\mathbf{q}}^{(\mathsf{g})}$ near the 2 meV are an order of magnitude smaller to meet the instability condition $2\chi_{\mathbf{q}}^{(\mathsf{g})}=\hbar\Omega_{P}=40$ meV. In particular, the axial wavevector $\mathbf{q}^{*}$ for the global maximum of $|\tilde{\mathsf{g}}(\mathbf{q};\mathbf{k})|^{2}$ remains invisible.

The above analysis for the el-ph Hamiltonian Eq. 4 apparently fails to find an instability at the axial candidate wavevector $\mathbf{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.

The prevailing issue in the physics of underdoped cuprates is the conundrum of the pseudogap Lee et al. (2006). 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 Ref. Melikyan and Norman (2014):

[TABLE]

for $\mathbf{k}$ and $\psi_{\mathbf{k}}$ defined in the first quadrant of the BZ (see Fig. 2c). $Z_{\mathbf{k}}$ is the quasiparticle weight tailored to continuously decrease from 1 at the nodal point to zero at the BZ faces. $\psi_{max/min}$ denotes the largest/smallest angle $\psi_{\mathbf{k}}$ for the Fermi surface momenta $(k_{F},\pi)$ and $(\pi,k_{F})$, respectively. The $\mathbf{k}$-dependence in Eq. 9d is analogously carried over to the other segments of the Fermi surfaces. We emphasize that Eq. 9d 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 Eq. 9a, 9d leads to the renormalized static susceptibility

[TABLE]

where the contributions from the incoherent part of $G(\mathbf{k},\omega)$ are neglected (see Supplementary Material).

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

In a second step we incorporate – besides the anisotropic weight factors $Z_{\mathbf{k}}$ – also the correlation induced renormalization of the band dispersion by replacing $\varepsilon_{\mathbf{k}}\Rightarrow\varepsilon^{R}_{\mathbf{k}}$ 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-2212 Norman (2007). Compared to the bare dispersion, $\varepsilon^{R}_{\mathbf{k}}$ has an almost three times narrower bandwidth and the nodal Fermi velocity $v_{F}$ is similarly reduced, while the shape of the Fermi surface remains almost preserved. The susceptibility $\chi^{(\mathsf{g}Z)}_{\mathbf{q}}$ with the renormalized $\varepsilon^{R}_{\mathbf{k}}$ is shown in Fig. 5b. We observe that the axial peaks are retained at wavevectors $\mathbf{q}^{(2)}_{CO}$ with $q^{*}<q^{(2)}_{CO}<q^{(1)}_{CO}$, whereas the peak at the $(\pi,\pi)$ point loses its strength. $q^{(2)}_{CO}$ is about $16\%$ larger than $q^{*}$ and follows the same doping dependence as $q^{*}$. Furthermore, the strength of the peak at $\mathbf{q}^{(2)}_{CO}$ has tripled. This increase is naturally tied to the downward renormalization of $v_{F}$. 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 $\mathbf{q}^{(2)}_{CO}$ 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 cuprates Li et al. (2018), report a quantitatively larger quasiparticle broadening ($\sim 10$ meV near the gap node Yamasaki et al. (2007)), 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 meV Harris et al. (2006), 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 experiments Pan et al. (2001); Fischer et al. (2007) 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.

Without aiming at a quantitative description we nevertheless translate the obtained result into an estimate. The absolute units in Fig. 5b indicate that $\chi^{\mathsf{g}Z}_{\mathbf{q}^{(2)}_{CO}}\ll\hbar\Omega_{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 Ref. Grüner (1994), the renormalized $B_{1g}$ phonon frequency $\tilde{\omega}_{\mathbf{q}}$ follows from

[TABLE]

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

Several factors may act to enhance the el-ph coupling in underdoped cuprates. We recall that a dispersion $\varepsilon^{R}_{\mathbf{k}}$ for an optimally doped material was used for the evaluation of $\chi^{(\mathsf{g}Z)}_{\mathbf{q}}$. On the underdoped side, the nodal Fermi velocity in Bi-2212 drops by as much as 50% Vishik et al. (2010). Such a drop necessarily enhances $\chi^{(\mathsf{g}Z)}_{\mathbf{q}^{(2)}_{CO}}$ and the corresponding phonon softening is estimated to rise to about $2.12\%$ reflecting somewhat enhanced CDW correlations at wavevectors $\mathbf{q}^{(2)}_{CO}$ 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 $\mathbf{q}^{(2)}_{CO}$. 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 subtle Matt et al. (2015); Kivelson and Lederer (2019). 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 lengths Ghiringhelli et al. (2012); Chang et al. (2012); Blackburn et al. (2013), and is most likely nucleated by defects Campi et al. (2015); Wu et al. (2015); Atkinson et al. (2018). A strong magnetic field Gerber et al. (2015); Jang et al. (2016) or uniaxial pressure Kim et al. (2018) 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.

IV Discussion and Concluding Remarks

Our proposed mechanism therefore appears compatible with experimental observations. The momentum-dependence of the el-ph coupling to the $B_{1g}$ bond-buckling phonon and the specific variation of the oxygen orbital content on the Fermi surface select an incommensurate, axial wavevector $\mathbf{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 $\mathbf{q}^{*}$ to move the el-ph systems at least towards a charge ordered state with an axial wavevector $\mathbf{q}^{(2)}_{CO}$ near $\mathbf{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 patches Tabis et al. (2017). 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 coupling Johnston et al. (2009). 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 Ref. Ahmadova et al. (2020), 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 $\mathsf{g}(\mathbf{q};\mathbf{k})$ achieves its largest values in the near nodal regions.

First, the global maximum value of $\mathsf{g}(\mathbf{q};\mathbf{k})$ on the doping dependent Fermi surfaces drops with hole doping, i.e. upon leaving the underdoped region towards the overdoped side (see Sec. E in Supplemental Material). Simultaneously, the nodal Fermi velocity almost doubles Vishik et al. (2010). 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 $\mathbf{q}^{*}$, also the charge on the Cu ion will be sizably modulated, as indeed is detected by nuclear magnetic resonance Wu et al. (2011). Although a prominent $d$-wave character for the charge modulations on the oxygen $p$-orbitals has been reported in early resonant x-ray scattering Comin et al. (2015) and STM experiments Fujita et al. (2014), this initial conclusion has recently been disputed. Instead, the x-ray data in Ref. McMahon et al. (2019) 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.

V Acknowledgement

We acknowledge helpful conversations with T. P. Devereaux, P. J. Hirschfeld, T. Kopp, L. Chioncel, S. Johnston, and C. Morice. WAA acknowledges support by the Natural Sciences and Engineering Research Council (NSERC) of Canada. APK acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)- project-ID-107745057-TRR 80.