Charge Order in the Pseudogap Phase of Cuprate Superconductors

In a multiorbital model of the cuprate high-temperature superconductors soft antiferromagnetic (AF) modes are assumed to reconstruct the Fermi surface to form nodal pockets. The subsequent charge ordering transition leads to a phase with a spatially modulated transfer of charge between neighboring oxygen p_x and p_y orbitals and also weak modulations of the charge density on the copper d_{x^2-y^2} orbitals. As a prime result of the AF Fermi surface reconstruction, the wavevectors of the charge modulations are oriented along the crystalline axes with a periodicity that agrees quantitatively with experiments. This resolves a discrepancy between experiments, which find axial order, and previous theoretical calculations, which find modulation wavevectors along the Brillouin zone (BZ) diagonal. The axial order is stabilized by hopping processes via the Cu4s orbital, which is commonly not included in model analyses of cuprate superconductors.

An intriguing feature of charge order in YBCO and BSCCO is that there appears to be a strong intra-unit cell transfer of charge between oxygen atoms in each CuO 2 plaquette, rather than the inter-unit cell charge transfer normally associated with charge-density waves. The most direct evidence for this comes from STM experiments [23,33,40], and further support is provided by x-ray scattering [28]. Roughly then, the charge ordered phase can be thought of as a finite-q modulation with a d x 2 −y 2 form factor describing the intra-unit cell charge transfer, and with relatively little transfer of charge between neighbouring unit cells. For this reason, the phase is sometimes called a "dCDW". Alternatively, because the charge order is a generalization of a q = (0, 0) nematic phase that breaks rotational but not translational symmetry, it has been labelled a "modulated nematic".
It is natural to ask whether the dCDW charge order identified in YBCO and BSCCO is related to stripe order (see [44,45] for further discussion). Stripe order is well established in La 2−x Ba x CuO 4 and dynamical stripes are inferred in La 2−x Sr x CuO 4 [46]. Stripes in La 2−x Ba x CuO 4 are characterized by a static or quasistatic spin modulation whose period is double that of a concomitant charge modulation. [47] The doping dependence of the modulation wavevector is opposite to what one would expect for a Fermi surface instability, and suggests instead a strong coupling picture in which holes and spins segregate into one dimensional stripes [47]. The dCDW described above has some similarities to this stripe order: both compete with superconductivity, and both have maximal intensity near a hole doping p = 1/8 in all cuprates for which the doping dependence has been measured. [48,49,50] On the other hand, there are also significant differences. First, the local spins in YBCO are dynamic, rather than (quasi)static. Models of fluctuating stripes have been proposed to describe this [46,51,52]; however, recent NMR experiments clearly show that the charge order is static in YBCO up to high temperatures. [53] This suggests that the intertwining of charge and spin textures that is key to stripe formation in the La-cuprates is not a factor in YBCO. Consistent with this we note that, while the doping dependences of the spin and charge modulation wavevectors are closely connected in La 2−x Ba x CuO 4 , they appear unconnected in YBCO [29]. Finally, recent x-ray experiments have shown that the structure factor for charge order in La 2−x Ba x CuO 4 has an extended-s symmetry, [45] consistent with multiorbital models of a magnetically driven stripe instability, [54] and in contrast to YBCO and BSCCO.
Whether these differences are due to small differences in the band structure that tip the balance towards particular phases, or point to larger differences between the cuprate families is not yet established. Here, we adopt the point of view that the mechanism driving charge order in YBa 2 Cu 3 O 6+x is distinct from that in the La-cuprates. The majority of previous theoretical work along these lines is based on one-band effective models of a single CuO 2 plane; in such models, the analogue of intra-unit cell charge redistribution is bond order, namely, an anisotropic renormalization of the electronic effective mass along the x and y axes. Several theories have argued that bond order follows from AF exchange interactions; a vital role for the charge instabilities is thereby ascribed to "hot spot" regions of the Fermi surface where scattering from AF spin fluctuations is especially strong [55,56,43,57,42,58,59,60]. Alternative one-band [61,62,63] and three-band [64] model calculations with generic interactions have found similar charge instabilities. With the exception of Ref. [59], which additionally found current-carrying stripes, these models universally obtained a charge density with a d x 2 −y 2 form factor and an ordering wavevector q * along the BZ diagonal. While the form factor is compatible with experiments [23,33,28], the magnitude of q * is typically too small by a factor of 2, and the direction of q * is rotated by 45 • relative to the experiments. The robustness of these discrepancies suggests that the underlying models lack an essential ingredient.
In this work, we show that the Fermi surface topology affects the emergent charge order in a fundamental way, and can explain the discrepancy between the observed and predicted values of q * . Starting from a simplified model of the Fermi surface in the pseudogap phase, we obtain a charge instability that quantitatively agrees with that found experimentally. The implication of this work is that charge order emerges from the pseudogap phase, rather than contributing to it directly.
Experimentally, the pseudogap is characterized by a partial depletion of the density of states around the Fermi level. Photoemission experiments have revealed that this depletion occurs near the BZ boundary at (±π, 0) and (0, ±π), and that spectral weight near these points is pushed away from the Fermi energy [65]. Early proposals ascribed this to nearly AF spin fluctuations that partially nest these regions of the Fermi surface and shift spectral weight to higher energy [2,3]. This picture remains physically appealing because the underdoped regime lies near the AF insulating phase of the parent compounds, and it is supported by quantum Monte Carlo (QMC) [66] and cluster dynamical mean-field theory [16] (cDMFT) calculations that draw a link between spin fluctuations and the pseudogap.
With this in mind, we adopt a model that we believe contains the essential ingredients to understand charge order in the cuprates. The basic element of this model is a two dimensional CuO 2 plane, and we retain both Cu and O orbitals, as well as the short-range Coulomb interactions between them. While it is likely that similar ordering instabilities to the ones described here may be found in one-band models, by choosing a multiorbital model we are able to obtain details of the intra-unit cell charge redistribution, which can be directly probed by STM, NMR, and x-ray experiments. We examine specifically the charge instabilities generated by the short-range Coulomb forces, which have been shown to be attractive in the relevant charge ordering channel [67,64]. We note that spin fluctuations are also attractive in this channel, at least in one-band models; these will modify T co , but should not alter the relationship between q * and the Fermi surface structure. For simplicity, then, we focus on the charge fluctuations only and omit spin dynamics.
Pseudogap physics in the three-orbital model derives from the large local Coulomb interaction, U d ∼ 10 eV, on the Cu sites, which is the source of strong correlation physics in the cuprates. Despite the achievements of state-of-the-art computational methods, there is still a paucity of tools available that can capture both the short-range strong correlation physics and the long-range physics of incommensurate charge order. Numerical methods like QMC and cDMFT, which have proved capable of verifying the pseudogap structures in the density of states, suffer from finite-size effects that render the charge instability inaccessible, and cluster methods further face the difficulty of treating the nonlocal Coulomb interactions that drive the charge order [68]. On the other hand, weak-coupling methods, which capture long-range physics, find that spin fluctuations introduce only weak pseudogap-like spectral features [3,4,5,6].
For these reasons, we follow a partially phenomenological approach. To leading order, the effect of U d is to suppress double occupancy of the Cu d x 2 −y 2 -orbital and thereby create local moments; as a consequence, the itinerant electrons reside primarily on the oxygen sites, although the Fermi surface does nevertheless have some Cu character due to the hybridization of Cu and O orbitals. Above the superconducting transition, the spin spectrum measured by neutron scattering [69] is centred at (π, π), indicating dynamical AF correlations. Our subsequent diagrammatic analysis is based on two strong simplifying assumptions: first, the moments are assumed quasistatic on electronic timescales, and second, the AF correlation length ξ AF is larger than the charge order correlation length ξ co , which is estimated from experiments to be ∼ 50Å [49]. In essence, this implies that the local moments can be treated, as if they are ordered antiferromagnetically. (A similar ansatz was made in [70], where static spin textures of classical magnetic moments on the copper sites were shown to induce charge order.) Both of these assumptions are not satisfied throughout most of the doping range where charge order is observed experimentally (ξ AF ∼ 20Å in YBa 2 Cu 3 O 6.5 [71]); however, by making these assumptions we are able to explore a different physical regime than previous weak coupling calculations. We also note that there is evidence in YBa 2 Cu 3 O 6+x that charge order survives inside a static magnetic phase that exists at very low doping. [49] Ultimately, however, one should think of the above approximations as a simple phenomenology for strong correlation physics on the copper sites, which is justified after the fact by the surprising accuracy with which we predict certain properties of the charge ordered phase. The effect of the ordered local moments is to create a pseudogap-like shift of spectral weight away from the Fermi level at (±π, 0) and (0, ±π) and to reconstruct the Fermi surface to form hole and electron pockets. In this pseudogapped state we find that residual Coulomb interactions between the quasiparticles can drive a d x 2 −y 2 -like charge redistribution between Op x and Op y orbitals, accompanied by a weaker periodic modulation of the Cu charge density. The obtained charge pattern with an ordering wavevector q * along the BZ axis is indeed consistent with what has been observed experimentally. This charge order induces a second Fermi surface reconstruction which generates diamond-shaped electron pockets. The existence of such pockets was earlier inferred from quantum-oscillation experiments.
We introduce the model in section 2 and describe briefly the calculations for the charge susceptibility, with details left for the appendices. The results of our calculations are discussed in section 3, with an emphasis on comparisons to experiments. The main implication is that the proposed model calculation, while still not a complete description of the microscopic physics underlying charge order, provides a route to understand the experiments, and suggests that charge order and pseudogap features are in fact distinct phenomena. A short summary is contained in section 4.

Model and Calculations
The goal is to model charge order in YBa 2 Cu 3 O 6+x , and to this end we employ a multiband description of the CuO 2 planes due to Andersen et al. (ALJP) that was derived specifically for YBa 2 Cu 3 O 7 [72]. In an extension to the Emery model [73], which is based on the Cu3d x 2 −y 2 and two σ-bonded oxygen orbitals, Op x and Op y , ALJP included also the Cu4s orbital. The latter resides well above the Fermi energy, approximately 6.5 eV above the Cud orbital, and has a large overlap with the Op orbitals. Downfolding this orbital leads to an effective three-band model (see Appendix is an array of electron creation operators for the d, p x , and p y orbitals. Parameters t pd and t pp denote hopping amplitudes, s x,y = sin(k x,y /2),˜ x,y (k) = p + 4t i pp s 2 x,y , and d and p are orbital energies. The tilde denotes renormalization by hopping through the Cu4s orbital. In particular,t pp = t d pp + t i pp where the superscripts indicate direct (d) and indirect (i; through the Cu4s orbital) hopping between Op orbitals.
For the reasons outlined above we introduce AF moments on the Cud orbitals by adding a staggered spin-dependent potential M (r j ) to the Hamiltonian and thereby obtain a pseudogap-like reconstruction of the Fermi surface. It is natural to think of this potential as the auxilliary field that appears when the Coulomb interaction U dnjd↑njd↓ on the Cud orbitals is removed by a Hubbard-Stratonovich transformation. In this transformation, the quartic interaction term is replaced by an interaction between the electrons and a spin-polarizing time-dependent auxiliary field M (r j , t). As mentioned in section 1, we make two assumptions in order to isolate the physics of interest: first, that the field is static, and second that M (r j ) has long range AF order. Under these assumptions, an additional term, −M j e iQ·r j (n j↑ −n j↓ ) with Q = (π, π), is added to the Hamiltonian. On physical grounds, we expect this potential to be quite large: within a saddle-point approximation, M = U d m, where m is the static AF moment on the Cu sites. Given that U d ∼ 10 eV in the cuprates, even a modest value of m leads to M ∼ 1 eV. The Fermi surface reconstruction generated by M is illustrated in Fig. 1, where the local Cud moments open a gap along Fermi surface segments near the AF hot spots, i.e. those points where the Fermi surface intersects the magnetic BZ boundary.
Charge order is driven by interactions between quasiparticles in the reconstructed bands. It has been shown that, in one band models at least, the exchange of spin fluctuations may drive a charge ordering transition; here, we consider only short range Coulomb interactions. Electrons interact at short distances through intra-orbital U d and U p and nearest-neighbor V pd and V pp Coulomb repulsions. The corresponding interaction part of the Hamiltonian iŝ where j implies summation over unit cells, and δ is summed over nearest-neighbor orbitals of type Op x,y (for V pd ) or Op y (for V pp ). In our model, the charge instability is driven by V pp .
To study charge ordering tendencies, we calculate the charge susceptibility χ αβ (q) = −(∂n α /∂ β )(q), where n α denotes electron densities and α and β are orbital labels. The onset of charge order is signalled by a diverging susceptibility at a specific q * upon lowering the temperature. The interactions are treated in a generalized random-phase approximation (see Ref. [64] and Appendix B), which allows one to find the leading charge instability without any bias towards a particular ordering wavevector q * or orbital type.

Results
The main results of this calculation are summarized in Fig. 1. The Fermi surface for the ALJP bands is shown in Fig. 1(a), along with the wavevectors q 1 and q 2 at which the charge susceptibility first diverges upon cooling in the absence of staggered Cu moments. As in previous calculations [56,57,42,61,62,43,64,58], these wavevectors lie along the BZ diagonals and the charge instability primarily involves an intra-unit cell charge transfer between Op x and Op y orbitals. q 1 and q 2 connect points close to nearby hotspot regions of the Fermi surface. When M is finite but small, as in Fig. 1(b), the Fermi surface breaks up into hole pockets around (±π/2, ±π/2) and electron pockets centered at the "antinodal" points on the BZ boundary; the modulation wavevectors remain diagonal and connect these pockets.
While the directions of q 1 and q 2 are consistent with previous calculations, they conflict with experiments [20,21,30,34,31,29], which clearly indicate that the charge densities are modulated along the axial Cu-O bond directions. This discrepancy is resolved when the electron pockets are fully eliminated [ Fig. 1(c)] by a sufficiently large staggered potential M and the modulation wavevectors rotate to the axial direction. Furthermore, the magnitude of q * 1,2 agrees quantitatively with the experimental data of Blackburn et al. [29] as shown in Fig. 1(d) for the doping dependence of |q * 1,2 |. We emphasize that no fine tuning of the model parameters was done to obtain these results: the band parameters were taken from Ref. [72], and q * 1 and q * 2 depend only weakly on the size of M once it is large enough to remove the electron pockets.
The fact that the Fermi surface forms well defined pockets is an artefact of the assumption that the AF correlation length ξ AF is infinite [74], and indeed there is no experimental evidence for hole pockets of the type shown in Fig. 1(c). When the ξ AF is finite, however, the pockets become arcs, which is consistent with experiments. For our purposes, it is important to note that the portions of the Fermi surface connected by q * 1 and q * 2 remain well defined when ξ AF is finite [74], while the back sides of the pockets are wiped out. For this reason, we believe that the leading charge instability described here will also be the leading instability in models with short range AF correlations.
The charge modulation amplitudes on the different orbitals are determined from the eigenvector v χ j of the divergent eigenvalue of the 3 × 3 susceptibility matrix χ αβ (q * j ) (j = 1, 2) at the transition. The three components of v χ j give the relative (but not Thus, for q * 1 , the charge modulation amplitudes on the Cu and Op x sites are about 25% and 43% respectively of the amplitude on the Op y site. This distinction between Op x and Op y sites is consistent with the observation of anisotropic NMR linewidths in YBa 2 Cu 3 O 6.5 [53]: the linewidths of O(2) oxygen sites, which lie perpendicular to q * [29] are roughly 50% greater than for O(3) sites, which lie parallel to q * . As pointed out in [75], the ratio of the Op y to Op x modulation grows (shrinks) rapidly with increasing (decreasing) U d . Figure 1(e) illustrates the unidirectional charge modulations derived from v χ 1 . As v χ 1 directly tells, the charge modulations on the Op x and Op y orbitals are out of phase, so there is a significant intra-unit cell nematic-like charge transfer between them. The charge ordered phase is not purely nematic, however, as there are also modulations of the total charge per unit cell and of the Cu charge. This structure is consistent with the observation of nematic-like modulations of the oxygen orbitals by STM [23,32,33] and elastic RXS [28], and the observation of Cu charge modulations by NMR [19]. For a unit cell centered on a Cud orbital at r, the total charge modulation is δn tot (r) = δn Cu (r) + 1 2 δ δn p (r + δ), where r + δ are the locations of the four neighboring oxygen atoms; the nematic modulation is defined by δn nem (r) = 1 2 δ (−1) δy δn p (r + δ). Figure 1(f) clearly shows that all three types of modulation are present. These different symmetries must in fact mix because χ αβ (q) is not invariant under fourfold rotations when q = 0.
To understand the role of the Cu4s orbital, we compare our results to those for the Emery model, which does not include it.t pp = −1 eV is chosen for both models, so that the only difference between them is that the diagonal matrix elements of H(k) are unrenormalized in the Emery model. As shown in Fig. 2, this changes the Fermi-surface shape and the underlying band structure only quantitatively, with a noticeable increase of the Fermi-surface curvature. Indeed, the incommensurate peak positions q * j in the charge susceptibility shift only by about 5% between the two models for M = 1.5 eV. Surprising and important, however, is that the leading instability in the Emery model is to a q = 0 nematic phase, and that the incommensurate phase is subleading. We have traced this difference to the oxygen spectral weight distribution along the Fermi surface, which is strongly anisotropic in the Emery model, but nearly isotropic in the ALJP model (see Appendix A). Thus, the Cu4s orbital stabilizes the ALJP model against q = 0 nematic order.
To discuss the Fermi surface reconstruction from charge order we show in Fig. 3(a) the four original Fermi surface hole pockets centered at (±π/2, ±π/2), which we label (±1, ±1); these are the "nodal" pockets. Charge order along a direction q * j scatters quasiparticles through ±q * j and generates replica Fermi-surface pockets. Red contours mark those first-order replicas, generated by shifting the (−1, 1) pocket by ±q * 1 and the (1, −1) pocket by ±q * 2 , that touch the (1, 1) nodal pocket. Where original and replica pockets touch, the bands hybridize and a gap opens. Importantly, at any doping q * 1,2 are such that replica and original pockets precisely touch without crossing. We include also a second-order replica (blue dotted) by shifting the (−1, −1) pocket by q * 1 + q * 2 . This replica appears only when the order is bi-directional, and it hybridizes with two of the first order replicas and the original (1, 1) nodal hole pocket to form a diamond-shaped electron pocket shown as the grey region on the front side of the (1, 1) pocket [closest to the origin] in Fig. 3(a). It was argued empirically [36] that electron pockets of this diamond type could explain observed magneto-oscillations in YBa 2 Cu 3 O 6.5 . Yet, the interpretation is complicated because, in addition to a central frequency of F expt ∼ 530 T [76,77,78], a pair of side frequencies is observed [79]. The latter have been attributed to bilayer splitting into bonding and antibonding bands [79,36]. For the ALJP model, we find that the electron pocket has an area A 1 = 0.50/a 2 0 (a 0 is the lattice constant) which gives an oscillation frequency F 1 = ( /2πe)A 1 = 340 T, slightly less than F expt . However, since A 1 represents only ∼ 1% of the BZ area, it is far more sensitive to the Fermi surface shape than is q * . We obtain, for example, F 1 = 1000 T using the Emery model with M = 1.5eV and the parameters in the caption of Fig. 2; this is a factor of 3 larger than the ALJP result, even though the incommensurate q * differs by only ∼ 10% between the two models. Obviously, fine tuning of the ALJP model, which is based on band structure calculations for YBa 2 Cu 3 O 7 , is needed to quantitatively match quantum oscillation experiments performed on YBa 2 Cu 3 O 6.5 .
One difference to the proposal in Ref. [36] is that we find four electron pockets attached to each nodal pocket, rather than one. In addition to the electron pocket discussed above, there is a second electron pocket with identical area (not shown) on the back side of the nodal pocket [closest to (π, π)]. Two further diamond-shaped electron pockets form at opposite ends of the each nodal pocket. One of these, with an area area A 2 = 0.10/a 2 0 and corresponding oscillation frequency F 2 = 65 T, is shown as a shaded pink region in Fig. 3(a). These additional electron pockets are an artefact of the assumed infinite AF correlation length. As we said previously, when ξ AF is finite, the spectral function is characterized by Fermi arcs that resemble the front side of the nodal pockets; the back and side electron pockets only emerge as ξ AF diverges. [74] To see the effect of charge order on the spectral function, we model bi-directional charge order as a perturbation of the Cud, Op x , and Op y site energies by Adding the corresponding potential term to the Hamiltonian, we calculate the spectral function A(k, ω) = α n |φ αn (k)| 2 δ(ω − E nk ) at the Fermi energy ω = ε F , where φ αn (k) are the energy eigenvectors indicating the projection of band n onto orbital α, and E nk are the energy eigenvalues. Figure 3(b) shows A(k, ε F ) without charge order (δ = 0). In Figs. 3(c) and (d) the modulation potential is increased to δ = 0.25 eV and δ = 0.5 eV, respectively. These selected values are exaggerated for presentation purposes. The main effect of charge order is to erode spectral weight along segments of the Fermi surface that touch replicas as in Fig. 3(a). In contrast, the spectral weight is almost unaffected by charge order along short arcs on the insides of the nodal pockets. Also, the diamond-shaped electron pockets shown in Fig. 3(a) are unobservable, even for the unphysically large value of δ used in Fig. 3(d).
In our model calculations, the charge instability is driven by the Coulomb repulsion V pp between electrons on neighboring oxygen atoms. In the doping window 0.1 < p < 0.14 the ordering wavevector q * j continuously decreases with p as in the x-ray diffraction experiments by Blackburn et al. [29] [ Fig. 1(d)]. In the same doping regime the calculated charge ordering temperature T co rises with increasing p. Because the calculation of T co is numerically intensive, we show instead in Fig. 4 the inverse of the critical interaction strength, namely V −1 pp , required to drive the charge ordering transition at fixed T = 110 K. This quantity is a useful proxy for T co : a large value of V −1 pp indicates that the system is very susceptible for charge ordering, and should therefore have a large T co . In our calculations, the susceptibility towards charge order with growing hole density p concomitantly increases with the increasing size of the nodal hole-Fermi pockets.
Experimentally, the variation of T co with hole doping remains inconclusive. RXS data indicate that T co decreases with increasing p [24], but this trend is at variance with earlier x-ray data and with the field-tuned T co (H) observed by NMR [19,34]. From the latter data a maximum T co around p = 0.12 was inferred [34], and a similar dome-shaped p-dependence was determined for the Fermi-surface reconstruction from Hall measurements [19]. Recent x-ray experiments on YBa 2 Cu 3 O 6+x [49] also find a dome-shaped dependence of T co on p, peaked at p ∼ 0.10.
The evolution of T co in model calculations likely depends on the detailed doping dependence of both the effective interaction strength in the charge ordering channel, which has not been considered here, and the Fermi surface, which is the central topic of this work. A further complication is the role of disorder, which is unavoidable in doped cuprates and should influence the spatial lock-in of any charge-density wave. The issue of how T co evolves with doping is an open question that needs to be resolved.

Conclusions
In this work, we have described a model calculation that provides a route to understand the doping dependence of the charge-ordering wavevectors q * in cuprate superconductors. The essential model ingredients are a realistic multiorbital description of the CuO 2 planes, the assumption and the ansatz that strong correlation effects on the Cud x 2 −y 2 orbitals can be modeled by antiferromagnetically correlated moments, and the inclusion of short range Coulomb forces that drive the charge-ordering instability. While the model analysis is still incomplete, e.g. inelastic spin-scattering processes and the spin dynamics are neglected, it nonetheless provides an important result: quantitatively correct charge-ordering wavevectors q * are obtained, if the charge order is presumed to emerge from the pseudogap phase, rather than to generate the pseudogap itself.
Also a subtle but important role played by multiorbital physics is highlighted. While the three-orbital Emery model and the four-orbital ALJP model have similar Fermi surfaces, the leading instability in the Emery model is to a q = 0 nematic phase, while the ALJP model correctly reproduces the structure seen experimentally. This distinction is traced to subtle differences in the orbital composition of the conduction band.
A number of questions necessarily remains open, in particular the relationship between the charge order and the pseudogap and also the possible connection to the emergence of spontaneous loop currents await further clarification. Notably, the dependence of T co on p is a challenging question that demands an improved treatment of the pseudogap phase beyond the initial steps presented in this work. where s x = sin(k x /2) and s y = sin(k y /2), and wherẽ is an array of electron annihilation operators for the Cu3d x 2 −y 2 , Op x , Op y , and Cu4s orbitals, respectively. The spin index is suppressed in Eqs. (A.1) and (A.2). We can integrate out the 4s orbital in the usual downfolding procedure [72,80]. Writing the four-band Hamiltonian matrix in a block structure, where the subscript notation i × j denotes the size of each block, we solve the equationsof-motion for the Green's function in the subspace of d x 2 −y 2 , p x , and p y orbitals: From the structure of G(k, ω) at ω = ε F , an effective three-band Hamiltonian matrix is generated where t d pp is the direct hopping between p x and p y orbitals, and is the indirect hopping amplitude, through the 4s orbital, between p orbitals. Importantly, we note that ε F < s , so that Based on the signs of the orbital lobes, we would expect t d pp > 0; however, Andersen et al. proposed that t d pp is negligible compared to the indirect contribution, and that t pp ∼ −1 eV. Throughout this work, we adopt the values of t pd , t d pp , t i pp , and d − p given by ALJP [72] and listed in Table A1. Figure A1 shows the spectral functions at the Fermi energies projected onto both Cu and Op x orbitals. For comparison, results are also shown for the Emery model.

Appendix A.2. Slater Antiferromagnetism
We add a staggered magnetic field at the copper sites to the Hamiltonian to generate local moments on the Cud orbitals. Then, with the spin index included, the Hamiltonian isĤ where Q = (π, π), In the state with staggered copper moments, the wavevector k is restricted to the first antiferromagnetic (AF) BZ, labelled I in Fig. A2. Hence, k + Q belongs to the second AF BZ, labelled II in Fig. A2.

. Diagrammatic perturbation theory
We calculate the nematic susceptibility by summing the ladder and bubble diagrams shown in Fig. B1. This is analogous to what was done in Ref. [64], and we describe here how that calculation has been extended to the AF case. In Fig. B1, the wavevectors k and k + Q are constrained to the first and second AFBZs, respectively, pictured in Fig. A2, while q is unconstrained. In this notation k is conserved along each propagator and the indices 1 , 2 , . . . = 0, 1 label the AFBZ to which the creation or annihilation operators at the ends of the lines belong. For example, the line end labeled k 1 θ 1 in Fig. B1 has a corresponding annihilation operator c θ k 1 + 1 Qσ , where σ is the electron spin. Fig. B1(a) shows the bare interaction vertex V ρ (k, k , q) between charges, which includes both direct (first term) and exchange (second term) diagrams, which is Figure B1. Diagrams evaluated in the calculation of the charge susceptibility. (a) Effective interaction in the charge channel, including both Hartree (first term) and exchange (second term) contributions. The wavevector k is restricted to the first AFBZ, and the greek labels, denote the orbital type (d, p x , p y ). The labels 1 , etc. indicate the AFBZ of the corresponding electron creation or annihilation operator, which has momentum k + 1 Q. (b) Diagrams summed in the calculation of the charge susceptibility χ αβ (q).
Appendix B.2. Origin of the q = 0 instability in the Emery model A comparison between the ALJP and Emery models is made in Fig. A1. In both models, the Cu spectral weight is large and uniformly distributed along the Fermi surface.
The Op x spectral weight is comparatively weak, but because the charge instability involves primarily oxygen atoms, the details of the Op x spectral weight distribution are important. Notably, the Op x spectral weight is highly anisotropic in the Emery model and more isotropic in the ALJP model. (The Op y spectral function A py (k, ε F ) is obtained by rotating A px (k, ε F ) by 90 • .) As a consequence, the matrix element of the bare susceptibility is strongly reduced in the Emery model (the superscript 0 indicates the susceptibility in the noninteracting limit). As we show below, this matrix element tends to stabilize the system against nematic order. We focus on the nonmagnetic case where approximate analytic expressions are easily obtained. Within a simplified random phase approximation in which all interactions except V pp are ignored, we have at q = 0 which has a diverging eigenvalue when 1 + 8V pp χ 0 xy − χ 0 xx χ 0 yy = 0. (B.9) (The factor of 8 arises because of a sum over spin and over the four neighboring oxygen sites for each Op orbital.) From this equation, it is clear that χ 0 xx and χ 0 yy drive the nematic transition while χ 0 xy opposes it. Thus, it appears that the strong anisotropy of oxygen spectral weight in the Emery model is the principal difference between the Emery and ALJP models which makes the former unstable to a q = 0 nematic instability.